Modeling Linear Polarization of the Didymos–Dimorphos System before and after the DART Impact

The polarization observations of the Didymos system were done using both VLT and NOT. With VLT, the Focal reducer/low dispersion spectrograph 2 (FORS2) instrument was used in both imaging polarimetric (IPOL) and spectropolarimetric (PMOS) modes, and the Alhambra Faint Object Spectrograph and Camera (ALFOSC) instrument was used in imaging polarimetric mode with the NOT. Observations were done over a time span of 153 days, starting from 2022 August 19 and ending on 2023 January 19. The DART spacecraft impacted Dimorphos on 2022 September 26, so our observations span from 39 days pre-impact to 114 days post-impact.

The observations are published in depth in Bagnulo et al. (2023) and Gray et al. (2024); please consult those for details. All the observational data used in this article are available in Gray et al. (2024). Four filters were used (B, V, R, I) with the VLT observations, and two (B, R) were used with the NOT. The numbers of observations per filter are 40 (B), 14 (V), 43 (R), and 15 (I). Since the B and R filters have the majority of the observations, we decided to use only these in our analysis. Please see a graphical presentation of the observations with different filters on the timeline in Figure 1.

Zoom In Zoom Out Reset image size

Regarding phase angle, the angle between the Sun and the telescope as seen from the target, the observations span between 655 and 7634. This means that the observations cover the so-called negative polarization branch, the inversion angle where the phase curve turns from negative to positive, and the start of the positive polarization branch for typical asteroid polarization phase curves.

We use the so-called trigonometric function to model the continuous polarization phase curve, i.e., the DoLP (P) as a function of the phase angle (Lumme & Muinonen 1993; Penttilä et al. 2005). The model is of the form

$$\begin{eqnarray}&&{P}_{m}(\alpha ;b,{c}_{1},{c}_{2},{\alpha }_{0})=b\,{\sin }^{{c}_{1}}\alpha \,{\cos }^{{c}_{2}}\tfrac{1}{2}\alpha \,\sin \left(\alpha -{\alpha }_{0}\right),\end{eqnarray} \tag{ 1 }$$

where b, c₁, c₂, and α₀ are parameters to be fitted. The parameter α₀ is the inversion angle, b controls the amplitude of the curve, and c₁ and c₂ control the shape. One benefit of this model over other empirical polarization phase curve models is that the trigonometric function is designed to meet the physical constraints of zero DoLP at both α = 0° and 180°. Therefore, it is a more reliable model when extrapolating values outside the observed phase angles.

For statistical analysis between polarization phase curve models with different numbers of free parameters, we use the Bayesian information criterion (BIC). BIC can be defined with

$$\begin{eqnarray}&&\mathrm{BIC}=n\mathrm{log}\left(\displaystyle \frac{\mathrm{SSE}}{n}\right)+\mathrm{log}\left(n\right)p,\end{eqnarray} \tag{ 2 }$$

where n is the number of observations, p the number of fitted parameters, and SSE the sum of squared errors of the model. Smaller values of BIC indicate better models. Small SSE will result in small BIC, but large values of p, i.e., many free parameters in the model, penalize by increasing the BIC value. The trigonometric model itself is fitted to data using nonlinear least-squares minimization.

We will derive physical properties of the regolith and the ejecta particles in the Didymos system using geometric optics and radiative transfer light-scattering simulations (see the next section). These simulations need a shape model for the particles in the regolith and the ejecta, while the other parameters (size and the complex refractive index of the material) are varied and then inferred by comparing the simulation results and the observations. We choose to model the particles with so-called random Voronoi polyhedra. The modeled particles are formed through an algorithmic procedure. First, a spherical volume is filled with a number of seed points, which are randomly distributed into the volume. Second, a Voronoi division is computed around the seed points. A Voronoi cell or particle for a specific seed point is defined by the volume where any point is closer to that specific seed point than any other seed point. The definition leads to polyhedral particles. The particles close to the surface of the original spherical volume might have no closing faces, so only particles from the middle parts of the spherical volume are included in the final sample.

All the sampled particles are normalized to have a volume of $\tfrac{4}{3}\pi$ so that their volume-equivalent sphere radius is 1 and can later be scaled to any value. An example of Voronoi particle shapes is shown in Figure 2.

Zoom In Zoom Out Reset image size

We use light-scattering simulations to derive the DoLP of the modeled Didymos system to be compared to observations. Reflected light in visual wavelengths for Didymos comes either from the particles in the very surface of the Didymos and Dimorphos (pre-impact) or partly also from the particles in the ejecta cloud around the system (post-impact). Therefore, we model the system as an optically thick cloud of particles. In this work, the geometric shape of the ejecta cloud and the binary geometry of the bodies are not taken into account, but a spherical and optically thick cloud of particles is used. Further, the particles are assumed to be optically large, that is, have sizes much larger than the wavelength of light. In practice, this means particles with a diameter of tens of micrometers and above. The reason for using large particles is given in Section 3.2.1.

The simulations are done using two codes, SIRIS4 (Muinonen et al. 2009; Lindqvist et al. 2018) and radiative transfer and coherent backscattering (RT-CB; Muinonen 2004) codes; see, e.g., Penttilä et al. (2020) for previous work with a similar work flow. The single-scattering properties, i.e., the single-scattering albedo and the scattering matrix, of an individual grain are simulated using SIRIS4. The input parameters are the shape of the particle, where a collection of 30 random Voronoi shapes are used; size in terms of size parameter x; and complex refractive index m. The size parameter is defined as $x=\tfrac{2\pi r}{\lambda }$, where r is the particle radius (volume-equivalent sphere radius) and λ is the wavelength. Complex m is m = n + ik, where n is the real refractive index and k is the imaginary part.

The output of the single-scattering simulations is averaged over the 30 shape samples and fed to RT-CB simulation of a practically infinitely large and optically thick spherical volume of particles. Since the input scattering matrices come from a geometric optics simulation, they do not include phase information of the wave, so the CB part of the code is not used. The output of the RT simulation is the scattering matrix M of the system as a function of phase angle α. From the scattering matrix, the DoLP can be computed as $P(\alpha )=-\tfrac{{\left[{\boldsymbol{M}}\right]}_{21}}{{\left[{\boldsymbol{M}}\right]}_{11}}$.

We run the simulations in a regular grid with input parameters x, n, and k. The real part n has values 1.25, 1.5, 1.75, and 2, but k and x have values equally spaced in log₁₀-space: ${\mathrm{log}}_{10}(x)=1,1.2,1.4,\ldots ,3.6$, and ${\mathrm{log}}_{10}(k)=-6,-5.5,-5,\ldots ,0$. For B and R filters, the size parameters map to particle diameters from about 1.4 to 700 μm. The simulation results include the whole scattering matrix of the system as a function of phase angle, but we only need certain elements to compare to our polarization observations. As a diagnostic feature in polarization, we select the maximum value ${P}_{\max }$ of the polarization phase curve P(α). We can extract the maximum value from both the simulations and the observations fitted with the trigonometric polarization phase curve model.

The simulations show interesting trends that allow us to simplify many result features. Within this RT framework, features such as the integrated scattered intensity and the polarization maximum depend on x and k together as ${\mathrm{log}}_{10}(x)+{\mathrm{log}}_{10}(k)$, that is, the size and the imaginary part of the refractive index have similar effects, especially when combined in log-scale (Muinonen et al. 1996). The real part n, however, has only a minor effect. We fixed n = 1.6, a typical value for olivine and asteroid material having abundant olivine, and were able to successfully fit a smooth regression spline to ${P}_{\max }\left({\mathrm{log}}_{10}\left(x\right)+{\mathrm{log}}_{10}\left(k\right)\right)$ data.

The polarization observations of the Didymos system in the B and R filters are shown in Figure 3, together with the best-fitted models. With the continuous phase curve model, we can inspect the observations both as a function of phase angle and as a function of time and compare between different subsets, such as pre- and post-impact observations. As noted already in Bagnulo et al. (2023) and Gray et al. (2024), there is a drop in DoLP at the time of the DART impact, and the effect lasts at least some tens of days after the impact. In what follows, we verify this statistically and model the ${P}_{\max }$ data for pre- and post-impact cases and for B and R filters.

Zoom In Zoom Out Reset image size

The differences in the DoLP phase curves are verified by comparing BIC values for different models. By fitting one model P_m (α) to all data, we can first conclude that parameter c₂ in the model can be almost 0, meaning that the part ${\cos }^{{c}_{2}}\left(\tfrac{\alpha }{2}\right)$ can be approximately 1. This is typical for the trigonometric model we have used, since this part describes behavior close to maximum, and often there is a lack of data supporting a confident estimation of c₂. The fitted parameter values and the BIC value for this model and all the other models mentioned in this section are listed in Appendix A.1. This model is referred as "m0."

Second, we compare models where the filters B and R are fitted either separately (m1) or jointly having shared c₁ and α₀ parameters (m2). The BIC analysis shows that the best model is the one where both filters have different amplitudes (parameter b) but shared shape of the curve (c₁ and α₀ parameters).

Third, we check that the three different observing instruments (NOT/ALFOSC and VLT in IPOL and PMOS modes) are not producing statistically different data. We compare BIC values between one model (m0), three models with different amplitudes but shared shape parameters (m3), and three completely different models (m4). The results indicate that the one joint model between all the instruments performs best, so there is no difference between the instruments.

Lastly, we get to the interesting part when comparing models that separate pre- and post-impact data (m5–m8). Here the result is that the model where pre- and post-impact data have fitted b and α₀ parameters separately (m7) outperforms all the other combinations, including the one where they are fitted with a single model. Therefore, we can conclude that the pre- and post-impact polarization is (statistically) different, and that difference shows in the P_m model in different magnitude (b) and inversion angle (α₀) parameters for pre- and post-impact models, but the shape parameter c₁ is better shared.

In Figure 3 we show the best models, i.e., model m7 with separate and joined parameters resulting in four different models for (i) filter B and pre-impact, (ii) filter R and pre-impact, (iii) filter B and post-impact, and (iv) filter R and post-impact; see Appendix A.1. The models are shown together with the observations as a function of both phase angle and time. The time-dependent model can be constructed by interpolating the phase angle versus time (α vs. t) dependence of the Didymos system, as seen from Earth, and then plotting P_m (α(t)).

The polarization properties of the final fitted models, that is, the angles of minimum and maximum polarization ${\alpha }_{\min }$ and ${\alpha }_{\max }$, the inversion angle α₀, the values of minimum and maximum polarization ${P}_{\min }$ and ${P}_{\max }$, and the slope of polarization at inversion h, are shown in Tables 1 and 2. These values compress our observational data. We find especially the ${P}_{\max }$ values interesting; large phase angles are present in both pre- and post-impact observations, the model selection suggests differences between both colors (${{\rm{\Delta }}}_{c}={P}_{\max ,\mathrm{filter}\ R}-{P}_{\max ,\mathrm{filter}\ B}$) and times (${{\rm{\Delta }}}_{t}={P}_{\max ,\mathrm{post}-\mathrm{impact}}-{P}_{\max ,\mathrm{pre}-\mathrm{impact}}$), the differences are the largest in large phase angles, and the RT simulation without the CB component is accurate with large phase angles, but not so much with small phase angles, where CB would have a significant effect. For all these reasons, we will base our analysis on the Didymos system particles on these four differences shown in Table 2.

Table 1. Properties of the Final Phase Curve Models for the DoLP Observations, Apart from Maximum Polarization, Which Is Shown in Table 2

	α0 (deg)	${\alpha }_{\min }$ (deg)	${\alpha }_{\max }$ (deg)	${P}_{\min }$ (%)		Slope (%/deg)
				B	R	B	R
Pre	19.68	8.37	101.31	−0.55	−0.44	0.090	0.073
Post	22.56	9.60	102.96	−0.65	−0.55	0.094	0.078

Download table as:  ASCIITypeset image

Table 2. The Maximum DoLP Values in the Models Fitted to Observations, and Their Differences between Colors (Filters B and R; Δ_c ) and Time (Pre-impact and Post-impact; Δ_t )

	${P}_{\max }$ (%)
	B	R	Δc
Pre	11.21	9.03	−2.18
Post	10.50	8.77	−1.73
Δt	−0.71	−0.26

Download table as:  ASCIITypeset image

3.2.1. Polarization Mechanism with Optically Large Particles

The polarization properties of the Didymos system show signs of optically large particles dominating the polarization behavior. There are two indications of this. First, the angle of polarization maximum is shifted from the 90° value, typical for Rayleigh scattering from optically small particles, toward larger values above 100° (Muñoz et al. 2021; Escobar-Cerezo et al. 2018; Frattin et al. 2022). Second, the observed behavior of the DoLP maximum value decreasing with increasing wavelength (filters B and R) supports large particles (Dabrowska et al. 2015). In the geometric optics domain, the DoLP is decreased when there are more reflections from the vacuum−particle boundaries, since this is where the polarization state is changed and therefore randomized toward zero DoLP. The absorption of the light happens in the volume of the particles. When the same enclosing space is hosting a material volume that is divided into larger or smaller particles, the one with smaller particles will have more boundaries between material and vacuum, therefore introducing more randomized DoLP (Sorensen 2001).

The mechanism of decreasing DoLP can be resulting from either decreasing size parameter of particles or decreasing imaginary part k of the material refractive index, since the latter also increases the boundary reflections by decreasing absorption in the material. The size parameter can decrease if either the particle size is decreased or the wavelength is increased. The same mechanism that decreases polarization will simultaneously increase brightness for the same reason, increased surface reflections versus volume absorption. In our simulations (Section 2.4) we can verify this in all cases. Other verification comes from, e.g., the Amsterdam-Granada Light Scattering Database ⁷ measurements of olivine (Muñoz et al. 2000) and the Catalogue of Asteroid Polarization Curves; ⁸ see Figure 4. In Granada measurements, the "L" and "XL" size fractions of olivine are sieved to contain only large particles, 20–65 μm for "L" and 65–125 μm for "XL" in diameter, and maximum polarization is decreased from 7.6% to 7.5% (at 442 nm) and from 9.1% to 8.8% (at 633 nm) when moving from "XL" to smaller "L." In the asteroid polarization catalog it is clearly seen how brighter objects such as S-type asteroids, with material absorption k being smaller than with C-type asteroids, have much more neutral polarization. From Figure 4, one can further observer how polarization observations confirm the Didymos system being an S-type object, as concluded also by de León et al. (2010), Bagnulo et al. (2023), Gray et al. (2024), and Lin et al. (2023).

Zoom In Zoom Out Reset image size

3.2.2. Derived Properties of Particles in the Didymos System

We will derive the properties of the particles dominating the polarization signal of the Didymos system by looking at the ${P}_{\max }$ values (8.77%–11.21%) and their differences between the colors (−1.73 and −2.18 percentage points (p.p.)). We can inspect these values with respect to our simulation results. In Figure 5, the simulated ${P}_{\max }$ and Δ_c values are visualized with density plots having x and k as the axis, in log-scale. If we take the observed ${P}_{\max }$ values roughly between 8% and 12%, we can find an area of suitable (x,k) value pairs. As seen in Figure 6, the area is a wide line from the upper left, that is, ${\mathrm{log}}_{10}(k)$ around −2 and ${\mathrm{log}}_{10}(x)$ at 1, to lower right, ${\mathrm{log}}_{10}(k)$ around −4.4 and ${\mathrm{log}}_{10}(x)$ at 3.5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can add a constraint from the observed difference Δ_c in DoLP, although it is more complex to interpret than ${P}_{\max }$. The simulated Δ_c values are computed from simulated ${P}_{\max }$ values between different size parameters owing to different wavelengths but assuming constant k. This assumption is not completely valid—the imaginary part of the refractive index for olivine or ordinary chondrite meteor material linked with S-type asteroids is typically lower on R-filter than on B-filter wavelengths (see, e.g., Lucey 1998; Muñoz et al. 2000; Database of Optical Constants for Cosmic Dust ⁹ ). This means that the material is brighter with the R filter and the DoLP is lower. For this reason, our observed Δ_c of −1.73 and −2.18 is at least partly due to different k values between the filters and can offer only a lower negative limit to the change owing to the size parameter changing with wavelength. This will exclude the area of fast ${P}_{\max }$ change above the stripe proposed by ${P}_{\max }$ values in Figure 6. In conclusion, based on our simulations, we believe that the Didymos system particle size and absorption properties should reside in the area constrained by suitable ${P}_{\max }$ values.

The multiple scattering mechanism in optically large particles and the interplay between surface reflections and refractions (neutralizing polarization, not dependent on k) and volume absorption (decreasing surface interactions, dependent on k) lead to the relation where there is no single (x,k) pair that can produce the observations, but a collection of (x,k) value pairs. Larger values of k will counteract the smaller sizes; smaller k, the larger sizes.

Olivine is a very abundant mineral in asteroids and also the major component affecting the scattering properties of S-type asteroids. This is seen in observed S-type spectra that closely follow laboratory-measured spectra of pure olivine. It is therefore natural to assume that Didymos particles mainly consist of olivine and that the optical constant of the particles should match that seen in olivine and in ordinary chondrite meteorites. This would indicate values of ${\mathrm{log}}_{10}(k)$ around −4, typically between −6 and −3 (see, e.g., Muñoz et al. 2000; Database of Optical Constants for Cosmic Dust; see footnote 9). Turning this into particle size according to constraints from our simulations, the dominating particle size would be around ${\mathrm{log}}_{10}(x)=3.1$, ranging from 2.1 up to about 3.5. In visual wavelengths at the middle of B and R filters, ${\mathrm{log}}_{10}(x)=3.1$ corresponds to particle diameters of ∼200 μm.

3.2.3. Particles before and after DART Impact in the Didymos System

After finding a plausible particle size (about 200 μm) and imaginary part of refractive index (about 0.0001), we can turn to the observed change in DoLP of the Didymos system with the DART impact. We fitted the pre- and post-impact DoLP observations in B and R filters and observed a change of −0.71 p.p. (B filter) and −0.26 p.p. (R filter) (Table 2). With a fixed wavelength at the given filter, a decrease in ${P}_{\max }$ can be due to either the size of particles decreasing or the imaginary part of the material k decreasing. We can fix the initial particle properties to our plausible pair of particle diameter d = 200 μm and k = 0.0001 and see what kind of change in either d or k at the DART impact would explain the observed Δ_t .

Since both x and k have an effect, it is best to fix one and see how the other one influences ${P}_{\max }$. Let us first fix k = 0.0001 and study what kind of change in d would cause the observed Δ_t. This can be seen in Figure 7; Δ_t = −0.71 in the B filter can be explained by particle diameter decreasing from about 230 to 210 μm, and the same change of about 215 to 205 μm in the R filter. Hence, a 10–20 μm decrease in the dominant particle size between pre- and post-impact observations can explain the observed DoPL.

Zoom In Zoom Out Reset image size

On the other hand, we can fix d and see what kind of change in k is needed. Now d is fixed to 200 μm. From the right panel of Figure 7 we see that a decrease between 0.02 and 0.04 in ${\mathrm{log}}_{10}(k)$ can explain the observed Δ_t with both the B filter (${\mathrm{log}}_{10}(k):-3.94\to -3.97$) and the R filter (${\mathrm{log}}_{10}(k):-3.97\to -3.99$). In k this change is about −0.00001.

Table A1 summarizes the different phase function models fitted to the DoLP observations.