© ESO, 2017

1. Introduction

Remote sensing of comets with ground-based radar can be accomplished if their distance to Earth is less than <0.5 AU (Harmon et al. 2004). By transmitting a monochromatic electromagnetic wave and measuring the backscattered signal, it is possible to obtain information about the nucleus, its size, rotation, surface roughness, and its environment. At smaller distances, <0.1 AU, Doppler imaging of nuclei can be conducted (Harmon et al. 2010). Typically, a circularly polarized pulse is transmitted and both opposite-circular (OC) and same-circular (SC) polarization echoes are measured. The OC Doppler spectra of comets with a significant population of large grains show narrowband peaks from the nucleus, and a broadband signal is expected from moving particles in the coma. Using this technique, Doppler spectra measured for some of the comets show the presence of the so-called large-grain coma around cometary nuclei (Harmon et al. 2011, 2006, 2010; Nolan et al. 2006). Same-circular echoes are usually much weaker. The circular polarization ratios SC/OC calculated for these spectra indicate the presence of scatterers that have radii larger than λ/ 2π, which is the Rayleigh-Mie criterion. Model-based estimations of the dust ejection from nuclei can constrain the material of these cm-sized scatterers, which are most likely icy grains (Kelley et al. 2013).

Modeling the circular polarization ratio also allows an estimation of the particle size distribution (Harmon et al. 2011). A more reliable study is possible with rigorous simulations of electromagnetic scattering by such particles. The Mie theory can be applied (e.g., Epifani et al. 2001) for approximate modeling, assuming spherical grains. However, this approach would not be as rigorous, because of uncertainties, as electromagnetic scattering from spherical and irregular particles can be very different, in particular at exact backscattering. The other option is to use an aggregate model of clusters of spheres (e.g., Kelley et al. 2013). This seems to be reasonable for cometary dust particles that have sizes on the order of 1−10 μm. We note, however, that computations where clusters of spheres are used do not provide satisfactory fits for the photometric and polarimetric measurements of light scattered by cometary dust (Lasue et al. 2009; Das et al. 2011). In the case of the cm-sized chunks, it is also doubtful that such large objects can have a structure of ballistic clusters with μm-sized constituent grains.

Another approach is to model cometary dust particles with irregularly shaped agglomerates, which resemble the morphology of dust sampled in the vicinity of comets. One of these models uses agglomerated debris particles and is capable of reproducing astronomical observations of comets in the visible (e.g., Zubko et al. 2015, 2016). However, this particle morphology was designed to reproduce features of cometary dust particles on a scale of up to a few microns; whereas the applicability of this shape model on the decimeter scale is not justified. In the latter case, it is more reasonable to expect a resemblance of shape with what appears in a whole cometary nucleus.

Therefore, we need model shapes that are random, irregular and compact at the same time. The so-called Gaussian Random Field particles (Grynko & Shkuratov 2002, 2003, 2007; Shkuratov & Grynko 2005) seem to be a good first approximation to decimeter-sized pieces of a comet.

In this paper, we analyze the radar measurements of the circular polarization ratio in cometary comae, using our computational results of electromagnetic backscattering from ensembles of irregular particles of different size. For numerical simulations of electromagnetic scattering we use the discrete dipole approximation (DDA) method, which allows arbitrary shapes of scatterers. Section 2 describes the details of the radar observations of comets. Furthermore, we present our numerical model in Sect. 3 and discuss the results of simulations in Sect. 4.

2. Radar observations of comets

Table 1 shows a few of the comet radar observation epochs that we consider in our work. The data comes from the Arecibo and Goldstone observatories, and the listed comets have shown a detectable response from the grain coma. In the table, A_s and G_x indicate S-band (λ = 12.6 cm) for Arecibo and X-band (λ = 3.5 cm) for Goldstone, respectively.

A typical observation was done for comet Hartley in October 2010. Radar observations were made at Arecibo (Harmon et al. 2011). A 750 kW, circularly polarized wave was transmitted and an echo in both OC and SC polarizations was received. Comet Hartley shows a two-component echo. These two components are distinguishable and the coma does not seem to cloud the nucleus. This shows the coma to be optically thin (Harmon et al. 2004), i.e., its optical thickness τ is much smaller than 1. In other words, it implies no multiple scattering between large grains in the coma. Overall, this is quite consistent with low spatial concentration of large chunks detected in situ in Comet 103P/Hartley 2 (Kelley et al. 2013). The observed circular polarization ratio, μ_c, is 0.34 which is an intermediate value in the μ_c range (0.10−0.59) of the observed cometary nuclei (Harmon et al. 2011). Other radar observed comets shown in Table 1 also show a response from the grain coma. Their μ_c values and wavelengths of observation are listed in Table 2.

3. Numerical model

The DDA is a versatile numerical method used to compute the absorption and scattering of electromagnetic radiation by particles of arbitrary shapes and composition (Purcell & Pennypacker 1973). In this method the scatterer is divided into a discrete set of dipoles on a cubic lattice. The dipoles interact with each other and the incident field, which form a set of linear equations that can then be solved to obtain dipole polarizations. All the other scattering quantities can be computed from these polarizations. To be realistic in the simulations, the shape of the scatterers should be close to those found in nature. The shape of cometary dust particles is irregular at different scales with respect to the wavelength. Relatively large particles may have sharp edges and surface roughness smaller than or comparable to the wavelength of incident radiation. This plays an important role in the formation of the linear and circular polarization response. Flexibility of the DDA method to simulate particles of irregular shapes effectively makes it a good choice for our study. We use the Amsterdam DDA (ADDA) parallel code (Yurkin et al. 2007).

The model particles in our simulations are Gaussian random field particles (GRF; Grynko & Shkuratov 2002, 2003, 2007; Grynko et al. 2013). The generation technique includes the creation of a 3D random field with Gaussian statistics in computer memory. A 3D field can be considered as a 4D single-valued random topography with the Gaussian statistics of heights and that of a Gaussian correlation function. This topography is then dissected by a 3D hyperplane at a certain level parallel to the average plane of the topography. With this procedure an ensemble of 3D random shapes of different sizes is formed. Alternatively, each value of the random field can be compared to a certain threshold constant. The point is attributed then to the material if it is larger than the constant and to the empty space in the opposite case. In the end, we convert this set of grid points representing the particle volume to the faceted description of the particle surface (Fig. 1).

Fig. 1 Examples of randomly shaped GRF particles used in the simulations.

Generally, the shapes of the generated particles are random and irregular, but they can also be ellipsoid-like or have very complicated porous geometries (Grynko & Shkuratov 2008). Tuning the free parameters, we chose realistic compact shapes like those shown in Fig. 1.

In theory and modeling of electromagnetic scattering, instead of the absolute size of the scattering particle, the size relative to the incident wavelength, known as the size parameter, is used (1)where r is the radius of the particle and λ is the wavelength. In the case of non-spherical scatterers, the assessment of size parameters is not obvious and a volume equivalent size parameter is often used. Here, we define r as the radius of the sphere circumscribing the particle. As suggested by the radar observations the grain coma includes a significant population at decimeter sizes or even larger (Harmon et al. 2011). Our study includes ensembles of GRF particles, 100 samples for each size, ranging from 6.5 cm to 100.8 cm (diameter of the circumscribing sphere) with a step of 12.6 cm. This means a size parameter varying from 1.5 to 25. Larger sizes become a difficult problem due to quickly growing computational costs.

We also consider two types of compositions expressed through the complex refractive index m: pure water ice with m = 1.78 + 0.003i at λ = 12.6 cm (Warren 1984) and dirty ice that is contaminated with 30% silicate inclusions. We assume olivine inclusions with refractive index m = 2.49 + 0.0124i (Campbell & Ulrichs 1969; Virkki et al. 2014) as an inhomogeneity diffused randomly in the scatterer. The choice of the composition stems from the prior knowledge of the constituents of interplanetary dust.

In DDA, the volume of the scatterer is divided into a set of dipoles. For irregular particles the cells must satisfy the condition kd | m | < 1, where k is the wave number, d is the size of the dipoles, and m is the complex refractive index of the material (Zubko et al. 2010). For GRF particles our model showed convergence in intensity and linear polarization for 20 dipoles per wavelength, which corresponds to approximately kd | m | = 0.5.

The scattering response of the particles is calculated for all scattering angles averaged over 180 azimuthal planes. In the output, the ADDA code gives the angular dependence of the 16 elements S_ij of the square scattering matrix (Bohren & Huffman 2003). For our study we are interested in the first element S₁₁, which is the intensity, and the element S₄₄, which characterizes the circular depolarization property of the scatterer. So, the radar measurable circular polarization ratio, μ_c, in terms of scattering matrix is defined as (Bohren & Huffman 2003) (2)where θ is the scattering angle. At backscattering the formula becomes (3)Furthermore, it is relevant to know the size distribution of the dust particles in the grain coma. The distribution follows a power law r^{− n} (Mazets et al. 1986; Green et al. 2004), meaning the probability of occurrence is inversely proportional to the size of the particle raised to an appropriate index. As mentioned above, we average the S₁₁ and S₄₄ values for 100 particle samples of each size. These average values, and of each individual size, weighted by a power-law distribution function are summed. This is done for a range of power indices ranging from 2 to 6.5 in steps of 0.5. For a power index n and size range a_i ranging from x_min to x_max the final values of S₁₁ and S₄₄ are obtained as With these values we then calculate the polarization ratios μ_c.

Fig. 2 Intensity S11 (a) and ratio S44/S11 (b) as functions of the scattering angle computed for ensembles of spherical and irregular “icy” particles with the size distribution power index n = 3.

To analyze the relative contribution of particles with different size parameters to scattering we plot in Fig. 3a their average scattering cross sections σ_s weighted by the distribution function with power index n = 3 and perform a normalization according to (6)For particles with X < 4 the step is 0.5 in size parameter units. Those particles that are smaller than the wavelength appear to be the most efficient scatterers for the considered refractive indices. However, this regards only the intensity of the scattered light. The μ_c values that these particles produce at backscattering are very low, nearly zero (see Fig. 3b), and therefore a wide range of sizes should be taken into account.

It is also necessary to make sure that multiple scattering can be neglected. Harmon et al. (2004) estimated the optical depth τ of the grain coma for a line of sight passing through its central point, implying the decrease in the grain number density with the distance from the nucleus as 1 /r². The obtained upper limit of τ ~ 10^-4 indicates that single scattering dominates, whereas multiple scattering appears important only when τ is greater than 1 (Chandrasekhar 1960). Such a low value is also in agreement with the images of the nucleus environment of comet 103P/Hartley 2 obtained in situ by the Deep Impact probe (Kelley et al. 2013). An additional argument for a thin coma is the ratio of the absorption and backscatter cross sections σ_a/σ, estimated in Harmon et al. (2004) as low as ~ 0.1.

As the circular polarization ratio is sensitive to the shape of scatterers, the corresponding scattering matrix elements computed for ideal spheres and irregular particles appear to be quite different (Virkki & Muinonen 2016). To illustrate this point in our specific case we compared the scattering angle dependencies of intesity S₁₁ and ratio S₄₄/S₁₁ computed for two sets of “icy” particles, spherical and GRF, with sizes distributed according to power law (Fig. 2). A difference is clearly seen in both plots for angles >90°. In particular, this regards the ratio S₄₄/S₁₁ at exact backscattering.

Fig. 3 a) Normalized weights (average scattering cross sections σs weighted by the size distribution function) of particles with different materials as a function of size and with size distribution power index n = 3. b) Unweighted circular depolarization ratio μc as a function of size.

4. Results and discussion

In Fig. 4 we present the circular polarization ratios computed for various size distributions and compare them with those measured for several comets. As the probability distribution of the cometary dust varies as the power law of constituting sizes, variation of power indices governs the relative contribution of particles of different sizes to the scattering. The plot shows data for power indices from 2 to 6.5. The exact values of the radar measured μ_c values can be found in Table 2.

Fig. 4 Simulated μc values vs. power index compared with measurements for seven observed comets.

With these results we can apply our simulations to constrain the properties of the grain coma chunks. We can see that the distribution law strongly influences the μ_c value. With reasonable choices of the power index n in the range from 2 to 6 the value of μ_c can be tuned so that the simulations reproduce the observed values, except for comet IRAS-Araki-Alcock (IAA) as it has very low value of μ_c close to zero. For the rest of the comets we can constrain the particle size distribution laws finding power indices that produce matching μ_c values. It is clear that n can range from ≈3 for comet Halley up to 5 for comet Tuttle.

A different material, dirty water ice contaminated with small olivine grains, results in smaller circular depolarization ratios. However, the difference mostly does not exceed the accuracy limits of the radar measurements. Only at small n, i.e., if large chunks have larger number densities, does it become significant. At n > 3 the size distribution law appears to be a stronger factor that influences the μ_c value. For comets with μ_c < 0.2 the difference in the radar response cannot only be explained by the difference in chemical composition of the chunks.

In situ observations of the large chunks near the nucleus of comet 103P/Hartley 2 done by the EPOXY mission give us another interesting opportunity to apply our results. At close approach of the spacecraft to the nucleus of the comet many individual chunks were seen in the images taken by the spacecraft cameras (Kelley et al. 2013). From the dynamics, approximate lifetime, and brightness of the chunks, their size range has been roughly estimated between 1 and 15 cm for pure water ice and 10 and 150 cm for dirty ice. This is in agreement with our simulations. Comparison of our results with the radar measured circular depolarization ratio of 0.34 gives us a constraint for the power index for the coma of this comet in the range of n≈ 3.3−4.3. Assuming dirty ice results in a similar tendency of increasing particle sizes through decreasing power index of the size distribution law.

Based on the data it can be inferred that the comets showing lower μ_c values did not have many particles of appreciable size to contribute to the observed response or have a different size distribution. We note that the observations for comet Hayakutake (Fig. 4) were performed at λ = 3.5 cm.

Comet IAA (Fig. 4), showed the coma echo only ≈ 1% depolarized (Fig. 4). It is the lowest for any solar system radar echo. This is quite consistent with the real cutoff for maximum grain radii not being much larger than λ/ 2π (Harmon et al. 2004).

On the other hand, comets Halley and Hayakutake showed a non-negligible μ_c value, but the significant echo from comet C/2001 A2 (Fig. 4) gave information about the proportionality of the depolarization and particle size (Nolan et al. 2006). The nucleus of C/2001 A2 had split before radar observations (Sekanina et al. 2002) and its quite probable that the strong depolarization was the response from these big boulders that formed a part of the coma. The same can be said for comet 103P/Hartley 2 and 8P/Tuttle (Fig. 4). A high response for depolarization and hence the μ_C value may represent the presence of boulders that are larger than the wavelength. Based on our model results obtained with sizes ranging up to 1 m, the power index of the size distribution for comet 103P/Hartley 2 can be estimated in the range 3.5−4 and even higher for comet 8P/Tuttle.

The effect of composition is also interesting. As can be seen in Fig. 4, the effect of impurity tends to lower the depolarization. As an example, for comets 103P/Hartley, C/2001 A2, and Hayakutake the pure ice computations result in a size distribution with power index around 4. Considering the presence of absorbing impurities the power index shifts towards lower values due to the decreased μ_c. As the size distribution assumes inverse power law dependence, for a μ_c value that fits in the range of the larger power index, the response from smaller particles remains more dominant. In this region μ_c is very small, which is a property of particles much smaller than the wavelength, and the difference in chemical composition does not play a role.

5. Conclusions

The principal results of our numerical simulation of radar response from comets can be summarized as follows:

• 1.

  Large chunks detected by the Deep Impact space probe in thevicinity of the nucleus ofcomet 103P/Hartley 2 may have apredominantly water ice composition, as was suggested in(A’Hearn et al. 2011). Under theassumption of a power-law size distribution of large chunks, thepower index n in comet 103P/Hartley 2 can be constrained to the range n ≈ 3.3–4.3.

• 2.

  Except for comet IRAS-Araki-Alcock, we have evidence for large chunks in the vicinity of all considered comets that return detectable radar echo from the coma. However, the observed difference in the circular polarization ratio μ_c between comets cannot only be explained with the difference in chemical composition of the chunks. On the contrary, varying solely the size distribution can reproduce the whole range of μ_c measured in comets. At that, the power index n may range from n ≈ 2 up to 5. We also note that the population of the coma and the size distribution function can vary with the distance from the Sun.

• 3.

  The chemical composition of large chunks, pure water ice or water ice contaminated with olivine, noticeably affects the circular polarization ratio μ_c only if it takes on relatively large values, μ_c > 0.2. Below that threshold μ_c is mostly governed by the size distribution.

Acknowledgments

The authors gratefully acknowledge the computing time granted by the Paderborn Center for Parallel Computing (PC2). The authors are also thankful to the reviewer Anne Virkki for her thoughtful comments on the manuscript.

References

All Tables

All Figures

	Fig. 1 Examples of randomly shaped GRF particles used in the simulations.
In the text

	Fig. 2 Intensity S11 (a) and ratio S44/S11 (b) as functions of the scattering angle computed for ensembles of spherical and irregular “icy” particles with the size distribution power index n = 3.
In the text

	Fig. 3 a) Normalized weights (average scattering cross sections σs weighted by the size distribution function) of particles with different materials as a function of size and with size distribution power index n = 3. b) Unweighted circular depolarization ratio μc as a function of size.
In the text

	Fig. 4 Simulated μc values vs. power index compared with measurements for seven observed comets.
In the text