Gas kinetics and dust dynamics in low-density comet comae

Icarus 210 (2010) 455–471

Contents lists available at ScienceDirect

Icarus journal homepage: www.elsevier.com/locate/icarus

Gas kinetics and dust dynamics in low-density comet comae Björn J.R. Davidsson a,*, Samuel Gulkis b, Claudia Alexander b, Paul von Allmen b, Lucas Kamp b, Seungwon Lee b, Johan Warell a a b

Department of Physics and Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

a r t i c l e

i n f o

Article history: Received 30 October 2009 Revised 18 June 2010 Accepted 18 June 2010 Available online 25 June 2010 Keywords: Comets, Coma Comets, Dust Spectroscopy

a b s t r a c t Extensive regions of low-density cometary comae are characterized by important deviations from the Maxwell–Boltzmann velocity distribution, i.e. breakdown of thermodynamic equilibrium. The consequences of this on the shapes of emission and absorption lines, and for the acceleration of solid bodies due to gas drag, have rarely been investigated. These problems are studied here to aid in the development of future coma models, and in preparation for observations of Comet 67P/Churyumov–Gerasimenko from the ESA Rosetta spacecraft. Two topics in particular, related to Rosetta, are preparation for in situ observations of water, carbon monoxide, ammonia, and methanol emission lines by the mm/sub-mm spectrometer MIRO, as well as gas drag forces on dust grains and on the Rosetta spacecraft itself. Direct Simulation Monte Carlo (DSMC) modeling of H2O/CO mixtures in spherically symmetric geometries at various heliocentric distances are used to study the evolution of the (generally non-Maxwellian) velocity distribution function throughout the coma. Such distribution functions are then used to calculate Doppler broadening profiles and drag forces. It is found that deviation from thermodynamic equilibrium indeed is commonplace, and already at 2.5 AU from the Sun the entire comet coma displays manifestations of such breakdown, e.g., non-equal partitioning of energy between kinetic and rotational modes, causing substantial differences between translational and rotational temperatures. We exemplify how deviations from thermodynamic equilibrium affect the properties of Doppler broadened line profiles. Upper limits on the size of liftable dust grains as well as terminal grain velocities are presented. Furthermore, it is demonstrated that the drag-to-gravity force ratio is likely to decrease with decreasing cometocentric distance, which may be of relevance both for Rosetta and for the lander probe Philae. Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction The European Space Agency cornerstone mission Rosetta was launched on March 2, 2004 and is currently en route to its primary target, Jupiter Family Comet 67P/Churyumov–Gerasimenko. Initiation of near-nucleus operations is scheduled for August 22, 2014, at a heliocentric distance of rh = 3.25 AU. Rosetta carries a number of in situ and remote-sensing instruments for characterization of the gas and dust properties of the comet coma (see, e.g., Glassmeier et al., 2007), including the millimeter–submillimeter spectrometer MIRO, built and operated by Jet Propulsion Laboratory (Gulkis et al., 2009). The MIRO instrument has broadband continuum

* Corresponding author. Fax: +46 (0)18 4715999. E-mail addresses: [email protected] (B.J.R. Davidsson), samuel.gulk [email protected] (S. Gulkis), [email protected] (C. Alexander), paul.a. [email protected] (P.von Allmen), [email protected] (L. Kamp), seung [email protected] (S. Lee), [email protected] (J. Warell). 0019-1035/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.icarus.2010.06.022

channels centered on 190 GHz (1.6 mm) and 562 GHz (0.5 mm), as well as a 44 kHz-resolution spectrometer connected to the submillimeter channel, capable of measuring rotational emission lines 17 18 of H16 2 O; H2 O; H2 O, CO, CH3OH, and NH3. In order to produce synthetic MIRO observations (needed for science operations planning and for facilitating the reconstruction of 3D gas properties of the comet coma from MIRO line of sight observations), substantial modeling efforts have recently been performed. The core of the MIRO spectral simulator is a code which employs the accelerated Monte Carlo method to solve the radiative transfer equation at the relevant wavelengths, using rate equations which account for non-equilibrium intermolecular and electronmolecule collisional excitation, absorption, spontaneous and stimulated emission as well as infrared pumping by solar radiation. The radiative transfer model, described by Lee et al. (2009), requires input regarding the molecular number density n, translational temperature Ttr, gas expansion velocity W as well as the distribution of random molecular velocities throughout the coma. In one set

456

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

of models, Lee et al. (2009) assumed a spherically symmetric single-species water coma without dust, with constant gas temperature and expansion velocity, leading to a R2 number density dependence on cometocentric distance R. Such simple coma models and the corresponding radiative transfer solutions are valuable for code validation purposes. The current paper describes a more realistic gas coma model, where density, temperature, expansion velocity and the local distribution of random velocities develop in accordance with the Boltzmann equation. These gas models are also used by Lee et al. (2009) in a second set of radiative transfer models, thereby yielding synthetic MIRO observations that are likely to resemble measured data. Here, the kinetic Direct Simulation Monte Carlo (DSMC) technique (Bird, 1994) is used to calculate macroscopic parameters (i.e., {n, Ttr, Trot, W}, where Trot is the rotational temperature) as well as microscopic properties (i.e., the complete molecular velocity distribution function g) of comet comae consisting of pure H2O or H2O/CO mixtures. In particular, conditions at large heliocentric distances (2.5–3.5 AU) are considered. Although the primary reason for the modeling is to provide input to the MIRO spectral simulator, additional calculations with relevance for other aspects of the Rosetta mission are performed as well. The dynamics of dust grains entrained in the gas is considered in order to place constraints on dust velocities and the largest grain size that can be lifted off the nucleus surface. Such calculations are not only scientifically motivated but are of interest for spacecraft safety reasons. Furthermore, gas drag forces acting on the Rosetta spacecraft itself are calculated and compared with the nucleus gravitational force. The drag-to-gravity force ratio has important implications for spacecraft dynamics and fuel budget. The DSMC technique is used here due to its capability of accurately describing gas kinetics even when the molecular velocity distribution function deviates strongly from the ideal Maxwell– Boltzmann distribution. Such deviations from thermodynamic equilibrium unavoidably arise near surfaces which sublimate into vacuum, forming so-called Knudsen layers (e.g., Cercignani, 2000; Crifo et al., 2002; Davidsson, 2008). At sufficiently high gas production rates, the gas equilibrates above the Knudsen layer and the exterior coma is in the hydrodynamical regime (except for very large distances, where breakdown takes place again). However, at sufficiently low gas production rates, thermodynamic equilibrium is never achieved and the entire comet coma is an extended Knudsen layer. This causes severe difficulties for hydrodynamical coma models (based on Euler or Navier–Stokes equations), since hydrodynamics itself assumes that thermodynamic equilibrium prevails. Hydrodynamical coma models therefore always face problems concerning inner boundary conditions, and may not be applicable at all for comets at large heliocentric distances. By using DSMC, such problems are avoided completely. Furthermore, we emphasize that standard gas drag force formulae also assume a Maxwell–Boltzmann distribution for the gas, and may therefore not be reliable. Therefore, we here calculate drag forces numerically, directly from DSMC molecular distribution functions. The dust grain dynamics and drag forces on Rosetta presented here may therefore be more accurate than calculations assuming that thermodynamic equilibrium prevails. To our knowledge, this is the first time that coma formation is studied over hot dust mantles (normally, the much cooler conditions relevant for exposed sublimating surface ice are studied). The initial water velocity distribution function is for the first time calculated self-consistently with the assumed surface conditions (regarding porosity of the medium and the subsurface temperature structure), whileas semi-Maxwellians typically are assumed but their correctness not validated. As mentioned previously, drag coefficients suitable for non-equilibrium gases are employed for the first time, whileas the rule normally is to use standard expres-

sions from hydrodynamics even when thermodynamic breakdown occurs and such expression no longer necessarily are valid. In spite of these improvements, the coma model also contains significant simplifications, such as a spherically symmetric geometry and an assumption that molecular dynamics is unaffected by molecular absorption of solar light, molecular emission of radiation and the presence of dust (these assumptions lead to isentropic gas expansion). Such conditions are highly idealized and will not be encountered at 67P/Churyumov–Gerasimenko except perhaps under special conditions. However, models with intermediate level of complexity are needed for comparison between simple approaches (like the set of n / R2 models in the Lee et al. (2009) paper) and highly sophisticated models (e.g., Crifo et al., 2002, 2003; Zakharov et al., 2008a,b; Tenishev et al., 2008). Furthermore, the currently considered model is interesting from an engineering point of view as it combines relatively low computational costs with a reasonable level of physical accuracy. The paper is structured as follows. Section 2 describes the theory and models used in this work, and Section 3 summarizes the results. Specifically, Sections 3.1–3.3 describe the gas comae at 1.3 AU, 2.5 AU and 3.5 AU, respectively, Section 3.4 deals with Doppler broadening of molecular emission lines, Section 3.5 treats dust dynamics, and Section 3.6 presents drag forces on Rosetta. Finally, the results are discussed in Section 4.

2. Theory and models The current work relies on five different models, designed to treat specific aspects of the highly complex comet sublimation and coma formation problem. I. A thermophysical nucleus model, considering an icy interior covered by a porous dust mantle, used to calculate gas production rates. II. A random-walk Monte Carlo code for studying molecular diffusion inside porous dust mantles. It is used for calculating probabilities for successful molecular migration through the dust mantle versus recondensation (input to the thermophysical model I). It also calculates the complete velocity distribution of molecules emerging from the dust mantle (input for the DSMC model III). III. A DSMC model which calculates the evolution of the H2O and CO velocity and rotational energy distribution functions throughout the coma, as well as the associated macroscopic parameters {n, Ttr, Trot, W}. IV. An integrator using numerical molecular velocity distribution functions to calculate drag forces on spherical and flat targets as function of cometocentric distance, target size, and target velocity. V. An equation of motion integrator used to calculate dynamics of solid targets under influence of nucleus gravity and gas drag forces. The thermophysical model (I) consists of the spatially onedimensional differential equation for heat conduction, solved for a domain composed of two layers with different densities q, specific heat capacities c and heat conductivities j. An upper layer with q = 980 kg m3, c = 830 J kg1, j = 0.22 J s1 m2 K1 and pffiffiffiffiffiffiffiffiffi resulting thermal inertia I ¼ qcj ¼ 420 J m2 K1 s1=2 acts as an ice-free dust mantle with 70% porosity. A lower layer with q = 925 kg m3, c = 1570 J kg1, j = 3.1 J s1 m2 K1 and I ¼ 2100 J m2 K1 s1=2 represents water ice, considered compact for simplicity. The upper boundary condition balances energy fluxes due to time-dependent solar illumination, thermal reradiation into space and heat conduction into the interior. The lower boundary

457

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

300

Temperature [K]

280

Dust mantle surface Mid−mantle (10cm depth) Ice/dust interface (20 cm depth) Icy interior (30 cm depth)

260 240 220 200 180 0

2

4

6

8

10

12

Local time [h] 300 Midnight Sunrise Noon Sunset

Temperature [K]

280 260 240 220 200 180 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

Depth [m] Fig. 1. The upper panel shows the physical nucleus temperature versus time at four different depths during a full nucleus rotation. Equatorial conditions on a nucleus having its spin axis perpendicular to the orbital plane and rh = 1.3 AU are considered. The lower panel shows the temperature versus depth for four different instances during nucleus rotation.

condition at a large depth consists of a vanishing temperature gradient. A sink term at the ice/dust interface governs energy consumption due to sublimation. The net sublimation rate includes a correction for recondensation of molecules unable to diffuse through the mantle. The molecular escape probability pesc depends on the dust mantle thickness, porosity, and radius of the largest cavities of the mantle (here taken as 2 mm), and is calculated using a Monte Carlo molecular random-walk approach (model II). Trialand-error calculations, using 10 latitudinal slabs to estimate the total water production rate of a spherical nucleus with radius Rnuc = 1700 m, the equal-volume spherical equivalent of 67P/Churyumov–Gerasimenko (Lamy et al., 2007), resulted in a preferred dust mantle thickness of 20 cm and a corresponding pesc = 2.8%. Under such conditions, the total model water production rate at rh = 1.3 AU amounts to Q H2 O ¼ 1:18  1028 molec s1 , similar to that of Comet 67P/Churyumov–Gerasimenko near perihelion (Osip et al., 1992; Schleicher, 2006). This dust mantle thickness was kept for all three considered heliocentric distances (rh = 1.3 AU, 2.5 AU, and 3.5 AU) to produce a self-consistent decline of comet activity with diminishing solar flux. Fig. 1 exemplifies the model calculation at rh = 1.3 AU by showing temperature versus time at various depths (upper panel), and temperature versus depth at various instances of time (lower panel). Note that substantial diurnal temperature variations take place in the mantle, but that the ice/dust interface is virtually isothermal due to sublimation at the interface and the relatively high thermal inertia of the ice. For a given latitude slab (here, the equator on a nucleus with its spin pole perpendicular to the orbital plane is considered), the gas production rate is therefore quasi-constant during day and night. For more details on model I, see Appendix A. The random-walk Monte Carlo model (II) tracks individual test particles from ice release to escape into the coma, or recondensation at the ice/dust interface. Every time a test particle interacts with the dust mantle, a new direction of motion is drawn randomly from an isotropic angular distribution, a new velocity is drawn randomly from a Maxwellian distribution evaluated at the local temperature, and the distance to the next interaction point is

selected. The dust mantle temperature distribution is taken from noon temperature-profiles calculated with model I (e.g., lower panel of Fig. 1). Therefore, the molecular velocity distribution naturally adjusts to the temperature profile of the dust mantle. Molecular flight distances are drawn randomly from a pre-calculated distribution of free paths which has been sampled by letting a test particle bounce around inside a computer-generated medium consisting of randomly positioned equal-sized spheres, here with a global porosity of 70%. Test particles encountering ice either scatter or condense according to a temperature-dependent recondensation coefficient. Model II is used to calculate both pesc and the distribution function of escaping molecules (i.e., the distribution of molecular speeds for a large number of angles to the surface normal). Generally, the combination of a temperature gradient and the physical obstacle presented by the dust mantle, generates an initial distribution function that deviates from the ideal semi-Maxwellian distribution that most authors use by default, primarily by being more collimated. Note that models I and II only consider water sublimation. Pure water comae are calculated for all considered heliocentric distances, rh = 1.3 AU, 2.5 AU and 3.5 AU. At 3.5 AU additional models are calculated for H2O/CO mixtures. It is then assumed that the CO net production rate is 1 or 65 that of water. A semi-Maxwellian initial distribution function is selected for CO molecules, evaluated for the dust mantle surface temperature (assuming that the last CO interaction with the dust mantle occurred near the surface). For more details on model II, see Appendix B. Coma simulations are made with a standard DSMC code1 written by Bird (1994) and modified substantially by the authors (model III) for comet modeling applications. It considers a spherically symmetric geometry and uses noon production rates and distribution functions calculated by models I and II as initial conditions. The DSMC code tracks the positions and velocity vectors of individual test particles, each representing a certain number of H2O or CO mol1 The code DSMC1.FOR is distributed by Bird (1994) and available at URL http:// www.aeromech.usyd.edu.au/dsmc_gab/.

458

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

ecules. Velocity vectors are reoriented during molecular collisions according to a Variable Hard Sphere (VHS) model which realistically mimics temperature-dependent viscosity of the gases by having velocity-dependent collision probabilities. Exchange between kinetic and internal energy modes during inelastic collisions are governed by the Larsen–Borgnakke model. Here, water molecules have three degrees of rotational freedom while carbon monoxide molecules have two. The velocity vectors (and rotational energies) of test particles are sampled repeatedly for a number of cells throughout the considered volume, to build the local velocity and rotational energy distribution functions. These distribution functions can then be used to calculate the corresponding macroscopic properties of the cell, e.g., {n, Ttr, Trot, W}. It has been shown (Wagner, 1992) that Bird’s DSMC algorithm converges towards the Boltzmann equation, i.e., it reproduces the underlying physical processes described by that equation. A DSMC solution is therefore equivalent, within numerical errors, to a solution to the Boltzmann equation. Hence, DSMC is more general than, e.g. hydrodynamical modeling, since the governing equations of hydrodynamics are based on moments of a simplified version of the Boltzmann equation (i.e., the limit where the collision integral vanishes). Systematic differences between DSMC solutions and hydrodynamical calculations are therefore seen when the collision integral is not negligible (i.e., deviations from thermodynamic equilibrium are strong). The DSMC technique (model III) is described in more detail in Appendix C. If the complete velocity distribution function is known, it can be used to calculate drag forces on targets with arbitrary geometry by considering the momentum transfer to the target during collisions with individual molecules in the gas. Model IV uses distribution functions from model III (DSMC) for numerical integration of drag forces for targets shaped as spheres and plates. The former shape is here used to represent coma dust grains, while the latter shape is used to represent Rosetta (since the primary drag arises at the flat solar panels). Target velocity relative to the nucleus is easily accounted for by recalculating the gas velocity distribution function to a system co-moving at the target speed. Since individual distribution functions are available for the species H2O and CO, drag forces from both gases are calculated and added. Model IV assumes that all molecules are specularly reflected. This may over- or underestimate the drag force if there are significant temperature differences between the target and the gas, since specular reflections do not lead to molecular speed modifications by thermalization with the target. However, allowance for grain temperature (and its variation with grain size and solar exposure time) would complicate the problem substantially without necessarily leading to more realistic results (e.g., optical grain properties, responsible for their temperature, remains a guess-work). Due to the negligence of the uncertain effects of diffuse scattering, we estimate that our drag forces are off by less than 40% near the nucleus, where most of the acceleration takes place. In model V, the equation of motion is solved for grains being subjected to nucleus gravity and the gas drag force obtained from model IV. This provides grain position and speed versus time, including the terminal velocity of dust grains. Comparison between the drag force in the immediate vicinity of the nucleus and the local gravitational force yields an estimate of the largest grains which can be lifted off the nucleus. More details are given in Appendix D.

temperature is Tsurf = 295 K, and the local water production rate is Z = 4.2482  1020 molec m2 s1 (corresponding to a total production of Q = 1.5428  1028 molec s1 if Z is applied everywhere on a Rnuc = 1700 m sphere). The near-nucleus molecular mean free path is 0.5 m, growing rapidly with height. Fig. 2 shows the number density, temperatures, and gas expansion velocity obtained in the DSMC simulations. In addition to Ttr and Trot, Fig. 2 also shows the temperature components Tx = Ty and Tz, which may require some explanation (see also Eq. (C.13)). Temperature is a measure of the kinetic energy related to random molecular motion in a gas (i.e., a measure of the velocity field seen in a frame co-moving at W). Since velocity vectors can be decomposed into three orthogonal components, one can also define a b temperature for each orthogonal direction {Tx, Ty, Tz}, where ^z ¼ R. Equality between these components (taking the same value as Ttr) shows that the gas is isotropic in a frame moving at the expansion (or bulk drift) speed W, and is a strong sign of thermodynamic equilibrium. However, differences between the temperature components show that the velocity distribution function is nonisotropic in the co-moving frame and thereby non-Maxwellian. Difference between the translational and rotational temperatures further illustrates that the gas is not in equilibrium, since the total kinetic energy in a volume differs from the total internal energy. The middle panel of Fig. 2 shows a strong difference between temperature components in the lower 15 m of the coma, and is de facto a visible manifestation of the Knudsen layer. In the region above the Knudsen layer the gas is in thermodynamic equilibrium, which means that the DSMC solution can be cross-checked against hydrodynamical calculations. This can be done by solving the Euler or Navier–Stokes equations, but since the geometry is trivial (spherically symmetric expansion), much simpler options are available, which does not require the solution of differential equations. For isentropic expansion, temperature and number density are related by

3. Results

where the Mach number M is

3.1. The gas coma at perihelion

M¼W

We here summarize our results for rh = 1.3 AU, the perihelion distance of Comet 67P/Churyumov–Gerasimenko. The ice/dust interface temperature is Tice = 205 K, the dust mantle surface

and indices 1 and 2 here refer to two different cometocentric distances (or more generally, two different cross sections through which gas is flowing).

T tr n1c ¼ C;

ð1Þ

where c = (f + 5)/(f + 3) is the specific heat ratio, f is the number of rotational degrees of freedom, and C is a constant. The number den_ by sity and bulk drift velocity are related to the net mass flux m

_ ¼ mnW; m

ð2Þ

and the net mass flux changes with height as

_ ¼m _0 m

 2 Rnuc ; R

ð3Þ

_ 0 is the net mass flux at the nucleus surface (equal to the where m gas production rate Z since all returning coma molecules are diffusely scattered at the dust mantle surface). Combining Eqs. (2) and (3) yields,

nW ¼

  _ 0 Rnuc 2 m : m R

ð4Þ

Finally, temperature and bulk drift speed are related by

! cþ1  2 2 2ðc1Þ R2 M1 1 þ 12 ðc  1ÞM1 ¼ ; R1 M2 1 þ 1 ðc  1ÞM22 2 rffiffiffiffiffiffiffiffiffiffiffi m ; ckT tr

ð5Þ

ð6Þ

459

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

−3

n [molec m ]

18

10 17 10 16 10 15 10 14 10 13 10 12 10 −1 10

0

10

1

10

2

10

3

10

4

10

5

6

10

10

T [K]

300

Txy T z Ttr T

200 100 0 −1 10

rot

0

10

1

10

2

10

3

10

4

10

5

10

6

10

W [m s−1]

1000 800 600 400 200 −1 10

0

10

1

10

2

10

3

10

4

10

5

10

6

10

Height above surface [m] Fig. 2. Water coma conditions at rh = 1.3 AU. The upper panel shows the number density calculated with DSMC (thick solid line) compared to a n / R2 law fitted to the number densities at large cometocentric distances (thin solid line). The middle panel shows the translational temperature Ttr, its three components {Tx, Ty, Tz}, and the rotational temperature Trot according to the legend. The lower panel shows the gas expansion velocity calculated with DSMC (solid line). In all panels, squares show the solution according to isentropic hydrodynamical relations, and vertical lines indicate the sonic point.

Eqs. (1), (4) and (5) thereby constitute three relations that can be solved for the three unknowns {n(R2), Ttr(R2), W(R2)}, as long as the values of {n, Ttr, W} are provided for a single point R1 in the _ 0 ¼ Z is already coma (this condition will fix C and M1 while m known). Here, the sonic point is chosen, i.e. the point where the temperature and expansion velocity are such that the Mach number equals unity according to Eq. (6). The sonic point (M1 = 1) is located at R1 = Rnuc + 14.001 m, where {n, T, W} = {1.1327  1018 molec m3, 221.02 K, 368.95 m s1}, all according to the DSMC solution. The sonic point is marked with a vertical line in Fig. 2 and the analytical hydrodynamic solution is marked by squares (extending far beyond the point where DSMC calculations had to be interrupted in order not to suffer excessive computer CPU time). In the region where DSMC and hydrodynamical solutions overlap they are in perfect mutual agreement, illustrating the validity of the DSMC code implementation. Note that the hydrodynamical solution cannot be extended closer to the nucleus than the sonic point. The right hand side of Eq. (5) cannot take values below unity when M1 = 1, for any value of M2, i.e., hydrodynamical solutions for R2 < R1 do not exist. To enter into the R2 < R1 region (obviously containing gas), the molecular velocity distribution function must necessarily become non-Maxwellian, which invalidates hydrodynamics as such and all relations following from it, e.g., Eq. (5). Stating it differently – solid or porous media which inject molecules into space (in a realistic manner) will inevitably create a non-equilibrium gas that needs time to reach thermodynamical equilibrium through molecular collisions. For this reason, the lower part of the comet coma constitutes a kinetic Knudsen layer, which separates the nucleus from the region of the coma where thermodynamical equilibrium prevails. Note particularly, that the sonic point coincides with the height where differences between Tx,y, Tz, Ttr and Trot first become negligible. Of course, the Boltzmann equation is valid at all distances from the surface, the only important thing being that the collision integral

in that equation takes non-zero values for R < R1 and is zero for R P R1. Since a hydrodynamic solution never can extend entirely to the nucleus surface, it also means that measured downstream (hydrodynamic) coma conditions never can be used to infer the nucleus surface conditions (e.g., the surface temperature), unless a proper kinetic model is used to connect the two. According to the hydrodynamical solution, the expansion velocity W increases gradually from W  200 m s1 near the nucleus, to a terminal value just below W  1 km s1, reached a few 100 km from the nucleus. When W is constant, the gas density falls as n / R2. A R2 law has been fitted to the R > 1000 km part of the hydrodynamic n(R) solution (thin solid line in the upper panel of Fig. 2), and extrapolated towards the nucleus. Note how such an extrapolation underestimates the near-nucleus coma density by a factor 4.4. In order to investigate whether the isentropic solution is reliable at large heights (where hydrodynamics breaks down anew), a DSMC simulation was performed in the R = 20 km–300 km region. At the inner boundary, the local isentropic solution was used as input. Weak signs of non-equilibrium appear in the distant flow, but at the outer boundary at 300 km the difference between the DSMC and isentropic solutions are merely 0.6% in n, 3 K in Ttr and 6 m s1 in W. Since the terminal temperature and expansion velocity nearly have been reached, we therefore consider the isentropic solution, plotted up to a distance of 1000 km in Fig. 2 to be representative of the real flow. Finally, we consider the effect of using a numerically sampled initial velocity distribution (model II). In the ideal case, molecules would cross the surface according to a Maxwellian transmission distribution function (see, e.g. Huebner and Markiewicz, 2000). If this distribution is integrated over the azimuth angle, a speed bin dV centered on V, and a bin dh centered on the angle h from the surface normal, the number of molecules entering this {V, h}-bin (m2 s1) would be given by

460

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

NMTDF ¼ 4Zb4 V 3 expðb2 V 2 Þ cos h sin h dV dh:

tion of the population at a slightly larger depth (since many molecules are scattered back towards the underlying ice, while a continuously decreasing population migrates towards the surface). The consequence is that the emerging flow becomes more collimated than NMTDF. The right panel of Fig. 3 compares the translational temperature and expansion velocity calculated from the numerical initial conditions (thin curves), with a calculation based on the ideal transmission distribution function (thick curves). Due to the higher degree of collimation for the sampled distribution with respect to NMTDF, the initial expansion velocity is higher by 11 m s1. Since the level of random motions is smaller, the initial translational temperature is 18 K lower. Since both flows have the same net mass flux, the larger expansion velocity for the numerical initial conditions lead to a 5% lower initial number density. More realistic initial conditions therefore lead to a somewhat thinner, faster and cooler medium, compared to the classical semi-Maxwellian initial condition. At larger distances from the surface (2.4 km), initial semi-Maxwellian conditions leads to a gas that is 4 K warmer, 18 m s1 faster, and 3% thinner than the more realistic model. The reason for this difference is that the more realistic model includes some molecules arriving directly from a slightly cooler interior, i.e., the total energy is somewhat lower than for a transmission distribution function evaluated at Tsurf = 295 K. Note that the difference between the two approaches here is minimized since the considered dust mantle is rather thick (20 cm). For a substantially thinner mantle (e.g. 2 cm), a large fraction of molecules would originate directly from the cool ice, while the rest would scatter one or several times and acquire an enhanced velocity by interacting with the hot dust. Such a mixed population can never be well represented by a Maxwellian transmission distribution function evaluated for a single temperature, and usage of numerically calculated initial distribution functions would be mandatory.

ð7Þ

Note that the transmission distribution function describes the properties of molecules that simultaneously cross a surface. As the molecules spread into a volume above the surface, slower molecules accumulate while faster molecules become underrepresented with respect to the transmission distribution function. The population in the volume is described by the semi-Maxwellian velocity distribution function, i.e. a truncated non-drifting Maxwellian with velocity vectors confined to the outward direction. The left panel of Fig. 3 shows a contour plot of NMTDF (thick lines) evaluated for the surface temperature Tsurf = 295 K, using bins of width dV = 10 m s1 and dh = 1°. The left panel also shows the corresponding sampled distribution (thin lines), which is more collimated than NMTDF, i.e., peaks for a lower h-value. The difference can be understood as follows. Davidsson and Skorov (2004) illustrated that an isothermal, and in all points sublimating, porous medium produces a transmission distribution function that (after appropriate integration) is identical to NMTDF. However, if temperature decreases with depth the sublimation rate is comparably large very close to the surface, which statistically enhances the probability for molecules to escape at large h-angles. The corresponding flow is therefore less collimated than NMTDF. Conversely, if temperature increases with depth the flow becomes more collimated than NMTDF. This is because most molecules are produced at comparably large depths, from where they escape predominantly along the surface normal. The currently considered medium (a non-sublimating thick dust mantle) behaves similarly to the latter case. This can be understood by realizing that, in both cases, the number density of molecules increases rapidly with depth. For the sublimating medium, this is because the production rate of molecules increases with depth. For the non-sublimating medium, it is because the population at a shallow depth constitutes a frac-

600 250

80

60

50

40

30

400 200

Expansion velocity [m s−1]

500

Translational temperature [K]

Angle from outward normal θ [deg]

70

300 20

10

200

400

600

800

1000

Molecular speed V [m s−1]

150

0

2

10

200

10

Height [m]

Fig. 3. Comparison between the usage of numerically sampled initial conditions and using an idealized semi-Maxwellian. The left panel shows contour plots of the numerically sampled transmission distribution function (thin lines), compared to the Maxwellian transmission distribution function evaluated for Tsurf = 295 K (thick lines), both integrated over azimuth and {V, h} bins as described in the text. The right panel shows the translational temperature (solid) and expansion velocity (dashed) obtained for two different initial conditions: model II (thin curves) and a Maxwellian transmission distribution function (thick curves).

461

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

3.2. The gas coma at 2.5 AU The ice/dust interface temperature at 2.5 AU is Tice = 182 K, the dust mantle surface temperature is Tsurf = 201 K, and the local water production rate is Z ¼ 1:0173  1019 molec m2 s1 (corresponding to a modeled total production of Q = 3.6945  1026 molec s1 if Z is applied everywhere on a Rnuc = 1700 m sphere). No measurements of the water production rate are available for Comet 67P/Churyumov–Gerasimenko at such large heliocentric distances, but Schleicher (2006) report Q = 9.5 (±0.8)  1026 molec s1 at rh  1.9 AU, which can be scaled to Q  5.5  1026 molec s1 at rh = 2.5 AU, assuming a Q / r 2 h relation. This shows that the applied thermophysical model (for which the dust mantle thickness was selected to reproduce the perihelion production rate) yields roughly the correct Q = Q(rh) behavior. The near-nucleus molecular mean free path is 9.2 m. Fig. 4 shows the DSMC solution for number density, temperatures and expansion velocity in the lower 300 km of the water coma at 2.5 AU from the Sun. The temperature evolution with height is completely different compared to rh = 1.3 AU conditions. The rotational temperature Trot is several tens of Kelvin higher than the translational temperature Ttr at all times. The temperature components Tx,y and Tz merge temporarily around R  500 m, but then depart and differ several tens of Kelvin at large distances. Evidently, molecular collisions occur too infrequently to allow efficient energy exchange between kinetic and internal modes, and the gas is incapable of maintaining a Maxwellian distribution function against the disruptive process of expansion. The transverse temperature component Tx,y tends to zero, showing that molecules travel radially  at large distances. However, Tz tends towards a nonzero value T 1 z  24 K , indicating a remaining velocity dispersion for the molecules along the radial. Had the number density been sufficiently high to allow frequent collisions (i.e., if hydrodynamics had prevailed), the molecules would have shared their kinetic energy equally between themselves, so that all molecules eventually

travelled radially at exactly W (i.e., no random motions, and 1 T1 z ¼ 0 K). This would have yielded T tr ¼ 0 K, as seen for the hydrodynamical solution at 1.3 AU (see Fig. 2). Instead, the translational temperature at 2.5 AU levels off at T 1 tr ¼ 9 K, a phenomenon called freeze-out (e.g., Sone and Sugimoto, 1993). Note that the terminal expansion velocity of W1 = 750 m s1 is some 250 m s1 lower than at 1.3 AU, which primarily is a consequence of the cooler dust mantle at 2.5 AU, giving molecules less initial kinetic and rotational energy. 3.3. The gas coma at 3.5 AU The ice/dust interface temperature at 3.5 AU is Tice = 156 K, the dust mantle surface temperature is Tsurf = 166 K, and the local water production rate is Z ¼ 4:2611  1016 molec m2 s1 (corresponding to a total production of Q = 1.5475  1024 molec s1 if Z is applied everywhere on a Rnuc = 1700 m sphere). Nothing is known about the real water production rate of Comet 67P/Churyumov–Gerasimenko at this large distance, but we rely on a model that produces roughly the correct production rates at 1.3 AU and 2.5 AU. Note, however, that if exposed surface ice is present on Comet 67P/Churyumov–Gerasimenko, the large-distance gas production behavior may deviate substantially from that predicted by the current model, which assumes the presence of a global dust mantle. The assumed carbon monoxide production rate is either zero, 1 or 65 that of water (the latter two rates being QCO = 1.5475  1024 molec s1, and QCO = 1.0059  1026 molec s1). Such dominance of CO over H2O in comet comae at large heliocentric distances (rh J 3 AU) was measured for, e.g., Comet C/1995 O1 (Hale-Bopp) and is due to the fact that the condensation temperature of CO is substantially lower than that of water (Biver et al., 2002). An upper limit on the carbon monoxide production rate of QCO 6 9.2  1026 molec s1 at 3.0 AU for Comet 67P/Churyumov– Gerasimenko was established by Bockelée-Morvan et al. (2004).

17

n [molec m−3 ]

10

Actual profile n ∝ R−2

16

10

15

10

14

10

13

10

12

10

11

10

1

10

2

10

3

10

4

10

5

10

200 Txy Tz Ttr Trot

T [K]

150 100 50 0

1

10

2

10

3

10

4

10

5

10

−1 W [m s ]

800 600 400 200

1

10

2

10

3

10 Height above surface [m]

4

10

5

10

Fig. 4. Water coma conditions at rh = 2.5 AU (see the caption of Fig. 2 for details). Note that no hydrodynamic model can be calculated for rh = 2.5 AU since the gas never achieves thermodynamic equilibrium.

462

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

Fig. 5 shows the behavior of {n, Ttr, Trot, W} for the pure-water coma as solid curves. The rotational temperature is virtually constant throughout the inner coma, illustrating that very few collisions take place (the molecular mean free path near the nucleus is 2.2 km, i.e. comparable to Rnuc). In reality, Trot would decrease due to spontaneous emission, a process that is counteracted by IR excitation through absorption of solar radiation. In the current isentropic calculations such processes are neglected. The decrease of Trot in Figs. 2 and 4 illustrate the importance of energy transfer between rotational and kinetic modes through collisions, which always takes place in parallel with radiative processes when gas densities are sufficiently high. The lack of Trot decrease in Fig. 5 for a pure-water coma at 3.5 AU is not expected in real observations, but here serves the purpuse of showing that the gas essentially is collisionless. For more realistic calculations of Trot for the currently considered coma models we refer to Lee et al. (2009). They calculate the effect of absorption, spontaneous and stimulated emission, as well as infrared pumping on Trot for a number of different H2O rotational transition lines, showing that they vary substantially, both with respect to each other and as function of distance to the nucleus. The problem of rotational excitation of water molecules in comet atmospheres has been discussed and modeled thoroughly by e.g. Crovisier (1984), Bensch and Bergin (2004), Zakharov et al. (2007). The translational temperature freezes out at T 1 tr ¼ 22 K, and the terminal expansion velocity is W1 = 460 m s1, roughly half of the perihelion value. Note that the differences between {Tx, Ty} and Tz at 3.5 AU (not plotted here) are even more extreme than at 2.5 AU (middle panel of Fig. 4). Adding a carbon monoxide production at the same rate as that for water (dashed-dotted curves in Fig. 5) has virtually no effect on the water dynamics since the number of H2O–CO collisions are too few. However, when carbon monoxide dominates the coma (pro-

duction rate ratio H2O:CO = 1:65), the water feels the presence of the CO (dotted curves in Fig. 5). The water rotational temperature now cools substantially with height, as internal energy is transfered to kinetic energy. The translational temperature freeze-out is lowered to T 1 tr ¼ 12 K. The near-nucleus water expansion velocity drops by about 80 m s1, as the water molecules have to diffuse through the heavier and thereby slower carbon monoxide. The reduced expansion velocity leads to a corresponding 40% increase in the near-nucleus water number density (to keep the same mass flux). However, the terminal water expansion velocity is somewhat higher when substantial amounts of CO is present, compared to pure water, due to the previously mentioned exchange between kinetic and internal energy modes. The behavior of carbon monoxide at rh = 3.5 AU is seen in Fig. 6. Thick curves in Fig. 6 denote water for reference and are identical to those in Fig. 5, while thin curves show CO. The translational temperatures do not differ substantially between the species, particularly not when H2O:CO = 1:65. In the high-QCO case, molecules collide frequently enough to cause a roughly equal kinetic energy partition between the species. The rotational temperature is quasi-constant for both species when H2O:CO = 1:1, but falls significantly when H2O:CO = 1:65. For the CO-dominated coma, the terminal rotational temperature is T 1 rot ¼ 94 K for CO, while the water rotational temperature is some 20 K cooler. The reason why Ttr becomes roughly equal for the two species while Trot does not, is that only every fifth collision here is considered inelastic (see Appendix C). Energy transfer between kinetic and internal modes is too slow to allow the internal energies of H2O and CO to equilibrate. The drift speed of CO is almost 100 m s1 lower than that for H2O when their production rates are comparable, due to the molecular mass difference. When the carbon monoxide production rate increases, CO–CO collisions are more efficient in turning inter-

14

−3 n [molec m ]

2.5

x 10

H2O:CO=1:0 H2O:CO=1:1 H2O:CO=1:65

2 1.5 1 0.5 2

10

3

10

4

10

5

10

T [K]

150 100 50 0

Ttr Trot 2

10

3

10

4

10

5

10

W [m s−1]

500 400 300 200

2

10

3

10

4

10

5

10

Height above surface [m] Fig. 5. Water coma conditions at rh = 3.5 AU. The upper panel shows the water number density n versus height above the nucleus surface. The middle panel shows water translational temperature Ttr (thin curves) and water rotational temperature Trot (thick curves). The lower panel shows the water expansion velocity W. In all panels, line styles identify H2O:CO production ratios according to the legend in the upper panel.

463

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

200 H2O (1:1)

T [K] tr

150

CO (1:1) H2O (1:65)

100

CO (1:65)

50 0

2

10

3

4

10

10

5

10

200

100

T

rot

[K]

150

50 0

2

10

3

4

10

10

5

10

−1

W [m s ]

500 400 300 200 2

10

3

4

10

10

5

10

Height above surface [m] Fig. 6. Carbon monoxide and water coma conditions at rh = 3.5 AU. The upper panel shows translational temperature Ttr, the middle panel shows rotational temperature, and the lower panel shows expansion velocity. In all panels, line styles identify H2O:CO production ratios according to the legend in the upper panel.

nal energy into kinetic energy, hence W increases significantly. At this high CO production rate, the water is forced to flow virtually at the same expansion velocity as the carbon monoxide. Finally, it is noted that the entire 3.5 AU coma is severely out of equilibrium for both water and carbon monoxide, as {Tx, Ty} and Tz differ by many tens of Kelvin.

3.4. Doppler broadening of molecular emission lines One way to quantify the level of deviation of a velocity distribution function with respect to an ideal Maxwellian, is to calculate the difference between the temperature components Tx = Ty and Tz. Another method is to calculate the Doppler broadening of emission lines due to the internal velocity field of the gas and study variations in the line profile, e.g., with changing viewing geometry. We here present examples of the latter technique, considering isolated parcels of gas taken from the simulations at different heliocentric and cometocentric distances. We emphasize that these examples primarily are of academic interest. The intention is to illustrate properties of the numerically sampled velocity distribution functions, not to evaluate whether MIRO actually could detect deviations from thermodynamic equilibrium. Such an assessment would require that radiative transfer calculations are performed for gas along an entire line of sight. Furthermore, instrumental effects (e.g., spectral and spatial resolution, signal to noise ratio) must be considered. This is out of scope in the current investigation. However, variations in the line FWHM versus viewing angle are compared with the MIRO spectral resolution to determine whether the detection of deviations from thermodynamic equilibrium is excluded solely based on spectral resolution considerations. If the rest frequency (in a cometocentric system) of a molecular emission line is f0, molecular motion will create a Doppler broadening profile U = U(Df) as a function of frequency f relative to f0 (Df = f  f0), that will be partially responsible for the observed line

shape. Since the DSMC modeling provides the full molecular velocity distribution function, U(Df) can be evaluated numerically. Fig. 7 shows examples (squares) of numerically evaluated U(Df) profiles for four viewing geometries (a is the angle between the b and the line of sight b outward normal R L), using the velocity distribution function sampled within a parcel of gas at a height of 20 m above the nucleus surface, at rh = 2.5 AU. At a = 180° (i.e., gas expansion towards the observer), the line is blue-shifted, while it is red-shifted for a = 0° (i.e., gas expansion away from the observer). The considered line is the 1(1, 0)–1(0, 1) transition of H16 2 O with f0 = 556.936002 GHz. If the molecular distribution function is Maxwellian, and the gas b has translational temperature Ttr and bulk drift velocity W along R, the Doppler broadening profile is given by

0 !2 1 bb 1 1 Wð R LÞf 0 A; UðDf Þ ¼ pffiffiffiffi exp @ 2 Df  c pD D

ð8Þ

where

D¼

rffiffiffiffiffiffiffiffiffiffiffi 2kT tr f0 ; m c

ð9Þ

and c is the speed of light in vacuum. Best-fit solutions {Ttr, W}fit (in a v2-minimization sense) using Eq. (8) have been fitted to the numerically sampled data and are shown in Fig. 7 as solid curves. The fitted bulk velocity of the gas is Wfit = 246–249 m s1 and corresponds rather well to the actual bulk velocity of the gas in the parcel, W = 250.5 m s1. However, the fitted translational temperature varies with a and takes values in the range 133 6 T fit tr 6 148 K (135° 6 a 6 180°), which is substantially smaller than the actual translational temperature for the parcel, Ttr = 153.5 K. As a approaches 90°, the fitted temperature exceeds 160 K. The a-dependence of T fit tr , which would not be seen for a Maxwellian gas, is thus a manifestation of the breakdown of thermodynamic equilibrium.

464

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

Φ(Δ f)× 106

α=180°: {T,W}fit={133.2 K, 249 ms−1}

α=135°: {T,W}fit={147.9 K, 246 ms−1}

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

−2

−1

0

1

0

2

Φ(Δ f)× 106

α=45°: {T,W}fit={147.9 K, 246 ms−1} 1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

−2

−1

0

1

−1

0

1

2

α=0°: {T,W}fit={133.2 K, 249 ms−1}

1

0

−2

0

2

Δ f=f−f0 [MHz]

−2

−1

0

1

2

Δ f=f−f0 [MHz]

Fig. 7. Doppler broadening line profiles calculated for the H16 2 O 1(1, 0)–1(0, 1) line at f0 = 556.936002 GHz from numerical molecular distribution functions (squares). Solid curves show Eq. (8) fitted to the data (best-fit values {Ttr, W}fit given for each panel). Heliocentric and cometocentric distances are rh = 2.5 AU and R  Rnuc + 20 m, b and line of sight b respectively. The angle between gas expansion direction R L is denoted by a (as given for each panel).

forces calculated in such a manner are seen in the upper panel of Fig. 8 (solid and dashed-dotted curves), considering rh = 1.3 AU, a = 1 lm, q = 1000 kg m3 and two different dust grain velocities along the outward radial direction, vd = 0 m s1 and 300 m s1. Note that vd = 0 m s1 grains will experience a net outward force at any cometocentric distance, while a vd = 300 m s1 grain will feel a strong ‘‘air resistance” when vd > W (i.e., within 2 m of the surface according to the lower panel of Fig. 2), resulting in a net force towards the nucleus. Since the gas at R J 20 m is in thermodynamical equilibrium, the accuracy of the numerically calculated drag force can be evaluated by comparing with an analytical formula that depends on {n, Ttr, W} of the gas, the velocity vd and cross section A of the target. The formula assumes that the surrounding gas has a Maxwellian velocity distribution function. If the target is shaped as a sphere (a is the grain radius and A ¼ pa2 is the cross section), then the drag force is given by

It is emphasized that the same procedure repeated at rh = 1.3 AU, at a height of 1 km above the nucleus yields {Ttr, W}fit = {134.4 ± 0.2 K, 676 ± 0 m s1} for all considered a values, which agrees well with the actual values {Ttr, W} = {134.1 K, 676 m s1}. As long as the gas is Maxwellian, the method recovers the expected values. Examples of {Ttr, W}fit solutions for other heliocentric and cometocentric distances are collected in Table 1. At rh = 3.5 AU and R = 100 km, the fitted translational temperature changes by 34 K as a decreases from 180° to 135°. The corresponding decrease in the FWHM of the line profile amounts to nearly 230 kHz, or 5 MIRO resolution elements. This result shows that the finite spectral resolution of MIRO does not exclude that deviations from thermodynamic equilibrium can be detected. A more detailed assessment is beyond the scope of this work.

3.5. Dust dynamics

8 < F drag ¼ 12 AC D nðW  v d Þ2 ; 2 2 4 2 : C D ¼ ð2W þ1pÞffiffiffiexp ðW Þ þ ð4W þ4W 1Þerf ðWÞ : pW3 2W4

The drag force on spherical dust grains are here evaluated by calculating the momentum transfer to the grain from specularly reflected molecules, using the numerical molecular velocity distribution function sampled during DSMC modeling. Examples of drag

ð10Þ

Table 1 Columns 2–5 shows the best fit {Ttr (K), W (m s1)} parameters to numerical Doppler broadening profiles (for the cases specified in column 1), if Eq. (8) is used. Generally, that method does not manage to recover the actual values of {Ttr, W} (sixth column) since Eq. (8) assumes that the gas has a Maxwell–Boltzmann velocity distribution, which not necessarily is the case. The change of the line FWHM (measured in number of MIRO 44 kHz spectral resolution elements) when going from a = 180° to 135° is given in column 7. These numbers increase by a factor 2 when going to a = 90°. {rh, R}

{Ttr, W}a=180°

{Ttr, W}a=135°

Translational temperature and expansion velocity from the 1(1, 0)–1(0, 1) line of {1.3 AU, 1 km} {134.5, 676} {134.6, 676} {2.5 AU, 20 m} {133.2, 249} {147.9, 246} {2.5 AU, 100 km} {23.9, 745} {13.9, 745} {3.5 AU, 100 km} {67.7, 436} {33.6, 436} (H2O:CO = 1:0)

{Ttr, W}a=45°

{Ttr, W}a=0°

{Ttr, W}actual

DFWHM

{134.2, 676} {133.2, 249} {23.9, 745} {67.7, 436}

{134.1, 676} {153.5, 250.5} {11.3, 749} {22.4, 462}

0 1.3 2.5 5.2

H16 2 O {134.6, 676} {147.9, 246} {13.9, 745} {33.6, 436}

465

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471 −14

4

x 10

vd=0 m s−1 vd=300 m s−1

Drag force [N]

3 2 1 0 −1 −1 10

0

10

1

2

10

10

3

10

4

10

5

10

1000

−1 Velocity [m s ]

800

Gas Dust (a=1μm, ρ=1000 kg m−3)

600 400 200 0 −1 10

0

10

1

2

10

10

3

10

4

10

5

10

Height above surface [m] Fig. 8. The upper panel shows drag forces calculated from numerical molecular distribution functions, considering rh = 1.3 AU, a = 1 lm, q = 1000 kg m3, at two different grain velocities vd according to the legend (solid and dashed-dotted curves). The analytical drag formula valid under hydrodynamic conditions, Eq. (10), is also shown as squares (vd = 0 m s1) and rings (vd = 300 m s1). The lower panel shows the grain velocity versus height when q = 1000 kg m3 (thick solid line), as well as for q = 500 kg m3 and q = 1500 kg m3 (dashed-dotted curves). Note that the isentropic gas solution was used for the R > 4.1 km region. The lower panel also shows the gas expansion velocity (thin solid curve) for comparison.

Here, erf(x) is the error function and W is given by

W ¼ ðW  v d Þb;

ð11Þ

where

b¼

rffiffiffiffiffiffiffiffiffiffiffi m ; 2kT tr

ð12Þ

and m is the molecular mass and k is the Boltzmann constant. If the target is shaped as a square plate (s is the side of the square and A ¼ s2 ), then CD in Eq. (10) changes to

!     pffiffiffiffi 4 exp W2 1 erf ðWÞ : C D ¼ pffiffiffiffi þ p 1þ W p 2W2

ð13Þ

Eqs. (10)–(13) are discussed by, e.g., Gombosi (1994). To illustrate the quality of the numerically calculated drag forces, Fig. 8 also shows Eq. (10) evaluated for vd = 0 m s1 (squares) and vd = 300 m s1 (circles). The agreement is excellent at both small and large grain speeds. A numerical routine to calculate drag forces for flat targets has also been verified by comparing with Eq. (13). Even if Eqs. (10)–(13) formally loose validity within the Knudsen layer, it is useful to investigate the discrepancy between those formulae (evaluated for {n, Ttr, W} of the non-equilibrium gas) and the numerically calculated drag force. For the current simulations, the analytical expressions still perform very well, with discrepancies on the order of 5% (this is true also at 3.5 AU). Usage of Eqs. (10)–(13) therefore seem possible, even near the comet surface. However, it remains to be demonstrated that Eqs. (10)–(13) are applicable also for more extreme conditions, e.g., when the dust mantle is very thin. In that case, a large fraction of the molecules may escape directly from the ice at comparably low speed, while the rest interacts with the hot dust and come off at comparably high speed. The resulting mixed initial velocity distribution func-

tion is likely to be even more out of equilibrium than in the currently considered simulations, and there is no guarantee that the analytical formulae continue to perform well. Model V uses the numerical drag forces Fdrag = Fdrag(R, a, vd) and nucleus gravity (assuming qnuc = 500 kg m3) when integrating the equation of motion. The results for grains with a = 1 lm and three different densities (q = 500, 1000, or 1500 kg m3) are seen in the lower panel of Fig. 8. The terminal velocity of micron-sized q = 1000 kg m3 grains at rh = 1.3 AU is obtained as v 1 d ¼ 150 m s1 , which is almost an order of magnitude lower than the terminal gas velocity. Table 2 shows terminal grain velocities for different heliocentric distances, H2O:CO production rate ratios, grain densities, and grain sizes. Note that these velocities are lower than obtained when assuming the presence of exposed surface ice under similar illumination conditions. The drag force acting on vd = 0 m s1 grains in the DSMC cell closest to the nucleus can be used to estimate the size of the largest grain that can be lifted off the nucleus. Table 3 shows the maximum liftable grain size for various heliocentric distances, H2O:CO production rate ratios, and dust grain densities. Near perihelion, grains a few centimeters in radius can be lifted, while pure water comae at 2.5 AU and 3.5 AU shifts the maximum radius to the millimeter and micrometer range, respectively. Strong CO production at 3.5 AU can increase the maximum grain size in the coma to a few tens of a millimeter. Again, these results depend on assumptions regarding the location of sublimating ice (e.g., exposed or present underneath a dust mantle). 3.6. Drag forces on the rosetta spacecraft The gas drag force on the Rosetta spacecraft has been calculated by Fertig (2009), for various heliocentric distances and gas production rates, assuming a constant gas expansion velocity of

466

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

Table 2 Terminal dust grain velocities for various heliocentric distances and H2O:CO production ratios, for three different dust densities for each considered grain radius a. Presence of a global nucleus dust mantle is assumed, having a thickness fine-tuned to reproduce perihelion water production rates observed for Comet 67P/Churyumov–Gerasimenko. a = 1 lm

Model (rh, H2O:CO)

a = 10 lm

Terminal dust grain velocities (m s1) for (q = 500, 1000, 1500 kg m3) 1.3 AU, 1:0 200 150 120 70 2.5 AU, 1:0 31 22 18 10 3.5 AU, 1:0 1.6 0.9 0.6 – 3.5 AU, 1:1 2.6 1.7 1.3 – 3.5 AU, 1:65 17 12 9.7 5.3

a = 0.1 mm 50 7 – – 3.7

Table 3 Maximum liftable dust grain size for various heliocentric distances, H2O:CO production ratios and grain densities q, assuming a global nucleus dust mantle with a thickness fine-tuned to reproduce perihelion water production rates observed for Comet 67P/Churyumov–Gerasimenko. It is assumed that Rnuc = 1700 m and qnuc = 500 kg m3. Model (rh, H2O:CO)

q = 500 kg m3

Maximum liftable grain size 1.3 AU, 1:0 5.9 cm 2.5 AU, 1:0 1.1 mm 3.5 AU, 1:0 4.0 lm 3.5 AU, 1:1 8.9 lm 3.5 AU, 1:65 0.3 mm

q = 1000 kg m3

q = 1500 kg m3

3.0 cm 0.6 mm 2.0 lm 4.5 lm 0.2 mm

2.0 cm 0.4 mm 1.3 lm 3.0 lm 0.1 mm

W = 800 m s1 throughout the coma. The spacecraft is treated as a stationary object with mass mS/C = 1500 kg and surface area A ¼ 70 m2 , while the nucleus has Rnuc = 2400 m and either qnuc = 200 kg m3 or qnuc = 1000 kg m3. Fig. 9 shows the drag-togravity ratio Fdrag/Fgrav calculated by Fertig (2009) as thin curves. Note that Fdrag/Fgrav becomes independent on cometocentric distance when a constant W is assumed, since the coma gas number density changes as R2, just as gravity.

41 6 – – 3.0

23 3 – – 1.5

a = 1 mm 16 2 – – 0.8

13 2 – – 0.4

7 0.5 – – –

5 – – – –

4 – – – –

Drag forces on Rosetta have also been estimated using the current coma model. Specifically, numerical molecular distribution functions sampled during DSMC simulations have been used to calculate Fdrag for a target shaped as a plate with area 70 m2 (it has been verified that the numerical integration reproduces Eq. (13) when the coma is in the hydrodynamical regime). Rather high densities qnuc for the comparably small R = 1700 m nucleus have been assumed to match the same nucleus masses used by Fertig (2009). Calculated Fdrag/Fgrav ratios have been increased by factors 1.5, 1.4 and 2.0 at 1.3 AU, 2.5 AU, and 3.5 AU, respectively, to compensate for somewhat lower gas production rates in current modeling compared to those assumed by Fertig (2009). The results are seen in Fig. 9 as thick curves. The results are qualitatively different from those of Fertig (2009) due to a decrease of Fdrag/Fgrav with decreasing cometocentric distance. The drag force compared to the local gravity therefore becomes weaker near the nucleus (at rh = 3.5 AU, Fdrag/Fgrav is 63% higher at R = 100 km compared to R  Rnuc). Such a behavior of Fdrag/Fgrav could be relevant, i.e. for the dynamics of the Rosetta lander probe Philae during descent towards the nucleus surface. The quantitative differences between the current calculations and those of Fertig (2009) are not large, but the present coma mod-

Fdrag/Fgrav

10 Low mass x 1.5 Fertig (low Mnuc) 5

1.3 AU

drag

/Fgrav

0 1 10

F

High mass x 1.5 Fertig (high Mnuc)

0.2 0.15 0.1

2

10

3

10

4

10

5

10

Low mass x 1.4 Fertig (low Mnuc) High mass x 1.4 Fertig (high Mnuc)

2.5 AU

0.05 0 1 10

2

10

3

10

4

10

5

10

0.1

0.05

3.5 AU, H2O:CO=1:65

High mass x 2.0 Fertig (high Mnuc)

F

drag

/F

grav

Low mass x 2.0 Fertig (low Mnuc)

0 1 10

2

10

3

10

4

10

5

10

Height above surface [m] Fig. 9. The drag-to-gravity force ratio Fdrag/Fgrav according to Fertig (2009) (thin lines) and the present investigation (thick curves), at three different heliocentric distances as indicated in the panels. The ‘‘low” and ‘‘high” nucleus mass cases refer to nucleus bulk densities of qnuc = 200 kg m3 and qnuc = 1000 kg m3, respectively, for an assumed spherical nucleus of radius Rnuc = 2400 m. Note that our force ratios have been boosted by factors 1.4–2.0 (as indicated in the panels), in order to compensate as well as possible for differences in assumptions regarding the gas production rate in the two investigations.

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

el consistently yields lower Fdrag/Fgrav ratios. Particularly at large heliocentric distances, where the assumed high gas expansion velocity of W = 800 m s1 is increasingly at odds with simulations, the current Fdrag/Fgrav ratios are typically 2–3 times lower than those of Fertig (2009).

4. Discussion A number of gas kinetic models of comet water and carbon monoxide comae at various heliocentric distances have been presented. Although certain aspects of the modeling is idealized and thereby unrealistic (e.g., spherically symmetric outgassing, negligence of ionization, dissociation and chemical reactions), other processes are treated rather realistically. Special attention has been given to the thermophysics of the nucleus, the gas diffusion through the dust mantle, the evolution of the coma gas under non-equilibrium conditions, and the evaluation of drag forces in such comae. It has been demonstrated that the velocity distribution function deviates from the ideal Maxwellian, especially at rh P 2.5 AU. Macroscopic manifestations of such deviations include differences between parallel and perpendicular temperature components (with respect to the gas drift vector), as well as variations in the FWHM of Doppler broadened line profiles with viewing geometry, in addition to the normal change in the line central wavelength. Such variations were demonstrated for isolated parcels of gas. In reality, isolated parcels of gas are not observed, and the empirical line shape is instead determined by the gas properties along an entire line of sight. To reconstruct coma properties from Rosetta/MIRO line profiles by inverse modeling will therefore be very difficult. Not only will variations in {n, Ttr, Trot, W} along the line of sight play a rôle, but local deviations from the Maxwell– Boltzmann velocity distribution may complicate the problem further. It is therefore extremely important that the forward problem can be modeled, in order to gain experience on line shape behavior in various conditions, and to identify suitable approaches for solving the inverse problem. The merger of the current gas kinetic models with a radiative transfer model presented by Lee et al. (2009) is a crucial first step in that direction. An important application of such a unified model is to develop observational strategies that facilitates the solution to the inverse problem. It is quite clear that snapshot observations will not be sufficient for 3D reconstruction – e.g., one may have to observe the same region of the coma from a variety of geometries, in order to collect enough information to arrive at unambiguous solutions. To our knowledge, we have presented the first simulations of drag forces on dust grains and the Rosetta spacecraft, that do not assume that the velocity distribution of the gas is Maxwellian. Currently, such an increase in modeling realism may have little value. For example, we do not know the real nucleus mass of Comet 67P/ Churyumov–Gerasimenko or its gas production rate at rh J 2 AU and such uncertainties will always dominate over any second-order effect due to differences in modeling approach. At this point, it is not possible to say whether the drag-to-gravity ratios presented in Fig. 9 will correspond well to the actual conditions, or be off by an order of magnitude or more. However, once Rosetta in situ observations commences, many of the currently free parameters will be constrained accurately. Then, one may find that dust grain and spacecraft dynamics cannot be properly understood unless the assumption of Maxwellian gas is abandoned, and we believe it is important to prepare for such scenarios also at an early stage. Direct quantitative comparisons between our work and previously published DSMC coma models are difficult due to differences in assumptions about total and local gas production rates as well as nucleus size or shape. For the 1.3 AU perihelion case, the closest

467

match is provided by Tenishev et al. (2008) who studied a slightly larger nucleus (Rnuc = 1.98 km) with a three times smaller total production rate Q (however, the noon local production rate Z is only 20% lower than ours). The most interesting difference between the models concerns the terminal expansion velocity for water (along the comet–Sun axis), W1  750 m s1 according to Tenishev et al. (2008) compared to W1  1000 m s1 in the current work. A similar behavior is seen for our 2.5 AU case and their 2.7 AU case, for which local noon productions rates Z differ by only 10% (although their total production rate Q is 4.6 times smaller than ours). Here, the Tenishev et al. (2008) terminal velocity is W1  600 m s1 while we obtain W1  750 m s1. Note that W1 is not very sensitive to Q (for example, reducing Q a factor 40 when going from 1.3 AU to 2.5 AU only reduces W1 by 250 m s1). Hence, it is likely that our consistently larger terminal velocities for similar Q, compared to Tenishev et al. (2008), is due to our higher nucleus surface temperatures (295 K versus 196 K at perihelion and 201 K versus 186 K around 2.5 AU). The relatively high initial velocities of molecules exiting a hot dust mantle, compared to the slower molecules originating from exposed surface ice, may therefore result in measurably different terminal gas expansion velocities. Crifo et al. (2002, 2003) also present models with Q within a factor 2 of our 1.3 AU production rate, but only for nuclei that are substantially larger than ours (Rnuc = 5.0 km or 24.54 km). Large differences in nucleus size makes meaningful comparisons between models difficult, as illustrated by cross-checking our 2.5 AU model against model ‘‘205” by Zakharov et al. (2008a) which has a Q-value only 11% below ours but considers a substantially larger nucleus (Rnuc = 6.95 km). In our model, the gas obtains the expansion velocity W = 600 m s1 at a height of 2 km (where we have n  3.4  1015 molec m3), while model 205 reaches the same value at a height of 7 km along the comet–Sun axis (at that point, their n is 5 times smaller than ours). At 3.5 AU our pure water case has Q H2 O  1:5  1024 molec s1 . The models presented by Crifo et al. (2002) and Zakharov et al. (2008a) all have production rates that are at least an order of magnitude higher. Crifo et al. (2003) performed calculations for Q being a factor 3 smaller than ours, but only show the number density within 2 km of the irregularly shaped nucleus, which is substantially smaller (nucleus radius 0.5 km) than in the current study. The only DSMC calculation we are aware of, performed for a similar Q H2 O and nucleus size, showing all parameters of importance up to a large distance, was made by Tenishev et al. (2008). Their 3.25 AU case for water appears to match our simulation rather closely. Tenishev et al. (2008) also considers presence of carbon monoxide, but in proportion H2O:CO = 1:0.05 while we have studied H2O:CO = 1:1 or 1:65. The calculations regarding CO are therefore not directly comparable. As far as we know, only Crifo et al. (2004) have performed DSMC simulations for CO-dominated comae, considering a 2 km nucleus at 3 AU, assuming a carbon monoxide production rate of Q = 1027 molec s1 (i.e., one order of magnitude larger than our highest carbon monoxide production rate). At a distance of 20 km from the nucleus, they obtain Ttr  30 K and an expansion velocity W somewhat above 700 m s1. We obtain roughly the same temperature, but a much lower expansion velocity (W  470 m s1), due to our smaller QCO. Empirical estimates of water expansion velocities are available from e.g. analysis of the 18 cm hydroxyl (OH) line profile observed by radio telescope observations (Bockelée-Morvan et al., 1990; Tseng et al., 2007). These estimates are restricted to fairly high production rates (Q J 1028 molec s1) and depend somewhat on assumptions regarding, e.g. coma symmetry, OH velocity with respect to its parent, and the telescope beam size. Generally, beams tend to be large (5  103–5  105 km on target) compared to the currently simulated regions (3  102–103 km). Such large beams

468

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

may include gas which has been accelerated by photochemical heating (see e.g. Fig. 1 in Combi et al. (1999), where the plateau corresponding to the terminal velocity of pristine gas is followed by an expansion velocity increase when photodissociation sets in). Such photolytic heating is ignored in the current work due to its limited importance in the regions under consideration. Empirical expansion velocities (Vobs) would therefore tend to be higher than the pre-photolytic terminal velocities (V1), which should be compared to the W1-values presented in the current work. Expansion velocity measurements performed for conditions similar to our 1.3 AU model (1.3 6 rh 6 1.5 AU, 1.0  1028 6 Q 6 2.0  1028 molec s1) presented by Tseng et al. (2007) show substantial scatter, taking values in the range 0.5 6 Vobs 6 1.1 km s1. If photolytic heating is insignificant in these measurements (V1  Vobs), then our W1  1 km s1 value is consistent with the higher observed values. However, our terminal expansion velocity is several 0.1 km s1 higher than the lowest Vobs, a discrepancy that would worsen if photolytic heating indeed is a factor. Low empirical expansion velocities could indicate that outgassing models which assume presence of exposed surface ice are more realistic than the currently considered model, where sublimation takes place under a hot dust mantle. Cooler surfaces yield lower expansion velocities. However, spacecraft flybys show that at least some comets, e.g. Comet 19P/Borrelly (Soderblom et al., 2003) and Comet 9P/Tempel 1 (Groussin et al., 2007; Davidsson et al., 2000) have surfaces completely dominated by hot dust mantles. Furthermore, the reduction of the terminal expansion velocity corresponding to cooler surface conditions, e.g. W1 = 0.75 km s1 (Tenishev et al., 2008), may not be sufficient to explain some very low Vobs measurements. Bockelée-Morvan et al. (1990) suggested a number of factors apart from low surface temperatures that could explain low empirical expansion velocities. For example, the assumed OH ejection velocity could be too high, coma asymmetry could be a problem, heavy molecules could be more common than previously assumed. An additional, and perhaps more important, factor is that most models ignore feedback effects that the dust has on gas. Many models pre-calculate the gas environment, then study the dynamics of dust grains in such gas. However, in reality a potentially large amount of dust (comparable to the gas mass) must be accelerated from rest to the dust terminal velocity. In this process, the gas looses momentum and energy which should result in a lower expansion velocity. Multifluid hydrodynamics calculations where such feedback of dust on the gas is considered (e.g., Probstein, 1968; Gombosi et al., 1985) confirm that expansion velocities could be lowered measurably. For this mechanism to be effective it may be necessary that gas is not heated during diffuse scattering, i.e., dust grains may need to remain comparably cold (Tdust  Tsurf) until the gas–dust coupling essentially is over. The 1 lm grains at 1.3 AU decouple dynamically from the gas about 5 km above the surface, a distance they cross in about a minute. Hence, if substantial grain heating can be avoided within the first few minutes after entering the coma, the net effect of dust may indeed be to slow down the gas flow. The observed dispersion in measured expansion velocities could simply be a consequence of variations in telescope beam size. Relatively large beams would sample gas that has been accelerated by photolytic heating and yield larger Vobs, while smaller beams may recover expansion velocities before substantial photodissociation has occurred. However, differences in gas-to-dust mass ratio and surface coverage of mantles versus ice among comets could also contribute to the dispersion. The measurements by MIRO and other instruments on the Rosetta spacecraft will definitively shed light on these problems, and push model development and data analysis in the right direction.

Acknowledgments Davidsson most gratefully thanks the US Rosetta Project Manager Dr. Claudia Alexander and the Rosetta/MIRO Principal Investigator Dr. Samuel Gulkis for the invitation to work at the Jet Propulsion Laboratory in Pasadena (CA) for six months as a visiting scientist. A part of this research was carried out at the Jet Propulsion Laboratory (JPL), California Institute of Technology, under a contract with the National Aeronautics and Space Administration. This work was funded by the NASA-Rosetta project. Appendix A The governing equation for model I (thermophysical nucleus modeling) is

qc

@T @2T ¼ j 2 þ S; @t @x

ðA:1Þ

solved on a domain x 2 [0, xmax] with one set of material parameters {q, c, j} for 0 6 x < xd (dust mantle) and another for xd 6 x 6 xmax (icy interior), as specified in Section 2. The term S governs energy consumption due to net sublimation of ice and take non-zero values only at x 2 [xd, xd + Dx] according to

S ¼ Z HK ðTÞas ðTÞLðTÞ

pesc ; PDx

ðA:2Þ

where the Hertz–Knudsen formula (Fanale and Salvail, 1984; Cercignani, 2000) is

Z HK ¼ 3:56  1012  rffiffiffiffiffiffiffiffiffiffiffiffi 6141:667 m  exp ðkg m2 s1 Þ; T 2pkT

ðA:3Þ

the sublimation coefficient (Kossacki et al., 1999) is

asub ðTÞ ¼ 0:573 þ 0:427    T  150:5 p ;  tanh 4:353 tan p  122:5 2

ðA:4Þ

the latent heat (Tancredi et al., 1994) is 1

LðTÞ ¼ 2:888  106  1116T ðJ kg Þ;

ðA:5Þ

and the escape probability pesc (see Section 2) is corrected for multiple scattering by

PðTÞ ¼ 1  ð1  acond ðTÞÞð1  pesc Þ;

ðA:6Þ

where the condensation coefficient is taken as acond = asub. The upper boundary condition to Eq. (A.1) is given by

    S ð1  AÞ 1 P cosðd Þ cosðl Þ þ sinðd Þ sinðl Þ cos x t  co co co co 2 r 2h   @T ¼ erSB T 4  j  ; @x x¼0

ðA:7Þ 2

1

where the solar constant is S ¼ 1367 J m s , the albedo is A = 0.04, co-declination is dco = 90° (i.e., spin axis assumed perpendicular to orbital plane), co-latitude lco = 90° (i.e., equatorial conditions), unless lco is evaluated for many latitudinal slabs when estimating the total gas production of a sphere, x ¼ 2Pp, where the rotational period is taken as P = 12.3 h (Lamy et al., 2003, 2007), the emissivity is e = 1.0, and rSB is the Stefan–Boltzmann constant. The lower boundary condition is taken as

 @T  ¼ 0: @x x¼xmax

ðA:8Þ

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

Eq. (A.1), along with boundary conditions, is solved by using the Finite Element Method with spatial grid cells of Dx = 5  103 m (xmax = 0.4 m yields 80 cells), time discretization Dt = 30 s and typically performing 10 nucleus rotations to reach steady-state. Appendix B The Monte Carlo model used to calculate the molecular escape probability pesc and the molecular distribution function of molecules which successfully has percolated through a porous dust mantle (model II) is a variant of the model presented by Davidsson and Skorov (2004). Test particles are released at the ice/dust interface (x = xd) with a direction of motion specified by the angle h with respect to the outward surface normal and azimuth angle u according to



h ¼ arccosð1  2n1 Þ;

u ¼ 2pn2 ;

ðB:1Þ

where n1, n2 2 [0, 1] are normally distributed real random numbers. Molecular speeds V are drawn from the Maxwellian transmission distribution, integrated over all angles, thus obtaining the speed distribution,

g t ¼ 2b4 V 3 expðb2 V 2 Þ;

ðB:2Þ

where b, defined by Eq. (12) is evaluated for the local temperature. Selecting a speed for an individual test particle that statistically is consistent with Eq. (B.2) is made with the acceptance–rejection method (e.g., Bird, 1994). Here, V* is generated randomly, and if gt(V*) 6 n3, where n3 2 [0, 1] is another random real number, then V = V* (otherwise, another {V*, n3} pair is tested). The distance k to the next point of interaction with the dust mantle is drawn from a distribution of free paths K given by

ln K ¼

6 X

ci ðln kÞi ;

ðB:3Þ

The time-development of the vapor is modeled by alternating between a translation phase and a collision phase. During the translation phase (of duration Dt), the test particles are moved linearly in accordance with their velocity vectors (here, in the absence of a gravitational field, due to the low mass of a cometary nucleus). During the translation phase, no intermolecular collisions are considered, but if a boundary is encountered, this is accounted for (diffuse reflection at the lower boundary, escape through the outer boundary). During the collision phase, the test particles are considered stationary, while their velocity vectors and rotational energies are updated, as a result of collisions. In order to treat collisions, the considered volume is divided into Nc spatial cells. For each cell and collisional phase, the number of potential collision partners is calculated,

Npcp ¼

Ntp;i hNtp;j iðV rel rÞij F N Dt ; 2V c

Appendix C The variant of the DSMC technique utilized here (model III) to calculate the gas evolution throughout the coma is that of Bird (1994), summarized briefly in the following. In the DSMC framework, the vapor is modeled at a kinetic (i.e., molecular) level by considering a large number of Monte Carlo test particles, each representing a certain number FN of real molecules. Each test particle is characterized by its position coordinates {x, y, z}, its velocity vector {u, v, w}, its rotational energy Erot and species-related properties (mass, number of degrees of rotational freedom). Note that z = R, the height above the sublimating surface, is the only spatial coordinate of importance in the present spherically symmetric 1D problem.

ðC:1Þ

where Ntp,i is the instantaneous number of test particles of species i in the cell, hNtp,ji is the number of test particles of species j in the cell averaged over time, Vc is the cell volume, Vrel is the relative velocity, and r is the collision cross section ((Vrelr)ij is the largest value of the product Vrelr ever encountered for individual particle pairs {i, j} in the cell during the course of simulation). The 2Npcp test particles (generally an overestimate) are then paired up for potential collisions. Whether a collision between such a pair of potential colliders actually takes place, is subject to random selection, according to the probability function Vrelr/ (Vrelr)ij, evaluated for the pair in question. The collision cross section of a test particle, and the outcome of a collision (i.e., its effect on the velocity vectors and rotational energies of colliding particles), depend on the applied collision model. In this work, the Variable Hard Sphere (VHS) model is applied, using molecular diameters that are functions of relative velocity in order to mimic a viscosity dependence on temperature of the form l / Tx, where x is a constant. The molecular diameter d is given by,

i¼0

with coefficients c0 = 4.179, c1 = 0.24791, c2 = 0.41939, c3 = 1.5275  104, c4 = 4.9138  103, c5 = 9.3158  103, c6 = 1.6195  103. Eq. (B.3) is an analytical fit to data on the distribution of free paths sampled by having a test particle bouncing around within a computer-generated medium with porosity w = 0.7, consisting of randomly placed spheres. Fig. 2 of Davidsson and Skorov (2004) shows a plot of Eq. (B.3). No molecular recondensation takes place within the dust mantle. If migrating molecules hit the ice/dust interface they either recondense or scatter according to the condensation coefficient acond = asub (Eq. (A.4)). Molecules which escape the medium (i.e., reach x < 0) have their {V, h} values sampled to build the distribution function of molecules exiting the nucleus.

469

d ¼ dref

 x12 hV rel;ref i ; hV rel i

ðC:2Þ

where dref is the molecular diameter at a reference temperature Tref and hVreli is the mean relative velocity,

sffiffiffiffiffiffiffiffiffi 8kT ; hV rel i ¼ pmr

ðC:3Þ

where the reduced mass is mr = mimj/(mi + mj) for species i and j. The mean relative velocity at the reference temperature Tref is denoted hVrel,refi. The VHS reference diameter is given by

dref

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 15 mkT ref =p ; ¼ 2ð5  2xÞð7  2xÞlref

ðC:4Þ

where lref = l(Tref) is the viscosity at the reference temperature. Given an empirical x, and the values of {hVreli, l} at an arbitrary Tref, the VHS molecular diameter can therefore be calculated using Eqs. (C.2)–(C.4), and thereby, the collision cross section. For water we use x = 1.108, dref = 5.6159  1010 m and m = 2.99  1026 kg at Tref = 373.16 K (Davidsson and Skorov, 2004). For carbon monoxide we use x = 0.73, dref = 4.19  1010 m and m = 4.65  1026 kg for Tref = 273 K (Bird, 1994). For H2O–CO collisions we average the x and dref values for the individual species. The total cross section is given by r = p (di + dj)2/4, where di,j is the molecular diameter of species i or j. If a VHS collision is elastic, energy and momentum conservation leads to an (on average) isotropic reorientation of the relative velocity vector V rel (its magnitude jV rel j does not change).

470

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

Inelastic collisions are treated according to the phenomenological Larsen–Borgnakke model. For a particular collision, the relative translational energy Etrel ¼ mr V 2rel =2, is first added to the sum of internal energies for the two particles Eint ¼ Erot1 þ Erot2 , thereby obtaining the total energy, Etot = Etrel + Eint. When redistributing Etot between the translational and rotational modes (as a result of the collision), it is required that the post-collision values Etrel and Eint ¼ Etot  Etrel are selected from the distribution functions

    3x E g E / Etrel 2 exp  trel ; trel kT      f1 Eint exp  g E / Eint ; int kT

ðC:5Þ ðC:6Þ

where f = (fi + fj)/2 and {fi, fj} are the number of internal degrees of freedom for species i and j, while the T is proportional temperature

to Etot. The probability that a pair Etrel ; Eint fulfills these distribution functions simultaneously is given by the product of Eqs. (C.5) and (C.6), normalized to yield a peak value of unity,

  P1 Etrel ¼

 !32x f  x þ 12 Etrel 3 Etot x 2   f1 f  x þ 12 E  : 1  trel f1 Etot

ðC:7Þ

The practical problem of finding a single value Etrel , that statistically conforms with the probability function P1 ðEtrel Þ is here solved by applying the acceptance–rejection method. Here, a new translational energy Etrel is proposed by randomly selecting a value in the interval [0, Etot] (leading to a potential new rotational energy Eint ¼ Etot  Etrel ). Eq. (C.7) is evaluated for this Etrel , and in case R 6 P 1 ðEtrel Þ is fulfilled, the suggested value is accepted (R 2 ½0; 1 is a random number). Otherwise it is rejected, and the procedure is repeated. The new rotational energy Eint must be distributed between the two test particles. This is done by applying the acceptance–rejection method, to find Erot1 from the distribution

P2 ðErot1 Þ ¼ 2f2

  2f 1   f 1 Erot1 Erot 2 1  1 :  Eint Eint

ðC:8Þ

It is noted that P2 peaks for energy equipartition between the test particles. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The new relative speed is given by V rel ¼ 2Etrel =mr , and as was the case for elastic collisions, a simple isotropic reorientation of the relative velocity vector is performed. A finite relaxation rate for energy equipartition between internal and kinetic energy, is introduced by only considering a fraction Kc ¼ 15 of the collisions as inelastic, while the remainder are treated as elastic (Kc is the relaxation collision number). Macroscopic gas quantities are calculated per cell and species, with instantaneous values (for a cell containing Ntp test particles) obtained by first evaluating the following quantities:

fV ðxÞ ; V ðyÞ ; V ðzÞ g ¼

N tp 1 X fui ; Ntp i¼1

vi;

N tp n o 2 1 X ðxÞ ðyÞ ðzÞ V2 ; V2 ; V2 ¼ u ; Ntp i¼1 i

wi g;

v 2i ;

w2i ;

ðC:9Þ ðC:10Þ

Ntp;i F N ; Vc  m ðlÞ Tl ¼ V 2  ðV ðlÞ Þ2 ; k 1 T tr ¼ ðT x þ T y þ T z Þ; 3 2erot : T rot ¼ fi kN tp;i

n¼

ðC:11Þ ðC:12Þ ðC:13Þ ðC:14Þ

Finally, the expansion (or bulk drift) velocity is given by W = V(z). These quantities are prone to fluctuate violently with time, since Ntp generally is a very small, constantly changing number. Therefore, the flow properties are sampled a large number of times during the course of simulation, i.e., Eqs. (C.11)–(C.14) are ðlÞ based on accumulated values of fN tp ; V ðlÞ ; V 2 ; erot g. The original code of Bird (1994) has been modified in order to allow injection of test particles at the lower boundary in accordance with distribution functions sampled with model II. Also, diffuse scattering of coma molecules with the (hot) dust mantle surface has been implemented. Another important modification consists of sampling the entire distribution function within each cell (for each species), thereby documenting the cumulative number of test particles having a certain speed S and propagating at an b Speed bins of width angle h with respect to the outward normal R. 1 DS = 10 m s are used, with an upper cut-off at Smax = 1200 m s1 or Smax = 2000 m s1 (the maximum naturally occurring molecular speed depends on rh). The angular bins has widths of Dh of 1°–2°. If the sampled number of test particles is denoted by g0 (R, S, h), the distribution function g = g(R, S, h, u) is obtained by calculating the normalization factor

N¼

2pDSDh

XX S

!1 g 0 ðR; S; hÞ

;

ðC:15Þ

h

and evaluate

gðR; S; h; uÞ ¼

Ng 0 ðR; S; hÞ S2 sin h

:

ðC:16Þ

Cell sizes, time steps, and FN values are selected in accordance with general DSMC recommendations. The current simulations divides the computational region (1700 6 R 6 4100 m for the 1.3 AU case and 1700 6 R 6 3  105 m for all other cases) into 6000– 20,000 cells (selected so that a molecular mean free path k should be resolved by at least three cells at all points in the flow). The time step is typically Dt = 104–102 s, selected to make sure that test particles cross a cell in at least four time steps. Typically, 5.7  104–1.6  106 test particles are present inside the modeled region at any given moment, each representing 8.3  1021– 1.0  1024 molecules (aiming at having at least 10 test particles per cell). The sampling procedure has two phases – a relaxation phase and a cumulative sampling phase. During the relaxation phase, test particles gradually fill the initially empty calculational region, and a steady-state flow is achieved. The flow properties are sampled every 4–8 time step, and 1000 samples are performed between file dumps. Typically, 10–40 dumps are performed to monitor the gradual relaxation (e.g., at some point {n, Ttr, Trot, W} are no longer modified at a given height, apart from statistical scatter). At that point, cumulative sampling is initiated, during which the flow is sampled typically 6  104–2.1  105 times to achieve good statistics. Appendix D

and summing up the total rotational energy for all test particles in the cell, erot. The number density n, translational temperature components Tl (with l = x, y, z), translational temperature Ttr, and rotational temperature Trot for a species i are then given by

In order to calculate the drag force on a spherical or flat target b (model IV), the molecular distribution functions with velocity v d R g = g(R, S, h, u) obtained from model III sampling are first trans-

B.J.R. Davidsson et al. / Icarus 210 (2010) 455–471

b For formed into a coordinate frame which is co-moving at v d R. each {S, h, u} combination, the corresponding components {S, h, u} in the co-moving system are given by

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 2 2 > > < S ¼ S  2Sv d cos h þ v d ;  S cos hv d ; > S > h ¼ arccos : u ¼ u;

ðD:1Þ

and a distribution function for the co-moving frame is built as G(R, S, h, u) = g(R, S, h, u). In practice, G is integrated over NS speed bins of size DS, Nh angular bins of size Dh and Nu azimuthal bins of thickness Du, which may or may not be identical to the discretization used during model III sampling. The resulting quantity is denoted G0 ðR; S ; h ; u Þ ¼ GðR; S ; h ; u ÞS2 sinðh ÞDS Dh Du . Targets (spheres of radius a = 1 m or plates with square side a = 1 m) are built by dividing their surfaces into Na locally planar facets, with outward surface normals N i and surface area ai. The drag force due to momentum transfer by specularly reflected molecules is then calculated as,

F drag ðR;

v d Þ ¼ mnðRÞ

Ni X

ai N i

i¼1

 

Nu NS X Nh X X 

!

2 0 min 0; 2Ni  V jkl G ðR; Sj ; hk ; ul Þ ; ðD:2Þ

j¼1 k¼1 l¼1

where indices {j, k, l} point at S-, h-, and u-bins, respectively. Note that V jkl ¼ fSj sin hk cos ul ; Sj sin hk sin ul ; Sj cos hk g and that b due to symmetry. Scaling to other target sizes are readily F drag k R done as Fdrag(R, a, vd) = a2Fdrag(R, vd). Eq. (D.2) has to be evaluated for a large number of vd values for each height R (typically, 10 m s1 resolution for 0 6 vd 6 100 m s1 and 50 m s1 resolution for vd P 150 m s1). Model V is used to calculate the dynamics of dust grains in the comet coma. The governing equation is the equation of motion,

4 dv pqa3 d ¼ F drag ðR; a; 3 dt

vdÞ 

4pqa3 GMnuc 3ðR þ Rnuc Þ2

;

ðD:3Þ

with initial conditions R(t = 0) = Rnuc and vd(Rnuc) = 0 m s1. The solutions R = R(t) and vd = vd(t) are obtained by using the ordinary differential equation solver ‘‘ode45” in MATLAB. References Bensch, F., Bergin, E.A., 2004. The pure rotational line emission of ortho-water vapor in comets. I. Radiative transfer model. Astrophys. J. 61, 531–544. Bird, G.A., 1994. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Oxford University Press Inc., New York. Biver, N., and 22 colleagues, 2002. The 1995–2002 long-term monitoring of Comet C/1995 O1 (Hale-Bopp) at radio wavelengths. Earth Moon Planets 90, 5–14. Bockelée-Morvan, D., Crovisier, J., Gérard, E., 1990. Retrieving the coma gas expansion velocity in P/Halley, Wilson (1987 VII) and several other comets from the 18-cm line shapes. Astron. Astrophys. 238, 382–400. Bockelée-Morvan, D., Moreno, R., Biver, N., Crovisier, J., Crifo, J.-F., Fulle, M., Grewing, M., 2004. CO and dust productions in 67P/Churyumov–Gerasimenko at 3 AU post-perihelion. In: Colangeli, L., Epifani, E.M., Palumbo, P. (Eds.), The New Rosetta Targets. Observations, Simulations and Instrument Performances, Astrophysics and Space Science Library, vol. 311. Kluwer Academic Publishers, Dordrecht, pp. 25–36. Crovisier, J., 1984. The water molecule in comets: Flourescence mechanisms and thermodynamics of the inner coma. Astron. Astrophys. 130, 361–372. Cercignani, C., 2000. Rarefied Gas Dynamics. From Basic Concepts to Actual Calculations. Cambridge University Press, Cambridge. Combi, M.R., Kabin, K., DeZeeuw, D.L., Gombosi, T.I., Powell, K.G., 1999. Dust–gas interrelations in comets: Observations and theory. Earth Moon Planets 79, 275– 306. Crifo, J.F., Lukianov, G.A., Rodionov, A.V., Khanlarov, G.O., Zakharov, V.V., 2002. Comparison between Navier–Stokes and direct Monte-Carlo simulations of the cicumnuclear coma: I. Homogeneous, spherical source. Icarus 156, 249–268.

471

Crifo, J.F., Loukianov, G.A., Rodionov, A.V., Zakharov, V.V., 2003. Navier–Stokes and direct Monte Carlo simulations of the cicumnuclear coma: II. Homogeneous, aspherical sources. Icarus 163, 479–503. Crifo, J.-F., Lukyanov, G.A., Zakharov, V.V., Rodionov, A.V., 2004. Physical model of the coma of Comet 67P/Churyumov–Gerasimenko. In: Colangeli, L., Epifani, E.M., Palumbo, P. (Eds.), The New Rosetta Targets. Observations, Simulations and Instrument Performances, Astrophysics and Space Science Library, vol. 311. Kluwer Academic Publishers, Dordrecht, pp. 119–130. Davidsson, B.J.R., 2008. Comet Knudsen layers. Space Sci. Rev. 138, 207–223. Davidsson, B.J.R., Skorov, Y.V., 2004. A practical tool for simulating the presence of gas comae in thermophysical modeling of cometary nuclei. Icarus 168 (1), 163– 185. Davidsson, B.J.R., Gutiérrez, P.J., Rickman, H., 2000. Physical properties of morphological units on Comet 9P/Tempel 1 derived from near-IR Deep Impact spectra. Icarus 201, 335–357. Fanale, F.P., Salvail, J.R., 1984. An idealized short-period comet model: Surface insolation, H2O flux, dust flux and mantle evolution. Icarus 60, 476–511. Fertig, J., 2009. Observing scenarios during escort phase. Presentation during the 26th Rosetta Science Working Team Meeting, June 11–12, 2009, at ESTEC, Noordwijk, The Netherlands. Glassmeier, K.-H., Boehnhardt, H., Koschny, D., Kührt, E., Richter, I., 2007. The Rosetta mission: Flying towards the origin of the Solar System. Space Sci. Rev. 128, 1–21. Gombosi, T.I., 1994. Gaskinetic Theory. Cambridge University Press, Cambridge. Gombosi, T.I., Cravens, T.E., Nagy, A.F., 1985. Time-dependent dusty gasdynamical flow near cometary nuclei. Astrophys. J. 293, 328–341. Groussin, O., A’Hearn, M.F., Li, J.-Y., Thomas, P.C., Sunshine, J.M., Lisse, C.M., Meech, K.J., Farnham, T.L., Feaga, L.M., Delamere, W.A., 2007. Surface temperature of the nucleus of Comet 9P/Tempel 1. Icarus 187, 16–25. Gulkis, S., and 32 colleagues, 2009. MIRO: Microwave Instrument for Rosetta Orbiter. In: Schulz, R., Alexander, C., Boehnhardt, H., Glassmeier, K.-H. (Eds.), ROSETTA: ESA’s Mission to the Origin of the Solar System, Springer, New York, pp. 291–314. Huebner, W.F., Markiewicz, W.J., 2000. The temperature and bulk flow speed of a gas effusing or evaporating from a surface into a void after reestablishment of collisional equilibrium. Icarus 148, 594–596. Kossacki, K.J., Markiewicz, W.J., Skorov, Y., Kömle, N.I., 1999. Sublimation coefficient of water ice under simulated cometary-like conditions. Planet. Space Sci. 47 (12), 1521–1530. Lamy, P.L., Toth, I., Weaver, H., Jorda, L., Kaasalainen, M., 2003. The nucleus of Comet 67P/Churyumov–Gerasimenko, the new target of the Rosetta mission. Bull. Am. Astron. Soc. 35 (4), 970. Lamy, P.L., Toth, I., Davidsson, B.J.R., Groussin, O., Gutiérrez, P.J., Jorda, L., Kaasalainen, M., Lowry, S.C., 2007. A portrait of the nucleus of Comet 67P/ Churyumov–Gerasimenko. Space Sci. Rev. 128, 23–66. Lee, S., von Allmen, P., Kamp, L., Gulkis, S., Davidsson, B., 2009. Non-LTE radiative transfer calculations for a submillimeter water line in cometary comae observed by Rosetta/MIRO. Astron. Astrophys, submitted for publication. Osip, D.J., Schleicher, D.G., Millis, R.L., 1992. Comets: Groundbased observations of spacecraft mission candidates. Icarus 98, 115–124. Probstein, R.F., 1968. The dusty gasdynamics of comet heads. In: Problems of Hydrodynamics and Continuum Mechanics. Society for Industrial and Applied Mathematics, Philadelphia, pp. 568–583. Schleicher, D.G., 2006. Compositional and physical results for Rosetta’s new target Comet 67P/Churyumov–Gerasimenko from narrowband photometry and imaging. Icarus 181 (2), 442–457. Soderblom, L.A., Britt, D.T., Brown, R.H., Buratti, B.J., Kirk, R.L., Owen, T.C., Yelle, R.V., 2003. Short-wavelength infrared (1.3–2.6 lm) observations of the nucleus of Comet 19P/Borrelly. Icarus 167, 100–112. Sone, Y., Sugimoto, H., 1993. Kinetic theory analysis of steady evaporating flows from a spherical condensed phase into a vacuum. Phys. Fluids A 5 (6), 1491– 1511. Tancredi, G., Rickman, H., Greenberg, J.M., 1994. Thermochemistry of cometary nuclei. I: The Jupiter family case. Astron. Astrophys. 286, 659–682. Tenishev, V., Combi, M., Davidsson, B., 2008. A global kinetic model for cometary comae: The evolution of the coma of the Rosetta target Comet Churyumov– Gerasimenko throughout the mission. Astrophys. J. 685 (1), 659–677. Tseng, W.-L., Bockelée-Morvan, D., Crovisier, J., Colom, P., Ip, W.-H., 2007. Cometary water expansion velocity from OH line shapes. Astron. Astrophys. 467 (2), 729– 735. Wagner, W., 1992. A convergence proof for Bird’s direct simulation Monte Carlo method for the Boltzmann equation. J. Statist. Phys. 66 (3–4), 1011–1044. Zakharov, V.V., Bockelée-Morvan, D., Biver, N., Crovisier, J., Lecacheux, A., 2007. Radiative transfer simulation of water rotational excitation in comets. Comparison of the Monte Carlo and escape probability methods. Astron. Astrophys. 473, 303–310. Zakharov, V.V., Rodionov, A.V., Lukyanov, G.A., Crifo, J.F., 2008a. Navier–Stokes and direct Monte-Carlo simulations of the cicumnuclear coma III. Spherical, inhomogeneous sources. Icarus 194, 327–346. Zakharov, V.V., Bockelée-Morvan, D., Biver, N., Crovisier, J., Crifo, J.F., Gulkis, S., Rodionov, A.V., 2008b. Numerical simulations of water spectra obtained with the microwave instrument for the Rosetta orbiter (MIRO) from Comet 67P/ Churyumov–Gerasimenko. LPI Contributions 1405, 8144.