A NEW 3D MULTI-FLUID MODEL: A STUDY OF KINETIC EFFECTS AND VARIATIONS OF PHYSICAL CONDITIONS IN THE COMETARY COMA

Our multi-fluid model solves a set of hydrodynamic equations for each species, which are shown in Equations (1)–(3). They are the continuity equation, the momentum equation, and the pressure equation. ρ_s, ${{\boldsymbol{u}}}_{s}$, p_s are the mass density, the velocity vector, and the scalar pressure of the neutral species s. γ_s denotes the specific heat ratio of species s. The terms on the right-hand side of the equations represent the source terms, which are specified for each species and include collisions, chemical reactions, and the excess energies accompanied by photo-dissociations.

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{s}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{s}{{\boldsymbol{u}}}_{s})=\displaystyle \frac{\delta {\rho }_{s}}{\delta t}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{s}{{\boldsymbol{u}}}_{s}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{s}{{\boldsymbol{u}}}_{s}{{\boldsymbol{u}}}_{s}+{p}_{s}{ \mathcal I })=\displaystyle \frac{\delta {\rho }_{s}{{\boldsymbol{u}}}_{s}}{\delta t}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {p}_{s}}{\partial t}+{\rm{\nabla }}\cdot ({p}_{s}{{\boldsymbol{u}}}_{s})+({\gamma }_{s}-1){p}_{s}({\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{s})=\displaystyle \frac{\delta {p}_{s}}{\delta t}.\end{eqnarray} \tag{ 3 }$$

The source terms are described in Equations (4)–(7). Those source terms with chemical reaction frequencies ${\nu }_{s\to t}$ are related to photochemical reactions. The photochemical reactions, reaction rates, and the corresponding excess energies used in the model are shown in Table 1. In the pressure source term, the excess energies are partitioned under the restriction that the momenta of the daughter species should be conserved and thus are inversely proportional to the mass. The source terms involving momentum transfer coefficients ${\overline{v}}_{s,t}$ between species s and t are terms accounting for collisions. The collision frequency is linearly proportional to the density of each involved gas species, and the relative speed of the colliding gas molecules or atoms, which is calculated from the thermal and bulk velocities of both species. The cross sections for the modeled collisions are also listed in Table 2. We note here that cross sections of self-collisions are not included, since fluid approaches assume an approximately Maxwellian distribution for each fluid, implying plenty of self-collisions in the gas of the species. To account for the infrared cooling effect of H₂O, a cooling function L is added to the pressure source term of H₂O. L is expressed as $L=\tfrac{8.5\times {10}^{-19}{T}_{{{\rm{H}}}_{2}{\rm{O}}}^{2}{n}_{{{\rm{H}}}_{2}{\rm{O}}}^{2}}{{n}_{{{\rm{H}}}_{2}{\rm{O}}}+2.7\times {10}^{7}{T}_{{{\rm{H}}}_{2}{\rm{O}}}}$ erg cm⁻³ s⁻¹, which was first proposed by Shimizu (1976), and works well for the cases in this paper. But it needs additional adjustment involving optical depth for production rates higher than 10³⁰ s⁻¹, since the inner coma under such circumstances is optically thick, preventing efficient radiative cooling. The adjusted cooling rate can be found in Gombosi et al. (1986), which does not make much difference for the cases shown here compared with more complex approaches that account for the non-LTE effects between kinetic and rotational temperature (Combi 1996; Tenishev et al. 2008).

$$\begin{eqnarray}&&\displaystyle \frac{\delta {\rho }_{s}}{\delta t}=\displaystyle \sum _{n={\rm{neutrals}}}{\nu }_{n\to s}{n}_{n}{m}_{s}-\displaystyle \sum _{t={\rm{other}}\,{\rm{species}}}{\nu }_{s\to t}{n}_{s}{m}_{s}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\displaystyle \frac{\delta {\rho }_{s}{{\boldsymbol{u}}}_{s}}{\delta t} & = & \displaystyle \sum _{n={\rm{neutrals}}}{\nu }_{n\to s}{n}_{n}({{\boldsymbol{u}}}_{n}-{{\boldsymbol{u}}}_{s}){m}_{s}+\displaystyle \frac{\delta {\rho }_{s}}{\delta t}{{\boldsymbol{u}}}_{s}\\ & & +\displaystyle \sum _{t={\rm{other}}\,{\rm{species}}}{\overline{\nu }}_{s,t}{n}_{s}{m}_{s}({{\boldsymbol{u}}}_{t}-{{\boldsymbol{u}}}_{s})\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\overline{\nu }}_{s,t}=\displaystyle \frac{{m}_{t}}{{m}_{t}+{m}_{s}}{n}_{t}{\sigma }_{s,t}\sqrt{\displaystyle \frac{8k}{\pi }\left(\displaystyle \frac{{T}_{t}}{{m}_{t}}+\displaystyle \frac{{T}_{s}}{{m}_{s}}\right)+{({{\boldsymbol{u}}}_{t}-{{\boldsymbol{u}}}_{n})}^{2}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}\displaystyle \frac{\delta {p}_{s}}{\delta t} & = & ({\gamma }_{s}-1)\displaystyle \frac{{m}_{s}}{{m}_{n}}{\nu }_{n\to s}{n}_{n}{\rm{\Delta }}{E}_{n\to s}\\ & & +\displaystyle \frac{1}{2}{\nu }_{n\to s}{p}_{n}+\displaystyle \frac{1}{3}{\nu }_{n\to s}{n}_{n}{({{\boldsymbol{u}}}_{s}-{{\boldsymbol{u}}}_{t})}^{2}{m}_{s}\\ & & -\displaystyle \sum _{t={\rm{other}}\,{\rm{species}}}{\nu }_{s\to t}{p}_{s}+\displaystyle \sum _{t={\rm{other}}\,{\rm{species}}}{\overline{\nu }}_{s,t}\displaystyle \frac{{m}_{t}{m}_{s}}{{m}_{t}+{m}_{s}}{n}_{s}\\ & & \times \left[\displaystyle \frac{2}{3}{({{\boldsymbol{u}}}_{s}-{{\boldsymbol{u}}}_{t})}^{2}+\displaystyle \frac{2k({T}_{t}-{T}_{s})}{{m}_{t}}\right].\end{eqnarray} \tag{ 7 }$$

We note here that collisions in our model conserve momentum and energy. The sums of the collisional sources of momentum and energy of two colliding fluids are zero. However, the pressure sources rather than the energy sources are presented here, because the collisional terms in pressure equations are more concise and clear.

Table 1.  Chemical Reactions and Corresponding Parameters

Wavelength (Å)	Reaction	Reaction Rate (s−1)	Excess Energy (eV)
1357–1860	H2O + hν $\to$ H + OH	8.0 × 10−6	1.7
	H2O + hν $\to$ H2 + O	8.4 × 10−8	1.7
1216	H2O + hν $\to$ H + OH	8.0 × 10−6	4.5
	H2O + hν $\to$ H2 + O	2.8 × 10−7	1.7
984–1357 (excluding 1216)	H2O + hν $\to$ H + OH	3.6 × 10−7	4.5
	H2O + hν $\to$ H2 + O	4.8 × 10−8	1.7
<984	H2O + hν $\to$ ionization products	7.0 × 10−7
2160	OH + hν $\to$ H + O	4.5 × 10−6	0.36
2450		5.0 × 10−7	0.67
1400–1800		1.4 × 10−6	3.2
1216	OH (12Δ)	3.0 × 10−7	3.8
	OH (${{\rm{B}}}^{2}{{\rm{\Sigma }}}^{+}$)	5.0 × 10−8	1.6
	OH (2 ${}^{2}{\rm{\Pi }}-{3}^{2}{\rm{\Pi }}$)	5.0 × 10−8	3.8
<1200	OH (${{\rm{D}}}^{2}{{\rm{\Sigma }}}^{-}$)	1.0 × 10−8	2.7
	H2 + hν $\to$ H + H	1.1 × 10−7	1.8

Download table as:  ASCIITypeset image

Our model is based on the BATS-R-US (Block-Adaptive Tree Solar wind Roe-type Upwind Scheme) code (Powell et al. 1999; Tóth et al. 2012), which is capable of solving the magnetohydrodynamics (MHD) and hydrodynamics equations efficiently on adaptive grids. The simulations are performed on a 3D spherical grid, the radius of which is 1.4 × 10⁶ km for the comet 67P case. A spherical body with a radius of 2 km is placed at the origin to model the comet nucleus. We note here that a realistic nucleus radius of comets with a production rate higher than 10³⁰ s⁻¹ should be more than 20 km. However, as our model shows for comets with high production rates, a realistic radius does not change the results at cometocentric distances larger than 50 km. The resolution in the radial direction is about 370 m near the nucleus and about 1.2 × 10⁵ km near the outer boundary of the computing domain. The resolutions in the polar and azimuthal directions are about 07 and 14. This grid allows studying the coma in different length scales but without too heavy of a computational burden. In addition to the adaptive grid, two features in the BATS-R-US code also greatly improve the efficiency and accuracy of the model. The first is the steady-state mode, where the time steps are different in every grid cell limited by the local stability condition only, so that the time to the convergence is reduced. The second is the point-implicit scheme (Tóth et al. 2012), which facilitates the calculation of the stiff source terms, which are mainly related to the photo-chemistry in the coma without spatial derivatives involved. For example, during a single computational time step, one minor species may have a density so tiny that it may be comparable to the incremental density added by the photochemical reactions. An explicit scheme cannot compute such terms efficiently and accurately, but the point-implicit scheme can handle them well.

To compare our results with the DSMC results (Tenishev et al. 2008), we run the model to simulate comet 67P at four heliocentric distances with the same or close set of parameters as possible between fluid and kinetic models. The inner boundary condition on the surface of the nucleus for the only parent species, water, is fixed. H₂O flux F and H₂O temperature T are set as functions of the solar zenith angle, which is described by Tenishev et al. (2008). Following Huebner & Markiewicz (2000) and Bieler et al. (2015), the magnitude of the H₂O velocity u on the surface is $0.8257\sqrt{8{kT}/(\pi m)}$, where k is the Boltzmann constant and m is the molecular mass of H₂O. The velocity is normal to the surface. The number density is calculated by $F/| u|$. The pressure at the boundary is set to be $(F/u)k(0.9049T)$. For all other species, the photo-dissociation products of water, a floating boundary condition is set for all variables, in which a zero gradient is imposed. At the outer boundary, floating boundary conditions are also applied for all variables.

In Section 3.2, we study the effects of production rates on the coma morphology, so the inner boundary conditions are slightly modified. The flux is set to be uniform on the surface and the temperature is fixed to 180 K. The heliocentric distances are fixed to 1.0 au for all cases and we only vary the neutral gas production rate from 10²⁷ to 10³⁰ s⁻¹ in these simulations.

To compare with remote observations of several comets, we run more cases with varying heliocentric distances and different production rates, and obtain the expansion speed of H₂O at about 10⁵ km from the nucleus. The distance is chosen mainly because the expansion speeds reported by Tseng et al. (2007) were measured approximately at that distance.

2D cuts of densities, speeds, and temperatures of H₂O and H from the multi-fluid model results at 1.3 au are shown in Figure 1. The Sun is in the direction of negative x-axis. The effect of solar illumination can be seen from the H₂O results. Density, speed, and temperature of H₂O are higher on the dayside than the nightside. But this is not true for H. Because of more collisions taking place on the dayside, the dayside speed and temperature of H are lower. A similar picture to H applies to other daughter species.

Zoom In Zoom Out Reset image size

In this section, we will juxtapose our model results and the DSMC solution, and display the similarities and the differences between them. Specifically, the 1D profiles of velocities, temperatures, and densities of modeled species extracted along the comet-Sun line are compared. Such comparison are made at four heliocentric distances: 1.3, 2.0, 2.7, and 3.3 au with production rates of 5 × 10²⁷, 8 × 10²⁶, 8 × 10²⁵, and 1 × 10²⁴ s⁻¹, respectively.

3.1.1. Velocity

Figure 2 shows the speed of each species at four heliocentric distances, with the left column displaying our fluid model results and the right column the DSMC results from Tenishev et al. (2008). The following figures of temperatures and densities have the same format. We can spot three groups of lines in the four cases in both columns. H and H₂ behave as one group, while O and OH are another one. Each group has similar masses, and gains energy from the photo-dissociation. The group of H and H₂ has the highest speeds, the group of O and OH has the second highest. H₂O almost stays level after a short distance of acceleration, with collisions with daughter species as its only source for acceleration. The speeds of H and H₂ are decoupled from H₂O before O and OH diverge from H₂O at a larger distance. The decoupling distance decreases with lower production rates. This trend can be readily explained by fewer collisions in a thinner coma and more excess energy translated to the kinetic energy. Therefore, as production rate decreases, H and H₂ reach the plateau of the bulk speed at about 5 km s⁻¹ in a shorter distance, resulting in a steeper slope in the velocity profile.The new fluid and DSMC kinetic results are quite similar in all respects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.2. Temperature

Figure 3 displays temperatures for all species at the selected four heliocentric distances. The temperature behavior can also be divided into the three groups mentioned in the previous section. H and H₂ have the highest temperatures. OH and O are intermediate. Without any excess energy input, H₂O is the coldest. The temperature of H₂O decreases due to the adiabatic expansion at distances less than 100 km, then remains at the level of around 10 K. For cases with larger production rates, the temperatures of H, H₂, O, and OH first drop slightly before increasing to high levels. As the drops are mainly caused by the collisions with the cold H₂O, the lack of collisions for the cases with low production rates leads to the disappearance of the dips.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also notice that temperatures of all daughter species decline at large distances. This may be for two possible reasons. The first is the cooling effect caused by adiabatic expansion. The second reason is that the relative abundance of the parent species to the daughter species decreases farther out. As a result, the percentage of newly born daughter species with a high photo-dissociation temperature drops in the population of the daughter species. The bulk temperature is thus decreased. Later we will show the relative abundance of the parent species to the daughter species also decreases faster at smaller heliocentric distances due to the shorter photo-lifetime of the parent species, which explains why in the 1.3 au case the temperatures of daughter species drop most quickly.

The most obvious difference between our fluid model and the DSMC model is the decreasing rate of the water temperature. In the fluid model results, the water temperature decreases to the minimum temperature of 10 K within 20 km for all four cases. But the slopes in the DSMC results are flatter, especially in the cases of 1.3 and 2.0 au with the minimum temperature reached at 100 km. This may be caused by the fluid model overestimating the self-collisions, so adiabatic expansion happens more efficiently than in the DSMC which yields less frequent collisions.

It is interesting to see in the DSMC results that the H₂O temperature jumps to about 100 K near 10⁶ km in the 1.3 au case. According to Tenishev et al. (2008), this is caused by the selection effect. Most of the slower H₂O particles have already been destructed by photo-dissociation. We can also see that in the velocity plot, close to 10⁶ km the water velocity goes up significantly. Due to the nature of the fluid approximation, our model is not able to include such purely kinetic effects related to distribution functions.

3.1.3. Density

Figure 4 presents the densities of all species versus the cometocentric distances. The densities decrease almost linearly on log–log plots. The two columns look similar. H₂O decreases the fastest, since it gets photo-dissociated at a high rate without any new supply. In the cases at 1.3 and 2.0 au, we can see H and OH have the same density near the nucleus and diverge after some distance. The explanation is that they are produced by the same chemical reaction and they share the same velocity in the collisional region. Outside the collisional region, freshly produced fast H dominates compared with slower OH. As a result, the H density declines faster. In other cases, while even in the vicinity of the nucleus is collisionless, the H density is lower than OH all the way out. The same reasoning also applies to H₂ and O.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To study the effects of the production rate on coma dynamics, we run the model at 1.0 au but with four different production rates: 10²⁷, 10²⁸, 10²⁹ and 10³⁰ s⁻¹. Figure 5 shows, for all four cases, the speed of water, the mean speed of the three heavy species (H₂O, OH, O), and the mean speed of all five neutral species. The 10²⁷ and 10²⁸ s⁻¹ cases are very similar in the speed profile, suggesting the role of collisional heating is negligible with a low production rate or a low density. The water terminal speed in the 10²⁷ and 10²⁸ s⁻¹ cases is 0.7 km s⁻¹, which is mainly determined by the initial temperature. As the production rate increases, the collisional heating effect kicks in. In the 10²⁹ s⁻¹ case, the water speed at 10⁵ km is 0.9 km s⁻¹ and it rises to about 1.4 km s⁻¹ in the 10³⁰ s⁻¹ case. This result also agrees well with Bockelee-Morvan & Crovisier (1987). In the other two plots, we find the mean speed accelerates beyond 10⁴ km because of the contributions by the fast species. As a result, the increase of the production rate has more impact on water speed than the bulk speeds of heavy species and all combined species.

Zoom In Zoom Out Reset image size

Figure 6 shows the temperature of H₂O, the mean temperatures of three heavy species (H₂O, OH, O), and the mean temperatures of all five neutral species for the four cases. In the H₂O temperature profile, there is a peak beyond 10⁴ km in the cases of 10²⁹ and 10³⁰ s⁻¹ production rates. The general trend is similar to that in Bockelee-Morvan & Crovisier (1987). But the specific peak values are slightly different between our model and theirs. In the 10³⁰ s⁻¹ case, our peak near 140 K is higher than the peak at 110 K in their model. It may be because their model is single-fluid and treats photochemical heating and radiative cooling in a different way to ours. In the 10²⁹ s⁻¹ case, both models have roughly the same peak near 45 K. The temperatures in the other two plots all have dips within 1000 km but then the uptrend remains beyond 1000 km. Though the temperatures of the four cases increase at different rates, they are close to each other around 10⁵ km. This suggests that, when the cometocentric distance is large enough, the average temperature of all gases may not vary much with the production rate. At large cometocentric distances all species become collisionally decoupled and the secondary species trend to the photo-production temperatures, which are independent of the production rates. Figures 5 and 6 illustrate the radiative cooling effect on the speed and the temperature of water. Within 1000 km, in the 10³⁰ s⁻¹ case, because of the most significant cooling effect caused by the highest water density, the speed and temperature are dampened for a larger distance than those in other cases.

Zoom In Zoom Out Reset image size

Figure 7 shows the mean molecular mass versus the cometocentric distance for varying production rates. Since the comet in this study is kept at a heliocentric distance of 1.0 au and has H₂O as its only parent species, the four cases yield very similar curves. The difference in the speed profiles is the only factor capable of altering the relative abundances and thus the mean molecular mass. We can expect, if all species have fixed speed profiles for each individual species, that the four cases should have the same mean molecular mass profile. We also notice that the mean molecular masses at 10⁵ km almost converge to a single value of about 15.4 amu.

Zoom In Zoom Out Reset image size

In this section, we characterize the cometary H₂O expansion speeds extracted from our multi-fluid gas coma model at a cometocentric distance of 10⁵ km. The model is run with several selected production rates and heliocentric distances. The results are listed in Table 3. One can readily see from the table that a larger production rate and a smaller heliocentric distance (i.e., a higher photon flux) lead to a higher expansion speed as expected (Combi 1987).

Tseng et al. (2007) derived the H₂O expansion speeds from the 18 cm line shapes of the OH radicals observed by radio telescopes in over 30 comets, which is reproduced in Figure 8 and serves as benchmark for our model. Our model results, which are denoted by solid diamonds, are superimposed on the observations. We note here the x-axis in Figure 8 represents the OH production rate, which is often obtained by multiplying a factor of 0.86, the photo-dissociation branching ratio of H₂O to OH, with the H₂O production rate (Harris et al. 2002; Combi et al. 2004, pp. 523–552). For simplicity, we assume the OH and H₂O production rates are the same. Beyond 0.6 au, the diamonds in the 10²⁸ and 10²⁹ s⁻¹ cases are close to each other, while the spread between the 10²⁹ and 10³⁰ s⁻¹ cases is significantly larger. In the figure where the heliocentric distance is smaller than 0.6 au, the increase in production rate results in a more evenly increase in water speeds than at other distances. This reflects the nonlinear effect of the production rate and the heliocentric distance on the expansion speed, which is similar to the threshold effect mentioned by Tseng et al. (2007). In addition, the model is also applied to comet C/1995 O1 (Hale-Bopp) at the heliocentric distance of 1.0 au. The radiative cooling effect is neglected within 10⁴ km because of the large high production rate of 8 × 10³⁰ s⁻¹. Our model yields similar results to that of the single fluid model in Combi et al. (1999) and Combi (2002), which matched the observations in Biver et al. (2002). It is because the production rate is so large that most heavy species are coupled within a cometocentric distance of 10⁵ km and thus the single fluid assumption is still valid. In any case, our multi-fluid model well reproduces the observed variation in coma outflow speeds with different comet gas production rates and heliocentric distances.

Zoom In Zoom Out Reset image size