5862
OPTICS LETTERS / Vol. 39, No. 20 / October 15, 2014
Early-time diffusion in pulse propagation through dilute random media Elizabeth Bleszynski,* Marek Bleszynski, and Thomas Jaroszewicz Monopole Research, Thousand Oaks, California 91360, USA *Corresponding author: [email protected] Received July 7, 2014; accepted August 26, 2014; posted September 5, 2014 (Doc. ID 215211); published October 9, 2014 Propagation of short infrared/optical pulses in dilute random media (e.g., atmospheric clouds, fog, dust, or aerosols) consisting of large, compared to the wavelength, scatterers is analyzed. A rigorous approach based on analytic complex-contour integration of numerically determined cut and pole singularities of the radiative transport equation solution in the Fourier space is presented. It is found that the intensity of a propagating pulse, in addition to the coherent (“ballistic”) contribution and a long late-time diffusive tail, also exhibits a sharply rising early-time component that (i) can be attributed to the small-angle diffractive part of the scattering cross-section on medium particles, (ii) is attenuated proportionally to the nondiffractive rather than total cross-section, and (iii) can be extracted by high-pass filtering of the received pulse. © 2014 Optical Society of America OCIS codes: (030.5620) Radiative transfer; (290.1090) Aerosol and cloud effects; (290.1990) Diffusion; (290.2558) Forward scattering; (290.4210) Multiple scattering; (070.2615) Frequency filtering. http://dx.doi.org/10.1364/OL.39.005862
Imaging and communication with optical or infrared pulsed signals through obscuring random media (e.g., atmospheric clouds, fog, dust, or aerosols) is a longstanding and challenging problem, both experimentally and theoretically. From the experimental and applications perspective, a coherently detected pulse field preserves its time profile, but is strongly attenuated, at the rate proportional to the total cross-section of the wave on an individual medium scatterer. The incoherently detected field intensity [or, more generally, the mutual coherence function (MCF)], although attenuated at a lower rate (proportional to the absorption cross-section), develops a long diffusive tail that causes loss of resolution in imaging and loss of bandwidth in communication. From the theoretical point of view, pulse propagation in a statistically homogeneous, dilute, random medium can be adequately described by the radiative transport equation (RTE) [1–3]. However, obtaining its analytic solution for the time-dependent intensity (except for the case of isotropic scatterers [4,5]) proved to be a challenging problem. Consequently, most of the analysis and understanding of these transient phenomena have been based on Monte Carlo (MC) simulations (e.g., [6–9]) and some approximate methods (e.g., [10–12]). Only recently analytic solutions have been derived for arbitrary anisotropic scatterers [13–15] and applied in the context of biological media. This Letter presents first results of our analysis of propagation of short optical pulses in an atmospheric cloud-type medium. In this regime, the scatterers are relatively large compared to the wavelength (a typical average droplet radius [16] being a 5 μm ≈ 7.9λ0 for λ0 0.633 μm), the absorption is negligible, and the differential scattering cross-section on a medium constituent clearly exhibits a diffractive peak of angular width Δθ ∼ λ0 ∕2πa ≡ 1∕k0 a. We find, by means of the analytic solution technique that we briefly describe in this Letter, that the evolved pulse, in addition to the coherent (ballistic) contribution and a long diffusive tail, develops a distinctive 0146-9592/14/205862-04$15.00/0
early-time component characterized by (i) a rise time much shorter than the usual diffusion time scale, and (ii) attenuation significantly lower than that of the coherent intensity. Further, we show that, because of property (i), the early-time pulse component can be isolated by means of a “frequency-gating” technique implementable as high-pass filtering. In contrast to time-gating methods used widely for isolating an early-time part of the received pulse [17–19], the proposed gating technique does not require information on the expected arrival time of the signal and is, therefore, well suited in communication and in range imaging through a scattering medium. The early-time pulse structure can be interpreted as due to a diffusion-type process driven by a twodimensional angular random walk with steps of the order of the diffractive peak width, Δθ. Time delay in a single step is proportional to lt ∕cθ2 , where lt is the coherent attenuation length and θ the angle between the path segment and the observation direction. In N ∼ R∕lt steps needed to reach the distance R, the total delay becomes Δt ∼ Nlt ∕chθ2 i and, because of the angular random walk, hθ2 i ∼ NΔθ2 . Therefore, Δt ∼ Δθ2 R2 ∕clt , which is a relation characteristic of diffusion, but with a small coefficient Δθ2 ∼ 1∕k0 a2 . This justifies referring to the early-time pulse component as “early-time diffusion,” distinct from the usual latetime diffusion associated with the large-angle isotropic scattering. The above interpretation indicates that the appearance of the early-time diffusion component is due to a distinctive diffractive peak characteristic of scatterers large in comparison to the wavelength. Early time enhancement of evolved pulses has been noticed before in MC simulations [7] for very large (a ∼ 50 μm) strongly diffractive scatterers (Fig. 1) typical of sprays in combustion processes [20]. However, small length and time scales characterizing dense sprays (lt ∼ 1 cm, lt ∕c ∼ 3.3 ps) [21] make it difficult to experimentally separate the expected early-time diffusion © 2014 Optical Society of America
October 15, 2014 / Vol. 39, No. 20 / OPTICS LETTERS
5863
~ with the two-field Green function GΩ; P; sˆ μt − iΩ∕v0 − P · sˆ−1 . Here v0 ≈ c is the speed of coherent wave in the effective medium, μt 1∕lt is the attenuation coefficient, and the scattering function is given by Σˆs · sˆ 0 n0 σ s ˆs · sˆ 0 . The RTE (2) can be solved by projectingR the specific intensity onto partial-waves, I~ l Ω; P≔ d2 sˆ∕ ~ 4πpl Pˆ · sˆIΩ; P; sˆ , with normalized Legendre polynop mials pl x≔ 2l 1Pl x. Truncating the expansion at l < Lmax ∼ k0 a yields the matrix equation ~ ~ IΩ; P − GΩ; PΣIΩ; P GΩ; P0 ; Fig. 1. Small-angle Mie and Henyey–Greenstein (H–G) differential cross-sections at λ0 0.633 R 1 μm, normalized to total scatd cos θ σ s cos θ, for various tering cross-section σ s 2π −1 scatterer sizes.
signal from other scattering mechanisms. Observing early-time diffusion in biomedical imaging [22] is also difficult due to smaller scatterers [23], and thus insignificant diffraction [see Fig. 1, Mie scattering for red blood cells (a 2.78 μm) and Henyey–Greenstein [24] model for most biological tissues] combined with even higher scatterer concentration and smaller length and time scales. In our problem, on the other hand, droplets in an atmospheric cloud are large enough (a 5 μm) to give rise to considerable diffraction and the scatterer number density is much lower, typically n0 109 m−3 , resulting in the coherent attenuation length of lt 1∕n0 σ t ≈ 6 m. The corresponding attenuation time scale lt ∕c ≈ 20 ns is orders-of-magnitude longer; hence resolving the early-time signal should prove feasible [Eq. (7) and Fig. 4]. In the analysis described here, we adopt the customary approximate treatment of electromagnetic field components as scalar fields. We assume a coherent point-like isotropic source located at the origin in infinite space, emitting an intensity pulse St of duration T p , such that its physical length cT p is large compared to both the wavelength λ0 and the scatterer size a. Those two conditions are necessary in order to be able to reduce the Bethe–Salpeter equation (BSE) to its approximate RTE form: the first to simplify the BSE Green function and the second to neglect the dependence of the scattering amplitudes on the relative frequency (denoted below by Ω). The evolved field distribution is then described by the specific intensity (radiance) I S t; R; sˆ [1,2], where the unit vector sˆ ; is the energy flux direction at the time t and at the point R. The Fourier-space radiance I~ S Ω; P; sˆ≔
Z
dt d3 R eiΩt e−iP·R I S t; R; sˆ
(1)
~ ~ is given by I~ S Ω; P; sˆ SΩ IΩ; P; sˆ, where S~ is the Fourier transform of the source and I~ satisfies the RTE for an impulse (∼δt) source, ~ ~ IΩ; P; sˆ − GΩ; P; sˆ
Z
~ ~ d2 sˆ 0 Σˆs · sˆ 0 IΩ; P; sˆ 0 GΩ; P; sˆ (2)
(3)
where the vector I~ consists of Lmax compo~ nents R 1 I l , components of 2the vector Σ are Σl 2π −1 dx Pl x Σx, and the Lmax elements of the matrix G are R
Gl;l0 Ω; P
d2 sˆ ∕4πpl Pˆ · sˆ pl0 Pˆ · sˆGΩ; P; sˆ;
4
the r.h.s. of Eq. (3) is the l R0 column of G. The field intensity, It; R d2 sˆ ∕4πIt; R; sˆ, is now spherically symmetric and is given by Z It; R
∞
0
Z dP
∞
−∞
dΩ
P sinPR −iΩt ~ e I 0 Ω; P; 4π 3 R
(5)
where I~ 0 Ω; P I − GΩ; PΣ−1 GΩ; P0;0 is the l 0 partial wave of the solution of Eq. (3). While computing the Fourier-space solution I~ 0 Ω; P is straightforward, evaluation of the Fourier integral (5) (as discussed, e.g., in [14]) is nontrivial, especially for large scatterers and in the “early time,” near-light-cone region (v0 t ≳ R ≫ lt and v0 t − R ≪ lt ). To overcome those difficulties, we devised an approach in which we locate (partly analytically and partly numerically) singularities of I~ 0 Ω; P, perform contour integration over Ω, and integrate the result over real P [25]. The procedure relies on the fact that the function I~ 0 Ω; P has two types of singularities [26]: (i) Branch points at Ω v0 P − iμt and the associated cut in GΩ; P. The resultant cut contribution to the integral (5) is straightforward to calculate, and it decays with R at the same rate as the coherent intensity component (i.e., it becomes negligible for R ≫ lt ). (ii) Poles Ωj P, due to the matrix I − GΩ; PΣ not being invertible. Since the determinant of the matrix I − GΣ is holomorphic on the cut complex Ω plane, it may only have isolated zeros, resulting in a meromorphic function I~ 0 Ω; P, having only isolated poles. As P increases, the poles travel across the complex Ω plane until they reach the cut and disappear, one by one, from the “physical” Riemann sheet. The pole contribution to the intensity for a given P is, therefore, the sum of terms 1 It;Rj 2 2π R
Z 0
Pj
dPP sinPReImΩj Pt Imfαj Pe−iReΩj Pt g; (6)
where Ωj P is the jth pole trajectory, αj P is the pole residue, and P j is the value above which the pole
5864
OPTICS LETTERS / Vol. 39, No. 20 / October 15, 2014
disappears. For jPj ≥ maxfP j g ∼ k0 a2 μt , there are no more poles left on the physical sheet. This behavior of singularities can be understood by observing that, for P ≫ μt and k0 a2 ≫ 1, Eq. (2) becomes equivalent to the Schrödinger equation with an imaginary potential proportional to P −1 and the Fourier transform of Σ in the transverse components of sˆ . In the large P (thus weak coupling) regime, the discrete modes disappear and only the continuous spectrum remains. In Fig. 2 we plot, for propagation distances R1 12lt and R2 16lt , the dimensionless intensity J S v0 t∕lt ; tR∕lt l3t ∕v0 I S t; R, expressed as a function of dimensionless variables v0 t∕lt and R∕lt , and computed as the sum of cut and pole contributions. We traced the pole locations, and evaluated their residues, by means of eigenvalue decomposition of the matrix I − GΣ. We assumed a coherent source emitting a Gaussianmodulated field of carrier frequency c∕λ0 , with λ0 0.633 μm, resulting in intensity St exp−t2 ∕2T 2p ∕ p 2π T p , with T p 60 ps. The time-dependent intensity exhibits a distinctive early-time diffusion behavior, characterized by a narrow spike followed by a much broader late-diffusion maximum and a long diffusive tail. As the propagation distance increases, the broad diffusion shoulder starts overlapping the trailing edge of the early-time diffusion peak, but its sharply rising leading edge remains. Theoretical estimates and a number of computations for various scatterer sizes indicate that the rise time of the early-diffusion structure is approximately given by ΔtED ≈
κ ED R2 ; k0 a2 v0 lt
(7)
with a numerical factor κ ED ≲ 0.1. In the considered example, ΔtED ≈ 100 ps for R 12lt and ΔtED ≈ 180 ps for R 16lt ; these time scales are significantly smaller than the late-diffusion times (tens of nanoseconds). We note that the chosen pulse length T p < ΔtED allows us to take full advantage of the achievable range resolution cΔtED . Similar estimates indicate that the early-time diffusion intensity is attenuated, approximately, at the rate μED n0 σ t − σ d ;
proportional to the nondiffractive cross-section, rather than the total one (the diffractive cross-section σ d is the differential cross-section integrated over the diffractive forward-scattering cone). In our example the attenuation is reduced relative to the coherent attenuation by σ t − σ d ∕σ t ≈ 0.65. Both features, (7) and (8), can be explained in terms of the evolution of the poles appearing in Eq. (6). The two leading poles in our example problem are shown in Fig. 3. The j 0 pole emerges from the cut at P 0 ≈ 1674μt and, as long as P ≳ 10μt , the propagation velocity Re Ω0 P∕P remains very close to v0 , while the attenuation coefficient −Im Ω0 P∕v0 decreases from μt to about 60% of that value, resulting in an early, narrow, and relatively weakly attenuated contribution to the intensity. The j 1 pole, which emerges from the cut at a much smaller P value and moves up to Ω 0, is responsible for most of the late-time diffusion. The origin of the early-time diffusion component suggests that it should be possible to isolate it by selecting the region of large P values in the Fourier-space intensity; this amounts to high-pass filtering of the received timedomain intensity. In our example, we used a Gaussian ~ filter ΦΩ 1 − exp−T 2f Ω2 ∕2 of width v0 T f 0.015lt 9 cm, chosen such that v0 ΔtED < v0 T f ≪ lt . The original and filtered pulses, J S and Φ ∘ J S are plotted in Fig. 4 for distances R1 12lt and R2 16lt . The R2 intensity was scaled by the factor expμED R2 − R1 R2 ∕R1 2 ≈ 24. The close similarity between the R1 and the scaled R2 filtered intensities validates the reduced attenuation prediction of Eq. (8). The time shift in their maxima confirms the estimate of Eq. (7). The plots show that, already for the propagation distance of 12lt , the magnitude of the filtered early diffusive signal exceeds that of the coherent intensity by the factor ∼100. It is also evident that frequency gating does suppress the smooth late-time diffusion contribution while preserving the sharp leading edge of the early-time diffusion signal. To summarize, we expounded main properties (small rise time and reduced attenuation) of the early-time
(8)
Fig. 2. Intensity distributions showing early- and late-time diffusion behavior, for propagation distances R1 12lt and R2 16lt , with the time scale indicated on top. The curves are visually indistinguishable from the MC simulation results.
Fig. 3. Evolution, on the complex Ω plane, of two leading poles Ωj P, j 0 and j 1, of the integrand I~ 0 Ω; P of the Fourier transform (5). Each pole is either located on the imaginary axis or is accompanied by its mirror reflection in that axis. The thick horizontal line represents the cut and the arrows mark positions of poles at several values of P.
October 15, 2014 / Vol. 39, No. 20 / OPTICS LETTERS
Fig. 4. Received (top) and filtered (bottom) intensities for R1 12lt and R2 16lt . Coherent contributions, multiplied by factors 102 and 104 , are shown for reference.
diffusion component in the evolved short pulse intensity, and its relation to small-angle diffractive scattering; we also devised a filtering technique for extracting it from the received intensity. Our analysis indicates that, for cloud-type media, the expected duration of the early diffusion signal (∼100 ps) falls within resolution capabilities of current opto-electronic devices, and, at the same time, provides range resolution of a few centimeters. Our results may have important implications in highresolution range imaging as well as communication through obscuring (atmospheric clouds, fog, dust, or aerosols) media. Further developments may include (i) generalizations to other sources and detectors (limited detector field-of-view, collimated beams, source and detector arrays, partially coherent transmitted fields), (ii) generalization to transmission through a slab geometry, and (iii) alternative signal processing techniques in lowintensity measurements, such as single-photon counting [27] or first-photon-imaging [28] methods. This research is supported by the U. S. Air Force Office of Scientific Research under grant no. FA 9550-12-1-0121.
5865
References and Notes 1. S. Chandrasekhar, Radiative Transfer (Dover, 1960). 2. A. Ishimaru, Wave Propagation and Scattering in Random Media (Wiley, 1999). 3. L. Ryzhik, G. Papanicolaou, and J. B. Keller, Wave Motion 24, 327 (1996). 4. K. M. Case, Ann. Phys. 9, 1 (1960). 5. J. Palmeri, J. Stat. Phys. 58, 885 (1990). 6. E. A. Bucher, Appl. Opt. 12, 2391 (1973). 7. C. Rozé, T. Girasole, L. Méès, G. Gréhan, L. Hespel, and A. Delfour, Opt. Commun. 220, 237 (2003). 8. C. Calba, L. Méès, C. Rozé, and T. Girasole, J. Opt. Soc. Am. A 25, 1541 (2008). 9. B. Wu, Z. Hajjarian, and M. Kavehrad, Appl. Opt. 47, 3168 (2008). 10. A. Ishimaru and S. T. Hong, Radio Sci. 10, 637 (1975). 11. S. T. Hong and A. Ishimaru, Radio Science 11, 551 (1976). 12. W. Cai, M. Lax, and R. R. Alfano, Phys. Rev. E 63, 016606 (2000). 13. A. Liemert and A. Kienle, Phys. Rev. A 83, 015804 (2011). 14. A. Liemert and A. Kienle, Biomed. Opt. Express 3, 543 (2012). 15. A. Liemert and A. Kienle, Phys. Rev. E 86, 036603 (2012). 16. We assume droplet radii spread according to the gamma distribution [29] with the shape parameter ν 10. 17. K. M. Yoo and R. R. Alfano, Opt. Lett. 15, 320 (1990). 18. K. M. Yoo, Q. Xing, and R. R. Alfano, Opt. Lett. 16, 1019 (1991). 19. J. B. Schmidt, Z. D. Schaefer, T. R. Meyer, S. Roy, S. A. Danczyk, and J. R. Gord, Appl. Opt. 48, B137 (2009). 20. M. Linne, Prog. Energy Combust. Sci. 39, 403 (2013). 21. M. A. Linne, M. Paciaroni, E. Berrocal, and D. Sedarsky, Proc. Combust. Inst. 32, 2147 (2009). 22. L. V. Wang and H.-i. Wu, Biomedical Optics: Principles and Imaging (Wiley, 2007). 23. W.-F. Cheong, S. A. Prahl, and A. J. Welch, IEEE J. Quantum Electron. 26, 2166 (1990). 24. L. G. Henyey and J. L. Greenstein, Astrophys. J. 93, 70 (1941). 25. Our technique utilizes the partial-wave representation of the Green function G, rather than the operator G−1 − Σ, as in Refs. [14,15]; we intend to discuss the relation between the two approaches in a separate work. 26. Our cut and pole singularities are analogues of the continuous and discrete spectra in the Case’s method [4]. 27. A. McCarthy, R. J. Collins, N. J. Krichel, V. Fernández, A. M. Wallace, and G. S. Buller, Appl. Opt. 48, 6241 (2009). 28. A. Kirmani, D. Venkatraman, D. Shin, A. Colaço, F. N. C. Wong, J. H. Shapiro, and V. K. Goyal, Science 343, 58 (2014). 29. D. Deirmendjian, Appl. Opt. 3, 187 (1964).
OPTICS LETTERS / Vol. 39, No. 20 / October 15, 2014
Early-time diffusion in pulse propagation through dilute random media Elizabeth Bleszynski,* Marek Bleszynski, and Thomas Jaroszewicz Monopole Research, Thousand Oaks, California 91360, USA *Corresponding author: [email protected] Received July 7, 2014; accepted August 26, 2014; posted September 5, 2014 (Doc. ID 215211); published October 9, 2014 Propagation of short infrared/optical pulses in dilute random media (e.g., atmospheric clouds, fog, dust, or aerosols) consisting of large, compared to the wavelength, scatterers is analyzed. A rigorous approach based on analytic complex-contour integration of numerically determined cut and pole singularities of the radiative transport equation solution in the Fourier space is presented. It is found that the intensity of a propagating pulse, in addition to the coherent (“ballistic”) contribution and a long late-time diffusive tail, also exhibits a sharply rising early-time component that (i) can be attributed to the small-angle diffractive part of the scattering cross-section on medium particles, (ii) is attenuated proportionally to the nondiffractive rather than total cross-section, and (iii) can be extracted by high-pass filtering of the received pulse. © 2014 Optical Society of America OCIS codes: (030.5620) Radiative transfer; (290.1090) Aerosol and cloud effects; (290.1990) Diffusion; (290.2558) Forward scattering; (290.4210) Multiple scattering; (070.2615) Frequency filtering. http://dx.doi.org/10.1364/OL.39.005862
Imaging and communication with optical or infrared pulsed signals through obscuring random media (e.g., atmospheric clouds, fog, dust, or aerosols) is a longstanding and challenging problem, both experimentally and theoretically. From the experimental and applications perspective, a coherently detected pulse field preserves its time profile, but is strongly attenuated, at the rate proportional to the total cross-section of the wave on an individual medium scatterer. The incoherently detected field intensity [or, more generally, the mutual coherence function (MCF)], although attenuated at a lower rate (proportional to the absorption cross-section), develops a long diffusive tail that causes loss of resolution in imaging and loss of bandwidth in communication. From the theoretical point of view, pulse propagation in a statistically homogeneous, dilute, random medium can be adequately described by the radiative transport equation (RTE) [1–3]. However, obtaining its analytic solution for the time-dependent intensity (except for the case of isotropic scatterers [4,5]) proved to be a challenging problem. Consequently, most of the analysis and understanding of these transient phenomena have been based on Monte Carlo (MC) simulations (e.g., [6–9]) and some approximate methods (e.g., [10–12]). Only recently analytic solutions have been derived for arbitrary anisotropic scatterers [13–15] and applied in the context of biological media. This Letter presents first results of our analysis of propagation of short optical pulses in an atmospheric cloud-type medium. In this regime, the scatterers are relatively large compared to the wavelength (a typical average droplet radius [16] being a 5 μm ≈ 7.9λ0 for λ0 0.633 μm), the absorption is negligible, and the differential scattering cross-section on a medium constituent clearly exhibits a diffractive peak of angular width Δθ ∼ λ0 ∕2πa ≡ 1∕k0 a. We find, by means of the analytic solution technique that we briefly describe in this Letter, that the evolved pulse, in addition to the coherent (ballistic) contribution and a long diffusive tail, develops a distinctive 0146-9592/14/205862-04$15.00/0
early-time component characterized by (i) a rise time much shorter than the usual diffusion time scale, and (ii) attenuation significantly lower than that of the coherent intensity. Further, we show that, because of property (i), the early-time pulse component can be isolated by means of a “frequency-gating” technique implementable as high-pass filtering. In contrast to time-gating methods used widely for isolating an early-time part of the received pulse [17–19], the proposed gating technique does not require information on the expected arrival time of the signal and is, therefore, well suited in communication and in range imaging through a scattering medium. The early-time pulse structure can be interpreted as due to a diffusion-type process driven by a twodimensional angular random walk with steps of the order of the diffractive peak width, Δθ. Time delay in a single step is proportional to lt ∕cθ2 , where lt is the coherent attenuation length and θ the angle between the path segment and the observation direction. In N ∼ R∕lt steps needed to reach the distance R, the total delay becomes Δt ∼ Nlt ∕chθ2 i and, because of the angular random walk, hθ2 i ∼ NΔθ2 . Therefore, Δt ∼ Δθ2 R2 ∕clt , which is a relation characteristic of diffusion, but with a small coefficient Δθ2 ∼ 1∕k0 a2 . This justifies referring to the early-time pulse component as “early-time diffusion,” distinct from the usual latetime diffusion associated with the large-angle isotropic scattering. The above interpretation indicates that the appearance of the early-time diffusion component is due to a distinctive diffractive peak characteristic of scatterers large in comparison to the wavelength. Early time enhancement of evolved pulses has been noticed before in MC simulations [7] for very large (a ∼ 50 μm) strongly diffractive scatterers (Fig. 1) typical of sprays in combustion processes [20]. However, small length and time scales characterizing dense sprays (lt ∼ 1 cm, lt ∕c ∼ 3.3 ps) [21] make it difficult to experimentally separate the expected early-time diffusion © 2014 Optical Society of America
October 15, 2014 / Vol. 39, No. 20 / OPTICS LETTERS
5863
~ with the two-field Green function GΩ; P; sˆ μt − iΩ∕v0 − P · sˆ−1 . Here v0 ≈ c is the speed of coherent wave in the effective medium, μt 1∕lt is the attenuation coefficient, and the scattering function is given by Σˆs · sˆ 0 n0 σ s ˆs · sˆ 0 . The RTE (2) can be solved by projectingR the specific intensity onto partial-waves, I~ l Ω; P≔ d2 sˆ∕ ~ 4πpl Pˆ · sˆIΩ; P; sˆ , with normalized Legendre polynop mials pl x≔ 2l 1Pl x. Truncating the expansion at l < Lmax ∼ k0 a yields the matrix equation ~ ~ IΩ; P − GΩ; PΣIΩ; P GΩ; P0 ; Fig. 1. Small-angle Mie and Henyey–Greenstein (H–G) differential cross-sections at λ0 0.633 R 1 μm, normalized to total scatd cos θ σ s cos θ, for various tering cross-section σ s 2π −1 scatterer sizes.
signal from other scattering mechanisms. Observing early-time diffusion in biomedical imaging [22] is also difficult due to smaller scatterers [23], and thus insignificant diffraction [see Fig. 1, Mie scattering for red blood cells (a 2.78 μm) and Henyey–Greenstein [24] model for most biological tissues] combined with even higher scatterer concentration and smaller length and time scales. In our problem, on the other hand, droplets in an atmospheric cloud are large enough (a 5 μm) to give rise to considerable diffraction and the scatterer number density is much lower, typically n0 109 m−3 , resulting in the coherent attenuation length of lt 1∕n0 σ t ≈ 6 m. The corresponding attenuation time scale lt ∕c ≈ 20 ns is orders-of-magnitude longer; hence resolving the early-time signal should prove feasible [Eq. (7) and Fig. 4]. In the analysis described here, we adopt the customary approximate treatment of electromagnetic field components as scalar fields. We assume a coherent point-like isotropic source located at the origin in infinite space, emitting an intensity pulse St of duration T p , such that its physical length cT p is large compared to both the wavelength λ0 and the scatterer size a. Those two conditions are necessary in order to be able to reduce the Bethe–Salpeter equation (BSE) to its approximate RTE form: the first to simplify the BSE Green function and the second to neglect the dependence of the scattering amplitudes on the relative frequency (denoted below by Ω). The evolved field distribution is then described by the specific intensity (radiance) I S t; R; sˆ [1,2], where the unit vector sˆ ; is the energy flux direction at the time t and at the point R. The Fourier-space radiance I~ S Ω; P; sˆ≔
Z
dt d3 R eiΩt e−iP·R I S t; R; sˆ
(1)
~ ~ is given by I~ S Ω; P; sˆ SΩ IΩ; P; sˆ, where S~ is the Fourier transform of the source and I~ satisfies the RTE for an impulse (∼δt) source, ~ ~ IΩ; P; sˆ − GΩ; P; sˆ
Z
~ ~ d2 sˆ 0 Σˆs · sˆ 0 IΩ; P; sˆ 0 GΩ; P; sˆ (2)
(3)
where the vector I~ consists of Lmax compo~ nents R 1 I l , components of 2the vector Σ are Σl 2π −1 dx Pl x Σx, and the Lmax elements of the matrix G are R
Gl;l0 Ω; P
d2 sˆ ∕4πpl Pˆ · sˆ pl0 Pˆ · sˆGΩ; P; sˆ;
4
the r.h.s. of Eq. (3) is the l R0 column of G. The field intensity, It; R d2 sˆ ∕4πIt; R; sˆ, is now spherically symmetric and is given by Z It; R
∞
0
Z dP
∞
−∞
dΩ
P sinPR −iΩt ~ e I 0 Ω; P; 4π 3 R
(5)
where I~ 0 Ω; P I − GΩ; PΣ−1 GΩ; P0;0 is the l 0 partial wave of the solution of Eq. (3). While computing the Fourier-space solution I~ 0 Ω; P is straightforward, evaluation of the Fourier integral (5) (as discussed, e.g., in [14]) is nontrivial, especially for large scatterers and in the “early time,” near-light-cone region (v0 t ≳ R ≫ lt and v0 t − R ≪ lt ). To overcome those difficulties, we devised an approach in which we locate (partly analytically and partly numerically) singularities of I~ 0 Ω; P, perform contour integration over Ω, and integrate the result over real P [25]. The procedure relies on the fact that the function I~ 0 Ω; P has two types of singularities [26]: (i) Branch points at Ω v0 P − iμt and the associated cut in GΩ; P. The resultant cut contribution to the integral (5) is straightforward to calculate, and it decays with R at the same rate as the coherent intensity component (i.e., it becomes negligible for R ≫ lt ). (ii) Poles Ωj P, due to the matrix I − GΩ; PΣ not being invertible. Since the determinant of the matrix I − GΣ is holomorphic on the cut complex Ω plane, it may only have isolated zeros, resulting in a meromorphic function I~ 0 Ω; P, having only isolated poles. As P increases, the poles travel across the complex Ω plane until they reach the cut and disappear, one by one, from the “physical” Riemann sheet. The pole contribution to the intensity for a given P is, therefore, the sum of terms 1 It;Rj 2 2π R
Z 0
Pj
dPP sinPReImΩj Pt Imfαj Pe−iReΩj Pt g; (6)
where Ωj P is the jth pole trajectory, αj P is the pole residue, and P j is the value above which the pole
5864
OPTICS LETTERS / Vol. 39, No. 20 / October 15, 2014
disappears. For jPj ≥ maxfP j g ∼ k0 a2 μt , there are no more poles left on the physical sheet. This behavior of singularities can be understood by observing that, for P ≫ μt and k0 a2 ≫ 1, Eq. (2) becomes equivalent to the Schrödinger equation with an imaginary potential proportional to P −1 and the Fourier transform of Σ in the transverse components of sˆ . In the large P (thus weak coupling) regime, the discrete modes disappear and only the continuous spectrum remains. In Fig. 2 we plot, for propagation distances R1 12lt and R2 16lt , the dimensionless intensity J S v0 t∕lt ; tR∕lt l3t ∕v0 I S t; R, expressed as a function of dimensionless variables v0 t∕lt and R∕lt , and computed as the sum of cut and pole contributions. We traced the pole locations, and evaluated their residues, by means of eigenvalue decomposition of the matrix I − GΣ. We assumed a coherent source emitting a Gaussianmodulated field of carrier frequency c∕λ0 , with λ0 0.633 μm, resulting in intensity St exp−t2 ∕2T 2p ∕ p 2π T p , with T p 60 ps. The time-dependent intensity exhibits a distinctive early-time diffusion behavior, characterized by a narrow spike followed by a much broader late-diffusion maximum and a long diffusive tail. As the propagation distance increases, the broad diffusion shoulder starts overlapping the trailing edge of the early-time diffusion peak, but its sharply rising leading edge remains. Theoretical estimates and a number of computations for various scatterer sizes indicate that the rise time of the early-diffusion structure is approximately given by ΔtED ≈
κ ED R2 ; k0 a2 v0 lt
(7)
with a numerical factor κ ED ≲ 0.1. In the considered example, ΔtED ≈ 100 ps for R 12lt and ΔtED ≈ 180 ps for R 16lt ; these time scales are significantly smaller than the late-diffusion times (tens of nanoseconds). We note that the chosen pulse length T p < ΔtED allows us to take full advantage of the achievable range resolution cΔtED . Similar estimates indicate that the early-time diffusion intensity is attenuated, approximately, at the rate μED n0 σ t − σ d ;
proportional to the nondiffractive cross-section, rather than the total one (the diffractive cross-section σ d is the differential cross-section integrated over the diffractive forward-scattering cone). In our example the attenuation is reduced relative to the coherent attenuation by σ t − σ d ∕σ t ≈ 0.65. Both features, (7) and (8), can be explained in terms of the evolution of the poles appearing in Eq. (6). The two leading poles in our example problem are shown in Fig. 3. The j 0 pole emerges from the cut at P 0 ≈ 1674μt and, as long as P ≳ 10μt , the propagation velocity Re Ω0 P∕P remains very close to v0 , while the attenuation coefficient −Im Ω0 P∕v0 decreases from μt to about 60% of that value, resulting in an early, narrow, and relatively weakly attenuated contribution to the intensity. The j 1 pole, which emerges from the cut at a much smaller P value and moves up to Ω 0, is responsible for most of the late-time diffusion. The origin of the early-time diffusion component suggests that it should be possible to isolate it by selecting the region of large P values in the Fourier-space intensity; this amounts to high-pass filtering of the received timedomain intensity. In our example, we used a Gaussian ~ filter ΦΩ 1 − exp−T 2f Ω2 ∕2 of width v0 T f 0.015lt 9 cm, chosen such that v0 ΔtED < v0 T f ≪ lt . The original and filtered pulses, J S and Φ ∘ J S are plotted in Fig. 4 for distances R1 12lt and R2 16lt . The R2 intensity was scaled by the factor expμED R2 − R1 R2 ∕R1 2 ≈ 24. The close similarity between the R1 and the scaled R2 filtered intensities validates the reduced attenuation prediction of Eq. (8). The time shift in their maxima confirms the estimate of Eq. (7). The plots show that, already for the propagation distance of 12lt , the magnitude of the filtered early diffusive signal exceeds that of the coherent intensity by the factor ∼100. It is also evident that frequency gating does suppress the smooth late-time diffusion contribution while preserving the sharp leading edge of the early-time diffusion signal. To summarize, we expounded main properties (small rise time and reduced attenuation) of the early-time
(8)
Fig. 2. Intensity distributions showing early- and late-time diffusion behavior, for propagation distances R1 12lt and R2 16lt , with the time scale indicated on top. The curves are visually indistinguishable from the MC simulation results.
Fig. 3. Evolution, on the complex Ω plane, of two leading poles Ωj P, j 0 and j 1, of the integrand I~ 0 Ω; P of the Fourier transform (5). Each pole is either located on the imaginary axis or is accompanied by its mirror reflection in that axis. The thick horizontal line represents the cut and the arrows mark positions of poles at several values of P.
October 15, 2014 / Vol. 39, No. 20 / OPTICS LETTERS
Fig. 4. Received (top) and filtered (bottom) intensities for R1 12lt and R2 16lt . Coherent contributions, multiplied by factors 102 and 104 , are shown for reference.
diffusion component in the evolved short pulse intensity, and its relation to small-angle diffractive scattering; we also devised a filtering technique for extracting it from the received intensity. Our analysis indicates that, for cloud-type media, the expected duration of the early diffusion signal (∼100 ps) falls within resolution capabilities of current opto-electronic devices, and, at the same time, provides range resolution of a few centimeters. Our results may have important implications in highresolution range imaging as well as communication through obscuring (atmospheric clouds, fog, dust, or aerosols) media. Further developments may include (i) generalizations to other sources and detectors (limited detector field-of-view, collimated beams, source and detector arrays, partially coherent transmitted fields), (ii) generalization to transmission through a slab geometry, and (iii) alternative signal processing techniques in lowintensity measurements, such as single-photon counting [27] or first-photon-imaging [28] methods. This research is supported by the U. S. Air Force Office of Scientific Research under grant no. FA 9550-12-1-0121.
5865
References and Notes 1. S. Chandrasekhar, Radiative Transfer (Dover, 1960). 2. A. Ishimaru, Wave Propagation and Scattering in Random Media (Wiley, 1999). 3. L. Ryzhik, G. Papanicolaou, and J. B. Keller, Wave Motion 24, 327 (1996). 4. K. M. Case, Ann. Phys. 9, 1 (1960). 5. J. Palmeri, J. Stat. Phys. 58, 885 (1990). 6. E. A. Bucher, Appl. Opt. 12, 2391 (1973). 7. C. Rozé, T. Girasole, L. Méès, G. Gréhan, L. Hespel, and A. Delfour, Opt. Commun. 220, 237 (2003). 8. C. Calba, L. Méès, C. Rozé, and T. Girasole, J. Opt. Soc. Am. A 25, 1541 (2008). 9. B. Wu, Z. Hajjarian, and M. Kavehrad, Appl. Opt. 47, 3168 (2008). 10. A. Ishimaru and S. T. Hong, Radio Sci. 10, 637 (1975). 11. S. T. Hong and A. Ishimaru, Radio Science 11, 551 (1976). 12. W. Cai, M. Lax, and R. R. Alfano, Phys. Rev. E 63, 016606 (2000). 13. A. Liemert and A. Kienle, Phys. Rev. A 83, 015804 (2011). 14. A. Liemert and A. Kienle, Biomed. Opt. Express 3, 543 (2012). 15. A. Liemert and A. Kienle, Phys. Rev. E 86, 036603 (2012). 16. We assume droplet radii spread according to the gamma distribution [29] with the shape parameter ν 10. 17. K. M. Yoo and R. R. Alfano, Opt. Lett. 15, 320 (1990). 18. K. M. Yoo, Q. Xing, and R. R. Alfano, Opt. Lett. 16, 1019 (1991). 19. J. B. Schmidt, Z. D. Schaefer, T. R. Meyer, S. Roy, S. A. Danczyk, and J. R. Gord, Appl. Opt. 48, B137 (2009). 20. M. Linne, Prog. Energy Combust. Sci. 39, 403 (2013). 21. M. A. Linne, M. Paciaroni, E. Berrocal, and D. Sedarsky, Proc. Combust. Inst. 32, 2147 (2009). 22. L. V. Wang and H.-i. Wu, Biomedical Optics: Principles and Imaging (Wiley, 2007). 23. W.-F. Cheong, S. A. Prahl, and A. J. Welch, IEEE J. Quantum Electron. 26, 2166 (1990). 24. L. G. Henyey and J. L. Greenstein, Astrophys. J. 93, 70 (1941). 25. Our technique utilizes the partial-wave representation of the Green function G, rather than the operator G−1 − Σ, as in Refs. [14,15]; we intend to discuss the relation between the two approaches in a separate work. 26. Our cut and pole singularities are analogues of the continuous and discrete spectra in the Case’s method [4]. 27. A. McCarthy, R. J. Collins, N. J. Krichel, V. Fernández, A. M. Wallace, and G. S. Buller, Appl. Opt. 48, 6241 (2009). 28. A. Kirmani, D. Venkatraman, D. Shin, A. Colaço, F. N. C. Wong, J. H. Shapiro, and V. K. Goyal, Science 343, 58 (2014). 29. D. Deirmendjian, Appl. Opt. 3, 187 (1964).
