3D-PSTD simulation and polarization analysis of a light pulse transmitted through a scattering medium

Fabrice Devaux; Eric Lantz

doi:10.1364/OE.21.024969

1. Introduction

In a recent paper, we have presented results of ultrafast compensation of turbidity of ex-vivo biological tissues by type 2 three-wave-mixing phase conjugation (TWMPC), where images transmitted through biological tissues with thicknesses up to 5 mm were restored [1]. These results indicate that the scattering process in such biological samples is more or less independent of the polarization state of the light and that a polarization change of the phase conjugated wave with respect to the incident wave does not preclude the image restoration process. Indeed, the phase conjugated wave (i.e. the idler wave) that retraces the scattering path is not polarized as the scattered light exiting from the biological tissues (i.e. the signal), because in a type 2 three-wave-mixing process the signal and the idler wave are polarized in orthogonal directions. The hypothesis of polarization insensitivity, first confirmed by a simple numerical model based on the Monte Carlo method, deserves further study with a more efficient numerical model. In fact, we have to implement a numerical model able to simulate accurately, in time and in space, electromagnetic phenomena such as the propagation of light in a scattering medium and the non linear optical process of TWMPC. Although the Monte Carlo offers a flexible and accurate method approach to model the variation in time of the state of polarization of the light transmitted through scattering media [2], non linear optical phenomena and coherent effects of light propagation are not accessible with this method.

The most popular method to simulate transient electromagnetic wave propagation is the finite-difference time domain (FDTD) method, where the spatial derivatives in Maxwell’s equations are approximated by finite differences [3]. Even if FDTD has been developed to treat a wide range of problems involving interaction of electromagnetic waves with all kind of materials, the numerical dispersion characteristics associated with FDTD can be problematic especially when large-scale problems are considered. In these cases, a space sampling density of at least few tens of cells per wavelength is necessary to ensure that FDTD methods produce acceptable results. Therefore, 3D sampling of a scattering medium with realistic dimensions (with dimensions of few tens of micrometers in each direction) will require a very large number of sampling points and very powerful computing resources. An alternative solution has been proposed by Liu [4]. This method, called the pseudo-spectral time-domain (PSTD) method, differs from the FDTD by the fact that spatial derivatives are calculated using Fourier transforms. This spatial differential process converges with infinite order of accuracy for grid-sampling densities of two or more points per wavelength, provided that the medium optical properties are sampled in accordance with the Nyquist theorem. In consequence, it allows the study of various problems on larger scales, more efficiently (with a factor 4^D to 8^D where D is the dimensionality of the problem) and with a better accuracy than FDTD methods [4–6]. In particular, the problem of light scattering was addressed with 2D or 3D-PSTD algorithms in order to calculate the scattered far fields [7] and the particle single-scattering properties with sizes up to 200 wavelengths [8] or for comparison with the Mie theory [9] and the discrete dipole approximation [10]. The problem of time reversal was also studied in two dimensions using the PSTD [11, 12]. In this case, the phase conjugation was obtained by a simple inversion of the magnetic field and not with the equations governing the non-linear effect producing phase conjugation. In any case, the polarization state and the time of flight of the scattered electromagnetic fields were not discussed.

In this paper, we propose a 3D-PSTD algorithm to model the propagation of a light pulse through scattering media with realistic dimensions. The variation in time and the state of polarization of the transmitted light are analyzed for different characteristics of the scattering medium and different polarization states of the incident light. To validate our algorithm, we compare the results obtained from our simulations with a Monte carlo numerical method and with some previously reported results. We show that, while the incident light is strongly depolarized on the whole, the light that reaches the output face of the scattering medium is much less depolarized. This simple fact could explain the results reported in [1].

2. Principle of the tri-dimensional pseudo-spectral time domain method

In PSTD algorithms, Maxwell’s curl equations are calculated with discrete Fourier transforms in order to solve the spatial derivatives on an unstaggered, collocated grid [13]. The Fast-Fourier-Transform (FFT) is used to implement these spatial derivatives and limitations of FFT, due to the periodic boundary conditions, are eliminated by using absorbing boundary conditions formulated for perfect matched layer (PML) in nonconductive media [4, 14, 15].

2.1. Tri-dimensional pseudo-spectral time domain scheme

We consider here the propagation along the z axis of a transverse, pulsed electromagnetic wave in a linear, lossless, non dispersive, non conductive and inhomogeneous medium. Because a non dispersive medium is considered, we assume that light is monochromatic or, more precisely, that its Gaussian time envelope is Fourier transform. According to the Maxwell’s equations, the components along the (x, y, z) axes of the electric displacement field D⃗, the electric field E⃗ and the magnetic field B⃗ are calculated using the time-stepping iterations [4]. We give here the details of the time-variation calculations for the x components of the different fields. The other components can be deduced by a simple circular permutation of the x, y, z indices. The equations governing electromagnetic fields in the medium are given by :

{\begin{array}{l} \frac{\partial D_{x}}{\partial t} = & \frac{\partial D_{x y}}{\partial t} + \frac{\partial D_{x z}}{\partial t} = & \frac{1}{μ_{0}} [\frac{\partial B_{z}}{\partial y} - \frac{\partial B_{y}}{\partial z}] \\ E_{x} = & \frac{D_{x}}{ε} \\ \frac{\partial B_{x}}{\partial t} = & \frac{\partial B_{x y}}{\partial t} + \frac{\partial B_{x z}}{\partial t} = & - \frac{\partial E_{z}}{\partial y} + \frac{\partial E_{y}}{\partial z} \end{array}

D_x and B_x are split respectively into two terms [4] where the second indices y and z are related to the y, z derivatives so that $\frac{\partial D_{x y}}{\partial t} = \frac{1}{μ_{0}} \frac{\partial B_{z}}{\partial y}$ , $\frac{\partial D_{x z}}{\partial t} = - \frac{1}{μ_{0}} \frac{\partial B_{y}}{\partial z}$ and $\frac{\partial B_{x y}}{\partial t} = - \frac{\partial E_{z}}{\partial y}$ , $\frac{\partial B_{x z}}{\partial t} = \frac{\partial E_{y}}{\partial z}$ . As we use the FFT to approximate the spatial derivatives, the y derivative of a general field component, for example A_z that is known at all grid points (j, k, l), can be computed as :

{\frac{\partial A_{z}}{\partial y}}_{j k l} \approx {F_{y}^{- 1} (- 2 i π ν_{y} F_{y} (A_{z}))}_{j k l},

where ν_y, F_y and

F_{y}^{- 1}

denote, respectively, the spatial frequency, the forward FFT and the inverse FFT along the y direction. The use of the approximation given by Eq. (2) to calculate the spatial derivatives in Eq. (1) yields backward Euler time-stepping relations of the following form:

{\begin{matrix} {D_{x} |}_{j k l}^{n + 1} = & \frac{{D_{x y} |}_{j k l}^{n + \frac{1}{2}} + \frac{Δ t}{μ_{0}} F_{y}^{- 1} (- 2 i π ν_{y} F_{y} ({B_{z} |}_{j k l}^{n}))}{1 + Δ t {γ_{y} |}_{j k l}} + \frac{{D_{x z} |}_{j k l}^{n + \frac{1}{2}} - \frac{Δ t}{μ_{0}} F_{z}^{- 1} (- 2 i π ν_{z} F_{z} ({B_{y} |}_{j k l}^{n}))}{1 + Δ t {γ_{z} |}_{j k l}} \\ {E_{x} |}_{j k l}^{n + 1} = & \frac{{D_{x} |}_{j k l}^{n + 1}}{ɛ_{j k l}} \\ {B_{x} |}_{j k l}^{n + 1} & \frac{{B_{x y} |}_{j k l}^{n - \frac{1}{2}} - Δ t F_{y}^{- 1} (- 2 i π ν_{y} F_{y} ({E_{z} |}_{j k l}^{n}))}{1 + Δ t {γ_{y} |}_{j k l}} + \frac{{B_{x z} |}_{j k l}^{n + \frac{1}{2}} + Δ t F_{z}^{- 1} (- 2 i π ν_{z} F_{z} ({E_{y} |}_{j k l}^{n}))}{1 + Δ t {γ_{z} |}_{j k l}} \end{matrix}

where (j, k, l) are cartesian coordinates indices of the spatial sampling point (jΔx, kΔy, lΔz) in the sampled volume and n is the index of the temporal sampling point nΔt. γ_v|_jkl and ε_jkl are, respectively, the boundary absorbing layer function along the v direction [4, 15] and the permittivity of the medium at the spatial sampling point (jΔx, kΔy, lΔz). The absorbing layers are obtained by setting the term γ_v|_jkl to zero inside the region of interest and to a value increasing linearly inside the absorbing domains [4]. In these equations, we assume that the permeability and the permittivity of vacuum are equal to 1 (μ₀ = 1, ε₀ = 1) so that the speed of light in vacuum is normalized : c = 1. According to the Nyquist sampling theorem, the spatial sampling step is fixed to a value corresponding to two or more samples per wavelength in a single direction [4] : Δx = Δy = Δz = 0.3λ. With c = 1, the time step is fixed to

Δ t = \frac{Δ x}{16 c} = \frac{Δ x}{16}

, which is smaller than the maximum time step

\frac{2 Δ x}{π c \sqrt{3}}

corresponding to the stability limit of the 3D-PSTD algorithm [13]. The time-stepping relations for the other components of the electromagnetic fields can be easily deduced from Eq. (3) by a circular permutation of the x, y, z indices. The source terms are designed such as the transverse electromagnetic wave emitted by the source propagates along the z direction and added, at the intermediate time-step

n + \frac{1}{2}

, to the terms

{D_{x z} |}_{j k l}^{n}

and

{D_{y z} |}_{j k l}^{n}

using the relations :

{\begin{matrix} {D_{x z} |}_{j k l}^{n + \frac{1}{2}} = & {D_{x z} |}_{j k l}^{n} + {S_{x} |}_{j k l}^{n} \\ {D_{y z} |}_{j k l}^{n + \frac{1}{2}} = & {D_{y z} |}_{j k l}^{n} + {S_{y} |}_{j k l}^{n} \end{matrix}

${S_{x} |}_{j k l}^{n}$ and ${S_{y} |}_{j k l}^{n}$ are the transverse amplitude components of the pulsed source term (delayed with a time t₀ and with a time width σ_t) at the time step nΔt and at the spatial sampling point (jΔx, kΔy, lΔz). These components are given by:

{\begin{matrix} {S_{x} |}_{j k l}^{n} = & {S_{0} |}_{j k l} cos ψ e^{- i (ω n Δ t + φ_{x})} e^{- \frac{{(n Δ t - t_{0})}^{2}}{2 σ_{t}^{2}}} \\ {S_{y} |}_{j k l}^{n} = & {S_{0} |}_{j k l} sin ψ e^{- i (ω n Δ t + φ_{y})} e^{- \frac{{(n Δ t - t_{0})}^{2}}{2 σ_{t}^{2}}} \end{matrix}

ψ, φ_x and φ_y are the parameters used to define the polarization state of the electromagnetic wave emitted by the source. The spatial amplitude S₀|_jkl is designed with a Gaussian shape in the (x, y) transverse plane and with the optimized three-cells normalized pattern [ $\frac{1}{4}$ , $\frac{1}{2}$ , $\frac{1}{4}$ ] along the z axis in order to suppress the aliasing errors [16]. Indeed, this optimized three-cells normalized pattern, according to the properties of Pascal’s triangle, has a spatial frequency spectrum which is a discrete decreasing function becoming null at the cutoff spatial frequency. Finally, the propagation of the wave in the increasing z direction is ensured by adding the source terms simultaneously on the magnetic field so that $\vec{B} = \sqrt{\frac{μ_{0}}{ε}} {\vec{e}}_{z} \times \vec{D}$ . It is obtained as follows:

{\begin{matrix} {B_{x z} |}_{j k l}^{n + \frac{1}{2}} = & {B_{x z} |}_{j k l}^{n} - \sqrt{\frac{μ_{0}}{ε_{j k l}}} {S_{y} |}_{j k l}^{n} \\ {B_{y z} |}_{j k l}^{n + \frac{1}{2}} = & {B_{x z} |}_{j k l}^{n} + \sqrt{\frac{μ_{0}}{ε_{j k l}}} {S_{x} |}_{j k l}^{n} \end{matrix}

2.2. Characteristics of the sampled volume and of the scattering medium

The considered volume is sampled with n_x × n_y × n_z points and sampling steps Δx = Δy = Δz = 0.3λ. Figure 1(a) shows a slice along the xz plane of the sampled volume. This volume is made up of several domains. Absorbing layers, with widths of $\frac{n_{x}}{8}$ , $\frac{n_{y}}{8}$ and $\frac{n_{z}}{8}$ pixels, are defined at the boundaries of the volume along each spatial dimension. The scattering medium is delimited by these absorbing layers along the x and y dimensions. Along the propagation axis (z axis), the input face, corresponding to the center of the three-cells source, and the output face, where measurements are made, are distant respectively of three and two pixels from the absorbing layers. Finally, in some cases the scattering medium will be surrounded by reflecting layers along the x and y directions; this point is discussed in section 3.1. The scattering medium is constituted by an homogeneous dielectric material (with a relative permittivity ε_m) randomly filled with dielectric spheres, with a radius r_s, a relative permittivity ε_s > ε_m and a volume concentration β_s. Figure 1(b) shows a typical distribution, in a (18λ)³ volume, of the dielectric spheres with a radius r_s = 0.9λ and a volume concentration β_s ≃ 5%. In this figure, we can observe that some dielectric spheres overlap and form aggregates with sizes greater than that of a sphere. Moreover, the spatial sampling step introduces some distortions in the spherical shape of the spheres. Although aggregates and the imperfect shape of the spheres introduce some small uncertainties, we do not observe significative variations in our numerical results when different realizations of the scattering medium with the same characteristics are performed. Moreover, numerical results are consistent with the estimation of the scattering coefficients calculated from the Mie theory [17] for spherical particles. In the presented results, the total volume is sampled with 128 × 128 × 256 points and the volume of the scattering medium is V_s ≃ (30λ)² × 56λ.

Fig. 1 (a) Section along the xz plane of the sampled volume. (b) Example of a scattering medium modelized by dielectric spheres randomly embedded in a homogeneous dielectric medium.

Download Full Size | PDF

3. Numerical results

First, we investigate the temporal properties of the scattered light exiting from inhomogeneous media designed with different values of β_s. In order to validate our numerical model, the temporal profile of the output light intensity is compared with the time-resolved transmittance function of a semi-infinite scattering medium given by Patterson et al. [18]:

\begin{matrix} T (t) = & {(4 π D ν)}^{- \frac{1}{2}} t^{- \frac{3}{2}} \times e^{- μ_{a} v t} \\ \times {(d - l_{s}^{*}) e^{- \frac{{(d - l_{s}^{*})}^{2}}{4 D v t}} - (d + l_{s}^{*}) e^{- \frac{{(d + l_{s}^{*})}^{2}}{4 D v t}} \\ + (3 d - l_{s}^{*}) e^{- \frac{{(3 d - l_{s}^{*})}^{2}}{4 D v t}} - (3 d + l_{s}^{*}) e^{- \frac{{(3 d + l_{s}^{*})}^{2}}{4 D v t}}} \end{matrix}

$D = \frac{1}{3 [μ_{a} + (1 - g) μ_{s}]}$ is the diffusion coefficient where μ_s, μ_a and g = 〈cos θ〉 are, respectively, the scattering coefficient, the linear absorption coefficient and the anisotropy coefficient. Because non-absorptive particles are considered, calculations are performed with μ_a = 0. $l_{s}^{*} = \frac{1}{(1 - g) μ_{s}}$ is the reduced mean free path, d is the thickness of the scattering medium, t is the time and v is the light velocity in the medium. We also compare the results of the PSTD algorithm with Monte Carlo simulations. In the Monte-Carlo method for photon transport, the propagation distance between two scattering events and the scattering direction cosines are randomly generated using the inverse distribution laws method [19].

Then, we analyze the time dependence of the polarization state of the scattered light for different polarization states of the incident light.

3.1. Time-resolved propagation of a short pulse in scattering media

A short pulse polarized along the x axis and with a 27 fs duration ( $\frac{1}{e^{2}}$ full width) is emitted by the source. The source has a very narrow circular gaussian shape in the transverse plane with a full width 6Δx = 1.8λ. For a radius r_s = 0.9λ and a refractive index of the non-absorptive dielectric spheres $n_{s} = \sqrt{ε_{s}} = 1.34$ , the size parameter of the spheres is $k r_{s} = \frac{2 π n_{s} r_{s}}{λ} = 7.6$ , which implies that light scattering occurs in the Mie regime, with an anisotropy coefficient g = 0.85. Although the scattering properties of a medium mostly depend on the relative size of the spheres with respect to the wavelength, we perform calculations by considering a wavelength λ = 1μm corresponding approximately to the wavelength used in our experiments reported in [1]. Considering that spheres are embedded in a lossless homogeneous medium with a refractive index $n_{m} = \sqrt{ε_{m}} = 1$ and with volume concentrations of 5%, 7%, 9% and 11%, the scattering coefficients are calculated using the Mie theory. Table 1 gives the values of these coefficients with respect to the volume concentrations β_s. μ_s and $μ_{s}^{*}$ are given in μm⁻¹. $N_{s} = μ_{s} d = \frac{d}{l_{s}}$ and $N_{s}^{*} = μ_{s}^{*} d = \frac{d}{l_{s}^{*}}$ are, respectively, the number of scattering events and the number of reduced scattering events calculated for a thickness d = 56λ = 56μm of the scattering medium.

Table 1. Scattering coefficients with respect to the volume concentrations β_s of the dielectric spheres.

View Table | View all tables in this article

Figures 2(a)–2(d) show, at different propagation times, single-frame excerpts from the video recording in the xz plane the intensity of the pulse propagating in a scattering medium with a volume concentration of 7% ( Media 1). It corresponds to a medium with a thickness twice as great as the reduced mean free path ( $d = 2 l_{s}^{*}$ ). The propagation axis is graduated in optical pathlength units : δ = 〈n_s〉z, where 〈n_s〉 = n_m +(n_s − n_m)β_s is the average value of the refractive index of the scattering medium. Typically, with β_s = 7%, 〈n_s〉 = 1.024 and the average speed of light in this medium is $v = \frac{c}{〈 n_{s} 〉} = 0.977 c$ . In Fig. 2(a), the white contours show the location and the shape of the scattering particles and in Figs. 2(b)–2(d) the white dotted lines represent the boundaries of the scattering medium. In Fig. 2(d), we can observe the increasing temporal and spatial spreading due to scattering of the pulse during its propagation. Light scattering produces a speckle pattern which is closely related to the location of the spheres in the medium. In the presented speckle patterns like in Fig. 2(b), some peak intensities much larger than the mean intensity make difficult viewing the temporal profile of the pulse during its propagation. Then, in order to resolve with a better accuracy the time spreading of the temporal shape of the pulse intensity, we plot in Fig. 2(e) the normalized pulse intensity profiles along the z axis, at the same times considered in Fig. 2(a)–2(d), where the pulse intensity is integrated in the (x, y) transverse plane. In the video we clearly see that the speckle patterns shown in Figs. 2(a)–2(d) are strongly correlated to the particular realization of the particles in the medium and we can observe the propagation of the backscattered light. On the other hand, we verified that the time shape of the pulse intensity, integrated in the (x, y) plane, is substantially the same whatever the realization.

Fig. 2 (a)–(d) Single-frame excerpts from the video ( Media 1) recording, in the xz plane and at the times 0, 70, 120 and 170 fs, the propagation of the 27 fs pulse in a scattering medium with β_s = 7%. (a) The white contours show the locations and the shapes of the particles in the considered xz plane of the medium. (b)–(d) The white dotted lines represent the boundaries of the scattering medium. (e) Corresponding profiles along the z axis of the pulse intensity integrated in the (x, y) transverse plane.

Download Full Size | PDF

In order to validate our numerical model, we compare systematically the temporal shape of the transmitted light given by our PSTD algorithm with the analytical expression of the time-resolved transmittance of a scattering medium given by Eq. (7) and with Monte-Carlo simulations [19]. The curves giving the results of the Monte Carlo simulations correspond to the histogram of the optical pathlengths calculated over 10000 trajectories. First, simulations are performed with scattering media only delimited by the absorbing layers. With the PTSD algorithm, the durations of the transmitted pulses are significantly shorter than the results given by the two other methods. We explain this discrepancy by the small transverse dimensions of the scattering mediium : indeed, a large part of the scattered light is lost in the absorbing layers. Therefore, in order to avoid these losses, we have added reflective layers at the boundaries of the transverse plane (see Fig. 1).

Figure 3 shows the different temporal shapes of the output light (green curves), obtained with the reflective layers, with respect to the thickness expressed in number of reduced mean free paths. Each curve is compared with the corresponding time-resolved transmittance curve (red curves) convolved with the input pulse (blue curves) and with the Monte-Carlo simulations (cyan curves). The x-axis is graduated in units of optical pathlength and the vertical dotted lines correspond to the optical pathlength corresponding to a straight propagation in the scattering medium : 〈n_s〉d. Each curve is normalized by its peak value. Though Monte-Carlo simulations are calculated over 10000 trajectories, some oscillations remain in the presented results, but these oscillations are progressively smoothed out when the number of trajectories increases. We can observe that the results of the PSTD algorithm fit with a very good agreement the Monte-Carlo simulations but not the transmittance function. This discrepancy can be explained by the fact that the diffusion model used to calculate the transmission function is valid for scattering media with a large number of reduced mean free paths, which is not the case here. In future works, scattering media with a larger numbers of $l_{s}^{*}$ should be investigated in order to demonstrate a possible better agreement between the PSTD and the diffusion theory.

Fig. 3 Comparison of the time-shapes of the transmitted pulse given by the PSTD algorithm (green curves), Monte-Carlo simulations (cyan curves) and analytical transmittance convolved by the input pulse shape (red curves) for different values of $N_{s}^{*}$ . The x-axis is graduated in optical pathlength units. The blue curves represent the input pulse and the vertical dotted lines correspond to the optical length of the scattering medium.

Download Full Size | PDF

In addition to the temporal analysis of our results, we used the method proposed in [20] to determine the diffusion properties of a scattering medium by measuring the speckle contrast resulting from the transmission of a femtosecond pulse. The contrast of a speckle is shown to be $C = \frac{1}{\sqrt{N}}$ , where N is the number of independent speckle patterns incoherently added. When a scattering medium is illuminated by a Fourier-transform short pulse, with a time duration τ_l, the transmitted pulse is stretched to approximately the average one-way traversal time $τ_{s} = \frac{d^{2}}{2 D v} > τ_{l}$ [21]. Then, the transmitted pulse is no more Fourier-transform and it can be considered as the superposition of N Fourier-transform pulses with a frequency bandwidth limited by $~ \frac{1}{τ_{s}}$ . It leads to the addition of independent speckle patterns whose contrast is related to this time by :

C ~ \sqrt{\frac{τ_{l}}{τ_{s}}}

From the time-integrated speckle patterns obtained with our PSTD algorithm, we have calculated the contrast $\frac{σ_{I}}{〈 I 〉}$ , where 〈I〉 is the mean intensity and σ_I is the standard deviation of the intensity fluctuations. Then, we have compared these contrast values with the contrast calculated with Eq. (8) as well as with the contrasts calculated from measurements of the time widths of the transmitted pulses. Figure 4 shows the typical output speckle pattern obtained with an input pulse, with a duration τ_l = 27 fs ( $\frac{1}{e^{2}}$ full width), transmitted through a medium with a thickness $d = 3 l_{s}^{*}$ . Table 2 summarizes the values of the contrasts calculated with the different methods. τ_s and τ′_l are, respectively, the average one-way traversal time (in fs) and the time full width of the transmitted pulse (in fs measured at $\frac{1}{e^{2}}$ . We can observe that all the methods give values of the contrast in good agreement even if the calculated one-way traversal times are significantly greater than the durations of the output pulses. The same measurements have been performed for increasing time durations of the incident pulse (from 27 to 253 fs). A scattering medium with a thickness $d = 2 l_{s}^{*}$ that corresponds to an average one-way traversal time τ_s = 577 fs is considered. As expected, the contrast of the speckle patterns observed at the output of the scattering medium increases from 0.24 to 0.52 with the time width of the input pulse and the contrasts calculated with the different methods give results in good agreement.

Fig. 4 Normalized intensity of the output speckle pattern obtained with a $3 l_{s}^{*}$ scattering medium. The contrast is C = 0.21.

Download Full Size | PDF

Table 2. Comparison of the contrasts of the speckle patterns calculated with different methods and for different values of $N_{s}^{*}$ .

View Table | View all tables in this article

The good fit between the results obtained with the PSTD and the Monte-Carlo simulations and the good agreement between the contrast values calculated by the different methods tend to prove that our PSTD algorithm correctly simulates the propagation of a light pulse through an inhomogeneous medium.

3.2. Time-resolved polarization analysis of scattered light

In this section, we study the space and time-resolved polarization state of the transmitted light with respect to the properties of the scattering medium and with respect to the polarization state of the input pulse. The time variation of the local polarization state of the transmitted light is characterized by calculating the components of the Stokes vector in the the output plane, at the time step nΔt and at the spatial sampling point (jΔx, kΔy) :

\begin{array}{l} {I |}_{j k}^{n} = & {E_{x} |}_{j k}^{n} {E_{x}^{*} |}_{j k}^{n} + {E_{y} |}_{j k}^{n} {E_{y}^{*} |}_{j k}^{n} \\ {Q |}_{j k}^{n} = & {E_{x} |}_{j k}^{n} {E_{x}^{*} |}_{j k}^{n} - {E_{y} |}_{j k}^{n} {E_{y}^{*} |}_{j k}^{n} \\ {U |}_{j k}^{n} = & {E_{x} |}_{j k}^{n} {E_{y}^{*} |}_{j k}^{n} + {E_{x}^{*} |}_{j k}^{n} {E_{y} |}_{j k}^{n} \\ {V |}_{j k}^{n} = & i ({E_{x} |}_{j k}^{n} {E_{y}^{*} |}_{j k}^{n} - {E_{x}^{*} |}_{j k}^{n} {E_{y} |}_{j k}^{n}) \end{array}

Figures 5(a)–5(d) show pictures of the components of the normalized local Stokes vector. These results are obtained with a $2 l_{s}^{*}$ scattering medium. 〈〉_t denotes the time integration of the considered variables. Correlations between the intensity speckle pattern (Fig. 5(a)) and the local value of the Q component (Fig. 5(b)) can be observed and show that the shinning speckle grains are related to the coherent part of the transmitted light. Figure 5(e) presents the corresponding degree of polarization defined as :

D O P = \frac{\sqrt{Q^{2} + U^{2} + V^{2}}}{I}

Fig. 5 Time-integrated components of the local Stokes vector and the local degree of polarization (DOP) of the light transmitted through a $2 l_{s}^{*}$ scattering medium. Stokes parameters are normalized by the peak intensity of the speckle pattern.

Download Full Size | PDF

For the different values of $N_{s}^{*}$ , we calculated the average values and the standard deviations of the Stokes parameters and of the DOP. Figure 6 and Table 3 summarize the statistical properties of these variables. 〈〉_xyt denotes the integration in time and space of the considered variables. As the input pulse is linearly polarized in the x direction, its normalized Stokes vector reads I_in = 1, Q_in = 1, U_in = 0 and V_in = 0 (DOP_in = 1). When $N_{s}^{*}$ increases, we can observe that the average values of the U and V components remain null with a large standard deviation and the average values of the Q parameter and of the DOP decrease. More surprisingly, we find that, although the rate of depolarized light increases with $N_{s}^{*}$ , a significant part of the transmitted light remains linearly polarized along the x axis even when the thickness of the scattering medium is larger than $l_{s}^{*}$ .

Fig. 6 Histograms of the time-integrated local Stokes parameters and of the DOP for different values of $N_{s}^{*}$ .

Download Full Size | PDF

Table 3. Average values and standard deviations of the space and time-integrated local Stokes parameter Q and DOP.

View Table | View all tables in this article

Now, let’s us examine the evolution in time of the state of polarization of the light transmitted through the scattering medium. We have plotted the variation with time of the instantaneous values of the spatially integrated Stokes parameters for the different scattering media (Fig. 7). The x-axis (i.e. time axis) is still graduated in optical pathlengths. The first observation concerns the peak of the Q parameter. When $N_{s}^{*}$ increases, its value decreases relatively to the peak intensity but its position is less delayed than the peak intensity position. The rising edges of the Q and I variations are confounded and the time width of the Q component (τ_Q ∼ 100 fs) is much larger than the width of the input pulse (τ_l = 27 fs). Moreover, while the time width of the intensity profiles increases with $N_{s}^{*}$ , the time width of the Q component remains constant as well as its temporal shape. We can also observe that the light is completely depolarized when the optical pathlength of the scattered light is greater than twice the optical thickness of the medium. These results confirm that the early transmitted light (corresponding to the so-called ballistic photons) is preferentially vertically polarized as well as the weakly scattered light corresponding to the so-called snake photons. Moreover, it shows that the time arrival of the weakly scattered light is not determined by the scattering coefficient of the medium. Future works should be done in order to characterize more precisely this property with respect to the thickness of the scattering medium, the size and the shape of the scattering particles.

Fig. 7 Time variation of the space-integrated components of the Stokes vector for the different media. x-axes are graduated in optical pathlength units.

Download Full Size | PDF

We studied also the variation of the DOP during time. Figure 8(a) represents the transient regime of the space and time-integrated 〈DOP〉_xyt obtained with the different scattering media. The curves exhibit maxima very close to 1 at the early times and decrease up to the final values given in Table 3. Figures 8(b)–8(d) show the local time-integrated DOP at different times for the $2 l_{s}^{*}$ scattering medium. In Fig. 8(b), that corresponds to the early time of the scattered light, the local DOP is close to 1 in an area located at the center of the transverse plane and with transverse dimensions close to the size of the source. This corresponds to the light that travels in a straight line across the medium. At the time corresponding to the peak value (Fig. 8(c)), the DOP is close to 1 in the whole transverse section. For increasing times, the local DOP exhibits some increasing fluctuations until it reaches its final state depicted by the Fig. 8(d).

Fig. 8 (a) Time variation of the time and space integrated DOP for the different media. (b)–(d) Time integrated local DOP at different time for the $2 l_{s}^{*}$ medium. The vertical black dotted line represents the optical pathlength through the homogeneous medium.

Download Full Size | PDF

The last numerical experiments concern the influence of the polarization state of light on the scattering process. Linear and circular polarizations of the light emitted by the source are considered. In both cases, the time-resolved variation of the Stokes components and of the DOP of the transmitted light are studied for different scattering media with increasing values of β_s. We have modelized scattering media made of dielectric spheres, with a radius r_s = 0.6λ and a refractive index n_s = 1.34, embedded in a medium with a refractive index n_m = 1. At λ = 1μm, the anisotropy coefficient is g = 0.81 and the size parameter of a sphere is $k r_{s} = \frac{2 π n_{s} r_{s}}{λ} = 5.0$ . With the volume concentrations of 7, 9 10 and 12%, the scattering media are characterized by the coefficients $μ_{s}^{*} = 0.043$ , 0.053, 0.063 and 0.073 μm⁻¹, and the thickness of the scattering media corresponds, respectively, to a number of reduced scattering events $N_{s}^{*} = 2.4$ , 3, 3.6 and 4.1.

In these numerical experiments, we removed the reflective layers at the boundaries of the scattering medium in order to prevent the right ↔ left changes of a circular polarization due to the reflections. Figures 9(a) and 9(b) show, for the two values of β_s, the time-resolved variations of the space-integrated Stokes components (〈I〉_xy, 〈Q〉_xy) and (〈I〉_xy, 〈V〉_xy) of the transmitted light initially linearly and circularly polarized. These components are normalized by the peak values of the intensity. In both figures, we can observe that the time variations of the intensities are confounded. It shows that the intensity temporal shape of the transmitted light is independent of the incident polarization state [2]. However, we can observe on the same figures, that the component Q of the linearly polarized light decreases a little faster than the component V of the circularly polarized light. This result is consistent with previous works showing that, in the Mie regime, linearly polarized light is more rapidly depolarized than circularly polarized light [2, 22]. The exponential time decay of the Stokes components Q and V is also clearly exhibited, with a decay constant of about 0.25μm⁻¹ with the time axis graduated in optical pathlength units (corresponding to a time decay constant of about 0.075 fs⁻¹). With the x-axis graduated in number of scattering events unit (c(t − t₀) + 〈n_s〉d) $μ_{s}^{*}$ , we can observe that the light is depolarized if it has experienced an optical path longer than an optical thickness of the medium of approximately one reduced mean free path.

Fig. 9 Comparison of the time variations of the space integrated Stokes parameters (I, Q) or (I, V), respectively for linearly (blue lines) or circularly (red lines) polarized incident pulses for scattering media with volume concentrations (a) β_s = 7% and (b) β_s = 9%. (c) Time variation of the space integrated linear and circular degrees of polarization when the incident light is, respectively, linearly and circularly polarized with β_s =7, 9, 10 and 12 %.

Download Full Size | PDF

Because the mean values of the Stokes components (U, V) and (Q, U) are null, respectively for the linear and the circular polarization state of the incident light, we calculated the corresponding instantaneous time variation of the space integrated degree of linear polarization ${〈 {DOP}_{L} 〉}_{x y} = {〈 \frac{Q}{I} 〉}_{x y}$ and the degree of circular polarization ${〈 {DOP}_{C} 〉}_{x y} = {〈 \frac{V}{I} 〉}_{x y}$ . Figure 9(c) depicts the time variation of the DOP_L and the DOP_C for scattering media with β_s =7, 9, 10 and 12 %. The x-axis origin of the curves corresponds to the optical thickness of the media. In agreement with the Q and V variations,the DOP_L decreases a little faster than the DOP_C. Since the FWHM of the curves decreases approximately from 19 μm (i.e. 63 fs) to 16 μm (i.e. 53 fs) when β_s increases, the FWHM of the DOP curves expressed in optical pathlength unit divided by the reduced mean free paths increases from 0.8 to 1.1, if we take into account the reduced mean free paths calculated for the different values of β_s. It confirms that the exiting scattered light has undergone, on average, one reduced scattering event before being depolarized. In Fig. 9(b) weak oscillations (they would be almost invisible with a linear scale) can be observed in the time variation of the Stokes parameters V for an incident wave circularly polarized (red dashed curve) between optical pathlengths of 80 and 90 μm. These oscillations, traducing a fluctuation of the instantaneous polarization state of light, are reproducible but with different shapes for different realizations of the scattering medium. Similar oscillations can be observed in Fig. 7 for the Stokes parameters U and V.

4. Conclusion

For the first time to our knowledge, a 3D-PSTD algorithm has been implemented to solve the full-wave Maxwell’s equations in order to resolve, in time and space, the propagation of electromagnetic waves in inhomogeneous media with realistic dimensions in the three dimensions of space. The polarization state of the light exiting the scattering media was also analyzed with respect to the characteristics of the scattering medium and to the polarization state of the incident light. The numerical results match with a very good agreement the Monte-Carlo simulations and are consistent with previous works. We show that, with this numerical model, it is possible to resolve and fully characterize all the properties of an electromagnetic wave propagating in a complex medium. In addition, we confirm that the ballistic part and the weakly scattered light emerging from a scattering medium maintain a polarization state very close to the polarization state of the incident light even when the thickness of the scattering medium is much greater than one reduced mean free path. This explains in part our surprising experimental results reported in [1], where the object image is retrieved by the phase conjugated light despite the rotation of the polarization due to the nonlinear process of phase-conjugation. These results encourage us to model other problems related to phenomena of light scattering as the polarization dependence property of the backscattered light [23] or the propagation of light in inhomogeneous media with birefringence or optical activity properties. In future works, we will model time reversal of scattered light by phase conjugation, obtained by three or four wave mixing, in order to more fully understand and explain our experimental results or works performed by other groups [24].

Acknowledgments

This work has been supported by the Agence Nationale de la Recherche (ICLM, project ANR-2011-BS04-017-03) and the Conseil Régional de Franche-Comté. Computations have been performed on the supercomputer facilities of the Mésocentre de calcul de Franche–Comté.

References and links

1. F. Devaux and E. Lantz, “Real time suppression of turbidity of biological tissues in motion by three-wave mixing phase conjugation,” J. of Biomed. Opt. 18, 111405 (2013). [CrossRef]

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

3. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Transactions on Antennas and Propagation 14, 302–307 (1966). [CrossRef]

4. Q. H. Liu, “The PSTD algorithm: a time-domain method requiring only two cells per wavelength,” Microw. Opt. Technol. Lett. 15, 158–165 (1997). [CrossRef]

5. Z. Tang and Q. H. Liu, “The 2.5D FDTD and Fourier PSTD methods and applications,” Microw. Opt. Technol. Lett. 36, 430–436 (2003). [CrossRef]

6. T. W. Lee and S. C. Hagness, “Pseudospectral time-domain methods for modeling optical wave propagation in second-order nonlinear materials,” J. Opt. Soc. Am. B 21, 330–342 (2004). [CrossRef]

7. X. Liu and Y. Chen, “Applications of transformed-space non-uniform PSTD (TSNU-PSTD) in scattering analysis with the use of the non-uniform FFT,” Microw. Opt. Technol. Lett. 38, 16–21 (2003). [CrossRef]

8. C. Liu, R. L. Panetta, and P. Yang, “Application of the pseudo-spectral time domain method to compute particle single-scattering properties for size parameters up to 200,” J. Quant. Spectrosc. Transfer 113, 1728–1740 (2012). [CrossRef]

9. S. H. Tseng, A. Taflove, D. Maitland, and V. Backman, “Pseudospectral time simulations of mutiple light scattering in three-dimensional macroscopic random media,” Radio Science 41, RS4009 (2006). [CrossRef]

10. C. Liu, L. Bi, R. L. Panetta, P. Yang, and M. A. Yurkin, “Comparison between the pseudospectral time domain method and the discrete dipole approximation for light scattering simulations,” Opt. Express 20, 16763–16776 (2012). [CrossRef]

11. S. H. Tseng and C. Yang, “2-D PSTD simulation of optical phase conjugation for turbidity suppression,” Opt. Express 15, 1605–1616 (2007). [CrossRef]

12. S. H. Tseng, “PSTD simulation of optical phase conjugation of light propagating long optical paths,” Opt. Express 17, 5490–5495 (2009). [CrossRef] [PubMed]

13. Q. H. liu and G. Zhao, “Review of PSTD methods for transient electromagnetics,” Int. J. Numer. Model. 17, 299–323 (2004). [CrossRef]

14. J. -P. Berenger, “A perfect matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1994). [CrossRef]

15. W. C. Chew and W. H. Weedon, “A 3D perfectly matched medium from modified Maxwell’s equations,” Microw. Opt. Technol. Lett. 7, 599–604 (1994). [CrossRef]

16. Z. Li, “The optimal spatially-smoothed source patterns for the pseudospectral time-domain method,” IEEE Transactions on Antennas and Propagation 58, 227–229 (2010). [CrossRef]

17. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley & Sons, 1983 Chap. 4, pp. 82–129).

18. M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. 28, 2331–2336 (1989). [CrossRef] [PubMed]

19. L. H. Wang, S. L. Jacques, and L. Q. Zheng, “Monte Carlo modeling of photon transport in multi-layered tissues,” Computer Methods and Programs in Biomedicine 47, 131–146 (1995). [CrossRef]

20. N. Curry, P. Bondareff, M. Leclerq, N. K. Van Hulst, R. Sapienza, S. Gigan, and S. Grésillon, “Direct determination of diffusion properties of random media from speckle contrast,” Opt. Lett. 36, 3332–3334 (2011). [CrossRef] [PubMed]

21. R. Landauer and M. Büttiker, “Diffusive traversal time: Effective area in magnetically induced interference,” Phys. Rev. B 36, 6255–6210 (1987). [CrossRef]

22. D. Bicout, C. Brosseau, A. S. martinez, and J. M. Schmitt, “Depolarization of multiply scattered waves by spherical difFusers: Influence of the size parameter,” Phys. Rev. E 49, 1767–1770 (1994). [CrossRef]

23. F. C. MacKintosh, J. X. Zhu, D. J. Pine, and D. A. Weitz, “Polarization memory of multiply scattered light,” Phys. Rev. B 40, 9342–9345 (1989). [CrossRef]

24. M. Cui, E. J. McDowell, and C. Yang, “Observation of polarization-gate based reconstruction quality improvement during the process of turbidity suppression by optical phase conjugation,” Appl. Phys. Lett. 95, 123702 (2009). [CrossRef] [PubMed]

β_s (%)	μ_s (μm⁻¹)	$μ_{s}^{*}$ (μm⁻¹)	N_s	$N_{s}^{*}$

5	0.18	0.027	10.1	1.5
7	0.24	0.036	13.4	2.0
9	0.30	0.045	16.7	2.5
11	0.35	0.053	19.7	3

$N_{s}^{*}$	$\frac{σ_{I}}{〈 I 〉}$	τ_s (fs)	$\sqrt{\frac{τ_{l}}{τ_{s}}}$	τ′_l (fs)	$\sqrt{\frac{τ_{l}}{τ_{l}^{'}}}$

1.5	0.29	433	0.25	311	0.30
2	0.23	577	0.22	449	0.24
2.5	0.22	724	0.19	512	0.23
3	0.21	863	0.18	586	0.21

$N_{s}^{*}$	${\bar{〈 Q 〉}}_{x y t}$	σ_Q	${\bar{〈 D O P 〉}}_{x y t}$	σ_DOP

1.5	0.44	0.17	0.49	0.14
2	0.31	0.17	0.40	0.13
2.5	0.25	0.17	0.35	0.13
3	0.20	0.17	0.31	0.12

β_s (%)	μ_s (μm⁻¹)	$μ_{s}^{*}$ (μm⁻¹)	N_s	$N_{s}^{*}$

5	0.18	0.027	10.1	1.5
7	0.24	0.036	13.4	2.0
9	0.30	0.045	16.7	2.5
11	0.35	0.053	19.7	3

$N_{s}^{*}$	$\frac{σ_{I}}{〈 I 〉}$	τ_s (fs)	$\sqrt{\frac{τ_{l}}{τ_{s}}}$	τ′_l (fs)	$\sqrt{\frac{τ_{l}}{τ_{l}^{'}}}$

1.5	0.29	433	0.25	311	0.30
2	0.23	577	0.22	449	0.24
2.5	0.22	724	0.19	512	0.23
3	0.21	863	0.18	586	0.21

3D-PSTD simulation and polarization analysis of a light pulse transmitted through a scattering medium

Abstract

1. Introduction

2. Principle of the tri-dimensional pseudo-spectral time domain method

2.1. Tri-dimensional pseudo-spectral time domain scheme

2.2. Characteristics of the sampled volume and of the scattering medium

3. Numerical results

3.1. Time-resolved propagation of a short pulse in scattering media

3.2. Time-resolved polarization analysis of scattered light

4. Conclusion

Acknowledgments

References and links

Supplementary Material (1)

Cited By

Figures (9)

Tables (3)

Equations (10)

Optics Express