Probing Coronal Magnetic Fields with Sungrazing Comets: H i Lyα from Pickup Ions

To compute the rate of Lyα scattering we precompute a lookup table of the number of photons scattered per second per H atom as a function of the flow speed in the radial direction and thermal width, then scale by the dilution factor of the solar radiation at the heliocentric height of the comet. For the integration over frequency and the illuminating disk, we use a calculation by S. Cranmer (2019, private communication). We use a Lyα disk intensity appropriate for solar minimum and increase it by a factor of about 1.8 for solar maximum conditions. For first generation neutrals we include the opacity in the Lyα line using the outflow speed of the first generation neutrals, V₁, for the Doppler width, and we reduce the scattering rate accordingly. V₁ is determined by the energy liberated during photodissociation of H₂O and OH, and ranges from about 8 to 24 km s⁻¹ for the H atoms (Shimizu 1991). For the opacity calculation, we assume that B is radial.

The ionization rate is given by

$$\begin{eqnarray}&&{q}_{\mathrm{ion}}={n}_{p}{q}_{\mathrm{ex}}\,+\,{n}_{e}{q}_{i}\,+\,{q}_{\mathrm{phot}},\end{eqnarray} \tag{ 1 }$$

where q_ex is the charge transfer rate coefficient, q_i is the ionization rate coefficient for electron impact, q_phot is the photoionization rate, and n_p and n_e are the coronal proton and electron densities (we use the approximation n_p = 0.83n_e valid for solar coronal plasma). The charge transfer rate is computed with the cross sections of Schultz et al. (2008) with the average velocity. The effective speed for charge transfer between first or second generation neutrals and coronal protons or PUIs is given by the sum in quadrature of the relative bulk speed of the two components and the random (thermal) speeds of the two components. For instance, the effective speed for charge transfer between first generation neutrals and coronal protons is given by

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}=\sqrt{{V}_{\mathrm{rel}}^{2}+{V}_{T}^{2}+{V}_{1}^{2}},\end{eqnarray} \tag{ 2 }$$

where V_T is the thermal speed of the coronal protons (computed as the most probable speed of protons at T_cor) and V_rel is the bulk speed of the coronal gas relative to the comet. V₁ is the outflow speed of the first generation neutrals, which is small compared to the other terms. We ignore charge transfer between second generation neutrals and coronal protons, because that does not change the number of second generation neutrals or their speed distribution. For charge transfer between third generation neutrals and coronal protons, the effective speed for charge transfer becomes

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}=\sqrt{{V}_{\mathrm{rel}}^{2}+{V}_{T}^{2}+{V}_{\perp }^{2}}.\end{eqnarray} \tag{ 3 }$$

For the third generation, the value of V_rel uses the parallel velocity component of the PUIs and the solar wind speed.

For collisional ionization, we use the ionization rate from Scholz & Walters (1991) at an electron temperature of log T = 6.2. It varies by only 3% over the range of log T = 6.0–6.6. For photoionization, we scale the rate at 6.8 R_⊙ from Raymond et al. (1998) with the dilution factor as a function of heliocentric distance. This is the value appropriate for solar minimum, and we increase it by a factor of 2 for solar maximum conditions.

The model uses a Cartesian grid, with the comet moving along the X-axis and the magnetic field in the XZ plane. The models can then be rotated to give Doppler shifts and foreshortening for comparison with observations. Neutrals are generated at the comet and followed as they scatter Lyα and undergo charge transfer until they become ionized or leave the simulation box. We use analytic descriptions of the distributions of first and second generation neutrals. PUIs are produced by ionization of first generation neutrals (or by ionization of third generation neutrals). Third generation neutrals are produced by charge transfer between PUIs and first or second generation neutrals at each grid cell, and assumed to move in the direction of the magnetic field at ${V}_{\parallel }$ and to spread out with a speed V_⊥.

We approximate the first generation distribution as a spherically symmetric outflow moving with the comet and expanding at a constant speed, with an exponential cutoff due to ionization. Its density is

$$\begin{eqnarray}&&{n}_{1}(r)=\displaystyle \frac{\dot{N}}{4\pi {r}^{2}{V}_{1}}{e}^{-\tfrac{{q}_{\mathrm{ion}}}{{V}_{1}}r},\end{eqnarray} \tag{ 4 }$$

where V₁ is the initial outflow speed due mostly to the 8–24 km s⁻¹ speeds of the H atoms ejected during photodissociation. At high densities, the H atoms share some of their momentum with O atoms, so V₁ can be smaller.

The second generation neutrals form a conical structure whose axis lies between the comet trajectory and the solar wind direction. The angle between the comet trajectory and the second generation axis is given by

$$\begin{eqnarray}&&\sin (\beta )={V}_{\mathrm{wind},\perp }/{V}_{\mathrm{rel}},\end{eqnarray} \tag{ 5 }$$

where V_wind,⊥ is the solar wind speed component perpendicular to the comet trajectory. We approximate the second generation neutrals as a cone with a density falling off with distance from the axis based on a Gaussian distribution with the coronal thermal speed V_T. The tail fades away on a length scale ${V}_{\mathrm{rel}}/({q}_{i}+{q}_{\mathrm{phot}})$.

The PUIs, and therefore the third generation neutrals, form a structure that moves along the B field with ${V}_{\parallel }={V}_{\mathrm{com}}\sin \theta$. Since ${V}_{\parallel }$ is less than V_com and has a component perpendicular to the comet trajectory, the PUI structure appears to trail behind the comet at an angle, displaced to the side opposite to that of the second generation neutrals. Because the cloud of first generation neutrals is fairly small owing to the relatively small V₁ and the PUIs are tightly coupled to the field lines on which they form, the main PUI structure is narrow. The third generation neutrals, however, can cross field lines, and they have a randomly directed velocity ${V}_{\perp }$, so the third generation neutrals form a wider structure around the PUI field lines.

The comet moves across the model grid at one cell per time step. The first and second generation densities are analytic functions that move with the comet. At each time step the number of PUIs generated at each point is computed from the density of first generation neutrals and the total (ionization plus charge transfer) ionization rate, and the PUIs are advected along the magnetic field at ${V}_{\parallel }$ and allowed to diffuse along it at ${V}_{\perp }/\sqrt{3}$. They are assumed to be unable to move across field lines. The densities of PUIs and first and second generation neutrals are used to compute the formation rate of third generation neutrals, and they also advect along the field lines, but they are allowed to diffuse in three dimensions with a speed V_⊥.

Once the densities of the different components are known, the Lyα emissivity at each grid point is computed using the central velocity and line width to compute the scattering rate considering Doppler dimming and distance from the Sun. At each grid point the first component emission is reduced by ${e}^{-\tau }$, where τ is the optical depth between the grid point and the Sun. The second and third generation neutrals are Doppler shifted away from the first component, and they generally have much larger line widths, so they are much less absorbed. Since the velocity centroid of each component is known, the average velocity of the scattered photons is the weighted average of the velocities of the three components. For comparison with observations, we rotate the model to the angle between the comet path and the LOS, then integrate to produce a 2D projection.

Figure 4 shows a model for the parameters R = 8.0 R_⊙, V_wind = 200 km s⁻¹, and n_e = 10⁴ cm⁻³, with θ = 23°. The intensity images show a very bright spherical first component, a conical second component, and a third component that juts off to one side.

Zoom In Zoom Out Reset image size

We have computed approximate models to show the overall nature of the spatial and velocity structure of the Lyα emission from a comet. We assume that the comet does not perturb the corona in a major way, such as inducing a bow shock like those measured near larger comets (Gombosi et al. 1996). One criterion is that the amount of mass produced by the comet must be small compared with the amount of coronal mass it interacts with. The former is just the outgassing rate. The amount of coronal mass swept up per second is

$$\begin{eqnarray}&&\dot{M}\,=\,4\pi {r}_{i}^{2}{\mu }_{c}{n}_{p}{V}_{\mathrm{rel}},\end{eqnarray} \tag{ 6 }$$

where r_i is the interaction length scale for ionization of cometary neutrals, V₁/q_i, and μ_c is the mean weight of the coronal nuclei (Jones et al. 2018). This gives a $\dot{N}\,\lt \,{10}^{35}{V}_{200}/{n}_{p}$ for V₁ = 10 km s⁻¹, where V₂₀₀ is the comet speed in units of 200 km s⁻¹. This is consistent with the parameter ranges of interest, but the next level of approximation would be very challenging. We defer it to a future investigation because it would require at least an MHD code of the level of sophistication of the ones used by Gombosi et al. (1996) and Jia et al. (2014).

Our assumption that the PUIs can be described by a shell in velocity space moving along the magnetic field at ${V}_{\parallel }$ is probably reasonable for perpendicular geometries, but if θ < 45° streaming instabilities are likely to set in that would drive the PUI distribution to a Maxwellian and couple it to the coronal plasma. We also assume a spherical shell in velocity space rather than the bispherical shell predicted by theory, and we neglect the loss of energy from PUIs to Alfvén waves, which depends on the Alfvén speed and θ.

Our method of advancing time does not allow second generation neutrals to get ahead of the comet, which can happen if V_com and V_wind are small compared to V_T. This is not too serious a problem for most parameters, because V_com is large close to the Sun, and V_wind becomes larger farther away from the Sun.

We also neglect the possible production of fourth and higher generation neutrals by subsequent charge transfer events. The probability of such subsequent charge transfers increases with density, and therefore with $\dot{N}$, so the assumption that the comet does not disturb the corona may break down before the neglect of the fourth and higher generations becomes a serious problem. Similarly, we neglect the loss of first and second generation neutrals to charge transfer with PUIs, which can become important at high $\dot{N}$.

To address these limitations, we would probably need an MHD model that includes a neutral fluid, as in Jia et al. (2014), along with a PUI fluid. The model from such a code would serve as the basis for a Monte Carlo calculation of the generation and transport of the neutrals and a calculation of the scattering of Lyα.

UVCS observed comet C/2002 S2 at 6 R_⊙ from 23:17 until 00:16 UT on 2002 September 18/19. We assume the comet position and velocity vector from the orbit as given in Giordano et al. (2015), and we take the magnetic field to be radial. That ought to be a good approximation between the top of the closed field arcades at the streamer cusps, which are up to about 4 R_⊙ (Strachan et al. 2002), and the Alfvén radius at 10–20 R_⊙ (Zhao & Hoeksema 2010; Deforest et al. 2014; Tasnim et al. 2018). That assumption is critical because the angle between the comet trajectory and the field determines the perpendicular and parallel velocities of the PUIs and third generation neutrals. Moreover, the direction of the field determines the LOS components of the second and third generation neutrals. The solar wind velocity is not independently known, but to the extent that the comet orbit is reliable, we can take the agreement of the predicted and observed redshifts of the third generation emission as some support for the assumption that the field is radial.

We assume a solar wind velocity parallel to the magnetic field with a magnitude of 175 km s⁻¹. This is much larger than the value of 75 km s⁻¹ inferred from the shape of the Lyα intensity image by Giordano et al. (2015), but we find that it is required to match the observed blueshift of 111–145 km s⁻¹. This large velocity implies strong Doppler dimming, so a much larger mass loss rate, 1 × 10²⁹ H per second, is needed to match the Lyα brightness. We have no independent measurement of the wind velocity at this height, but Cho et al. (2018) give a range of about 130–430 km s⁻¹ at 6 R_⊙ based on white light observations. The density chiefly controls falloff of the second generation intensity; the simulation with a coronal proton density of 10⁴ cm⁻³, in agreement with Giordano et al. (2015), agrees reasonably well with the data (Figure 9).

The assumed thermal velocity of 130 km s⁻¹ enters the Doppler dimming and the opening angle of the second generation emission. This corresponds to a kinetic temperature of log T = 6.0, however, the effective temperature including nonthermal (wave) motions could be larger. Frazin et al. (2003) report a 1/e width of Lyα in a streamer at 5.1 R_⊙ (their largest height) of about 155 km s⁻¹.

Another poorly known parameter is the outflow speed from the comet. We assume 10 km s⁻¹, though that is somewhat larger than is generally assumed. The speed of the H atoms produced by photodissociation of H₂O ranges from about 8 to 24 km s⁻¹, but some fraction of that speed is lost to interactions with O atoms. With smaller outflow speeds our simulation code does not resolve the region where first generation neutrals and PUIs coexist. A smaller outflow speed would reduce the first generation intensity by increasing the optical depth of the neutral cloud formed by dissociation of water. It would also increase the number of third generation neutrals because the PUIs would be formed in a region of higher neutral density.

Finally, there is some uncertainty in the orbit. This will mainly affect the projection of the cloud of Lyα onto the plane of the sky. In particular, a 10°–15° error in the ascending node could significantly change the projected angle between the second and third generation clouds.

Comparison of Figure 6 with Figures 1 and 2 shows that the models get one essential aspect right. The redshifted third generation emission lies to the south and the blueshifted second generation to the north, and their velocities agree with the observations. However, the models do not really resemble the data. The redshifted third generation emission appears as a short stub at a large angle to the main body of the comet tail in the model, while in the observations it is longer, and it makes a small enough angle to the main comet tail that it blends into the second generation emission. Because the predicted third generation emission is so spatially confined, it creates the strong peak for about 100 s when the comet crosses the slit, but the sharp peak is not resolved in the observations because of the 120 s integration time of each exposure.

Figure 9 compares the Lyα intensity predicted by the model for first, second, and third generation neutrals with the values observed by Giordano et al. (2015), and they match the observed intensities quite well. The model parameters are different than those derived by Giordano et al. (2015), who assumed a slower wind speed. The slower wind speed gives a good match to the opening angle of the Lyα tail, which then requires a higher density to match the falloff with time. It now seems likely that the opening angle of the tail is determined by the separation between second and third generation tails, and the second generation tail is narrower than we had thought, consistent with a faster solar wind. The faster wind speed also implies more severe Doppler dimming of the second generation emission, so our fit requires a higher outgassing rate of 10²⁹ s⁻¹, as opposed to 10²⁸ s⁻¹ from the earlier Monte Carlo simulation. The model in Figure 9 assumes a wind speed of 100 km s⁻¹, and a somewhat larger speed is needed to fully match the observed blueshift.

Overall, it appears that the basic idea of the third generation emission is correct, in that it explains the otherwise baffling red- and blueshifted emission from the south and north sides of the comet tail. Part of the discrepancy may be that the observations missed some of the first and third generation emission, because the comet was already in the UVCS slit during the first exposure. On the other hand, there is clearly something wrong with either the model parameters or with the basic physics in the model.

The parameter most closely related to the appearance of the third generation emission as a stub nearly at a right angle to the second generation emission is the angle between the comet trajectory and the magnetic field. When this angle is small, the PUIs move along the field at nearly the speed of the comet, which then means that they almost move with the comet and occupy a small interval relative to the comet before they are ionized. As described above, the angle is set by the comet trajectory and the assumption that the field is radial. The extremely elongated orbit and the small perihelion distance of the Kreutz sungrazers imply that the trajectory is not far from radial at 6 R_⊙, and the angle is 257. The magnetic field is expected to be close to radial in a steady model, but it could be perturbed by a CME. LASCO C2 difference images show a 450 km s⁻¹ west limb CME at 08:30 UT and continuous low level activity in the SW quadrant, so it is possible that the field in 10°–20° from radial. We note that the LASCO images are dominated by plasma near the plane of the sky, while the comet was about 4 R_⊙ away from the plane of the sky, so the LASCO movies may not accurately reflect conditions near the comet.

UVCS observed the comet at 8 R_⊙ from 21:31 until 22:22 UT on 2002 September 18. The UVCS images at 8 R_⊙ in Figures 1 and 2 show a split tail as well as red- and blueshifts on the southern and northern sides of the tail. The blueshifted northern side is very faint, which is the expected result of Doppler dimming of the second generation neutrals for a solar wind speed around 200 km s⁻¹. The split tail is also expected as a result of the angle between the second and third generation populations. However, the observed brightness increases gradually over the course of about 15 minutes (Figure 5 of Giordano et al. 2015), while the models predict a bright, short-lived third generation peak followed by a slowly fading second generation tail. The width of the tail during even the first UVCS exposure suggests that the observations may have missed the passage of the nucleus and first generation cloud across the slit, but that still would not explain the gradual brightening.

The models can easily match a third generation tail brighter than the second generation tail with a wind velocity high enough to give strong Doppler dimming. However, the models predict a large angle between the two tails and a third generation tail that is bright for only a short time. This is at least partly due to our approximate treatment of the third generation neutrals as discussed below.

There are three apparent aspects of the physics in the model that might be problematic; the assumption that the PUIs do not interact with the ambient medium, the assumption that we can ignore the effects of the comet on the ambient plasma and the assumption that we can treat each species as a single fluid when computing the Doppler dimming.

We have assumed that the PUIs form with parallel and perpendicular velocities determined by the angle between the comet trajectory and the B field, and that ${V}_{\parallel }$ is conserved while V_⊥ is scattered into a shell centered on ${V}_{\parallel }$. This assumption receives some support from the AIA observations of Comet Lovejoy (Raymond et al. 2014), though the AIA observation indicated that the shell in velocity space was quickly scattered into a Maxwellian. However, for the modest angle between the orbit and the radial B field at 6 R_⊙, ${V}_{\parallel }$ is more than twice ${V}_{\perp }$, so that the PUIs may be subject to violent streaming instabilities. Those instabilities can transfer momentum between the wind and the PUIs and they can heat both populations. They can also create large amplitude fluctuations in the magnetic field, so that a range of ${V}_{\parallel }$ and ${V}_{\perp }$ is present. Those effects will broaden the third generation tail, push it closer to the second generation tail, and make it longer. However, if the relative speed of the comet and the wind is less than the Alfvén speed the turbulence may be weak, while if the interaction is too strong it will reduce the redshift vanish, in contradiction to the observations.

We have also assumed that the mass flow from the comet does not greatly perturb the wind. That assumption is likely to be good for the relatively small $\dot{N}$ found by Giordano et al. (2015), but it is somewhat questionable for the higher value indicated by our model. The interaction can decelerate the wind, which would make it difficult to explain the strong blueshift observed. On the other hand, we do not independently know the wind speed, and we might be observing a flow that has been decelerated. If the perturbation is significant, a bow shock can form in the mass-loaded wind (Gombosi et al. 1996; Jia et al. 2014), and the changes in velocity and temperature can be large. Dramatic shifts of more than 100 km s⁻¹ in the velocity centroid and line width can be seen at outgassing rates above 10³⁰ s⁻¹ in Comet Lovejoy (Raymond et al. 2018).

The largest part of the discrepancy may be related to the perturbation of the wind by the comet. Behar et al. (2018a, 2018b) used analytic and numerical models to study the interaction between the outflowing gas of comet 67P/Churyumov–Gerasimenko and the solar wind. In that case, the Larmor radii of the ions are comparable to the size of the interaction region, and significant electric fields arise. In the case of C2/2002 S2, even the Larmor radius of the O⁺ ions is only a percent or so of the size of the interaction region, but the formation of third generation neutrals is strongly biased toward small distances from the comet because of the rapid density falloff. Behar et al. (2018b) showed that cometary ions are deflected in one direction and solar wind ions in the other, which would tend to produce the split tail we observe and separate the second and third generation neutrals. In the case of comet C/2002 S2, the main effect may be that the comet deflects the magnetic field. That both increases the angle between the field and the wind and deflects the wind away from its radial flow. The former increases the perpendicular component of the PUI velocity, reduces the parallel component, and produces a more elongated tail in the direction of the comet motion, producing a redshifted component more like that observed. The latter changes the angle between the solar wind and the LOS, decreasing the wind speed needed to match the observed LOS component. That helps with severe Doppler dimming expected from the observed blueshift if an unperturbed radial flow is assumed, and it therefore, reduces the required outgassing rate.

Finally, we have computed the Doppler dimming under the assumption that each component can be characterized by a Maxwellian velocity distribution or a shell in velocity space with a given width and centroid. This is a reasonable approach if there is no correlation between velocity and position, but it can break down for a noninteracting cloud of expanding neutrals. In particular, the third generation neutrals have a substantial velocity component toward the Sun and a somewhat narrower velocity width. Atoms that have relatively small velocities toward the Sun will stretch farther along the tail of the comet, and they will also experience relatively weak Doppler dimming. This coupling of the velocity and Lyα emission may explain the gradual intensity rise in the 8 R_⊙ observation, since the H atoms farther behind the comet scatter Lyα photons from the disk very effectively, while H atoms that remain close to the comet experience stronger Doppler dimming.