A model for full-field optical coherence tomography in scattering media

TL;DR

This paper presents a comprehensive model for full-field optical coherence tomography (FF-OCT) in scattering media, accounting for partial coherence and advanced scattering theory to improve understanding of imaging features and quantitative measurements.

Contribution

It introduces a novel FF-OCT model that incorporates partial coherence and goes beyond the Born approximation, enhancing the analysis of imaging and signal decay in scattering media.

Findings

The model describes how partial coherence affects the optical transfer function.

It explains the decay of FF-OCT signals with depth in scattering media.

Conditions are derived for exponential decay, enabling quantitative measurements using Beer-Lambert law.

Abstract

We develop a model of full-field optical coherence tomography (FF-OCT) that includes a description of partial temporal and spatial coherence, together with a mean-field scattering theory going beyond the Born approximation. Based on explicit expressions of the FF-OCT signal, we discuss essential features of FF-OCT imaging, such as the influence of partial coherence on the optical transfer function, and on the decay of the signal with depth that is captured by the model. We derive the conditions under which the spatially averaged signal exhibits a pure exponential decay with depth, providing a clear frame for the use of the Beer-Lambert law for quantitative measurements of the extinction length in scattering media.

Equations

$$\delta E_s({\mathbf{r}}_D)=k_0^2\int ⟨G({\mathbf{r}}_D,{\mathbf{r}}^{\prime })⟩\Delta \epsilon ({\mathbf{r}}^{\prime })⟨E({\mathbf{r}}^{\prime })⟩d^3{r}^{\prime },$$

$$⟨E(\mathbf{r})⟩=A\int ⟨G(\mathbf{r},{\mathbf{r}}_S)⟩E_0({\mathbf{r}}_S)d^2{\rho }_S$$

$$\delta E_s({\mathbf{r}}_D)=A{k}_0^2\int d^3{r}^{\prime }\int d^2{\rho }_S⟨G_T({\mathbf{r}}_D,{\mathbf{r}}^{\prime })⟩\Delta \epsilon ({\mathbf{r}}^{\prime })⟨G_T({\mathbf{r}}^{\prime },{\mathbf{r}}_S)⟩E_0({\mathbf{r}}_S).$$

$$⟨E_s({\mathbf{r}}_D)⟩=A\int ⟨G_R({\mathbf{r}}_D,{\mathbf{r}}_S)⟩E_0({\mathbf{r}}_S)d^2{\rho }_S.$$

$$S({\mathbf{r}}_D,\omega )=\overline{E_{\mathrm{ref}}^{\ast }({\mathbf{r}}_D)E_s({\mathbf{r}}_D)},$$

$$E_0(\mathbf{r})=\int d^{2}q{e}_0(\mathbf{q})\exp [i\mathbf{q}\cdot \mathbf{\rho }-i{k}_z(q)z],$$

$$S({\mathbf{r}}_D,\omega )=⟨S({\mathbf{r}}_D,\omega )⟩+\delta S({\mathbf{r}}_D,\omega )$$

$$⟨S({\mathbf{r}}_D,\omega )⟩$$

$$\delta S({\mathbf{r}}_D,\omega )$$

$$W_0({\mathbf{\rho }}_S-{\mathbf{\rho }}_D)=\overline{E_0^{\ast }({\mathbf{\rho }}_D,z=0)E_0({\mathbf{\rho }}_S,z=0)},$$

$$\displaystyle\langle G_{T}({\bf r},{\bf r}^{\prime})\rangle=\frac{i}{2\pi}\int d^{2}q\frac{g({\bf q})}{k_{z}(q)}\exp\big{[}i{\bf{q}}\cdot(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})+ik_{z}(q)\lvert z-z^{\prime}\rvert\big{]}\,.$$

$$\widetilde{\delta S}(\mathbf{q},\omega )=-4{\pi }^{2}A\int d{z}^{\prime }\exp [2i{k}_0(z_0-z^{\prime })]\widetilde{\Delta \epsilon }(\mathbf{q},z^{\prime })\exp (-|z^{\prime }|/{\ell }_e)$$

$$\times \int d^2{q}^{\prime }g({\mathbf{q}}^{\prime })g(\mathbf{q}+{\mathbf{q}}^{\prime }){\tilde{W}}_0({\mathbf{q}}^{\prime }),$$

$${\int }_{{\omega }_0-\Delta \omega /2}^{{\omega }_0+\Delta \omega /2}\exp [2i{k}_0(z_0-z^{\prime })]d\omega =\exp [2i{\bar{k}}_0(z_0-z^{\prime })]\mathrm{sinc}[(z_0-z^{\prime })\Delta \omega /c]\Delta \omega ,$$

$$\widetilde{\delta S}(\mathbf{q})=-4{\pi }^{3}Ac\widetilde{\Delta \epsilon }(\mathbf{q},z_0)\exp (-z_0/{\ell }_e)\int d^2{q}^{\prime }g({\mathbf{q}}^{\prime })g(\mathbf{q}+{\mathbf{q}}^{\prime }){\tilde{W}}_0({\mathbf{q}}^{\prime }).$$

$$\widetilde{\delta S}(\mathbf{q})=\tilde{H}(\mathbf{q}){\widetilde{\Delta \epsilon }}_{z_0}(\mathbf{q}),$$

$${\tilde{H}}_c(\mathbf{q})\sim I_{0}g(0)g(\mathbf{q}).$$

$${\tilde{H}}_i(\mathbf{q})\sim I_0^{\prime }\int d^2{q}^{\prime }g({\mathbf{q}}^{\prime })g(\mathbf{q}+{\mathbf{q}}^{\prime }).$$

$$k_z(q)$$

$$k_z^{\mathrm{eff}}(q)$$

$$\widetilde{\delta S}(\mathbf{q}=0)=-4{\pi }^{3}Ac\widetilde{\Delta \epsilon }(\mathbf{q}=0,z_0)\exp (-z_0/{\ell }_e)\int d^2{q}^{\prime }{g}^2({\mathbf{q}}^{\prime }){\tilde{W}}_0({\mathbf{q}}^{\prime })\exp [-q^{\prime 2}{z}_0/(2{\bar{k}}_0^2{\ell }_e)],$$

$${\nabla }^2⟨E(\mathbf{r})⟩+k_0^2{\epsilon }_{\mathrm{eff}}⟨E(\mathbf{r})⟩=0,$$

$${\nabla }^2⟨G(\mathbf{r},{\mathbf{r}}^{\prime })⟩+k_0^2{\epsilon }_{\mathrm{eff}}⟨G(\mathbf{r},{\mathbf{r}}^{\prime })⟩=-\delta (\mathbf{r}-{\mathbf{r}}^{\prime }),$$

$${\nabla }^{2}E(\mathbf{r})+k_0^2\epsilon (\mathbf{r})E(\mathbf{r})=0.$$

$${\nabla }^2\delta E_s(\mathbf{r})+k_0^2{\epsilon }_{\mathrm{eff}}\delta E_s(\mathbf{r})=-k_0^2[\epsilon (\mathbf{r})-{\epsilon }_{\mathrm{eff}}]E(\mathbf{r}).$$

$$\delta E_s(\mathbf{r})=k_0^2\int ⟨G(\mathbf{r},{\mathbf{r}}^{\prime })⟩[\epsilon ({\mathbf{r}}^{\prime })-{\epsilon }_{\mathrm{eff}}]E({\mathbf{r}}^{\prime })d^3{r}^{\prime }.$$

$$\delta E_s(\mathbf{r})=k_0^2\int ⟨G(\mathbf{r},{\mathbf{r}}^{\prime })⟩[\epsilon ({\mathbf{r}}^{\prime })-{\epsilon }_{\mathrm{eff}}]⟨E({\mathbf{r}}^{\prime })⟩d^3{r}^{\prime }.$$

$$E_0({\mathbf{r}}_S)=\exp (-i{k}_0{z}_S)\int d^{2}q{e}_0(\mathbf{q})\exp (i\mathbf{q}\cdot {\mathbf{\rho }}_S)\exp [i{q}^2{z}_S/(2k_0)],$$

$$E_{ref}^{\ast }({\mathbf{r}}_D)=\exp [-i{k}_0(z_D-2z_0)]\int d^{2}q{e}_0^{\ast }(\mathbf{q})\exp (-i\mathbf{q}\cdot {\mathbf{\rho }}_D)\exp [i{q}^2(z_D-2z_0)/(2k_0)].$$

$$\widetilde{\delta S}(\mathbf{q},\omega )=-4{\pi }^{2}A\int d{z}^{\prime }\exp [2i{k}_0(z_0-z^{\prime })]\widetilde{\Delta \epsilon }(\mathbf{q},z^{\prime })\exp (-|z^{\prime }|/{\ell }_e)$$

$$\times \exp [i{q}^2(z^{\prime }-z_D)/(2k_0)]\exp [-q^2|z^{\prime }|/(4k_0^2{\ell }_e)]\int d^2{q}^{\prime }g({\mathbf{q}}^{\prime })g(\mathbf{q}+{\mathbf{q}}^{\prime }){\tilde{W}}_0({\mathbf{q}}^{\prime })$$

$$\times \exp [iq{q}^{\prime }(z^{\prime }-z_D)/k_0]\exp [i{q}^{\prime 2}(z^{\prime }-z_0)/k_0]\exp [-(q{q}^{\prime }+q^{\prime 2})|z^{\prime }|/(2k_0^2{\ell }_e)].$$

$$\widetilde{\delta S}(\mathbf{q}=0,\omega )=-4{\pi }^{2}A\int d{z}^{\prime }\exp [2i{k}_0(z_0-z^{\prime })]\widetilde{\Delta \epsilon }(\mathbf{q}=0,z^{\prime })\exp (-|z^{\prime }|/{\ell }_e)$$

Full text

I Introduction

Since its initial development huang1991optical , Optical Coherence Tomography (OCT) has proven its ability to image inside scattering materials with micrometer resolution in three-dimensions schmitt1999optical . The possibility to access cellular structures in tissues at millimeter depth has been a breakthrough in biomedical optics, which stimulated a fast and broad dissemination of the technique.

The OCT setup is essentially a low coherence Michelson interferometer, in which one arm collects the light backscattered from the sample, while the other arm produces a reference beam reflected on a mirror. The main feature of an OCT setup is the ability to decouple the depth (longitudinal) and transverse resolutions fercher2003optical . Depth resolution is produced by temporal coherence gating, and is controlled by the spectral width of the incident light. Transverse resolution is controlled by the numerical aperture (NA) of the microscope objective in the sample arm. In scanning OCT (S-OCT) huang1991optical ; nassif2004vivo ; huber2006fourier ; hoeling2000optical ; zhang2005adaptive , a point-by-point image is formed by three-dimensional scanning of a focal spot. OCT systems recording en face images in planes perpendicular to the optical axis have also been developed, using spatially coherent illumination as in wide-field OCT (WD-OCT)  bourquin2001optical ; bordenave2002wide ; laubscher2002video , or spatially incoherent illumination as in full-field OCT (FF-OCT)  vabre2002thermal ; fercher2000thermal .

Since OCT is expected to collect the singly backscattered photons, the signal is substantially affected by multiple scattering, whose contribution overcomes the signal at depths larger than the scattering mean free path schmitt1997model . Several approaches have been followed to decrease the mutiple scattering contribution by reducing the weight of long light paths, including spatial filtering through confocal detection webb1996confocal , time gating hee1993femtosecond , or polarization gating schmitt1992use ; macdonald2017numerical . Other strategies address an inverse problem to correct a posteriori for multiple scattering, and increase resolution and penetration depth. Approaches based on interferometric synthetic aperture microscopy marks2009partially , or computational adaptive optics adie2012computational , have proven to be successful. Interestingly, it has been shown that the effect of multiple scattering depends on the degree of spatial coherence of the illuminating beam marks2009partially . It was also demonstrated that aberrations in the sample arm do not influence the transverse resolution in FF-OCT using spatially incoherent light xiao2016full . Recently, a method based on the measurement of the reflection matrix has demonstrated an efficient discrimination between singly and multiply scattered light, with an unprecedented increase in the OCT working depth badon2016smart , thus pushing the limits of optical microscopy in highly scattering media badon2017multiple .

Despite the success of OCT, from basic to clinical studies, a comprehensive theoretical model of the OCT signal, that handles a description of partial coherence together with a realistic scattering model (beyond the Born approximation), is still missing. In this paper we develop such a model, and use it to discuss different aspects of FF-OCT imaging. We study the influence of partial coherence on the optical transfer function. We also address the question of the decay of the signal with depth, that is captured by a mean-field scattering approach. Depending on the degree of spatial coherence and on the numerical aperture of the illumination/collection optics, we derive the conditions under which the spatially averaged FF-OCT signal exhibits a pure exponential decay with depth. The provides a clear frame for the use of the Beer-Lambert law for quantitative measurements of the extinction length in scattering media.

II Scattering model

We consider a sample made of large-scale inhomogeneities immersed in a scattering medium made of randomly distributed scattering centers, as represented in Fig. 1(a).The sample is characterized by a dielectric function $\varepsilon({\bf r})$ that we write as $\varepsilon({\bf r})=\varepsilon_{b}+\varepsilon_{d}({\bf r})+\delta\varepsilon({\bf r})$, where $\varepsilon_{b}$ denotes a uniform background, $\varepsilon_{d}({\bf r})$ is a large-scale deterministic dielectric function varying on a scale $L_{\varepsilon}\gtrsim\lambda$, with $\lambda$ the wavelength in vacuum, and $\delta\varepsilon({\bf r})$ is a real random variable describing a disordered distribution of small-scale scattering centers. We assume that $\delta\varepsilon({\bf r})$ satisfies $\langle\delta\varepsilon({\bf r})\rangle=0$, the brackets denoting a statistical average over an ensemble of realizations of the disordered scattering background, and that the correlation function $\langle\delta\varepsilon({\bf r})\delta\varepsilon({\bf r}^{\prime})\rangle$ is of the form $f(|{\bf r}-{\bf r}^{\prime}|/\ell_{\varepsilon})$, where $f$ is a positive decaying function with range close to unity. This defines $\ell_{\varepsilon}$ as the microscopic length scale of the disordered scattering medium, and we assume $\ell_{\varepsilon}\lesssim\lambda$. Note that this description is very general. For example by setting $\varepsilon_{d}({\bf r})=0$ we would describe a purely scattering medium, as that represented in Fig. 1(b) (OCT is sometimes used to measure the extinction length $\ell_{e}$ in such materials). By setting $\delta\varepsilon({\bf r})=0$ we would describe large-scale objects in a uniform background.

In the statistical approach, in which the scattering background is treated as a random medium, the scattered field is a random variable. In a given realization of the sample (deterministic large-scale objects immersed in one statistical realization of the scattering background), the total field can be written $E=\langle E\rangle+\delta E$. In OCT, one measures the backscattered field $E_{s}=E-E_{0}$, $E_{0}$ being the incident field, and we can also write $E_{s}=\langle E_{s}\rangle+\delta E_{s}$. In the regime of weak-scattering characterized by the condition $|\delta E_{s}|\ll|\langle E\rangle|$, multiple scattering in the background scattering medium can be described in the mean-field approximation. This approximation is based on well-established result in multiple scattering theory sheng2007introduction , that are summarized in Appendix A. In the mean-field approach, the fluctuating part $\delta E_{s}$ of the scattered field at point ${\bf r}_{D}$ in the detector plane, and at frequency $\omega$, reads

[TABLE]

where $k_{0}=\omega/c$ is the wavenumber in free space, $\langle G({\bf r}_{D},{\bf r}^{\prime})\rangle$ is the average Green’s function connecting a point ${\bf r}^{\prime}$ in the sample to the detection point ${\bf r}_{D}$, and the integration is over the volume of the sample. Here the average Green’s function $\langle G\rangle$ and the average field $\langle E\rangle$ are defined as solutions of a propagation equation in an effective medium, described by an effective dielectric function $\varepsilon_{\mathrm{eff}}$, representing the average contribution of the random scattering background (the detailed derivation of Eq. (1) is given in Appendix A). It is important to note that this effective medium description is a rigorous result of multiple scattering theory sheng2007introduction ; akkermans2007mesoscopic . Noting that $\Delta\varepsilon({\bf r}^{\prime})=\varepsilon({\bf r}^{\prime})-\varepsilon_{\mathrm{eff}}$ is the local dielectric contrast between a heterogeneity in the sample and the effective medium, a clear physical meaning can be given to Eq. (1): It expresses the fluctuating part $\delta E_{s}$ of the scattered field as the result of single scattering in the effective medium, a result valid in the weak-scattering regime $|\delta E_{s}|\ll|\langle E\rangle|$. The mean-field approach differs from the Born approximation by the renormalization of the background medium into an effective medium, that accounts in particular for the decay of the average Green’s function due to scattering, as we shall see. The average field is connected to the incident field $E_{0}$ in the source plane by the relation

[TABLE]

where $\langle G({\bf r},{\bf r}_{S})\rangle$ is the average Green’s function connecting a point ${\bf r}_{S}=(\boldsymbol{\rho}_{S},z_{S})$ in the source plane to an arbitrary point ${\bf r}$, $A$ is a constant that we do not need to specify, and the integration is along the source plane. Inserting Eq. (2) into Eq. (1) leads to

[TABLE]

In this expression, for the sake of clarity, we have denoted by $\langle G_{T}\rangle$ the Green’s function that accounts for propagation through the optics in the sample arm, transmission at the effective medium surface, and propagation inside the effective medium. The average scattered field at the detector is readily deduced from Eq. (2), and reads

[TABLE]

Here, the average Green’s function $\langle G_{R}\rangle$ accounts for propagation through the optics in the sample arm, and reflection at the effective medium surface. For practical calculations, the average Green’s function in transmission or reflection can be approximated using a simple model, as we shall see below.

III FF-OCT signal

An OCT setup is based on a Michelson interferometer, as represented schematically in Fig. 2. Starting from the source, the beam is divided by a beam-splitter to travel along two arms. The sample arm collects the field backscattered from the sample. In the second arm, a reference beam is produced by reflection on a mirror. The detector collects the intensity resulting from the interference between the sample and reference beams. In FF-OCT, a full-field illumination is used in combination with an array of detectors (in practice a CCD camera) to record the signal at multiple transverse locations in parallel, and produce an en face image. This means that both arms contain a microscope objective, not represented in Fig. 2 for simplicity. The state of coherence of the light source plays a crucial role in FF-OCT. Beyond the longitudinal sectioning given by the temporal coherence length (or equivalently the spectral bandwidth), spatial coherence influences the lateral resolution, as well as the sensitivity to aberrations in the sample arm xiao2016full . In order to understand precisely the role of temporal and spatial coherence, we need to build a model of the FF-OCT signal that includes the scattering model introduced in the previous section, and the main features of the interferometric and broadband detection used in OCT.

Under illumination by statistically stationary light, the intensity recorded by a Michelson interferometer is characterized by three terms: The (time) average intensity of the reference field, the average intensity of the scattered field, and the cross-correlation between them (interference term). Assuming $|E_{s}|\ll|E_{0}|$, the signal carrying information on the sample is the interference term (the intensity of the scattered field is negligible, and the intensity of the reference field only contributes as a background signal) davis2007autocorrelation . In the frequency domain, the OCT signal measured at a given point ${\bf r}_{D}$ in the detector plane is

[TABLE]

where $E_{\mathrm{ref}}$ and $E_{s}$ are, respectively, the complex amplitude of the reference and backscattered field at frequency $\omega$, the superscript $*$ denotes the complex conjugate, and the overline means a time averaging over the fluctuations of the partially coherent source.

The incident field $E_{0}$ at a point ${\bf r}=(\boldsymbol{\rho},z)$ can be written in the form of a plane-wave expansion

[TABLE]

with $k_{z}(q)=(k_{0}^{2}-q^{2})^{1/2}$. The $z$-axis is chosen such that the plane $z=0$ coincides with the sample plane (the sample filling the half-space $z<0$), and the axis $z>0$ follows the optical axis in each arm (see Fig.1(d)). Assuming a weakly focused beam (the influence the numerical aperture is discussed in section V), we can perform a zeroth-order paraxial approximation $k_{z}(q)\simeq k_{0}$, resulting in $E_{0}({\bf r})=E_{0}(\boldsymbol{\rho},z=0)\exp(-ik_{0}z)$. The field in the reference arm is assumed to coincide with the incident field longitudinally shifted by the mirror displacement. Under the same paraxial approximation, we can write $E_{\mathrm{ref}}^{*}({\bf r}_{D})=E^{*}_{0}(\boldsymbol{\rho}_{D},z=0)\exp[-ik_{0}(z_{D}-2z_{0})]$ where $z_{0}$ is the position of the mirror along the $z$-axis in the reference arm and $z_{D}$ is the position of the detection plane along the $z$-axis in the detection arm. From Eqs. (3), (4), (5) and the expressions of $E_{0}$ and $E_{\mathrm{ref}}^{*}$ above , we can write the FF-OCT signal as

[TABLE]

with

[TABLE]

In these expressions we have introduced the cross-spectral density

[TABLE]

that characterizes the state of coherence of the incident field in the sample plane $z=0$ mandel1995optical . Here we assume that the incoming intensity is uniform over the sample surface, and use a homogeneous Schell model such that $W_{0}$ is a function of $\boldsymbol{\rho}_{S}-\boldsymbol{\rho}_{D}$ only marks2009partially .

It is clear that $\langle S({\bf r}_{D},\omega)\rangle$ characterizes the effective medium, and does not carry information on $\Delta\varepsilon({\bf r}^{\prime})$. In order to build the simplest model, we assume weak scattering, meaning that $k_{0}\ell_{e}\gg 1$ with $\ell_{e}$ the extinction mean free path characterizing the decay of the average field (the exinction coefficient $\mu_{e}=1/\ell_{e}$ can also be used equivalently). The effective medium, as seen by the average field, is described by a dielectric function $\varepsilon_{\mathrm{eff}}=\varepsilon_{b}+i/(k_{0}\ell_{e})$ (this expression being valid to first order in the small parameter $1/(k_{0}\ell_{e})$) sheng2007introduction ; akkermans2007mesoscopic . In this limit, the reflected average Green’s function $\langle G_{R}\rangle$ simply describes reflection at the surface of the effective medium, and contributes as a uniform background. Assuming a low index mismatch with $\varepsilon_{b}\simeq 1$, the contribution of $\langle S({\bf r}_{D},\omega)\rangle$ becomes negligible. In the following we focus on the contribution $\delta S({\bf r}_{D},\omega)$ that carry the relevant information on the image formation process. It is useful to introduce the plane-wave expansion of the average Green’s function:

[TABLE]

We take $k_{z}(q)=(k_{0}^{2}-q^{2})^{1/2}$ for propagation outside the sample (for $z,z^{\prime}>0$), and $k_{z}(q)=k_{z}^{\mathrm{eff}}(q)=(\varepsilon_{\mathrm{eff}}k_{0}^{2}-q^{2})^{1/2}$ for propagation inside the sample (for $z,z^{\prime}<0$). The filter $g({\bf q})$ plays the role of a pupil function that limits the transverse wavevector ${\bf q}$ within a region bounded by the numerial aperture $NA$ of the objectives (the simplest model for $g({\bf q})$ is a disk with radius $k_{0}\,NA$) sentenac2018tutorial . For $g({\bf q})=1$ (infinite pupil), we recover the Weyl expansion of the free-space Green’s function. To get an expression of the FF-OCT signal relevant for an analysis in terms of optical transfer function, we insert Eq. (11) into Eq. (9), and perform again a zeroth-order paraxial approximation, using $k_{z}(q)\simeq k_{0}$ outside the sample, and $k_{z}^{\mathrm{eff}}(q)\simeq k_{0}+i/(2\ell_{e})$ inside the sample. The result is easily written in terms of the 2D Fourier transform of the signal, $\widetilde{\delta S}({\bf q},\omega)=\int\delta S({\bf r}_{D},\omega)\exp(-i{\bf q}\cdot\boldsymbol{\rho}_{D})d^{2}\rho_{D}/{4\pi^{2}}$, and reads

[TABLE]

where $\widetilde{\Delta\varepsilon}({\bf q},z^{\prime})$ is the 2D Fourier transform of $\Delta\varepsilon(\boldsymbol{\rho}^{\prime},z^{\prime})$, and $\widetilde{W}_{0}({\bf q}^{\prime})$ is the 2D Fourier transform of $W_{0}({\boldsymbol{\rho}}_{S}-{\boldsymbol{\rho}}_{D})$. Note that $\widetilde{W}_{0}({\bf q}^{\prime})\propto{I}_{0}({\bf q}^{\prime})$, with $I_{0}({\bf q}^{\prime})$ the angular distribution of the incident intensity barabanenkov1968radiation ; mandel1995optical . Also note that this expression does not rely on any assumption regarding the optics in the illumination, reference and detection arms, and does not depend on the definition of a focal plane for either illumination of detection (these features drive the precise form of the Green’s function $g({\bf q})$). In practice, in order to control the degree of spatial coherence in the sample plane $z=0$, one could choose to control the intensity distribution in the source plane $z=z_{S}$ and conjugate this plane with the plane $z=0$, but this practical choice does not influence the general form of Eq. (12).

A feature of OCT is to integrate the signal over a broad spectral range. For a statistically stationary source, the broadband signal is obtained by integrating the different frequency components over the source bandwidth $\Delta\omega$. Assuming that $W_{0}$ and $\Delta\varepsilon$ have a weak dependence on $\omega$, the spectral integration gives

[TABLE]

where $\bar{k}_{0}=\omega_{0}/c$, with $\omega_{0}$ the central frequency of the source. The $\mathrm{sinc}$ function, considered as a function of $z^{\prime}$, is centered at $z_{0}$ with a width $\ell_{\omega}=c/\Delta\omega$ corresponding to the temporal coherence length of the incident light, and is responsible for the longitudinal sectioning in OCT. When $\ell_{\omega}\ll\ell_{e}$, we can use the approximation $\mathrm{sinc}[(z_{0}-z^{\prime})/\ell_{\omega}]\simeq\pi\ell_{\omega}\delta(z_{0}-z^{\prime})$. We endup with a closed-form expression of the 2D Fourier transform of the broadband OCT signal $\widetilde{\delta S}({\bf q})=\int_{\Delta\omega}\widetilde{\delta S}({\bf q},\omega)d\omega$, that depends on the position $z_{0}$ of the mirror and on the source bandwith $\Delta\omega$:

[TABLE]

Expression (14) of the the FF-OCT signal implicitly involves different length scales, whose interplay is critical in the analysis of the signal: The length scale $\ell_{\varepsilon}$ characterizing the microscopic random inhomogeneities in the sample, the length scale $L_{\varepsilon}$ characterizing large-scale deterministic fluctuations of the dielectric function, the temporal coherence length $\ell_{\omega}$ and the spatial coherence length $\ell_{c}$ of the incident light. Expression (14) is similar to that previously derived in Ref. marks2009partially . The main difference is that our derivation includes a scattering model based on a mean-field approach, that only relies on the weak-scattering assumption $|\delta E_{s}|\ll|\langle E\rangle|$, always satisfied when the scattered field remains much smaller than the incident field. This approach goes beyond a first-order or second-order Born approximation, by accounting rigorously for the propagation of the average field in a renormalized effective medium. An interesting consequence is that the extinction of the signal with depth, due to scattering and absorption, emerges explicitly.

IV Spatial coherence and response function

Expression (14) is an interesting starting point for the analysis of the role of spatial coherence on the image formation in FF-OCT. To proceed, let us rewrite it in the compact form

[TABLE]

where $\widetilde{H}({\bf q})\sim\int d^{2}q^{\prime}g({\bf q}^{\prime})g({\bf q}+{\bf q}^{\prime})\widetilde{W}_{0}({\bf q}^{\prime})$ is the FF-OCT transfer function (we omit constant prefactors for simplicity), and $\widetilde{\Delta\varepsilon}_{z_{0}}({\bf q})=\widetilde{\Delta\varepsilon}({\bf q},z_{0})\exp(-z_{0}/\ell_{e})$ is the weighted 2D Fourier transform of the dielectric contrast. The width of the cross-spectral density $W_{0}(\boldsymbol{\rho})$ defines the spatial coherence length $\ell_{c}$ of the incident field in the sample plane $z=0$. When $\ell_{c}\gg L_{\varepsilon}\gg\ell_{\varepsilon}$, the incident light can be considered fully coherent, which corresponds to $\widetilde{W}_{0}({\bf q}^{\prime})\sim I_{0}\delta({\bf q}^{\prime})$, $I_{0}$ being proportionnal to the incident intensity. The transfer function for spatially coherent illumination becomes

[TABLE]

The coherent transfer function $\widetilde{H}_{c}({\bf q})$ posseses the same spatial frequency content as $g({\bf q})$, and covers the spatial frequency range $q\leq\bar{k}_{0}NA$, leading to a transverse spatial resolution limit $\lambda/(2NA)$, where $\lambda=2\pi/\bar{k}_{0}$ is the central wavelength of the incident field. The regime of spatially incoherent illumination corresponds formally to $W_{0}(\boldsymbol{\rho})\sim\delta(\boldsymbol{\rho})$, or equivalently $\widetilde{W}_{0}({\bf q}^{\prime})=I_{0}^{\prime}$. The transfer function for spatially incoherent illumination reads as

[TABLE]

In practice, since $\ell_{c}$ cannot be made smaller than $\lambda/2$ under far-field illumination, the condition $\ell_{c}\ll\ell_{\varepsilon}\ll L_{\varepsilon}$ that rigorously corresponds to $W_{0}(\boldsymbol{\rho})\sim\delta(\boldsymbol{\rho})$ is out of reach. Nevertheless, provided that FF-OCT is used to image the large-scale inhomonegeities at the scale $L_{\varepsilon}$, and not to resolve the small-scale random scattering centers, the condition of incoherent illumination can be relaxed to be $\ell_{c}\simeq\ell_{\varepsilon}\ll L_{\varepsilon}$. The integral term in Eq. (17) shows that the spatial frequency range encompassed by $\widetilde{H}_{i}({\bf q})$ is broader by a factor of two compared to coherent illumination, which is a usual result in optical microscopy, leading to a theoretical resolution limit $\lambda/(4NA)$. The results above, for coherent and incoherent illumination, are identical to those previsouly established in Ref. marks2009partially , and in agreement with the qualitative analysis of coherent and incoherent FF-OCT presented in Ref. sentenac2018tutorial .

Expressions (16) and (17) may also provide a theoretical frame to study the influence of aberrations in the sample arm on the image quality, recently discussed in Ref. xiao2016full . A precise study is beyond the scope of the present work. Qualitatively, we can simply note that weak aberrations that would tend to squeeze the Green’s function $g({\bf q})$ in Fourier space (without changing the cutoff frequency) would narrow the shape of $\widetilde{H}_{c}({\bf q})$, thus degrading the image quality. For incoherent illumination, the convolution product in Eq. (17) reduces the sensitivity of the transfer function to change in shapes due to weak aberrations, resulting in a better protection of the image quality.

V Depth dependence of the spatially averaged signal

The decay of the signal with depth is a feature of OCT. For paraxial illumination and detection, and in a medium with $\widetilde{\Delta\varepsilon}({\bf q},z_{0})$ independent of $z_{0}$, the signal follows an exponential decay $\exp(-z_{0}/\ell_{e})$, as described by Eq. (14). This exponential decay can be used for the measurement of the extinction length in weakly scattering materials plamann2019OCT . At higher numerical aperture, a correction to a pure exponential decay is expected, that may also depend on the degree of spatial coherence of the incident field. To address this question, we make use of a second-order paraxial approximation

[TABLE]

in the plane-wave expansions of the incident field [Eq. (6)] and of the average Green’s function [Eq. (11)]. Note that, as a result of the mean-field approach, the paraxial approximation in the scattering medium amounts to replacing $k_{0}$ by $\sqrt{\varepsilon_{\mathrm{eff}}}k_{0}$. This strongly influences the dependence of the FF-OCT on the numerical aperture, as we shall see. Following again the steps leading to Eq. (14), we end up with an expression of the FF-OCT signal that explicitely accounts for the angular aperture of the illumination and detection beams. In order to discuss the depth dependence of the spatially integrated signal, we consider $\widetilde{\delta S}({\bf q}=0)$, which can be cast in the following form (details of the derivation are given in Appendix B):

[TABLE]

where $\bar{k}_{0}=\omega_{0}/c$ is the central wavenumber of the polychromatic incident field. This expression is well suited to a discussion of the depth dependence of the integrated FF-OCT signal. Note that in practice, an exponential fit to the signal is often used to estimate the extinction mean free path $\ell_{e}$ in statistically homogeneous and isotropic scattering media. From Eq. (20), it is clear that at finite $NA$ a deviation from a pure exponential decay $\exp(-z_{0}/\ell_{e})$ may be observed due to the dependence on $z_{0}$ of the integral over $q^{\prime}$. This integral contains several cutoffs. First, the last exponential term gives a depth-dependent cutoff $q^{\prime}_{z_{0}}\simeq\bar{k}_{0}(\ell_{e}/z_{0})^{1/2}$. For $z_{0}\gg\ell_{e}$, the integral is restricted to vanishingly small $q^{\prime}$, and a pure exponential decay $\exp(-z_{0}/\ell_{e})$ is always expected in the tail of the signal versus depth. Second, the angular spectrum of the Green’s function $g({\bf q}^{\prime})$ produces a cutoff $q^{\prime}_{NA}\simeq\bar{k}_{0}NA$, showing that for $NA\to 0$ a pure exponential decay is observed, in agreement with Eq. (14). Third, the Fourier transform of the cross spectral density $\widetilde{W}_{0}({\bf q}^{\prime})$ introduces a cutoff $q^{\prime}_{c}\simeq 2\pi/\ell_{c}$ that depends on the degree of spatial coherence of the incident light in the sample plane. For $\ell_{c}\to\infty$ (spatially coherent illumination), a pure exponential decay is also observed. More precisely, to prevent the $\exp[-q^{\prime 2}z_{0}/(2\bar{k}_{0}^{2}\ell_{e})]$ term in the integral to play a role, we need either $NA\ll(\ell_{e}/z_{0})^{1/2}$ (low numerical aperture regime), or $\ell_{c}\gg\lambda(z_{0}/\ell_{e})^{1/2}$ (coherent illumination).

This qualitative analysis is confirmed by numerical calculations of the spatially integrated signal using Eq. (20), as shown in Fig. 3. A non-exponential decay is observed in the regime $NA\simeq(\ell_{e}/z_{0})^{1/2}\simeq\lambda/\ell_{c}$. The deviation from a pure exponential decay is enhanced at high $NA$ and low spatial coherence. Although the curves are displayed in the regime $z_{0}\simeq\ell_{e}$, we have verified that an exponential decay $\exp(-z_{0}/\ell_{e})$ is recovered in any case when $z_{0}\gg\ell_{e}$. Numerical calculations based on ab initio simulations could extend the analysis beyond the weak-scattering regime and second-order paraxial approximations used in this model. This is left for future work.

VI Conclusion

In summary, we have developed a model of FF-OCT that includes a mean-field scattering theory, in addition to a precise description of temporal and spatial coherence. The model describes the decay of the signal with depth due to scattering and absorption, an essential feature in OCT imaging. It also permits an analysis of the interplay between different length scales characterizing the scattering medium and the degree of coherence of the incident field. Based on explicit expressions of the FF-OCT signal, we have discussed several features of FF-OCT imaging. We have analyzed the influence of spatial coherence on the optical transfer function, and discussed the particular cases of fully coherent and incoherent illuminations. We have also studied the depth dependence of the signal integrated over the transverse directions, that not only limits the penetration depth in OCT, but is also used to measure the extinction length in scattering materials. For spatially incoherent illumination, and/or with high numerical aperture of the illumination/detection optics, deviations from a pure exponential decay $\exp(-z_{0}/\ell_{e})$ can be observed. A pure exponential decay is always recovered in the tail of the signal versus depth when $z_{0}\gg\ell_{e}$. Our analysis provides a clear frame for the use of the Beer-Lambert law for quantitative measurements of the extinction length in scattering media. The model also provides a framework for a precise analysis of the role of aberrations generated by the optics in the sample arm, or by the scattering medum itself, and for the development of advanced inverse reconstruction procedures going beyond the first and second-order Born approximations.

Funding

This work was supported by the European Union’s Seventh Framework Program under HELMHOLTZ grant agreement No 610110, and by LABEX WIFI (Laboratory of Excellence ANR-10-LABX-24) within the French Program “Investments for the Future” under reference ANR-10-IDEX-0001-02 PSL*∗*.

Acknowledgments

We thank A. Aubry, V. Barolle, C. Boccara, R. Bocheux, K. Irsch and K. Plamann for stimulating discussions.

Appendix A Elements of multiple scattering theory

In order to build a mean-field expression of the scattered field, we first consider a purely scattering medium with dielectric function $\varepsilon_{b}+\delta\varepsilon({\bf r})$, as represented in Fig. 1(b). We assume that the random variable $\delta\varepsilon({\bf r})$, that describes a spatial distribution of small-scale scattering centers, is statistically homogeneous and isotropic. An important result of multiple scattering theory, states that the average field obeys a propagation equation in an effective homogeneous medium, as represented in Fig. 1(c). We take this result as our starting point (for a derivation, see for example Refs. sheng2007introduction ; akkermans2007mesoscopic ). More precisely, the average field obeys

[TABLE]

where $\varepsilon_{\mathrm{eff}}$ is the effective dielectric function that describes the average contributio of the random scattering medium. Equivalently, the average Green’s function satisfies the Dyson equation

[TABLE]

with an outgoing wave condition at infinity. This effective medium approach is strictly valid provided that $\ell_{\varepsilon}\ll\lambda$, meaning that non-locality in the effective dielectric function can be neglected, $\ell_{\varepsilon}$ being the microscopic clength scale associted to $\delta\varepsilon({\bf r})$ sheng2007introduction . Interestingly, the imaginary part of the effective dielectric constant describes the attenuation of the average field by scattering (and absorption), and defines the extinction mean free path $\ell_{e}$ such that $\mathrm{Im}\varepsilon_{\mathrm{eff}}=(k_{0}\ell_{e})^{-1}$. In absence of absorption, the extinction mean free path coincides with the scattering mean free path $\ell_{s}$.

In the general situation represented in Fig. 1(a), in which large-scale deterministic objects are superimposed to the random scattering background, the dielectric function is $\varepsilon({\bf r})=\varepsilon_{b}+\varepsilon_{d}({\bf r})+\delta\varepsilon({\bf r})$. The total field obeys

[TABLE]

In order to obtain the equation satisfied by the fluctuating scattered field $\delta E_{s}=E-\langle E\rangle$, we subtract Eq. (21) to Eq. (23), which leads to

[TABLE]

Using the average Green’s function $\langle G({\bf r},{\bf r}^{\prime})\rangle$ defined in Eq. (22), the solution to Eq. (24) is shown to obey the following integral equation:

[TABLE]

Assuming that scattering is sufficiently weak for the condition $|\delta E_{s}|\ll|\langle E\rangle|$ to be valid, the above expression can be approximated by

[TABLE]

Equation (26) is a mean-field approximation to Eq. (25). Physically, it describes the scattered field in one realization of the medium as resulting from a single scattering process in the effective medium.

Appendix B FF-OCT signal in the second-order paraxial approximation

In this appendix we develop the steps leading to Eq. (20). We start by performing the expansion (18) in Eq. (6), resulting in the following expressions of the incident and reference fields:

[TABLE]

Next we use expansions (18) and (19) in the expressions (11) of the average Green’s function. Here the expansion to second order in $q$ is performed only in the exponential term, keeping a zeroth-order expansion $k_{z}(q)\simeq k_{0}$ in the denominator that plays a minor role in the overall $z$ dependence. Using the resulting expressions of $E_{0}$, $E_{ref}^{*}$ and $\langle G_{T}\rangle$ into Eqs. (3) and (5), we obtain

[TABLE]

This expression extends Eq. (12) to incident and detection beams with non-negligible numerical apertures. We now focus on the expression of the signal integrated over the transverse direction, which is obtained by taking ${\bf q}=0$:

[TABLE]

It can be verified by a numerical evaluation that the integral over $q^{\prime}$ weakly depends on frequency over a bandwidth $\Delta\omega$ corresponding to a depth resolution on the order of one micrometer, which corresponds to $\ell_{\omega}\simeq 1\mu m$. Making use of Eq. (13) in the limit $\ell_{\omega}\ll\ell_{e}$, the signal integrated over frequencies becomes

[TABLE]

where $\bar{k}_{0}=\omega_{0}/c$ is the central wavenumber of the incident field. This expression is Eq. (20) of the main text.