We present experimental scattering matrices of the JSC Mars-1, MMS-2, and MGS-1 simulants at 488 and 640 nm. The analogs were processed so that narrow size distributions representative of Martian dust aerosols during different dust cycles were obtained. We find that the forward peak of the phase function depends on particle size as it becomes narrower with increasing size, whereas the side- and backscattering directions depend on both composition and size so that increasing size and decreasing absorption produce a flatter curve. The position and maximum of the degree of linear polarization varies based on particle size and composition, and the negative polarization branch is more prominent for wavelength-scale particles diminishing with increasing size. The linear depolarization is strongly affected by size and composition. Finally, we compare sky-brightness curves measured by the Navcam and Hazcam engineering cameras on board the Mars Science Laboratory rover to the measured phase functions. The observations show a narrower peak at the forward direction and a flatter curve toward the side- and backscattering directions with an increasing dust load in the atmosphere, similar to what can be seen for the measured phase functions of the analogs with increasing particle size. In the case of the analogs, the flattening of the curve can be caused by an increase in multiple scattering within a particle by wavelength-scale surface roughness and/or internal inclusions. For the observed sky brightnesses, particle aggregation and multiple scattering among particles in denser dust conditions play a major role.