, , 1–8 ()

\begin{matrix}c \\ ◯ \end{matrix} Kluwer Academic Publishers, Boston. Manufactured in The Netherlands.

On unsteady oscillations developing in the stream effusing from a cometary nucleus due to a water molecule excitation ^{*}

ANDREY L. ITKIN itkin@momipao.msk.su International Institute for Problems of the Asteroid Hazard Apt.37, Taganskaya 44, Moscow, 109147, Russia

DOMENIQUE BOCKELEE-MORVAN

Observatory Medon, CNRS, France

Editor:

Abstract. We present the results of the modelling of 1-D unsteady water vapor outflow from a cometary surface, taking into account the interaction between solar radiation and the H_{2}O molecules. The vibrational and rotational molecular excitation model of [1], which includes infrared rovibrational pumping by solar radiation, thermal excitation by collisions and radiation trapping in the rotational and rovibrational lines, is used. The outflow of water vapor is described by time-dependent Euler equations with source terms representing the interaction between molecules and radiation. The coupled gasdynamic and relaxation kinetic equations are solved using factorization over the physical processes. The model is used to study the transient regime following a sudden onset of gas production from the nucleus surface.

Keywords: inner coma, water molecules, rotational excitation, simulation

1. Introduction

Radiative processes are expected to play an important role in the formation of the inner cometary comae, due to their strong effect on the gas temperature. The inner coma is characterized by a low gas temperature, due to the effect of expansion cooling. For H_{2}O molecules, the region of maximum solar radiation absorption lies in the IR. For these reasons, the vibrational degrees of freedom are not excited, but only the rotational and rovibrational ones. A physical model of the water molecule excitation by solar light has been developed in [1], [2] and confirmed by the observational results of [3]. Owing to the small number of populated rovibrational levels, and to the relatively high density of water molecules in the coma, a certain number of radiative lines are optically thick [1], [2]. While the number of significantly populated levels is not large, it is nonetheless necessary to consider a large total number of levels when modelling the radiative excitation. Accordingly, the

* This work was supported in part by the Russian Ministry of Science and Technical Policy, the Russian Foundation for Basic Researches, Grant N 97-03-46003 and by CNRS. The support of the French Embassy in Russia is also appreciated (A.L.I.). The authors thank Dr. J.F. Crifo for fruitful discussions.

use of exact radiative transfer equations would be a formidable task and approximate analytical methods have been developed to allow for optical depth effects. One of them is the escape probability method [4] originally proposed for a spherically symmetrical atmosphere. This method allows one to take into account the optical thickness effects by retaining the usual system of balance kinetic equations, and simply multiplying the Einstein coefficients by a proper factor. This factor is strongly nonlinear with respect to level populations.

Previous modelling of water vapor flow in the inner coma with allowance for radiative effects has been restricted to 1-D stationary conditions, in view of the large number of governing equations, and of their strongly non-linear nature [1]. In contrast, the present work is based on a time-dependent formalism, and allows the investigation of transient and/or unsteady effects in order to reveal how rotational nonequilibrium of the water molecules affects the gasdynamic parameters of the stream.

2. Gasdynamic model

The present approach is a rather classical way of treating flows with nonequilibrium physical and chemical processes. Water vapor issuing from the nucleus surface under solar illumination is assumed to be a mixture of ideal gases, each of which consists of water molecules in a certain internal state. In accordance with [1], we will here consider only the nonequilibrium excitation of the rovibrational sublevels which belong to the ground vibrational state of the molecule. We assume that the translational relaxation is fast enough that it is possib= le to introduce the total translational temperature of the mixture T, and we assume that the vibrational excitation of the water molecule is negligible.

2.1. Governing equations

For the sake of simplicity we assume the precedingly defined mixture to be an inviscid fluid – an acceptable assumption for the inner coma of active comets. We also assume that all components of the mixture have the same flow velocity and kinetic temperature. Let us introduce a spherical coordinate system (r,t) with the origin at the center of the comet nucleus. Under the preceding assumptions, the equations which govern the flow of the mixture can be written in conservative form [5]:

\frac{∂Ur^{2}}{∂t} + \frac{∂Er^{2}}{∂r} = ρr^{2}H (1)

where

U = | \begin{matrix}ρα_{i} \\ ρu \\ ρh_{0} − P \end{matrix} | , E = | \begin{matrix}ρuα_{i} \\ ρu^{2} + P \\ ρuh_{0} \end{matrix} | , H = | \begin{matrix}f_{m,i} \\ 2P/ρr \\ f_{e} \end{matrix} | (2)

Here t is the time, r the radial coordinate, u, ρ, P, h_{0} are the velocity, density, static pressure and total entalphy of the whole mixture, h_{0} = h + u^{2}/2, where h is the specific entalphy o= f the mixture, α_{i} the mass fraction of the i-th component of the mixture. Thus, ρα_{i} is the relative (mass) population of the i-th rovibrational sublevel of the ground vibrational state. Hereinafter for α_{i} we shall also use the term “population” because the connection between α_{i} and a true population of the corresponding level is transparent enough.

Equations (1) are, respectively, the budget equations for mass, momentum and energ. The subscript i labels the components of the mixture i = 1, N, where N is the total number of the components. The source term f_{m,i} describes the variation of population of the i-th sublevel due to collisional relaxation and radiative exchanges, and f_{e} is the variation of the total energy of the mixture due to radiation.

Summing the equations of the mass budget over all components yields the mass budget equation for the whole mixture, which using the normalization equations ∑_{i} α_{i} = 1 and ∑_{i} f_{m,i} = 0, takes the form

\frac{∂ρr^{2}}{∂t} + \frac{∂ρur^{2}}{∂r} = 0 (3)

In order to close system (1) two equations are necessary yet - the caloric and thermodynamic equations of state, which permit one to express the translational temperature of the mixture T and the specific entalphy h through the gasdynamic variables introduced above, and to find out an explicit form of the source terms. The thermodynamic equation of state is written in the usual form

P = ρRT ∑_{i} \frac{α_{i}}{μ_{i}} ≡ \frac{ρRT}{μ} ∑_{i} α_{i} = \frac{ρRT}{μ} (4)

where the above equalities hold because all μ_{i} are equal to one another, i.e. μ_{i} ≡ μ for any i, and the mass populations are normalized to unity. An expression of the specific entalphy of the mixture can be easily obtained from the specific internal energy of the mixture e, which can be introduced as follows

e = ∑_{i} α_{i}e_{i} = ∑_{i} α_{i} (a\frac{RT}{μ_{i}} + \frac{RE_{i}}{k_{B}μ_{i}} ) = \frac{RT}{μ} (a + \frac{1}{k_{B}T} ∑_{i} α_{i}E_{i}) (5)

Here e_{i} is the specific internal energy of the i-th component of the mixture, R is the universal gas constant, k_{B} is the Boltzmann constant, E_{i} is the energy of the i-th rovibrational sublevel counted from the energy of a zero level, a is the number of equilibrium degrees of freedom divided by two (in our case only the translational degrees of freedom are in equilibrium so a = 3/2). Hereinafter we shall use the value of E_{i} expressed in m^{−1}, designating this variable as E_{i}. Thus, E_{i} = h_{p}cE_{i}, where h_{p} is Planck’s constant, c is the speed of the light. Finally, taking into account, that h = E + P/ρ from (5) and (4) we obtain

h = \frac{RT}{μ} (a + 1 + \frac{ch_{p}}{k_{B}T} ∑_{i} α_{i}E_{i}) (6)

2.2. Boundary conditions

In this work we do not consider the effects connected with the existence of a Knudsen layer near the nucleus surface, therefore the boundary conditions are set directly on its surface. It is assumed that in the initial moment when the calculation is started the values of the total density, velocity and temperature of the flow are given as well as the values of N − 1 mass populations (since the population of the N-th component of the mixture can be found from the normalization condition). We suppos that the initial density, temperature and populations at the surface do not vary with time . As far as the gas velocity is concerned, here two situations are possible [6]. If the initial velocity of the flow at a certain time is greater than, or equal to the local speed of sound, then it can be preset at the surface of the nucleus. In the case of a subsonic velocity, however, it cannot be preset, since small perturbations can propagate downstream and thus influence the initial value of the velocity at the nucleus surface.

3. Kinetic model

To describe kinetics of radiative transitions model [1] is utilized. Within this model as the main processes affecting the energy redistribution between the components of the mixture and the external radiation we take into account infrared rovibrational pumping by the solar radiation flux, thermal excitation by collisions and radiation trapping in the rotational and rovibrational lines.

3.1. Kinetic equations

The system of kinetic equations describing the above mentioned processes is written as follows

\frac{dα_{i}}{dt} = f_{m,i} = Ω_{c,i} + Ω_{r,i},  Ω_{c,i} = C_{c} (α_{ie} − α_{i})  (7)
Ω_{r,i} = ∑_{i′} A_{i′i}α_{i′} − α_{i} ∑_{i′} A_{ii′} + (∑_{i′} B_{i′i}J(λ_{i′i})α_{i′} − α_{i} ∑_{i′} B_{ii′}J(λ_{ii′}))

In (7) Ω_{c,i} and Ω_{r,i} are the collisional and radiative fluxes for the i-th component of the mixture, A_{ii”} and B_{ii”} are the Einstein coefficients for spontaneous and induced radiation, λ_{ii”} is the wavelength of the transition, J is the intensity of the radiative flux averaged over all the directions and integrated over the whole line, C_{c} is the total collision rate for rotational transitions, α_{ie} is the equilibrium populations at the given temperature. When writing (7) a model of the strong collisions is used allowing to write the collisional term in the BGK form. The coefficients A_{ij} are equal to zero if E_{i} < E_{j}.

In accordance with the accepted form of the kinetic equations the source term f_{e} can be represented as follows: f_{e} = ∑_{i} E_{i}Ω_{r,i}.

3.2. Constants and parameters of the model

To close the system (7), it is necessary to know the Einstein coefficients, the intensity of radiation and the parameters present in the collisional term. In accordance with [1] we shall consider the system of 32 ortho and 32 para rotational sublevels of the ground vibrational state of the water molecule. In view of the existence of selection rules, the initial system (7) of 64 ordinary differential equations is factorized into two independent systems of 32 equations.

In the equilibrium all the ortho and para states are distributed with statistical weights 0.75 and 0.25, respectively. Thus, the equilibrium populations can be determined from the Boltzmann distribution, having the form

α_{ie}^{ν} = \frac{p^{ν}}{Q^{ν}} w_{i}^{ν}exp (−\frac{E_{i}^{ν}}{k_{B}T} ) ,  Q^{ν} = p^{ν} ∑_{i=0}^{N/2} w_{i}^{ν}exp (−\frac{E_{i}^{ν}}{k_{B}T} )  (8)

where index ν accepts the values ”” and ”p” for ortho and para states, respectively; p^{o} = 0.75, p^{p} = 0.25.

For the total collisional rate C_{c} the usual expression of the kinetic theory of gases [7] is used

C_{c} = \frac{ρσN_{a}}{μ} \sqrt{\frac{16RT}{πμ}} (9)

where the collisional cross-section of the water molecule is assumed to be constant σ = 5 · 10^{−14} cm^{−2}, N_{a} is the Avogadro constant.

To compute the mean intensity of radiation, the escape probability approximation is used [4], whereby the profile-averaged mean radiation intensity at frequency ν_{ij} of the transition i → j is given by

J_{ij} = (1 − β_{ij})S_{ij} + β_{ij}^{c}I(ν_{i}j) (10)

where S_{ij} is the source function, β_{ij} is the average over 4π solid angle of β_{ij}(n) which is the probability that an emitted photon escapes the line trapping region in direction n without suffering any absorbtion, β_{ij}^{c}I(ν_{ij}) describes the amount of external radiation I(ν_{ij}) - here the solar field - participating to the excitation, β_{ij}^{c} is β_{ij}(n) averaged over the solid angle of the excitation source. With the assumption of complete frequency distribution we have J_{ij} = J_{ji}. It is readily shown that including radiation trapping effects - i.e. the first term of (10) into the rate equations is equivalent to multiply the Einstein A_{ij} coefficient by the photon escape probability β_{ij}. One can easily show that the opacity effects will reduce the optically thin solar excitation B_{ij} by a factor β_{ij}^{c}.

The calculation of β_{ij} and β_{ij}^{c} in case of spherical geometry and for positive optical depth was made in [8] and gives

β_{ij} = \frac{2}{3τ_{ij}} − \{\frac{τ_{ij}}{3} [K_{2}(τ_{ij}/2) − K_{1}(τ_{ij}/2)] − K_{1}(τ_{ij}/2)\} exp(−τ_{ij}/2)
τ_{ij} = \frac{A_{ji}c^{3}g_{j}}{8πν_{ij}^{3}} [\frac{α_{i}}{g_{i}} − \frac{α_{j}}{g_{j}} ] \frac{ρ(r)Rr}{μu}  (11)

Here τ_{ij} is the Sobolev optical depth, K_{2} and K_{1} are modified Bessel functions of the second kind, g_{i} and g_{j} are the degeneracies in the lower and upper levels of the transition. Thus, right part (7) is essentially nonlinear function of populations α.

It is necessary to emphasize that expression (10), derived in [8], is rigorous, only if the flow velocity u does not depend upon r. If u ∼ r, (10) is transformed to the following expression

β_{ij} = \frac{1 − exp(−τ_{ij})}{τ_{ij}} (12)

Values of β_{ij} obtained from (10) and (12) coincide at τ_{ij} = 0, but at large values of τ_{ij} their ration is equal to 2/3. In case of an arbitrary dependence of u upon r the values of β_{ij} depend upon r. However, it is possible to expect that an error in determination of β_{ij} due to disregard of the real dependence of u upon r does not significantly influence, at least, the qualitative character of the populations’ evolution.

At calculation of the phonon escape probability β_{ij}^{c} we assume it to be related to β_{ij} through the dilute factor

β_{ij}^{c} = β_{ij}/r_{h}^{2} (13)

where r_{h} is the heliocentric distance of the comet. Generally speaking, this assumption is strictly valid only if u ∼ r.

For negative optical depth β_{ij} is no more a probability but a maser amplification factor. Since population inversions may occur in water molecules for rotational transitions - notably for 6_{16} − 5_{23} transition - the opportunity of massing must be included in the model. We will assume that the dependence of the amplification factor on optical depth behaves like (12) because the exponential behaviour of (12) ensures a rapid increase of β_{ij} with negative optical depths and appears realistic for including maser amplification.

Energies of orto and para levels, considered in the model, and Einstein coefficients A_{i′i} and B_{i′i}J(λ_{i′i}) for optically thin lines are given in [9]

4. Numerical method

To solve numerically the system of equations (1), (3), (4), (6) a finite-difference method is used which has been originally proposed in [10] and being a modification of a widely known Godunov’s method. The method is besed upon a piece–linear approximation of flow parameters in each cross-section t = const and possesses a second order of accuracy on smooth solutions with preservation of a proper monotony. This method was specially developed to handle flows with strong discontinuities (for instance, shock waves) in a general way, i.e., without the need to prescribe a priori the localization of these discontinuities. From other similar methods it is distingvished by use of a predictor-corrector procedure that makes the algorithm the rather simple and universal one. In addition, a method of factorization over physical processes is used, consisting in the two steps - gasdynamic and kinetic one. Details of the scheme are given in [9].

5. Results and discussion

For the numerical calculations we considered a weak comet with nucleus radius R_{n} = 6 km and located at the distance 2 AU from the Sun. In the first series of calculations the surface temperature was taken equal to T_{s} = 140K and the vapor mass density on the external boundary of the Knudsen layer adjacent to the surface was taken to the saturation value ρ_{s} = 6 · 10^{−14}g/cm^{3}. The initial velocity of the vapo= r was determined in accordance with the special algorithm described in section 2.2. Finally, the populations of rotational sublevels were assumeed to be in local equilibrium at the given temperature and at all times.

At the initial time the following values of parameters were set at all radial grid points outside of the nucleus surface: ρ_{i} = ρ_{s}/100, P_{i} = P_{s}/100, T_{i} = T_{s}, (α_{ik} = α_{sk} = α_{ek}, k = 1...64) =, where i labels the radial points, s the values at the surface, and e the equilibrium values. The low initial values of the density and temperature in the coma in comparison with the values of at the surface were used for modeling a possible sudden increase in surface sublimation.

Our simulation shows (see [9]) that in such a case nonequilibrium processes take place, induced by the interaction of water molecules with the radiative flux from the Sun. Thus, the populations of some levels become smaller than their equilibrium values, but for other ones the inverse situation is realized. The reason is that for the present conditions the radiative processes prevail over the collisional relaxation, and the steady distribution of populations is largely determined by the external radiation. It turns out that with time the absolute values of the populations do not practically change at the bottom levels, while at the top levels they decrease. An existence of inversion for certain levels in the population distribution could also been noted that does not result in the generation of radiation from them, because the computed factor of amplification everywhere is less than zero.

When the picture of the flow is close to the stationary one the values of populations for some levels smoothly fall down with the growth of the distance from nucleus, for the other levels — on the contrary, increase, and for the third ones — have a wavy character.

In the second series of calculations all the initial and boundary conditions were chosen as in the previous case, except that the temperature of the surface was taken equal to T_{s} = 170K and the heliocentric distance to 1AU. The corresponding density of the saturated water vapor at the given temperature is ρ_{s} = 5.4 · 10^{−11}g/cm^{3}, whereby Q ≈ 6.8 · 10^{28} molecule s^{−1}.

An increase of the initial vapor density results in more frequent collisions, which reduce the time of the collisional relaxation τ so, that it becomes less the the specific gasdynamic time. As the result in each computational point over the space coordinate and at all times the equilibrium distribution over the energy of rotational sublevels has time to be established.

On the other hand, an intensive exchange by energy between the internal and translational degrees of freedom results in a sharp difference of the flow parameters from the adiabatic ones. As shown in [9] during the flow evolution non-stationary oscillations of the temperature are observed in the stream. To be convinced that these oscillations have physical nature and are not connected with the possible oscillations of the finite-difference sheme of the second order, the similar calculations were executed with the help of the sheme having first order of accuracy, which revealed the same qualitative results. However, for infinitely large times a stationary regime of an adiabatic flow with exponent equal to 1.29 should be obtained. Thus, the preceding oscillations should vanish with time.

We shall notice that the heliocentric distance of the comet r_{h} appears only in formula (13). Thus, a variation of r_{h} affects only the rate of the radiative processes. At high vapor initial temperatures and densities where at each time an equilibrium distribution of populations has time to be established, such a variation of the kinetic coefficients hardly influences the parameters of the flow. Thus, under such conditions it is possible to use the present results to discuss the coma behaviour at differing heliocentric distances.

References

1. Bockelee-Morvan, D., ”A model for the excitation of water in comets,” Astron.Astrophys., 181, 1987, pp.169–181.

2. Crovisier, J., ”The water molecule in comets: Fluorescence mechanisms and thermodynamics of the inner coma,” Astron.Astrophys., 130, 1984, pp.361–372.

3. Weaver, H.A. and Mumma, M.J., ”Infrared investigations of water in comet P/Halley,” Astron.Astrophys., 187, 1987, pp.411–418.

4. Rybicky, G.B., In W. Kalkofen, editor, ”Methods in Radiative Transfer”, Cambridge University Press, Cambridge, 1984, p.21.

5. Sternin, L.E., ”Foundations of gasdynamics of two-phase flows in nozzles”, Mashinostroenie, Moscow, 1974. in Russian.

6. Crifo, J.F., Itkin, A.L., and Rodionov, A.V., ”The near-nucleus coma formed by interacting dusty gas effusing from a cometary nucleus. I.,” Icarus, 116, 1995, pp.77–112.

7. Ferziger, J.H., and Kaper, H.G., Mathematical theory of transport processes in Gases. North-Holland, Amsterdam, 1972.

8. Litvak, M.M., and Rodrigues Kuiper, E.N., ”Cometary NH: Ultraviolet and submillimeter emission,” Astrophys.J., 253, 1982, pp.622–633.

9. Itkin, A.L., Bockelee-Morvan, D., and Rodionov A.V., Space Researches, 35(1), 1997, pp.18– 36.

10. Rodionov, A.V., ”Methods of increasing the accuracy of godunov’s scheme,” USSR Comput. Maths. Math. Phys., 27(6), 1987, pp.164–169.