## Abstract

Transient (time-dependent) polarized radiative transfer in a scattering medium exposed to an external collimated beam illumination is conducted based on the time-dependent polarized radiative transfer theory. The transient term, which persists the nanosecond order time and cannot be ignored for the time-dependent radiative transfer problems induced by a short-pulsed beam, is considered as well as the polarization effect of the radiation. A discontinuous finite element method (DFEM) is developed for the transient vector radiative transfer problem and the derivation of the discrete form of the governing equation is presented. The correctness of the developed DFEM is first verified by comparing the DFEM solutions with the results from the literature. The DFEM is then applied to study the transient polarized radiative transfer induced by a pulsed beam. The time-dependent Stokes vector components are calculated, plotted and analyzed as functions of the axis coordinate and discrete direction. Effects of the diffuse/specular boundary and the incident beam polarization state with respect to the Stokes vector components are further analyzed for cases of different boundary reflection modes and incident beam illuminations.

© 2017 Optical Society of America

## 1. Introduction

Radiative transfer of polarized light in a participating medium has drawn much attention on atmospheric radiation [1] and optical applications [2]. Research work on the fields of optical imaging [3, 4], astronomy and planetary science [5], and target detection in turbid media [6] can also benefit from exploiting information associated with the polarized radiation. In particular, it is noted that the vector radiative transfer in scattering media induced by a short-pulsed beam has significance in the fields of biomedical diagnosis and fast remote sensing both in theoretical research and engineering applications. The polarization effect of the radiation can be accurately described by the vector radiative transfer equation (VRTE). The polarized radiation contains four Stokes vector components depending on both the scattering zenith angle and azimuth angle, which makes the solution of the VRTE quite complex, thereby limiting the applications of polarized radiation information. In recent decades, researchers have developed several numerical methods to investigate the vector radiative transfer process in scattering media. Garcia and Siewert [7, 8] reported a solution for all four Stokes vector components using a generalized spherical harmonics method and the *F _{N}* method respectively. Evans and Stephens [9] solved the VRTE for a plane-parallel medium with the doubling and adding method and listed the Legendre series coefficients for the scattering phase matrix. A partial list of the numerical methods for VRTE includes the successive order of scattering (SOS) method [10], discrete ordinates method (DOM) [11–13], discontinuous finite element method (DFEM) [14], the matrix operator method [15] and the Monte Carlo method (MCM) [16, 17]. Kokhanovsky et al. [18] compared several different numerical methods for solving vector atmospheric radiative transfer and, benchmarked results with a 1° viewing zenith angle grid for different viewing azimuth angles. Very recently, the International Polarized Radiative Transfer (IPRT) of the International Radiation Commission (IRC) [19] initiated a model competition and the results for the one-dimensional plane-parallel model were summarized.

Nowadays, the transient effect induced by ultrafast laser beam interactions with participating media have attracted significant curiosity due to its widespread applications in microstructure fabrication, biomedical and laser writing areas [20, 21]. Transients persist for a short time, the order of a nanosecond for typical length scales, after which the radiation energy reaches a steady state inside the medium as the radiation propagates at the speed of light. For the quasi-steady problems, the transient term associated with time in the complete radiative transfer equation (RTE) can be omitted. However, for nanosecond pulses, the time for the radiation to reach a steady state is of the same order with that of the beam duration, and thus the non-steady terms in the RTE cannot be ignored [20]. The time-dependent radiative transfer theory should be adopted, in which the transient radiative transfer equation (TRTE) for time-resolved radiation distributions needs to be solved. For the scalar radiative transfer problems, some analytical and numerical results are available [22–32]. Su and Olson [22] gave an analytical solution for non-equilibrium radiative transfer in an infinite medium by applying the Fourier transform with respect to the spatial variable and the Laplace transformation with respect to the temporal variable. Oliveira and associates [23] combined the spectral method and the Laplace transform method to get the semi-analytical results for time-dependent radiative transfer problem in a slab-geometry medium. Tan and Hsu [24] developed a time-dependent integral formulation for modeling transient radiative transfer in participating media and some other researchers studied the transient scalar radiative transfer problems by using various numerical methods [25–32]. For the polarized radiative transfer problems, owing to the complexity of the transient vector radiative transfer equation (TVRTE), only a few studies [33–37] take both the polarization and transient effects into account. The propagation mechanism and, the energy distributions of the external polarized beam illumination in the scattering medium have not been fully investigated. The time-resolved Stokes vector components need investigation, and the effects of the pulse duration, reflection mode of the medium boundary and incident polarization state on the behavior of time-dependent signals need to be examined in detail.

The discontinuous finite element method (DFEM) first proposed by Reed and Hill [38] is an accurate numerical procedure that combines the salient features of both the finite element method (FEM) and finite volume method (FVM) [39]. The DFEM was first applied to solve the radiative transfer equation by Cui and Li [40]. Afterward, Liu and Hsu [41] extended the DFEM to solve the scalar radiative transfer equation within a graded index medium. The previous work shows the potential of DFEM for solving radiative transfer problems and our purpose in this paper is to apply the DFEM to solve the TVRTE in a one-dimensional parallel slab containing a Mie scattering atmosphere and to carry forward research on the problem by revealing the physical behavior and characteristics of transient polarized radiative transfer in the scattering medium.

The outline of this paper is as follows. In the following section, the mathematical formulations of the TVRTE are introduced and the DFEM discretization of the TVRTE is given. In Section 3, the correctness of the DFEM is verified by comparing the DFEM results with those in reported literature. Afterward, the DFEM is applied to solve the TVRTE in the scattering atmosphere for a short-pulsed beam. Different boundary reflection modes are investigated and the time-resolved Stokes vector components are analyzed. Finally, some conclusions are drawn in the last section.

## 2. Governing equation and DFEM discretization

We consider a polarized collimated beam illumination of radiation as shown in Fig. 1. The pulsed beam has a peak Stokes vector of **I**_{0} = (*I*_{0}, *Q*_{0}, *U*_{0}, *V*_{0})^{T}, is incident on the left boundary at a polar angle of *θ*_{0} and azimuthal angle of *φ*_{0}. The duration of the beam, which is of a nanosecond order, is denoted by *t _{p}*. Based on the incoherent addition principle of Stokes parameters and, considering the polarization and transient effects, the radiative transfer equation can be written as

*c*

_{0}is the light speed in the vacuum,

*z*is the axis coordinate,

*t*is the time,

**I**= (

*I*,

*Q*,

*U*,

*V*)

^{T}is the Stokes vector,

**Ω**=

*μ*

**i**+

*η*

**j**+

*ξ*

**k**= sin

*θ*cos

*φ*

**i**+ sin

*θ*sin

*φ*

**j**+ cos

*θ*

**k**is the unit direction vector,

**β**is the extinction coefficient matrix, and

**S**is the source term written as

**Z**is the scattering phase matrix for scattering from an incoming direction

**Ω**' to an outgoing direction

**Ω**.

Considering diffuse and specular reflection, the boundary conditions are given by

**Ω**” denotes the corresponding incident direction of the current reflected beam of direction

**Ω**,

**n**

_{w}denotes the unit outward normal vector of the boundary.

**R**

_{s}and

**R**

_{d}are the specular and diffuse reflection matrix [9].

The short-pulsed collimated beam considered in this work is a square pulse arriving at the left boundary at *t* = 0. The radiation incident on the boundary at *z* = 0 can be written as

*H*(

*t*) is the Heaviside step function.

For the medium exposed to an external collimated beam, the Stokes vector inside the medium is usually divided into two parts

where**I**

_{c}denotes the collimated part (Stokes vector of the direct radiation) and

**I**

_{d}denotes the diffuse part (Stokes vector of the diffuse radiation).

The collimated part **I**_{c} within the medium decreases exponentially according to Beer's law and can be derived mathematically as

The diffuse part **I**_{d} inside the medium can be obtained by solving Eq. (1) after modifying the right-hand source term as

Without any ambiguity, the subscript *d* in **I**_{d} is omitted in the following expressions. As the results in terms of dimensionless variables are more general in the transient radiative transfer field [24, 29–32], we use the dimensionless time *t** = *c*_{0}*t*/*L* to analyze the propagation characteristics of the beam propagation and radiation distribution in this paper. The radiative transfer equation for the diffuse radiation can be written as

By using the central time differencing, Eq. (8) can be written as

Substitute Eq. (10) into Eq. (8) we can get the TVRTE at each time step

*k*denotes the

*k*th time step, $\tilde{\beta}$ and ${\tilde{S}}_{k}$ are defined as

For the DFEM application to the TVRTE, the one-dimensional parallel slab is divided into *N _{elem}* non-overlapped uniform adjacent elements, and the angular space is discretized into

*M*=

*N*×

_{θ}*N*directions. The subscript

_{φ}*k*is omitted in the following expression for simplification and Eq. (10) can be written in discrete form as

*n*th discrete direction. For each element

*e*, Eq. (13) is weighted by

*W*and integrated over

*K*using Gauss divergence theorem

The symbol *∂e* in Eq. (14) denotes the boundary of the element *e* and, **n*** _{∂e}* denotes the unit outward normal vector of the boundary, [

**I**] denotes the discontinuous radiation on the boundary associated with two adjacent elements. By applying the local Lax-Friedrichs numerical scheme, a classical time-front, space-center difference scheme based on finite difference first proposed by Lax in 1954 [42], the discontinuous radiation on the boundary can be expressed as

*e*, the elements, element boundaries and the radiation values on the boundaries in DFEM discretization are sketched in Fig. 2.

By substituting Eq. (16) into Eq. (14) we can obtain the general DFEM formulation of TVRTE as

By using the shape function on each element *K*, an approximated solution of **I*** ^{n}* over an element

*e*can be assumed in the form as

*i*and

*ϕ*is the shape function. By substituting Eq. (19) into Eq. (18) and taking the shape function

_{i}*ϕ*as the weight function, we can obtain the DFEM discretization for the TVRTE on element

_{j}*e*in matrix formwhere the matrix ${{\rm K}}^{n}$ and ${H}^{n}$are defined as

Up to now, the DFEM discretization for the TVRTE has been completed. As a result, the Stokes vector components at each time step can be obtained by calculating all the component matrices in Eq. (21) and solving Eq. (20) by Gaussian elimination on each discrete element.

## 3. Results and discussions

As the source term in Eq. (21)b) contains the radiative information of the other directions, global iterations are necessary to update the source term. In this paper the maximum relative error 10^{−6} of incident radiation is taken as the stop criterion for the global iteration. A Matlab code is constructed to calculate the matrix terms in Eq. (21) and solve the transient polarized radiative transfer problems. The correctness of the developed code is verified by comparing the DFEM results for a transient scalar and a quasi-steady polarized radiative transfer problem with the reported results. Then the DFEM is applied to study the transient polarized radiative transfer in the one-dimensional atmosphere investigated by Garcia and Siewert [8]. The time-dependent Stokes vector component distributions are plotted and, effects of the boundary reflection mode and the incident beam's polarization state are further investigated.

#### 3.1 Validation of the discontinuous finite element method

To validate the correctness of the procedure code developed in this paper, we first compare the DFEM results for a scalar radiative transfer problem against the FVM results [26]. As shown in Fig. 1, a short-pulse beam source with a pulse width of *t _{p}** = 1.0 is vertically exposed to the left boundary of a purely scattering medium. The incident beam is unpolarized with a uniform irradiance of (1, 0, 0, 0)

^{T}. By setting all off-diagonal elements of the phase matrix

**Z**to be zero, polarization can be ignored and the DFEM procedure developed in this paper is applied to the transient scalar radiative transfer problem. The computational domain is divided into

*N*= 20 uniform elements, and the angular space is divided into

_{elem}*N*×

_{θ}*N*= 40 × 60 directions using the piecewise constant approximation (PCA) uniform angular discretization according to our recent work on solving the steady vector radiative transfer equation [14], and the time step

_{φ}*Δt** = 0.05 is found to be sufficient to obtain stable and satisfactory results. The reflectance on the left boundary and the transmittance on the right boundary are plotted in Fig. 3 for different optical thickness, namely

*τ*=

*βL*= 0.1, and 1.0. It can be seen from Fig. 3 that the DFEM results agree well with the FVM results, which shows the correctness and accuracy of the DFEM to handle the time-dependent term in the radiative transfer equation.

A polarized radiative transfer is further studied to examine the correctness of the DFEM for the polarization effect. The quasi-steady solutions of vector radiative equilibrium in a one-dimensional scattering atmosphere with a diffuse boundary have been reported [8]. Thus, the present transient solutions at long enough time can be compared with the quasi-steady solutions to examine the validity of the present scheme. Consider the one-dimensional medium exposed to a uniform external illumination as shown in Fig. 1, the left boundary is transparent and the right boundary is a reflector. Garcia and Siewert [8] carried out the quasi-steady vector radiative transfer solutions for a cold atmosphere with optical thickness *τ* = 1.0 and a single scattering albedo *ω* = 0.99. The collimated beam illumination with linear polarization state and a flux intensity **I**_{0} = *π* × (1, 0, 0, 0)^{T}, is incident on the left boundary in the direction *μ*_{0} = cos*θ*_{0} = 0.2, *φ*_{0} = *π*/2. The right boundary is diffused with a reflectivity *ρ*_{d} = 0.1. The scattering phase matrix is that of Mie scattering at a wavelength of 0.951 μm from a gamma distribution of particles and the phase matrix can be written as

*P*

_{1},

*P*

_{2},

*P*

_{3}and

*P*

_{4}are the four unit elements. The Legendre coefficients for the four elements are presented by Evens and Stephens [9], and the four unique elements of the scattering phase matrix varying with scattering angles are presented in Fig. 4. In this paper, we consider a continuous external beam incident on the left boundary at time

*t** = 0 and apply the transient DFEM model to solve this problem.

Figure 5 plots the diffuse Stokes vector components (*I*_{*}, *Q*_{*}, *U*_{*}, *V*_{*}) at position *z*/*L* = 0.5 for the viewing azimuth angle 0° to compare with the quasi-steady results [8] and all the following results are diffuse Stokes vector components for the viewing azimuth angle 0°. As it takes dimensionless time *t** = 0.5 for external light photons to reach position *z*/*L* = 0.5, the time-resolved Stokes components after *t** = 0.5 (*t** = 1.0, 2.0, 5.0, 10.0, 15.0) are plotted to compare with the quasi-steady results. It can be seen from Fig. 5 that the distribution trend varies with time and the peak values of the Stokes component lines at a different time appear at different direction cosines. The curves of the transient solutions gradually approach that of the quasi-steady solutions and the transient solutions at a large *t**, we say *t** ≥ 10, completely overlap with the quasi-steady solutions, which means the energy distributions inside the medium will reach a steady state for *t** ≥ 10. The DFEM solutions approach and overlap with the quasi-steady solutions demonstrates the correctness of the transient DFEM model for solving vector radiative transfer equation in the scattering atmosphere.

#### 3.2 Transient vector radiative transfer in an atmosphere exposed to a continuous collimated beam illumination

The atmosphere exposed with a linear polarized beam illumination is considered in this section, the computation parameters are same with the polarized validation case in section 3.1 and thus not repeated here. Figure 6 plots the time-resolved Stokes vector component results at three different positions *z*/*L* = 0, 0.5, 1.0. For the position *z*/*L* = 0, the incident laser just begins to enter the medium and has not yet been scattered, all of the diffuse radiation from inside the medium with a propagation direction cosine *μ* < 0 will be refracted into the outside environment at the left transparent boundary, and all the four Stokes components remain zero for the direction *μ* > 0. The Stokes vector component values at the position *z*/*L* = 0 begin to increase from *t** = 0, have apparent changes at *t** = 0-2.0 and gradually become stable at *t** = 10.0 as shown in Fig. 6(a). These distributions at position *z*/*L* = 0.5 plotted in Fig. 6(b) have non-zero values for all discrete directions due to the scattering of the atmosphere. For the Stokes vector components at *z*/*L* = 1.0 plotted in Fig. 6(c), as the right boundary is diffused reflection, only the Stokes vector component *I _{*}* (radiative intensity) is reflected, the values for Stokes vector components

*Q*,

_{*}*U*and

_{*}*V*remain zero for the directions

_{*}*μ*< 0. The

*I*values for the directions

_{*}*μ*< 0 are the same at a given time step because the radiation energy is reflected uniformly in the hemisphere space at the diffuse boundary.

The polar angle space [0, *π*] is discretized into *N _{θ}* = 40 sub-angles. The discrete angles are ${\theta}_{i}=(i-1/2)\pi /{N}_{\theta}$, and the discrete directions

*θ*

_{5},

*θ*

_{10}

*θ*

_{15}and their symmetrical directions about

*μ*= 0, i.e.,

*θ*

_{26},

*θ*

_{31},

*θ*

_{36}are chosen for further analysis. The direction cosines and the scattering angles with the incident beam [incident polar angle

*θ*

_{0}= arccos (0.2) = 78.46°] are listed in Table 1. The reflected and transmitted Stokes components for the chosen directions with direction cosines

*μ*= ± 0.4187, ± 0.7343, and ± 0.9382 are plotted in Fig. 7.

For the reflected Stokes vector components plotted in Fig. 7 (a), the values of *I _{*}* and

*Q*at

_{*}*z*/

*L*= 0 gradually increase from time

*t** = 0 due to the multiple scattering effects of the light bundles induced by the incident beam illumination. For the three chosen discrete directions

*μ*= −0.4187, −0.7343, and −0.9382, it can be seen from Fig. 4 that the scattering angle is smallest for the direction

*μ*= −0.4187, the element

*P*

_{1}in the scattering phase matrix is biggest for

*μ*= −0.4187, and then most radiation will be scattered into

*μ*= −0.4187. This leads to a result that the intensity of the Stokes vector component

*I*is strongest for the direction

_{*}*μ*= −0.4187 and weakest for the direction

*μ*= −0.9382. A similar situation can be found for the Stokes vector component

*Q*as it is transformed from

_{*}*I*, which can be seen from the scattering phase matrix Eq. (22). In contrast, the value of the Stokes vector component

_{*}*U*is largest for

_{*}*μ*= −0.9382 and smallest for the direction

*μ*= −0.4187. It can be also found from Fig. 6(a) that the Stokes vector component curves do not cross with each other for the selected discrete directions, as the

*I*and

_{*}*Q*curves with largest (smallest for

_{*}*U*and

_{*}*V*curves) steady values increase/decrease at the fastest rate.

_{*}The transmitted Stokes vector components at *z*/*L* = 1.0 are shown in Fig. 7(b). It can be seen that the values of all the Stokes vector components begin to appear at *t** = 1.0 when the first light photon reaches the right boundary. The intensity of Stokes vector component *I _{*}* is strongest for the direction

*μ*= 0.4187 due to the smallest scattering angle and biggest

*P*

_{1}value in the scattering phase matrix. The

*I*and

_{*}*Q*curves gradually increase while

_{*}*U*curves gradually decrease, and the

_{*}*I*,

_{*}*Q*and

_{*}*U*curves change at different rates before they reach constant asymptotic values. It can be found that the

_{*}*V*curves for

_{*}*μ*= 0.7343 and

*μ*= −0.9382 increase at first and then decrease before reaching stable values, while the

*V*curve for

_{*}*μ*= 0.4187 decreases at first and then increases to a constant value. As a result, the transmitted Stokes vector components curves for different directions cross each other.

Figure 8 plots the radiative flux *q* distributions for Stokes vector components *I _{*}* (${q}_{{I}_{*}}$) and

*Q*(${q}_{{Q}_{*}}$) at four selected times

_{*}*t** = 0.5, 1.0, 2.0, 10.0. For

*t** = 0.5, the incident beam illumination just propagates to the location

*z*/

*L*= 0.5 and thus, the radiative fluxes only have values in the area from

*z*/

*L*= 0 to

*z*/

*L*= 0.5, consistent with the curves for

*t** = 0.5 in Fig. 8. The radiative flux of Stokes vector component

*I*is larger for the region from

_{*}*z*/

*L*= 0 to

*z*/

*L*= 0.2, but smaller for the region from

*z*/

*L*= 0.55 to

*z*/

*L*= 1.0 at earlier times. This is mainly due to the multiple scattering of the energy bundles inside the medium. For the locations near the left boundary (

*z*/

*L*= 0), there are fewer bundles inside the medium for an earlier time because of the shorter incident beam duration, and fewer bundles are propagating in the backward directions with

*μ*< 0. At later times, the number of bundles with

*μ*< 0 increases naturally and, the radiative flux of Stokes vector component

*I*decreases. In contrast, for the locations near to the right boundary (

_{*}*z*/

*L*= 1.0), most of the bundles propagate in the forward directions with

*μ*> 0, the radiative flux values of Stokes vector component

*I*are positive and increase with time as more bundles propagate to the locations near to the right boundary. Similar trends can be found for the radiative flux of Stokes vector component

_{*}*Q*. At early times, the radiative flux of Stokes vector components

_{*}*Q*is bigger for the region from

_{*}*z*/

*L*= 0 to

*z*/

*L*= 0.15, while smaller for the region from

*z*/

*L*= 0.42 to

*z*/

*L*= 1.0.

#### 3.3 Transient vector radiative transfer in an atmosphere exposed to a short-pulsed beam illumination

In this case, we consider the transient vector radiative transfer in the atmosphere exposed to a short-pulsed beam illumination. The dimensionless duration of the beam illumination is *t _{p}** = 0.5 and the other computational parameters are the same with those in section 3.1. Figure 9 plots the Stokes vector component distributions as a function of time and zenith angle cosine. As the incident beam is of short duration, the source term

**I**

*in Eq. (7) will disappear after time ${t}^{*}=1.5$, when the last incident collimated beam propagates through the medium, resulting in the rapid decrease of the diffuse radiation and in the intensity of Stokes vector component*

_{c}*I*, as manifested in Fig. 9. The Stokes vector components

_{*}*I*,

_{*}*Q*and

_{*}*U*last longer for the directions with

_{*}*μ*= 0 compared with that for other directions. The changing trend for

*V*curve is much complicated and distinctive in spite of its small magnitude. For the location

_{*}*z*/

*L*= 0, the

*V*value around direction

_{*}*μ*= −0.1 during time period

*t** = [0.4, 2.0] is positive and relatively large. There exists an ellipsoidal area where the

*V*value is negative for location

_{*}*z*/

*L*= 0. Similar areas with teardrop and triangle shapes can be found for locations

*z*/

*L*= 0.5 and

*z*/

*L*= 1.0 respectively.

Figure 10 plots the time-resolved Stokes vector components at locations *z*/*L* = 0 and *z*/*L* = 1.0 for the directions *μ* = ± 0.4187, ± 0.7343, and ± 0.9382. For the reflected Stokes vector components at *z*/*L* = 0 plotted in Fig. 10(a), the *I _{*}* and

*Q*curves climb gradually before they decrease quickly from

_{*}*t** = 0.5 when the incident beam illumination disappears, leading to a sharp peak at

*t** = 0.5. It is noted that the Stokes vector component

*I*for

_{*}*μ*= −0.9382 increases again from

*t** ≈2.0, when the reflected radiation by the right boundary travels back to the left boundary. The phenomenon is not so apparent for

*μ*= −0.4187 and

*μ*= −0.7343 since there are not so many reflected bundles travelling in the directions

*μ*= −0.4187 and

*μ*= −0.7343. The

*U*curves exhibit an opposite trend compared to those of

_{*}*I*and

_{*}*Q*curves, and the

_{*}*V*curves have smooth peaks instead of sharp ones.

_{*}For the transmitted Stokes vector components at *z*/*L* = 1.0 plotted in Fig. 10(b), all the curves start to grow from *t** = 1.0, when the incident beam first reaches the right boundary. For the direction *μ* = 0.9382, the *I _{*}* and

*Q*curves increase most rapidly and reach their highest values at

_{*}*t** = 1.5. For the direction with smaller direction cosines, i.e., the direction having bigger intersection angle with the positive

*z*-axis direction, the

*I*and

_{*}*Q*curves climb more slowly, reaching a lower peak value at a later time. The transmitted

_{*}*U*curve for

_{*}*μ*= 0.7343 has the lowest value, and the

*U*curve for

_{*}*μ*= 0.4187 reaches the valley peak at the latest time. The

*V*curves for

_{*}*μ*= 0.7343 and

*μ*= 0.9382 increase at first and then decrease, while the

*V*curve for

_{*}*μ*= 0.4187 first decreases and then increases before

*t** = 3.0. All of the Stokes vector components approach zero after

*t** = 3.0 due to the continual decline of the diffuse radiation.

The radiative flux distributions at different times are plotted in Fig. 11 to analyze the energy distributions along the optical thickness (*z*-axis) direction. For the time *t** ≤ 1.0, the radiative flux curves of *I _{*}* increase continuously before reaching a peak at locations

*z*/

*L*= 0.2, 0.3 and 0.45 for time

*t** = 0.5, 0.75 and 1.0, respectively. The wavefront of the incident beam propagates to the locations

*z*/

*L*= 0.5, 0.75 and 1.0 for the times

*t** = 0.5, 0.75 and 1.0 respectively, leading to the result that the radiative flux curves of

*I*is only appreciable before

_{*}*z*/

*L*= 0.5, 0.75 and 1.0 for the time

*t** = 0.5, 0.75 and 1.0. When

*t** ≥ 1.5, the incident short-pulse beam has propagated through the right boundary of the medium, the collimated source term in Eq. (7) disappears and the radiative flux drops quickly. It can be seen from Fig. 10 that the radiative flux values of

*I*for

_{*}*t** = 2.0, 5.8, and 8.0 are much smaller than those for

*t** ≤ 1.0, and the radiative fluxes of

*I*are linearly distributed along the

_{*}*z*-axis. A similar phenomenon can be found for the radiative flux curves of

*Q*, but it is noted that the radiative flux curve of

_{*}*Q*at

_{*}*t** = 2.0 is not linearly distributed, which demonstrates that the Stokes vector component

*Q*responds slower than

_{*}*I*does to the disappearance of the incident beam.

_{*}#### 3.4 Transient vector radiative transfer in an atmosphere with a specular boundary and exposed to a short-pulsed beam illumination

The atmosphere-ocean system is vital for remote sensing applications. The surface of the ocean can be seen as a specular boundary for the one-dimensional atmosphere-ocean system and thus we consider the vector radiative transfer in the atmosphere with a specular reflection boundary in this case. The computational parameters are the same with those in section 3.3 except for the reflection mode of the right boundary, which is modeled as a Fresnel boundary with a reflectivity *ρ*_{s} = 0.1 and a complex refractive index *m* = 3.724-2.212*i* for water at 27 °C [10]. The Stokes vector component distributions at locations *z*/*L* = 0, 0.5, and 1.0 are plotted in Fig. 12. Obvious differences can be found compared with the Stokes vector components for the medium with a diffuse boundary. For the Stokes vector components at location *z*/*L* = 0, as shown in Fig. 12(a), the *V _{*}* distribution for the medium with a specular boundary is most different from that for the medium with diffuse boundary, the area of relatively large

*V*values for the medium with specular boundary starts from

_{*}*t** = 2.5, which is much later than that for the medium with diffuse boundary. This is due to the reason that the polarization state of the reflected bundles by the specular right boundary is dramatically changed after the specular reflection event so that the reflected bundles have appreciable

*V*values and contribute to the

_{*}*V*distribution. From Fig. 12(b), the influence of the reflected bundles to all the four Stokes vector components at location

_{*}*z*/

*L*= 0.5 is significant, and this influence is most obvious for the Stokes vector components at the specular reflection boundary

*z*/

*L*= 1.0 plotted in Fig. 12(c).

Figure 13 compares the Stokes vector components at *z*/*L* = 0, *μ* = −0.4187 and those at *z*/*L* = 1.0, *μ* = 0.4187 for media with different reflection boundaries. It can be seen from Fig. 13(a) that the reflected Stokes vector component curves for the media with different reflection boundaries overlap exactly with each other before *t** = 2.0. This is due to the fact that the left boundary first begins to receive the bundles reflected by the right boundary at *t** = 2.0, the reflected bundles and the reflection mode of the right boundary do not have any influence on the Stokes vector components distributions before *t** = 2.0. The Stokes vector component curves for media with different reflection boundaries begin to separate after time *t** = 2.0 when the radiation reflected by the right boundary travels back to the left boundary. The reflection mode of the right boundary has the most obvious effects on the Stokes vector component *V _{*}*, which coincides with the analysis results from Fig. 12. A similar situation and corresponding explanation can be applied to the transmitted Stokes vector component curves plotted in Fig. 13(b). From the results for the medium with specular boundary, we reach the conclusion that the reflection mode of the boundary is an important parameter for the transient vector radiative transfer and it changes the Stokes vector component

*V*in a significant magnitude.

_{*}#### 3.5 Transient vector radiative transfer in an atmosphere exposed to a short-pulsed beam of circular polarization

An incident short-pulsed beam illumination with circular polarization (CP) and flux intensity **I**_{0} = *π* × (*I*_{0}, *Q*_{0}, *U*_{0}, *V*_{0})^{T} = *π* × (1, 0, 0, 1)^{T} is investigated and compared with the above-mentioned linear polarization (LP) beam to study the effects of the polarization state on the time-resolved Stokes vector components distributions. The other computational parameters are the same with those in section 3.3.

Figure 14 plots the four Stokes vector components for the linear and circular polarization beams. The *I _{*}*, and

*Q*curves for different incident beams coincide as the

_{*}*I*value for different incident beams are the same and the component

_{0}*Q*is transformed from the component

*I*. The

*U*curves for different incident beams have similar trends but different peak values, while the

_{*}*V*curves for different incident beam are different, in both trends and peak magnitudes. For the Stokes vector components at location

_{*}*z*/

*L*= 0, as shown in Fig. 14(a), the component

*U*for the circular polarization beam (CPB) is slightly higher than that for the linear polarization beam (LPB). The component

_{*}*V*for the CPB increases linearly during the beam duration

_{*}*t** = 0 ~0.5 because the

*V*value in the incident CPB equals to 1 and the diffuse bundles increase continually during the beam duration. As

_{0}*V*

_{0}= 1.0 in the incident CPB while

*V*

_{0}= 0 in the incident LPB, the Stokes vector component

*V*for the LPB only accumulates from the scattering bundles whose polarization state is changed by the scattering event, leading to the fact that the

_{*}*V*curve for LPB is much smaller than that for CPB. The cases of Stokes vector components distributions at location

_{*}*z*/

*L*= 1.0 plotted in Fig. 14(b) are similar to those for location

*z*/

*L*= 0.

Figure 15 plots the radiative fluxes of the Stokes vector components for the atmosphere with a diffuse boundary and a circular polarized beam illumination. The radiative fluxes for all the four Stokes vector components have similar trends. Due to the fact that *V*_{0} = 1.0 in the incident beam, the radiative flux for the Stokes vector component *V _{*}* has a relatively large value compared to the other three components due to the contribution of the incident beam.

## 4. Conclusions

Transient polarized radiative transfer in a scattering atmosphere exposed to an external beam illumination is numerically solved by a discontinuous finite element method (DFEM). The transient solutions at large time are compared with the quasi-steady solutions to examine the validity of the present scheme. The transient results approach and finally coincide with the steady results, which demonstrates the correctness of the developed DFEM method. The time-dependent Stokes vector components distributions as a function of time and discrete directions for a Mie-scattering atmosphere exposed to a collimated beam are presented and analyzed. The following remarkable conclusions can be drawn.

- (I) Stokes vector components in different directions are compared together and it is found that the elements distributions with the scattering angle in the scattering phase matrix have a significant influence to the Stokes vector component distributions.
- (II) Transient polarized radiative transfer in an atmosphere with a diffuse boundary and a short-pulsed beam is studied. For the atmosphere with a linear-polarized beam, the radiative flux distributions before the beam propagates through the medium are much larger and have different trends compared with those after the beam propagates through the atmosphere.
- (III) Transient polarized radiative transfer in an atmosphere with a specular boundary and a short-pulsed beam is studied. Results show that the specular boundary has a significant influence on the Stokes vector component distributions, especially for the transmitted Stokes vector components at location
*z*/*L*= 1.0 and for the Stokes vector component*V*at all locations._{*} - (IV) Effects of the polarization state of the incident beam are investigated. Results show that the Stokes vector component
*V*is the one most different between the linear and circular cases._{*}

## Funding

National Natural Science Foundation of China (NSFC) (51422602).

## References and links

**1. **S. Chandrasekhar, *Radiative transfer* (Dover, 1960).

**2. **D. S. Kliger and J. W. Lewis, *Polarized light in optics and spectroscopy* (Elsevier, 2012).

**3. **M. Moscoso, J. B. Keller, and G. Papanicolaou, “Depolarization and blurring of optical images by biological tissue,” J. Opt. Soc. Am. A **18**(4), 948–960 (2001). [CrossRef] [PubMed]

**4. **A. J. Brown, “Spectral bluing induced by small particles under the Mie and Rayleigh regimes,” Icarus **239**, 85–95 (2014). [CrossRef]

**5. **A. J. Brown, T. I. Michaels, S. Byrne, W. Sun, T. N. Titus, A. Colaprete, M. J. Wolff, G. Videen, and C. J. Grund, “The case for a modern multi-wavelength, polarization-sensitive LIDAR in orbit around Mars,” J. Quant. Spectrosc. Radiat. Transf. **153**, 131–143 (2015). [CrossRef]

**6. **S. A. Kartazayeva, X. Ni, and R. R. Alfano, “Backscattering target detection in a turbid medium by use of circularly and linearly polarized light,” Opt. Lett. **30**(10), 1168–1170 (2005). [CrossRef] [PubMed]

**7. **R. D. M. Garcia and C. E. Siewert, “A generalized spherical harmonics solution for radiative transfer models that include polarization effects,” J. Quant. Spectrosc. Radiat. Transf. **36**(5), 401–423 (1986). [CrossRef]

**8. **R. D. M. Garcia and C. E. Siewert, “The *F _{N}* method for radiative transfer models that include polarization effects,” J. Quant. Spectrosc. Radiat. Transf.

**41**(2), 117–145 (1989). [CrossRef]

**9. **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]

**10. **J. Lenoble, M. Herman, J. L. Deuzé, B. Lafrance, R. Santer, and D. Tanré, “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]

**11. **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]

**12. **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 media with different refractive indices,” J. Quant. Spectrosc. Radiat. Transf. **111**(4), 616–633 (2010). [CrossRef]

**13. **D. Cohen, S. Stamnes, T. Tanikawa, E. R. Sommersten, J. J. Stamnes, J. K. Lotsberg, and K. Stamnes, “Comparison of discrete ordinate and Monte Carlo simulations of polarized radiative transfer in two coupled slabs with different refractive indices,” Opt. Express **21**(8), 9592–9614 (2013). [CrossRef] [PubMed]

**14. **C. H. Wang, H. L. Yi, and H. P. Tan, “Discontinuous finite element method for vector radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. **189**, 383–397 (2017). [CrossRef]

**15. **X. Q. He, Y. Bai, Q. Zhu, and F. Gong, “A vector radiative transfer model of coupled ocean–atmosphere system using matrix-operator method for rough sea-surface,” J. Quant. Spectrosc. Radiat. Transf. **111**(10), 1426–1448 (2010). [CrossRef]

**16. **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]

**17. **T. Yun, N. Zeng, W. Li, D. Li, X. Jiang, and H. Ma, “Monte Carlo simulation of polarized photon scattering in anisotropic media,” Opt. Express **17**(19), 16590–16602 (2009). [CrossRef] [PubMed]

**18. **A. A. Kokhanovsky, V. P. Budak, C. Cornet, M. Z. Duan, C. Emde, I. L. Katsev, D. A. Klyukov, S. V. Korkin, L. C-Labonnote, B. Mayer, Q. Min, T. Nakajima, Y. Ota, A. S. Prikhach, V. V. Rozanov, T. Yokota, and E. P. Zege, “Benchmark results in vector atmospheric radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. **111**(12), 1931–1946 (2010). [CrossRef]

**19. **C. Emde, V. Barlakas, C. Cornet, F. Evans, S. Korkin, Y. Ota, L. C. Labonnote, A. Lyapustin, A. Macke, B. Mayer, and M. Wendisch, “IPRT polarized radiative transfer model intercomparison project – Phase A,” J. Quant. Spectrosc. Radiat. Transf. **164**, 8–36 (2015). [CrossRef]

**20. **S. Kumar and K. Mitra, “Microscale aspects of thermal radiation transport and laser applications,” Adv. Heat Transf. **33**, 187–294 (1999). [CrossRef]

**21. **M. Malinauskas, A. Žukauskas, S. Hasegawa, Y. Hayasaki, V. Mizeikis, R. Buividas, and S. Juodkazis, “Ultrafast laser processing of materials: from science to industry,” Light Sci. Appl. **5**(8), e16133 (2016). [CrossRef]

**22. **B. Su and G. L. Olson, “An analytical benchmark for non-equilibrium radiative transfer in an isotropically scattering medium,” Ann. Nucl. Energy **24**(13), 1035–1055 (1997). [CrossRef]

**23. **J. V. P. De Oliveira, A. V. Cardona, M. T. Vilhena, and R. C. Barros, “A semi-analytical numerical method for time-dependent radiative transfer problems in slab geometry with coherent isotropic scattering,” J. Quant. Spectrosc. Radiat. Transf. **73**(1), 55–62 (2002). [CrossRef]

**24. **Z. M. Tan and P. F. Hsu, “An integral formulation of transient radiative transfer,” ASME J. Heat Transfer **123**(3), 466–475 (2001). [CrossRef]

**25. **K. M. Katika and L. Pilon, “Modified method of characteristics in transient radiation transfer,” J. Quant. Spectrosc. Radiat. Transf. **98**(2), 220–237 (2006). [CrossRef]

**26. **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]

**27. **Q. Cheng, H. C. Zhou, Z. F. Huang, Y. L. Yu, and D. X. Huang, “The solution of transient radiative transfer with collimated incident serial pulse in a plane-parallel medium by the DRESOR method,” ASME J. Heat Transfer **130**(10), 102701 (2008). [CrossRef]

**28. **Z. X. Guo, J. Aber, B. A. Garetz, and S. Kumar, “Monte Carlo simulation and experiments of pulsed radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. **73**(2), 159–168 (2002). [CrossRef]

**29. **H. L. Yi, C. H. Wang, and H. P. Tan, “Transient radiative transfer in a complex refracting medium by a modified Monte Carlo simulation,” Int. J. Heat Mass Transfer **79**, 437–449 (2014). [CrossRef]

**30. **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), 3799–3806 (2010). [CrossRef]

**31. **P. Rath, S. C. Mishra, P. Mahanta, U. K. Saha, and K. Mitra, “Discrete transfer method applied to transient radiative transfer problems in participating medium,” Numer. Heat Transf. A **44**(2), 183–197 (2003). [CrossRef]

**32. **C. H. Wang, Y. Zhang, H. L. Yi, and M. Xie, “Analysis of transient radiative transfer induced by an incident short-pulsed laser in a graded-index medium with Fresnel boundaries,” Appl. Opt. **56**(7), 1861–1871 (2017). [CrossRef] [PubMed]

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

**34. **X. 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]

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

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

**37. **H. Yi, X. Ben, and H. Tan, “Transient radiative transfer in a scattering slab considering polarization,” Opt. Express **21**(22), 26693–26713 (2013). [CrossRef] [PubMed]

**38. **W. H. Reed and T. Hill, “Triangular mesh methods for the neutron transport equation,” Los. Alamos Rep. LA-UR-73-479 (1973).

**39. **J. S. Hesthaven and T. Warburton, *Nodal discontinuous Galerkin methods: algorithms, analysis, and applications* (Springer Science & Business Media, 2007).

**40. **X. Cui and B. Q. Li, “Discontinuous finite element solution of 2-D radiative transfer with and without axisymmetry,” J. Quant. Spectrosc. Radiat. Transf. **96**(3), 383–407 (2005). [CrossRef]

**41. **L. H. Liu and L. J. Liu, “Discontinuous finite element method for radiative heat transfer in semitransparent graded index medium,” J. Quant. Spectrosc. Radiat. Transf. **105**(3), 377–387 (2007). [CrossRef]

**42. **P. D. Lax, “Weak solutions of nonlinear hyperbolic equations and their numerical computation,” Commun. Pure Appl. Math. **7**(1), 159–193 (1954). [CrossRef]