Transient radiative transfer in a scattering slab considering polarization.

Transient radiative transfer in a scattering slab considering polarization Hongliang Yi,1 Xun Ben,1 and Heping Tan1,2,* 1

School of Energy Science and Engineering, Harbin Institute of Technology, 92 West Dazhi Street, Harbin 150001, China 2 [email protected] *[email protected]

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 OCIS codes: (010.5620) Radiative transfer; (140.7090) Ultrafast lasers; (290.5855) Scattering, polarization.

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). 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). 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). 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). 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). 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). 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). 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). 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).

#195986 - $15.00 USD Received 19 Aug 2013; revised 11 Oct 2013; accepted 14 Oct 2013; published 29 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026693 | OPTICS EXPRESS 26693

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). 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). 15. J. C. Chai, P. F. Hsu, and Y. C. Lam, “Three-dimensional transient radiative transfer modeling using the finitevolume method,” J. Quant. Spectrosc. Radiat. Transf. 86(3), 299–313 (2004). 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). 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). 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). 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). 20. R. D. M. Garcia and C. E. Siewert, “The FN method for radiative transfer models that in include polarization effects,” J. Quant. Spectrosc. Radiat. Transf. 41(2), 117–145 (1989). 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). 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). 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). 24. F. Xu, R. A. West, and A. B. Davis, “A hybrid method for modeling polarized radiative transfer in a sphericalshell planetary atmosphere,” J. Quant. Spectrosc. Radiat. Transf. 117, 59–70 (2013). 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). 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). 27. M. Xu, “Electric field Monte Carlo simulation of polarized light propagation in turbid media,” Opt. Express 12(26), 6530–6539 (2004). 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). 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). 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). 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). 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). 33. M. Sakami and A. Dogariu, “Polarized light-pulse transport through scattering media,” J. Opt. Soc. Am. A 23(3), 664–670 (2006). 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). 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). 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). 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).


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). 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).

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], FN method [20], adding-doubling method [21], discrete-ordinate (DO) method [22], successive-orderscattering (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: 1 ∂S(r, Ω ) ⋅ + Ω ⋅ ∇S(r, Ω ) = −κeS(r, Ω) + κa I b (r) +  Z(r, Ω' → Ω)S(r, Ω')d Ω, 4π ∂t c

(1)

where S is the Stokes vector; Ω is the discrete direction vector; r is the spatial coordinate vector; κe is the extinction matrix; κa is the absorption vector; Ib is the black body intensity; c is the velocity of light in the medium; t is the time; Z is the scattering phase matrix. In Eq. (1), the Stokes vector S is expressed by S = ( I , Q , U ,V ) , T

(2)

where I is the intensity, Q is the linear polarization aligned parallel or perpendicular to the zaxis, 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 κaIb describing the thermal radiation emitted by the medium in Eq. (1) is ignored for the cold medium considered in the present work. 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 tp, the number of pulses is p, the time interval between two


successive pulses is td, and the intensity of the pulse is I0. Each pulse emits N photons within the pulse width of tp, and every photon starts its trajectory at the time of t = ( p − 1 + Rt )t p + ( p − 1) td ,

(3)

Laser Intensity

where Rt is the uniform random number.

tp

td I0

0

Time

Fig. 1. Schematic of square pulse chain

The distance a photon travels without being absorbed, scattered, reflected or refracted is determined by Eq. (4a): L=−

ln(1 − RL )

κ

,

(4a)

where L is the free distance, κ is the extinction coefficient and RL is a pseudo random number with an uniform distribution. After travelling the free distance, the new coordinate position ( x ', y ', z ' ) of the photon is obtained by:  x ' = x + L sin θ cos ϕ ,   y ' = y + L sin θ sin ϕ ,  z ' = z + L cos θ , 

(4b)

where (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 tL, given by Eq. (5a): tL =

L nL , = c c0

(5a)

where n is the refractive index of the medium, and c0 and c are the speed of light in vacuum and medium, respectively. Then the new time is determined by Eq. (5b): t ' = t + tL .

(5b)

In this paper, the time t is normalized to t* as Eq. (6):

t * = c0κ 0t ,

(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 Si and the resulting scattering Stokes vector Ss determine a plane AOB, defined as the scattering plane. The incident Stokes vector Si and the 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 Ss (scattering Stokes vector) and the zaxis, 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: Ss = L(π − i2 )M( Θ)L( −i1 )Si ,

(7)

where the matrix M is the single scattering Mueller matrix defined in Van de Hulst [36]; L is a rotation matrix; i1 and i2 are the angles of rotation into and out of the scattering plane, respectively. The phase matrix Z in Eq. (1) is equal to L(π−i2)M(Θ)L(−i1). z

Ss

C

Si

i2 B

θ2

Θ

θ1

i1 A

O

y

ϕ2 ϕ1 ϕ2

ϕ1

x

Fig. 2. Geometry of Scattering plane and Meridian planes. The photon’s direction of propagation before and after scattering is (θ1,φ1) and (θ2,φ2) respectively.

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 M 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:  M1 M M( Θ) =  2  0   0

M2

0

M1

0

0 0

M3 −M 4

0  0  , M4   M3 

(8)

where Θ 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:


0 0 0 1  0 cos 2φ sin 2φ 0  . (9) L(φ ) =   0 − sin 2φ cos 2φ 0    0 0 1 0 The physical meaning of Eq. (7), expression of Stokes vector after scattering, can be well understood. First, the incident Stokes vector Si is rotated by the angle of i1 counterclockwise from the meridian plane AOC to the scattering plane AOB. Then, the incident Stokes vector Si interacts with the particle and is scattered into a new direction. Finally, rotate the scattering Stokes vector by the angle of π − i2 clockwise from the scattering plane AOB to the meridian plane BOC and obtain the scattering Stokes vector Ss. 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 i1 and Θ are obtained, the angles i2, θ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 i1 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: S* = S I = (1, Q I ,U I ,V I ) . T

(10)

The normalization of Stokes vector must be carried out for every state-changing event such as scattering, reflection and transmission. 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): Sref = RSi ,

(11a)

Stra = TSi ,

(11b)

where R and T are reflection matrix and transmission matrix, respectively. Note that Sref and Stra also need to be normalized as done in Eq. (10). For the semitransparent Lambertian surface, the reflection matrix and transmission matrix are given as follows: ρ 0 R = 0  0

0 0 0 0 0 0  , 0 0 0  0 0 0

(12a)


1 − ρ  0 T=  0   0

0 0 0 0 0 0  , 0 0 0  0 0 0

(12b)

where ρ 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 Rref with the ρ to judge whether the reflection takes place or not at the boundary. If Rref < ρ 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 = arc cos

(

)

1 − Rθ , ϕ r = 2π Rϕ ,

(13)

where 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):  r2 + r⊥2  2 2 1 r − r⊥ R(θi ) =   2 0   0

 0 0  0 0 , 2 Re{r⊥ r* } 2 Im{r r⊥* }   2 Im{r⊥ r* } 2 Re{r⊥ r*} 

r2 − r⊥2 r2 + r⊥2 0 0

(14a)

in which r⊥ =

ni cos(θi ) − nt2 − ni2 sin 2 (θi ) ni cos(θi ) + nt2 − ni2 sin 2 (θi )

, r =

nt2 cos(θi ) − ni nt2 − ni2 sin 2 (θi ) nt2 cos(θi ) + ni nt2 − ni2 sin 2 (θi )

, (14b)

where the subscripts ⊥ and denote the components of perpendicular and parallel to the incident plane, respectively; Ref{ } and Im{ } indicate the real and imaginary parts, respectively; ni and nt 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2 + t⊥2  2 2 1 nt cos(θ t )  t − t⊥ T(θi ) = 2 ni cos(θi )  0   0

  , 2 Re{t⊥ t* } 2 Im{tt⊥* }   2 Im{t⊥ t*} 2 Re{t⊥ t* }

t2 − t⊥2 t2 + t⊥2

0 0

0 0

0 0

(14c)

in which t⊥ =

2ni cos(θi ) ni cos(θi ) + n − n sin (θi ) 2 t

2 i

2

, t =

2ni nt cos(θi ) n cos(θi ) + ni nt2 − ni2 sin 2 (θi ) 2 t

. (14d)

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 , nt sin θ t = ni sin θi , ϕ r = ϕ t = ϕi ,

(15)


where the subscripts 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 Rref 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):

ρ=

r2 + r⊥2 2

(16)

.

From the Eq. (16) we can see that the reflection coefficient of Fresnel surface for scalar radiative transfer has nothing to do with the incident photon. While for the vector radiative transfer, the reflection coefficient is associated with the incident photon. According to the definition, the reflection coefficient is equal to the ratio of the intensity for reflected wave to that for incident wave, denoted by ρ ' , written as Eq. (17a):

ρ'=

I' , I

(17a)

where I ' is the intensity for reflected wave and I is the intensity for incident wave. We can get I ' from the Eq. (11a) and Eq. (14a), and then ρ ' is obtained, expressed by Eq. (17b):

ρ'=

1 2 2 ( r + r⊥ ) + 2QI ( r2 − r⊥2 ) = ρ + 2QI ( r2 − r⊥2 ) . 2

(17b)

We can see if 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 Di,t,θ,φ is defined as follows: Nt

Di ,t ,θ ,ϕ =

S q =1

Δt

q

N tp

,

(18)

in which N indicates the number of photons emitted within the pulse width of tp; Nt denotes the total number of the photons passing through the cross-section i along (θ, φ) direction in the solid angle dΩ from t to t + Δt moment; Sq is the normalized Stokes vector carried by photon 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: Si ,t ,θ ,ϕ =

I 0 cos θ 0 cos θ dΩ

⋅ Di ,t ,θ ,ϕ ,

(19)

in which I0 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 tp into NΔt time steps, and each time step Δt = tp / NΔt. We only need to solve the Stokes vector response caused by the photons emitted at the time zero within the time step of Δt according to the Eq. (18) and Eq. (19) letting tp = Δt, denoted by S’, 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 Si,t,θ,φ can be calculated by Eq. (20): Si ,t ,θ ,ϕ =

1 min ( N Δt , int ( t Δt ) )

min ( N Δt ,int ( t Δt ) )

 j =1

Si' ,( int ( t Δt )− j )Δt ,θ ,ϕ .

(20)

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 steadystate 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 tp* = 1.0 is loading on the A1− 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 S0 = (1, 0, 0, 0)T. The medium with refractive index n = 1.0 is isotropic scattering and non-absorbing (ω = 1.0). Both A1 and A2 are assumed Lambertian boundaries. According to [16] the dimensionless reflectance and transmittance are defined as Eq. (21) R (t ) =

I

2π

A1− ,t ,θ ,ϕ

cos θ d Ω

I 0 cos θ 0

, T (t ) =

I

2π

A2− , t ,θ ,ϕ

cos θ d Ω

I 0 cos θ 0

.

(21)


A2− A2+

τ

A1+ A1−

θ S0

Fig. 3. Schematic of the geometry of the model.

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 × 106 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 × 108, while once the TSS principle is considered in the MC method, only 1.0 × 106 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 × 108 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. ×10

−2

1.0

I-Reflectance

0.8

6

MC_TSS N=1.0×10 FVM Ref. [16]

4.0 3.0

I-Transmittance

5.0

τ = 0.1

2.0 1.0 0.0

0

1

2

3

t* ×10

0

0.8

τ = 2.0

0.6 0.4 0.2 2

4

6 t*

8

1

2

(b)

4

−1

×10

6


2.0 1.5

τ = 2.0

1.0 0.5 0.0

10

3

t*

2.5 I-Transmittance

I-Reflectance

6


1.0

0

0.2

3.0

1.2

0.0

τ = 0.1

(a)

−1

1.4

0.4

0.0

4

6


0.6

0

2

4

6

8

10

t*

Fig. 4. Comparison of transmittance and reflectance signals with respect to dimensionless time for different optical thickness: (a) τ = 0.1; (b) τ = 2.0.


5.0 ×10

−2

I-Reflectance

4.0

6

MC_TSS N=1.0×10 8 MC N=1.0×10

3.0 2.0

τ = 0.1

1.0 0.0

0

1

2

3

4

t*

Fig. 5. Comparison of computational accuracy between MCM and MCM_TSS. Table 1. Comparison of computation time by MCM and MCM_TSS MC_TSS, N = 1.0 × 106 2s 8s

τ 0.1 2.0

MC, N = 1.0 × 108 81s 293s

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 atmospherewater 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 S0 = (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 × 109 and the angular discretization is taken as Nθ × Nφ = 60 × 20 for the comparison. To satisfy the steady-state working condition, we set tp* = ∞ in the transient vector MC model.

Atmosphere Interface

Ocean

Bottom

Fig. 6. Schematic of the coupled atmosphere-water system with an oblique incident unpolarized beam.


−1

−2

1.0 ×10

3.5

by ρ' by ρ DOM [44]

0.8 0.6

×10

by ρ' by ρ DOM [44]

3.0 2.5

I

Q

2.0

0.4

1.5 1.0

0.2

ϕ = 90°

0.0 -1.0

-0.5

0.0

ϕ = 90°

0.5

0.5

0.0 -1.0

1.0

-0.5

cosθ

(a)

-2

6.0

×10

by ρ by ρ' DOM [44]

3.0 2.0 ϕ = 90°

1.0 0.0 -1.0

-0.5

0.0 cosθ

(c)

1.0

×10

4.0 3.0 V

U

4.0

0.5

(b)

−4

5.0

5.0

0.0 cosθ

2.0

ϕ = 90°

by ρ by ρ' DOM [44]

1.0 0.5

1.0

0.0 -1.0

-0.5

0.0

0.5

1.0

cosθ

(d)

Fig. 7. Stokes vector elements just above the atmosphere–water interface for a collimated unpolarized incident beam.

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 ρ ' 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 ρ ' agree well with those in [44] by DOM, and compared with the results by ρ, the results by ρ ' are a little better. Thus in the following simulation, only the results by ρ ' 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.


1

1

0.080 0.074 0.068 0.062 0.056 0.050 0.044 0.038 0.032 0.026 0.020 0.014 0.008 0.002

0

-0.5

-1

0

5

10

t*

15

20

0.5

cos(polar angle)

cos(polar angle)

0.5

0.033 0.030 0.027 0.024 0.021 0.018 0.015 0.012 0.009 0.006 0.003 0.000

0

-0.5

-1

25

0

5

10

1

0.049 0.045 0.041 0.037 0.033 0.029 0.025 0.021 0.017 0.013 0.009 0.005 0.001

20

25

0

-0.5

0

5

10

t*

(c)

15

20

25

0.5

cos(polar angle)

cos(polar angle)

15

1

0.5

-1

t*

(b)

(a)

4.4E-04 4.0E-04 3.6E-04 3.2E-04 2.8E-04 2.4E-04 2.0E-04 1.6E-04 1.2E-04 8.0E-05 4.0E-05 0.0E+00

0

-0.5

-1

0

5

10

t*

15

20

25

(d)

Fig. 8. Distributions of Stokes vector just above the atmosphere–water interface varying with time and direction.

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 A1− surface, and the pulse width is tp* = 1. The pulse is a polarized beam of irradiance S0 = (1,0, 0.5, 0.5,0.5)T. The optical thickness τL with respect to the layer thickness 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 A1 and A2 surfaces of the medium are Fresnel ones. The number of the photons is taken as 1.0 × 1010, 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 timeresolved results, defined as t* = nτx/μ, where τx is optical thickness with respect to the normal distance x from the surface A1. 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 τx = 1, t* = 2 for μ = 0.72 and τx = 1, t* = 4 for μ = 0.72 and τx = 2, t* = 1.8 for μ = 0.8 and τx = 1, and t* = 2.8 for μ = 0.8 when the photons emitted at t* = 1 (at the end of the pulse) just arrive at the surface A2 from A1. For the case of n = 1, the two characteristic times may be used for A2: t* = 1 for μ = 1 when the photons emitted at t* = 0 are just received by A2, and t* = 2 for μ = 1 when the photons emitted at t* = 1 (at the end of the pulse) are just received by A2.


1

1

0.0245 0.0225 0.0205 0.0185 0.0165 0.0145 0.0125 0.0105 0.0085 0.0065 0.0045 0.0025 0.0005

0

-0.5

-1

0

2

4

6

8

0.5

cos(polar angle)

cos(polar angle)

0.5

0.0028 0.0024 0.0020 0.0016 0.0012 0.0008 0.0004 0.0000 -0.0004 -0.0008 -0.0012 -0.0016 -0.0020

0

-0.5

-1

10

0

2

4

t*

(a)

10

1

0.001 0.000 -0.001 -0.002 -0.003 -0.004 -0.005 -0.006 -0.007 -0.008 -0.009 -0.010 -0.011

0

-0.5

0

2

4

6

t*

8

0.5

cos(polar angle)

0.5

cos(polar angle)

8

(b)

1

-1

6

t*

0.0028 0.0020 0.0012 0.0004 -0.0004 -0.0012 -0.0020 -0.0028 -0.0036 -0.0044 -0.0052

0

-0.5

-1

10

0

2

4

6

8

10

t*

(c)

(d)

Fig. 9. Distributions of Stokes vector on the surface A1+ varying with time and direction for n = 1.44. 1

1

0.061 0.056 0.051 0.046 0.041 0.036 0.031 0.026 0.021 0.016 0.011 0.006 0.001

0

-0.5

-1

0

2

4

6

8

0.5

cos(polar angle)

cos(polar angle)

0.5

0.0100 0.0085 0.0070 0.0055 0.0040 0.0025 0.0010 -0.0005 -0.0020 -0.0035 -0.0050

0

-0.5

-1

10

0

2

4

t*

10

1

-0.002 -0.005 -0.008 -0.011 -0.014 -0.017 -0.020 -0.023 -0.026 -0.029 -0.032 -0.035 -0.038

0

-0.5

cos(polar angle)

0.5

0.5

cos(polar angle)

8

(b)

1

-1

6

t*

(a)

0

-0.5

-1

0

2

4

t*

(c)

6

8

10

0.026 0.023 0.020 0.017 0.014 0.011 0.008 0.005 0.002 -0.001 -0.004 -0.007

0

2

4

6

8

10

t*

(d)

Fig. 10. Distributions of Stokes vector on the surface A2+ varying with time and direction for n = 1.44.

Figures 9 and 10 demonstrate contour maps of the Stokes vector on the surfaces A1+ and A2 with respect to time and directions, respectively. If the scattering is not considered, the shortest time t* reaching the surface A2 is equal to nτL/μ from the surface A1, and if considering scattering, the time is shorten to nτL. Therefore, as shown in Fig. 9, from the beginning, the surface A1+ 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 A2+. In addition, the temporal range of peaks (valleys) is determined by the pulse width. At the surface A2+, 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.51. 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).


Single-scattering solutions to radiative transfer in infinite turbid media.

One-dimensional transient radiative transfer by lattice Boltzmann method.

Ultrafast radiative heat transfer.

Interference-exact radiative transfer equation.

Polarization recovery through scattering media.

Comparative analysis of radiative transfer approaches for calculation of diffuse reflectance of plane-parallel light-scattering layers.

Polarized radiative transfer in an arbitrary multilayer semitransparent medium.

Approximate solution to vector radiative transfer in gradient-index medium.

Radiative and nonradioactive electron transfer in donor-acceptor phenanthrimidazoles.

Nanoantenna enhanced emission of light-harvesting complex 2: the role of resonance, polarization, and radiative and non-radiative rates.

Near-field radiative heat transfer between metamaterial thin films.

Quantitative analysis of light scattering in polarization-resolved nonlinear microscopy.

Combining the absorptive and radiative loss in metasurfaces for multi-spectral shaping of the electromagnetic scattering.

Decoupling scattering and absorption of turbid samples using a simple empirical relation between coefficients of the Kubelka-Munk and radiative transfer theories.

Inelastic electron scattering from a helical potential: transverse polarization and the structure factor in the single scattering approximation.

Polarization-controlled optimal scatter suppression in transient absorption spectroscopy.

Geometric phase coded metasurface: from polarization dependent directive electromagnetic wave scattering to diffusion-like scattering.

Carrier-induced transient defect mechanism for non-radiative recombination in InGaN light-emitting devices.

Polarization-controlled directional scattering for nanoscopic position sensing.

A path to practical Solar Pumped Lasers via Radiative Energy Transfer.

Comparative analysis of radiative transfer approaches for calculation of plane transmittance and diffuse attenuation coefficient of plane-parallel light scattering layers.

A Note on Radiative Heat Transfer to Peristaltic Flow of Sisko Fluid.

Energy streamlines in near-field radiative heat transfer between hyperbolic metamaterials.

Influence of Orientation and Radiative Heat Transfer on Aluminum Foams in Buoyancy-Induced Convection.