Variance reduction method for particle transport equation in spherical geometry
X Blanc C Bordin G Kluth G Samba
math.NA

Dust coagulation regulated by turbulent clustering in protoplanetary disks
    Takashi Ishihara1,Naoki Kobayashi2,Kei Enohata2,Masayuki Umemura3, Kenji Shiraishi4
    December 30, 2023
=======================================================================================




  This article is devoted to the design of importance sampling method for the Monte Carlo simulation of a linear
  transport equation. This model is of great importance in the simulation of inertial confinement fusion
  experiments. Our method is restricted to a spherically symmetric idealized design : an outer sphere emitting 
  radiation towards an inner sphere, which in practice should be thought of as the hohlraum and the fusion
  capsule, respectively. We compute the importance function as the solution of the corresponding stationary adjoint problem. Doing
  so, we have an important reduction of the variance (by a factor 50 to 100), with a moderate increase of
  computational cost (by a factor 2 to 8). 




§ INTRODUCTION

In inertial confinement fusion (ICF) experiments, a small ball of hydrogen (the target) is submitted to intense
radiation by laser beams. These laser beams are either pointed directly to the target (direct drive approach), or pointed
to gold walls of a hohlraum in which the target is located (indirect drive approach, see Figure <ref>). These gold walls heat up, emitting X-rays toward the
target. The outer layers of the target are heated up, hence ablated. By momentum conservation, the inner part of
the target implodes (this is usually called the rocket effect). Hence, the pressure and temperature of the hydrogen
inside the target increase, hopefully reaching the thermodynamical conditions for nuclear fusion. This process is
summarized in Figure <ref>.





The numerical simulation of such an experiment involves many physical phenomena such as hydrodynamics, radiation
transfer, neutronics, etc... In the present article, we focus on the simulation of radiation, that is, the
transmission of the (X-ray) energy to the target. A simplified model for this is the grey radiative transfer
equation:

    {∂_t u + ·∇ u + κ_t u      = κ_s ∫_S^2  u(t,,') 
    k(,',)dΩ' + Q() 
    
    	 u(t=0,,)    = g(,),
    .

where the solution u is the radiation intensity and depends on the time t, the position ∈ℝ^3, the direction of propagation ∈ S^2. The term Q()
represents a source of radiation. In the present case, Q() is a modelling of the emission of X-rays by the
hohlraum walls. Furthermore, κ_t is the total 
cross-section. It satisfies κ_t = κ_a +
κ_s, where κ_a≥ 0 is the absorption
cross-section and κ_s≥ 0 the scattering cross-section. The kernel k(,',) is a
probability density with respect to ' and , that is, k≥ 0 and ∫
k(,',)d' = ∫
k(,',)d=1. Note that we have assumed here that we use units such that the speed of light
is c=1. 


Equation (<ref>) may be simulated using a Monte Carlo method. If so, the probability distribution k
may be interpreted as the probability density associated to the new direction propagation 
for a particle having a shock with initial direction '. 
In Monte Carlo simulations of such situations, variance reduction
methods are important to reduce the statistical noise. Indeed, as the target implodes, hydrodynamic instabilities
develop, which are a source of energy loss. Should this loss be too important, the experiment would be
compromised. Thus, it is important to have a precise numerical description of these instabilities. In the case of
Monte Carlo simulations, this implies a statistical noise as small as possible (at least smaller than the amplitude of
the instabilites). A small variance is particularly important on the target boundary. 

A widely used reduction variance technique in such a situation is the importance sampling method. It may be
  summarized as follows:
  
  
  *  Calculate the importance function (in our case, the solution of the adjoint equation);
  
  *  Use the importance function to modify the transport equation, and apply a Monte Carlo method;
  
  * Calculate the forward intensity from <ref> and <ref>
  
 



Importance sampling is a well-known reduction variance method, which has
  been applied to transport problems in many situations. We refer for instance to the textbooks <cit.>
  <cit.> for a general presentation. The key-point in such a method is the way one computes the
  importance function. If it is solution to the adjoint problem, then one achieves a zero-variance method. However,
solving the adjoint problem is at least as difficult as solving the direct problem at hand. Therefore, many methods
using approximations of the adjoint solution have been developped. This is the spirit of the exponential transform
(see <cit.> and <cit.>). In some situations, a diffusion approximation is used for this
calculation, as for instance in <cit.>. In other situations, discrete ordinates approximation is
preferred <cit.>. The method which is the closest to the one presented here is probably
<cit.>, in which the adjoint equation is formulated as an integral equation, and solved using a
space discretization. An importance difference is, however, that when solving the adjoint problem, the scattering
is neglected in <cit.>. Here, we use the same kind of method, but taking advantage of the
radially symmetric geometry, we are able to take scattering effects into account.



The article is organized as follows: in Section <ref>, we give a rapid presentation of the
Monte Carlo method applied to transport equations, then of the importance sampling method. This method is based on
the computation of an importance function, which is the subject of Section <ref>. In
Section <ref>, we present some numerical experiments, while the appendices contain some
technical result which we do not want to detail in the main body of the article. 




§ MONTE CARLO METHOD FOR TRANSPORT EQUATIONS





 §.§ Natural method



We give in this subsection a short overview of the application of Monte Carlo method applied to transport
equation. More details and mathematical justifications are given in <cit.>. We
concentrate here on practical aspects. 

Considering equation (<ref>), we define the following quantities (we assume here that Q and g
are integrable functions):

    α = ∫_S^2∫_𝒟 g(,)d d,   g(,) = 1/αg(,),


    β = ∫_S^2∫_𝒟 Q()d d = 4π∫_𝒟 Q()d,
        Q() = 1/βQ()

Here, 𝒟 is the spatial domain. Hence, g and Q are probability measures on the
phase space 𝒟× S^2. Of course, we assume that both g and Q are non-negative, which is
physically relevant. 

We first deal with the case Q = 0, and then extend it to the general case. 
We define N independent realizations (X_i(t),Ω_i(t)) of the jump Markov process (X(t),Ω(t)) as
follows:


  * (X_i(0),Ω_i(0)) are drawn independently of each other, following the law g(,)d
  d. 
 To each of
  them is assigned a weight w_i(0)=1/N. 

  * Between jumps, (X_i,Ω_i) follows the characteristics of Equation (<ref>), that is,
  
    {
      
           Ẋ_i(t) = Ω_i(t), 

           Ω̇_i(t) = 0,
      
.
  
which is equivalent to the fact that Ω_i is constant[Note that in curvilinear coordinates, this is not the
  case. For instance, in Section <ref> below, spherical coordinates are used, hence the
  direction μ_i is a non-trivial function of t.] and 

    X_i(t) = X_i(t_0) + (t-t_0)Ω_i.

Moreover, the weight w_i(t) is assumed to satisfy the equation ẇ_i(t) + (κ_t - κ_s)w_i(t) = 0,
that is,

    w_i(t) = w_i(t_0) e^-(κ_t - κ_s)(t-t_0).


  * Time jumps are defined by a Poisson process 𝒩(t) of intensity one as follows: if the process
  𝒩(tκ_s) has a jump at time t, then Ω_i(t) has a jump, and the conditional law of
  Ω_i(t^+) knowing Ω_i(t^-) is given by k(,_i(t^-),)d


It can be proved that such a strategy gives a good approximation of the solution u(t,,) to
(<ref>) in the following sense <cit.>: assume that 

    η_N(t,d,d) = α∑_i=1^N w_i(t) δ_X_i(t),Ω_i(t)(d,d),

then this measure converges narrowly to u(t,,)d d, as N→+∞. 



It remains to include the influence of the source Q(x). For this purpose, we split the time interval into time
steps of equal size Δ t. What follows can easily be generalized to non-constant time steps, but this not our
purpose here. At each time step, we generate M more realizations of another jump Markov process, independently of
the initial ones, as follows:


  * At time mΔ t, we draw M independent couples (X_j^m, Ω_j^m) according to the law
  Q()d d.

  The weight of each particle is w_j^m(mΔ t) = Δ t/M

  * Each of these random variables follow the same evolution as in the preceding case, with positions,
  velocities, weights X_j^m(t),Ω_j^m(t),w_j^m(t), respectively.


Finally, the measure η_N defined by (<ref>) is replaced by (here, we assume that nΔ t ≤ t <
(n+1)Δ t)

    η_N,Δ t(t,d,d) = α∑_i=1^N w_i(t) δ_X_i(t),Ω_i(t)(d,d) +
      β∑_m=0^n ∑_j=1^M w_j^m(t) δ_X_j^m(t),Ω_j^m(t)(d,dΩ).

Here again, this measure narrowly converges to u(t,,)d d, as n,N,M→+∞, (with
nΔ t→ t) where u solves
(<ref>), according to <cit.>. 

As we already pointed out, the time step Δ t may be non-uniform, and, moreover, the number M of particles
generated at each time step may depend on the time step m. 



Finally, we point out that we did not take care about boundary conditions. The boundary conditions we aim at
imposing are either free boundary condition, or imposed incoming flux. In the first case, if a process hits the
boundary, it simply vanishes. In the second case, we write the boundary condition as a source consisting of a Dirac mass supported by the boundary,
therefore including it into the source Q. 



 §.§ Importance sampling



The method of importance sampling is widely used in many applications of Monte Carlo methods, and in particular in
the case of the simulation of transport equations (see for instance
<cit.>). A nice account of this method in the present context can also be
found in <cit.>.

The idea of importance sampling is to introduce an importance function, which we call I(t,,), and which
is assumed to be positive. Instead of applying a Monte Carlo method to compute u, we are going to apply it to the
function

  ũ(t,,) = I(t,,) u(t,,).

A simple computation gives the equation satisfied by ũ:

    ∂_t ũ + ·∇ũ + κ̃_t(t,,)ũ =  ∫_S^2κ̃_s(t,,') ũ (t,,')k̃(t,,',)d' + Q̃(t,,),

with

    κ̃_t(t,,) = κ_t - ∂_t I(t,,)/I(t,,) -
      ·∇ I(t,,)/I(t,,),


    κ̃_s(t,,) = κ_s/I(t,,)∫_S^2k(,,”)I(t,,”)d”,


    k̃(t,,',) = ( ∫_S^2 k(,',”)I(t,,”) d”)^-1 k(,',)I(t,,),

and

    Q̃(t,,) = I(t,,) Q().



Equation (<ref>) has a similar structure as (<ref>), and a Monte Carlo method can easily
be designed to compute an approximation of its solution. 
It should be noted that, although the coefficients κ_t and κ_s where assumed to be constant in
(<ref>), the new coefficients defined by (<ref>) and (<ref>) do
depend on t,  and  since I does. Similarly, k̃ depends on t although k does not, and
Q̃ is now a function of t,,. This is not a problem for Monte Carlo simulations, the only point
is that equation (<ref>) should be modified as follows:

    w̃_i(t) = w̃_i(t_0) exp( -∫_t_0^t (κ̃_t(s,X_i(s),Ω_i) - κ̃_s(s,X_i(s),Ω_i) )ds).

Note that the definition (<ref>) of κ̃_t does not imply that
  κ̃_t≥ 0. This may be a problem when dealing with Monte Carlo simulations of
  (<ref>). However, if I is solution to the
  adjoint equation (see (<ref>) below), then κ̃_t = κ̃_s ≥ 0. In the present work, I is not
an exact solution of (<ref>), but, as it is pointed out in Remark <ref> below, we still have
κ̃_t≥0 and κ̃_s≥ 0.
The importance function I should be chosen in such a way that the variance of the computed approximation of
ũ has a smaller variance than the approximation of u computed with the method described above.





 §.§ Adjoint equation



It is known (see <cit.>) that, in order to have zero variance, the importance function should be
solution to the adjoint equation:

    -∂_t I -∇ I + κ_t I = κ_s ∫_S^2 I(t,,') k(x,,')d'.

A rigorous proof of the above fact may be found in <cit.>, but let us give a simple argument which
indicates that this is indeed the case. We assume that the spatial domain D is a ring between R_0 and
R_1>0:

  𝒟 = { ∈^3,   R_0< || < R_1}.

We consider equation (<ref>) in this domain, with initial condition g=0 and boundary conditions 

  u(t,,) =
  
    1    if  || = R_1,   ·n() <0, 

    0    if  || = R_0,   ·n() <0.
  

Here, n() is the outer normal unit vector to the domain 𝒟 at point . Actually, for the
domain we are studying, n() = /|| if ||=R_1 and n() = -/|| if ||=R_0.
We assume that we are interested in computing the flux on the target || = R_0, that is,

    F(T) = ∫_0^T ∫_||=R_0∫_· n()>0 |· n()| u(t,,)
      d d dt.

This can be done using the above Monte Carlo method. An estimator of the quantity F(T) is then given by the
following:

    F = ∑_j,  |X_j|=R_0,  Ω_j· n(X_j)>0 w_j,

provided that, at each time step, particles are created with |X_j|= R_1, and Ω_j drawn according to
Lambert cosine law (see <cit.>), with initial weights equal to Δ t/M, where M
is the number of particles created at each time step.

Let us now make precise the equation which the importance function I solves: we assume that (<ref>) is satisfied, and that
the following boundary conditions are imposed

    I(t,,) =
      
        0    if  || = R_1,   · n() >0, 
    
        1    if  || = R_0,   · n() >0.


Now, consider a Monte Carlo method applied to ũ, that is, to equation (<ref>). Here, Q̃ = 0, but the boundary condition on ũ is different:

  ũ(t,,) =
  
    I(t,,)    if  || = R_1,   ·n() <0, 

    0    if  || = R_0,   ·n() <0.
  

Hence, the boundary data for ũ amounts to sampling the distribution I(t,,). In particular, a
good choice for the initial weights in such a case is the following:

  w̃_j(t) = Δt/M∫_||=R_1∫_·n()<0 |Ω·n()| I(t,,)dd,

with velocities drawn according to the law |Ω· n()|I(t,,)1_{· n() <0}/∫_· n() <0 |Ω· n()|I(t,,)d. 

First, we note that, since I solves (<ref>), we have κ̃_s = κ̃_t. Hence, the
weight of a particle does not change between shocks. Second,
we point out that all particles go to the target. In order to see this, we assume that it is not the
case. Then, there exists a particle which exits the computation domain through the outer boundary
||=R_1. Denote by  its position and  its direction when it exits. Then, · n()>0, || = R_1, and u(t,,)
I(t,x,) ≠ 0. But this is impossible since the boundary condition satisfied by I is I(t,,Ω) =
0 for such values or  and .

Hence, the estimated flux F̃ is now deterministic:

    F̃ = ∑_j,  |X_j|=R_0,  Ω_j· n(X_j)<0w̃_j = ∫_0^T
      ∫_||=R_1∫_· n()<0|· n() |
      I(t,,)d  d dt.

Finally, multiplying (<ref>) by I, integrating,
and using (<ref>), a simple integration by parts proves that

  ∫_0^T ∫_||=R_1 ∫_·n()<0 |Ω·n()| I(t,,)dd=
  ∫_0^T ∫_||=R_0 ∫_·n()>0 |Ω·n()| u(t,,)dd= F(T).

Hence, if one is able to compute I solution to (<ref>) with boundary conditions
(<ref>), and to compute the integral on the right-hand side of (<ref>), then
we have an exact evaluation of the quantity F(T).



§ COMPUTATION OF THE IMPORTANCE FUNCTION



As it was stated in the previous section, computing the solution to the adjoint equation allows to have an
importance function such that the result of importance sampling computations has zero variance. However, solving this equation is at least as difficult as solving (<ref>). Therefore, a tractable
approximation of the solution to this equation should be sought in order to be used as an importance function.
It is in general
not possible to compute it exactly, but we will see that in particular situations, simplified expressions may
be 
derived which give good approximation of the solution to the adjoint problem. First, in
Section <ref>, we review the work of <cit.>, in which an analytic expression was
derived for the importance function in dimension one. Then, in Section <ref>, we extend this
analysis to the spherically symmetric case. In such a case, an analytic expression is no longer valid, but a
numerical solution may be computed with the characteristics method, if κ_t, κ_s and k do not vary
in space. If they do vary in space, we do not know for now how to generalize the calculations of Section <ref>. As mentionned in
Section <ref> below, a possible way to address this question is to compute numerically the
solution. Doing so, one should be careful to have a good balance between precision (which allows for a significant
improvement of the variance) and computational cost.



 §.§ Analytic expression : one-dimensional planar case



In his Phd thesis <cit.>, J.-M Depinay developed a method to compute an importance function in case
of stationary transport equation in slab geometry. This Section is not directly related to what we
  do in the spherically symmetric case. It is only a simple example in which an explicit computation of the
  importance function I is available. In Section <ref>, we generalize (to some extent) this
  approach. In the planar case, the unknown u of (<ref>) is
assumed to depend only on one space variable x, where = (x,y,z). Moreover, the structure of the equation
implies (see  <cit.>) that u depends
on  only through the scalar product μ = · (1, 0, 0). Therefore, using these
notations, Equation (<ref>) reads

∂_t u + μ∂_x u + κ_t u  =κ_s ∫_-1^1  k(x,μ',μ)u(x,μ') dμ' + Q(x).

Looking for stationary solutions, and assuming that the source is zero (Q = 0), this reduces to

    μ∂_x u + κ_t u  =κ_s ∫_-1^1  k(x,μ',μ)u(x,μ') dμ'.

Hence, the corresponding adjoint problem reads

    -μ∂_x I + κ_t I  =κ_s ∫_-1^1  k(x,μ,μ')I(x,μ') dμ'.

If one assumes that I is of the form

    I(x,μ) = exp(Kx)Φ_K(μ),

then (<ref>) implies that Φ_K(μ) satisfies the equation 

    Φ_K(μ) = κ_s/κ_t-Kμ∫_-1^1 k(x,μ,μ')ϕ_K(μ')dμ',

where the parameter K is chosen such that 

    ∫_-1^1 Φ_K(μ)dμ = 1.

It can be proved that such (<ref>) always has a unique solution K (see <cit.>),
if k∈ L^∞. Moreover, assuming that k is constant (that is, k=1/2), then (<ref>) reduces to

  Φ_K(μ) = 1/2κ_s/κ_t-Kμ.

Hence, K is the unique solution of equation (<ref>), which reads

  1/2 ∫_-1^1 κ_s/κ_t - Kμdμ= 1.

Let us insert ũ(x,v) = u(x,v) I(x,v) into Equation (<ref>). We have the following modified equation for ũ

    μ∂_x ũ + κ̃_t ũ = ∫_-1^1 κ̃_s(x,μ')k̃(x,μ',μ)ũ(x,μ') dμ'.

with the modified parameters

    κ̃_s    = κ_t - Kμ, 
    κ̃_t    = κ_t - Kμ, 
    
    k̃(μ',μ)    = 1/2Φ_K(μ)/Φ_K(μ')κ_s/κ_t-Kμ'.

These expressions are (<ref>), (<ref>), (<ref>), adapted to the
particular case of 1D planar (slab) geometry.

Let us point out two important things here: first, the above importance function I(x,μ) is a solution to the
adjoint equation, but does not in general satisfy appropriate boundary conditions. Therefore, it might result in
poor variance reduction in some situation. However, the tests in <cit.> indicate very good efficiency
of the method for a case in which a detector is placed far away from an emitting source. Second, this kind of
solution is related to those exhibited in <cit.> (see also <cit.>). Such solutions are
eigenvectors of the (adjoint) transport operator, corresponding to the largest possible eigenvalues. This is why
they play an important role here. 

Finally, although the above derivation is done with a stationary transport equation, an implicit time scheme leads,
at each time step, to solving a stationary transport equation. Therefore, although we do not use an
implicit time scheme to solve our transport equation (this would imply additional difficulties that go beyond the
scope of the present work, see <cit.>, the review paper <cit.> and the references therein) the use of the corresponding adjoint
solution in the time-dependent case may prove efficient. This is the strategy we are going to apply in the
spherically symmetric case. 


 §.§ The spherical case


Now, we want to compute an importance function in the case of spherical geometry. In such a case, the unknown u
is assumed to depend on  only through r = ||, which in turn implies that it depends on  only through
μ = ·/||. As a consequence, Equation (<ref>) becomes <cit.>:

    ∂_t u + μ∂_r u + 1-μ^2/r∂_μ u + κ_t u = κ_s∫_-1^1  k(r,μ',μ)u(t,r,μ') dμ' + Q(r).

Here, μ = cosθ, where θ is the angle formed by the radial direction  and the direction .


[scale=1]

(2,0) arc (0:90:2) ;

(0,0) – (2,0) ;

(0,0) – (0,2) ;

(0,0) node[below]0 ;

[dashed] (0,0) – (1.414213562,1.41421362) ;
	
(1.414213562,1.41421362) node[right] ;

(1.414213562,1.41421362) – (2.5,2.5) ;
(1.414213562,1.41421362) node +;

(1.6,1.6) arc (0:90:0.2) node[above right] θ;

[->] (1.414213562,1.41421362) – (1.414213562,2.5) ;
(1.414213562,2.5) node[left] ;



As it has been done in the general case, we consider an importance sampling function I(t,r,μ) and define
ũ(t,r,μ)=u(t,r,μ)I(t,r,μ) in Equation (<ref>). Thus, we have the following equation for
ũ: 

    ∂_t ũ+ μ∂_r ũ + 1-μ^2/r∂_μũ + κ̃_tũ = ∫_-1^1κ̃_s(x,μ') k̃(x,μ',μ)ũ(x,μ') dμ' + Q̃,

with 










Q̃(t,r,μ) = Q(r)I(t,r,μ),   


    κ̃_s(μ) = κ_s1/I(t,r,μ)∫_-1^1 I(t,r,μ”) k(r,μ,μ”) dμ”,


    κ̃_t(μ) = κ_t - ( ∂_t I + μ∂_r I + 1-μ^2/r∂_μ I )1/I(t,r,μ) ,


    k̃(μ',μ) = I(t,r,μ)k(μ',μ)/∫_-1^1 I(t,r,μ”) k(r,μ',μ”) dμ”


The setting we are going to use is the following: we want to reproduce the geometry of an ICF experiment, with a
good statistical convergence on the boundary of the target (or at the ablation front, which is even better). In
order to do so, we assume that the computation domain is

    D = { r,    R_0≤ r ≤ R_1},

where R_0 is the radius of the target (or of the ablation front), and R_1>R_0 is the outer boundary of the
domain. In all the following, R_1 is assumed to be fixed, whereas R_0 = R_0(t) may be a function of time,
reflecting the dynamics of the implosion. An incoming flux is imposed on the outer boundary, while the quantity we
want to compute is the outgoing flux at r=R_0. Therefore, Equation (<ref>) is set with the boundary
conditions:

    u(r,μ) =
      
        1     if    r=R_1,  μ <0, 
    
        0     if    r=R_0,  μ>0.


Equivalently, equation (<ref>) is set with the boundary conditions:

    ũ(r,μ) =
      
        I(r,μ)     if    r=R_1,  μ <0, 
    
        0     if    r=R_0,  μ>0.


Finally, we assume that k is constant, although this is not a limitation in our strategy:

  k(r,μ',μ) = 1/2.






  §.§.§ Solution of the adjoint equation

To find I, we are going to solve the following adjoint equation 

    -μ∂_r I - 1-μ^2/r∂_μ I + κ_t I = κ_s ⟨ I ⟩ ,

where

    ⟨ I ⟩ = 1/2∫_-1^1 I(r,μ)dμ.

We use a stationary approximation for I, although equation (<ref>) is not stationary.
This proves sufficient in the tests below, but this approximation will need to be assessed in
  the presence of material motion, as we point out in the conclusion below. The reason for
this is twofold: first, some of the numerical tests we are going to use are in fact stationary, and second, even in
the case of a non-stationary situation, the only dependence on time in the model is that of R_0, which is assumed
in fact to be constant in each time step of the simulation. Therefore, at each time step, the stationary importance
function should give a good variance reduction.

In order to have a zero variance on the inner boundary R_0, the boundary conditions for I should be the
following:

    {
    I(R_0,μ)     = 1   if μ <0 
    
    I(R_1,μ)    = 0   if μ > 0,
    .

Considering S(r) = κ_s⟨ I ⟩ as a source, Equation (<ref>) reads :

    -μ∂_r I - 1-μ^2/r∂_μ I + κ_t I = S(r).

with the same boundary conditions. We are going to use the method of characteristics to solve (<ref>)-(<ref>). In order to do so,
we change variables, setting x = rμ and y = r√(1-μ^2). The domain is thus (see Figure <ref>):

    { (x,y)∈[-R_1 ; R_1]×[0 ; R_1 ]  | R_0^2 ≤ x^2+y^2 ≤R_1^2
    } = { (x,y)∈ℝ×ℝ^+  | R_0^2 ≤ x^2+y^2 ≤R_1^2}.

We denote by J the new unknown, that is, I as a function of the new variables:

  I(r,μ) = J(rμ,r√(1-μ^2)),   and  J(x,y) = I(√(x^2+y^2), x/√(x^2+y^2)).

Equation (<ref>) becomes 

    -∂_x J + κ_t J = S(√(x^2+y^2)).

This is equivalent to 

    -∂_x(J(x,y)e^-κ_t x)  = S(√(x^2+y^2))e^-κ_t x

with the boundary conditions

    {    J(x,y) = 1   if  x < 0,    x^2+y^2=R_0^2, 
        J(x,y) = 0   if  x > 0,    x^2+y^2=R_1^2.
    .

To solve (<ref>)-(<ref>), we split the domain into three different parts
(see Figure <ref>): 

    𝒟_1     = { (x,y)   |   x<0,   y<R_0,    R_0< √(x^2+y^2)<R_1 }, 
    𝒟_2     = { (x,y)  |   x>0,   y<R_0,    R_0< √(x^2+y^2)<R_1}, 
    𝒟_3     = { (x,y)  |    y>R_0,    R_0< √(x^2+y^2)<R_1}.

 
In each domain, we integrate the equation along the characteristics and apply the boundary conditions. Thus, we
have, in
domain 𝒟_1, 

J(x,y)= exp( κ_t ( x+√(R_0^2-y^2)) ) + ∫_x^-√(R_0^2-y^2) S(√(s^2+y^2) ) exp(κ_t (x-s))ds

In domain 𝒟_2,

J(x,y)= ∫_x^√(R_1^2-y^2) S(√(s^2+y^2) ) exp(κ_t (x-s))ds,

and in domain 𝒟_3,

J(x,y)= ∫_x^√(R_1^2-y^2) S(√(s^2+y^2) ) exp(κ_t (x-s))ds.

Collecting all these results, the importance function I(r,μ) reads

    I(r,μ) =  1_{μ <
        -√(1-R_0^2/r^2)}exp(κ_t(rμ + √(R_0^2-r^2+r^2μ^2)) ) 
     + ∫_rμ^R(r,μ) S(√(s^2+r^2-r^2μ^2)) exp(κ_t (rμ-s))ds,

where R(r,μ) = - √(R_0^2-r^2+r^2μ^2) 1_{μ < -√(1-R_0^2/r^2)}
+ √(R_1^2-r^2+r^2μ^2) 1_{μ > -√(1-R_0^2/r^2)}.
In the above formulae (and in the sequel), we use the notation 1 to indicate the step
  function: for any m∈,
  
    1_{μ< m} =
    
      1     if  μ< m, 

      0      if  μ≥m.
    
  





  §.§.§ Integral equation on S

Recalling that S(r) = κ_s⟨ I ⟩, and integrating
(<ref>) with respect to μ, one obtains an integral equation on ϕ
defined by

  ϕ(r) = r S(r) = κ_s r ⟨I⟩.

This equation reads

    ϕ(r) = κ_s/2∫_-1^-√(1-R_0^2/r^2) rexp(κ_t(rμ + √(R_0^2-r^2+r^2μ^2))
    )dμ
    + κ_s/2∫_-1^1 ∫_rμ^R(r,μ) rϕ(√(s^2+r^2-r^2μ^2))/√(s^2+r^2-r^2μ^2)exp(κ_t (rμ-s))dsdμ.

It is possible to compute exactly the integrals with respect to μ, using the exponential integral function (see
<cit.>, and Appendix <ref> below). This gives

    ϕ( r ) = κ_s/4   ( 1/κ_t[ exp(κ_t θ) ]_-r+R_0^-√(r^2-R_0^2) + (R_0^2 - r^2 )[ κ_t ( κ_t θ) - exp( κ_t θ)/θ]_-r+R_0^-√(r^2-R_0^2)) 
        + κ_s/2∫_R_0^R_1ϕ(r') [ ( κ_t (-√(r^2-R_0^2) - √(r'^2-R_0^2)) ) - ( κ_t(-| r -r' | ))   ] dr'.

Here, (x) is the exponential integral defined by

    (x) = ∫_-x^∞exp( -t)/tdt.




  §.§.§ Numerical computation of I


A priori, it is not possible to solve this integral equation exactly. However, it is possible to solve it numerically. For this purpose, we introduce a mesh to discretize the space variable r. In all the numerical
examples we are going to give, we use a uniform mesh, although this is not essential. Thus, we use the following
notations: let N_r be a positive integer, and define Δ r = R_1-R_0/N_r. Note however that in some
of the cases treated below, R_0 depends on t. In such cases, we use a mesh independent of t, with Δ r = R_1/N_r.


    ∀ 0≤ j ≤ N_r-1,    M_j = [r_j-1/2,r_j+1/2],    r_j-1/2 = R_0 + jΔ r,   
      r_j+1/2 = R_0 +  (j+1)Δ r, 
     r_j = R_0+(j+1/2)Δ r.

We then use a piecewise constant approximation of ϕ, defining

  ϕ(r) = ∑_j=0^N_r-1 ϕ_j 1_M_j(r),   hence   ∫_M_j ϕ= |M_j|ϕ_j =
  Δr ϕ_j.

Inserting this into (<ref>), and using a piecewise constant approximation for all the
functions appearing in the integrals, we infer

  ϕ_j = b_j + κ_s/2 ∑_i=0^N_r-1 ϕ_i Δr[ (-κ_t √(r_j^2 - R_0^2) -
      κ_t√(r_i^2 - R_0^2)) - (-κ_t |r_j-r_i| )],

where

    b_j = κ_s/4( 1/κ_t[ exp(κ_t θ) ]_-r_j+R_0^-√(r_j^2-R_0^2) + (R_0^2 - r_j^2 )[ κ_t ( κ_t θ) - exp( κ_t θ)/θ]_-r_j+R_0^-√(r_j^2-R_0^2))

Hence, we are lead to the following linear system satisfied by ϕ:

    (I_d - A )ϕ = b,
 
where I_d is the identity matrix and A is defined by:

    A_ij = κ_s/2Δ r[
     ( κ_t (-√(r_i^2-R_0^2) - √(r_j^2-R_0^2)) ) - ( κ_t(-| r_i -r_j | ))],
 
and the right-hand side b is defined by (<ref>). This formula is valid only in the case i≠
j. If i=j, the singularity of  at the origin does not allow for the use of (<ref>). In order
to compute them, we note that H=1 is the unique solution of the system

    {    -μ∂_r H - 1-μ^2/r∂_μ H + κ_t H = κ_s ⟨ H ⟩ + κ_t-κ_s 
       H(R_0,μ)  = 1  , μ<0
       H(R_1,μ)  = 1  , μ>0
    .
 
Now, applying the method of characteristics as above to this system, we have an equation for H similar to
(<ref>), with an additional term due to the boundary condition at r=R_0:

    H(r,μ) = exp( κ_t( rμ + √( r^2μ^2 - r^2 + R_0^2)))1_{μ < -√(1-R_0^2/r^2)}
    
    +  exp( κ_t( rμ + √( r^2μ^2 - r^2 + R_1^2)))1_{μ > -√(1-R_0^2/r^2)}
     + ∫_rμ^R(r,μ)(Ψ(√(s^2+r^2-r^2μ^2))/√(s^2+r^2-r^2μ^2)+ κ_t-κ_s)exp(κ_t (rμ-s) )ds.

Moreover, integrating with respect to μ, we also have a relation similar to (<ref>):

    Θ( r ):=κ_sr⟨ H⟩ = κ_s/4( 1/κ_t[ exp(κ_t θ) ]_-r+R_0^-√(r^2-R_0^2) + (R_0^2 - r^2 )[ κ_t ( κ_t θ) - exp( κ_t θ)/θ]_-r+R_0^-√(r^2-R_0^2)) 
    
    + κ_s/4( 1/κ_t[ exp(κ_t θ) ]^r+R_1_-√(r^2-R_0^2)+√(R_1^2-R_0^2) + (R_1^2 - r^2 )[ κ_t ( κ_t θ) - exp( κ_t θ)/θ]^r+R_1_-√(r^2-R_0^2)+√(R_1^2-R_0^2)) 
    
    + κ_s/2∫_R_0^R_1(Ψ(r') + (κ_t - κ_s )r') [ ( κ_t (-√(r^2-R_0^2) - √(r'^2-R_0^2)) ) - ( κ_t(-| r -r' | ))   ] dr'

Using the fact that H=1, and assuming a piecewise constant approximation of Θ, we may assume

    Θ_j = κ_s r_j.

Inserting this into (<ref>), we find

κ_s( r_j - ∑_i=0^N_r-1A_jir_i) = b_j + c_j + ∑_i=0^N_r-1 A_jid_i ,

where the coefficients A_ij and b_j are defined as above by (<ref>) and
(<ref>). Here, the coefficients c_j are given by 

c_j  = κ_s/4  ( 1/κ_t [ exp(κ_t θ)
  ]^r_j+R_1_-√(r_j^2-R_0^2)+√(R_1^2-R_0^2) + (R_1^2 - r_j^2 )[ κ_t
    ( κ_t θ) - exp( κ_t θ)/θ
  ]^r_j+R_1_-√(r_j^2-R_0^2)+√(R_1^2-R_0^2)),

and 

d_j  = (κ_t - κ_s )r_j.

Hence, the vector (Θ_j) satisfies the equation

(I_d - A )Θ= b + c + Ad,

Hence, using (<ref>),

  r_jκ_s = b_j + c_j + κ_t ∑_i=0^N_r-1 A_jir_i.

Finally,

    A_jj = r_jκ_s - b_j - c_j/κ_t r_j - ∑_i≠ j A_jir_i/r_j.

We have an expression of the diagonal coefficients A_jj in terms of the off-diagonal ones.




After solving (<ref>),
we use ϕ to define an approximation of the importance function I(r,μ):

    I(r,μ) = 1_{μ < -√(1-R_0^2/r^2)}exp( κ_t( rμ + √( r^2μ^2 - r^2 + R_0^2))) 
     +  ∑_i=0^N_r-1ϕ_i∫_rμ^R(r,μ)1_M_i(√(s^2+r^2-r^2μ^2)) exp(κ_t (rμ-s) )/√(s^2+r^2-r^2μ^2)ds.




Since κ̃_s depends on μ, we need a mesh in propagation direction μ. We use a piecewise constant approximation
of I on each cell (both in space and direction). This implies discontinuities at the boundaries of the cells, so
we have to change the value of ũ when going from a cell to its neighbour. Indeed, between a cell M_1 and
a cell M_2, the continuity of u implies that 

ũ_M_1(r,μ)/I_M_1(r,μ) = ũ_M_2(r,μ)/I_M_2(r,μ).

In order to take this into account, we multiply the weight of a particle going from M_1 to M_2 by
I_M_2(r,μ)/I_M_1(r,μ). 


  
A simple computation proves that, if I is an exact solution of (<ref>), then κ̃_t =
κ̃_s ≥ 0. Here, we do not have an exact solution, but one can still prove, with the same computation,
that κ̃_s ≥ 0 and κ̃_t ≥ 0. 


Proof: If I is solution to (<ref>), then I≥ 0. Hence, by definition, we have κ̃_s
≥ 0, according to (<ref>). Now, considering κ̃_t, we have, using
(<ref>) and (<ref>),

  κ̃_t = κ_t -1/I (μ∂_r I + 1-μ^2/r ∂_μI) = κ_t -
  1/I ( κ_t I - κ_s ⟨I ⟩) = κ_s⟨I ⟩/I ≥0.

Next, we consider the case in which (<ref>) is replaced by (<ref>), where S is no more equal to
κ_s ⟨ I ⟩. However, S(r)=ϕ(r)/r, where ϕ is numerically computed by solving
(<ref>), where A is defined by (<ref>) and (<ref>), and b by
(<ref>). With these definitions and the fact that the function  is negative and decreasing on
^-, we infer that I_d -A is an M-matrix, and that b_j≥ 0, for all j. Hence, S≥ 0, from which we
deduce again that I≥ 0. The proof of κ̃_s≥ 0 is exactly the same as above. Turning to κ̃_s
≥ 0, we point out that S≥ 0, hence the same computation as above gives

  κ̃_t = κ_t -1/I (μ∂_r I + 1-μ^2/r ∂_μI) = κ_t -
  1/I ( κ_t I - S) = S/I ≥0.





§ NUMERICAL RESULTS


The law of large numbers states that for independent and identically distributed random variables, the sample
average converges to the expected value when the number of random variables increases. (In the
  present context, a random variable is synonymous to a Monte Carlo particle.) In addition, the central
limit theorem implies that the rate of convergence is 1/√(N), where N is the
number of random variables. In our case, N is equal to the number of particles used in the simulation. 
For Monte Carlo methods, the computational cost is proportional to the number N of random variables. So, to
compare the performance of different methods, we have to take the calculation time T into account. Knowing that
T≈ Nt, where t is the time to generate one realization of a random variable, we define the figure of
merit (F.O.M) of a method by 

    F.O.M = 1/σ^2 t = 1/Σ_N^2T,

with σ^2 is the variance of X_i, where ( X_i)_i ∈ℕ is used in the simulation sequence of independent random variables. An unbiased estimator for σ^2 is 

    σ_N^2 = 1/N-1∑_i=1^N(X_i - X̅_N )^2   with  X̅_N = 1/N∑_i=1^NX_i.

Moreover Σ_N^2 is the variance of X̅_N. An unbiased estimator for Σ_N^2 is σ_N^2/N. In steady cases below, we apply the above formulas. In unsteady cases, we
  apply them to time-integrated values.




We are now going to give some numerical results obtained with the method developed so far. First, we provide two
verification cases, which indicate that our implementation of the Monte Carlo method is correct. In these cases, we
have analytical solutions, allowing to assess the statistical convergence of the method. Second, we provide
variance reduction tests, in which we do not have any analytical solution, so we only study variance reduction when
importance sampling is applied. The first case is stationary, and the second one is unsteady, with data in
agreement with FCI simulations. 



 §.§ A stationary verification test case



This test is borrowed from <cit.> and <cit.>, and is used as a verification procedure for
our Monte Carlo code (without using the importance sampling method). We solve Equation (<ref>) with a
point source located at r=R_ source (actually, this is a point source only if R_0=0, but we
  nevertheless use this denomination even if R_0>0.):

    Q(r) = 1/4π r^2δ(r-R_ source),

The domain is (0,1), that is, (<ref>) with R_0=0 and R_1=1. The cross sections are such that
κ_t = 1. We test κ_s = 0.3 and κ_s = 0.9. The boundary conditions correspond to zero incoming
flux:

  u(R_1,μ) = 0,   ∀μ<0. 

We use 10^5 particles, and compute the zero-moment of the intensity ψ as a function of r, on a mesh with 1000 identical
cells, that is, Δ r= 10^-3. Recall that

    ψ(r) = 2π∫_-1^1 u(r,μ)dμ,

so that an unbiased estimator of the average value of ψ on the cell M is given by 

  ψ_j = ∑_X_j∈M w_j  ⟶_N→+∞  ∫_M ψ,

where j is the index of a particle, X_j its position and w_j its weight. 

A semi-analytical solution is derived in <cit.>. This solution involves an integral which is
computed numerically. In Figure <ref> and Figure <ref>, we compare the result of
our Monte-Carlo code (without importance sampling) with this analytical solution. 



This test shows a good agreement between the result of our code and the analytical solution. The statistical noise
is more important in the cells near the origin. This can be explained by the fact that these cells are small, so
very few particles are present in them. The singularity at r=R_ source is well reproduced. 




 §.§ An unsteady verification test case



In this test, we assume that R_1 = 1, and that the target is locate at r=R_0(t), where we have set

    R_0(t) = α + β t,   α = 0.37625,   β = -0.027625

These data are borrowed from physically relevant cases of inertial confinement fusion (after
adimensionalization). The incoming flux imposed at r=R_1 is equal to 1 (with Lambert cosine law) between
  times 0 and T_ max, and 0 between T_ max and T=10. We tune T_ max so as to have an
  exact value for the flux on the target integrated between 0 and T, as is explained below.





In the case κ_s=0 and κ_t=0, we have an analytic expression for the solution, and the flux at the
boundary of the target, integrated in time from 0 to T_ max, is equal to

    F = ∫_0^T_ max∫_-1^μ(t) |μ| dμ dt,

where the time T_ max is made precise below. 

The direction  μ(t) is, for a particle generated at time t and reaching the target, the largest possible
propagation direction. In order to compute it, we compute the trajectory of the corresponding particle, defined by ṙ =
μ and μ̇= 1-μ^2/r, hence




  r(t+t')^2 = R_1^2 +2t'R_1μ(t) +t'^2,   μ(t+t') = R_1μ(t)+t'/√(R_1^2 +2t'R_1μ(t) +t'^2).

This trajectory crosses the inner boundary if and only if the equation r(t+t') = R_0(t+t') has a solution. This
equation reads R_0(t+t')^2 =  R_1^2 +2t'R_1μ(t) +t'^2, which is a second-degree equation in t'. The maximum
direction μ(t) corresponds to the case when the discriminant is 0. Computing it, we find

    μ(t) = 1/R_1[β R_0(t) - √((1-β^2 )(R_1^2 - R_0(t)^2 ))].

Hence, (<ref>) also reads

    F = 1/2∫_0^T_ max( 1-μ(t)^2) dt = 1/2∫_0^T_ max( 2
       -β^2 -R_0(t)^2/R_1^2 -
        2βR_0(t)/R_1√((1-β^2 )(1 - R_0(t)^2/R_1^2)))dt

In (<ref>), the time T_ max is such that a particle generated at t<T_ max reaches the target
with initial direction μ∈ [-1,μ(t)], where μ(t) is given by (<ref>). Indeed, for such value of T_ max<T, formula (<ref>) is valid, whereas if T_ max is
chosen to be larger, some particles, generated between T_ max and T, never reach the target. Therefore the
value of μ(t) is no more given by (<ref>), and formula (<ref>) should be modified. In order to avoid
technical difficulties associated to this new value of μ(t), we
restrict the time integral to [0,T_ max]. A simple computation shows that

T_max = T - R_0(T)β- √(R_0(T)^2β^2 - R_0(T)^2+R_1^2).

Using the above values of α and β, we find T_ max = 9.00777122797, for which we find

  F = 0.25172763696.

This computation is carried out in the following conditions: Δ r = 10^-2, Δ t = 10^-3. The
results are displayed in Figure <ref>, showing statistical convergence to the exact value
as the number of particles grows. Note that the abscissa in Figure <ref> is the number of
particles generated at each time step. Therefore, the total number of particles in the simulation is equal to
N× T/Δ t = N× 10^4. 




In the case κ_s=0 and κ_t=1, we still have an exact expression of the flux on the inner ball R_0(t):

F = ∫_0^T_max ∫_-1^μ(t) |μ|exp[-(κ_t -κ_s)(τ(t,μ) -t )] dμdt

with

τ(t,μ) = 1/1-β^2 (t+αβ- R_1μ-√((R_1μ-t-αβ)^2 - ( 1-β^2 )(t^2+R_1^2-α^2-2R_1μt ))).

This integral is not explicit, so we applied a numerical integration method to compute it. The result is

  F =  0.11385526445,

up to an estimated error of 10^-6. The parameters for this test are as follows: 100 cells in Δ r =
10^-2,  Δ t= 10^-3. 

Figure <ref> gives the results of this test, showing statistical convergence as the
number of particles grows. 





 §.§ Stationary test case: variance reduction



We now consider a stationary test case again. Here, we solve the transport equation in the domain (R_0,R_1), with
R_0 = 0.1 fixed, and R_1 = 1. An incoming flux is imposed on the outer boundary r=R_1, of value 1. The
discretization corresponds to Δ r = 10^-2 and Δμ = 2× 10^-3. Note that the discretization
in μ is only used when the importance sampling method is applied.






  §.§.§ Case κ_s=0.9, κ_t=1



  





Table <ref> and table <ref> show the results for
this test. The last column gives the proportion of particles reaching the inner sphere r=R_0. We note that very
few of them reach the inner sphere without the importance sampling method. On the contrary, an important proportion
(almost 90%) reach it when the importance sampling is applied. The figure of merit (F.O.M) is computed according to
formula (<ref>). When applying the importance sampling method, we provide two execution times, and therefore
two values for the F.O.M. The first one includes the computation of the importance function, which is not
meaningful from a statistical viewpoint, although it is from a computational cost viewpoint. On the contrary, the second value (Time2, and F.O.M2), do not include it,
and therefore give a clear meaning to the statistical efficiency of the method. As expected, when the number of
particles grows, these two F.O.M are very close to each other. In the case of a small number of particles, the
importance sampling method is less efficient (the F.O.M is increased only by a factor 20) because the calculation of the importance function is too expensive compared to the Monte
Carlo method. When a large number of particles is used, however, the method is much more efficient, and we see that
the F.O.M is more than 100 times better than without importance sampling. Figure
<ref> shows an important variance reduction, for any number of particles. 



  §.§.§ Case κ_s=0.1, κ_t=1






Table <ref> and table <ref> show the results for
this test. Here again, the last column gives the proportion of particles reaching the inner sphere r=R_0. We note that very
few of them reach the inner sphere without the importance sampling method. On the contrary, an important proportion
(about 90%) reach it when the importance sampling is applied. The figure of merit (F.O.M) is computed according to
formula (<ref>). 




Here again, we provide two execution times and two F.O.M, the first one including the computation
  of the importance function, the second one excluding it.
In the case of a small number of particles, the
importance sampling method is less efficient (the F.O.M is increased by a factor 50) because the calculation of the importance function is too expensive compared to the Monte
Carlo method. When a large number of particles is used, however, the method is much more efficient, and we see that
the F.O.M is more than 150 times better than without importance sampling. Figure
<ref> shows an important variance reduction, for any number of particles. 











 §.§ Unsteady test case



We use now a test case in which the inner sphere has a radius which depends on time, according to the same law as
in Subsection <ref>, that is,

  R_0(t)=0.37625-0.027625t.

As mentioned above, this value of R_0(t) has been derived from three-dimensional simulations of ICF. The incoming
flux imposed on the outer sphere R_1 = 1 is constant in time, so we impose a number of particles generated at
each time step. We fix the values of κ_s=0.9 and κ_t=1, although other values give the same kind of
results. Here again, we use a uniform mesh in r
with Δ r= 10^-2, a time step Δ t = 10^-3 and a direction discretization Δμ = 2×
10^-3. The output value of the code is the flux at the moving boundary R_0(t), integrated from time t=0 to
time t=10, which is the final time of the simulation. 

The results are presented in Table <ref> et Table
<ref>, in which the first column N is the number of particles generated at each
time step. Therefore, the total number of particles generated in the simulation is N× T/Δ t = N×
10^4. As above, the computation time Time1 includes the computation of the importance function, which is done at
each time step, whereas the computation time Time2 does not. The last column in each table represents the
proportion of particles reaching the target. 



[1]N is the number of particles generated at each time step.



















§ CONCLUSION


We have presented in this paper a new method of variance reduction based on importance sampling for the transport
equation in spherical geometry. The importance function is computed as the solution of the adjoint equation, which
is solved numerically. In order to do so, we use an integral equation derived by Siewert and Thomas
<cit.>, and solve this equation numerically in order to find the first moment (with respect to the
direction μ of
the importance function). Once this is computed, we apply the method of characteristics to compute the importance
function. Contrary to what has been done in <cit.> in a similar context, we do not have an analytical
expression for the importance function. However, it should be noted that the importance function used in
<cit.> does not satisfy the correct boundary conditions. Therefore, it is adapted only if boundary
conditions are not of too much importance in the computation at hand. We expect that this is case in ICF
  experiments, in which photons can move a long distance before absorption or scattering. In our method, we are not
  limited by such considerations. 

Numerical tests indicate that the method is efficient, including situations close to the case of inertial
confinement fusion. 






In future works, we plan to test this method in situations closer to experiments. In this respect,
  several issues need to be considered:
  
  
  * the absorption and scattering coefficients are not constant, contrary to the assumptions we have made
    here. These heterogeneities make the calculation of I (Section <ref>) much more
    difficult. For instance, the generalization of (<ref>) might lead to complicated expressions, thereby
    impeding the numerical efficiency of the method. One way (among others) to circumvent this difficulty would
    then be to use a numerical approximation for the computation of I itself. 
  
  * The problem is by nature frequency dependent. Although using a grey importance function in such a
    simulation is possible, one should bear in mind that it might prove insufficient. Hence, including a dependence
    of I upon the frequency will probably be an important question to be dealt with. 
  
  * The problem we studied here does not take into account interaction with matter. Although we have used
    transient simulations in which this effect is partly represented by the movement of the detector, the
    hydrodynamics of the plasma imply a much richer interaction. This will imply new issues to be considered. 
  
  * In relation with the preceding point, a mesh used for an ICF simulation is in general highly
    heterogeneous. This is an additional problem to be considered. 
  






§ INTEGRAL EQUATION ON Φ



We give in this Appendix the details of the derivation of the integral equation satisfied by ϕ. We define

  μ_d(r) = -√(1-R_0^2/r^2),

and

  ϕ(r) = r S(r) = κ_s r ⟨I⟩.

Equation (<ref>) is equivalent to 


    I(r,μ)    = exp( κ_t( rμ + √( r^2μ^2 - r^2 + R_0^2)))1_{μ < μ_d(r)} + ∫_rμ^R(r,μ) S(√(s^2+r^2-r^2μ^2))exp(κ_t (rμ-s) )ds 
        = exp( κ_t( rμ + √( r^2μ^2 - r^2 + R_0^2)))1_{μ < μ_d(r)}
        + 1_{μ < μ_d(r)}∫_rμ^-√(r^2μ^2-r^2+R_0^2) S( √(s^2 + r^2-r^2μ^2))exp(κ_t ( rμ-s))  ds 
        + 1_{μ_d(r)<μ < 0}∫_rμ^√(r^2μ^2-r^2+R_1^2) S( √(s^2 + r^2-r^2μ^2))exp(κ_t( rμ-s)) ds 
        + 1_{μ > 0}∫_rμ^√(r^2μ^2-r^2+R_1^2) S( √(s^2 + r^2-r^2μ^2))exp(κ_t( rμ-s))  ds,

That is,

    I(r,μ)    = exp( κ_t( rμ + √( r^2μ^2 - r^2 + R_0^2)))1_{μ < μ_d(r)}
        + 1_{μ < μ_d(r)}∫_rμ^-√(r^2μ^2-r^2+R_0^2) S( √(s^2 + r^2-r^2μ^2))exp(κ_t( rμ-s ))  ds 
        + 1_{μ_d(r)<μ < 0}( ∫_rμ^0 S( √(s^2 + r^2-r^2μ^2))exp(κ_t( rμ-s))ds
    . 
       .+ ∫_0^√(r^2μ^2-r^2+R_1^2) S( √(s^2 + r^2-r^2μ^2))exp(κ_t( rμ-s )) ds)    
        + 1_{μ > 0}∫_rμ^√(r^2μ^2-r^2+R_1^2) S( √(s^2 + r^2-r^2μ^2))exp(κ_t ( rμ-s ))  ds.

Hence, changing variables r' = √(s^2+r^2-r^2μ^2), we find

    I(r,μ)    = exp( κ_t( rμ + √( r^2μ^2 - r^2 + R_0^2)))1_{μ < μ_d(r)}
        + 1_{μ < μ_d(r)}∫_R_0^r r'S( r')exp( κ_t( rμ+√(r^2μ^2-r^2+r'^2)) ) dr'/√(r^2μ^2-r^2+r'^2)
        + 1_{μ_d(r)<μ < 0}( ∫_r√(1-μ^2)^r r'S( r')exp( κ_t( rμ+√(r^2μ^2-r^2+r'^2)) ) dr'/√(r^2μ^2-r^2+r'^2).
        + .∫_r√(1-μ^2)^R_1 r'S( r')exp( κ_t( rμ-√(r^2μ^2-r^2+r'^2)) ) dr'/√(r^2μ^2-r^2+r'^2)) 
        + 1_{μ > 0}∫_r^R_1 r'S( r')exp( κ_t( rμ-√(r^2μ^2-r^2+r'^2)) ) dr'/√(r^2μ^2-r^2+r'^2).

Next, we integrate with respect to μ:

    ∫_-1^1I(r,μ)dμ   = ∫_-1^1exp( κ_t( rμ + √( r^2μ^2 - r^2 + R_0^2)))1_{μ < μ_d(r)} dμ
        + ∫_-1^11_{μ < μ_d(r)}∫_R_0^r r'S( r')exp( κ_t( rμ+√(r^2μ^2-r^2+r'^2)) ) dr'/√(r^2μ^2-r^2+r'^2) dμ
        + ∫_-1^11_{μ_d(r)<μ < 0}( ∫_r√(1-μ^2)^r r'S( r')exp( κ_t( rμ+√(r^2μ^2-r^2+r'^2)) ) dr'/√(r^2μ^2-r^2+r'^2).
        + .∫_r√(1-μ^2)^R_1 r'S( r')exp( κ_t( rμ-√(r^2μ^2-r^2+r'^2)) ) dr'/√(r^2μ^2-r^2+r'^2)) dμ
        + ∫_-1^11_{μ > 0}∫_r^R_1 r'S( r')exp( κ_t( rμ-√(r^2μ^2-r^2+r'^2)) ) dr'/√(r^2μ^2-r^2+r'^2)dμ.

We deal with each term seperately.


  -Let 
    I_0(r) = ∫_-1^μ_d(r)exp( κ_t( rμ + √( r^2μ^2 - r^2 + R_0^2))) dμ.


Setting

    θ = rμ + √( r^2μ^2 - r^2 + R_0^2),
 we obtain

    I_0(r) = 1/2r( [1/κ_texp(κ_t θ) ]_-r+R_0^-√(r^2-R_0^2) + ( R_0^2-r^2 ) [ κ_tEi(κ_tθ) - exp( κ_tθ)/θ]_-r+R_0^-√(r^2-R_0^2)).



  -Let

    I_1(r)    = ∫_-1^11_{μ < μ_d(r)}∫_R_0^r r'S(r')exp( κ_t( rμ+√(r^2μ^2-r^2+r'^2)) ) dr'/√(r^2μ^2-r^2+r'^2) dμ
       =  ∫_R_0^r r'S(r') ( ∫_-1^μ_d(r)exp( κ_t( rμ + √( r^2μ^2 - r^2 + r'^2))) dμ/√(r'^2-r^2 +r^2μ^2)) dr'


We change variables as follows:

    ν = rμ + √( r^2μ^2 - r^2 + r'^2),
 and get

    I_1(r)     = ∫_R_0^r r'S(r') ( ∫_-r+r'^-√(r^2-R_0^2)+√(r'^2-R_0^2)exp( κ_t ν) dν/rν)  dr' 
        = 1/r∫_R_0^r r'S(r') [ Ei( κ_t ν) ]_-r+r'^-√(r^2-R_0^2)+√(r'^2-R_0^2) dr',

whence 

    I_1(r) = 1/r∫_R_0^r r'S(r') [ Ei( κ_t ( -√(r^2-R_0^2)+√(r'^2-R_0^2))  ) - Ei( κ_t (-r+r' )) ] dr'.



  -Let

    I_2(r)    = ∫_-1^11_{μ_d(r)<μ < 0}∫_r√(1-μ^2)^r r'S( r')exp( κ_t( rμ+√(r^2μ^2-r^2+r'^2)) ) dr'/√(r^2μ^2-r^2+r'^2)dμ
       =  ∫_R_0^r r'S(r') ( ∫_μ_d(r)^-√(1-r'^2/r^2)exp( κ_t( rμ + √( r^2μ^2 - r^2 + r'^2))) dμ/√(r'^2-r^2 +r^2μ^2)) dr'


Changing variables according to 

    ν = rμ + √( r^2μ^2 - r^2 + r'^2),
 we find

    I_2(r)     = ∫_R_0^r r'S(r') ( ∫^-√(r^2-r'^2)_-√(r^2-R_0^2) +√(r'^2-R_0^2)exp( κ_t ν) dν/rν)  dr' 
        = 1/r∫_R_0^r r'S(r') [ Ei( κ_t ν) ]^-√(r^2-r'^2)_-√(r^2-R_0^2)+√(r'^2-R_0^2) dr',

thus

    I_2(r) = 1/r∫_R_0^r r'S(r') [Ei( κ_t (-√(r^2-r'^2))) -  Ei( κ_t ( -√(r^2-R_0^2)+√(r'^2-R_0^2))  ) ] dr'.



  -Let

    I_3(r)    = ∫_-1^11_{μ_d(r)<μ < 0}∫_r√(1-μ^2)^R_1 r'S( r')exp( κ_t( rμ-√(r^2μ^2-r^2+r'^2)) ) dr'/√(r^2μ^2-r^2+r'^2)
       =  ∫_R_0^R_1 r'S(r') ( ∫_μ_d(r)^-√(1-r'^2/r^2)1_{ r'<r }exp( κ_t( rμ - √( r^2μ^2 - r^2 + r'^2))) dμ/√(r'^2-r^2 +r^2μ^2)) dr'

 
Changing variables by setting

    ν = rμ - √( r^2μ^2 - r^2 + r'^2),
 we infer

    I_3(r)     = ∫_R_0^R_1 r'S(r') ( -∫^-√(|r^2-r'^2|)_-√(r^2-R_0^2) - √(r'^2-R_0^2)exp( κ_t ν) dν/rν)  dr' 
        = 1/r∫_R_0^R_1 r'S(r') [ Ei( κ_t ν) ]_-√(|r^2-r'^2|)^-√(r^2-R_0^2)-√(r'^2-R_0^2) dr'.

Hence,

    I_3(r) = 1/r∫_R_0^R_1 r'S(r') [Ei( κ_t ( -√(r^2-R_0^2)-√(r'^2-R_0^2))  ) - Ei( κ_t (-√(|r^2+r'^2|)))  ] dr'.



  -Let

    I_4(r)    = ∫_-1^11_{μ > 0}∫_r^R_1 r'S(r') exp( κ_t( rμ - √( r^2μ^2 - r^2 + r'^2))) dr'/√(r'^2-r^2 +r^2μ^2)dμ
       =  ∫_r^R_1 r'S(r') ( ∫_0^1exp( κ_t( rμ - √( r^2μ^2 - r^2 + r'^2))) dμ/√(r'^2-r^2 +r^2μ^2)) dr'


We change variables as follows:

    ν = rμ - √( r^2μ^2 - r^2 + r'^2),
 and get

    I_4(r)     = ∫_r^R_1 r'S(r') ( -∫^r-r'_-√(r'^2-r^2)exp( κ_t ν) dν/rν)  dr' 
        = 1/r∫_r^R_1 r'S(r') [ Ei( κ_t ν) ]_r-r'^-√(r'^2-r^2) dr',

from which we deduce

    I_4(r) = 1/r∫_r^R_1 r'S(r') [ Ei( κ_t ( -√(r'^2-r^2))  ) - Ei( κ_t (r-r' )) ] dr'.

Collecting all the above results, we have

    ∫_-1^1 I(r,μ) dμ    = I_0(r) + I_1(r) + I_2(r) + I_3(r) + I_4(r) 
        = I_0(r) + 1/r∫_R_0^R_1 r'S(r') [ Ei( κ_t (-√(r^2-R_0^2) -√(r'^2-R_0^2))  ) - Ei( κ_t (-|r-r'| )) ] dr'.

Recalling the definition of ϕ, we get

    ϕ(r)     = κ_s r/2( I_0(r) + 1/r∫_R_0^R_1ϕ(r') [ Ei( κ_t (-√(r^2-R_0^2) -√(r'^2-R_0^2))  ) - Ei( κ_t (-| r-r'| )) ] dr' ),

hence

    ϕ(r) = κ_s/4   ( 1/κ_t[ exp(κ_t θ) ]_-r+R_0^-√(r^2-R_0^2) + (R_0^2 - r^2 )[ κ_t ( κ_t θ) - exp( κ_t θ)/θ]_-r+R_0^-√(r^2-R_0^2)) 
        + κ_s/2∫_R_0^R_1ϕ(r') [ ( κ_t (-√(r^2-R_0^2) - √(r'^2-R_0^2)) ) - ( κ_t(-| r -r' | ))   ] dr'.




plain



