The inverse scattering problem for orthotropic media in polarization-sensitive optical coherence tomography

In this paper we provide for a first time, to our knowledge, a mathematical model for imaging an anisotropic, orthotropic medium with polarization-sensitive optical coherence tomography. The imaging problem is formulated as an inverse scattering problem in three dimensions for reconstructing the electrical susceptibility of the medium using Maxwell’s equations. Our reconstruction method is based on the second-order Born-approximation of the electric field.


Introduction
Optical coherence tomography (OCT) is an imaging technique producing highresolution images of the inner structure of biological tissues. Standard  broadband, continuous wave light for illumination and the images are obtained by measuring the time delay and the intensity of the backscattered light from the sample. For a detailed description of OCT systems we refer to the books (Bouma and Tearney 2002;Drexler and Fujimoto 2015) and for a mathematical modeling we refer to Elbau et al. (2015).
Apart form standard OCT, there exist also functional OCT techniques such as the polarization-sensitive OCT (PS-OCT) which considers the differences in the polarization state of light to determine the optical properties of the sample. PS-OCT is based on polarization-sensitive low coherence interferometry established by Hee et al. (1992) and then first applied to produce two-dimensional OCT images (de Boer et al. 1997(de Boer et al. , 1998. In this work, we consider the basic scheme of a PS-OCT system which consists of a Michelson interferometer with the addition of polarizers and quarter-wave plates (QWP).
More precisely, a linear polarizer is added after the source and the linear (horizontal or vertical) polarized light is split into two identical parts by a polarization-insensitive beam splitter (BS). In the reference arm the light is reflected by a mirror and in the sample arm the light is incident on the medium. At the BS, the back-reflected beam and the backscattered light from the sample, in an arbitrary polarization state, are recombined. The recombined light passes through a polarizing BS which splits the output signal into its horizontal and vertical components to be measured at two different photo detectors. See Fig. 1 for an illustration of this setup.
To describe the change in the polarization state of the light due to its propagation into the sample we adopt the analysis based on the theory of electromagnetic waves scattered by anisotropic inhomogeneous media (Colton and Kress 2013;Wolf and Foley 1989). We assume that the dielectric medium is linear and anisotropic. In addition, we impose the property that the medium is invariant under reflection by the x 1 −x 2 plane. A medium with this property is called orthotropic in the mathematical com-

BS Source
Polarizer D x D y PBS QWP1 Mirror QWP2 Sample Fig. 1 Schematic representation of the light travelling in a time-domain PS-OCT system. In the reference and sample arms are placed quarter-wave plates (QWP) at specific orientations munity (Cakoni and Colton 2014) or monoclinic in the material science community (Torquato 2002).
The medium is also considered as weakly scattering and we present the solution in the accuracy of the second-order Born-approximation. As we are going to see later, we consider higher-order approximation in order to be able to recover all the material parameters. We describe the change in the polarization state of the light by the Jones matrix formalism which is applicable since OCT detects the coherent part of the electric field of the backscattered light (Jiao and Wang 2002). As in standard OCT systems, the backscattered light is detected in the far field.
In the medical community, the sample is usually described by a general retarder and the change in the polarization state of the light returning from the sample can be modelled by a Jones matrix (Hitzenberger et al. 2001;Jiao and Wang 2002). However, even though the produced images are satisfactory they are mainly used qualitatively. The usage of these images comes only secondarily to quantify the optical parameters by image processing techniques.
In this work we are interested in the quantitative description of PS-OCT. To do so, we have first to describe mathematically the system properly. Thus, we represent the polarized scattered field as solution to the full-wave Maxwell's equations. This has not yet been applied to PS-OCT, since for isotropic media, the Born-approximation decouples the effects of the optical properties of the sample from the polarization state of the scattered field. However, this analysis for anisotropic media provides enough information to reconstruct the electric susceptibility of the medium. The scattered field satisfies then an integral equation of Lippmann-Schwinger type. Under the far-field approximation and the assumption of a homogeneous background medium we obtain a system of integral equations for the unknown optical parameters.
In the mathematical literature, the scattering problem by anisotropic objects has been widely considered over the last decades (Beker and Umashankar 1989;Geng et al. 2003;Graglia and Uslenghi 1984;Papadakis et al. 1990). Recently, the connection between the inverse problem to reconstruct the refractive index and the interior transmission problem has been investigated (Cakoni et al. 2010;Cakoni and Haddar 2007). For the specific case of an orthotropic medium we refer the reader to the book (Cakoni and Colton 2014) and to Colton et al. (1997) and Potthast (1999) for results concerning the uniqueness and existence of solutions of the inverse problem.
The paper is organized as follows: In Sect. 2, we derive the integral representation of the scattered field for an orthotropic medium in the accuracy of the second-order Bornapproximation in the far-field zone. In Sect. 3, we describe mathematically the standard PS-OCT system using the Jones matrix formalism and we derive the expression for the cross-spectral density. The system of equations for all the components of the susceptibility is presented in the last section using two incident illuminations.

Notation
In this paper, we use the following conventions: -Let f : R → C be integrable, then the one-dimensional Fourier-transform is defined byf -Let f : R → C be integrable, then the one-dimensional inverse Fourier-transform is defined byf -Let f : R 3 → C be integrable, then the three-dimensional Fourier-transform is defined byf

The direct scattering problem
In absence of external charges and currents, the propagation of electromagnetic waves in a non-magnetic medium is mathematically described by Maxwell's equations relating the electric and magnetic fields E : R × R 3 → R 3 and H : R × R 3 → R 3 and the electric displacement D : where c is the speed of light. Maxwell's equations are not sufficient to uniquely determine the fields D, E and H. Therefore additional material parameters have to be specified: Definition 1 An anisotropic medium is called linear dielectric if there exists a function, called the electric susceptibility, A linear dielectric medium is called orthotropic (Cakoni and Colton 2014;Colton et al. 1997) if it admits the special symmetric form Application of the Fourier transform to Maxwell's equations (1) and taking into account (2), it follows that the Fourier-transform E of E satisfies the vector Helmholtz equation 2. For every x ∈ R 3 the function can be extended to a square integrable, holomorphic function on the upper half plane 3. E solves the Lippmann-Schwinger integral equation where is the fundamental solution of the scalar Helmholtz equation.
The integral operator G is strongly singular and we address its properties in the last section.
Proof From the initial condition (7) it follows for every solution E of Maxwell's equations (1) which fulfills (7) that the inverse Fourier-transform of g satisfieš Thus, the second assertion is a direct consequence of the Paley-Wiener theorem (Papoulis 1962).
To prove the first part, note that the electric field E is uniquely defined by (4) together with the assumption that the function g can be for every x ∈ R 3 extended to a square integrable, holomorphic function on the upper half plane.
Finally, the solution of Eq. (4) can be written as the solution of the integral Eq. (8), see Cakoni and Colton (2014) and Potthast (2000).

Born and far-field approximation
We assume that the medium is weakly scattering, meaning thatχ is sufficiently small (Chew 1990;Colton and Kress 2013) such that the incident field E i is significantly

Definition 3
The first order Born-approximation of the solution E of the Lippmann-Schwinger equation (8) is defined by The second order Born-approximation is defined by Inserting (9) into (10) gives or in coordinate writing where now G is the Green tensor of Maxwell's equations (Haddar 2004;Hazard and Lenoir 1996) The physical meaning of the second order Born-approximation is that at a point x the total field E 2 contains all single and double scattering events.
In an OCT setup, the measurements are performed in a distance much bigger compared to the size of the sample. Thus, setting x = ρϑ, ρ > 0 and ϑ ∈ S 2 , we can replace the above expression by its asymptotic behaviour for ρ → ∞, uniformly in ϑ, see for instance Elbau et al. (2015, Equation (4.1)), resulting to Here we have introduced the operator defined for functions f : R × R 3 → R 3 .

Polarized-sensitive OCT
We describe the standard PS-OCT system in the context of a Michelson interferometer first presented by Hee et al. (1992). The detector array is given by D = R 2 ×{d} with d > 0 sufficiently large. Moreover, we specify the CIF function to be E 0 as defined in Example 1 and we assume that We describe now the change in the polarization state of the light through the PS-OCT system. The effect of the polarization-insensitive beam splitter (BS) is not considered in this work since it only reduces the intensity of the beam by a constant factor. For simplicity, we place the sample and the mirror around the origin and the detector at the BS, for more details see Elbau et al. (2015, Section 3.3). The BS splits the light into two identical beams entering both arms of the interferometer.
Reference arm: The light (at some negative time) passes through a zero-order quarter-wave plate (QWP) oriented at angle φ 1 to the incident linear polarization. It is reflected by a perfect mirror placed in x 3 = l and then passes through the QWP again, at time t = 0, see the right picture in Fig. 2. We formulate this process as a linear operator to be specified later. Then, the reference field E l takes the form Sample arm: The incoming light passes through a QWP (oriented at a different angle φ 2 ) at some time t < 0, placed in the plane given by the equation x 3 = l Q . This process results to a field that until t = 0 does not interact with the medium, see the left picture in Fig. 2. Detector: The electric field E which is obtained by illuminating the sample with the incident field E 0,inc is combined with the reference field E l . We assume here that the backscattered light does not pass through the QWP again. At every point on the detector surface D we measure the two intensities (Elbau et al. 2015) We assume that we do not measure the incident fields at the detector, meaning E 0 (t, ξ ) = E 0,inc (t, ξ ) = 0 for t > 0 and ξ ∈ D and recalling (16) we obtain The two scattering problems in PS-OCT. On the left picture the incoming light in the sample arm passes through a QWP and is incident on the medium. On the right picture, in the reference arm, the light is back-reflected by a mirror (passing twice a QWP) We use Plancherel's theorem, and since E ∈ R 3 it follows that E(−ω, ·) = E(ω, ·). Thus, the above formula can be rewritten as

Jones calculus
Here we describe the operators J l and J , introduced in (15) and (17), respectively. We consider the fields in the frequency domain. Then, for positive frequencies we can apply the Jones matrix method (keeping also the zero third component of the fields) in order to model the effect of the QWP's on the polarization state of light. We assume that the properties of the QWP's are frequency independent and that the light is totally transmitted through the plate surface.

Definition 4
We define where is the rotated Jones matrix for a QWP with the fast axis oriented at angle φ (Gerrard and Burch 1975).
The above definition summarizes what we described before: In the reference arm, the incoming field passes through the QWP (at angle φ 1 ) is reflected by the mirror and then passes through the QWP again. The field travels additionally the distance 2(x 3 − l). In the sample arm, the field passes only through the QWP at angle φ 2 .
Now we can define our approximated data. We approximate in (19) the termÊ j − E 0,inc j byÊ 2 j −Ê 0,inc j and for the termÊ l j −Ê 0 j we consider (16) and (21).

Definition 5 We call
the second order approximated measurement data of OCT.

The inverse problem of recovering the susceptibility
The problem we address here, is to recoverχ from the knowledge of I 2 (l, ξ ) for l ∈ R, ξ ∈ D. First, we show that the measurements provide us with expressions which depend onχ non-linearly.
Proof We consider Eq. (13) where now E i is replaced by E 0,inc for ω > 0. Then, we We apply the inverse Fourier transform with respect to l in (22), to obtain which forω > 0,f = 0 and η j = 0, using that E and f are real valued, results to This identity together with (24), results asymptotically to (23).
We observe here that we want to reconstruct four four-dimensional functions from two three-dimensional measurement data. Thus, we have to consider some additional assumptions on the medium in order to cancel out the lack of dimensions and handle the non-linearity of (23) with respect toχ .
Assumption 1 Specific illumination: The support of the initial pulse is small enough such that the optical parameters in this spectrum can be assumed constant with respect to frequency.
Medium: The susceptibility can be decomposed into two parts, a background susceptibility which is constant and assumed to be known and a part that counts for the local variations of the susceptibility and can be seen as deviation from the constant value.
In the following, we consider this type of media, which is typical for biological tissues, and we assume in addition that the behavior of the homogeneous medium ( = 0) is known. Then, as a consequence, also the measured data from PS-OCT are known, let us call them I 0 , and we can assume the following form for the measurements for some known functions M j .
Proposition 2 Let the assumptions of Proposition 1 and the additional Assumption 1 hold. We define v = ω c (ϑ + e 3 ), ϑ ∈ S 2 + . Then, the spatial Fourier transform of the matrix-valued function ψ : R 3 → C 3×3 , satisfies the equations The operators K and K † are defined by for functions f : R 3 → C 3×3 , with kernels for α = z, y.
Proof We substituteχ , considering Assumption 1, and (26) in (23) and we equate the first order terms ψ and M to obtain In order to analyse the left hand side of the above equation we consider the definition (14) and the analytic form (12). Then, we rewrite (30) as where m j is given by (28). Taking the Fourier transform of ψ with respect to space, we get This equation form(v) := m(ω, ϑ), using the definitions of the integral operators (29) admits the compact form (27).
Regarding the integral operators appearing in (29), we prove the following property.
Proof We consider the following decomposition The above expression in compact form reads for the operators The operator F : L 2 (Ω) → L 2 (S 2 ) is a modification of the usual far-field operator with smooth kernel, thus compact. The operators G : L 2 (Ω) → L 2 (Ω) and G 1 : (L 2 (Ω)) 3×3 → (L 2 (Ω)) 3×3 are also compact due to their weakly singular kernels, see for instance (Colton and Kress 2013;Potthast 2000), and the operator G 0 : (L 2 (Ω)) 3×3 → (L 2 (Ω)) 3×3 is bounded (Colton et al. 2007). Thus K : (L 2 (Ω)) 3×3 → (L 2 (S 2 )) 3×3 is also compact. The same arguments hold also for the operator K † . Now, we are in position to formulate the inverse problem: Recover from the expressions the matrix-valued function ψ : Ω → C 3×3 , assuming that we have measurements for every incident polarization. Let us now specify the polarization vectors η and p. We choose two different incident polarization vectors q (1) = e 1 and q (2) = e 2 , and using the formulas (21) we obtain the vectors Remark 1 To find, for instance, the form of the incident wave p (1)f (ω) e −i ω c x 3 , for ω > 0, in the time domain we have to extend it for negative frequencies and consider its inverse Fourier transform. Then, we have If the small spectrum is centered around a frequency ν, we approximatef (ω) δ(ω − ν), for ω > 0, to obtain We see that E (1) describes also a circularly polarized wave with a phase shift.
If we neglect the zeroth third components, we observe that η (1) , η (2) ∈ R 2 and p (1) , p (2) ∈ C 2 form a basis in R 2 and C 2 , respectively. The following result shows that measurements for additional polarization vectors q do not provide any further information.
Proposition 3 Let ϑ ∈ S 2 + be fixed and q = q (1) , q (2) . Then, the Eq. (27) is equivalent to the system of equations where Y :=ψ(v) j , k, j = 1, 2, and P ϑ denotes the orthogonal projection in direction ϑ. The upper index on the data counts for the different incident polarizations.
We see that Proposition 3, for for k, j = 1, 2 and two different polarization vectors q = e 1 and q = e 2 uniquely determine the projections [ P ϑ Y p (k) ] j for k, j ∈ {1, 2}.
Moreover, measurements for additional polarizations q do not provide any further informations so that at every detector point, corresponding to a direction ϑ ∈ S 2 + , only the four elements [ P ϑ Y p (k) ] j , k, j = 1, 2, of the projection influence the measurements.
Remark 2 In contrast to standard OCT where three polarization vectors were needed, see Elbau et al. (2015, Proposition 11), and to first order Born-approximation where Y =ψ, as we are going to see in the following, the above measurements due to the special form of Y allow for reconstructing all the unknowns functions ψ i j .
Proof In order to reformulate equations (33), first we consider an arbitrary vector p and we split the expression P ϑ Y p into the sum omitting for simplicity the v dependence of the unknownψ. The first term on the right hand side admits the decomposition where we observe the independence onψ 33 . To analyse the other two terms, we consider (29) and define the operators acting now on the components of the matrixvalued function f : for k, j = 1, 2, 3. Since we are interested only in the first two components of P ϑ Y p and the calculations are rather lengthy we are going to omit the third component in the following expressions. The second term on the right hand side of (37) reads The only term whereψ 33 appears is the last one (as expected), namely We can combine now all the above formulas to obtain and y = ψ 11 ,ψ 12 ,ψ 22 ,ψ 33 . Then, the system of Eq. (33), considering (32) reads We observe that in all equations the coefficient in front of the operator M is the same, which is the only operator applying on the fourth component of y. In addition, from (38), we see that ϑ 2 M 14 = ϑ 1 M 24 . Thus, in order to eliminate y 4 we reformulate the above system as follows: we subtract from Eq. (39a) the Eq. (39c), from Eq. (39b) the Eq. (39d) and from ϑ 2 · (39a) the equation ϑ 1 · (39b), resulting to I( p (1) ) + χ 0 L( p (1) ) + χ 0 M) y] 1 −ϑ 1 [(I( p (1) ) + χ 0 L( p (1) ) + χ 0 M) y] 2 = ϑ 2 b (1) (1) 2 .
Recall that ϑ ∈ S 2 + , meaning ϑ 3 > 0. Then, if in addition we impose that ϑ 1 = ϑ 2 for all ϑ ∈ S 2 + , the matrixĨ is invertible withĨ −1 = det(Ĩ) −1 adj(Ĩ). Then, Eq. (40) can be written in the form which is the Fredholm integral equation of the second kind (35), for C :=Ĩ −1 N , Once (41) is solved for y 1 , y 2 and y 3 we can choose one of the four equations from the system (39) resulting to a Fredholm integral equation of the first kind for the unknown y 4 now: for some known function b, depending onỹ andb. This is Eq. (36) for C := M 3 . The compactness of the integral operators C and C follows from the compactness of the operators K and K † , see Lemma 2, since they are operators that act on the components of the matrix-valued function.
Remark 3 Equation (36) reflects the ill-posedness of the inverse problem, due to the compactness of the integral operator.

Conclusions
In this work we have formulated the inverse problem of recovering the electric susceptibility of a non-magnetic, inhomogeneous orthotropic medium, placed in a polarized-sensitive optical coherence tomograph (PS-OCT) as a system of Fredholm integral equations (both of first and second kind). Under the assumptions of a nondispersive, weakly scattering medium with small background variations we have shown that we can reconstruct all the coefficients of the matrix-valued susceptibility, given the data for two different incident polarization vectors. This paper can be seen, on one hand, as a first attempt to model mathematically PS-OCT and on the other hand, as a theoretical basis for an upcoming paper where the numerical validation of the proposed method will be examined.