Issue 
A&A
Volume 625, May 2019



Article Number  A125  
Number of page(s)  15  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201834997  
Published online  24 May 2019 
Nonconservative Rayleigh scattering
A perturbation approach
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, Bd de l’Observatoire, CS34229, 06304 Nice Cedex 4, France
email: frisch@oca.eu
Received:
30
December
2018
Accepted:
18
March
2019
Context. The continuous spectrum of stellar and planetary atmospheres can be linearly polarized by Rayleigh or Thomson scattering. The polarization rate depends on the ratio κ_{c}/(κ_{c} + σ_{c}), κ_{c} and σ_{c} being the absorption coefficients due to photoionizations and scattering processes, respectively. The scattering process is conservative if κ_{c} = 0, and in this case the centertolimb variation of the polarization rate follows Chandrasekhar’s law. Deviations from this law appear if the scattering is nonconservative, that is, if photons have a probability ϵ = κ_{c}/(κ_{c} + σ_{c}) of being destroyed at each scattering.
Aims. Nonconservative Rayleigh scattering is addressed here with a perturbation point of view, using ϵ, assumed to be a constant, as an expansion parameter. The goal is to obtain a perturbation expansion of the polarized radiation field that can be used to measure of the effects of a nonzero ϵ on the polarization rate of the emergent radiation and to check the accuracy of numerical codes.
Methods. The expansion method is an application to Rayleigh scattering of a general perturbation approach developed for scalar monochromatic transport equations. The introduction of a space variable, rescaled by a factor √ϵ, transforms the radiative transfer equation into a new equation from which one can extract simpler equations to describe the field in the interior of the medium and in boundary layers.
Results. The perturbation method is applied to a planeparallel slab with no incident radiation and an unpolarized primary source of photons. The interior and boundary layer fields are expanded in powers of √ϵ. The expansion of the interior radiation field shows that it is unpolarized at leading order, with an intensity i_{0}(τ̃) satisfying a diffusion equation, and that the polarization appears at order ϵ. The emergent radiation is calculated up to and including order ϵ. The leading term yields the polarization rate of the Chandrasekhar’s law. The following one, of order √ϵ, accurately predicts the decrease of the polarization rate for values of ϵ up to 10^{−3} and shows that it varies roughly as (1 − μ) for any unpolarized primary source. Methods for testing the accuracy of numerical schemes are proposed. The perturbation method is also applied to a slab with an incident radiation field and a polarized primary source of photons.
Key words: radiative transfer / polarization / scattering / methods: analytical / stars: atmospheres
© H. Frisch 2019
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
In stellar or planetary atmospheres, the continuous spectra are formed by bound–free and free–free transitions and by scattering on particles with a size smaller than the wavelength, mainly Rayleigh scattering on atoms and molecules in latetype stars, and Thomson scattering on electrons in early type stars. These two scattering mechanisms, which are actually described by the same phase matrix, henceforth simply referred to as the Rayleigh phase matrix, produce a certain amount of linear polarization. The ratio κ_{c}/(κ_{c} + σ_{c}), with σ_{c} being the scattering coefficient and κ_{c} the photoionization coefficient, defines a destruction probability per scattering, here denoted ϵ. The difference 1 − ϵ is the single scattering albedo. The sum κ_{c} + σ_{c} is usually referred to as the absorption coefficient, and κ_{c} as a true absorption or destruction coefficient (e.g., Fluri & Stenflo 1999; Trujillo Bueno & Shchukina 2009; Kostogryz & Berdyugina 2015).
If ϵ = 0, that is, if there is no destruction of photons, the scattering is said to be conservative. Chandrasekhar (1946) showed that conservative Rayleigh scattering in a semiinfinite medium can be solved exactly. He showed in particular that the emergent radiation field for the Milne problem, which describes conservative scattering in a semiinfinite medium with a given total flux F, no incident radiation, and diffuse reflection problems, which describe conservative scattering in a semiinfinite medium with an incident radiation and no primary source, can be expressed in terms of two functions H_{l}(μ) and H_{r}(μ); here μ is the direction cosine of the radiation direction with respect to the normal to the medium. These functions have an explicit expression, albeit slightly complicated. They are tabulated in Chandrasekhar (1960, p. 248) and shown in Fig. 1. The assumption that ϵ = 0 is usually not verified in a stellar or planetary atmosphere, but this exact result, which predicts in particular that the polarization created by Rayleigh scattering cannot exceed 11.7%, is very useful for checking the accuracy of numerical codes (e.g., Fluri & Stenflo 1999).
Fig. 1. Functions H_{l}(μ) and H_{r}(μ). Solid line: H_{l}(μ); dashed line: H_{r}(μ). The functions H_{l}(μ) and H_{r}(μ) are defined for all positive values of μ. For μ → ∞, H_{l}(μ) goes to infinity as and H_{r}(μ) to . The curves have been plotted with the data given in Chandrasekhar (1960). 
If the scattering is nonconservative, that is, ϵ is nonzero, there is no exact expression for the emergent radiation field. For a uniform ϵ, there are two exact results of moderate interest. The law, which holds in the scalar case for a uniform ϵ in a semiinfinite medium with a uniform primary source has a Rayleigh scattering version (Ivanov 1990). This generalized law concerns the surface value of only one quantity, which moreover combines intensity and polarization. This inconvenience can be overcome if some polarization is artificially introduced into the primary source of photons (Frisch 1999; Štěpàn & Bommier 2007). Siewert et al. (1981) were able to construct an exact solution for a matrix H(μ), which, loosely speaking, replaces the functions H_{l}(μ) and H_{r}(μ) and provides the emergent radiation field. The construction of the matrix H(μ) by Siewert et al. (1981) is a real tour de force. However the expression that has been obtained is too unwieldy for practical use. The matrix H(μ) is usually calculated numerically via the solution of a nonlinear integral equation (de Rooij et al. 1989; Ivanov et al. 1996).
Here we propose an expression of the polarized emergent radiation field for nonconservative Rayleigh scattering. It is an approximate one, but its accuracy is fully under control and it can therefore be used for testing the accuracy of numerical schemes. The downside is that ϵ has to be uniform. The construction of this expression is based on a perturbation expansion of the radiative transfer equation in powers of ϵ. For a given ϵ, the accuracy is determined by the number of terms of the series being retained. For ϵ = 0, one recovers the exact Chandrasekhar results for conservative scattering.
The perturbation method is a generalization to Rayleigh scattering of a method developed for unpolarized monochromatic scattering and monokinetic neutron scattering (Larsen & Keller 1974; Frisch & Frisch 1977; Bardos et al. 1983). In a dilute and extended medium, where photons can suffer a large number of scatterings before being destroyed by a bound–free or freetransition, the destruction probability ϵ defines a scale of variation of the radiation field, usually referred to as the thermalization length. For monochromatic scattering, this length, l_{eff}, is of order , for unpolarized (scalar) as well as linearly polarized radiation. In the limit of small ϵ, the radiation field in the interior of the medium can be described in terms of a rescaled space variable . Close to a boundary the variations of the radiation field, in the direction perpendicular to the boundary itself, have to be described on optical depths of order unity. This is achieved by stretching the rescaled space variable in the direction perpendicular to the boundary by the factor l_{eff}. The stretching transforms the boundary layer into a plane parallel semiinfinite medium in which the field is written as the sum of the slowly varying interior field and a rapidly varying boundary layer field. The sum of the two fields must satisfy the boundary condition on the medium, such as a zero incident radiation, and the boundary layer field must go to zero in the interior of the medium.
The introduction of a slow and rapid variable, to adopt the terminology of timedependent problems, allows one to decompose the radiative transfer problem into two equations, one describing the variations of the radiation field in the interior of the medium and the other one acting in boundary layers. For the presentation of the method we use the example of a plane parallel medium. The small parameter of the problem being , we introduce the small parameter . The slow variable is , with τ the optical depth, and the rapid variable. Its domain of variation is [0, ∞], whereas the domain of τ can be finite.
Assuming for each field a perturbation expansion in powers of η, collecting terms with the same powers of η, one obtains two sets of coupled equations, simpler than the original transfer equation, which can be solved explicitly. One thus obtains a perturbation expansion of the emergent radiation field and of the polarization rate in powers of η. The method is presented in detail with an unpolarized internal source and no incident radiation. The case of a polarized primary source and the presence of an incident radiation is considered at the end of the article. We also indicate how this perturbation method could be generalized to multidimensional Rayleigh scattering.
The outline of the article is as follows. In Sects. 2 and 3 we present the radiative transfer equation and the main ingredients of the perturbation analysis, respectively. Section 4 is devoted to the interior field. It is shown in particular that the leading term of the interior field expansion satisfies a diffusion equation. In Sect. 5, which is devoted to the boundary layer, we show that it can be treated as a diffuse reflection problem with conservative scattering. In Sect. 6 we establish boundary conditions for the interior diffusion equation and solve the latter explicitly for a constant primary source of photons. The emergent radiation field is calculated to order η^{2} = ϵ in Sect. 7 and compared to the results of a numerical solution of the radiative transfer equation. We then apply the perturbation method to other physical situations: in Sect. 8 to the case of a polarized primary source and in Sect. 9 to the case of a nonzero incident radiation. The Appendix A is devoted to some technical details.
2. The radiative transfer equation
We first recall the astrophysical context (Sect. 2.1). We then introduce a representation of the radiation field in terms of the irreducible spherical tensors for polarimetry, well adapted to the perturbation analysis (Sect. 2.2).
2.1. The astrophysical context
In a planeparallel atmosphere, a cylindrically symmetric linearly polarized radiation field can be defined by its Stokes parameters I and Q. A linearly polarized continuous radiation field can thus be represented by a frequency dependent twocomponent vector I_{λ}(z, μ), where z is the space variable in the direction of the normal to the atmosphere and μ is the cosine of the direction of the ray with respect to the outward normal. It satisfies the radiative transfer equation,
The pure absorption coefficient κ_{c}(λ, z) comes from bound–free and free–free transitions and the scattering coefficient σ_{c}(λ, z) from Rayleigh and/or Thomson scattering. They have the same phase matrix, here denoted P(μ, μ′) and referred to as the Rayleigh phase matrix. We recall that σ_{c}(λ, z) is independent of frequency for Thomson scattering and varies as 1/λ^{4} for Rayleigh scattering. The primary source term given by the Planck function is unpolarized.
Introducing the monochromatic optical depth,
and the parameter
we can rewrite the radiative transfer equation as
In stellar atmospheres, ϵ_{λ}(τ_{λ}) generally increases towards the interior together with the density and the temperature. Here for the purpose of the perturbation analysis, we assume a uniform ϵ. Since there is no coupling between frequencies, the wavelength λ can be omitted.
Henceforth, we consider the radiative transfer equation
with ϵ being a constant and Q^{*}(τ) an unpolarized primary source. We assume that the medium is a planeparallel slab with a total optical thickness T and that there is no incident field on the boundaries, that is,
The factor ϵ in front of Q^{*}(τ) ensures that Stokes I is of order unity in the interior of the medium.
For the perturbation analysis, it is algebraically simpler to use what we refer to here as the (KQ) representation introduced in Frisch (2007), described below.
2.2. The (KQ) representation of the radiation field
The (KQ) representation makes use of a factorization of the Rayleigh phase matrix involving spherical irreducible tensors for polarimetry introduced by Landi Degl’Innocenti (1984). If a polarized radiation field is represented by its Stokes parameters, the elements of the full 4 × 4 Rayleigh phase matrix may be written
where Ω′ and Ω are the directions of the incident and scattered beams and the spherical irreducible tensors for polarimetry. The index i refers to the Stokes parameters, with i = 0, 1, 2, 3 corresponding to I,Q,U,V, respectively. The indices K and Q take the values K = 0, 1, 2 with −K ≤ Q ≤ K. Explicit expressions of all the can be found in Landi Degl’Innocenti & Landolfi (2004, Table 5.6, p. 211). They depend on the reference direction chosen to reckon the linear polarization. Here the reference direction is chosen parallel to the atmospheric limb (γ = 0 in Landi Degl’Innocenti & Landolfi 2004, Fig. 5.14, p. 196). For Rayleigh scattering, Stokes V is decoupled from the other Stokes parameters and if the linearly polarized field is cylindrically symmetric it can be represented by Stokes I and Stokes Q. Equation (7) reduces then to
with K = 0, 2, i, j = 0, 1, and
It is easy to verify that the matrix P(μ, μ′) can be written as
where the symbol T refers to the transposed matrix and
The factorization of P(μ, μ′) in Eq. (10) leads us to write the Stokes vector I(τ, μ)=(I(τ, μ),Q(τ, μ)) as
If the primary source is unpolarized, that is,
it can be written as
It is easy to verify that the radiative transfer equation for ℐ(τ, μ) has the form
The boundary conditions for ℐ(τ, μ) are the same as those given in Eq. (6) for I(τ, μ). The source vector is given by
where
The phase matrix Ψ(μ) is defined by
By construction Ψ(μ) is a symmetric matrix and an even function of μ. In explicit form,
Its integral over all directions is
We stress that the second element of the diagonal is smaller than unity. This property plays an important role in the perturbation analysis.
In component notation, Eq. (12) becomes
that is,
where and are the two components of the vector ℐ(τ, μ). The polarization is carried by the component .
Because the source term 𝒮(τ) depends only on optical depth, it satisfies a convolution type integral equation. For a semiinfinite medium, it may be written as
where the matrix kernel is defined by
The kernel decreases exponentially as τ goes to infinity and its integral over τ is
The integral equation for 𝒮(τ) is the starting point for various efficient iterative numerical methods of solution, which are of course not limited to a uniform ϵ (for references see Nagendra 2014).
We note that Eq. (21) can be generalized to describe the four Stokes parameters, transforming them into a ninecomponent vector, or a sixcomponent vector, if one considers linear polarization only (Frisch 2007), and also to handle multidimensional problems (Anusha & Nagendra 2011). The (KQ) decompositions are very convenient for numerical work. since the spherical tensors for polarimetry are absorbing the dependence of source terms on the ray directions.
3. The interior and boundary layer variables
For monochromatic scattering, the thermalization length l_{eff} scales as . Different techniques lead to this result: for example, random walks arguments, estimation of the roots of the socalled characteristic equation (see e.g., Sobolev 1963, p. 296; Ivanov 1973, p. 101), and asymptotic analysis of the integral equation of the source function (Frisch & Frisch 1977). This scaling law is not modified by Rayleigh scattering, since Stokes Q is entirely under the control of Stokes I, which behaves as the scalar intensity (e.g., FaurobertScholl et al. 1997).
To ensure that the scale of variation of the radiation field in the interior of the medium is controlled by the value of ϵ, we assume an effectively thick medium, i.e., T ≫ τ_{eff}, where τ_{eff} is the thermalization length in optical depth units. We also assume 𝒬^{*}(τ) is a slowly varying function of the optical depth.
Following standard asymptotic approaches, we introduce a rescaled optical depth , and a rescaled optical thickness defined by
where
We assume that inside the medium
To describe the variations of the radiation field in boundary layers, the variable must be stretched in such a way that the scale of variation of the radiation field is of order unity. This stretching must be performed in the direction perpendicular to the boundary. Here, we are dealing with a 1D medium with two boundaries, one at τ = 0 and one at τ = T. For the boundary at τ = 0, the stretched variable s is defined as
For the boundary at τ = T, the variable s should be defined as . The introduction of a stretched variable to describe a boundary layer is a classical technique for both space and timedependent problems (e.g., Cole 1968; Larsen & Keller 1974; Frisch & Frisch 1977; Bender & Orszag 1978; Bardos et al. 1983).
We write the radiation field in the boundary layer as
For simplicity we refer to as the boundary layer field.
The boundary layer and interior fields are coupled by two conditions. At the surface of the medium, say at τ = 0, in the direction of incoming rays, one has
where ℐ^{ext}(μ) is a given incident radiation. There is a similar condition at τ = T.
In the interior, one must recover the interior solution. This implies that the boundary layer solution goes to zero for s → ∞, that is,
We now consider the interior and boundary layer fields separately. It should be pointed out that the socalled interior region becomes closer to the surface when ϵ increases.
4. The interior radiation field
When performing the asymptotic analysis, it is convenient to write the radiative transfer equation as
where the operator 𝒪 is defined by
with V(μ)=(V_{1}(_{μ}),_{V}_{2}(μ)) an arbitrary 2D vector.
The interior field, now simply written , satisfies the rescaled radiative transfer equation
We look for solutions of the form
The notation 𝒪(η^{n}), with n being a number, signifies of order of η^{n} ^{1}. Regrouping all the terms with the same power of η, we obtain a hierarchy of equations:
and so on.
Equation (37) has nontrivial solutions, that is, solutions which are not identically zero, as we now show. We first remark on the fact that the integral term in Eq. (34) is independent of μ. Hence the vector should also be independent of μ. We deduce from the integral of Ψ(μ) given in Eq. (20) that the second component of is necessarily zero. The solutions of Eq. (37) have thus the form
where , a function of the optical depth only, is an arbitrary function of . The ensemble of such solutions is known as the kernel of the operator 𝒪.
We now consider Eq. (38). Since the kernel of the operator 𝒪 is not empty, the operator 𝒪 is not trivially invertible. Equations such as Eqs. (38) and (39) require a solvability condition. Known as the Fredholm alternative^{2}, which states that a necessary condition for a solution to exist is that the righthand side is orthogonal to the kernel of the adjoint operator. Here the operator 𝒪 is selfadjoint, therefore it has the same kernel as 𝒪. The solvability condition for Eq. (38) is thus
where the symbol ⋅ stands for the scalar product. This equation is satisfied since the μdependent term is an odd function of μ. Therefore Eq. (38) has a solution and it is easy to verify that it has the form
where the function is arbitrary at this stage.
We now consider Eq. (39). We infer from Eq. (37) that the first term on the righthand side is zero. Using Eq. (42) we can write
The solvability condition for this equation is
Performing the integration and the scalar product, one finds the solvability condition:
This equation provides a diffusion equation for the isotropic and unpolarized leading term of the interior solution. It is identical to the diffusion equation for scalar monochromatic scattering (Larsen & Keller 1974; Frisch & Frisch 1977; Bardos et al. 1983). The boundary conditions for this equation are established in Sect. 6, devoted to the asymptotic matching of the interior and boundary layer solutions. The solution of Eq. (43) has the form
The function is at this stage an arbitrary function of . A nonzero second component, denoted , is needed because the first component contains an anisotropic term, , which is an even function of μ. Inserting Eq. (46) into Eq. (43) and using Eqs. (19) and (20) one finds
The function provides the leading term of the component .
Hence, in the interior of the medium, Stokes I is of order unity and is given by , while Stokes Q is of order ϵ, given by
Absorbing the factor ϵ into the rescaled optical depth , Eq. (48) can be rewritten as Q(τ, μ)≃(1 − μ^{2})(∂^{2}I/∂τ^{2})/3.
Now considering the functions and , by continuing the hierarchy of equations to order η^{3}, one can establish that satisfies the homogeneous diffusion equation
The boundary conditions for Eq. (49) are derived in Sect. 6. The function appears in the emergent radiation in terms of order η^{2} (see Sect. 7) and the function in terms of order η^{3}. The asymptotic expansion not being pushed beyond the order η^{2}, the function has not been investigated.
5. The boundary layer field
According to Eq. (33), the boundary layer field satisfies the radiative transfer equation
The primary source term does not appear in this equation because it is incorporated into the equation for the interior. We stress that the domain of the variable s is [0, ∞[.
The boundary condition at s = 0, written in Eq. (31), reduces to
if ℐ^{ext}(μ) = 0. In this case, the role of the incident field is played by the opposite of the interior solution. One looks for solutions which go to zero at infinity (see Eq. (32)).
We now assume that has an expansion of the form
Inserting Eq. (52) into Eq. (50), we see that for k = 0, 1 satisfies the equation
The boundary conditions at s = 0 and at infinity must be satisfied at all orders in η, that is, one must have
Equations (53) and (54) also hold for k = 2 because (see Eq. (70)).
Equation (53), together with Eq. (54), defines a diffuse reflection problem, i.e., conservative scattering in a semiinfinite medium illuminated from outside by a given radiation field. This problem has been investigated in detail by Chandrasekhar (1960) and has also motivated many studies in the field of mathematics (e.g., Golse et al. 1988) in the framework of linear and nonlinear Boltzmann equations. The exact expressions of the emergent radiation field obtained by Chandrasekhar and an exact expression of the radiation field at infinity established below are essential ingredients for the perturbation analysis. All the results for the boundary layer are now given for the components I_{l} and I_{r} of the polarized field, related to the Stokes parameters I and Q by
As shown in Chandrasekhar (1960, p. 256), the two components of the emergent radiation field may be written as
Here and , μ ∈ [ − 1, 0] are the components of the incident field, and q and c are known constants (see Table A.1).
In the scalar case, the radiation field at infinity becomes isotropic, converging exponentially fast to a constant (e.g., Bardos et al. 1983; Frisch & Bardos 1981). For Rayleigh scattering, the radiation field also becomes isotropic and in addition unpolarized, that is, the I_{l} and I_{r} components become equal and independent of μ. They are given by
Equation (58), which is as far as we know a new result, can be demonstrated using different techniques.
One of them is based on the secondorder moment,
The total flux, defined by
is zero for diffuse reflection. The flux is constant because diffuse reflection describes a conservative problem and it is zero at infinity because the radiation field becomes isotropic, hence it is zero everywhere. The multiplication of Eq. (53) by μ followed by its integration over μ shows that the first derivative of J^{(2)}(s) is zero, being equal to F(s)/4π. Hence J^{(2)}(s) is a constant and one may write
The value of can thus be deduced from Eqs. (56) and (57).
It is possible to understand Eq. (58) by considering the solution of the Milne problem by Chandrasekhar (1946). The components of emergent radiation field may be written as
where F is the total radiative flux. Since F is a constant and the incident radiation at the surface is zero, one has the identity
We now consider diffuse reflection. First we assume that the incident radiation is unpolarized and independent of μ, that is, I_{l}(0,μ) = I_{r}(0,μ) = I^{inc}, μ ∈ [ − 1, 0]. In this case, the solution of Eq. (53) is isotropic, unpolarized, and independent of the optical depth, that is, I_{l}(s,μ) = I_{r}(s,μ) = I^{inc} for all values of s and μ ∈ [ − 1, +1]. Since Eq. (63) is an identity, we can set . Thus, if the incident radiation is independent of μ,
Equations (56) and (57) show that the functions H_{l}(μ) and H_{r}(μ) should be associated respectively to the l and r components of the radiation field. The two terms inside the integral describe the transformation of the incident radiation due to the scatterings inside the semiinfinite medium. This process does not depend on the incident field. Hence, Eq. (64) holds for any incident radiation field and leads to Eq. (58).
In the scalar case, the argument given here for Rayleigh scattering leads to
(e.g., Frisch & Bardos 1981).
6. Matching the interior and boundary layer solutions
The boundary conditions needed to solve the diffusion equations satisfied by and (see Eqs. (45) and (49)) are still missing. These conditions are deduced in this section from the constraint that the boundary layer field goes to zero at each order in η for s → ∞, the value at infinity being given by Eq. (58), where we set
Here and are the (lr) components of the interior field. The (lr) and (KQ) components are related by
These equations are readily derived from Eqs. (22) and (55). One can verify that in the direction normal to the surface, that is, for μ = 1, the polarization is zero since I_{l} = I_{r}.
At zeroth order, and . Hence, according to Eq. (67), . Since i(0)_{0}, the constants q and c, and the functions H_{l}(μ) and H_{l}(μ) are positive, the condition that the integral in Eq. (58) is zero leads to the boundary condition
for the interior solution, and to
for the boundary layer solution (see Eq. (54)). The solution of Eq. (53) is thus identically zero:
At first order, we deduce from Eqs. (42) and (67)
The condition that the integral in Eq. (58) is zero leads to
where the constant L^{p} is given by
The subscript p stands for polarization. It can be shown that
where is the common asymptotic value for τ → ∞ of the generalized Hopf functions q_{1}(τ) and q_{2}(τ) (Domke 1971). The constant L^{p} plays for Rayleigh scattering the role played in the scalar case by the constant value at infinity of the Hopf function q(τ). For the Hopf function q(∞) ≃ 0.7104 (e.g., Kourganoff 1952; Landi Degl’Innocenti & Landolfi 2004, p. 670).
The equation for with the boundary condition can be solved by the method of variation of parameters for any choice of the inhomogeneous term . Here, to deal with an explicit and simple expression of , we assume that is a constant q^{*}. Therefore is symmetric with respect to the middle of the slab. Using the boundary conditions i_{0}(0) = 0 and , we readily obtain
We recall that is the rescaled optical thickness of the slab. It is much larger than unity, since the total thickness of the slab has been assumed to be much larger than the characteristic scale τ_{eff}. Therefore at ,
The same values hold at the other boundary. The function increases from zero at the surface to a value close to q^{*} in the middle of the slab and decreases to zero towards the other boundary.
We note here that
for any choice of . This result is readily derived from the diffusion equation for written in Eq. (45) and from the boundary condition i_{0}(0) = 0.
Using the boundary condition in Eq. (72) and the symmetry of with respect to the middle of the slab, one finds for the solution of Eq. (49)
Therefore at ,
For an unpolarized and constant primary source (q^{*}, 0),
The function decreases from the boundary at to the middle of the slab, where it becomes exponentially small and increases again towards the other surface. We now have all the elements to determine the emergent radiation field.
7. The emergent radiation field and polarization rate
The emergent radiation field is given by the sum
where is the surface value of the interior solution and the contribution from the boundary layer. To calculate this latter term, we must replace in Eqs. (56) and (57) the incident field by the opposite of the interior solution taken at τ = 0, at each order in η. First we must express the interior field in terms of its (lr) components with the help of Eq. (67).
7.1. Surface value of the interior solution
At zeroth order, the interior solution is unpolarized and , hence .
At first order, the interior solution is unpolarized. Equations (71) and (72) lead to
The μdependent term contributes to the emergent radiation field, but not the constant term. Indeed, as pointed out in Sect. 5, if the incident field in Eqs. (56) and (57) is unpolarized and independent of μ, the emergent field is equal to the incident field. The latter being the opposite of the interior field, the interior and boundary layer contributions cancel each other in Eq. (81).
At second order, Eq. (46) combined with Eq. (67) and the boundary condition i_{0}(0) = 0 leads to
The terms q^{*}(0) and c_{2}(0) are independent of μ and appear in both Eqs. (83) and (84), and therefore will not contribute to the emergent radiation field. Keeping the relevant terms and using the expression of given in Eq. (47), we use for the interior field the expressions
The secondorder derivative has an explicit expression given in Eq. (77). We prefer not to use it, to keep track of the origin of the term.
7.2. Calculation of the emergent radiation
We now introduce Eqs. (82), (85), and (86) into Eqs. (56) and (57) and perform the integrations over μ′. We introduce the notation
At first order in η, the contribution of the boundary layer to the emergent radiation may be written as
with
The integrals involving the functions f_{lr}(μ, μ′) and f_{rl}(μ, μ′) are easy to express in terms of the moments of the functions H_{l}(μ) and H_{r}(μ). In the Appendix A we show how to calculate the other integrals. We thus obtain
At second order in η, the boundary layer contribution may be written as
where
After some algebra described in Appendix A, we obtain
where
and
The constants are moments of H_{l}(μ) and H_{r}(μ) defined by
Their numerical values for m = 0, 1, 2, 3 are listed in Table A.1.
For μ = 1, the emergent field should be unpolarized at all orders in η. Using the relations (Chandrasekhar 1960)
one easily verifies that
Adding the interior solution and the boundary layer contribution, one finds
where
and
Similar expressions hold for the emergent intensity at the surface τ = T, the derivatives of i_{0}(τ) at τ = 0, being replaced by their values at τ = T. The terms independent of H_{l}(μ) and H_{r}(μ) in Eqs. (91) and (96) disappear when adding the interior and boundary layer contributions.
All the terms in Eqs. (103)–(106) are of order unity, therefore in Eq. (102) the first term is truly of order η and the second one of order η^{2}. The perturbation expansion is consistent with the well known result that the radiation field at the surface of a semiinfinite or effectively thick medium is of order if the primary source is of order ϵ.
We also note that the term of order η has the same centertolimb variation as the emergent radiation of the Milne problem, the factor 3F/8π being replaced by one half the surface value of the derivative of .
7.3. The polarization rate
We now examine the effect of the destruction rate ϵ on the polarization rate, defined as in Chandrasekhar (1960), by
If we keep only the term of order η, we recover Chandrasekhar’s law for conservative scattering, namely
A nonzero value of ϵ leads to a decrease in the polarization rate. The effect is contained in the terms of order η^{2} given in Eqs. (105) and (106). To simplify the notation, we introduce
and the functions
With these notations, Eq. (107) becomes
The function has an explicit expression in terms of , leading to (see Eq. (79)).
The polarization decrease due to a nonzero ϵ is given by the difference p(μ)−p_{0}(μ). Subtracting p_{0}(μ)=h^{−}(μ)/h^{+}(μ) from Eq. (112), one obtains
If , Eqs. (76), (77), and (111) lead to . In this case
We now compare the predictions of the perturbation analysis with the results of a direct numerical simulation of the full radiative transfer equation.
7.4. Perturbation analysis versus numerical simulations
The polarized radiative transfer equation for the vector ℐ(τ, μ) written in Eq. (15) is solved numerically with a PALI (Polarized Accelerated Lambda Iteration) code (Nagendra et al. 1998, 1999) for several values of ϵ. The atmospheric model is a slab with a constant and unpolarized primary source term 𝒬^{*}(τ)=Q^{*}(τ)=(1, 0). For each choice of ϵ, the optical thickness of the slab satisfies to ensure that the thermalization is reached at the middle of the slab. The τgrid is logarithmic with 5 points per decade starting at τ = 10^{−4} and the μ gridpoints have a GaussLegendre distribution with 43 points between μ = 7.638 10^{−4} and μ = 0.9992.
Figure 2 is devoted to the polarization rate p(μ) as defined in Eq. (107). Solid lines show the numerical results obtained with the PALI code and dashed lines the perturbation values predicted by Eq. (112). The effects of a finite value of ϵ appear only above ϵ ≥ 10^{−4} and the predictions of the perturbation analysis are in excellent agreement with the numerical results for ϵ ≤ 10^{−3}. For larger values of ϵ, the polarization rate is underestimated by the perturbation analysis. For ϵ = 10^{−2}, the error on the polarization rate is about 30% for μ around 0.1. For larger values of μ, the asymptotic analysis becomes unreliable because the position of the change of sign of Stokes Q is not correctly evaluated. For ϵ = 10^{−1}, the expansion parameter is too large for the asymptotic expansion to be valid. Also in this figure, we note the fast decrease of the polarization rate for ϵ between 10^{−2} and 10^{−1}.
Fig. 2. Polarization rate in percent. The atmospheric model is an optically thick slab with an unpolarized and uniform primary source and no incident radiation. Dashed lines: perturbation expansion. Solid lines: numerical results calculated with a PALI code based on Nagendra et al. (1999). From top to bottom, Chandrasekhar’s law, log_{10}ϵ = −4, −3, −2, −1. For log_{10}ϵ = −1, only the numerical solution is shown. The Chandrasekhar limit is graphically indistinguishable from the ϵ = 10^{−4} numerical results. 
The asymptotic results for ϵ around 10^{−2} could probably be improved by continuing the asymptotic expansion to order η^{3} = ϵ^{3/2}. However, at this order, the boundary layer term is not a solution of Eq. (53). There is an additional inhomogeneous term depending on , making the construction of the solution more complicated.
The change of sign of the polarization rate observed in Fig. 2 for ϵ = 10^{−2} and ϵ = 10^{−1} is due to a change of sign of (see Eq. (22)). The radiative transfer equation for ℐ(τ, μ) shows that
where is the second component of 𝒮(τ). Remarking that is much larger than , Eq. (16) leads to
The righthand side (without the factor (1 − ϵ)), usually denoted (Trujillo Bueno 2001), measures the anisotropy of the radiation field intensity. It is positive if the intensity of the “vertical” rays (those with ) is larger than the intensity of the “horizontal” rays (those with ). At the surface, is positive. For a uniform primary source, the change of sign occurs at optical depths less than 10 and approaches the surface when ϵ increases. It becomes observable via the polarization rate if it occurs at optical depths smaller than unity.
Fig. 3. Deviation from Chandrasekhar’s law, p_{0}(μ)−p(μ), scaled by the expansion parameter . The dashed and dasheddotted lines show the contributions of and , and the first and second components of ℐ_{2}(0, μ), respectively. The figure is drawn ϵ = 10^{−2}, i.e., η = 10^{−1}. 
In Fig. 3, we show the ratio (p_{0}(μ)−p(μ))/η given by Eq. (114), for a uniform and unpolarized primary source term. This ratio is almost independent of η. The solid line corresponds to Eq. (114). The contributions of the first and second component of ℐ_{2}(τ, μ), the secondorder term in the interior solution, are shown separately with the dashed and the dasheddotted line, respectively. We can observe that they yield similar contributions, the second component q_{2}(0), carrying the polarization, having a somewhat larger contribution, and that they have a similar μdependence. Figure 3 suggests that Eq. (114) can be approximated by
For an unpolarized but τdependent primary source, the polarization rate p(μ) has a similar variation, since the τdependence affects only the constant factor g(0) in Eq. (113).
Figure 4 shows the centertolimb variation of the emergent value of Stokes I for different values of ϵ.
Fig. 4. Centertolimb variation of Stokes I. Same atmospheric model as for the other figures. Dashed lines: perturbation expansion. Solid lines: numerical results calculated with a PALI code. From bottom to top, Chandrasekhar’s law, log_{10}ϵ = −3, −2, −1. For log_{10}ϵ = −1, only the numerical solution is shown. 
One can observe that the perturbation expansion provides a very good fit to this variation for values of ϵ up to 10^{−2}. For 10^{−1}, it becomes invalid. For ϵ = 1, the centertolimb variation is equal to one for all values of μ, since the thermal source of our model is uniform.
8. Polarized and anisotropic source term
We now assume that the primary source term Q^{*}(τ, μ) is polarized and μdependent, but is still a slowly varying function of τ. In general, the decomposition Q^{*}(τ, μ)=A(μ)𝒬^{*}(τ), necessary to define a radiation field ℐ(τ, μ) solution of a radiative transfer equation in which the source term depends only on the optical depth (see Eq. (15)), will not be satisfied. As pointed out in Ivanov et al. (1997) and Frisch (2007), the technique to recover a decomposition similar to Eq. (12) is to introduce a diffuse field I^{d}(τ, μ) defined by
where I^{p}(τ, μ) is solution of the transfer equation
Simple algebra shows that the diffuse field satisfies Eq. (5) with a primary source term
which may be written as
The components q^{*}(τ) and of the vector can be calculated with Eq. (120). The (KQ) components of the diffuse field satisfy Eq. (15), the source term being given by Eq. (16).
The leading term and the function remain unchanged and the emergent radiation at order η is still given by Eqs. (103) and (104). Changes appear at order η^{2}. The second component of becomes
Hence, the square bracket in Eq. (85), giving the surface value of the lcomponent, contains the additional term
For the rcomponent, the additional term in the square bracket of Eq. (86) is
These terms contribute to the emergent intensity at order η^{2}. For the lcomponent, one must add
inside the square bracket of Eq. (105). For the rcomponent,
must be added inside the square bracket of Eq. (106). Using Eq. (100), one can verify that . The polarization remains at zero in the direction of the symmetry axis, at it should.
The polarization rate given in Eq. (112) is of course modified by these new terms. The square brackets in the numerator and in the denominator will contain the functions and , respectively, defined by
9. External incident radiation field
Here we assume that an external polarized radiation field is incident on the surface at τ = 0 and that there is no primary source of photons inside the slab. One expects the emergent radiation to be given at leading order by Eqs. (56) and (57), where the external field I^{ext}(μ) plays the role of the incident radiation. We use the asymptotic expansion described in the preceding sections to establish this result and to find a correction term of order . All the fields are represented by their (l, r) components.
The interior diffusion equation is not modified by the presence of an external radiation field, the boundary layers at τ = 0 and τ = T can also be described by diffuse reflection problems, but changes appear in the boundary conditions at τ = 0 for the interior and boundary layer equations.
The boundary condition at τ = 0 for the total field (interior plus boundary layer) is now
For the zeroth and firstorder terms, this yields
Higher orders also satisfy Eq. (130). Here we consider only the zeroth and firstorder terms.
At zerothorder in η, , where is the solution of Eq. (45) with . To calculate the emergent radiation, one must add the contribution from the interior solution, given by [I^{int}(0, μ)]_{0}, and the contribution from the boundary field given by Eqs. (56) and (57) in which the incident radiation is
The contribution from the term (i_{0}(0),i_{0}(0))/2 is canceled by the interior solution. One thus recovers Eqs. (56) and (57) in which the incident field is simply the exterior field. These equations may be written as
where
The functions depending on μ and μ′ are defined in Eq. (87).
The value of i_{0}(0) does not enter into the leading term of the emergent radiation but to evaluate the effects of a nonzero ϵ, one needs the solution of the homogeneous diffusion equation for . The boundary condition at is readily derived from the constraint that goes to zero at infinity. Inserting Eq. (131) into Eq. (58), one finds
The righthand side comes from the first term on the righthand side of Eq. (131) and the lefthand side from the second term. We recall that the value at infinity is equal to the incident field, if the latter is unpolarized and independent of μ. At , the boundary condition is still . The solution of the diffusion equation for is
with i_{0}(0) given in Eq. (134). The derivative of at τ = 0 can be approximated by
We now turn to the firstorder term. We can apply the results obtained in the preceding sections for a slab with no incident radiation, since the exterior field is zero at this order (see Eq. (130)). The components and are given by Eqs. (103) and (104). Adding the zerothorder contribution to the firstorder one, calculated with the approximate value of given in Eq. (136), and using the expression of i(0)_{0} given in Eq. (134), one obtains
with
and
If we assume that the external incident field is a beam in the direction μ_{0}, there is no need to integrate over μ′ in Eq. (137).
We recover here an approximation proposed in Eq. (71) of Domke (1973), which is also based on an analysis of the radiation field in the deep layers, albeit somewhat less systematic than ours. By comparison with an exact numerical calculation of Lenoble (1970) for ϵ = 10^{−2}, Domke has correctly evaluated this approximation to have an error of order ϵ.
10. Concluding remarks
This paper shows, apparently for the first time, how a standard asymptotic expansion method for scalar transport equations can be generalized to a transport equation in which the phase function is a matrix. Applied to nonconservative Rayleigh scattering, it leads to explicit expressions for the polarized emergent radiation field. These are given in Eqs. (102)–(106) for a medium with an unpolarized primary source of photons and no incident field and in Eq. (137) for a medium with an incident radiation but no primary source. These expressions are not exact results, but being constructed with a systematic expansion in powers of , with ϵ the destruction probability per scattering, their accuracy can be precisely defined and they can be used to check the accuracy of numerical schemes and to optimize the choice of discretization grids. They also provide ways of testing the very small ϵ regime, for which good accuracy is hard to reach.
Comparisons between the predictions of the perturbation expansion and numerical solutions calculated with a PALI code are carried out for a uniform primary source to investigate the domain of validity of the expansion. They show that the centertolimb variations of the emergent values of Stokes I and Stokes Q follow the Chandrasekhar’s law for the conservative scattering given in Eqs. (103) and (104) for small values of ϵ, up to 10^{−4}. For larger values, deviations from the Chandrasekhar’s law become significant. The polarization rate is correctly predicted for values of ϵ up to 10^{−3}, but underestimated for larger values (see Fig. 2). For Stokes I, the centertolimb variations are correctly predicted for ϵ up to 10^{−2} (see Fig. 4).
The perturbation expansions given in Eqs. (102)–(106) allow one to check separately the accuracy on the two components of the radiation field and to verify that the constant factor is independent of μ and close to for an unpolarized and uniform primary source . For this test, an adequately small value of ϵ should be chosen; a value certainly less than 10^{−4}. The same kind of test can be performed for the secondorder derivative which should be equal to −3q^{*}. To isolate the secondorder term in Eq. (102), one must subtract the firstorder term, calculated with a very small value of ϵ, say ϵ = 10^{−8}. In Eqs. (105) and (106), the derivative can be replaced by . The results given here concerning the domain of validity of the perturbation expansion are derived from tests made with a uniform primary source. They should hold for any spacedependent source provided it shows no significant variations over optical depths of order unity. Special asymptotic techniques have been developed for rapidly varying media (e.g., Bensoussan et al. 1979).
We indicate how the perturbation expansion method presented here for a onedimensional medium could be generalized to the multidimensional case. The basic idea of introducing a largescale description of the radiation field in the interior of the medium holds independently of the dimension of the medium. For multidimensional problems, the scaling factor has simply to be applied to all the components of the space variable. The stretching of the rescaled space variable in the direction perpendicular to the surface has to be carried out at each point on the surface of the medium, transforming locally the boundary layer into a semiinfinite plane parallel medium (e.g., Larsen & Keller 1974; Frisch 1988). This rescaling procedure has been applied up to now only to scalar transport problems, but there is nothing preventing it from being extended to multidimensional Rayleigh scattering.
For a multidimensional geometry, it is possible to construct a multipolar (KQ) decomposition for the three Stokes parameters I, Q, and U that are needed to describe the radiation field in the absence of cylindrical symmetry. This decomposition, proposed in Anusha & Nagendra (2011), could be used to construct an asymptotic expansion for the interior solution in much the same way as in the 1D case. The boundary layer problem remains a conservative diffuse reflection problem, but with an incident field without cylindrical symmetry. Chandrasekhar (1960, p. 259) shows how to treat this general case. For multidimensional scalar problems, Larsen & Keller (1974) show how to match the interior and boundary layer solutions. For a polarized radiation field, a detailed analysis of the interior and the boundary layer solution should show how to establish matching conditions.
Acknowledgments
The author is very grateful to Dr. Sampoorna (IIA, Bangalore) for having kindly provided a PALI (Polarized Accelerated Lambda Iteration) code for the numerical calculation of Rayleigh scattering, to Dr. V. Bommier, Prof. C. Bardos, and Prof. F. Golse for stimulating discussions, and to an anonymous referee for constructive remarks.
References
 Anusha, L. S., & Nagendra, K. N. 2011, ApJ, 726, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Bardos, C., Santos, R., & Sentis, R. 1983, Trans. Am. Math. Soc., 284, 617 [CrossRef] [MathSciNet] [Google Scholar]
 Bender, C. M., & Orszag, S. A. 1978, Advanced Mathematical Methods for Scientists and Engineers (New York: McGrawHill) [Google Scholar]
 Bensoussan, A., Lions, J. L., & Papanicolaou, G. 1979, Publ. RIMS Kyoto Univ., 15, 53 [CrossRef] [MathSciNet] [Google Scholar]
 Chandrasekhar, S. 1946, ApJ, 104, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1960, Radiative Transfer (New York: Dover Publications) [Google Scholar]
 Cole, J. 1968, Perturbation Methods in Applied Mathematics (Waltham, MA: Blaisdell) [Google Scholar]
 de Rooij, W. A., Bosma, P. B., & van Hooff, J. P. C. 1989, A&A, 226, 347 [NASA ADS] [Google Scholar]
 Domke, H. 1971, Sov. Astron., 15, 266 [NASA ADS] [Google Scholar]
 Domke, H. 1973, Sov. Astron., 17, 81 [NASA ADS] [Google Scholar]
 FaurobertScholl, M., & Frisch, H. 1989, A&A, 322, 896 [NASA ADS] [Google Scholar]
 FaurobertScholl, M., Frisch, H., & Nagendra, K. N. 1997, A&A, 219, 338 [NASA ADS] [Google Scholar]
 Fluri, D. M., & Stenflo, J. O. 1999, A&A, 341, 902 [NASA ADS] [Google Scholar]
 Frisch, H. 1988, in Saasfee 18th Advanced Course, eds. R. P. Kudritzki, H. W. Yorke, & H. Frisch, 339 [Google Scholar]
 Frisch, H. 1999, ASSL, 243, 97 [NASA ADS] [Google Scholar]
 Frisch, H. 2007, A&A, 476, 665 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Frisch, H., & Bardos, C. 1981, JQSRT, 26, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Frisch, U., & Frisch, H. 1977, MNRAS, 181, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Golse, F., Perthame, B., & Sulem, C. 1988, Arch. Ration. Mech. Anal., 103, 81 [CrossRef] [Google Scholar]
 Ivanov, V. V. 1973, Transfer of Radiation in Spectral Lines (Washington Government Printing Office, NBS Spec: Publ), 385 [Google Scholar]
 Ivanov, V. V. 1990, Sov. Astron., 34, 621; transl. from Astron. Zh. 67, 1233 [NASA ADS] [Google Scholar]
 Ivanov, V. V., Kasaurov, A. M., & Loskutov, V. M. 1996, A&A, 307, 332 [NASA ADS] [Google Scholar]
 Ivanov, V. V., Grachev, S. I., & Loskutov, V. M. 1997, A&A, 318, 315 [NASA ADS] [Google Scholar]
 Kostogryz, N. M., & Berdyugina, S. V. 2015, A&A, 575, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kourganoff, V. 1952, Basic Methods in Transfer Problems. Radiative Equilibrium and Neutron Diffusion (London: Oxford University Press) [Google Scholar]
 Landi Degl’Innocenti, E. 1984, Sol. Phys., 91, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Landi Degl’Innocenti, E., & Landolfi, M. 2004, Polarization in Spectral Lines (Dordrecht: Kluwer Academic Publishers) [Google Scholar]
 Larsen, E. W., & Keller, J. B. 1974, J. Math. Phys., 15, 75 [CrossRef] [Google Scholar]
 Lenoble, J. 1970, JQSRT, 10, 533 [NASA ADS] [CrossRef] [Google Scholar]
 Nagendra, K. N. 2014, in ASP Conf. Ser., eds. K. N. Nagendra, J. Stenflo, Q. Qu, & M. Sampoorna, 489, 179 [NASA ADS] [Google Scholar]
 Nagendra, K. N., Frisch, H., & FaurobertScholl, M. 1998, A&A, 332, 610 [NASA ADS] [Google Scholar]
 Nagendra, K. N., Paletou, F., Frisch, H., & FaurobertScholl, M. 1999, ASSL, 243, 127 [NASA ADS] [Google Scholar]
 Štěpàn, J., & Bommier, V. 2007, A&A, 468, 7975 [Google Scholar]
 Siewert, C. E., Kelley, C. T., & Garcia, R. D. M. 1981, J. Math. Anal. Appl., 84, 509 [CrossRef] [Google Scholar]
 Sobolev, V. V. 1963, A Treatise on Radiative Transfer (Princeton, NJ: Von Nostrand Company); russian original Moscow (1956) [Google Scholar]
 Trujillo Bueno, J. 2001, in ASP Conf. Ser., ed. M. Sigwarth, 236, 161 [NASA ADS] [Google Scholar]
 Trujillo Bueno, J., & Shchukina, N. 2009, ApJ, 694, 1364 [Google Scholar]
Appendix A: Calculation of the emergent radiation field
Moments of H_{l}(μ) and H_{r}(μ) and some numerical constants.
We give here some details about the derivation of Eqs. (103)–(106), giving the emergent radiation field for a medium with no incident radiation and an unpolarized primary source.
We first show how to calculate integrals with the form
Other integrals appearing in Eqs. (56) and (57) are readily expressed in terms of the moments of H_{l}(μ) and H_{r}(μ) defined in Eq. (99). We proceed recursively, starting from α = 1, that is, from
Inserting μ′ = μ′ + μ − μ, we can write
The functions H_{l}(μ) and H_{r}(μ) satisfy the nonlinear integral equations
Hence,
The constants and are the zerothorder moments of H_{l}(μ) and H_{r}(μ) corresponding to k = 0 (see Eq. (99)). The numerical values of the four first moments are given in Table A.1.
By taking the limit μ → ∞ in Eqs. (A.4) and (A.5), one obtains two simple and useful relations, namely
Combined with Eqs. (A.6) and (A.7), they lead to
We now consider
Inserting
we can write
Thus
Using the expressions of A_{l}(μ) and A_{r}(μ) given in Eqs. (A.9) and (A.10), we find
Finally we consider
Proceeding as above, we write
Hence
Inserting the expressions of B_{l, r}(μ) given in Eqs. (A.15) and (A.16), we obtain
We now turn to the calculation of the functions and defined in Eqs. (89) and (90). They may be written as
Combined with Eqs. (A.15) and (A.16), these expressions become
They can be simplified with help of the following relations:
Equation (A.26) is the fundamental relation between q and c. Equations (A.27) and (A.28) express the property that the emergent field is equal to the incident field if the latter is unpolarized and independent of μ, and Eq. (A.29) is obtained by combining Eqs. (A.27) and (A.28) and taking Eq. (A.26) into account. The combination of Eq. (A.8) with Eqs. (A.27) and (A.28) leads to two other useful relations, namely
Employing Eq. (A.30) to express in Eq. (A.24) the secondorder moments in terms of firstorder moments, we obtain
The identity in Eq. (A.26) leads to the simple result
The expression of given in Eq. (A.25) is readily simplified with the help of Eq. (A.29). One obtains
Equations (A.32) and (A.33) appear in the text in Eq. (91) and lead to Eqs. (103) and (104) for the emergent radiation field at first order in η.
We now consider the functions and defined in Eqs. (94) and (95). They can be written as
Inserting the expressions of A_{l}(μ) and C_{l}(μ) given in Eqs. (A.9) and (A.20), we obtain
Expressing in terms of firstorder moments with the help of Eq. (A.30) and using Eq. (A.27) to simplify the second square bracket, we obtain
For , replacing in Eq. (A.35) A_{r}(μ) and C_{r}(μ) by the expressions given in Eqs. (A.10) and (A.21), we obtain
Proceeding as above, we express in terms of firstorder moments with the help of Eq. (A.30) and then use Eq. (A.28) to simplify the first square bracket. We obtain
The expressions of and in Eqs. (A.37) and (A.39) are used to construct the emergent radiation field at order η^{2} given in Eqs. (105) and (106). The functions and in Eq. (96) are the opposite of the square brackets in Eqs. (A.37) and (A.39).
All Tables
All Figures
Fig. 1. Functions H_{l}(μ) and H_{r}(μ). Solid line: H_{l}(μ); dashed line: H_{r}(μ). The functions H_{l}(μ) and H_{r}(μ) are defined for all positive values of μ. For μ → ∞, H_{l}(μ) goes to infinity as and H_{r}(μ) to . The curves have been plotted with the data given in Chandrasekhar (1960). 

In the text 
Fig. 2. Polarization rate in percent. The atmospheric model is an optically thick slab with an unpolarized and uniform primary source and no incident radiation. Dashed lines: perturbation expansion. Solid lines: numerical results calculated with a PALI code based on Nagendra et al. (1999). From top to bottom, Chandrasekhar’s law, log_{10}ϵ = −4, −3, −2, −1. For log_{10}ϵ = −1, only the numerical solution is shown. The Chandrasekhar limit is graphically indistinguishable from the ϵ = 10^{−4} numerical results. 

In the text 
Fig. 3. Deviation from Chandrasekhar’s law, p_{0}(μ)−p(μ), scaled by the expansion parameter . The dashed and dasheddotted lines show the contributions of and , and the first and second components of ℐ_{2}(0, μ), respectively. The figure is drawn ϵ = 10^{−2}, i.e., η = 10^{−1}. 

In the text 
Fig. 4. Centertolimb variation of Stokes I. Same atmospheric model as for the other figures. Dashed lines: perturbation expansion. Solid lines: numerical results calculated with a PALI code. From bottom to top, Chandrasekhar’s law, log_{10}ϵ = −3, −2, −1. For log_{10}ϵ = −1, only the numerical solution is shown. 

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.