A&A 488, 1237-1247 (2008)
F. Daniel1 - J. Cernicharo1,2
1 - Departamento de AstrofÃsica Molecular e Infrarroja, Instituo de la Estructure de la Materia, CSIC, Serrano 121, 28006 Madrid, Spain
2 - Laboratoire d'Astrophysique de l'Observatoire de Grenoble, 414 rue de la Piscine, BP 52, 38041 Grenoble Cedex 9, France
Received 30 April 2008 / Accepted 13 June 2008
Context. The improvement in observational facilities requires refining the modelling of the geometrical structures of astrophysical objects. Nevertheless, for complex problems such as line overlap in molecules showing hyperfine structure, a detailed analysis still requires a large amount of computing time and thus, misinterpretation cannot be dismissed due to an undersampling of the whole space of parameters.
Aims. We extend the discussion of the implementation of the Gauss-Seidel algorithm in spherical geometry and include the case of hyperfine line overlap.
Methods. We first review the basics of the short characteristics method that is used to solve the radiative transfer equations. Details are given on the determination of the Lambda operator in spherical geometry. The Gauss-Seidel algorithm is then described and, by analogy to the plan-parallel case, we see how to introduce it in spherical geometry. Doing so requires some approximations in order to keep the algorithm competitive. Finally, line overlap effects are included.
Results. The convergence speed of the algorithm is compared to the usual Jacobi iterative schemes. The gain in the number of iterations is typically factors of 2 and 4 for the two implementations made of the Gauss-Seidel algorithm. This is obtained despite the introduction of approximations in the algorithm. A comparison of results obtained with and without line overlaps for N2H+, HCN, and HNC shows that the J=3-2 line intensities are significantly underestimated in models where line overlap is neglected.
Key words: line: formation - radiative transfer - methods: numerical - ISM: molecules
Our understanding of the evolution of the molecular complexity in astrophysical media relies both on the quality of the observations and on their interpretation. The increasing spatial resolution that will be provided by the new generation of ground- and space-based telescopes (Herschel, ALMA, JWST) will enable finer geometrical structures to be infered in the objects under study. Hence, the best way to interpret the origin of the molecular line emission and to derive physical and chemical conditions is to use non local radiative transfer calculations with realistic geometries adapted to the angular resolution. While the Sobolev approximation can be safely used when a large velocity gradient is present, interstellar clouds often show linewidths that are thermal in nature or that involve low-velocity fields, so the use of such an approximation is rather doubtful, and a widely used method in the field of molecular astrophysics is the Monte-Carlo technique (Bernes 1979). This method is particularly attractive due to its conceptual simplicity and because it is easily implemented to treat multidimensional geometries. Nevertheless, an inherent drawback of the method is that it suffers from numerical noise and presents a slow convergence rate when optically thick lines are modelled. Although line overlaps can be easily implemented (Gonzalez-Alfonso & Cernicharo 1993), computing time becomes very expensive when treating a case such as N2H+ where the two nitrogen atoms contribute to the hyperfine structure.
Since the 90s, various techniques have been introduced to solve the coupled set of radiative transfer and statistical equilibrium equations. These methods lead to substantial gains in term of calculation speed and accuracy of the solutions. Among them, the preconditioning of statistical equilibrium equations entails drastic improvement in the convergence properties of the iterative scheme for problems with large optical depths (see Rybicki & Hummer 1992,1991). A crucial issue is that the Lambda iteration scheme does not always converge to the true solution for large optical depths (see Mihalas 1978; Auer 1991). Another progress consisted in the introduction by Trujillo Bueno & Fabiani Bendicho (1995) of the Gauss-Seidel algorithm for the solution of radiative transfer problems.
The description of this method was first made for the case of one-dimensional plane-parallel stellar atmospheres (Paletou & Léger 2007; Trujillo Bueno & Fabiani Bendicho 1995) and was further generalised to the case of two and three-dimensional radiative transfer problems with multilevel atoms in Cartesian coordinates (Léger et al. 2007; Fabiani Bendicho et al. 1997) and to the case of scattering line polarization (Trujillo Bueno & Manso Sainz 1999). It was then described (Asensio Ramos & Trujillo Bueno 2006) and used (Cernicharo et al. 2006) in spherical geometry. Although the geometry of interstellar clouds cannot always be assumed to be as simple as spherical, it enables a first guess more appropriate than plan-parallel geometry in order to model centrally condensed clouds. In this article, we give more details on the implementation in spherical geometry, and we extend the discussion to the case of line overlap. We find that in spherical geometry, the requirements of the original Gauss-Seidel algorithm cannot be fulfilled without increasing the computational time (even in the case of isolated lines). To address this problem, we introduced several approximations into the algorithm. Interestingly, despite these approximations, the usual increase in reaching convergence in terms of number of iterations persists. Moreover, these approximations do not entail false convergence, in comparison to results obtained within the iteration scheme.
Several molecules used to trace the physical conditions in molecular clouds harbour a hyperfine structure that complicates the analysis of the observed emission (see e.g. Daniel et al. 2007). Indeed, the influence of the overlap is usually not considered when interpreting observations despite the fact that for a molecule with hyperfine structure, the line blending increases with the rotational quantum number. For high rotational transitions, all the hyperfine lines are usually overlapped in astrophysical objects. Previous works have shown that, for HCN, line overlaps play an important role in the pumping of molecular levels and have to be considered in order to interpret the relative intensities of hyperfine transitions (cf. Guilloteau & Baudry 1981; Gonzalez-Alfonso & Cernicharo 1993). One goal of this paper is to extend the discussion to other molecules, i.e. HNC and N2H+, and to derive the general radiative transfer effects that are common to the different species.
This article is organised as follows. In Sect. 2.1, we first overview the short characteristics (SC) method. In Sect. 2.2, the main equations used in the preconditioning of statistical equations are given. In Sect. 2.3, the implementation of the Gauss-Seidel algorithm in spherical geometry is discussed. In Sect. 3, a test case for H2CO gives the main properties of the various algorithms. Finally, in Sect. 4, the influence of introducing line overlap in the calculations is discussed for the N2H+, HCN, and HNC molecules.
|Figure 1: Transfer along a short characteristic.|
|Open with DEXTER|
To solve the transfer equation, we use the SC method introduced by Olson & Kunasz (1987). Here, we briefly overview the method and give the main equations in order to introduce quantities that will be used in the following discussion.
If one considers the radiation propagation along
a ray that is discretized (cf. Fig. 1),
the specific intensity between two consecutive points varies as (Olson & Kunasz 1987)
|Figure 2: Geometric quantities in spherical geometry.|
|Open with DEXTER|
In spherical geometry, the radiation field is fully characterised introducing two variables, r and . Numerically, one has to adopt a discretized space for these two variables. At a given radius ri and for an angle , the specific intensity is known within the SC scheme from the geometrical quantities represented in Fig. 2. In the case , we find
is defined according to
The quantities d1,i and d2,i are used to derive and . The angle is used to interpolate the incoming intensity. Moreover, the angle is used in the determination of when a radial velocity field is introduced. For the case , the above formulae also apply, using the angle and then inverting respectively with and d1,i with d2,i.
Following Olson et al. (1986), rather than the lambda iteration scheme (LI), we use an alternative scheme known as the accelerated (or approximated) lambda iteration (ALI) since this one presents better convergence properties. In that scheme, the specific intensity obeys
Several studies aiming at determining the most effective ALO to be used so that the convergence is reached in the fastest way have been conducted. Non diagonal Lambda operators were found to be efficient in reducing the number of iterations but at a computational cost per iteration higher than for a diagonal Lambda operator. Finally, the use of a diagonal ALO seems to lead to the fastest convergence rates (see e.g. Olson et al. 1986; Puls & Herrero 1988).
When we use the diagonal part of the full lambda operator for
we find that
instead of solving the usual statistical equilibrium equations (SEE), one has to solve the system (Rybicki & Hummer 1991)
The overlap between hyperfine lines is treated as
described by Rybicki & Hummer (1992). Here we give a brief overview of
the main equations.
In the case of overlapping lines, the source function without background continuum radiation
is given by
In the appendix, the way the Lambda operator is constructed is presented and the differences with the plan-parallel case are emphasised.
The Gauss-Seidel (GS) iterative scheme was introduced in radiative transfer by Trujillo Bueno & Fabiani Bendicho (1995). The superiority of this algorithm was then shown, by comparison to Jacobi-based ones, in order to rapidly reach convergence even for the complex case of multilevel atoms in two and three-dimensional Cartesian geometry (see Fabiani Bendicho et al. 1997). Recently, a detailed description of the implementation of the GS algorithm for multilevel atomic (or molecular) case was presented by Paletou & Léger (2007). So far, the GS algorithm has been extensively described in plan-parallel geometry. The description of the implementation in spherical geometry has been done Asensio Ramos & Trujillo Bueno (2006), and here we extend the discussion made in this article. Before, we will review the main attributes required by the algorithm.
At each iteration, the grid is first swept from the outermost layer until
the innermost one (say, from
), and during this process, part of the averaged radiation field is computed, i.e. for
Then, once the innermost layer is reached, the grid is swept in the other direction. This leads to the computation of the averaged radiation field that corresponds to
Within the Jacobi iterative scheme, all the grid is swept until reaching the outermost layer, and then statistical equilibrium equations are solved. Within the GS iterative scheme, once the averaged radiation field is known for a layer, the statistical equilibrium equations are
solved for it. At this stage, various updates have to be made to account for the newly derived populations.
Using parabolic interpolation, this implies correcting both the outward and inward radiation field respectively, for the
current and the next grid point. These corrections are made so that when deriving
the populations for a layer k, the populations considered are those obtained during the current iteration
for the layers i<k and those corresponding to the previous iteration for the layers ,
We now describe how the GS algorithm should be implemented in spherical geometry and how we have implemented it in practice. In what follows, we note , the intensity for and the intensity for .
Let us consider first that the grid has been swept from the outermost layer to the innermost one. This leads to the computation of . Once the innermost layer is reached, we can compute the intensities since the incoming intensities are given by (if there is a central source of continuum radiation, the intensities are an input parameter). We then know , and we can invert the statistical equilibrium equations for this layer. Before going to the point k=2, there are three corrections that in principle have to be done to conform with the requirements of the original GS algorithm:
|Figure 3: Dependence of the intensities.|
|Open with DEXTER|
Suppose we reach the point k where new populations are obtained. Let us now consider the corrections that are needed to have an exact GS iterative scheme. Since new populations have been obtained for the layer k, a requirement of the algorithm is that, at the point k+1, the populations should be determined with for and for . Now, consider the way depends on the various populations and so, through the dependance on the incoming intensities. For this purpose, we refer to Fig. 3, which shows the dependencies of intensities. We see that to obtain the intensities , it is necessary to know the inward intensities , as well as the outward intensities , which in turn depend on , etc. Finally, we see that fully accounting for the new population in the computation of would require first to perform the propagation along the inward direction, leading to the updated set . Then, the propagation would be done again along the outward direction so that the set would effectively be given considering the new populations . In practical applications, such a process would lead to computational times longer than those obtained within the Jacobi iterative scheme. Moreover, it should be noted that this problem does not arise in plan-parallel geometry: in that case, once the population is updated at point k, the computation can be performed at point k-1 and will take exactly the new populations into account (cf. Fig. 3), since the radiation field at the innermost boundary is an input parameter of the model.
To make the algorithm computationally competitive, we proceed as follows. By analogy with the plan-parallel case, a way to implement the Gauss-Seidel algorithm would be not to account for the influence of the new populations on the incoming intensities and when we compute the outgoing intensities . This would lead to performing
As explained above, within the usual GS scheme and for each iteration, the grid is first swept from k=N to k=1 and then from k=1 to k=N, and during this second step, statistical equilibrium equations are solved. As mentioned by Trujillo Bueno & Fabiani Bendicho (1995) and Trujillo Bueno & Manso Sainz (1999), a way to increase the convergence speed is to invert the statistical equilibrium equations during both the upward and downward sweeping of the grid.
During the downward propagation and when the iterative process has been initiated, we already know the intensities . Thus, when reaching point k, we directly compute the averaged radiation field for . First we compute and then using and as incoming intensities. Here we note that the intensities in fact involve a mixing of the populations obtained during the two last iterations. This is because the first correction described in the previous section is not performed and that the intensities are saved after the second correction. Finally, since the intensities are computed at this stage, there is only one correction needed so that, after inversion, we perform an update of using the populations , , and . During the downward propagation, we have to save the intensities calculated with the new populations since these intensities will be used during the upward propagation. The upward sweeping of the grid is similar to the one described in the previous section.
The numerical implementation was checked against the results obtained with the code described by Gonzalez-Alfonso & Cernicharo (1993). Various models were used in the comparison, covering different physical conditions, as well as low, moderate, and large optical depths. The differences in the derived level populations ranged from a tenth of a percent to a few percent when integrated optical depths were higher than 100. This corresponds to the expected accuracy of numerical radiative transfer codes as explained in van Zadelhoff et al. (2002).
To accelerate the convergence we used the ALI method described in Sect. 2.2. In spherical geometry, the diagonal of the full Lambda operator is given by Eqs. (B.10)-(B.12). However, determining the exact diagonal is time consuming: first, it requires the evaluation of the opacities that depend on all the internal shells crossed by the rays and second, it requires interpolating in angle the coefficients cm+l+l to obtain them at the values . In practice, we thus chose not to keep the corresponding terms and we approximated the diagonal of the Lambda Operator just keeping, for ,
These terms permit us to accelerate the rate of convergence notably by comparison to the LI scheme. Moreover, the approximation made in deriving the diagonal of the ALO does not induce false convergence. Both aspects can be seen in Fig. 4 where the maximum difference with respect to the converged solution is plotted. The model considered here corresponds to a sphere of uniform physical parameters with n(H 2) = 104 cm-3 and temperature T = 25 K. The sampling consists in 500 spatial grid points and 20 angles. This has been found to insure a convergence for the level populations better than 0.1% for the first 12 levels. The molecule considered is ortho-H2CO, and 14 energy levels were included in the calculation. The molecular parameters were taken from the references given in Table 1. Since no analytical solution exists for this problem, we used as reference populations those obtained with the LI scheme. The error relative to the converged solution for the considered grid, noted , is defined according to Eq. (20) of Auer et al. (1994).
Table 1: References for the molecular parameters used in this work.
Figure 4 also shows the rate of convergence of the different algorithms in terms of the iteration number. For the GS algorithms, two curves are shown: the grey one corresponds to the case where the incoming intensities are corrected (i.e. correction 1 of Sect. 2.3.1) and the black one corresponds to the case where this correction is omitted. We see that, for this model, the ALI scheme speeds up the convergence by a factor 1.3 in comparison to the LI scheme. The GS down and GS up-down enabled us to reduce the number of iterations by factors of 2.0 and 4.4 respectively by comparison to the ALI scheme. In Sect. 2.3, the various corrections to be done for these two algorithms were discussed. For each correction, the code calls a formal solver that takes the source functions of the lines as input and outputs the averaged radiation field and Lambda operator, as well as the intensities. Thus, considering the number of calls to the formal solver and in the different schemes, we have for each layer 2 calls in the ALI scheme, 4 calls in the GS down scheme, and 6 calls in the GS up-down scheme. The number of calls gives the relative time per iteration of each algorithm. In practice, the GS schemes implementation can be improved by performing some bookkeeping, in particular for the incoming intensities. These lead to times per iteration, in the GS down and GS up-down schemes, which are respectively 1.45 and 2.3 times greater than in the ALI scheme. These times were corrected for the time spent in the code in internal storage on the hard drive. Such a procedure is done in our implementation of the GS schemes while not for the ALI one and makes the GS scheme 20% slower than permitted by the algorithmic. Moreover, these longer times are compensated for by the reduced number of iterations required (see above). Thus, for the current implementation, the GS down and GS up-down respectively reduce the total calculation time by factors of 1.4 and 1.9. The numerical implementation of the GS scheme could still be improved to gain some computing time. The corrections in spherical geometry are formally similar to the one in plan-parallel geometry, and for this case, Paletou & Léger (2007) find that the time per iteration of the GS down scheme was only 1.3 longer than that of the ALI scheme. Our code has been written in a very modular and structured way. By decreasing the number of internal calls to routines and functions, we could still gain computing time and approach the 1.3 value of plan-parallel geometry.
Nevertheless, the GS schemes become more advantageous if, rather than parabolic interpolation, we use linear interpolation. In that case, there are fewer corrections to carry out so that the time per iteration of the GS down and GS up-down scheme are only 1.4 and 1.7 that of the ALI scheme. The drawback is that if using linear interpolation, one needs to use more grid points in order to converge to the true solution. Nevertheless, for the model considered here, we find that the differences in level populations are inferior to 1%. In Fig. 5, the line parameters of the most opaque lines are plotted for the two types of interpolations and the results are indistinguishable on the scale of the graph.
|Figure 4: Convergence properties of the various algorithms in terms of number of iterations.|
|Open with DEXTER|
We note that the various algorithms discussed here can still be improved if combined with the Ng acceleration technique (Ng 1974). For the various models tested, we found that it typically leads to a reduction by a factor of 3 in the number of iterations. Alternatively, the successive over-relaxation technique can be used Trujillo Bueno & Fabiani Bendicho (1995).
The numerical implementation of the overlap was checked according to the results reported by Gonzalez-Alfonso & Cernicharo (1993). For the molecules considered in what follows the molecular parameters are summarised in Table 1. Figure 6 shows the variation in excitation temperature versus line opacity for the three HCN hyperfine lines associated with the J=1-0 rotational transition. The physical parameters of the models were taken from Gonzalez-Alfonso & Cernicharo (1993). The line parameters are given with (solid lines) and without (dashed lines) introduction of the overlap. By comparison to the results reported in Figs. 1a-d of Gonzalez-Alfonso & Cernicharo (1993), we see that the qualitative behaviour of is correctly reproduced for all the hyperfine lines when the overlap is introduced. Nevertheless, a quantitative comparison of the line parameters, with and without overlap, shows differences in the order of 5-10%. These differences are due to the choice of grid parameters, especially the number of radial grid points, adopted to reach convergence.
|Figure 5: Excitation temperature versus opacity for a model considering ortho-H2CO.|
|Open with DEXTER|
|Figure 6: Line parameters obtained for the models described in Gonzalez-Alfonso & Cernicharo (1993).|
|Open with DEXTER|
Figure 7 shows the influence of the overlap on the modelization of the hyperfine lines associated with the 3 lowest N2H+ rotational transitions. This molecule is particularly interesting for infering the physical conditions of the densest parts of cold dark clouds. Indeed, it was found that this molecule does not suffer from strong depletion for H2 volume densities below 106 cm-3. The model parameters are identical to those reported in Daniel et al. (2007) for the source L63. In this figure, we see that the J=1-0 line remains mainly unchanged with differences in the line fluxes lower than 5%. The differences increase when considering higher rotational lines. For the J=3-2 line, including the overlap increases the line flux by 20-30%. This can be understood as follows. For a group of blended hyperfine lines, the average radiation field seen by each line is greater when considering the overlap, since more molecules can interact with the radiation. Consequently, the pumping of the upper level of each transition is more efficient, which results in an enhancement of the populations of the highest rotational levels. Therefore the overlap changes the physical parameters reported in Daniel et al. (2007). Mainly, H2 volume densities and/or temperatures were overestimated in this study.
|Figure 7: Effect of the overlap on the line profile of L63. The black line correspond to a model without overlap and the red lines correspond to models that include overlap.|
|Open with DEXTER|
|Figure 8: Effect of overlap on the N2H+, HCN, and HNC J=1-0, 2-1, and 3-2 lines. Red lines correspond to models without introduction of the overlap and black lines to models considering line overlap. The differences between integrated intensities are given for each line. The intensity scale is antenna temperature and are obtained using the standard values of HPBW and beam efficiency of the 30 m IRAM telescope.|
|Open with DEXTER|
The effect induced by line overlap obviously depends on individual hyperfine line opacities, as well as on the frequency offset between the lines. Figure 8 shows the emerging intensities for the three lowest rotational lines of HCN, HNC, and N2H+. The models consist of a uniform density core of diameter 1'. The density is n(H 2) = 105 cm-3, the temperature T=10 K, and the turbulence is set to = 0.15 km s-1. For each molecule, calculations are performed at abundances of 10-10, , and . First, it can be seen that, at the lowest abundance, the effect introduced by the overlap is weak for N2H+ and HCN. The effect is stronger for HNC due to closer hyperfine frequencies. Second, we find that the effect introduced by the overlap is approximatively proportional to the column density, as long as individual line opacities 1. In the present calculations, this is found for the J=3-2 lines of each species for the 2 lowest abundances. For these lines, we see that the integrated intensity ratios of the calculations with and without overlap are approximatively proportional to the abundance ratio. In a more optically thick regime, the increase in the integrated intensity ratio depends less on the column density and finally reaches a saturation regime where the intensity ratio starts to be affected by self-absorption effects. This can be seen in Fig. 8 for the J=2-1 and J=1-0 lines. Finally, the error introduced in the interpretation of observations when neglecting the line overlap can be especially important for HNC. Indeed, we see in Fig. 8 that, for conditions typical of dark clouds, including line overlap leads to an intensity increase of up to a factor 2. Moreover, for similar abundances, the overlap differentially affects HCN and HNC and thus, neglecting it may lead to overestimating the X(HNC)/X(HCN) ratio.
Another effect introduced by the overlap is to modify the relative intensities of the hyperfine components of a given rotational transition. Figure 9 shows the J=2-1 line of N2H+ for a model with abundance . The line obtained without overlap is represented multiplied by a factor of 1.54, which corresponds to the integrated intensity ratio obtained previously. In this figure the position in frequency of the different hyperfine components is also represented as are their line strengths. These are given in percent by reference to the value obtained by summation over all the hyperfine components. Considering the red satellite, we see that, when including the overlap, the intensity associated with the strongest transition (group 1 at 2.5 km s-1) is reduced by comparison to the group of weaker lines situated at 3 km s-1 (group 2). This effect on the profile for the N2H+ J=2-1 red satellite agrees with the observations reported in Daniel et al. (2007). Two effects yield to the current emerging profile. The main effect is introduced by the number of strong transitions that overlap in each group. In group 1, the line strengths of the two weakest transitions are low enough so that there is mostly one line that contributes to the emission. On the other hand, group 2 corresponds to 4 lines with equivalent line strengths. Thus, for this group of lines, the pumping of the upper levels is enhanced producing an increase in the emergent intensities. Another effect of the overlap is to create a flow of population that depends on the relative line strengths of the overlapping lines. Figure 10 shows excitation temperature ratios between the transitions of group 1 obtained with (solid lines) and without (dashed lines) introduction of the overlap. We see in this figure that the overlap enhances the excitation temperatures of the 2 weakest lines. This means that the levels JF1F = 223 and JF1F = 221 are more efficiently populated when considering the overlap due to the presence of the JF1F=212-101 transition that is close in frequency. This can be understood when considering the way the average radiation field is obtained when solving the radiative transfer. For the strongest line, the presence of weak hyperfine components close in frequency does not have a strong influence since the opacities/intensities in the wings of the intrinsic line profiles are low. The average radiation field is thus only slightly modified. On the other hand, for a weak line that is close to the frequency of a strong line, the average radiation field is increased. This will pump the corresponding energy levels and result in an increase in intensity for the weakest lines.
|Figure 9: Effect of overlap on relative line intensities. The black curve corresponds to a model that includes line overlap and the red line to a model without overlap, the latter being scaled by the integrated intensity ratio.|
|Open with DEXTER|
|Figure 10: Excitation temperatures ratios for overlapping lines with different line strengths.|
|Open with DEXTER|
We have presented technical details concerning the determination of the Lambda operator. Moreover, we gave details on the implementation of the Gauss-Seidel algorithm in spherical geometry. Finally, the numerical code was applied to the treatment of line overlap. The main conclusions are:
Part of this work was supported by the AYA2006-14876, ESP2004-00665, ESP2007-65812-C02-01, ESO No. 010103050002, Molecular Universe MRTN-CT-2004-51230, ASTROCAM S-0505 ESP-0237 projects. The authors thank F. Paletou and L. Léger for the help provided concerning the Gauss-Seidel algorithm. We also thank J. R. Pardo for his useful revision of the present manuscript.
As pointed out by Rybicki & Hummer (1991) the Lambda operator is in fact an approximation of the true Lambda operator since it depends on the discretization of the structure of the object. Thus, the Lambda operator is a square matrix whose dimension corresponds to the number of spatial grid points.
Within the SC scheme, the Lambda operator can be constructed easily. As a first step,
from Eqs. (1) and (2) one finds that along a characteristic
|Figure A.1: Notation adopted in plan-parallel geometry.|
|Open with DEXTER|
In plan-parallel geometry, the values for
are readily assigned to values of
the source functions for each layer.
From Fig. (A.1), we have in the case
point k = 0 in expression (A.2) corresponds to the layer n = 1. The point k = i-1
corresponds to the layer n = m. In the above expression, we thus first introduce the variable
n = k+1 and then a notation where all the quantities are indexed according
to the layer, i.e.
Equivalent transformations can be done in the case
Defining N so that n = N-k,
identifying i = N-mand introducing the transformation
In spherical geometry, the same procedure of the plan-parallel case can be followed to obtain the components of the Lambda operator. In this case, however, all the shells are not necessarily crossed for a given direction, and furthermore, some shells are crossed two times before reaching the current point.
|Figure B.1: Definition of geometrical quantities.|
|Open with DEXTER|
We consider a cloud discretized in N concentric shells. In the case , the transformations are similar to the second case examined in the previous section with two exceptions. First, in the SC frame, the coefficients bn+l+l depend on an angle defined in a coordinate system associated to the layer n+l. This dependence was not outlined in the plan-parallel case since all the coefficients along a characteristic are at a unique angle value, namely . In what follows, we introduce the notation , which refers to the angle to be taken into account in the evaluation of b+ln+l (see Fig. B.1). The superscript in this notation distinguishes acute (-) and obtuse (+) angles.
Secondly, depending on the value of ,
the source function corresponding to the point i+1 can be associated either to shell m or
to shell m-1. In this case, we obtain for
|Figure B.2: Definition of geometrical quantities.|
|Open with DEXTER|
Now, let us consider the case , and designate K the innermost shell that is intersected by the ray. This shell is such that we have , with the angle defined according to Eq. (5).
Manipulations of Eq. (A.2) lead to