Issue 
A&A
Volume 619, November 2018



Article Number  A181  
Number of page(s)  12  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201833566  
Published online  26 November 2018 
Solving nonLTE problems in rotational transitions using the Gauss–Seidel method and its implementation in the Atmospheric Radiative Transfer Simulator
^{1} Terahertz Technology Research Center, National Institute of Information and Communications Technology, 421, NukuiKitamachi, Koganei, Tokyo, 1848795, Japan
email: ykasai@nict.go.jp
^{2} Department of Environmental Chemistry and Engineering, Tokyo Institute of Technology, G117, Nagatsuda, 4259, Midoriku, YokohamaCity, Kanagawa Prefecture, 2268502, Japan
^{3} Max Planck Institute for Solar System Research, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
Received:
4
June
2018
Accepted:
26
August
2018
This article presents our implementation of a nonLTE solver in spherical symmetry for molecular rotational transition in static or expanding atmospheres. The new opensource code relies on the Gauss–Seidel Accelerated Lambda Iteration methodology that provides a rapid and accurate convergence of the nonLTE problems, which is now routinely used in astrophysical and planetary research. The nonLTE code is interfaced with the widely used package, the Atmospheric Radiative Transfer Simulator (ARTS), to facilitate spectral line simulations for various viewing geometries. In this paper we describe the numerical implementation, provide the first validation results for the populations against two other nonLTE codes, and then discuss the possible application. The quantitative comparisons are performed using an established orthowater nonLTE model applied to cases of optical thick and thin conditions of Ganymede’s atmosphere. The differences in populations expressed as excitation temperatures show very good agreement in both cases. Finally, we also apply this model to a sample of data from the Microwave Instrument for the Rosetta Orbiter (MIRO) instrument. The new nonLTE package is demonstrated to be fast and accurate, and we hope that it will be a useful addition to the planetary community. In addition, being open source and part of the ARTS, it will be further improved and developed.
Key words: radiative transfer / methods: numerical / submillimeter: planetary systems / planets and satellites: atmospheres
© ESO 2018
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
The past two decades have seen numerous advances in heterodyne radio technologies and their application to modern themes in atmospheric science of planets, small bodies, and astrophysics. These instruments, operating typically in in the submillimetre wave range up to frequencies of 5 THz, provide a very high frequency resolution such that photons from individual rotational transitions of molecules can be measured. These measurements are currently accurate enough with high signaltonoise ratio such that Doppler shifts in line frequency can be used for wind determination along the line of sight. Perhaps, one of the most important lines is the orthowater molecule transition (1_{10} → 1_{01}) at 556.936 GHz, which has become crucial for learning about activities of cold icy objects in the solar system linked to our origins (BockeleeMorvan 1987; Bensch & Bergin 2004; Zakharov et al. 2007). Several observatories use this transition for detection, for example the Submillimeter Wave Astronomical Satellite (SWAS; Neufeld et al. 2000; Bensch et al. 2007), the Odin satellite equipped with a 1.1 m submillimetre telescope (Biver et al. 2007), and the Heterodyne Instrument for the Far Infrared (HIFI) on the Herschel Space Observatory (Hartogh et al. 2009; de Graauw et al. 2010). Similarly, the analogue transitions of water isotopologues ( and HDO) can be used to learn about the origin of water (Hartogh et al. 2011b; BockeléeMorvan et al. 2012; Biver et al. 2012; Lis et al. 2013). Naturally, molecules other than water can also be easily studied in this wavelength range (HCN, CH_{3}OH, H_{2}S, and CS; Biver et al. 2012; de ValBorro et al. 2012b, 2013).
More recently, the Rosetta mission carried the Microwave Instrument for the Rosetta Orbiter (MIRO; Gulkis et al. 2007) to comet 67P/ChuryumovGerasimenko, which for more than two years studied the near nucleus coma from rotational transitions of H_{2}O, NH_{3}, CO, and CH_{3}OH (Gulkis et al. 2015; Lee et al. 2015; Biver et al. 2015). Nevertheless, other solar system objects can be studied to detect and monitor the permanent and the transient atmospheres (Lellouch 2008). For instance, we have learned about the tenuous water atmosphere of Ceres from Herschel/HIFI (Küppers et al. 2014), and Saturn’s largest moon Enceladus and its water torus (Hartogh et al. 2011a).
In the future the exploration of the Jovian system as an archetype for gas giants and the characterization of the Ganymede, Europa, Callisto, and Io surface–atmosphere relationship will be taken up by the European JUpiter Icy Moon Explorer (JUICE) mission (Grasset et al. 2013). One of the relevant measurements will be made possible by the Submillimetre Wave Instrument (SWI) to profile the H_{2}O, winds, and temperature from emission of water lines at 557, 1153, and 1163 GHz in the atmospheres of moons^{1}.
A wellknown difficulty in interpreting these rotational transitions in the tenuous atmosphere of comets, moons, and asteroids stems from the fact that even the rotational populations cannot be assumed to be in local thermodynamic equilibrium (LTE), as noted in all the works cited in the previous paragraph. Furthermore, in the case of the lowest rotational transition of orthoH_{2}O a significant optical depth can exist even in such harsh conditions as a cometary coma (Zakharov et al. 2007; Lee et al. 2015). Therefore, accurate nonLTE molecular excitation codes have to be developed. An additional problem for these smaller objects is that their models cannot rely on the planeparallel assumption typically used for nonLTE in the rovibrational bands on the Earth and Mars (LópezPuertas & Taylor 2001). The solution of multilevel radiative transfer under nonLTE conditions in spherical symmetry can be considered a classical problem in astrophysics (Kunasz & Hummer 1974a,b; see also Anusha et al. 2009), and it is a demanding computational task even in simple 1D atmospheres. The challenging task is the repeated calculation of the mean radiative field at each atmospheric level required to solve the system of statistical equilibrium equations (SEE; see Sect. 2). A popular technique to simplify the calculation for expanding cometary atmospheres is the socalled escape probability (EP) method (Rybicki 1985; BockeleeMorvan 1987; Zakharov et al. 2007; de ValBorro et al. 2012a; O’Rourke et al. 2013). The advantage of EP is that the problem is treated as local (with photon escape probabilities calculated without radiative transfer), and that it is a deterministic method. Another popular method relies on the Monte Carlo (MC) method to obtain the mean intensity (Hogerheijde & Van Der Tak 2000; Zakharov et al. 2007; Lee et al. 2011), which can cope in principle with more complex geometry, but the calculated populations show intrinsic random noise. While both methods show good agreement in the forward simulations of cometary lines Zakharov et al. (2007), the MC method has a clear disadvantage if a routine inverse problem relying on a nonLTE Jacobian matrix has to be solved (Lee et al. 2015).
In this paper we describe the implementation of a stateoftheart deterministic method for the solution of multilevel nonLTE problem in a 1D spherically symmetric atmosphere able to accurately treat optically thick transitions in both static and expanding atmospheres. The method falls under the class of the Accelerated Lambda Iteration method (ALI; described in the next section), but relies on Gauss–Seidel approach for population corrections during iterations (Trujillo Bueno & Fabiani Bendicho 1995; Asensio Ramos & Trujillo Bueno 2006). The aim is to incorporate the implemented nonLTE model into the public version of the Atmospheric Radiative Transfer Simulator (ARTS; Buehler et al. 2017; Eriksson et al. 2011); the present version is made available via the ARTSassociated python package ‘typhon’.
This paper is set up as follows. Section 2 describes the formalism of radiative transfer in spherical geometry, and the solution of the multilevel nonLTE problem through the concepts of lambda iteration (LI), multilevel ALI (MALI), and multilevel Gauss–Seidel (MUGA). In Sect. 3 we compare the performance of nonLTE methods, and also validate our results under optically thick and thin conditions. Here we also present some forward calculations for the main water transitions in the Ganymede atmosphere predicted by the direct simulations model of Marconi (2007) and an example of additional application of our nonLTE model and spectra calculation in the context of cometary atmospheres as measured by the Rosetta/MIRO instrument. Finally, in Sect. 4 we summarize our main results and provide an outlook for future developments and improvements of the nonLTE model and ARTS.
2. Method of solution nonLTE problem
2.1. Overview
Our main concern is to model radiation intensity associated with rotational transitions of optically active molecules that falls roughly into the submillimetre wave range. Due to a small energy separation between the internal rotational states of molecules, the inelastic collisions within a volume are usually sufficient to keep the rotational populations in Boltzmann distribution, i.e. LTE, for even relatively low densities. In the LTE regime, the population distribution depends uniquely on the kinetic temperature. This is one of the advantages of submillimetre measurements, the possibility to target specific rotational lines typically of the ground vibrational state that are often safely assumed to be in LTE. In contrast, the vibrational populations of a molecule as a result of larger energy spacing among levels depart from LTE at somewhat high densities in atmospheric conditions typical of the Earth (see LópezPuertas & Taylor 2001; Feofilov & Kutepov 2012, for discussion of rovibrational nonLTE). Nevertheless, in the conditions of cometary comae, icymoon atmospheres, and at high altitudes in planetary atmospheres the collisions are not usually frequent enough to keep even the rotational level populations according to the Boltzmann distribution and accurate nonLTE modelling needs to be applied.
Generally speaking a nonLTE problem consist of two related parts: (1) the solution of the SEE expressing all the processes populating and depopulating the different levels including radiation and (2) the solution of the radiative transfer equation (RTE) integrated over angle and frequency to produce the mean intensity at each atmospheric level (Kutepov et al. 1998), which in turn depends on level populations.
The most obvious method of solving this interdependent problem is lambda iteration (LI), which repeatedly solves the RTE for the current level populations obtained from the SEE until convergence. However, after many iterations the LI method may often yield an incorrect population distribution, especially in the case of optically thick atmospheres; it is a problem of “false convergence” (Lambert et al. 2015) and for atmospheric science applications see the extensive discussion in Kutepov et al. (1998). In order to compute nonLTE populations for atmospheres with high opacities, an approximate (or accelerated) lambda iteration (ALI) method was devised which was capable of avoiding “trapped” photons in the thick lines, providing accurate and much faster convergence on the populations (Rybicki & Hummer 1991, 1992). An even more effective method (in number of iterations) of computing the nonLTE populations relies on the Gauss–Seidel approach for populations corrections during iteration of the SEE and RTE (Bueno & Sainz 1999; Asensio Ramos & Trujillo Bueno 2006). This method will be discussed in detail in the following section.
In the following section we give a brief overview of the nonLTE solution in spherical symmetry. We have adopted a notation system where each variable that needs to be integrated over the frequency and solid angle have been subindexed with ν and p, respectively. The details of p are described in Sect. 2.2. For variables where lines have to be considered individually, the subindex i j is adopted to indicate a transition from state i to state j. Most of these variables are also dependent on altitude, but we omit this subindex for simplicity until the Gauss–Seidel method description (Sect. 2.4.3), where for the sake clarity the ν subindex is implicit. Finally, we adapt the braket notations to clearly indicate access to individual elements of the matrices built below.
2.2. Radiative transfer calculations
In this work the solution of the transfer equation is performed on a set of rays tangent to the shells (see Fig. 1). The specific intensity I is computed along rays of constant impact parameter, p (Anusha et al. 2009; Harper 1994). Each ray in Fig. 1 is related to the shell radius and has angles, θ, with intersected shells. These angles are used to integrate specific intensity over directions, Ω, at each altitude, and it is also required for projecting the radial expansion onto the rays to account for the Doppler shift. The number of shells with radius smaller than the planet radius, R, is defined as NC, and the number of the shells larger than R is the same as the number of the given discrete spatial grid m = 1, 2, …, NP. Therefore, the p grid has the length of NC + NP, as shown in Fig. 1. The incoming rays may include a background radiation term and are being propagated until the tangent layer (TH), or until an intersection with the surface, where TH ≤ R. The outgoing parts are continued from the tangent height or started from the surface until the outer layer.
Fig. 1. Schematic representation of the placement of ray paths (arrows) of each p for computing the radiation field in spherical geometry. Dotted semicircles represent spherical atmospheric layers at a spatial grid of m − 1, m, m + 1, and NP. The solid semicircle represents a surface of a body where the spatial grid m = 0. Dashed semicircles represent shells inside the surface and dotted straight lines represent the tangent lines to the shells (see text for details). 
The macroscopic radiative transfer equation can formally be expressed as
where I_{ν} is the specific intensity [J m^{−2} sr^{−1} Hz^{−1} s^{−1}] defined as the radiative energy passing through a surface normal to the path per unit time, frequency, and solid angle. The emission coefficient, j_{ν}, has units of [J m^{−2} sr^{−1} Hz^{−1} s^{−1}], while the absorption coefficient, α_{ν}, is in [m^{−1}]. The equation describes the change in I_{ν} as the radiation travels along a path, where the distance along the path is given by s. In the subsequent description we limit ourselves to photons along a ray due to specific spectral (bound–bound) transitions, and neglect any continuum sources/sinks of radiation and scattering at a given frequency, ν [Hz]. The microphysical process of emission and absorption can be then expressed as
and
where h is the Planck constant, n_{i} denotes the number density of molecules in an upper level, A_{ij} is Einstein’s spontaneous emission probability of a transition i → j, and ϕ_{ν} is the frequencydependent line emission profile. The local absorption coefficient depends on Einstein’s absorption probability B_{ji} and Einstein’s stimulated emission probability B_{ij}. The Einstein coefficients obey the wellknown relationship, g_{i}B_{ij} = g_{j}B_{ji}, where g_{i} and g_{j} are the degeneracy at level i and j, respectively, and . The parameters ϕ′_{ν} and ϕ″_{ν} describe the line profiles for absorption and stimulated emission, respectively, and n_{j} is the number density of molecules in the lower level. Hereafter we rely on the typical assumption of complete angular and frequency redistribution such that the absorption and the emission profiles are approximately equal, .
In the case of macroscopic velocity fields present in the atmosphere, the emission and absorption coefficients incur an angular dependence. Because the angles Ω at a spatial point m vary with the set of rays p, the Doppler shift for a transition i to j at each spatial point m and ray path p is
where υ_{p} is the radial wind speed and θ_{p} is the angle between the wind direction at altitude m and line of sight placed at p in Fig. 1 (this is the angle θ as in Fig. 1, with explicit dependence on the ray p).
It is convenient to write the RTE in another more compact form by dividing both sides by α_{ν},
This notation introduces very useful physical quantities that turn out to be important for understanding the propagation of radiation along a ray. First, the socalled source function, S ≡ j_{ν}/α_{ν}, and second, the optical depth, dτ_{ν} ≡ ∑α_{ν}ds, where the sum is performed over all active absorber contributions at the given frequency (nonoverlapping lines). There are a number of different solutions to this equation, each with its own strengths, for example the method of short characteristics (SCs; Olson & Kunasz 1987), the Feautrier method (Rybicki & Hummer 1991), and the discontinuous finite element (DFE) method (Castor et al. 1992). In this paper, we only deal with spherically symmetric medium, thus we implemented the very efficient methods of first and secondorder SCs (Olson & Kunasz 1987; see Sect. 2.4.3 and Appendix C for a detailed description).
2.3. Statistical equilibrium
The SEE expresses a balance among all the detailed microphysical processes influencing the individual population levels. In this paper we assume steadystate conditions such that physical processes influencing the population of levels are much faster than variations in the thermophysical conditions of the atmosphere. In order to calculate the population n_{i} of each of the i = 1, 2, …, NL levels the code accounts for radiative and collisional processes (with one or more collisional partners), such that at each spatial point we can write
where t is time [s], C_{ij} [s^{−1}] is the total collisional rate^{2} for transitions from level i to j, and J̄_{ij} is the integrated mean intensity defined as
The line shape function is calculated including the Doppler shift of frequency Eq. (4) in the case of an expanding atmosphere. The first and second term on the right side of Eq. (6) represent the radiative rates between level i and the upper/lower levels, respectively. The third term on the right side of Eq. (6) represents transition rates due to collisional processes.
The resulting system of equations can be formally written as
where A′ is a matrix of size NL × NL whose elements contain the collisional and radiative rates, and n is a vector containing the population of each level. Under the stated assumptions the set of rate equations is singular and requires that one of the rows of the matrix should be replaced with ones, and setting a single value of the righthand side zerovector to the total number of molecules. The resulting equation is
We note that the diagonal of A (i.e. ⟨iAi⟩) represents the sum over the destruction rate coefficient of level i, and ⟨iAj⟩ represents the formation rate coefficient from level i to j. As mentioned, one row of A contains ones, and f is a zerovector, except having the same row replaced by the total number density.
The J̄_{ij} in Eq. (7) is a nonlocal function of the emitting and absorbing conditions throughout the medium, and these conditions depend on local values of the populations. Therefore, solving Eq. (9) for the atomic and molecular population distributions is a nonlocal and, generally, a nonlinear problem. The nonLTE multilevel problem consists of finding the populations n satisfying Eq. (9) which use the value of n at each spatial grid to calculate the radiative rates.
2.4. NonLTE methods
Our nonLTE model adopts the multilevel Gauss–Seidel flavour of the ALI iteration method, which features a very rapid convergence to a true solution (Trujillo Bueno & Fabiani Bendicho 1995). To demonstrate this method’s performance we compare it to two other approaches, namely LI and pure ALI, which are discussed below as precursors to the Gauss–Seidel description.
2.4.1. Lambda iteration method
The LI method uses a simple backandforth iteration between the SEE and RTE until convergence is reached. The lambda operator, Λ, acting on a source function is a convenient way to express the calculation of specific intensity, even though it is not implied that such matrix has to be explicitly formed. In this formalism the specific intensity at given altitude can be written as
where S is the source function vector and Λ_{ν, p} is a matrix whose elements ⟨mΛ_{ν, p}m′⟩ gives the response of the radiation field at point m due to a unitperturbation in the source function at point m′ along the line of sight p. Then the integrated mean intensity J̄_{ij} at each altitude can be written as
where and S^{old} represent the integrated mean intensity and source function vectors obtained from previous iteration, respectively. The line shape function, ϕ_{ν, p}, is calculated at each set of rays p and altitude grid m.
Once the radiative rates are updated the new level populations are calculated again by solving the linear system,
where n^{new} is the calculated level populations, and A^{old} is the A but using integrated mean intensity obtained from the populations calculated by the previous iteration step, n^{old}. The explicit statistical equilibrium equations (Eq. (6)) for Λiteration can be expressed as
The iteration continues until the conditions for convergence are satisfied. The Λiteration is the most direct and simple approach to the problem; however, it is known to be slow to reach convergence for optically thick lines. This occurs because each cycle of the iteration corresponds to photons moving about one mean free path in the medium, hence many iterations are required to move them any substantial distance. Such trapped photons in the cores of optically thick lines provide small updates in the populations (due to radiative field) from one iteration to the next (Rybicki & Hummer 1991; Lambert et al. 2015).
2.4.2. Multilevel accelerated lambda iteration (MALI) method
The multilevel accelerated Λiteration, MALI, has much faster convergence rates than Λiteration because it effectively removes the trapped photons from consideration through a clever analytical step (e.g. Cannon 1973; Rybicki & Hummer 1991; Auer & Paletou 1994).
In ALI, we write the full Λoperator as , where is the approximate (local) version of the Λoperator. Substituting this expression into the formula for specific intensity (10), we obtain
where S^{old} and S^{new} represent the source function vectors obtained from the previous iteration and current iteration, respectively. From Eqs. (7), (11), and (14), the integrated mean intensity at each spatial point is written as
where is the integrated mean intensity obtained from the population calculated by previous iteration step. The parameter is integral over frequency and direction of ⟨mΛ_{ν, p}m⟩ for a transition between i and j at spatial point at m:
Because S^{new} in Eq. (15) is
we can express SEE at each spatial point including the approximate radiative transfer effect as follows:
The terms and in the radiative rates in Eq. (18) are preconditioning the system, and remain linear in populations. We should also note that when the problem converges, Eqs. (14), (15), and (18) become exact as goes to 0. The resulting linear system of equation is
and the iteration proceeds basically in the same was as in the Λiteration, the difference being the preconditioning of the radiative rates.
2.4.3. Multilevel Gauss–Seidel (MUGA) method
The Gauss–Seidel iterative method based on the idea of solving the Eq. (9) sequentially from the lowermost (uppermost) layer toward the uppermost (lowermost) layer by updating the current source function value and performing full angular integration of the specific intensity by using new level populations. This method was first presented by Trujillo Bueno & Fabiani Bendicho (1995), and Fabiani Bendicho et al. (1997) gives a suitable summary of its application to the multilevel problem in Cartesian coordinates. The first implementation of this method in spherical symmetry for problems in astrophysical applications are given in Asensio Ramos & Trujillo Bueno (2006).
As we set out to solve the system of Eq. (19) the mean intensity at each iteration step, J̄_{m,i j} can be corrected by calculating the radiation field with new level populations at the spatial point of m − 1. To update J̄_{m,i j} at spatial point m we separate the specific intensity into outwarddirected intensity (I_{m, p↑}) and inwarddirected intensity (I_{m, p↓}), expressing them as
and
where
These outward and inwarddirected intensities, except for the bottom and top layers, can be updated when Eqs. (20) and (21) are discretized by the secondorder short characteristics method (Olson & Kunasz 1987). The secondorder method interpolates the source function by a parabola going through S_{m − 1}, S_{m}, S_{m + 1} using the Lagrangian interpolation,
and
where the each term with superscript “new” indicates that it is calculated by using new level populations in iteration steps. The details of how to calculate the interpolation coefficients, λ, in Eqs. (23) and (24) are detailed in Appendix C. As described in Olson & Kunasz (1987), the for spatial point m is
Then by using the corrected and integrated mean intensity, , the statistical equilibrium equation can be written as
The expression at spatial grid m is the integrated mean intensity obtained by Eq. (7). The resulting linearized system of the equation at every altitude point is
where indicates using the corrected mean intensity calculated by radiative transfer taking account of perturbation of the population correction at each spatial grid.
As mentioned in Asensio Ramos & Trujillo Bueno (2006), it is important to clarify that the multilevel Gauss–Seidel method requires a specific order of the loops over transitions, angles (direction of the radiation), and spatial points when propagating the radiation. The most external loop is the one over directions, first going inwards and then outwards. The next is the loop over spatial points, followed by the loop over transitions. The most internal ones are the loops over angles and frequencies. The reason is that we need to evaluate the mean intensity for all the radiative transitions in order to be able to perform the population correction before advancing to the following point of the spatial grid.
3. Validation of the MUGA nonLTE method implementation
In this section we describe the validation and performance of our implementation of the nonLTE solver based on the MUGA method (Asensio Ramos & Trujillo Bueno 2006). This is accomplished by comparing calculated populations to other existing nonLTE models for similar applications. In this validation we use a MALI code implementation for rotational level populations, as used in Rezac et al. (2013; hereafter MPSALI), and served as a benchmark for implementing the EP method used for cometary observations (Marshall et al. 2017). In addition, we also use the widely applied and freely available Monte Carlo nonLTE code RATRAN^{3}, described in the paper Hogerheijde & Van Der Tak (2000).
3.1. Input data and computational details
The most frequent and scientifically important observations in the submillimetre wave range from space telescopes target the lowlying rotational transitions of water. For this purpose we work with the orthoH_{2}O nonLTE model typically used in cometary applications (Bensch & Bergin 2004; Hartogh et al. 2011b; Lee et al. 2015). This model considers seven orthoH_{2}O levels with nine transitions (see Fig. 2). We use the LAMDA^{4} molecular database (Schöier et al. 2005) for level energies, degeneracies, frequencies, and spontaneous emission rates for each transition (see Appendix D).
The collisional coefficient data, and their temperature dependences were obtained from Buffa et al. (2000). In this paper we only consider H_{2}O–H_{2}O collisions (i.e. single collisional partner). In future work, when absolute accuracy for the scientific investigations with a proper nonLTE model for Ganymede’s spectral line measurements may be required, additional collision partners shall be considered.
In this validation we first consider a static water atmosphere of Ganymede as simulated by a kinetic model described in Marconi (2007). We consider optically thick and optically thin conditions to test our implementation, which are represented as subsolar latitude (SSL) 10° and 60°, respectively. The associated temperature and number density profiles are shown in Fig. 3.
The extent of the atmosphere is taken up to 450 km from the surface. Above this altitude the water number density becomes too low to significantly contribute to relevant measurements (except perhaps near the subsolar point), but other physical effects not considered in our model may also start to play a significant role in the excitation. Another assumption is that the molecular velocity distribution is Maxwellian, which may not always be the case in reality if the source of water molecules is sputtering (Marconi 2007; Turc et al. 2014). The vertical profiles of number density and temperature are discretized into 100 layers. The frequency grid of each transition extends to plus and minus five times that of the full width at half maximum of the Doppler broadening at the highest temperature. In the presented calculations we distribute frequency steps in such a way that there are 200 grid points for each transition line. The numerical integration in Eq. (7) relies on the trapezoidal rule. Nonetheless, these topics are not of importance for demonstrating the performance of our code implementation.
Fig. 2. Energy level diagram of the orthoH_{2}O molecule showing all the levels considered in this work, shown along with all the transitions. This is complementary to the information in Table D.2. 
Fig. 3. Vertical profiles of temperature (black lines) and number density (blue lines) of Ganymede atmosphere at subsolar latitude 10° (solid lines) and 60° (dashed lines) obtained from the kinetic model of Ganymede’s atmosphere (Marconi 2007). 
In the calculations, we assume an equilibrium between the first atmospheric layer and the surface temperature as the lower boundary condition, for convenience. Based on our initial sensitivity tests, we ignore additional external sources such as solar radiance or cosmic background radiation in these examples. Therefore, the upper boundary condition has no incoming radiation field. When discussing convergence performance of the different approaches we refer to the percent difference changes in the populations between successive iterations,
The iterations proceed until the relative change of populations on each level and each spatial point reach 0.1% for all the test cases.
Finally, the level populations calculated by the different codes are compared in terms of excitation temperature, T_{ex}, defined for each transition as
The T_{ex} provide an intuitive understanding where in the atmosphere the LTE/nonLTE transitions occurs and whether the excitation is sub or superthermal. Nevertheless, T_{ex} only denotes a relative degree of excitation and does not directly translate into spectral emission intensity, which depends on the total number of excited molecules in the given state.
3.2. Results: populations comparison
The profiles of T_{ex} for the optically thick case, SSL 10°, are plotted in Fig. 4.
Fig. 4. Vertical profiles of excitation temperatures of the nine transitions calculated for the atmospheric model corresponding to SSL = 10°. The labels displayed in the upper right corner denote the frequency for the given transition in GHz. The 556.9 GHz is typically one of the most important transitions for cometary and planetary contexts discussed in this text. The kinetic temperature profile is also shown as dashdotted line. 
In this case the populations start to depart from LTE at around 100 km, while the strong 557 GHz transition is kept thermalized to an altitude above 150 km. In the nonLTE region all the excitation temperatures are noticeably lower than kinetic (subthermal), which is understood as a result of the spontaneous emission rate being greater than the collisional rate, and the absorption of upwelling radiation cannot compensate for the emission process. Although we do not run this model with the solar radiation pumping we have estimated that this effect would be negligible in these examples. A more complex issue is the excitation of the 557 GHz transition from the ground up to ≈180 km, where the upper level of this transition is populated strongly due to spontaneous decay (2_{21} → 1_{10}), which leads to a value of T_{ex} slightly higher than the kinetic temperature. For this transition the thermalization process is purely radiative due to opacity effects (exchange of radiation among atmospheric layers) up to 180 km where it cannot compete with the spontaneous emission.
In the optically thin case, the excitation temperatures depart from LTE from the ground as demonstrated in Fig. 5. The collisional frequency is too low to thermalize the levels, and the opacity is too low for radiative exchange among layers to play any role. It is only the upwelling radiation which contributes to a degree to excitation in the nearground altitudes. The combination of optically weak medium and a low kinetic temperature, specifically the 1162 GHz line has an excitation temperature higher than kinetic above ∼100 km. Nevertheless, even in the weak line limit, i.e. J̄_{i j} = 0, this line’s T_{ex} is around 100 K at 100 km altitude in this case.
Fig. 5. Vertical profiles of excitation temperatures of the nine transitions calculated for the atmospheric model corresponding to SSL = 60°. The labels displayed in the upper right corner denote the frequency for the given transition in GHz. The kinetic temperature profile is also shown as dashdotted line. 
Figures 6 and 7 show the absolute differences between MUGA, MPSALI, and RATRAN calculations in terms of T_{ex} for the two Ganymede atmospheric cases.
Overall, in the optically thick case the codes provide very consistent vertical profiles of the populations. The maximum differences of 0.5 K are reached around a narrow altitude region of 200 km for the 557 GHz transition with respect to the MPSALI code (the left panel in Fig. 6), but agree almost perfectly at other altitudes. For the higher lying transitions the MUGA code obtains about 0.5–1 K higher excitation temperatures above 100 km altitude. The largest deviation is just below 1.5 K for the 1162 GHz transition between the highest energy levels considered in this model. The likely explanation for these discrepancies lies in the radiative transfer calculations of the mean intensity that strongly shape the population distribution in this case. This conclusion is further supported by the results in the optically thin scenarios presented below. The comparisons with RATRAN also shows a very good agreement, although there appears to be a semiconstant offset (cold bias in RATRAN calculations) of roughly 1.6 K as shown in the right panel in Fig. 6. We should note that RATRAN is not designed specifically to handle lower boundary conditions, so we artificially imposed in the first 100 m of the atmosphere an infinite optical thickness. The code then treats this implicitly as a surface. Nevertheless, if we were to adjust the offset the absolute differences in value would not exceed ∼2 K. Another obvious and expected feature of the RATRAN calculations is the random scatter of the T_{ex} on the order of 0.5 K, as we set up this run with a high signaltonoise ratio. The signaltonoise ratio and initial number of photons per cell for the RATRAN calculations were 50 and 10000, respectively. This is typically not a problem (for spectral line simulation) if only forward calculations are required, and shows that RATRAN, with slight modification, can be applied to such cases.
Fig. 6. Left panel: profiles of absolute difference of the excitation temperatures calculated in this work (MUGA) and the MPSALI code and right panel: the same, but compared against calculations made by the RATRAN package. The atmospheric inputs correspond to the SSL = 10° conditions (see Fig. 4 for T_{ex} magnitude). Lines of different colours correspond to the transitions as shown in the label in units of GHz. 
Fig. 7. As in Fig. 6, but for atmospheric conditions of SSL = 60° (optically thin case; see Fig. 5 for T_{ex} magnitude). 
The optically thin example (Fig. 7) demonstrates that MUGA, MPSALI, and RATRAN generally provide the same altitude profiles of populations with differences smaller than 1 K, typically less than 0.5 K. We note that larger deviations between MUGA and MPSALI are below 50 km; this is where the lower boundary emissions have the greatest influence on the populations. Above that, the very optically thin conditions mean that the radiative transfer plays no role, and we do not see any vertical structure in the T_{ex}, as was the case for the optically thick scenario. This indicates that the MUGA and MPSALI have small differences in handling the mean intensity calculations and layertolayer radiative transfer. The RATRAN calculations are performed in the same manner as described above, with the first 100 m of atmosphere acting as an artificial surface to get the lower boundary condition effects. Judging from the right panel in Fig. 7, a small but constant negative bias remains (perhaps 0.5 K), although the random noise has increased to about 1 K. This is a wellknown feature of Monte Carlo radiative transfer codes, it is difficult to build statistics for layers which are more than 99% transparent to the photons. Nevertheless, even in the optically thin conditions the comparisons show a good agreement among the codes.
The MUGA implementation appears to provide accurate populations, and here we discuss its superiority in terms of how many iterations are required to reach subsequent population changes of less than 0.1%. We choose a realistic case and show the maximum change in populations versus iterations for the optically thick scenario of Ganymede atmosphere presented previously. The Fig. 8 contrasts the convergence evolution for MUGA, MALI, and LI methods.
Fig. 8. Convergence plots by different iteration methods, MUGA (solid line), MALI (dashed line), and LI (dotted line), for the sevenlevel H_{2}O model of Ganymede’s atmosphere at SSL = 10^{o }. The vertical axis shows the maximum values of the relative difference between populations at each iteration step and the population which solved the system. 
This figure only demonstrates the wellknown fact that MUGA provides about a factor of 2 faster convergence compared to MALI. In this example, which is not extremely optically thick, the MUGA provides only a factor 3 improvement over LI. However, this is not a constant factor (as in MALI) since for even larger optical depths the LI method could fail to converge altogether due to the trapped photons (see similar plots and discussion in Bueno & Sainz 1999; Asensio Ramos & Trujillo Bueno 2006; Paletou & Léger 2007). The computation time for solving the nonLTE problem depends linearly on the number of spatial grid points in atmosphere (NP), the number of frequency points, the number of rays (NC+NP), and the number of iterations as described in Asensio Ramos & Trujillo Bueno (2006).
3.3. Results: spectral line simulations with ARTS
In this section we demonstrate the working interface between the nonLTE code, implemented in the python package “typhon”, used to calculate the level populations for a set of input data and ARTS, the general code for outgoing spectra synthesis. In the following examples we use ARTS to calculate synthetic pencil beam spectra for nadir and limb geometries for Ganymede’s atmospheres at the 557 GHz water transition, most often used for detection in the submillimetre wavelengths. As was demonstrated in Sect. 3.2, the populations depart from LTE even in the relatively high density region near the subsolar point. Therefore, to show the importance of taking into account these effects we use both LTE and nonLTE line simulations (Fig. 9) in the Rayleigh–Jeans (RJ) brightness temperature. The receiver is placed at 500 km.
Fig. 9. Pencil beam simulated spectra of H_{2}O(1_{10} − 1_{01}) for limb (black lines) and nadir (red lines) geometries for Ganymede’s atmosphere with different SSL, 10° (panel a) and 60° (panel b). The receiver was positioned just outside the atmosphere at 500 km. The tangent heights of the line of sight for the limb geometries are 250 km and 50 km altitudes for SSL 10° and 60°, respectively. Dashed lines represent the LTE simulated spectra and solid lines represent the simulated spectra including the nonLTE model. 
It is evident that in both geometries the nonLTE effects are pronounced. The LTE/nonLTE differences are consistent with the excitation temperatures shown in Figs. 4 and 5. The subthermal excitation yields nonLTE limb emissions weaker than LTE, while absorption lines are deeper than for the LTE case, as more molecules are in the ground state.
Finally, we demonstrate that the presented MUGA nonLTE code solver in combination with ARTS spectra simulator can be already applied to real measurements. For this purpose we take the submillimetre observations of cometary (67P/ChuryumovGerasimenko) spectra obtained by the MIRO instruments on board the Rosetta spacecraft. The instrument and early science results with the measurements are presented in Gulkis et al. (2007) and Gulkis et al. (2015), respectively. The aim here is not to derive any scientific results, conclusions, or parameters of interest from the measurement, but rather to demonstrate the applicability of the code. Therefore, we selected an ad hoc sample calibrated spectrum of the orthowater at 557 GHz from the MIRO database^{5}, corresponding to the beginning of “nearnucleus” observations. There are two reasons for the choice of this period: (1) an adequate parametrized model for the temperature, density, and velocity fields have already been published (Lee et al. 2015) and (2) the coma density for this period is just large enough to provide a moderately optical thickness, but too small to completely saturate the line, which would require a much more detailed coma model.
The MIRO measured spectra taken at 20140808T03:58:32 in nadir geometry is shown in Fig. 10. The spacecraft altitude was 82 km from the comet centre at this time. To increase the signaltonoise ratio we averaged a total of ten minutes of observations, as in Lee et al. (2015), such that the time stamp above indicates the centre time.
Fig. 10. Top panel: comparison of the measurement (black line) and forward model (red line) spectra for the orthoH_{2}O (557 GHz) transition in nadir geometry (see text for details). Bottom panel: residual between the simulation and measurement spectra. 
The simulated MIRO spectra rely on the same sevenlevel nonLTE model as described for our Ganymede simulations, which is the same as used by Lee et al. (2015). The five parameters needed for the coma model (see Table 1 in Lee et al. 2015) were estimated manually by comparing the measured and simulated spectra until a satisfactory fit was found. The physical parameters are production rate, Q, and terminal expansion velocity, ν_{0}. The other two coefficients are used in a power law for temperature and another in hyperbolictangent function for velocity. The original paper provides retrieved values for Q and ν_{0}. Our estimate for the sample spectra yields a production rate (≈Q = 2 × 10^{26} sec^{−1}) and ν_{0} = 0.66 km s^{−1}, which is consistent with the overall results of Lee et al. (2015). We did not perform an inverse problem to find the optimal parameters, which is beyond the scope of this work, and fitting other observations (with large optical depth) may even require a more complex coma model beyond the five parameters (Lee et al. 2015). Nevertheless, we demonstrate for a suitable set of molecular and atmospheric inputs the new implemented MUGA code for rotational transitions coupled with ARTS can be used as an effective tool for scientific investigations of relevant measurements in nonLTE conditions.
4. Conclusions
We presented an implementation of an accurate and efficient numerical method of the solution for the nonLTE problem in spherical symmetry. This nonLTE code is aimed at applications where populations of pure rotational transitions of a given vibrational level are required, for example interpretation of submillimetre observations of cometary comae or icy moon atmospheres.
The adopted nonLTE method relies on the stateoftheart approach developed in the context of stellar astrophysics, the so called Gauss–Seidel ALI method (Trujillo Bueno & Fabiani Bendicho 1995; Asensio Ramos & Trujillo Bueno 2006). This method has two advantages: it is deterministic, so that there is no random scatter in populations typical of Monte Carlo based RT solvers such as RATRAN, and b) it explicitly treats the RT among atmospheric layers, not relying on the simplifying assumption needed, for example, in the escape probability approach. The MUGA method also features a very rapid convergence to the solution, nearly a factor of two improvement relative to the popular ALI method for optically thick conditions. In addition, the ALI and MUGA approaches avoid the problem of extremely slow convergence which ultimately leads to an incorrect solution typical of the pure LI approach in optically thick media.
The nonLTE package is implemented in the Python language, with the goal of making this code a part of the ARTS radiative transfer package. The present communication with ARTS is via its Python interface, the equally opensourced typhon. ARTS uses the relative level populations passed via this interface to calculate emission and absorption coefficients^{6} for spectra calculation in any available geometry accounting for beam convolution if required in the simulation. The presented nonLTE code^{7} is available under MIT license and ARTS^{8} is available via GPL license.
In order to verify and validate our implementation of the MUGA method we compared the calculated nonLTE populations under optically thin and thick conditions in the context of Ganymede’s water atmosphere against other codes routinely applied to the nonLTE conditions of comets and other astrophysical sources, such as the RATRAN and ALI codes (Rezac et al. 2013; Marshall et al. 2017). Ganymede’s atmosphere is assumed to be static (no expansion velocity). Under the optically thick conditions we demonstrate very good agreements between the available codes; expressed as excitation temperatures for each transition, the differences are below 0.5 K for the main H_{2}O line (557 GHz) and below 1.5 K for any of the nine lines applied in this nonLTE model. For the optically thin conditions the nonLTE calculations become sensitive to the lower boundary conditions; however, even here the codes show very good agreement, T_{ex} less than 1 K for all transitions at all altitudes. The RATRAN 1D Monte Carlo approach shows a little larger discrepancy (random noise and small offset); however, we acknowledge that it is not designed explicitly to treat this kind of application, with static atmosphere with important contribution from the surface boundary condition. Nevertheless, it still compares well with the MUGA/ALI results modelling the Ganymede nonLTE atmosphere.
The spectra calculations performed by the ARTS package relying on the interface to the typhon package nonLTE calculations of vertical profile of populations are also presented. The nadir and limb observational geometries were compared, and showed very good agreement in the LTE tests case and also for the full nonLTE calculations. We have also applied the typhon and ARTS codes to simulations of expanding cometary atmospheres under moderately optically thick conditions, and compared them to the nadir measured line profile at 557 GHz obtained by the MIRO/Rosetta instrument. The coma atmosphere was estimated using the parametrized profiles described in Lee et al. (2015).
The new nonLTE package makes a great addition to the ARTS radiative transfer package already established and adapted in many application in atmospheric and planetary science. We demonstrate the accuracy and efficiency of the nonLTE method for the common examples of atmospheres where they are observable. Although we focus only on H_{2}O transitions in this work, any other molecule can be treated since the simulations are fully determined by the atmospheric, spectroscopic, molecular, and collisional data inputs. In the future, the nonLTE package can be further improved, for example to treat line overlapping in the case of hyperfine lines. In the near future, we will also implement the contribution of solar radiation in the excitation of the rotational populations via so called gfactors, as described in the work of Zakharov et al. (2007). The future applications of this code can be both the synthetic simulations of future missions (JUICE) or an existing dataset on cometary comas (MIRO).
see Appendix B for details.
see Appendix A for additional details.
Acknowledgments
We acknowledge the Principal Investigator M. Hofstadter (NASA/JPL, Pasadena, USA) and the ESA Planetary Science Archive for making available the MIRO data used in this paper. The spectroscopic and collisional data for the nonLTE simulations of orthowater are freely available from the LAMDA database (http://home.strw.leidenuniv.nl/~moldata/). The calibrated MIRO spectra (20140808T03:58:32) used in this work are also available from the ESA archive (ftp://psa.esac.esa.int/pub/mirror/INTERNATIONALROSETTAMISSION/MIRO/). The atmospheric profiles have been extracted from the paper of Marconi (2007) and are available from the author upon request.
References
 Anusha, L. S., Nagendra, K. N., Paletou, F., & Léger, L. 2009, ApJ, 704, 661 [NASA ADS] [CrossRef] [Google Scholar]
 Asensio Ramos, A., & Trujillo Bueno, J. 2006, EAS Pub. Ser., 18, 25 [CrossRef] [Google Scholar]
 Auer, L. H., & Paletou, F. 1994, A&A, 284, 675 [NASA ADS] [Google Scholar]
 Bensch, F., & Bergin, E. A. 2004, ApJ, 615, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Bensch, F., Melnick, G. J., Neufeld, D. A., et al. 2007, Icarus, 191, 267 [CrossRef] [Google Scholar]
 Biver, N., BockeléeMorvan, D., Crovisier, J., et al. 2007, Planet. Space Sci., 55, 1058 [NASA ADS] [CrossRef] [Google Scholar]
 Biver, N., Crovisier, J., BockeléeMorvan, D., et al. 2012, A&A, 539, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Biver, N., Hofstadter, M., Gulkis, S., et al. 2015, A&A, 583, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BockeleeMorvan, D. 1987, A&A, 181, 169 [NASA ADS] [Google Scholar]
 BockeléeMorvan, D., Biver, N., Swinyard, B., et al. 2012, A&A, 544, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Buehler, S. A., Mendrok, J., & Eriksson, P. 2017, Geoscientific Model Development Discussions, 1, 436 [Google Scholar]
 Bueno, J. T., & Sainz, R. M. 1999, ApJ, 1, 436 [NASA ADS] [CrossRef] [Google Scholar]
 Buffa, G., Tarrini, O., Scappini, F., & CecchiPestellini, C. 2000, ApJS, 128, 597 [NASA ADS] [CrossRef] [Google Scholar]
 Cannon, C. J. 1973, ApJ, 185, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I., Dykema, P. G., & Klein, R. I. 1992, ApJ, 387, 561 [NASA ADS] [CrossRef] [Google Scholar]
 de Graauw, T., Helmich, F. P., Phillips, T. G., et al. 2010, A&A, 518, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de ValBorro, M., Hartogh, P., Jarchow, C., et al. 2012a, A&A, 545, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de ValBorro, M., Rezac, L., Hartogh, P., et al. 2012b, A&A, 546, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de ValBorro, M., Küppers, M., Hartogh, P., et al. 2013, A&A, 559, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eriksson, P., Buehler, S., Davis, C., Emde, C., & Lemke, O. 2011, J. Quant. Spectr. Rad. Transf., 112, 1551 [NASA ADS] [CrossRef] [Google Scholar]
 Fabiani Bendicho, P., Bueno, J. T., & Auer, L. 1997, A&A, 324, 161 [NASA ADS] [Google Scholar]
 Feofilov, A. G., & Kutepov, A. A. 2012, Surv. Geophys., 33, 1231 [NASA ADS] [CrossRef] [Google Scholar]
 Grasset, O., Dougherty, M. K., Coustenis, A., et al. 2013, Planet. Space Sci., 78, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Gulkis, S., & Allen, M. 2015, Science, 347, aaa0709 [Google Scholar]
 Gulkis, S., Frerking, M., Crovisier, J., et al. 2007, Space Sci. Rev., 128, 561 [NASA ADS] [CrossRef] [Google Scholar]
 Harper, G. M. 1994, MNRAS, 268, 894 [NASA ADS] [CrossRef] [Google Scholar]
 Hartogh, P., Lellouch, E., Crovisier, J., et al. 2009, Planet. Space Sci., 57, 1596 [NASA ADS] [CrossRef] [Google Scholar]
 Hartogh, P., Lellouch, E., Moreno, R., et al. 2011a, A&A, 532, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hartogh, P., Lis, D. C., BockeléeMorvan, D., et al. 2011b, Nature, 478, 218 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Hogerheijde, M. R., & Van Der Tak, F. F. S. 2000, A&A, 362, 697 [NASA ADS] [Google Scholar]
 Kunasz, P. B., & Hummer, D. G. 1974a, MNRAS, 166, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Kunasz, P. B., & Hummer, D. G. 1974b, MNRAS, 166, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Küppers, M., O’Rourke, L., BockeléeMorvan, D., et al. 2014, Nature, 505, 525 [NASA ADS] [CrossRef] [Google Scholar]
 Kutepov, A., Oelhaf, H., & Fischer, H. 1997, J. Quant. Spectr. Rad. Transf., 57, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Kutepov, A., Gusev, O., & Ogibalov, V. 1998, J. Quant. Spectr. Rad. Transf., 60, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Lambert, J., Paletou, F., Josselin, E., & Glorian, J. M. 2015, Eur. J. Phys., 37, 15063 [Google Scholar]
 Lee, S., von Allmen, P., Kamp, L., Gulkis, S., & Davidsson, B. 2011, Icarus, 215, 721 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, S., von Allmen, P., Allen, M., et al. 2015, A&A, 583, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lellouch, E. 2008, Ap&SS, 313, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Lis, D. C., Biver, N., BockeléeMorvan, D., et al. 2013, ApJ, 774, L3 [NASA ADS] [CrossRef] [Google Scholar]
 LópezPuertas, M., & Taylor, F. W. 2001, NonLTE Radiative Transfer in the Atmosphere (World Scientific) [CrossRef] [Google Scholar]
 Marconi, M. 2007, Icarus, 190, 155 [NASA ADS] [CrossRef] [Google Scholar]
 Marshall, D. W., Hartogh, P., Rezac, L., et al. 2017, A&A, 603, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Neufeld, D. A., Stauffer, J. R., Bergin, E. A., et al. 2000, ApJ, 539, L151 [NASA ADS] [CrossRef] [Google Scholar]
 Olson, G. L., & Kunasz, P. 1987, J. Quant. Spectr. Rad. Transf., 38, 325 [NASA ADS] [CrossRef] [Google Scholar]
 O’Rourke, L., Snodgrass, C., de ValBorro, M., et al. 2013, ApJ, 774, L13 [NASA ADS] [CrossRef] [Google Scholar]
 Pack, R. T. 1979, J. Chem. Phys., 70, 3424 [NASA ADS] [CrossRef] [Google Scholar]
 Paletou, F., & Léger, L. 2007, J. Quant. Spectr. Rad. Transf., 103, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Rezac, L., Kutepov, A. A., Faure, A., Hartogh, P., & Feofilov, A. G. 2013, A&A, 555, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rybicki, G. B. 1985, in Methods in Radiative Transfer, ed. W. Kalkofen (Netherlands: Springer), 21 [Google Scholar]
 Rybicki, G., & Hummer, D. 1991, A&A, 245, 171 [NASA ADS] [Google Scholar]
 Rybicki, G. B., & Hummer, D. G. 1992, A&A, 262, 209 [NASA ADS] [Google Scholar]
 Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Trujillo Bueno, J., & Fabiani Bendicho, P. 1995, ApJ, 455, 646 [NASA ADS] [CrossRef] [Google Scholar]
 Turc, L., Leclercq, L., Leblanc, F., Modolo, R., & Chaufray, J.Y. 2014, ICARUS, 229, 157 [Google Scholar]
 Zakharov, V., BockeléeMorvan, D., Biver, N., Crovisier, J., & Lecacheux, A. 2007, A&A, 473, 303 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: A note on the ARTS implementation
The main calculations of state level distributions are performed outside of the main ARTS program, in a module of its accompanying python package: typhon. In ARTS absorption can be treated for a myriad of species and types of absorption from linebyline absorption to collisioninduced absorption, and various types of empirical/analytical continua. In addition, scattering is considered, adding to complications at even the most fundamental level of how to define the radiative transfer equation.
To fit within this wider theoretical framework, while not being detrimental to the speed of the algorithm in general, the main difference between the internal formalism in this work and the practical implementation in the ARTS code is that j_{ν} is computed indirectly as a relative ratio while being summed up. That is, for a single absorption line we have
where
While this might appear to be a great slowdown and cumbersome workaround, it allows us to add scattering absorption to α_{ν} of Eq. (A.1) without changing the treatment of the individual line in Eq. (A.2). It also means that we do not need to touch j_{ν} for overlapping lines or for cases where we can treat some absorption lines with LTE formalism. The latter is not expected for rotational nonLTE, but can happen when vibrational nonLTE matters.
Appendix B: LTE populations where collisional processes dominate the transitions
The collisional rate, C_{ij}, described in Eq. (6) is equal to
where n_{col} is the number density of the collision partner (in cm^{−3}) and γ_{ul} is the downward collisional rate coefficient (in cm^{3} s^{−1}), usually temperature dependent. In the case of multiple collisional partners the individual collisional rates are summed up.
The rate coefficient is the Maxwellian average of the collision cross section, σ,
where E is the collision energy, k is the Boltzmann constant, and μ is the reduced mass of the system. The upward rate is obtained through detailed balance
Therefore, the population ratio when the thermal collisional processes dominate the transitions can be easily derived as
(B.4) Kutepov et al. (1998) nicely explain the population in the TE and LTE state. Because the rate constant of γ is related to the Lorentz halfwidths, Δν_{L}, of the transition (see Pack 1979),
the rate constant of γ_{ul} can be described as
where α_{L} is the Lorentz halfwidth at the reference temperature T_{0} and reference pressure p_{0}, p is air pressure, and c is the speed of light. This approach was used by Kutepov et al. (1997) for their study of rotational nonLTE for CO.
Appendix C: Computation definitions of λ in Eqs. (23) and (24) in first and secondorder SC methods
The λ parameters in the Eq. (23) are calculated by
for secondorder SC, and for Eq. (24)
The intensities for any angle at the bottom and top atmospheric layer is calculated by firstorder short characteristic method. The first order SC formula of Eqs. (C.1) are simplified to
and Eqs. (C.2) are simplified to
Appendix D: Input data of molecular parameters in this study
Table D.1 shows the level energies and degeneracies at each level of orthoH_{2}O rotational state we considered in this study. Table D.2 shows frequencies and spontaneous emission coefficients for each transition in this study.
Rotational energies and degeneracy of orthoH_{2}O rotational states considered.
Spontaneous emission rates and transition frequencies of the nine rotational transitions of orthoH_{2}O molecule considered.
All Tables
Rotational energies and degeneracy of orthoH_{2}O rotational states considered.
Spontaneous emission rates and transition frequencies of the nine rotational transitions of orthoH_{2}O molecule considered.
All Figures
Fig. 1. Schematic representation of the placement of ray paths (arrows) of each p for computing the radiation field in spherical geometry. Dotted semicircles represent spherical atmospheric layers at a spatial grid of m − 1, m, m + 1, and NP. The solid semicircle represents a surface of a body where the spatial grid m = 0. Dashed semicircles represent shells inside the surface and dotted straight lines represent the tangent lines to the shells (see text for details). 

In the text 
Fig. 2. Energy level diagram of the orthoH_{2}O molecule showing all the levels considered in this work, shown along with all the transitions. This is complementary to the information in Table D.2. 

In the text 
Fig. 3. Vertical profiles of temperature (black lines) and number density (blue lines) of Ganymede atmosphere at subsolar latitude 10° (solid lines) and 60° (dashed lines) obtained from the kinetic model of Ganymede’s atmosphere (Marconi 2007). 

In the text 
Fig. 4. Vertical profiles of excitation temperatures of the nine transitions calculated for the atmospheric model corresponding to SSL = 10°. The labels displayed in the upper right corner denote the frequency for the given transition in GHz. The 556.9 GHz is typically one of the most important transitions for cometary and planetary contexts discussed in this text. The kinetic temperature profile is also shown as dashdotted line. 

In the text 
Fig. 5. Vertical profiles of excitation temperatures of the nine transitions calculated for the atmospheric model corresponding to SSL = 60°. The labels displayed in the upper right corner denote the frequency for the given transition in GHz. The kinetic temperature profile is also shown as dashdotted line. 

In the text 
Fig. 6. Left panel: profiles of absolute difference of the excitation temperatures calculated in this work (MUGA) and the MPSALI code and right panel: the same, but compared against calculations made by the RATRAN package. The atmospheric inputs correspond to the SSL = 10° conditions (see Fig. 4 for T_{ex} magnitude). Lines of different colours correspond to the transitions as shown in the label in units of GHz. 

In the text 
Fig. 7. As in Fig. 6, but for atmospheric conditions of SSL = 60° (optically thin case; see Fig. 5 for T_{ex} magnitude). 

In the text 
Fig. 8. Convergence plots by different iteration methods, MUGA (solid line), MALI (dashed line), and LI (dotted line), for the sevenlevel H_{2}O model of Ganymede’s atmosphere at SSL = 10^{o }. The vertical axis shows the maximum values of the relative difference between populations at each iteration step and the population which solved the system. 

In the text 
Fig. 9. Pencil beam simulated spectra of H_{2}O(1_{10} − 1_{01}) for limb (black lines) and nadir (red lines) geometries for Ganymede’s atmosphere with different SSL, 10° (panel a) and 60° (panel b). The receiver was positioned just outside the atmosphere at 500 km. The tangent heights of the line of sight for the limb geometries are 250 km and 50 km altitudes for SSL 10° and 60°, respectively. Dashed lines represent the LTE simulated spectra and solid lines represent the simulated spectra including the nonLTE model. 

In the text 
Fig. 10. Top panel: comparison of the measurement (black line) and forward model (red line) spectra for the orthoH_{2}O (557 GHz) transition in nadir geometry (see text for details). Bottom panel: residual between the simulation and measurement spectra. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.