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

A tridimensional pseudo-spectral time domain (3D-PSTD) algorithm, that solves the full-wave Maxwell's equations by using Fourier transforms to calculate the spatial derivatives, has been applied to determine the time characteristics of the propagation of electromagnetic waves in inhomogeneous media. Since the 3D simulation gives access to the full-vector components of the electromagnetic fields, it allowed us to analyse the polarization state of the scattered light with respect to the characteristics of the scattering medium and the polarization state of the incident light. 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 fact is consistent with our recently reported experimental results, where a rotation of the polarization does not preclude the restoration of an image by phase conjugation.


3D-PSTD simulation and polarization
analysis of a light pulse transmitted through a scattering medium Fabrice Devaux* and Eric Lantz

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 gridsampling 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][5][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].

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

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 : (1) 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 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 : where ν y , F y and F −1 y 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: 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 (µ 0 = 1, ε 0 = 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 = ∆x 16c = ∆x 16 , which is smaller than the maximum time step 2∆x πc √ 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 + 1 2 , to the terms D xz | n jkl and D yz | n jkl using the relations : S x | n jkl and S y | n jkl are the transverse amplitude components of the pulsed source term (delayed with a time t 0 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: ψ, ϕ x and ϕ y are the parameters used to define the polarization state of the electromagnetic wave emitted by the source. The spatial amplitude S 0 | jkl is designed with a Gaussian shape in the (x, y) transverse plane and with the optimized three-cells normalized pattern [ 1 4 , 1 2 , 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 It is obtained as follows: 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 n x 8 , n y 8 and 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λ ) 3 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λ ) 2 × 56λ .

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]: 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 = 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.

Time-resolved propagation of a short pulse in scattering media
A short pulse polarized along the x axis and with a 27 f s duration ( 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 = √ ε s = 1.34, the size parameter of the spheres is kr s = 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 = √ ε 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 −1 .
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.    . It corresponds to a medium with a thickness twice as great as the reduced mean free path (d = 2l * 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 = c n s = 0.977c. 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.
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.
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 = 1 √ 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 = d 2 2Dv > τ 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 ∼ 1 τ s . It leads to the addition of independent speckle patterns whose contrast is related to this time by : From the time-integrated speckle patterns obtained with our PSTD algorithm, we have calculated the contrast σ 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 f s ( 1 e 2 full width), transmitted through a medium with a thickness d = 3l * 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 f s) and the time full width of the transmitted pulse (in f s measured at 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 f s). A scattering medium with a thickness d = 2l * s that corresponds to an average one-way traversal time τ s = 577 f s 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  Fig. 4. Normalized intensity of the output speckle pattern obtained with a 3l * s scattering medium. The contrast is C = 0. 21. pulse and the contrasts calculated with the different methods give results in good agreement.  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.

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) : Figures 5(a)-5(d) show pictures of the components of the normalized local Stokes vector. These results are obtained with a 2l * 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 :  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 . 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 f s) is much larger than the width of the input pulse (τ l = 27 f s). 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 socalled 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.
We studied also the variation of the DOP during time. Figure 8 given in Table 3. Figures 8(b)-8(d) show the local time-integrated DOP at different times for the 2l * 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).
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 kr s = 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 −1 , 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 −1 with the time axis graduated in optical pathlength units ( corresponding to a time decay constant of about 0.075 f s −1 ). With the x-axis graduated in number of scattering events unit (c(t − t 0 ) + 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.
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 xy = Q I xy and the degree of circular polarization DOP C xy = V I xy . 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 f s) to 16 µm (i.e. 53 f s) 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 .

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 or-der to more fully understand and explain our experimental results or works performed by other groups [24].