## Abstract

The characteristics of the transient and polarization must be considered for a complete and correct description of short-pulse laser transfer in a scattering medium. A Monte Carlo (MC) method combined with a time shift and superposition principle is developed to simulate transient vector (polarized) radiative transfer in a scattering medium. The transient vector radiative transfer matrix (TVRTM) is defined to describe the transient polarization behavior of short-pulse laser propagating in the scattering medium. According to the definition of reflectivity, a new criterion of reflection at Fresnel surface is presented. In order to improve the computational efficiency and accuracy, a time shift and superposition principle is applied to the MC model for transient vector radiative transfer. The results for transient scalar radiative transfer and steady-state vector radiative transfer are compared with those in published literatures, respectively, and an excellent agreement between them is observed, which validates the correctness of the present model. Finally, transient radiative transfer is simulated considering the polarization effect of short-pulse laser in a scattering medium, and the distributions of Stokes vector in angular and temporal space are presented.

© 2013 Optical Society of America

## 1. Introduction

Short-pulse laser has widely practical applications in many fields. In particular, the femtosecond lasers are now commonly used in the optical measurements and experiments [1–3]. It also finds applications in atmospheric science to estimate radiation budget, in biomedical science to diagnose tumors and in performing laser surgery. The interaction of short-pulse laser and medium has become the focus of current research. Although in some cases the polarization can be neglected, both the transient and polarization effects must be considered in a complete and correct description for short-pulse laser transporting in a scattering medium. Especially for short-pulse laser transfer in a strong scattering medium, significant errors will be introduced into transient results if the polarization is ignored. Unfortunately, may be because of the complexity of polarization problem, most of the studies on short-pulse laser transmitting in scattering media only consider transient radiative transfer and ignore the polarized effect. To solve the transient scalar radiative transfer, many numerical methods have been developed, such as the discrete ordinates (DO) method [4–8], the Monte Carlo (MC) method [9–11], the discrete transfer (DT) method [4, 12, 13], the finite volume (FV) method [4, 14–16] and the quadrature method [17, 18].

While since the importance of polarization information, the investigation of vector radiative transfer has drawn great attention in many research fields, but most of works focus on the steady state problem. So many numerical methods for steady-state vector radiative transfer have been proposed and developed, including MC method [19], *F _{N}* method [20], adding-doubling method [21], discrete-ordinate (DO) method [22], successive-order-scattering (SOS) method [23], Markov-chain method [24] and so on. Of the available methods, MC method [25–32] is so popular especially in recent years with the rapid development of computer technology.

Owing to its complexity, only a few works have been reported that take into account both the polarization and the transient aspects of the radiative transfer problem. Ishimaru et al. [1] solved the transient vector radiative transfer equation for a plane-parallel Mie scattering medium by using the DO method, and calculated the time-dependent degree of polarization and cross-polarization discrimination. Wang et al. [2] examined polarized light transmitting through randomly scattering media with a polystyrene-microsphere solution, and using MC technique calculated temporal profiles for the Stokes vectors and the degree of polarization. While in [2], X. D. Wang et al. only presented time-resolved results only for one receiving angle and did not give the angular distribution of the results. Sakami et al. [33] used the DO method to investigate the propagation of a polarized pulse in random media, and only the time- and angular-resolved degree of polarization was calculated. Ilyushin et al. [34] applied the small angle approximation of the radiative transfer theory to study the problem of the propagation of the short light pulse in the scattering medium considering polarized effect. While the results from [34] only included the time-resolved intensity and degree of polarization. Up to now, we have only found the above several papers reporting transient vector radiative transfer problem, while they did not show the complete distributions of Stokes vector in angular and temporal space. It’s necessary to carry on a further research on the problem for revealing the physical behavior and characteristics of transient vector radiative transfer in scattering media.

Our purpose in this paper is to develop a Monte Carlo method for solving the radiative transfer problem that considers both the polarization and the transient effect. Compared with other numerical methods, the MC method has two major advantages in solving the radiative transfer problem including the transient and polarization effects. In the MC method, the conversion between space and time is easy to achieve, and for this merit, the MC method is very suitable for solving transient radiative transfer. Besides, as a method based on ray tracing, the MC scheme has a unique advantage for solving radiative transfer problem with polarization caused by Fresnel surface. While a major disadvantage of the MC method is that it is relatively time consuming for obtaining satisfactory results, especially for the transient radiative transfer problem. If the number of beams is not big enough, the rapid fluctuations in the results by the MC method will appear for the case of large optical thickness [9–11]. In this paper, we introduce and develop the basic theory for MC simulation of transient vector radiative transfer, including stochastic process of photons, conversion between transferring distance and time, the determination of scattering directions using the modified reject method, and normalization of Stokes vector. In order to describe the transfer characteristics of short pulse completely, the transient vector radiative transfer matrix is defined. A time shift and superposition principle is first applied to the MC model to improve the computational efficiency and accuracy. The correctness of the models for transient scalar radiative transfer and steady-state vector radiative transfer are validated, respectively. Finally, we simulate the radiative transfer with the consideration of the polarization and transient effects of short-pulse laser in a Mie scattering medium layer, and illustrate the distributions of Stokes vector in angular and temporal space.

## 2. Theory

#### 2.1 TVRTE and Stokes vector

If considering a finite propagation velocity of photons, the transient vector radiative transfer equation (TVRTE) is used to describe the transfer behavior for polarized monochromatic radiation in participating media with randomly oriented particles, and can be written as:

**is the Stokes vector;**

*S***is the discrete direction vector;**

*Ω***r**is the spatial coordinate vector;

*κ**is the extinction matrix;*

_{e}

*κ**is the absorption vector;*

_{a}*I*is the black body intensity;

_{b}*c*is the velocity of light in the medium;

*t*is the time;

**is the scattering phase matrix.**

*Z*In Eq. (1), the Stokes vector ** S** is expressed by

*I*is the intensity,

*Q*is the linear polarization aligned parallel or perpendicular to the

*z*-axis,

*U*is the linear polarization aligned ± 45° to the

*z*-axis and

*V*is the circular polarization. In this paper, the MC technique is used to solve the transient polarized radiative transfer in a participating medium. The term of

*κ**describing the thermal radiation emitted by the medium in Eq. (1) is ignored for the cold medium considered in the present work.*

_{a}I_{b}#### 2.2 Propagation

In the MC simulation of radiative transfer in a participating medium, the photons which experience free propagation, scattering, reflection, transmission and other events are traced in the medium until they escape from the boundaries or are absorbed. The deterministic behavior of light is determined by the random behavior of a finite number of photons in the medium, and each photon randomly propagates along the possible path.

To analyze the transient effect of radiative transfer, we consider a square continuous pulse laser chain with an external collimated irradiation at the bottom surface of a slab. Just as shown in Fig. 1, the pulse width is *t _{p}*, the number of pulses is

*p*, the time interval between two successive pulses is

*t*, and the intensity of the pulse is

_{d}*I*

_{0}. Each pulse emits

*N*photons within the pulse width of

*t*, and every photon starts its trajectory at the time of

_{p}*R*is the uniform random number.

_{t}The distance a photon travels without being absorbed, scattered, reflected or refracted is determined by Eq. (4a):

where*L*is the free distance,

*κ*is the extinction coefficient and

*R*is a pseudo random number with an uniform distribution. After travelling the free distance, the new coordinate position$\left(x\text{'},y\text{'},z\text{'}\right)$of the photon is obtained by:

_{L}*x*,

*y*,

*z*) is the initial position of the photon,

*θ*and

*φ*are the polar and azimuth angles of the propagation direction, respectively. The time taken for the distance of

*L*is

*t*, given by Eq. (5a):where

_{L}*n*is the refractive index of the medium, and

*c*

_{0}and

*c*are the speed of light in vacuum and medium, respectively. Then the new time is determined by Eq. (5b):In this paper, the time

*t*is normalized to

*t** as Eq. (6):where

*κ*

_{0}is the reference extinction coefficient just equal to 1 m

^{−1}.

#### 2.3 Scattering

The scattering problem for polarized radiative transfer is much more complicated than that for scalar one. In the transport process, the photons will probably interact with the particles in the medium. The photons may be absorbed by the medium or scattered to other directions by the particles. A pseudo random number is compared with the scattering albedo *ω* of the medium to judge whether the photon is absorbed or scattered. If the pseudo random number *R _{ω}* >

*ω*the photon will be absorbed by the medium and stop tracking the photon, otherwise the photon will be scattered by the particle.

Chandrasekhar [35] first proposed the idea of the photon scattering from one location to the next in the meridian and scattering planes. Figure 2 shows the schematic depiction of the scattering problem. The photon travelling directions, before and after scattering, are represented as those pointed to by the lines *OA* and *OB* respectively on the unit sphere. The incident Stokes vector ** S_{i}** and the resulting scattering Stokes vector

**determine a plane**

*S*_{s}*AOB*, defined as the scattering plane. The incident Stokes vector

**and the**

*S*_{i}*z*-axis determine a plane

*AOC*, called as the meridian plane. After a photon is scattered by the particle, a new meridian plane is created by the new propagation direction of

**(scattering Stokes vector) and the**

*S*_{s}*z*-axis, and a transformation of the Stokes vector must be done. The polarization caused by the scattering can be properly interpreted with respect to the new meridian plane

*BOC*, as shown in Fig. 2. For the Stokes vector, two rotations into and out of the scattering plane are required in the MC method. Based on the above analysis, the resulting Stokes vector after scattering is:

**is the single scattering Mueller matrix defined in Van de Hulst [36];**

*M***is a rotation matrix;**

*L**i*

_{1}and

*i*

_{2}are the angles of rotation into and out of the scattering plane, respectively. The phase matrix

**in Eq. (1) is equal to**

*Z***(**

*L**π*−

*i*

_{2})

**(**

*M**Θ*)

**(−**

*L**i*

_{1}).

The matrix ** M** contains 4 × 4 elements for describing the transformation of the Stokes vector of the incident wave into that of the scattered wave. For a spherical-particle scattering the matrix

**can be simplified as one having eight non-zero elements [37], of which only four are independent. The form of**

*M***is given as below:**

*M**Θ*is the scattering angle between the incident and scattering directions, as shown in Fig. 2.

The matrix ** L** represents a rotation in the clockwise direction with respect to an observer looking into the direction of the photon propagation [38] and can be written as:

The physical meaning of Eq. (7), expression of Stokes vector after scattering, can be well understood. First, the incident Stokes vector ** S_{i}** is rotated by the angle of

*i*

_{1}counterclockwise from the meridian plane

*AOC*to the scattering plane

*AOB*. Then, the incident Stokes vector

**interacts with the particle and is scattered into a new direction. Finally, rotate the scattering Stokes vector by the angle of**

*S*_{i}*π*−

*i*

_{2}clockwise from the scattering plane

*AOB*to the meridian plane

*BOC*and obtain the scattering Stokes vector

**.**

*S*_{s}For the MC simulation of light scattering problem considering the polarization, the determination of scattering angles is one of the major difficulties. For a given incident direction (*θ*_{1},*φ*_{1}) as shown in Fig. 2, if the angles *i*_{1} and *Θ* are obtained, the angles *i*_{2}, *θ*_{2} and *φ*_{1}−*φ*_{2} can be calculated from the spherical laws of sines and cosines [39]. For gaining these angles, the cumulative distribution method [32] was used. In this paper we use the rejection method to get *i*_{1} and *Θ*, and the detailed information about the method can be seen in [40].

In the Stokes vector, the first component denotes the energy carried by the photon, and the other three components store the states of polarization. So in the MC simulation, in order to ensure the conservation of energy, the normalization must be applied to the Stokes vector in the propagation of the photon. That is to say the *I*-Stokes element of the photon must always keep one as it propagates through. The normalized Stokes vector is written as:

#### 2.4 Reflection and transmission

At the medium boundary or the interface between the two adjacent layers, the reflection or transmission of the incident photon will take place. According to the reflection characteristics, the medium interfaces or boundary mainly include Fresnel surface and Lambertian surface. After reflection or transmission, the resulting Stokes vector is obtained by Eq. (11a) or (11b):

where**and**

*R***are reflection matrix and transmission matrix, respectively. Note that**

*T*

*S**and*

_{ref}

*S**also need to be normalized as done in Eq. (10).*

_{tra}For the semitransparent Lambertian surface, the reflection matrix and transmission matrix are given as follows:

*ρ*is the reflection coefficient of the Lambertian boundary which is independent to the incident light and just related to the properties of the surface.

We compare a uniform random number *R _{ref}* with the

*ρ*to judge whether the reflection takes place or not at the boundary. If

*R*<

_{ref}*ρ*the photon is reflected back to the medium on the incident side, otherwise it is refracted into the adjacent medium. The reflection direction is determined by the following equations:

*R*and

_{θ}*R*are uniform pseudo random numbers. Both the incident angle and reflection angle are measured with respect to the normal direction of the reflection plane.

_{φ}For the Fresnel surface, the reflection and transmission matrices have been presented in a few papers [41–43]. The reflection matrix is written as Eq. (14a):

*n*and

_{i}*n*are the refractive indexes of the media for the incident side and the transmission side, respectively. The transmission matrix is expressed as Eq. (14c):

_{t}The reflection and transmission behaviors are dominated by the laws of reflection and refraction. The angles after reflection and refraction are decided by the following equations:

*r*,

*i*and

*t*denote the reflection, incidence and transmission (refraction), respectively.

We also need compare the reflection coefficient *ρ* with a uniform pseudo random number *R _{ref}* to judge whether the reflection or transmission takes place when the photon hits the Fresnel boundary. For the scalar radiative transfer, the reflection coefficient is written as Eq. (16):

*I*is the intensity for incident wave. We can get $I\text{'}$from the Eq. (11a) and Eq. (14a), and then $\rho \text{'}$is obtained, expressed by Eq. (17b):

*Q*≡ 0, Eq. (17b) will degenerate to Eq. (16) just for scalar radiative transfer.

#### 2.5 Photon’s statistics and post-processing

For one-dimensional problem, in order to obtain the time-resolved Stokes vector distribution in angular space at a cross section of the medium, we need to record numbers of photons passing through the cross section in each direction within each time interval. In this paper, we present the concept of transient vector radiative transfer matrix (TVRTM) used to solve the Stokes vector. For one-dimensional problem, the TVRTM denoted by *D*_{i}_{,}_{t}_{,}_{θ}_{,}* _{φ}* is defined as follows:

*N*indicates the number of photons emitted within the pulse width of

*t*;

_{p}*N*denotes the total number of the photons passing through the cross-section

_{t}*i*along (

*θ, φ*) direction in the solid angle

*d*Ω from

*t*to

*t*+ Δ

*t*moment;

*S**is the normalized Stokes vector carried by photon*

_{q}*q*. The physical meaning of TVRTM is clear and it represents the transfer of intensity and polarized states of photons from the emission position to receiving position at the

*t*moment. Finally, we obtain the transient Stokes vector in terms of TVRTM, written by:

*I*

_{0}is the intensity of incident pulse laser, and cos

*θ*

_{0}is the incident-direction cosine.

#### 2.6 Time shift and superposition

As we know, the MC simulation of transient radiative transfer has a high computational cost. If considering the polarization in the transient problem, a very huge number of photons are required to gain reliable results, which results in a very low computational efficiency. In fact, in the current computational environment, in order to obtain the results with good accuracy, we can hardly afford the computational cost for running the MC simulation of transient vector radiative transfer. To improve the computational efficiency, we introduce the time shift and superposition (TSS) principle in the MC model. For one thing, the behaviors of the photons emitted at the same position in different time transferring in the medium are the same, so the time shift principle is suitable to the transient radiative transfer problem. For another thing, radiative transfer in the medium subjected to an external short-pulse laser irradiation is a typical linear problem, and for any linear system, the principle of superposition is applicable.

We divide the pulse width *t _{p}* into

*N*

_{Δ}

*time steps, and each time step Δ*

_{t}*t*=

*t*/

_{p}*N*

_{Δ}

*. We only need to solve the Stokes vector response caused by the photons emitted at the time zero within the time step of Δ*

_{t}*t*according to the Eq. (18) and Eq. (19) letting

*t*= Δ

_{p}*t*, denoted by

**, and for the photons emitted at any other moment within the same time step, from the time shift principle the corresponding Stokes vector response**

*S*^{’}**can be obtained on the base of the Stokes vector response for the photons emitted at the time zero. Finally, combined with the principle of superposition, the transient Stokes vector**

*S*^{’}

*S**can be calculated by Eq. (20):*

_{i,t,θ,φ}## 3. Results and discussion

In the section 2, we developed a transient vector MC model for polarized radiative transfer. So far, we have found only four papers [1, 2, 33, 34] addressed the related issues. Most works about transient radiative transfer are only for scalar one. In this section, we need to start with the validation of the correctness of the model. However, as some calculating parameters were not explicitly given, the results from [1, 2, 33, 34] cannot be used to compare with our results. Therefore, we can only compare the results for transient scalar radiative transfer and steady-state polarized radiative transfer against those available in the literatures, respectively. On the basis of the validation of the model, we have a deep investigation into the transient vector radiative transfer in the scattering slab. Note that the computing environment we used is Intel(R) Core i5-2320 CPU @ 3.00Hz.

#### 3.1 Transient radiative transfer in an isotropic scattering slab

In this case we present the results comparing against those obtained by FVM in [16] to verify the MC model for transient radiative transfer. A short-pulse laser source with pulse width of *t _{p}** = 1.0 is loading on the

*A*

_{1}

^{−}surface just as shown in Fig. 3. Photons emit continuously at an incident angle of

*θ*

_{0}= 0° for a pulse width from

*t** = 0 moment. The pulse is unpolarized with an irradiance of

*S*_{0}= (1, 0, 0, 0)

^{T}. The medium with refractive index

*n*= 1.0 is isotropic scattering and non-absorbing (

*ω*= 1.0). Both

*A*

_{1}and

*A*

_{2}are assumed Lambertian boundaries. According to [16] the dimensionless reflectance and transmittance are defined as Eq. (21)

Figure 4 plots the curves of the dimensionless reflectance and transmittance for different optical thickness. As shown in Fig. 4, an excellent agreement between the results obtained by the MC model of this paper and those by FVM [16] is achieved. Benefitting from the TSS principle, only the number of photons of 1.0 × 10^{6} is required for the comparison in the Fig. 4. If the TSS principle is not adopted in the MC model, a much bigger number of photons are needed for getting convergent results. As shown in Fig. 5, we can see that a rapid fluctuation still exists in the result obtained by MC method even with the number of photons of 1.0 × 10^{8}, while once the TSS principle is considered in the MC method, only 1.0 × 10^{6} photons are necessary for gaining the convergent results with good accuracy. In fact, to get smooth transient results for optical thickness of 0.1, more than 1.0 × 10^{8} photons are necessary in MC simulation without consideration of TSS. Thus we can conclude that using the TSS the number of photons required in MC method can be reduced by more than two orders of magnitude. In addition, Table 1 shows the comparison of computation time of two schemes, and we can see that the advantage of MC_TSS scheme in computational time is very obvious. If we consider the polarization in transient radiative transfer, the improvement in computational efficiency using MC_TSS will be more remarkable.

#### 3.2 Steady-state vector radiative transfer in an atmosphere-water system with Rayleigh scattering

In this case we present the results comparing against those obtained by the coupled atmosphere-water vector discrete ordinate radiative transfer (CAW-VDISORT) code in [44] to verify the MC steady-state polarized radiative transfer model developed in this paper. The simulation is performed for polarized radiative transfer in a plane–parallel medium composed of Rayleigh scattering atmosphere layer and water body, as shown in Fig. 6. The atmosphere-water interface is supposed to be Fresnel surface. A relative refractive index of water to the atmosphere is taken as 1.338. The optical depth is taken to be 0.15 for the atmosphere and 1.0 for the water. Both the atmosphere layer and water body are assumed to be non-absorbing (*ω* = 1.0), and the bottom of the water is assumed to be completely absorbing. The down radiation field at the top of the atmosphere is taken to be a collimated and unpolarized beam of irradiance *S*_{0} = (1,0,0,0)^{T} with an incident direction of (*θ*_{0}, *φ*_{0}) = (60°,0°) as shown in Fig. 6. The number of the photons is 3.0 × 10^{9} and the angular discretization is taken as *N _{θ}* ×

*N*= 60 × 20 for the comparison. To satisfy the steady-state working condition, we set

_{φ}*t** = ∞ in the transient vector MC model.

_{p}Figure 7 shows the curves for the four parameters *I*, *Q*, *U* and *V* of Stokes vector just above the water surface at *t** = 25 when the steady state has been reached. A negative cosine of a zenith angle indicates downward radiation, while a positive cosine indicates upward radiation. The results of Stokes vector only for azimuth angle of *φ* = 90° are presented. In the simulation, the reflection coefficients *ρ* from Eq. (20) and $\rho \text{'}$ from Eq. (21) are compared with the pseudo random number to judge whether the reflection takes place or not, respectively. From the Fig. 7 we can see that both the results by *ρ* and $\rho \text{'}$ agree well with those in [44] by DOM, and compared with the results by *ρ*, the results by $\rho \text{'}$ are a little better. Thus in the following simulation, only the results by $\rho \text{'}$ are presented. Figure 8 demonstrates the distributions of Stokes vector above the water surface varying with the time. We can see the values of Stokes vector change greatly before *t** = 10, and after that tend to stabilize. After *t** = 20 the polarized radiative transfer completely reach the steady state, and the distributions of Stokes vector don’t change again.

#### 3.3 Transient vector radiative transfer in Mie scattering medium with a short-pulse irradiation

In this case we investigate transient vector radiative transfer in a Mie scattering medium layer subjected to a short-pulse laser irradiation. As seen in Fig. 3, at the moment of *t** = 0, an oblique parallel beam of laser pulse (*μ*_{0} = cos*θ*_{0} = 0.5, *φ*_{0} = 0) penetrates into a scattering layer from the *A*_{1}^{−} surface, and the pulse width is *t _{p}** = 1. The pulse is a polarized beam of irradiance

*S*_{0}= (1,0, 0.5, 0.5,0.5)

^{T}. The optical thickness

*τ*with respect to the layer thickness

_{L}*L*is 1, and the single scattering albedo

*ω*is 1. The phase matrix is for Mie scattering at a wavelength of 0.951 μm from a gamma distribution of particles with effective radius of 0.2 μm and effective variance of 0.07. The Legendre coefficients for the four unique elements of the phase matrix have been obtained by Evens and Stephens [21]. Both

*A*

_{1}and

*A*

_{2}surfaces of the medium are Fresnel ones. The number of the photons is taken as 1.0 × 10

^{10}, and the angular discretization is taken as

*N*×

_{θ}*N*= 90 × 20. The results are plotted only for

_{φ}*φ*−

*φ*

_{0}=

*π*/ 2 (clockwise direction). The results for refractive index of

*n*= 1.44 is given in Figs. 9, 10, 11, 12, and 13, and those for

*n*= 1 is shown in Figs. 14 and 15.

In the following discussions, we use the concept of characteristic time for analyzing time-resolved results, defined as *t** = *nτ _{x}*/

*μ*, where

*τ*is optical thickness with respect to the normal distance

_{x}*x*from the surface

*A*

_{1}. From the law of refraction, we can get the critical angle for total reflection and the refraction angle at the internal surface. At the surface of the medium having

*n*= 1.44, the cosines of critical angle for total reflection and refraction angle are about 0.72 and 0.8, respectively. Thus for the case of

*n*= 1.44, there are several characteristic times including

*t** = 1.44 for

*μ*= 1 and

*τ*= 1,

_{x}*t** = 2 for

*μ*= 0.72 and

*τ*= 1,

_{x}*t** = 4 for

*μ*= 0.72 and

*τ*= 2,

_{x}*t** = 1.8 for

*μ*= 0.8 and

*τ*= 1, and

_{x}*t** = 2.8 for

*μ*= 0.8 when the photons emitted at

*t** = 1 (at the end of the pulse) just arrive at the surface

*A*

_{2}from

*A*

_{1}. For the case of

*n*= 1, the two characteristic times may be used for

*A*

_{2}:

*t** = 1 for

*μ*= 1 when the photons emitted at

*t** = 0 are just received by

*A*

_{2}, and

*t** = 2 for

*μ*= 1 when the photons emitted at

*t** = 1 (at the end of the pulse) are just received by

*A*

_{2}.

Figures 9 and 10 demonstrate contour maps of the Stokes vector on the surfaces *A*_{1}^{+} and *A*_{2}^{+} with respect to time and directions, respectively. If the scattering is not considered, the shortest time *t** reaching the surface *A*_{2} is equal to *nτ _{L}*/

*μ*from the surface

*A*

_{1}, and if considering scattering, the time is shorten to

*nτ*. Therefore, as shown in Fig. 9, from the beginning, the surface

_{L}*A*

_{1}

^{+}receives the photons and the Stokes vector is not equal to zero, although the values are very small. The values of Stokes vector increase with the time, and it is seen that from

*t** = 4, peaks (valleys) appear at the directions

*μ*= 0.72 and

*μ*= −0.72, and four obvious peak (valleys) areas corresponding to the four parameters of Stokes vector are gradually formed in a certain temporal and angular range of

*t**≥4 and |

*μ*|≤0.72, which is resulting from the total reflection at the surface

*A*

_{2}

^{+}. In addition, the temporal range of peaks (valleys) is determined by the pulse width.

At the surface *A*_{2}^{+}, some significant differences of the upward radiance (*μ* > 0) can be observed. As shown in Fig. 10, from *t** = 2, peaks (valleys) appear at the range of −0.5>*μ* ≥−0.72 for the downward radiance, while for the upward radiance (*μ* > 0), in a certain temporal range from *t** = 2, the peaks (valleys) appear in a bigger angular range of 0.5<*μ* ≤ 1. Much of the radiation from the total reflection with *μ* = 0.72 at the surface *A*_{1}^{+} arrives at the surface *A*_{2}^{+} at *t** = 2, and immediately contributes to the occurrence of the peaks (valleys) of upward Stokes vector at the surface *A*_{2}^{+} approximately in the range of 0.5<*μ* ≤ 1 not merely at *μ* = 0.72. In the following discussion, we analyze the temporal and angular distributions of Stokes vector, respectively.

Figures 11–14 show the curves of Stokes vector for three different moments and three different directions at *A*_{1}^{+} and *A*_{2}^{+} surfaces, respectively. From the Fig. 11, we can observe two sharp peaks having a symmetric distribution around *μ* = 0 at the critical angle for total reflection in each curve of Stokes parameters at *t** = 4 when the radiation from the total reflection happening at *A*_{1}^{+} just returns to the surface. While at *t** = 6, at the angle of total reflection, only abrupt changes are found, and two peaks appear at the other angle. At *t** = 4, only one peak is found at *μ* = 0. From the Fig. 11, it is also noticed that at the surface *A*_{1}^{+}, the Stokes parameters are all close to zero at *μ* = 1.

Figure 13 demonstrates the curves of Stokes parameters for the surface *A*_{2}^{+} at three different moments. As mentioned above, at *t** = 2 the radiation mainly from total reflection at the surface *A*_{1}^{+} just arrives at the surface *A*_{2}^{+}, so we analyze only the curves for *t** = 2 in the Fig. 13. Compared with the curves for *t** = 4 in the Fig. 11, an obvious difference can be observed in the curves for *t** = 2 in the Fig. 13. The sharp peaks appear at the angle of total reflection with *μ* = ± 0.72 in the curves for *t** = 4 in the Fig. 11, while in the curves for *t** = 2 in the Fig. 13, after abrupt rise or falling at *μ* = 0.72, approximately in the range of 1>*μ* ≥0.72, the values change gradually and maintain a high level. In addition, the Stokes parameters are approximately zero at *μ* = −1 at the surface *A*_{2}^{+}.

The curves of Stokes vector with respect to the time for three different directions are shown in Fig. 12 for *A*_{1}^{+} and Fig. 14 for *A*_{2}^{+}. From the Fig. 12, at *t** = 1 an abrupt change can be found in the curves for *A*_{1}^{+}, and approximately at *t** = 4 the radiance mainly from total reflection at *A*_{2}^{+} arrives at *A*_{1}^{+} and the values of Stokes vector increase rapidly with time, resulting in the occurrence of peaks. At the direction of *μ* = −0.602, owing to the gradually increasing effects of the backward scattering, the positive value of *V* component increases at *A*_{1}^{+} surface before *t** = 6. Due to the polarization effect caused by Fresnel reflection, the *V* component is converted into negative value at *A*_{2}^{+} surface. When the photons of circular polarization with negative value from total reflection arrive at *A*_{1}^{+} surface, the *V* component with positive value begins to fall approximately at *t** = 4.5, and therefore, a valley appears. As the photons leave the *A*_{1}^{+} surface gradually, the *V* component increases again till reaches the second peak, as shown in Fig. 12. After that the value of *V* component decreases again, which is resulting from the weakened pulse energy. From the Fig. 14, it is seen that at about *t** = 1.44, the surface *A*_{2}^{+} just receives photons with a small number, and at *t** = 1.8 the values of Stokes vector increase dramatically, causing peaks or valleys at *t** = 2.8.

In order to highlight the effect of refractive index on the transient vector radiative transfer, we present the curves of Stokes vector for the surface *A*_{2}^{+} for *n* = 1. Figure 15 shows the angular distributions of Stokes vector for three different moments, and Fig. 16 plots the curves of time-resolved Stokes vector for three different directions. At the boundary of the scattering medium with *n* = 1, no reflection will happen, which produces transient polarized transfer behavior of short pulse quite different from that for the case of *n* = 1.44 discussed above. From the Fig. 15, we can see that the value of Stokes vector is zero for *μ* < 0 as a result of no reflection and no scattering at *A*_{2}^{+}, and for the curves of *t** = 2.0, we can find an abrupt change in the values of four Stokes parameters at the incident direction with *μ* = 0.5. From the Fig. 16 it is observed that the surface *A*_{2}^{+} just receives photons at *t** = 1.0, and after a moment the values rise sharply with time until *t** = 2.0; after *t** = 2.0, the values continue to increase with a relative slow rate to the peaks and then fall rapidly.

In this simulation, the CPU time taken is about 12 hours for the case of *n* = 1, and for *n* = 1.44, it is about 37 hours. If the time shift and superposition principle is not applied in the MC model, it will take thousands of hours to get relatively satisfactory results for time-resolved polarized radiative transfer. From the above computation it is confirmed that the MCM_TSS is a strong numerical strategy for solving the transient vector radiative transfer.

With *ω* = 0.5 and *ω* = 1.0, effects of scattering albedo on the time-resolved Stokes vector are shown in Fig. 17. Except for the scattering albedo, other parameters used in Fig. 17 are taken the same values as Fig. 16. It can be observed in Fig. 17 that, for different value of scattering albedo, the peak value of Stokes vector appears at the same time. While for a bigger scattering albedo, the peak value is higher. It is owing to the fact that, for a higher value of albedo, more photons from incident direction are scattered to the surface *A*_{2}^{+}.

At the direction of *μ* = 0.707, the results of the time-resolved Stokes vector for different optical thickness are presented in Fig. 18. Except for the optical thickness, other parameters used in Fig. 18 are taken as the same values as Fig. 16. It can be seen from Fig. 18, with the increase in the optical thickness, the peak value of Stokes vector decreases significantly. Also, for the case with bigger optical thickness, owing to the stronger attenuation effect, the Stokes vector appears at a later time and the profile of the Stokes vector lasts longer.

## 4. Conclusion

In this work, we developed a MC model with a combination of time shift and superposition principle for solving transient vector radiative transfer in a scattering layer. According to the definition of reflectivity, we proposed a new criterion of reflection at Fresnel surface. We put forward the concept of transient vector radiative transfer matrix (TVRTM) used to describe the transient polarized behavior of radiative transfer. Using the TVRTM solved by the MC method, we obtained four time- and angular-resolved Stokes parameters. Time shift and superposition (TSS) principle was applied to the MC model and significantly improved the computational efficiency and accuracy in the simulation.

Adopting the model of transient polarized radiative transfer developed by MCM_TSS in the present work, we obtained the results for transient scalar radiative transfer and steady-state vector radiative transfer, respectively that were compared well with those available in reported literatures, which confirms the correctness and high efficiency of the present model.

Finally, we investigated the transient polarized behavior of radiative transfer in Mie scattering slab subjected to an inclined parallel incident short-pulse laser irradiation. In the simulation, the boundaries of the slab were assumed to be Fresnel surfaces, and two cases of *n* = 1.44 and *n* = 1 were taken into consideration. Temporal and angular distributions of Stokes vector were plotted and analyzed in detail. From the discussion of results some conclusions can be made. From the angular-resolved curves of Stokes vector, we find that in the curves for the same surface there are different directional distribution characteristics at different time as the surface receives different numbers of photons in the same direction at different time, and in these curves, there are sharp peaks or abrupt changes at the critical angle of total reflection. From the time-resolved curves of Stokes vector, we can observe some rapid changes in the curves at the characteristic times when the number of photons received by the surface has a sudden increase or decrease. For the medium having *n* = 1, due to no reflection at the boundary, some corresponding peaks or changes in value disappear in the curves, and the distribution characteristics with respect to time and direction become more simple compared with the case of *n*>1. In addition, the application of time shift and superposition principle to MCM for transient vector radiative transfer can improve the computational efficiency by hundreds of times.

## Acknowledgments

This work was supported by the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (Grant No. 51121004) and the National Natural Science Foundation of China (Grant No. 51176040).

## References and links

**1. **A. Ishimaru, S. Jaruwatanadilok, and Y. Kuga, “Polarized pulse waves in random discrete scatterers,” Appl. Opt. **40**(30), 5495–5502 (2001). [CrossRef] [PubMed]

**2. **X. D. Wang, L. V. Wang, C. W. Sun, and C. C. Yang, “Polarized light propagation through scattering media: time-resolved Monte Carlo simulations and experiments,” J. Biomed. Opt. **8**(4), 608–617 (2003). [CrossRef] [PubMed]

**3. **E. A. Sergeeva and A. I. Korytin, “Theoretical and experimental study of blurring of a femtosecond laser pulse in a turbid medium,” Radiophys. Quantum Electron. **51**(4), 301–314 (2008). [CrossRef]

**4. **S. C. Mishra, P. Chugh, P. Kumar, and K. Mitra, “Development and comparison of the DTM, the DOM and the FVM formulations for the short-pulse laser transport through a participating medium,” Int. J. Heat Mass Transfer **49**(11-12), 1820–1832 (2006). [CrossRef]

**5. **M. Sakami, K. Mitra, and P. F. Hsu, “Analysis of light pulse transport through two-dimensional scattering and absorbing media,” J. Quant. Spectrosc. Radiat. Transf. **73**(2-5), 169–179 (2002). [CrossRef]

**6. **J. M. Wang and C. Y. Wu, “Transient radiative transfer in a scattering slab with variable refractive index and diffuse substrate,” Int. J. Heat Mass Transfer **53**(19-20), 3799–3806 (2010). [CrossRef]

**7. **J. M. Wang and C. Y. Wu, “Second-order-accurate discrete ordinates solutions of transient radiative transfer in a scattering slab with variable refractive index,” Int. Commun. Heat Mass Transf. **38**(9), 1213–1218 (2011). [CrossRef]

**8. **M. Akamatsu and Z. X. Guo, “Ultrafast radiative heat transfer in three-dimensional highly-scattering media subjected to pulse train irradiation,” Numer. Heat Tranf. Anal. Appl. **59**, 653–671 (2011).

**9. **Z. X. Guo, S. Kumar, and K. C. San, “Multidimensional Monte Carlo simulation of short-pulse laser transport in scattering media,” J. Thermophys. Heat Transf. **14**, 504–511 (2000).

**10. **P. F. Hsu, “Effects of multiple scattering and reflective boundary on the transient radiative transfer process,” Int. J. Therm. Sci. **40**(6), 539–549 (2001). [CrossRef]

**11. **C. Y. Wu, “Monte Carlo simulation of transient radiative transfer in a medium with a variable refractive index,” Int. J. Heat Mass Transfer **52**(19-20), 4151–4159 (2009). [CrossRef]

**12. **P. Rath, C. S. Mishra, P. Mahanta, U. K. Saha, and K. Mitra, “Discrete transfer method applied to transient radiative transfer problems in participating medium,” Numer. Heat Tranf. Anal. Appl. **44**, 183–197 (2003).

**13. **R. Singh, S. C. Mishra, N. K. Roy, N. S. Shekhawat, and K. Mitra, “An insight into the modeling of short-pulse laser transport through a participating medium,” Numer Heat Tranf. B-Fundam. **52**(4), 373–385 (2007). [CrossRef]

**14. **J. C. Chai, “One-dimensional transient radiation heat transfer modeling using a finite-volume method,” Numer Heat Tranf. B-Fundam. **44**(2), 187–208 (2003). [CrossRef]

**15. **J. C. Chai, P. F. Hsu, and Y. C. Lam, “Three-dimensional transient radiative transfer modeling using the finite-volume method,” J. Quant. Spectrosc. Radiat. Transf. **86**(3), 299–313 (2004). [CrossRef]

**16. **S. C. Mishra, R. Muthukumaran, and S. Maruyama, “The finite volume method approach to the collapsed dimension method in analyzing steady/transient radiative transfer problems in participating media,” Int. Commun. Heat Mass Transf. **38**(3), 291–297 (2011). [CrossRef]

**17. **C. Y. Wu and S. H. Wu, “Integral equation formulation for transient radiative transfer in an anisotropically scattering medium,” Int. J. Heat Mass Transfer **43**(11), 2009–2020 (2000). [CrossRef]

**18. **S. H. Wu and C. Y. Wu, “Time-resolved spatial distribution of scattered radiative energy in a two-dimensional cylindrical medium with a large mean free path for scattering,” Int. J. Heat Mass Transfer **44**(14), 2611–2619 (2001). [CrossRef]

**19. **G. W. Kattawar and G. N. Plass, “Radiance and polarization of multiple scattered light from haze and clouds,” Appl. Opt. **7**(8), 1519–1527 (1968). [CrossRef] [PubMed]

**20. **R. D. M. Garcia and C. E. Siewert, “The F_{N} method for radiative transfer models that in include polarization effects,” J. Quant. Spectrosc. Radiat. Transf. **41**(2), 117–145 (1989). [CrossRef]

**21. **K. F. Evans and G. L. Stephens, “A new polarized atmospheric radiative transfer model,” J. Quant. Spectrosc. Radiat. Transf. **46**(5), 413–423 (1991). [CrossRef]

**22. **C. E. Siewert, “A discrete-ordinates solution for radiative-transfer models that include polarization effects,” J. Quant. Spectrosc. Radiat. Transf. **64**(3), 227–254 (2000). [CrossRef]

**23. **J. Lenoble, M. Herman, J. L. Deuze, B. Lafrance, R. Santer, and D. Tanre, “A successive order of scattering code for solving the vector equation of transfer in the earth’s atmosphere with aerosols,” J. Quant. Spectrosc. Radiat. Transf. **107**(3), 479–507 (2007). [CrossRef]

**24. **F. Xu, R. A. West, and A. B. Davis, “A hybrid method for modeling polarized radiative transfer in a spherical-shell planetary atmosphere,” J. Quant. Spectrosc. Radiat. Transf. **117**, 59–70 (2013). [CrossRef]

**25. **H. Ishimoto and K. Masuda, “A Monte Carlo approach for the calculation of polarized light: application to an incident narrow beam,” J. Quant. Spectrosc. Radiat. Transf. **72**(4), 467–483 (2002). [CrossRef]

**26. **X. D. Wang and L. V. Wang, “Propagation of polarized light in birefringent turbid media: A Monte Carlo study,” J. Biomed. Opt. **7**(3), 279–290 (2002). [CrossRef] [PubMed]

**27. **M. Xu, “Electric field Monte Carlo simulation of polarized light propagation in turbid media,” Opt. Express **12**(26), 6530–6539 (2004). [CrossRef] [PubMed]

**28. **R. Vaillon, B. T. Wong, and M. P. Mengüç, “Polarized radiative transfer in a particle-laden semi-transparent medium via a vector Monte Carlo method,” J. Quant. Spectrosc. Radiat. Transf. **84**(4), 383–394 (2004). [CrossRef]

**29. **C. Davis, C. Emde, and R. Harwood, “A 3-D polarized reversed monte carlo radiative transfer model for millimeter and submillimeter passive remote sensing in cloudy atmospheres,” IEEE Trans. Geosci. Rem. Sens. **43**(5), 1096–1101 (2005). [CrossRef]

**30. **J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express **13**(12), 4420–4438 (2005). [CrossRef] [PubMed]

**31. **J. N. Swamy, C. Crofcheck, and M. P. Mengüç, “A Monte Carlo ray tracing study of polarized light propagation in liquid foams,” J. Quant. Spectrosc. Radiat. Transf. **104**(2), 277–287 (2007). [CrossRef]

**32. **C. Cornet, L. C. Labonnote, and F. Szczap, “Three-dimensional polarized Monte Carlo atmospheric radiative transfer model (3DMCPOL): 3D effects on polarized visible reflectances of a cirrus cloud,” J. Quant. Spectrosc. Radiat. Transf. **111**(1), 174–186 (2010). [CrossRef]

**33. **M. Sakami and A. Dogariu, “Polarized light-pulse transport through scattering media,” J. Opt. Soc. Am. A **23**(3), 664–670 (2006). [CrossRef] [PubMed]

**34. **Y. A. Ilyushin and Y. P. Budak, “Analysis of the propagation of the femtosecond laser pulse in the scattering medium,” Comput. Phys. Commun. **182**(4), 940–945 (2011). [CrossRef]

**35. **S. Chandrasekhar, *Radiative Transfer* (Oxford University, 1950).

**36. **H. C. Van de Hulst, *Light Scattering by Small Particles* (Dover, 1981).

**37. **M. I. Mishchenko, L. D. Travis, and A. A. Lacis, *Scattering, Absorption, and Emission of Light by Small Particles* (Cambridge University, 2002).

**38. **H. H. Tynes, G. W. Kattawar, E. P. Zege, I. L. Katsev, A. S. Prikhach, and L. I. Chaikovskaya, “Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations,” Appl. Opt. **40**(3), 400–412 (2001). [CrossRef] [PubMed]

**39. **R. Green, *Spherical Astronomy* (Cambridge University, 1985).

**40. **B. A. Whitney, “Monte Carlo radiative transfer,” Arxiv preprint arXiv: 1104.4990 (2011).

**41. **K. Masuda and T. Takashima, “Computational accuracy of radiation emerging from the ocean surface in the model atmosphere–ocean system,” Pap. Meteorol. Geophys. **37**(1), 1–13 (1986). [CrossRef]

**42. **G. W. Kattawar and C. N. Adams, “Stokes vector calculations of the submarine light field in an atmosphere–ocean with scattering according to a Rayleigh phase matrix: effect of interface refractive index on radiance and polarization,” Limnol. Oceanogr. **34**(8), 1453–1472 (1989). [CrossRef]

**43. **P. W. Zhai, Y. X. Hu, J. Chowdhary, C. R. Trepte, P. L. Lucker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transf. **111**(7-8), 1025–1040 (2010). [CrossRef]

**44. **E. R. Sommersten, J. K. Lotsberg, K. Stamnes, and J. J. Stamnes, “Discrete ordinate and Monte Carlo simulations for polarized radiative transfer in a coupled system consisting of two medium with different refractive indices,” J. Quant. Spectrosc. Radiat. Transf. **111**(4), 616–633 (2010). [CrossRef]