© The Authors 2025

1 Introduction

The Stokes parameters of light incident on and scattered by particles are interrelated by a 4 × 4 Mueller matrix, the scattering matrix (van de Hulst 1957; Bohren & Huffman 1998). The scattering matrix contains complete information about the scattering process; in other words, how polarized incident light is transformed into polarized scattered light. In the scattering phase matrix here termed F, a normalization is introduced by separating the integrated scattering cross section for unpolarized incident light from the scattering matrix and normalizing the F₁₁-element, the scattering phase function, to 4π in the integration over the solid angle.

A commonly encountered scattering phase matrix is blockdiagonal, resulting from averaging over a particle ensemble. Natural particle ensembles can typically be assumed to consist of particles and their mirror particles, both oriented randomly in space. In such a case, the nonzero elements of the matrix are the diagonal elements F₁₁, F₂₂, F₃₃, and F₄₄ and the off-diagonal elements F₁₂ = F₂₁ and F₃₄ = -F₄₃.

The scattering matrices play a crucial role in scattering problems related to cosmic dust. We develop an empirical parameterization for the block-diagonal scattering phase matrix based on physics and practical aspects, striving toward a user-friendly formulation for inverse light-scattering problems and other applications. Before that, we offer a concise overview of related studies.

Baes et al. (2022) provide a recent review of empirical models for the scattering phase function. They further develop a two-term ultraspherical phase function, pointing out that the new phase function improves the scattering modeling, in particular, outside the spectral range of visible light.

The Henyey-Greenstein scattering phase function is by far the most commonly used empirical phase function (Henyey & Greenstein 1941). It allows for analytical manipulations and is parameterized by the asymmetry parameter, the average for the cosine of the scattering angle, characterizing the scattering into forward- and backward-directed hemispheres. Introducing two Henyey-Greenstein terms is a popular approach as it allows for realistic increases in the phase function toward both the forward-and backward-scattering directions.

Isotropic, Legendre polynomial, and Rayleigh phase functions are among the early choices for the phase function (e.g., Chandrasekhar 1960). These are, however, quite specific and have restrictions insofar as their application to wavelength-scale particle scattering is concerned.

Pollack & Cuzzi (1980) introduced a semiempirical scattering phase function for wavelength-scale particles based on three components. Mie theory is used for particles smaller than a given threshold in terms of the size parameter x = ka, where k = 2π/λ is the wave number for the wavelength, λ, and a describes the radial size. Above the threshold, the phase function is based on three components. The first component, diffraction, is modeled by using physical optics. The second component, external reflection, derives from geometric optics. The third component involves a single-parameter exponential expression. Pollack & Cuzzi (1980) obtained good matches for the semiempirical phase function and measured or computed phase functions. Bowell et al. (1989) extended the modeling with additional exponential terms and came up with an average scattering phase function for applications in multiple scattering by planetary regoliths.

Braak et al. (2001) devised an empirical phase matrix with the Henyey-Greenstein phase function and a scalable Rayleighlike degree of linear polarization. In order to allow for mathematically correct modeling of phase matrices, Domke (1974) introduced a series presentation in generalized spherical functions (see also de Haan et al. 1987). Due to the large number of coefficients in the series presentations and the nonlocal character of the basis functions, the generalized spherical functions are an inconvenient choice for the variation of matrix elements in, for example, inverse light-scattering problems.

Cubic splines have turned out to be highly useful in the presentation of Mie scattering phase matrix elements for arbitrary scattering angles, when the exact element values are available for a uniform set of scattering angles (e.g., angular resolution of a degree; Muinonen 2004). Furthermore, splines are an excellent choice for smoothening elements, either measured or computed, that may include random errors (e.g., Muñoz et al. 2012; Muinonen et al. 1996). For splines applied to wavelengthscale particles, some 15 nodes adaptively selected across the scattering-angle range can be sufficient for an accurate presentation. Consequently, using splines would result in coefficients on the order of 10² for presenting the six matrix elements. For controlled variation in the elements, splines are not optimal due to the large number of parameters, difficulties in enforcing symmetry relations (Hovenier & van der Mee 1983; Hovenier & van der Mee 2000), and local changes producing unwanted undulations in the angular patterns.

An earlier phenomenological scattering phase matrix has been developed for radiative transfer and coherent backscattering applications to atmosphereless planetary-system objects by Muinonen & Videen (2012). In that work, the phase matrix was modeled as a sum of two pure Mueller matrices with a double Henyey-Greenstein phase function. While the phenomenological matrix has been successfully used in the polarimetric modeling of icy satellites (Ejeta et al. 2013; Kiselev et al. 2022, 2024), it suffers from a number of shortcomings. First, the double Henyey-Greenstein phase function does not allow for satisfactory presentations of phase functions measured or computed for small particles. Second, there are limitations for the angular characteristics of matrix element ratios, making the modeling of measured and computed ratios challenging.

In Section 2, we introduce the ensemble-averaged scattering phase matrix together with the symmetry relations among its elements. We introduce the parameters and basis functions for the novel empirical phase matrix and describe the parameter constraints arising from the symmetry relations. The section further presents a number of additional descriptors for the matrix elements. In Section 3, we summarize parameter optimization methods and present Markov chain Monte Carlo (MCMC) methods for the evaluation of plausible parameter ranges. Section 4 focuses on the modeling of experimental and theoretical scattering phase matrices by using the empirical system. The set of examples serves to illustrate how the free parameters affect the shapes of the matrix elements. In Section 5, we close the study with conclusions and future prospects.

2 Scattering matrix

2.1 Ensemble-averaged scattering phase matrix

Consider a particle ensemble that includes equal numbers of nonspherical particles and their mirror particles, both in random orientation. The average 4 × 4 scattering phase matrix, F = F(θ), of the ensemble relates the Stokes parameters of incident light, I_i = (I_i, Q_i, U_i, V_i)^T, to those of scattered light, I_s = (I_s, Q_s, U_s, V_s)^T (θ is the scattering angle, a spherical polar angle with θ = 0 refering to the forward-scattering direction; Bohren & Huffman 1998), $${\mathbf{I}}_{\mathrm{s}}(\theta )=\frac{{\sigma }_{\mathrm{s}}}{4\pi r^2}\mathbf{F}(\theta )\cdot {\mathbf{I}}_{\mathrm{i}}.$$(1)

Here σ_s and r denote the scattering cross section for unpolarized incident light and the distance between the scatterer and the observer, respectively. The scattering phase matrix takes the block-diagonal form (Hovenier & van der Mee 2000) $$\mathbf{F}(\theta )=\left(\begin{array}{cccc}{a}_1(\theta )& b_1(\theta )& 0& 0\\ b_1(\theta )& a_2(\theta )& 0& 0\\ 0& 0& a_3(\theta )& b_2(\theta )\\ 0& 0& -b_2(\theta )& a_4(\theta )\end{array}\right),$$(2)

where we have the common normalization of $${\int }_{(4\pi )}\frac{d\mathrm{\Omega }}{4\pi }{a}_1(\theta )=\frac{1}{2}{\int }_0^{\pi }d\theta \sin \theta a_1(\theta )=1.$$(3)

We assessed the functions b₁, a₂, a₃, b₂, and a₄ relative to a₁ and, for convenience, defined $$\begin{array}{r}{a}_{j\mathrm{r}}(\theta )=\frac{a_j(\theta )}{a_1(\theta )},j=2,3,4,\\ b_{j\mathrm{r}}(\theta )=\frac{b_j(\theta )}{a_1(\theta )},j=1,2.\end{array}$$(4)

The scattering phase matrix is subject to symmetry relations (Hovenier & van der Mee 2000). For arbitrary scattering angles, $$\begin{gathered} |a_{j\mathrm{r}}(\theta )|\le 1,j=2,3,4, \\ |b_{j\mathrm{r}}(\theta )|\le 1,j=1,2, \\ [a_{3\mathrm{r}}(\theta )+a_{4\mathrm{r}}(\theta ){]}^2+4b_{2\mathrm{r}}^2(\theta )\le [1+a_{2\mathrm{r}}(\theta ){]}^2-4b_{1\mathrm{r}}^2(\theta ), \\ |a_{3\mathrm{r}}(\theta )-a_{4\mathrm{r}}(\theta )|\le 1-a_{2\mathrm{r}}(\theta ), \\ |a_{2\mathrm{r}}(\theta )-b_{1\mathrm{r}}(\theta )|\le 1-b_{1\mathrm{r}}(\theta ), \\ |a_{2\mathrm{r}}(\theta )+b_{1\mathrm{r}}(\theta )|\le 1+b_{1\mathrm{r}}(\theta ). \end{gathered}$$(5)

In the forward and backward directions, $$\begin{gathered} a_{2\mathrm{r}}(0)=a_{3\mathrm{r}}(0), \\ {a}_{2\mathrm{r}}(\pi )=-a_{3\mathrm{r}}(\pi ), \\ {a}_{4\mathrm{r}}(\pi )=1-2a_{2\mathrm{r}}(\pi ), \\ {b}_{1\mathrm{r}}(0)=b_{2\mathrm{r}}(0)=b_{1\mathrm{r}}(\pi )=b_{2\mathrm{r}}(\pi )=0. \end{gathered}$$(6)

2.2 Empirical system

We introduced a parameterized scattering phase matrix that captures the overall angular characteristics of typical ensemble-averaged experimental and theoretical phase matrix elements. The matrix elements are described as follows for arbitrary scattering angles, θ, through the dependence of $q=\sin \frac{1}{2}\theta$: $$\begin{gathered} a_1=\frac{D_1}{1+d_1{q}^{{\delta }_1}}+\frac{R_1}{1+r_1{p}^{{\rho }_1}}+S_1\exp (-s_1{p}^{{\sigma }_1}), \\ {a}_{j\mathrm{r}}=D_j{p}^{{\delta }_j}+R_j{q}^{{\rho }_j}+S_j\exp (-s_j{p}^{{\sigma }_j}),j=2,3,4, \\ {b}_{j\mathrm{r}}=A_j{p}^{{\beta }_j}{q}^{{\gamma }_j}+H_j{p}^{{\beta }_j}[1-\exp (-h_j{q}^{{\eta }_j})] \\ +K_j{q}^{{\gamma }_j}[1-\exp (-k_j{p}^{{\kappa }_j})],j=1,2, \end{gathered}$$(7)

where $$\begin{gathered} p=\cos \frac{1}{2}\theta =\sqrt{1-q^2}, \\ q=\sin \frac{1}{2}\theta . \end{gathered}$$(8)

We note that the system parameters are relevant only in the regime of parameter phase space, where the symmetry relations of Eqs. (5) and (6) are fulfilled. The parameterization is nonlinear and it is possible to arrive at equally satisfactory representations of the scattering matrix elements in different parts of the parameter space.

In the most general form in Eq. (7), there are 48 parameters describing six functions of θ. For a₁ as a rule of thumb, the parameter triplets D₁, d₁, δ₁ and R₁, r₁, ρ₁ characterize the angular dependencies in the forward-scattering (“D” for “drive”) and backward-scattering (“R” for “reverse”) regimes, whereas the triplet S, s₁, σ₁ describes backscattering characteristics for narrow to moderate angular ranges (“S” for “surge”). For a_jr(j = 2, 3,4), the doublet D_j, δ_j and the quintet R_j, ρ_j, S_j, S_j, σ_j refer to the forward- and backward-scattering regimes, respectively. For b_jr (j = 1,2), the triplet A_j, βj, γ_j refers to the overall angular characteristics, whereas the triplets H₁, h₁, η₁ and K₁, k₁, κ₁ refer to the forward- and backward-scattering regimes, respectively. Taking into account the condition of a₁ in Eq. (3) and the two conditions for a_jr (j = 2, 3,4) in Eq. (6), the number of free parameters reduces to 44. Furthermore, the triplets S_j, S_j, σ_j (j = 1,2, 3,4) and the parameters η_j and κ_j (j = 1,2) do not need to be always utilized, so the number of free parameters reduces to 28, which implies, on average, fewer than five free parameters for each of the six functions.

The a₁ model in Eq. (7) emerges from the Q-factor algebra (Sorensen et al. 2017) in the forward-scattering regime (see also van de Hulst 1957). Sorensen et al. (2017) show that the scattering phase functions, a₁, of varying small particle ensembles show a linear dependence on the logarithmic scale in the forward-scattering regime, when the phase functions are expressed as a function of $Q=k\sin \frac{1}{2}\theta$. In Eq. (7), a regularization of the linear dependence is introduced with the help of a rational function to remove the singularity in the forward-scattering direction. Furthermore, the wave number, k, does not enter the a₁ model (thus, $q=\sin \frac{1}{2}\theta$). It is clear that alternative basis functions can be introduced, when the particular forward-scattering characteristics are of interest (for example, the Fraunhofer diffraction pattern). The a₁ model includes an analogous dependency in the backward-scattering regime, suggested by the forward- and backward-propagating internal electromagnetic waves within small scattering particles, giving rise to similarities in forward- and backward-scattered waves (Muinonen et al. 2011). Finally, in order to account for a backward-scattering surge often detected for small particles (e.g., Lumme & Rahola 1998; Muñoz et al. 2012), an additional exponential-trigonometric term is introduced.

As was described above for a_jr (j = 2,3,4), there are formally seven parameters for each (the D_j, δ_j and R_j, ρ_j doublets and the S_j, S_j, σ_j triplets), amounting to 21 parameters in total. The symmetry relations in Eq. (6) introduce conditions that reduce the number of free parameters to 18: $$\begin{gathered} a_{3\mathrm{r}}(0)=a_{2\mathrm{r}}(0), \\ {a}_{3\mathrm{r}}(\pi )=-a_{2\mathrm{r}}(\pi ), \\ {a}_{4\mathrm{r}}(\pi )=1-2a_{2\mathrm{r}}(\pi ). \end{gathered}$$(9)

With the help of Eqs. (7) and (6), we express D₃, R₃, and R₄ with the help of D₂, R₂, S₂, s₂, S₃, s₃, and S₄: $$\begin{gathered} D_3=D_2+S_2\exp (-s_2)-S_3\exp (-s_3), \\ {R}_3=-R_2-S_2-S_3, \\ {R}_4=1-2(R_2+S_2)-S_4. \end{gathered}$$(10)

For the matrix element ratios b_1r and b_2r, there are nine free parameters for each. In addition to the argumentation above, the model is motivated by the trigonometric and linear-exponential models for polarimetric phase curves of atmosphereless Solar System objects (Lumme & Muinonen 1993; Penttilä et al. 2005; Muinonen et al. 2009). There are, however, differences: the present system allows for vanishing first derivatives in the backward- and forward-scattering directions and is thus applicable to scattering matrices of small particles.

2.3 Additional descriptors

As the exponents of the empirical system in Eq. (7) relate to the physics of scattering, they serve as useful descriptors of the angular dependencies. However, due to the complex parameter interplay, additional descriptors are called for to characterize the dependencies.

For the scattering phase function a₁, the norm is by definition N₁₁ = 1 and the asymmetry parameter g₁₁ describes the balance between scattering into forward and backward hemispheres (to distinguish between the descriptors, a numerical subscript of the matrix element in Eq. (2) was used), $$\begin{gathered} N_{11}={\int }_{(4\pi )}\frac{d\mathrm{\Omega }}{4\pi }{a}_1(\theta )=\frac{1}{2}{\int }_0^{\pi }d\theta \sin \theta a_1(\theta )=1, \\ {g}_{11}={\int }_{(4\pi )}\frac{d\mathrm{\Omega }}{4\pi }\cos \theta a_1(\theta )=\frac{1}{2}{\int }_0^{\pi }d\theta \sin \theta \cos \theta a_1(\theta ). \end{gathered}$$(11)

The predicted values in the forward and backward directions are of interest, as is the minimum value, a_1↓ and the corresponding scattering angle, θ_11↓: $$\begin{gathered} a_1(0)=D_1+\frac{R_1}{1+r_1}+S_1\exp (-s_1), \\ {a}_1(\pi )=\frac{D_1}{1+d_1}+R_1+S_1, \\ {a}_{1\downarrow }=a_1({\theta }_{11\downarrow }). \end{gathered}$$(12)

Note that a₁(0) and a₁(π) are independent of the exponents δ₁, ρ₁, and σ₁.

For the off-diagonal element b₁, N₁₂ and g₁₂ (Muinonen 2004), analogous to the parameters for a₁, describe the polarizing efficiency and asymmetry, respectively, $$\begin{gathered} N_{12}=-{\int }_{(4\pi )}\frac{d\mathrm{\Omega }}{4\pi }{b}_1(\theta )=-\frac{1}{2}{\int }_0^{\pi }d\theta \sin \theta b_1(\theta ), \\ {g}_{12}=-{\int }_{(4\pi )}\frac{d\mathrm{\Omega }}{4\pi }\cos \theta b_1(\theta )=-\frac{1}{2}{\int }_0^{\pi }d\theta \sin \theta \cos \theta b_1(\theta ). \end{gathered}$$(13)

Natural descriptors of −b_1r are the minimum and maximum degrees of polarization, −b_1r↓ and −b_1r↑, as well as the corresponding scattering angles, θ_12↓ and θ_12↑, together with the inversion angle, θ_12i, $$\begin{gathered} -b_{1\mathrm{r}\downarrow }=-b_{1\mathrm{r}}({\theta }_{12\downarrow }), \\ -b_{1\mathrm{r}\uparrow }=-b_{1\mathrm{r}}({\theta }_{12\uparrow }). \\ {b}_{1\mathrm{r}}({\theta }_{12\mathrm{i}})=0. \end{gathered}$$(14)

If there are two negative branches present, one near backward scattering and another near forward scattering, they can both be characterized similarly. For the off-diagonal element b_2r, the characterization is analogous.

For a_2r, a_3r, and a_4r, the values at forward- and backwardscattering directions are of interest, $$\begin{gathered} a_{2\mathrm{r}}(0)=a_{3\mathrm{r}}(0)=D_2+S_2\exp (-s_2), \\ {a}_{4\mathrm{r}}(0)=D_4+S_4\exp (-s_4), \\ {a}_{2\mathrm{r}}(\pi )=-a_{3\mathrm{r}}(\pi )=R_2+S_2, \\ {a}_{4\mathrm{r}}(\pi )=1-2a_{2\mathrm{r}}(\pi )=R_4+S_4=1-2R_2-2S_2. \end{gathered}$$(15)

The minimum and maximum values of these elements (not repeating the values in the forward- and backward-scattering directions) and the corresponding scattering angles can be used as descriptors, $$\begin{gathered} a_{j\mathrm{r}\downarrow }=a_{j\mathrm{r}}({\theta }_{jj\downarrow }), \\ {a}_{j\mathrm{r}\uparrow }=a_{j\mathrm{r}}({\theta }_{jj\uparrow }),j=2,3,4. \end{gathered}$$(16)

For the elements a_3r and a_4r, the inversion angles θ_33i and θ_44i (when relevant) can be utilized, $$a_{j\mathrm{r}}({\theta }_{jj\mathrm{i}})=0,j=3,4.$$(17)

3 Numerical methods

Least-squares optimization of the parameters was carried out by using the Nelder-Mead downhill simplex method (Press et al. 1992) and parameter uncertainties were evaluated with MCMC methods (cf., O’Hagan & Forster 2004). Note that, since the system is empirical in its nature, the parameters and their uncertainties must be treated with caution rather than related to the first principles of physics. Throughout the numerical methods, the sets of parameters to be incorporated in the analyses can be selected by the user.

Given the analytical expressions for the matrix elements, simple Matlab routines have been developed for the initial estimation of the parameters, complemented with a systematic mapping of the parameter space by using a Fortran implementation. While mapping the parameter space, symmetry relations are monitored and enforced. Least-squares optimization is then initialized in the proximity of the plausible solution, substantially improving the prospects for success in detecting the global and realistic extremum.

We followed Bayesian inference in defining a posteriori probability density function (p.d.f.), p_p(P), for the free parameters (unknowns), P (e.g., Muinonen et al. 2020), by assuming a constant a priori p.d.f. We utilized the Metropolis-Hastings algorithm with a symmetric Gaussian proposal p.d.f., p_t(P; P″) (the subscript “t” stands for transition). Here p_t is a product of uncorrelated monovariate Gaussian densities. The MCMC algorithm proceeds by the computation of the ratio a_r: $$a_r=\frac{p_{\mathrm{p}}({\mathbf{P}}^{\prime })p_{\mathrm{t}}({\mathbf{P}}_j;{\mathbf{P}}^{\prime })}{p_{\mathrm{p}}({\mathbf{P}}_j)p_{\mathrm{t}}({\mathbf{P}}^{\prime };{\mathbf{P}}_j)}.$$(18)

Here P_j and P″ denote the current and proposed parameters in a Markov chain, respectively. The proposed parameters, P″, are accepted or rejected with the help of a uniform random deviate, y ∈ [0,1]: $${\mathbf{P}}_{j+1}=\{\begin{array}{ll}{\mathbf{P}}^{\prime },& y\le a_r\\ {\mathbf{P}}_j,& y>a_r\end{array}.$$(19)

The proposed parameters are thus accepted with the probability of min(1, a_r). After a number of transitions in the burn-in phase, the Markov chain, in the case of success, converges on sampling the target p.d.f. p_p.

In our practical implementation, the MCMC part can first be carried out separately for all the matrix elements, completing the analysis for a₁ that is independent of the remaining matrix element ratios. Thereafter, the final MCMC simulation concerns simultaneously the five element ratios a_2r, a_3r, a_4r, b_1r, and b_2r. The root-mean-square (rms) values, σ_rms, of the least-squares fits are computed from $${\sigma }_{\mathrm{r}\mathrm{m}\mathrm{s}}=\sqrt{\frac{1}{N_{\mathrm{o}\mathrm{b}\mathrm{s}}}\sum _{j=1}^{N_{\mathrm{o}\mathrm{b}\mathrm{s}}}{[m_j-M({\theta }_j,\mathbf{P})]}^2},$$(20)

where N_obs denotes the number of observations, m_j and M(θ_j, P) stand for the observed and computed element values for a_1r, a_2r, a_3r, a_4r, −b_1r, and b_2r, respectively, and θ_j denote the scattering angles. Here the element a_1r is logarithmic and corresponds to $$a_{1\mathrm{r}}=\lg a_1.$$(21)

Furthermore, the restricted nature of the empirical model in Eq. (7) results in systematic differences between the observed scattering phase matrix elements and their model counterparts. We mitigated the systematics by introducing a so-called slip factor in the uncertainty model; that is, by increasing the formal standard deviation computed from the rms value of the leastsquares fit with a factor based on reducing the correlation from an expected high initial value to a targeted low final value. The procedure corresponds to adding a sufficient amount of fictitious white noise into the uncertainty model to mitigate the systematics. In what follows, we set the initial value of correlation to 0.9 and targeted the final value to be 0.1.

4 Results and discussion

We present example modeling results for polarimetric measurements obtained from the updated Granada-Amsterdam Light Scattering Database (Muñoz et al. 2025). Because measurements near the exact forward- and backward-scattering directions are typically subject to larger uncertainties, we limited the use of the measured data to the angular range between 10° and 173°. Furthermore, only those scattering angles for which all six Mueller matrix elements are available were retained. The exception was the numerical scattering simulations for water ice particles, in which case the full angular range of the data was used, as all matrix elements are available at all angles. The ice particle example further serves to illustrate that there are limitations in the use of the empirical scattering matrix system.

Across all samples, the empirical system demonstrates a consistent performance, reproducing the measured scattering behavior while adhering to symmetry constraints. Residual correlations were mitigated by introducing a sufficient amount of fictitious white noise during the data processing. For conciseness, a posteriori p.d.f.s of the parameters were characterized by the means and the standard deviations from the MCMC solutions and by the uncertainty envelopes for a number of example scattering matrices. The model parameters and their uncertainties expressed as standard deviations are summarized in Tables A.1 and A.2.

In order to supplement the model parameters, we incorporated key additional descriptors for the elements a₁ and −b_1r, which offer an enhanced physical insight into the scattering behavior. For the scattering phase function, a₁, the asymmetry parameter, g₁₁, is detailed. Regarding the degree of linear polarization, −b_1r, we provide the polarizing efficiency, N₁₂, the associated asymmetry parameter, g₁₂, and the maximum and minimum polarization degrees, −b_1r↑ and −b_1r↓. The remaining supplementary descriptors, which include extrema and inversion angles for a_2r, a_3r, a_4r, and b_2r, are compiled into Tables A.3 and A.4 as mean values derived from the MCMC analyses, accompanied by their respective standard deviations.

Fig. 1 Scattering phase matrix elements for the Allende meteorite sample: measurements (circles; Granada-Amsterdam Light Scattering Database) and modeling (solid line; present empirical system) at the wavelengths of 441.6 nm (blue) and 632.8 nm (red).

4.1 Allende meteorite

To demonstrate the capability of the parameterization, we begin with the measurement of a crushed sample of the Allende meteorite (Muñoz et al. 2000). The sample in the Granada-Amsterdam Light Scattering Database has a refractive index of 1.65 + i0.001, an effective particle size of r_eff = 0.8 μm, and an effective variance of v_eff = 3.3, indicating that the sample consists of small particles with a wide size distribution. The effect of the size distribution can be seen in Fig. 1, where only a slight wavelength dependency is shown, as well as in Tables A.1 through A.4. The Allende sample is a carbonaceous chondrite and shows strong forward scattering with high a₁(0) and gn, but a low polarizing efficiency, with shallow −b_1r extrema and low N₁₂. Figure 1 shows the measured data (36 data points in both the 441.6 nm and 632.8 nm datasets) alongside the modeling results. The rms values based on Eq. (20) for the elements lg a₁, −b_1r, a_2r, a_3r, b_2r, and a_4r are shown in Table A.1. In general, the modeling captures the features of the measured data well. The rms values are better than 1.2% for the 441.6 nm dataset and better than 2.0% for the 632.8 nm dataset. However, in the case of −b_1r, the analysis reveals slight differences at small scattering angles. It remains unclear whether these differences arise from limitations in the system or from uncertainties in the measurements themselves. Near the backward scattering direction, the linear polarization does not decrease as much as in the measured data. These differences most likely result from the small number of measurement data near the forward- and backward-scattering directions, causing the other directions to have excessive weight in the modeling.

Fig. 2 As in Fig. 1, but for the lunar surface analog (JSC-1A) at wavelengths of 488 nm (blue), 520 nm (green), and 647 nm (red).

4.2 Lunar analog

According to the rms values presented in Table A.2 for the lunar analog datasets (Escobar-Cerezo et al. 2018), modeling across all three wavelengths (488, 520, and 647 nm, each comprising 36 data points) generally yields a good agreement with the measurements. Nevertheless, the data corresponding to 488.0 nm pose certain challenges, exhibiting subtle differences in the angular dependence of the matrix elements (Fig. 2), which are likely to be attributable to noise within the measured data. The rms values are 0.41-4.6%, 0.49-2.7%, and 0.46-1.4% for the 488 nm, 520 nm, and 647 nm datasets, respectively (Table A.2). The database provides a refractive index for the sample akin to that of Allende, estimated to be 1.65 + i0.003, with r_eff = 10.50 1 μm and v_eff = 1.69, indicating a sample of large particles with a comparatively wide size distribution.

Overall, the findings demonstrate a satisfactory agreement for additional descriptors, including the minimum and maximum degrees of polarization, as well as the inversion angle. Symmetry constraints require the a_4r element to display a pronounced peak near the backward scattering direction, where experimental data are unavailable. Further investigation is warranted to ascertain whether this phenomenon arises as an artifact from the fitting process or constitutes an intrinsic characteristic of polarization. We consider it to be more likely the latter, as analogous surges result in theoretical studies. The lunar analog sample demonstrates stable and well-developed polarimetric behavior across the wavelengths examined. Its −b_1r↓ and −b_1r↓ values consistently exhibit robustness, and its a_2r↓ and a_4r↓ extrema remain relatively constant, implying scattering behavior insensitive to wavelength, which is likely to result from self-similarity due to a broad distribution of scatterer sizes.

Fig. 3 As in Fig. 1, but for the Martian surface analog (MMS-2L) at the wavelengths of 488 nm (blue) and 640 nm (red).

4.3 Martian analog

For the Martian analog, we present results at the wavelengths of 488 and 640 nm (Martikainen et al. 2024, 2025) with 32 and 35 points in the respective datasets. This red-brown-colored sample has re_eff = 16.2 μm, a small effective variance, v_eff = 0.12, and a refractive index of 1.5 + i0.0011 at 488 nm and 1.5 + i0.0035 at 640 nm. It offers an illustrative case of wavelength-dependent scattering behavior (Fig. 3).

The empirical system reproduces the measured scattering matrix elements overall with rms values of 0.45-3.4% for the 488 nm dataset and 0.39-2.6% for the 640 nm dataset (Table A.2). The 488 nm data exhibit slight deviations in the off-diagonal elements, −b_1r and b_2r, where the former shows a less pronounced negative branch compared to the measurements. Additionally, the a_4r element deviates slightly from the measured values, which could be attributed to noise in the measurements for b_2r.

Similar to the lunar analog, the Martian analog displays a significant surge in the a_4r element near the backward scattering direction, where measurement data are lacking. The sample displays a contrasting wavelength dependence. At 640 nm compared to 488 nm, MMS-2L shows a marked decrease in forwardscattering intensity and polarization amplitude, as is indicated by the decreasing trends in g₁₁, − b_1r↑, and the overall spread of Mueller matrix descriptors. The spectral sensitivity of its scattering properties distinguishes MMS-2L from the more stable silicate-rich materials in the dataset.

4.4 Feldspar

Feldspar, measured at 441.6 nm and 632.8 nm by Volten et al. (2001), represents the second case of wavelength-dependent measurements with 36 points in both datasets. The sample has r_eff = 1.0 μm and v_eff = 1.0, and the refractive index is estimated to have a real part of 1.5 to 1.6 and an imaginary part between 10⁻⁵ and 10⁻³. As is shown in Fig. 4, the modeling closely follows the measured data for both wavelengths. The rms values are small, at 0.24-1.0% for the 441.6 nm dataset and 0.28-0.88% for the 632.8 nm dataset (Table A.1). The differences in linear polarization, especially the amplitude and location of the negative polarization extremum, are well captured by the empirical system. However, at 632.8 nm there is a mild underestimation of the degree of linear polarization (−b_1r) near the backward scattering direction.

Feldspar exhibits strong forward scattering and a high polarizing efficiency, particularly at longer wavelengths. In Fig. 4, the effect of the particle size distribution is clearly visible, when the wavelength changes from 441.6 nm to 632.8 nm. This spectral change is not limited to amplitude: the spread between a_2r↓ and a_2r↑ is consistent with small particle size and narrow size distribution.

Fig. 4 As in Fig. 1, but for feldspar at the wavelengths of 441.6 nm (blue) and 632.8 nm (red).

4.5 Hematite and rutile

Hematite and rutile are samples featuring high refractive indices, small effective particle sizes, and narrow size distributions. For hematite, the size distribution is described by r_eff = 0.4 μm and v_eff = 0.6 and the complex refractive index is 3.0 + i[0.01-0.1]. The rutile sample has r_eff = 0.13 μm and v_eff = 0.5, with the real part of the refractive index 2.73 and the imaginary part near-zero. Both datasets correspond to a wavelength of 632.8 nm and include 36 data points (Shkuratov et al. 2004; Muñoz et al. 2006).

The polarization features are irregular for both samples, with a multilobed degree of linear polarization, −b_1r(θ). These features likely stem from the high refractive indices, optical anisotropy, or fine internal structures. They show pronounced matrix elements that are characterized by more complex angular dependencies, particularly in the off-diagonal matrix elements. The empirical system successfully reproduces these features (see Table A.1). For hematite, the rms values range across 0.161.2% and, for rutile, across 0.23-2.6%. Additionally, we present 10 000 MCMC results, visualized as uncertainty envelopes (hematite and rutile in Figs. 5 and 6, respectively). The results indicate that the model captures the dominant trends well and produces reliable outcomes. However, the lack of measurements near the backward-scattering and forward-scattering directions introduces uncertainties in these regions.

Fig. 5 Scattering phase matrix elements for the hematite sample: measurements (circles; Granada-Amsterdam Light Scattering Database) and modeling (solid line; empirical system). The uncertainty envelope (orange) depicts the minimum and maximum values from 10 000 MCMC sample solutions.

Fig. 6 As in Fig. 5, but for rutile.

4.6 Water ice

The scattering phase matrix for water ice, as presented in Muinonen & Markkanen (2023) with 181 points in the dataset, represents the most challenging scenario for the empirical system, containing more features than the measured data available. The fitting results, in conjunction with lg a_ι, a_3r, and a_4r, nevertheless demonstrate a reasonable fit, with rms values of 2.2-4.9% (Table A.2).

Conversely, −b_1r provides an adequate fit solely for scattering angles surpassing 50 degrees, while tending to average out prominent features below this threshold. The rms value is still excellent, at 1.1%.

In contrast, for b_2r, there are difficulties in reaching values as high as those anticipated by theoretical data near the backwardscattering direction, yet an otherwise satisfactory fit is achieved with an rms value of 2.2%. Meanwhile, a_2r corresponds well to the data near both forward and backward scattering, whereas the system falls short at intermediate scattering angles, with the resulting rms value as high as 7.6%.

Simulated water ice stands out with the most extreme forward scattering (highest g₁₁), a broad negative polarization structure, and a distinct backward-scattering peak in a_4r. Water ice shows unique behavior when compared with other samples, reinforcing its unique role as a reference for high-albedo surfaces such as those found on icy moons of giant planets or Kuiper Belt objects.

Fig. 7 Scattering phase matrix elements computed and ensemble-averaged for Gaussian-random-sphere particles composed of water ice (circles) and as described by the empirical system (solid line).

4.7 Comparative analysis with empirical descriptors

Figures 8 and 9 present a descriptor-based comparative analysis of scattering behavior across an example set of materials, derived by using the empirical scattering matrix system. The figures visualize physically meaningful quantities computed from the modeled matrix elements, such as the asymmetry parameter, gu, linear polarization extrema, −b_1r↓ and −b_1r↑, and various angular characteristics of a_2r, a_3r, and a_4r.

Figure 8 reveals a correlation between δ₁ and g₁₁ (top left). Across the sample set, higher values of δ₁, which correspond to narrower forward-scattering lobes in the model phase function, generally lead to larger g₁₁, consistent with enhanced forward-directed scattering. This trend is particularly evident for weakly absorbing silicate and icy materials such as feldspar and water ice. However, the correlation weakens for samples such as hematite and rutile, whereas, while δ₁ and g₁₁ are not directly linked in all cases, their joint behavior captures important distinctions in scattering.

The parameters δ₁ and ρ₁ (Fig. 8, top right) control the angular confinement of the forward- and backward-scattering lobes, respectively, in the model phase function a₁. The lack of a clear correlation between these descriptors indicates that the scattering mechanisms responsible for forward and backward scattering are largely different across the sample set. This separation highlights the flexibility of the empirical system in accommodating complex scattering behaviors present in planetary regolith analogs.

A secondary trend appears between the polarizing efficiency, N₁₂ - that is, the integrated strength of linear polarization across all angles - and the polarization asymmetry parameter, g₁₂ (Fig. 8, bottom left). This indicates that materials with higher overall linear polarization tend to exhibit polarization curves that are directionally biased. This relationship suggests that the angular structure and strength of linear polarization are connected. It is interesting that both N₁₂ and g₁₂ span positive and negative values. Based on N₁₂, this refers to the scatterers being net positively or negatively polarizing. For the former cases, g₁₂ is systematically positive, referring to the tendency of stronger positive polarization in the forward-scattering hemisphere and stronger negative polarization in the backward-scattering hemisphere. In contrast, for the latter cases, g₁₂ shows both negative (ice) and positive values (rutile). For ice, the negative value derives from strong positive polarization in the backward-scattering hemisphere and negative polarization in the forward-scattering hemisphere (Fig. 7). In comparison, for rutile, the angular ranges of positive and negative polarization are exchanged (Fig. 6). Notice that N₁₂ and g₁₂ relate to −b ₁ rather than −b_1r.

In Fig. 8 (bottom right), many samples with deep negative polarization minima (−b_1r↓) span a broad range of positive polarization maxima (−b_1r↑), indicating that, while both features often co-occur, their amplitudes can vary independently. This spread likely testifies to the combined influence of sample properties. The clustering in values of −b_1r↓ however, should be studied further in future work as it may suggest shared physical mechanisms underpinning both polarization extrema, with certain materials maximizing the angular polarization contrast more effectively than others.

Figure 9 extends the analysis to the angular behavior of the remaining matrix element ratios. It begins by comparing the minimum and maximum values of a_2r for those samples where both extrema are present, illustrating the angular modulation of linearly polarized incident light by the samples (top left). In Fig. 9 (top right), a_3r(π) is shown versus a_3r(0), illustrating the expected antisymmetric behavior of the element, consistent with theoretical symmetry relations for scattering by randomly oriented particles (Hovenier & van der Mee 2000). Finally (Fig. 9, bottom left and right), relationships in a_4r↓ vs. a_4r(0) and the extrema b_2r↕ of |b_2r| with the corresponding angle, θ_34r↕, suggest that strong backward-scattering features in polarization are consistent across silicate-rich materials.

A comparative look at multiple wavelengths reveals two cases with distinct spectral trends. Feldspar (441.6 nm vs. 632.8 nm) shows increased a₁(0) and −b_1r↑ at longer wavelengths, indicating enhanced forward scattering and polarization. It also shows a broader spread between minimum and maximum values for the red wavelength (632.8 nm) compared to the blue one (441.6 nm), indicating that the angular modulation of linear polarization retention becomes more pronounced with wavelength. This suggests that, in addition to amplitude scaling, feldspar’s scattering matrix undergoes subtle wavelengthdependent reshaping. In contrast, the Martian analog MMS-2L displays a suppression of both forward scattering and polarization at 640 nm compared to 488 nm. The shift is evident in both descriptor values and matrix element behavior pointing to decreased absorption at longer wavelengths. The spectral sensitivity of MMS-2L’s matrix elements further distinguishes it from silicate-rich samples.

Fig. 8 Additional descriptors derived from the parameterized scattering phase matrix elements a1 and b1 for selected materials: the asymmetry parameter, g11, versus the forward-scattering exponent, δ1 (top left), the backward-scattering exponent, ρ1, versus δ1 (top right), the polarization asymmetry parameter, g12, versus the polarizing efficiency, N12 (bottom left), and the minimum degree of linear polarization, -b1r↓, versus the corresponding maximum, -b1r↑ (bottom right).

Fig. 9 Angular characteristics of selected scattering phase matrix element ratios for the sample materials: the minimum versus maximum values of a2r (top left), a3r in the backward-scattering direction versus the same element ratio in the forward-scattering direction (top right), the minimum value of a4r versus the same ratio in the forward-scattering direction (bottom left), and the extreme value b2r↕ in the absolute sense versus the corresponding scattering angle, θ34↕ (bottom right).

5 Conclusions

We have provided a parameterized Mueller matrix for light scattering by cosmic dust particles. The resulting scattering phase matrix is analytical, user-friendly, and gives accurate representations of experimental and theoretical phase matrices. The empirical matrix comes with initialization, least-squares optimization, and MCMC schemes for the estimation of the free parameters. Hence, we have shown that it is possible to parameterize scattering phase matrices and simultaneously fulfil the symmetry relations among the matrix elements. We have illustrated the concept by applying the framework to 12 different cases of scattering phase matrices.

There are significant future prospects for the empirical scattering matrix system. First, we aim to analyze all the scattering matrices in the Granada-Amsterdam Light Scattering Database (Muñoz et al. 2012) and unveil possible trends in the parameters. Second, we expect to utilize the empirical system in inverse multiple scattering problems concerning the disk-integrated photometric and polarimetric phase curves of atmosphereless Solar System objects. Third, it may be possible to further improve the basis functions with the help of an artificial intelligence approach.

Acknowledgements

Research supported by the Research Council of Finland through grants nos. 336546 and 359893.

References

Appendix A Auxiliary tables

All Tables

All Figures

	Fig. 1 Scattering phase matrix elements for the Allende meteorite sample: measurements (circles; Granada-Amsterdam Light Scattering Database) and modeling (solid line; present empirical system) at the wavelengths of 441.6 nm (blue) and 632.8 nm (red).
In the text

	Fig. 2 As in Fig. 1, but for the lunar surface analog (JSC-1A) at wavelengths of 488 nm (blue), 520 nm (green), and 647 nm (red).
In the text

	Fig. 3 As in Fig. 1, but for the Martian surface analog (MMS-2L) at the wavelengths of 488 nm (blue) and 640 nm (red).
In the text

	Fig. 4 As in Fig. 1, but for feldspar at the wavelengths of 441.6 nm (blue) and 632.8 nm (red).
In the text

	Fig. 5 Scattering phase matrix elements for the hematite sample: measurements (circles; Granada-Amsterdam Light Scattering Database) and modeling (solid line; empirical system). The uncertainty envelope (orange) depicts the minimum and maximum values from 10 000 MCMC sample solutions.
In the text

	Fig. 6 As in Fig. 5, but for rutile.
In the text

	Fig. 7 Scattering phase matrix elements computed and ensemble-averaged for Gaussian-random-sphere particles composed of water ice (circles) and as described by the empirical system (solid line).
In the text

Fig. 8 Additional descriptors derived from the parameterized scattering phase matrix elements a1 and b1 for selected materials: the asymmetry parameter, g11, versus the forward-scattering exponent, δ1 (top left), the backward-scattering exponent, ρ1, versus δ1 (top right), the polarization asymmetry parameter, g12, versus the polarizing efficiency, N12 (bottom left), and the minimum degree of linear polarization, -b1r↓, versus the corresponding maximum, -b1r↑ (bottom right).

In the text

Fig. 9 Angular characteristics of selected scattering phase matrix element ratios for the sample materials: the minimum versus maximum values of a2r (top left), a3r in the backward-scattering direction versus the same element ratio in the forward-scattering direction (top right), the minimum value of a4r versus the same ratio in the forward-scattering direction (bottom left), and the extreme value b2r↕ in the absolute sense versus the corresponding scattering angle, θ34↕ (bottom right).

In the text