Implications of Philae Magnetometry Measurements at Comet 67P/Churyumov–Gerasimenko for the Nebular Field of the Outer Solar System

We assumed the interior of the comet to be filled with uniformly magnetized spheres of equal size, equal density, and equal magnetization with random orientations. The magnetic field outside a uniformly magnetized sphere is identical to that of a dipole with the same moment at the center of the sphere, so our model of 67P is equivalent to a collection of approximately uniformly spaced dipoles. Because any arbitrary magnetization pattern can be represented by a collection of properly positioned and oriented dipoles, our model geometry provides a robust way to estimate the intensity of the remanent magnetic field as a function of the length scale of coherent magnetization in the comet. In this model, using spheres with a radius (r_d) comparable to the radius of individual ferromagnetic grains (r_gr) in the comet produces an estimate for the remanent field of a comet with no signature of ADRM (i.e., the magnetic moments of individual grains have random mutual orientations). This is a lower bound for the NRM of the comet. Values of r_d > r_gr correspond to some degree of mutual alignment of the magnetic moments of individual ferromagnetic grains and produce higher estimates for the remanent magnetization of 67P.

Our simulation performs sphere packing using a hexagonal close-packed lattice to populate the bounding box of the comet shape and a modified ray-casting algorithm to discard spheres that have centers outside the comet (Figure 1; see Appendix A for details). This achieves a filling factor (fraction of volume filled by spheres) of $\phi =\pi /(3\sqrt{2})\approx 0.74048$. For simplicity, we did not require that the entirety of each sphere be inside the shape model. Since the size scales at which this inaccuracy is relevant are usually below the spatial resolution of the shape model, this does not significantly affect the results. The magnetic field is then calculated at a given location ${\boldsymbol{r}}$ outside 67P by summing the contributions from each dipole:

$$\begin{eqnarray}&&{\boldsymbol{B}}({\boldsymbol{r}})=\displaystyle \frac{{\mu }_{0}m}{4\pi }\displaystyle \sum _{i}\left(\displaystyle \frac{3{\hat{{\boldsymbol{d}}}}_{i}({\hat{{\boldsymbol{m}}}}_{i}\,\cdot \,{\boldsymbol{p}}{\hat{{\boldsymbol{d}}}}_{i})-{\hat{{\boldsymbol{m}}}}_{i}}{{d}_{i}^{3}}\right),\end{eqnarray} \tag{ 1 }$$

where μ₀ is the vacuum permeability, ${{\boldsymbol{r}}}_{{\boldsymbol{i}}}$ is the location of the ith dipole, and ${{\boldsymbol{m}}}_{{\boldsymbol{i}}}$ is the dipole's vector moment, which has magnitude m. The vector from the ith dipole to ${\boldsymbol{r}}$ is given by ${{\boldsymbol{d}}}_{{\boldsymbol{i}}}={\boldsymbol{r}}-{{\boldsymbol{r}}}_{{\boldsymbol{i}}}$. In the case that d_i < r_d, the contribution from that dipole is set equal to the field intensity at the surface of the uniformly magnetized sphere, d_i = r_d. This can occur as a result of requiring only sphere centers to be included in the shape model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The number of spherical magnetic domains is roughly proportional to r_d⁻³, where r_d is the sphere radius, approximately the scale of coherent magnetization. Because of this scaling, sphere packing and field calculations become computationally expensive at small magnetization scales. To address this problem, we incorporated an approximation based on a maximum interaction distance, D_max. Any point on the comet at a distance d_i > D_max for all sampled points in a trajectory was ignored. D_max was calculated as described below so that the combined magnetic field from the excluded region of the comet did not exceed a specified threshold B_min.

We seek to estimate D_max such that it it reduces computation time without introducing inaccuracies in the calculated field. We approach this by estimating an upper limit on the field of the ensemble of excluded dipoles. In particular, we estimate the field of the excluded dipoles as that of a single dipole with magnetic moment given by the vector sum of the magnetic moments of the excluded dipoles. From Heslop (2007), the probability density function for the magnitude, M_N, of the vector resulting from the summation of N randomly oriented vectors of magnitude m is

$$\begin{eqnarray}&&f({M}_{N})=\displaystyle \frac{3\sqrt{6}}{\sqrt{\pi }}\displaystyle \frac{{M}_{N}^{2}}{{m}^{3}{N}^{3/2}}{e}^{-3{M}_{N}^{2}/2{m}^{2}N}.\end{eqnarray} \tag{ 2 }$$

Integrating this to find the expected magnitude yields

$$\begin{eqnarray}&&\langle {M}_{N}\rangle =\sqrt{\displaystyle \frac{8}{3\pi }}m\sqrt{N}\approx m\sqrt{N},\end{eqnarray} \tag{ 3 }$$

the familiar scaling relationship for a random walk with step-size m. Taking the magnetization scale, r_d, as the sphere radius and assuming a specific magnetic moment, j, each dipole has a magnetic moment

$$\begin{eqnarray}&&m=\displaystyle \frac{4\pi {r}_{d}^{3}}{3}\rho j,\end{eqnarray} \tag{ 4 }$$

where ρ is the density of cometary material. Note that N is a priori unknown for the volume of the body with r > D_max. To save computation time, we estimate N without actually packing the body. We approximated N as the volume of the comet excluded (V_excl) divided by the dipole volume:

$$\begin{eqnarray}&&N\approx \displaystyle \frac{3{V}_{\mathrm{excl}}\phi }{4\pi {r}_{d}^{3}},\end{eqnarray} \tag{ 5 }$$

where ϕ ≈ 0.74048 is the filling factor for close-packing of equal spheres. To calculate the resulting field from the spheres in the excluded volume, we approximated the field from the excluded volume as that of a single net dipole at one location with the following conservative assumptions: (1) the dipole has an orientation that creates the maximum possible field at the measurement location (the measurement point lies along the axis of the magnetic moment), (2) the net dipole is at a distance D_max, and (3) the excluded volume is equal to the comet volume (V_excl = V_67P). The resulting magnetic field intensity of the excluded volume is found by combining Equations (4) and (5) with (3):

$$\begin{eqnarray}&&{B}_{N}=\displaystyle \frac{{\mu }_{0}j\rho }{2\pi {D}_{\max }^{3}}\sqrt{\displaystyle \frac{32{r}_{d}^{3}\phi }{9}}\sqrt{{V}_{67{\rm{P}}}}.\end{eqnarray} \tag{ 6 }$$

For a comet volume V_67P = 18.7 ± 1.2 km³ (Preusker et al. 2015), this yields the scaling relation

$$\begin{eqnarray}\begin{array}{rcl}{D}_{\max } & = & 147.27\,{\rm{m}}\times {\left(\displaystyle \frac{\rho }{535\mathrm{kg}{{\rm{m}}}^{-3}}\right)}^{1/3}{\left(\displaystyle \frac{{B}_{\min }}{100\mathrm{pT}}\right)}^{-1/3}\\ & & \times \ {\left(\displaystyle \frac{j}{1\times {10}^{-5}{\rm{A}}{{\rm{m}}}^{2}{\mathrm{kg}}^{-1}}\right)}^{1/3}{\left(\displaystyle \frac{{r}_{d}}{1{\rm{m}}}\right)}^{1/2}.\end{array}\end{eqnarray} \tag{ 7 }$$

Given the ROMAP upper limit of 0.9 nT on 67P's remanent field, Equation (7) shows that large portions of the comet can be ignored during many of our simulations. This dramatically improved the execution time of our simulations.

We performed simulations for r_d at 50, 5, and 1 m over the entire rebound trajectory. Even for a magnetization of j = 1 × 10⁻⁶ A m² kg⁻¹ field intensities approaching 0.9 nT were recorded (Figure 3). Magnetic domain sizes of r_d = 50, 5, and 1 m yielded maximum field intensities of 0.60, 0.53, and 0.53 nT, respectively.

Zoom In Zoom Out Reset image size

Figure 3 also illustrates the behavior of the magnetic field in different distance regimes. For a distance r ≫ r_d, the lander is in the far-field and the comet's magnetic field is approximately that of a single dipole, so B ∝ r⁻³. In Figure 3, the simulation with r_d = 1 m (dark blue) occupies this regime for altitudes h ≳ 10 m, indicated by the dashed black line.

Immediately after the bounce, and immediately before the collision, the simulated magnetic field intensities for the various values of r_d converge, indicating that these simulations have entered a near-field regime. Here, the field at the lander is dominated by the nearest dipole, located at a distance r ∼ r_d. In this regime, because the dipole's magnetic moment ${\text{}}m\propto {r}_{d}^{3}$, and the dipole field intensity B_dipole ∝ mr⁻³, the field intensity is independent of magnetization scale.

From these results, it is clear that, for the magnetization range of interest and r_d ≤ 1 m, the lander will not measure significant magnetic fields in the far-field regime (h ≫ r_d). We therefore conducted simulations at magnetization domain sizes smaller than 1 m on the subsections of the trajectory where we expected appreciable field intensity.

To probe smaller magnetization scales, we restricted our simulations to the time immediately after first touchdown and immediately prior to the collision. We simulated magnetic field measurements for magnetization scales of r_d = 0.25, 0.50, and 1 m at times when the lander was at an altitude h ≤ 10 m. For r_d = 10 cm, we limited the simulation to times when h ≤ 1 m.

The maximum field intensities predicted post-touchdown (Figure 4) for j = 10⁻⁶ A m² kg⁻¹ and magnetization scales of 1 m, 50 cm, 25 cm, and 10 cm were 0.57 nT, 0.38 nT, 0.10 nT, and 0.14 nT, respectively. The corresponding values prior to the collision (Figure 5) were 0.88 nT, 0.61 nT, 0.76 nT, and 0.85 nT, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The maximum field intensities from simulations at the same magnetization domain size (r_d) for the component of the trajectory immediately after touchdown, prior to the collision, and over the entire rebound trajectory agree to within an order of magnitude. The largest discrepancies are at the smallest magnetization scales, where the lander is almost always in the far-field regime, so that the field intensity is highly sensitive to minor changes in geometry or altitude. Even for the lowest assumed magnetization (1 × 10⁻⁶ ${\rm{A}}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$), the simulated field intensities approach the 0.9 nT threshold imposed by the Philae observations, suggesting a very low upper bound on 67P's magnetization.

From the results of each of our simulations, we calculated the maximum specific magnetic moment consistent with no detection of the comet's remanent field:

$$\begin{eqnarray*}&&{j}_{\max }=\displaystyle \frac{{B}_{\mathrm{obs}}}{{B}_{\max ,\mathrm{sim}}}{j}_{\mathrm{sim}},\end{eqnarray*}$$

where j_max is the maximum magnetization, j_sim is the input magnetization for the simulation, B_obs = 0.9 nT, and B_max,sim is the maximum value from the simulation. The maximum simulated field value in a simulation occurs at or near the point of lowest altitude, where the limited accuracies of the shape model (∼2 m) or reconstructed trajectory may play a role. However, this is not a significant limitation because the estimated maximum magnetizations are driven by the altitude of closest approach, and not its precise location. Because we know from the reconstruction of the trajectory (Heinisch et al. 2016) that ground contact occurred, the inaccuracy of the shape model is not significant.

The spherical domains in our simulations have a fixed magnetization, implying that they each are composed of the same ferromagnetic materials that are intimately mixed at spatial scales much smaller than the sphere radius. However, 67P is composed of both dust and non-magnetic ice. To determine the magnetization limit of the dusty material, we incorporated the dust-to-ice mass ratio, ξ:

$$\begin{eqnarray*}&&{j}_{\max }\,=\,\displaystyle \frac{{B}_{\max ,\mathrm{meas}}}{{B}_{\max ,\mathrm{sim}}}\displaystyle \frac{\xi +1}{\xi }{j}_{\mathrm{sim}}.\end{eqnarray*}$$

The dust-to-gas mass ratio measured for 67P is 4 ± 2 (Rotundi et al. 2015). Additionally, experiments on silica and ice pebble compression during collision coupled with simulations of pebble-pile formation for 67P suggest that the measured porosity requires 3 ≲ ξ ≲ 9 (Lorek et al. 2016). We therefore conservatively (i.e., higher inferred magnetization) assume ξ = 3. The results for each simulation are shown in Figure 6. For r_d = 0.5 m, the average of the maximum magnetizations from our simulations is 2.6 × 10⁻⁶ ${\rm{A}}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$. The analogous value from Auster et al. (2015), corrected with the updated bulk density , the new constraint on the measured field intensity (Heinisch et al. 2018), and ξ = 3, is 1.6 × 10⁻⁵ ${\rm{A}}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$. At domain sizes smaller than one meter, we obtained j_max ∼ 5 × 10⁻⁶ ${\rm{A}}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$, while at scales of one meter and larger, j_max ∼ 2 × 10⁻⁶ ${\rm{A}}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$. The values are in broad agreement with the low specific magnetic moment inferred for 67P by Auster et al. (2015).

Zoom In Zoom Out Reset image size

If 67P retains ADRM, then, over a given region, magnetic grains that were accreted during the formation of the comet would have nearly aligned magnetic moments. The alignment of these moments results in material that is approximately uniformly magnetized, producing a bulk specific magnetization comparable to that of the individual grains, j_bulk ∼ j_gr. The length scale (r_d in our model) of these regions is poorly constrained by theory and could, depending on the details of 67P's formation, range from several times the average grain size (r_d ∼ 5 cm) to the size of the entire comet (r_d ∼ 2 km). Conversely, random orientation of magnetic grains results in a specific magnetization for bulk material that is significantly less than the grain magnetization, j_bulk ∼ (r_gr/r_d)^3/2(ρ_gr/ρ_bulk)^1/2j_gr. If the maximum j_bulk permitted by the Philae measurements for a given value of r_d is below that expected for uniform magnetization (i.e., the magnetization of individual grains) then this suggests either that ADRM is not recorded on 67P or that ADRM operates with low efficiency. ADRM may have lower than perfect efficiency as a result of spatial or temporal variations in the aligning magnetic field or grain misalignment during accretion.

The Philae data do not provide a significant constraint on the magnetization of small (≲1 mm) grains on 67P, but samples of other solar system bodies can provide a range of plausible values. Studies of material returned from the JFC 81P/Wild ("Wild 2") by the Stardust mission show that the rocky component of JFCs is similar to primitive ordinary chondritic material (Brownlee et al. 2012). Meteorite measurements have shown that chondritic material can possess a remanent magnetization spanning from ∼4 × 10⁻⁵ to ∼1 × 10⁻¹ ${\rm{A}}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$ (e.g., Terho et al. 1993). Note, however, that many chondrites with moments >0.01 ${\rm{A}}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$ were likely remagnetized by collectors' hand magnets (Weiss et al. 2010). Chondrules and millimeter-sized bulk samples of the Allende CV carbonaceous chondrite tend toward the low end of this range, 10⁻⁵–10⁻⁴ ${\rm{A}}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$ (e.g., Fu et al. 2014a). Finally, pieces of the matrix of the Semarkona LL3.00 ordinary chondrite that have been completely demagnetized by alternating fields in the laboratory have a specific magnetization of ≳6 × 10⁻⁵ ${\rm{A}}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$ (Fu et al. 2014b).

In comparison to measurements of meteorite samples, 67P is strikingly non-magnetic—permissible values of the bulk magnetization of 67P on scales ≥10 cm are generally an order of magnitude less than even the lowest values reported from meteoritic materials. This suggests either that 67P does not contain regions of coherent magnetization on these scales from ADRM in the areas where Philae made ground contact, or that ADRM operates with a low alignment efficiency, at most ∼10% and possibly much lower.

The first stage in ADRM acquisition is the alignment of nebular grains to the background magnetic field. We review the requirements for grain alignment as described by Fu & Weiss (2012), adapting them where necessary for conditions in the outer solar nebula. The "compass-needle" mechanism of grain alignment described here is distinct from the theory of radiative grain alignment, which is more commonly discussed in the protoplanetary disk literature (e.g., Tazaki et al. 2017; Hull et al. 2018). In radiative alignment torque theory, a spinning paramagnetic grain acquires a magnetic moment parallel to its spin axis through the Barnett effect. Grain alignment then proceeds from the precession of a spinning grain's angular momentum vector about the direction of either the magnetic field or the anisotropic radiation field, depending on the relative strengths of the magnetic and radiative torques. Over ∼10–100 precession periods, the angular momentum vector becomes aligned with the precession vector (Lazarian & Hoang 2007).

The compass-needle mechanism of grain alignment discussed in this paper differs from radiative alignment in two ways. First, alignment of a grain occurs after the grain has been brought to near rotational rest by interaction with nebular gas. Second, the magnetic moments of these grains are produced by the presence of magnetized ferromagnetic inclusions, not grain rotation. As a result, these grains possess magnetic moments that are orders of magnitude larger than those produced through the Barnett effect. Additionally, the magnetization direction for these grains is independent from the shape of the grain. As a result, grains that are aligned to the nebular field by this mechanism should have no preferred orientation of their long axes and are therefore not expected to produce a polarized emission signature. Grain alignment in the compass-needle framework requires first that the torque exerted by the magnetic field be stronger than other randomizing torques on the ferromagnetic grain. Second, for the grain to remain aligned, the time between collisions, which disrupt the grain orientation, must be significantly greater than the time required to bring the grain back into alignment.

In these calculations, we adopt the minimum-mass solar nebula model from Hayashi (1981) with the gas temperature (T_g), gas density (ρ_g), and gas scale height (H_g) defined as

$$\begin{eqnarray}&&{T}_{g}=280\,{\rm{K}}\times {\left(\displaystyle \frac{a}{1\mathrm{au}}\right)}^{-1/2},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\rho }_{g}=1.36\,\times \,{10}^{-6}\,\mathrm{kg}\,{{\rm{m}}}^{-3}\times {\left(\displaystyle \frac{a}{1\mathrm{au}}\right)}^{-11/4},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{H}_{g}=0.0472\,\mathrm{au}\times {\left(\displaystyle \frac{a}{1\mathrm{au}}\right)}^{5/4},\end{eqnarray} \tag{ 10 }$$

where a is the heliocentric distance in the midplane of the nebula. The scale height H_g defines the exponential vertical profile of the gas at a given heliocentric distance: at a height z above the disk midplane, ${\rho }_{g}(z)\propto \exp (-{z}^{2}/{H}_{g}^{2})$. The H/He nebular gas has a mean molecular weight of μ = 2.34 u.

4.1.1. Grain Torques

Magnetic Torque—For grain alignment, the magnetic torque must exceed other randomizing torques on the grain. Depending on the grain size, the particular non-magnetic torque that dominates varies (Figure 8).

Zoom In Zoom Out Reset image size

The magnetic torque exerted by a field with intensity B on a spherical grain with a specific magnetization j_gr is

$$\begin{eqnarray}&&{\tau }_{B}=\displaystyle \frac{4\pi }{3}{r}_{\mathrm{gr}}^{3}{\rho }_{\mathrm{gr}}{j}_{\mathrm{gr}}B,\end{eqnarray} \tag{ 11 }$$

where r_gr and ρ_gr are the radius and density of the grain, respectively. We assume ρ_gr = 2000 kg m⁻³ in all our calculations. We adopt a conservative value of j_gr = 10⁻⁵ ${\rm{A}}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$ for the grain magnetization (see Section 3.3).

Brownian Motion—Very small grains are subject to randomization through rotational Brownian motion. Following Fu & Weiss (2012), we consider a grain to be misaligned when the thermal energy of rotation (k_BT/2, where k_B is the Boltzmann constant) exceeds the difference in magnetic potential energy for an aligned grain and one misaligned by an angle θ = π/10. Using the small-angle approximation, the threshold for alignment is then

$$\begin{eqnarray}&&\displaystyle \frac{4\pi }{6}{r}_{\mathrm{gr}}^{3}{\rho }_{\mathrm{gr}}{j}_{\mathrm{gr}}B{\left(\displaystyle \frac{\pi }{10}\right)}^{2}\gtrsim \displaystyle \frac{{k}_{{\rm{B}}}{T}_{g}}{2}.\end{eqnarray} \tag{ 12 }$$

Solving this inequality for the magnetic field intensity (B) and substituting into Equation (11), the magnetic torque required to align grains against Brownian randomization is ${\tau }_{{\rm{B}}}\gtrsim {\left(10/\pi \right)}^{2}{k}_{{\rm{B}}}{T}_{g}$.

Aerodynamic Drag—Grains in the solar nebula experience aerodynamic drag as a consequence of their motion relative to the nebular gas. The drag regime depends on the ratio of the mean free path of gas molecules to the size of the grain. This ratio is defined as the Knudsen number,

$$\begin{eqnarray}&&{\rm{Kn}}=\displaystyle \frac{1}{2{r}_{{\rm{gr}}}{n}_{g}\sigma },\end{eqnarray} \tag{ 13 }$$

where n_g = ρ_g/μ is the gas number density and σ = 2 × 10⁻¹⁵ cm² is the collision cross section of H₂ (Nakagawa et al. 1986). We define grains with Kn > 1/8 to be in the free molecular (FM) drag regime, where the mean free path is large compared to the grain size. All the grains of interest in this study fall into this regime. The net torque on a grain in an FM flow is determined by integrating the effects of gas molecule collisions on each surface element of the grain over the entire surface area. As in Fu & Weiss (2012), we assume the grains are prolate ellipsoids with an aspect ratio χ, defined as the ratio of the lengths of the semimajor and semiminor axes. For comparison to other torques, we consider r_gr as the semiminor axis, so that χr_gr is the semimajor axis. Expressions for the FM torque on an arbitrary convex figure of revolution, assuming a Maxwellian velocity distribution for the gas and thermal equilibrium between the gas and dust grain, are provided by Equation (3) from Ivanov & Ianshin (1980). We use their expression for M_x and assume an angle of π/4 between the long axis of the grain and the gas velocity vector (see Appendix B for details).

The relative velocity between the gas and dust grain has two components. First, the systemic drift velocity, caused by sub-Keplerian motion of the gas, and, second, velocity differences caused by turbulence in the gas. We determine the total systemic velocity from the radial (u) and transverse (w) velocity components, ${v}_{\mathrm{sys}}=\sqrt{{u}^{2}+{w}^{2}}$, where u and w are determined by solving Equations (15) and (20) in Weidenschilling (1977), as detailed in Appendix B.

The relative velocity induced by turbulence is

$$\begin{eqnarray}&&{v}_{{\rm{turb}}}={v}_{g}{\left[\displaystyle \frac{{{\rm{St}}}_{L}^{2}\left({{\rm{Re}}}^{1/2}-1\right)}{\left({{\rm{St}}}_{L}+1)({{\rm{St}}}_{L}{{\rm{Re}}}^{1/2}+1\right)}\right]}^{1/2},\end{eqnarray} \tag{ 14 }$$

St_L is the Stokes number relative to the largest eddy, Re is the flow Reynolds number, and v_g is the mean square of the turbulent gas velocity (Cuzzi & Hogan 2003). St_L is given by the ratio of the stopping time of a grain due to gas drag to the overturn time of the largest eddy, which is comparable to the orbital period:

$$\begin{eqnarray}&&{{\rm{St}}}_{L}=\displaystyle \frac{{r}_{{\rm{gr}}}{\rho }_{{\rm{gr}}}}{{v}_{{\rm{th}}}{\rho }_{g}}\displaystyle \frac{1}{{P}_{{\rm{orb}}}},\end{eqnarray} \tag{ 15 }$$

where ${v}_{\mathrm{th}}=\sqrt{8{k}_{{\rm{B}}}{T}_{g}/(\pi \mu )}$ is the thermal velocity of the gas (Nakagawa et al. 1986). The Reynolds number is given by

$$\begin{eqnarray}&&\mathrm{Re}\,=\,\displaystyle \frac{{{Lv}}_{g}{\rho }_{g}}{\eta }=\displaystyle \frac{\alpha {H}_{g}{\rho }_{g}\sigma }{\mu }\end{eqnarray} \tag{ 16 }$$

where η ∼ ρ_gv_th/(n_gσ) = v_th/(μσ) is the dynamic molecular viscosity of the nebular gas and $L\sim \sqrt{\alpha }{H}_{g}$ and ${v}_{g}\sim \sqrt{\alpha }{v}_{\mathrm{th}}$ (Cuzzi et al. 2001). The latter two expressions are derived from the assumption of an α-disk model, where angular momentum transport in the disk is assumed to be driven by a turbulent viscosity ∼αv_thH_g. Values for the dimensionless turbulence parameter α used here are 10⁻⁴ and 10⁻², corresponding to either a quiescent or a turbulent nebula. To find the total relative velocity between grains and the nebular gas, we make the conservative assumption that $v=| {v}_{\mathrm{sys}}| +| {v}_{\mathrm{turb}}|$ instead of adding the components in quadrature.

Accommodation Torque—Grains in the nebula can also experience torques as a consequence of surface heterogeneities. Differences in the efficiency of momentum transfer (or accommodation) during gas molecule collisions across the surface of a grain can result in a rotation rate higher than expected from random impacts of gas molecules (Purcell 1979). The magnitude of the effective torque, from Equation (14) of Purcell, is

$$\begin{eqnarray}&&{\tau }_{\mathrm{accom}}=\sqrt{\displaystyle \frac{2}{3}}{r}_{\mathrm{gr}}^{2}{n}_{g}{k}_{{\rm{B}}}{\rm{\Delta }}T\delta s,\end{eqnarray} \tag{ 17 }$$

where δ is the magnitude of the variation in the accommodation coefficient, s is the length scale over which the accommodation varies, and ${\rm{\Delta }}T\equiv {T}_{g}-\sqrt{{T}_{g}{T}_{s}}$ describes the difference between the grain temperature (T_s) and the gas temperature and is parameterized by η_a = ΔT/T_g.⁴ Following Purcell (1979) and Fu & Weiss (2012), we adopt s = 10⁻⁷ m, and the conservative values δ = 0.1 and η_a = 1.

Purcell Torque—Like the accommodation torque, the Purcell torque also arises as a consequence of variations in the properties of the grain surface. At certain sites on the grain, the binding energy of incoming hydrogen atoms is enhanced so that they become more tightly bound to the grain. The arrival of a second hydrogen atom can then result in recombination and the ejection of the resulting hydrogen molecule, causing a reaction force on the grain. The random distribution of these sites of hydrogen formation and ejection over the surface of the grain gives rise to a torque

$$\begin{eqnarray}&&{\tau }_{P}=\sqrt{\displaystyle \frac{32}{3\pi }}\gamma {r}_{\mathrm{gr}}^{3}{n}_{{\rm{H}}}\sqrt{\displaystyle \frac{{k}_{{\rm{B}}}{T}_{g}{E}_{\mathrm{ej}}}{\nu }},\end{eqnarray} \tag{ 18 }$$

where γ is the fraction of inbound hydrogen atoms that are ejected as molecules, n_H is the number density of hydrogen atoms, ν is the number of enhanced sites on the grain, and E_ej = 0.2 eV is the mean kinetic energy with which the molecule is ejected (Purcell 1979, Equation (19)). We adopt γ = 0.2 and ν = 100 × (r_gr/10⁻⁷ m)² from Fu & Weiss (2012) as values appropriate for the solar nebula. At the temperatures expected in the midplane region where 67P likely formed, the gas is dominated by molecular hydrogen, so that the abundance of atomic hydrogen is n_H/n_g ≲ 10⁻⁴ (e.g., Glassgold et al. 2009). We therefore take n_H/n_g = 10⁻⁴ as a conservative assumption in our calculations. Given the weak magnitude of the Purcell torque compared to other effects across the grain sizes of interest, the uncertainty in these parameters does not affect our results.

Radiative Torque—The final effect we consider is the radiative torque, affecting irregularly shaped grains in a uniform radiation field. From Equation (29) of Fu & Weiss (2012) the approximate magnitude of the torque is

$$\begin{eqnarray}&&{\tau }_{\mathrm{rad}}\approx \displaystyle \frac{1}{2}{r}_{\mathrm{gr}}^{2}{u}_{\mathrm{rad}}\bar{\lambda }{Q}_{{\rm{\Gamma }}},\end{eqnarray} \tag{ 19 }$$

where $\bar{\lambda }$ is the mean wavelength of the incident radiation, Q_Γ is an efficiency coefficient that depends on the grain size, and u_rad is the volumetric energy density of the radiation. In the interior of the disk, the radiation field is dominated by thermal emission from cool grains (e.g., Tazaki et al. 2017). Following Fu & Weiss, we therefore adopt $\bar{\lambda }\sim 100\,\mu {\rm{m}}$, Q_Γ ∼ 1 for ${r}_{\mathrm{gr}}\geqslant \bar{\lambda }$ and ${Q}_{{\rm{\Gamma }}}={({r}_{\mathrm{gr}}/\bar{\lambda })}^{3}$ for ${r}_{\mathrm{gr}}\lt \bar{\lambda }$, and a radiation density

$$\begin{eqnarray}&&{u}_{\mathrm{rad}}=\displaystyle \frac{8{\pi }^{5}{k}_{{\rm{B}}}^{4}}{15{c}^{3}{h}^{3}}{T}_{\mathrm{gr}}^{4},\end{eqnarray} \tag{ 20 }$$

where c is the speed of light and h is Planck's constant. Following Chiang & Goldreich (1997), we assume thermal equilibrium between the dust and gas in the disk midplane, T_gr ∼ T_g. This is a conservative choice, because it results in an overestimate of the radiation density by roughly an order of magnitude when compared to the more sophisticated model of Tazaki et al. (2017).

Considering all the torques above, we determine the magnetic field intensity required for the magnetic torque (Equation (11)) to be the strongest torque acting on grains of varying size (Figure 9). This provides a lower bound on the magnetic field intensity required to magnetically align the grains. For grains with r_gr ≲ 0.1 mm we find that alignment is determined by competition between Brownian motion and the magnetic torque, while for r_gr ≳ 0.1 mm the dominant non-magnetic torque is aerodynamic drag. We consider two cases for aerodynamic drag: a lower bound case with nearly circular grains and a quiescent nebula (χ = 1.1 and α = 10⁻⁴), and an upper bound case with prolate grains and a highly turbulent nebula (χ = 2 and α = 10⁻²).  The required field intensity is strongly dependent on the grain radius, but even in the higher drag scenario, the magnetic torque will control alignment of millimeter-sized grains for fields B ≳ 1 μT.

Zoom In Zoom Out Reset image size

4.1.2. Grain Alignment and Collision Timescales

Even when the magnetic torque dominates over other forces acting on nebular grains, alignment is disrupted by collisions with other grains. These collisions randomize the orientation and spin up the grains. When a grain is rapidly spinning (ω ∼ 10² s⁻¹), the rotational energy is much greater than the work done by the magnetic torque in one revolution, making direct realignment impossible and causing the grain to instead precess about the field direction (Fu & Weiss 2012). Realignment is therefore a two-stage process. First, on a timescale t_prec, the grain must shed its rotational energy through gas drag. Subsequently, the magnetic torque can reorient the grain, and the grain's magnetic moment oscillates about the field direction, with the gas drag now serving to damp the oscillation. This damping occurs on a timescale of t_drag. Grain alignment in the solar nebula requires that the time between collisions be much longer than the time required to de-spin and realign the grain: t_coll ≫ t_align = t_drag + t_prec. From Fu & Weiss (2012), the collision, drag, precession, and alignment timescales are

$$\begin{eqnarray}&&{t}_{\mathrm{coll}}=\displaystyle \frac{{\rho }_{\mathrm{gr}}}{3{v}_{\mathrm{rel}}{\rho }_{s}}{r}_{\mathrm{gr}},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{drag}}=\displaystyle \frac{4{\rho }_{\mathrm{gr}}}{5f{\rho }_{g}{v}_{\mathrm{th}}}{r}_{\mathrm{gr}},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray*}&&{t}_{\mathrm{prec}}=-\displaystyle \frac{{t}_{\mathrm{drag}}}{2}\mathrm{log}\left(\displaystyle \frac{5{{Bj}}_{\mathrm{gr}}}{{r}_{\mathrm{gr}}^{2}{\omega }_{0}^{2}}\right),\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{align}} & = & {t}_{\mathrm{drag}}+{t}_{\mathrm{prec}}\\ & = & \ {t}_{\mathrm{drag}}\left[1-\displaystyle \frac{1}{2}\mathrm{log}\left(\displaystyle \frac{5{{Bj}}_{\mathrm{gr}}}{{r}_{\mathrm{gr}}^{2}{\omega }_{0}^{2}}\right)\right],\end{array}\end{eqnarray} \tag{ 23 }$$

where v_rel is the relative velocity between grains in the nebula, ρ_s is the density of solids in the nebula, and ω₀ is the initial angular velocity of a grain after a collision. The factor f represents the efficiency of momentum transfer during gas molecule collisions and is assumed here to be 1. For conditions in the outer solar nebula, we find that t_prec ∼ 10 t_drag, so that t_align is roughly an order of magnitude greater than t_drag. This conclusion is insensitive to assumptions about the magnetic field strength and grain magnetization owing to the logarithmic dependence on these factors. Fu & Weiss (2012) show that in the solar nebula t_coll ≫ t_align for grains with r_gr ≲ 10 cm, suggesting that, if the magnetic torque dominates grain alignment, small grains will be aligned to the background field most of the time.

4.1.3. Upper Limit on Field Intensity

If the absence of ADRM on 67P is due to the failure of grain alignment in the solar nebula, then the calculated field intensity at which the magnetic torque dominates other forces (Figure 9) provides an upper bound on the nebular field intensity at the time and location of 67P's formation.

The location of 67P's formation is not well constrained. Comets are believed to have formed in a disk of icy planetesimals in the region of ∼15–45 au before being scattered by interactions with the giant planets to form the scattered disk and the Oort cloud (Fernandez & Ip 1981; Brasser & Morbidelli 2013; Johansen et al. 2014). Additional clues to the formation environment come from 67P's ice and volatile inventory. The release of volatiles from 67P is consistent with the decomposition of H₂O clathrates, suggesting that at least some of the water ice on 67P is not amorphous (Rubin et al. 2015; Luspay-Kuti et al. 2016; Mousis et al. 2016a, 2016b). Crystalline ice and clathrates are expected to form in the solar nebula through the vaporization of amorphous ice at temperatures of ∼150 K followed by condensation as the disk cools, requiring formation at distances ≲30 au (Willacy et al. 2015). Results from the Stardust mission, however, suggest substantial radial transport of material prior to or during the epoch of comet formation (Brownlee et al. 2012), so the presence of these ices may not directly reflect the formation location. Given these uncertainties, we calculate the minimum field intensity required to enforce grain alignment assuming a formation location of 15–45 au. For a given grain size, the results show that there is little change in calculated field intensities over this range of heliocentric distances.

However, the required field intensity is strongly influenced by the assumed size of the grains that agglomerated to form 67P (Figure 9). As discussed below (Section 4.2), this value is in turn highly dependent on the assumed mechanism by which the comet formed. Multiple lines of evidence, including 67P's thermal inertia, porosity, and surface features, suggest that 67P is composed largely of pebbles with a characteristic radius of 3–6 mm (Fulle et al. 2016; Poulet et al. 2016; Blum et al. 2017). Assuming, for the moment, that this approximates the size of nebular grains that accreted to form 67P, the absence of ADRM requires a magnetic field intensity B ≲ 1 μT for r_gr = 1 mm and B ≲ 3 μT for r_gr = 3 mm (Figure 10).

Zoom In Zoom Out Reset image size

The predicted alignment timescale for ferromagnetic grains (≲1 yr) is much shorter than the timescale for variation in the outer disk magnetic field found in simulations (∼10–100 yr) (e.g., Bai 2015). Absence of grain alignment on 67P therefore constrains the instantaneous magnetic field in the outer disk, not the weaker time-averaged net field (so that the orange and red curves in Figure 10 are the relevant curves to be compared to the Philae constraints). Our derived limits on field intensity are broadly consistent with theoretical predictions of the instantaneous magnetic field intensity required to drive stellar accretion at the assumed rate of ∼10⁻⁸ M_⊙ yr⁻¹ (red and orange lines in Figure 10).

Observations of (sub)millimeter thermal emission from dust grains in other protoplanetary disks offer a complementary way to measure magnetic field intensities. The polarization of this emission is related to the alignment of grains with sizes comparable to the wavelength, so the intensity and morphology of the observed polarization can indicate the mechanism that is controlling grain alignment. Recent observations of the protoplanetary disks around HL Tau and IM Lup show a pattern of polarization consistent with either self-scattering by unaligned dust grains, or radiative, not magnetic, alignment of dust grains (e.g., Kataoka et al. 2015, 2017; Yang et al. 2016; Stephens et al. 2017; Tazaki et al. 2017; Hull et al. 2018). Tazaki et al. (2017) estimate that a magnetic field intensity ≳13 μT is required for magnetic alignment to supersede radiative alignment for millimeter-sized grains at 50 au, consistent with, but less restrictive than, the upper bounds on the field intensity we infer from 67P. Tazaki et al. (2017), however, primarily considered non-magnetic grains or grains with super-paramagnetic inclusions. We consider grains with ferromagnetic inclusions, which can have much larger magnetic moments and can therefore become aligned at lower field intensities. Additionally, present-day observations of polarized emission from protoplanetary disks have typical resolution of ∼10 au, making it difficult to probe the inner regions of the disk where comet formation likely occurred. The 67P measurements therefore provide both a tighter constraint on the magnetic field and insight into a region of the disk that is difficult to access with current astronomical measurements.

Even under conditions where nebular grains are aligned to the background field, the mechanism by which these grains accrete to form a comet is critical for determining whether ADRM is possible. For grain alignment to be preserved and produce a magnetic record, the grains must be directly accreted by the comet while maintaining alignment.

There are two main hypothesized mechanisms for the formation of comets: hierarchical growth and pebble cloud collapse. If comets undergo hierarchical growth, forming gradually through collisional coagulation of increasingly large aggregates (e.g., Weidenschilling 1997; Davidsson et al. 2016), then ADRM is likely impossible or extremely limited in spatial scale. For plausible magnetic field intensities, the orientation of grains ≳1 cm in radius will be dominated by other factors, so that pairwise collisions of particles larger than this size will produce random mutual orientations of their magnetic moments. This scenario could explain the absence of ADRM on 67P at spatial scales ≥10 cm even in the presence of an intense background magnetic field.

Alternatively, recent modeling has demonstrated that the concentration of nebular solids through the streaming instability (Johansen et al. 2007), followed by the gentle gravitational collapse of the resulting cloud of millimeter- or centimeter-sized grains, can produce comet-like planetesimals matching the observed density, porosity, and tensile strength of comets (Blum et al. 2014; Wahlberg Jansson & Johansen 2014, 2017; Lorek et al. 2016). The detection of millimeter-sized pebbles on 67P, by direct imaging among other techniques, combined with the primitive nature of cometary material and the preservation of fluffy fractal dust aggregates, has been taken as evidence that 67P formed directly from the pebbles observed today (Fulle et al. 2016; Poulet et al. 2016; Blum et al. 2017; Fulle & Blum 2017). Blum et al. (2017) and Fulle & Blum (2017) argue that this evidence disfavors hierarchical growth models, which would compress fluffy fractal dust particles and would result in a range of grain sizes from small grains (∼10 μm) to larger aggregates, instead of the observed characteristic size of a few millimeters. In the pebble-pile scenario, because 67P directly incorporates millimeter-sized grains as it accretes, the absence of ADRM cannot be easily explained by the formation mechanism alone, but would also require the background field to be weak.

A possible factor limiting ADRM formation in the pebble pile is the repeated grain collisions in the pebble cloud preventing grain alignment (see Section 4.1.2). This is because the pebble cloud is much denser than the background nebula we had previously considered. To determine whether alignment can be maintained in these overdense regions, we recalculate the collision and alignment timescales (Equations (21) and (23)) for conditions at the start of cloud collapse. Following Lorek et al. (2016) and Wahlberg Jansson & Johansen (2014), the initial pebble cloud has uniform density, a mass M_67P ≈ 9.982 × 10¹² kg (Pätzold et al. 2016), and a radius equal to its Hill radius, R_H. We assume millimeter-sized grains, a formation distance of 20 au, a background field intensity B = 0.1 μT, and grain magnetization j_gr = 10⁻⁵ ${\rm{A}}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$. We take the velocity dispersion to be ${v}_{\mathrm{rel}}\sim \sqrt{2}{v}_{\mathrm{vir}}$, where the virial velocity is ${v}_{\mathrm{vir}}\sim \sqrt{{{GM}}_{67{\rm{P}}}/{R}_{{\rm{H}}}}$, and we conservatively estimate that the initial grain rotation rate after a collision is given by ω₀ = v_rel/r_gr. We find that the criterion for alignment is still satisfied in the initial pebble clouds: t_coll ≈ 19 yr ≫ 4 yr ≈ t_align.

As the cloud collapses, grain collisions continue. Initial models (Wahlberg Jansson & Johansen 2014) assumed a homogeneous particle size distribution in the cloud, leading to a steady increase in cloud density until solid density was reached. Under these conditions, the collision timescale rapidly decreases at the end of collapse, possibly randomizing grain orientations. More recent models (Wahlberg Jansson & Johansen 2017), which radially resolve the accreting planetesimals and include the effects of gas drag, suggest that comet-sized planetesimals may form through the rapid accretion of a core followed by the protracted accretion of a remnant low-density cloud. Because of the low density, collision rates during the end stage are low, so these grains may be able to maintain alignment to a background field. A caveat is that even these updated models lack important physics, including the rotation of the collapsing cloud. In summary, the absence of ADRM on 67P could be explained by a hierarchical growth model, but if comets form instead through the gentle gravitational collapse of a cloud of pebbles, then the formation mechanism alone may not explain the lack of magnetic alignment of 67P's constituent grains.

Significant post-formation restructuring would also provide a natural explanation for the observed lack of ADRM on 67P. There is some evidence that, rather than primordial objects, comets may be collisional rubble piles, forming from the re-accretion of material from a larger parent body (Morbidelli & Rickman 2015; Rickman et al. 2015). Others argue that 67P's high porosity, low tensile strength, abundance of supervolatiles, absence of aqueous alteration, and preservation of fractal dust particles require a primordial origin (Davidsson et al. 2016; Blum et al. 2017; Fulle & Blum 2017). Recent modeling of collisions, however, has shown that at least some primitive properties of cometary material can be preserved through significant collisional restructuring, including during disruption of the parent body and subsequent re-accretion (Jutzi et al. 2017; Schwartz et al. 2018). These scenarios do not predict global crushing or compaction on sub-meter scales, which would be difficult to reconcile with the retention of supervolatiles and high porosity, although these effects are well below the ∼10 m resolution limits of the models. Collisions may be capable of erasing ADRM on ≳10 m scales, but magnetization on ≲1 m scales may persist, even if the bulk structure of 67P is not primordial.

Outgassing, sublimation of volatiles, erosion, and sedimentation may all change the structure of 67P on smaller scales, resulting in the scrambling of grain orientations and loss of ADRM. Because our results are most sensitive to the measurements closest to the comet's surface, a regolith layer ∼1 m thick could mask ADRM if it were recorded on only ∼10 cm scales. Much thicker regolith layers are required, however, to mask magnetic coherence on greater length scales. The Philae data include low-altitude magnetic field measurements at multiple contact points on the nucleus. The initial touchdown location was in a region with an estimated 20–30 cm thick layer of granular regolith and meter-sized pits overlying more consolidated material (see discussion of "pitted plains" in Birch et al. 2017), while the final touchdown location is believed to be primordial, displaying a relatively dust-free landscape with nearby blocks and boulders and stratified terrain (Lucchetti et al. 2016; Poulet et al. 2016). The remanent field intensity is ≤0.9 nT at all contact points, regardless of the thickness of regolith or apparent degree of alteration. Birch et al. (2017) propose that observed comet nuclei display a range of alteration, with smoother terrain present on more altered comets. They argue that 67P displays more primordial morphology than most other JFCs. If the observed radial-layered structure is also primordial (Massironi et al. 2015), then this also places tight constraints on the degree of structural alteration.

Surface processes may be able to erase, or mask, ADRM on 67P if it was only ever acquired on small length scales. By contrast, if 67P, or its parent body, acquired large-scale ADRM during formation, a combination of collisions and surface processing would be required to remove it.

ADRM could also be removed by the remagnetization of ferromagnetic minerals after accretion. Processes that would be capable of remagnetizing 67P include heating to temperatures exceeding several hundred kelvin, aqueous alteration, or impact shocks exceeding ∼1 GPa (Gattacceca et al. 2007). However, these conditions are more extreme than those predicted by models of collisional rubble piles (e.g., Jutzi et al. 2017; Schwartz et al. 2018), and 67P shows no evidence of the extensive and widespread heating, aqueous alteration, or large shocks that would be required to demagnetize it over the entire region surveyed by Philae (e.g., Davidsson et al. 2016; Fulle & Blum 2017). We therefore consider this an unlikely explanation.