© H. Frisch 2019

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 late-type 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 photo-ionization 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 semi-infinite medium can be solved exactly. He showed in particular that the emergent radiation field for the Milne problem, which describes conservative scattering in a semi-infinite medium with a given total flux F, no incident radiation, and diffuse reflection problems, which describe conservative scattering in a semi-infinite 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).

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 semi-infinite 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 mono-kinetic 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 free-transition, 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 semi-infinite 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 time-dependent 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 multi-dimensional 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 η² = ϵ 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 plane-parallel 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 two-component 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,

(1)

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/λ⁴ for Rayleigh scattering. The primary source term given by the Planck function is unpolarized.

Introducing the monochromatic optical depth,

(2)

and the parameter

(3)

we can rewrite the radiative transfer equation as

(4)

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

(5)

with ϵ being a constant and Q^*(τ) an unpolarized primary source. We assume that the medium is a plane-parallel slab with a total optical thickness T and that there is no incident field on the boundaries, that is,

(6)

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

(7)

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

(8)

with K = 0, 2, i, j = 0, 1, and

(9)

It is easy to verify that the matrix P(μ, μ′) can be written as

(10)

where the symbol T refers to the transposed matrix and

(11)

The factorization of P(μ, μ′) in Eq. (10) leads us to write the Stokes vector I(τ, μ)=(I(τ, μ),Q(τ, μ)) as

(12)

If the primary source is unpolarized, that is,

(13)

it can be written as

(14)

It is easy to verify that the radiative transfer equation for ℐ(τ, μ) has the form

(15)

The boundary conditions for ℐ(τ, μ) are the same as those given in Eq. (6) for I(τ, μ). The source vector is given by

(16)

where

(17)

The phase matrix Ψ(μ) is defined by

(18)

By construction Ψ(μ) is a symmetric matrix and an even function of μ. In explicit form,

(19)

Its integral over all directions is

(20)

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

(21)

that is,

(22)

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 semi-infinite medium, it may be written as

(23)

where the matrix kernel is defined by

(24)

The kernel decreases exponentially as τ goes to infinity and its integral over τ is

(25)

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 nine-component vector, or a six-component vector, if one considers linear polarization only (Frisch 2007), and also to handle multi-dimensional 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 so-called 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., Faurobert-Scholl 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

(26)

where

(27)

We assume that inside the medium

(28)

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

(29)

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 time-dependent 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

(30)

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

(31)

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,

(32)

We now consider the interior and boundary layer fields separately. It should be pointed out that the so-called 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

(33)

where the operator 𝒪 is defined by

(34)

with V(μ)=(V₁(_μ),_V2(μ)) an arbitrary 2D vector.

The interior field, now simply written , satisfies the rescaled radiative transfer equation

(35)

We look for solutions of the form

(36)

The notation 𝒪(ηⁿ), 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:

(37)

(38)

(39)

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

(40)

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², which states that a necessary condition for a solution to exist is that the right-hand side is orthogonal to the kernel of the adjoint operator. Here the operator 𝒪 is self-adjoint, therefore it has the same kernel as 𝒪. The solvability condition for Eq. (38) is thus

(41)

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

(42)

where the function is arbitrary at this stage.

We now consider Eq. (39). We infer from Eq. (37) that the first term on the right-hand side is zero. Using Eq. (42) we can write

(43)

The solvability condition for this equation is

(44)

Performing the integration and the scalar product, one finds the solvability condition:

(45)

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

(46)

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

(47)

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

(48)

Absorbing the factor ϵ into the rescaled optical depth , Eq. (48) can be rewritten as Q(τ, μ)≃(1 − μ²)(∂²I/∂τ²)/3.

Now considering the functions and , by continuing the hierarchy of equations to order η³, one can establish that satisfies the homogeneous diffusion equation

(49)

The boundary conditions for Eq. (49) are derived in Sect. 6. The function appears in the emergent radiation in terms of order η² (see Sect. 7) and the function in terms of order η³. The asymptotic expansion not being pushed beyond the order η², the function has not been investigated.

5. The boundary layer field

According to Eq. (33), the boundary layer field satisfies the radiative transfer equation

(50)

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

(51)

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

(52)

Inserting Eq. (52) into Eq. (50), we see that for k = 0, 1 satisfies the equation

(53)

The boundary conditions at s = 0 and at infinity must be satisfied at all orders in η, that is, one must have

(54)

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 semi-infinite 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

(55)

As shown in Chandrasekhar (1960, p. 256), the two components of the emergent radiation field may be written as

(56)

(57)

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

(58)

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 second-order moment,

(59)

The total flux, defined by

(60)

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⁽²⁾(s) is zero, being equal to F(s)/4π. Hence J⁽²⁾(s) is a constant and one may write

(61)

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

(62)

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

(63)

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 μ,

(64)

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 semi-infinite 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

(65)

(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

(66)

Here and are the (lr) components of the interior field. The (lr) and (KQ) components are related by

(67)

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)₀, 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

(68)

for the interior solution, and to

(69)

for the boundary layer solution (see Eq. (54)). The solution of Eq. (53) is thus identically zero:

(70)

At first order, we deduce from Eqs. (42) and (67)

(71)

The condition that the integral in Eq. (58) is zero leads to

(72)

where the constant L^p is given by

(73)

The subscript p stands for polarization. It can be shown that

(74)

where is the common asymptotic value for τ → ∞ of the generalized Hopf functions q₁(τ) and q₂(τ) (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 and , we readily obtain

(75)

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 ,

(76)

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

(77)

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.

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)

(78)

Therefore at ,

(79)

For an unpolarized and constant primary source (q^*, 0),

(80)

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

(81)

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

(82)

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 leads to

(83)

(84)

The terms q^*(0) and c₂(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

(85)

(86)

The second-order 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

(87)

At first order in η, the contribution of the boundary layer to the emergent radiation may be written as

(88)

with

(89)

(90)

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

(91)

At second order in η, the boundary layer contribution may be written as

(92)

(93)

where

(94)

(95)

After some algebra described in Appendix A, we obtain

(96)

where

(97)

and

(98)

The constants are moments of H_l(μ) and H_r(μ) defined by

(99)

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)

(100)

one easily verifies that

(101)

Adding the interior solution and the boundary layer contribution, one finds

(102)

where

(103)

(104)

and

(105)

(106)

Similar expressions hold for the emergent intensity at the surface τ = T, the derivatives of i₀(τ) 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 η². The perturbation expansion is consistent with the well known result that the radiation field at the surface of a semi-infinite 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 center-to-limb 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

(107)

If we keep only the term of order η, we recover Chandrasekhar’s law for conservative scattering, namely

(108)

A nonzero value of ϵ leads to a decrease in the polarization rate. The effect is contained in the terms of order η² given in Eqs. (105) and (106). To simplify the notation, we introduce

(109)

(110)

and the functions

(111)

With these notations, Eq. (107) becomes

(112)

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₀(μ). Subtracting p₀(μ)=h⁻(μ)/h⁺(μ) from Eq. (112), one obtains

(113)

If , Eqs. (76), (77), and (111) lead to . In this case

(114)

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⁻⁴ and the μ grid-points have a Gauss-Legendre distribution with 43 points between μ = 7.638 10⁻⁴ 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⁻⁴ and the predictions of the perturbation analysis are in excellent agreement with the numerical results for ϵ ≤ 10⁻³. For larger values of ϵ, the polarization rate is underestimated by the perturbation analysis. For ϵ = 10⁻², 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⁻¹, 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⁻² and 10⁻¹.

The asymptotic results for ϵ around 10⁻² could probably be improved by continuing the asymptotic expansion to order η³ = ϵ^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⁻² and ϵ = 10⁻¹ is due to a change of sign of (see Eq. (22)). The radiative transfer equation for ℐ(τ, μ) shows that

(115)

where is the second component of 𝒮(τ). Remarking that is much larger than , Eq. (16) leads to

(116)

The right-hand 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.

In Fig. 3, we show the ratio (p₀(μ)−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 ℐ₂(τ, μ), the second-order term in the interior solution, are shown separately with the dashed and the dashed-dotted line, respectively. We can observe that they yield similar contributions, the second component q₂(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

(117)

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 center-to-limb variation of the emergent value of Stokes I for different values of ϵ.

One can observe that the perturbation expansion provides a very good fit to this variation for values of ϵ up to 10⁻². For 10⁻¹, it becomes invalid. For ϵ = 1, the center-to-limb 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

(118)

where I^p(τ, μ) is solution of the transfer equation

(119)

Simple algebra shows that the diffuse field satisfies Eq. (5) with a primary source term

(120)

which may be written as

(121)

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 η². The second component of becomes

(122)

Hence, the square bracket in Eq. (85), giving the surface value of the l-component, contains the additional term

(123)

For the r-component, the additional term in the square bracket of Eq. (86) is

(124)

These terms contribute to the emergent intensity at order η². For the l-component, one must add

(125)

inside the square bracket of Eq. (105). For the r-component,

(126)

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

(127)

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

(128)

For the zeroth- and first-order terms, this yields

(129)

(130)

Higher orders also satisfy Eq. (130). Here we consider only the zeroth- and first-order terms.

At zeroth-order 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, μ)]₀, and the contribution from the boundary field given by Eqs. (56) and (57) in which the incident radiation is

(131)

The contribution from the term (i₀(0),i₀(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

(132)

where

(133)

The functions depending on μ and μ′ are defined in Eq. (87).

The value of i₀(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

(134)

The right-hand side comes from the first term on the right-hand side of Eq. (131) and the left-hand 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

(135)

with i₀(0) given in Eq. (134). The derivative of at τ = 0 can be approximated by

(136)

We now turn to the first-order 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 zeroth-order contribution to the first-order one, calculated with the approximate value of given in Eq. (136), and using the expression of i(0)₀ given in Eq. (134), one obtains

(137)

with

(138)

and

(139)

If we assume that the external incident field is a beam in the direction μ₀, 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⁻², 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 center-to-limb 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⁻⁴. For larger values, deviations from the Chandrasekhar’s law become significant. The polarization rate is correctly predicted for values of ϵ up to 10⁻³, but underestimated for larger values (see Fig. 2). For Stokes I, the center-to-limb variations are correctly predicted for ϵ up to 10⁻² (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⁻⁴. The same kind of test can be performed for the second-order derivative which should be equal to −3q^*. To isolate the second-order term in Eq. (102), one must subtract the first-order term, calculated with a very small value of ϵ, say ϵ = 10⁻⁸. 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 space-dependent 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 one-dimensional medium could be generalized to the multi-dimensional case. The basic idea of introducing a large-scale description of the radiation field in the interior of the medium holds independently of the dimension of the medium. For multi-dimensional 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 semi-infinite 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 multi-dimensional Rayleigh scattering.

For a multi-dimensional 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 multi-dimensional 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

Appendix A: Calculation of the emergent radiation field

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

(A.1)

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

(A.2)

Inserting μ′ = μ′ + μ − μ, we can write

(A.3)

The functions H_l(μ) and H_r(μ) satisfy the nonlinear integral equations

(A.4)

(A.5)

Hence,

(A.6)

(A.7)

The constants and are the zeroth-order 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

(A.8)

Combined with Eqs. (A.6) and (A.7), they lead to

(A.9)

(A.10)

We now consider

(A.11)

Inserting

(A.12)

we can write

(A.13)

Thus

(A.14)

Using the expressions of A_l(μ) and A_r(μ) given in Eqs. (A.9) and (A.10), we find

(A.15)

(A.16)

Finally we consider

(A.17)

Proceeding as above, we write

(A.18)

Hence

(A.19)

Inserting the expressions of B_l, r(μ) given in Eqs. (A.15) and (A.16), we obtain

(A.20)

(A.21)

We now turn to the calculation of the functions and defined in Eqs. (89) and (90). They may be written as

(A.22)

(A.23)

Combined with Eqs. (A.15) and (A.16), these expressions become

(A.24)

(A.25)

They can be simplified with help of the following relations:

(A.26)

(A.27)

(A.28)

(A.29)

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

(A.30)

Employing Eq. (A.30) to express in Eq. (A.24) the second-order moments in terms of first-order moments, we obtain

(A.31)

The identity in Eq. (A.26) leads to the simple result

(A.32)

The expression of given in Eq. (A.25) is readily simplified with the help of Eq. (A.29). One obtains

(A.33)

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

(A.34)

(A.35)

Inserting the expressions of A_l(μ) and C_l(μ) given in Eqs. (A.9) and (A.20), we obtain

(A.36)

Expressing in terms of first-order moments with the help of Eq. (A.30) and using Eq. (A.27) to simplify the second square bracket, we obtain

(A.37)

For , replacing in Eq. (A.35) A_r(μ) and C_r(μ) by the expressions given in Eqs. (A.10) and (A.21), we obtain

(A.38)

Proceeding as above, we express in terms of first-order moments with the help of Eq. (A.30) and then use Eq. (A.28) to simplify the first square bracket. We obtain

(A.39)

The expressions of and in Eqs. (A.37) and (A.39) are used to construct the emergent radiation field at order η² 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