The inverse scattering problem for orthotropic media in Polarization-sensitive Optical Coherence Tomography

TL;DR

This paper develops a mathematical model for imaging orthotropic media using PS-OCT, formulating it as an inverse scattering problem and proposing a reconstruction method based on the second-order Born approximation.

Contribution

It introduces a new mathematical framework for reconstructing the electrical susceptibility of orthotropic media in PS-OCT using Maxwell's equations.

Findings

Mathematical model for anisotropic media in PS-OCT

Inverse scattering formulation for 3D reconstruction

Reconstruction method based on second-order Born approximation

Abstract

In this paper we provide a mathematical model for imaging an anisotropic, orthotropic medium with Polarization-Sensitive Optical Coherence Tomography (PS-OCT). 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.

Figures

Equations

$$\hat{f}(\omega )={\int }_{\text{mathbbm}R}f(t){\mathrm{e}}^{\mathrm{i}\omega t}\mathrm{d}t,\text{ for all }\omega \in \text{mathbbm}R.$$

$$\check{f}(t)=\frac{1}{2\pi }{\int }_{\text{mathbbm}R}f(\omega ){\mathrm{e}}^{-\mathrm{i}\omega t}\mathrm{d}\omega ,\text{ for all }t\in \text{mathbbm}R.$$

$$\tilde{f}(\mathbf{k})={\int }_{\text{mathbbm}{R}^3}f(\mathbf{x}){\mathrm{e}}^{-\mathrm{i}⟨\mathbf{k},\mathbf{x}⟩}\mathrm{d}\mathbf{x},\text{ for all }\mathbf{x}\in \text{mathbbm}{R}^3.$$

$$\mathrm{div}\mathbf{D}(t,\mathbf{x})$$

$$\mathrm{div}\mathbf{H}(t,\mathbf{x})$$

$$\mathrm{curl}\mathbf{E}(t,\mathbf{x})$$

$$\mathrm{curl}\mathbf{H}(t,\mathbf{x})$$

$$\mbox{\boldmath$\chi$}\in C^{\infty}_{\mathrm{c}}(\mathbbm{R}\times\mathbbm{R}^{3};\mathbbm{R}^{3\times 3}),\text{ with }\mbox{\boldmath$\chi$}(\tau,\bm{x})=0\text{ for all }\tau<0,\bm{x}\in\mathbbm{R}^{3},$$

$$\bm{D}(t,\bm{x})={\bm{E}}(t,\bm{x})+\int_{\mathbbm{R}}\mbox{\boldmath$\chi$}(\tau,\bm{x}){\bm{E}}(t-\tau,\bm{x})\,\mathrm{d}\tau\;.$$

$$\mbox{\boldmath$\chi$}=\begin{pmatrix}\chi_{11}&\chi_{12}&0\\ \chi_{12}&\chi_{22}&0\\ 0&0&\chi_{33}\end{pmatrix}.$$

$$\operatorname{\mathbf{curl}}\operatorname{\mathbf{curl}}\widehat{\bm{E}}(\omega,\bm{x})-\frac{\omega^{2}}{c^{2}}(\mathbbm{1}+\hat{\mbox{\boldmath$\chi$}}(\omega,\bm{x}))\widehat{\bm{E}}(\omega,\bm{x})=0,\quad\omega\in\mathbbm{R},\;\bm{x}\in\mathbbm{R}^{3}.$$

$$\mathrm{curl}\mathrm{curl}{\hat{\mathbf{E}}}^i(\omega ,\mathbf{x})-\frac{{\omega }^2}{c^2}{\hat{\mathbf{E}}}^i(\omega ,\mathbf{x})=0,\text{ and }\mathrm{div}{\hat{\mathbf{E}}}^i(\omega ,\mathbf{x})=0,\omega \in \text{mathbbm}R,\mathbf{x}\in \text{mathbbm}{R}^3,$$

$${\mathbf{E}}^0(t,\mathbf{x})=\mathbf{q}f(t+\frac{x_3}{c}),$$

$$\mathrm{supp}{\mathbf{E}}^0(t,\cdot )\cap \Omega =\varnothing \text{ for every }t\le 0.$$

$$\mathrm{curl}\mathrm{curl}{\mathbf{E}}^0$$

$$=-\frac{1}{c^2}{f}^{\prime \prime }(t+\frac{x_3}{c})\mathbf{q}=-\frac{1}{c^2}{\partial }_{tt}{\mathbf{E}}^0.$$

$$\mathbf{E}(t,\mathbf{x})={\mathbf{E}}^i(t,\mathbf{x}),t\le 0,\mathbf{x}\in \text{mathbbm}{R}^3.$$

$$g:\text{mathbbm}R$$

$$\omega$$

$$\text{mathbbm}H=\{\omega \in \text{mathbbm}C\mid \mathrm{\Im m}(\omega )>0\}.$$

$$\hat{\mathbf{E}}(\omega ,\mathbf{x})$$

$$\displaystyle=:\widehat{\bm{E}}^{i}(\omega,\bm{x})+\bm{\mathcal{G}}[\hat{\mbox{\boldmath$\chi$}}\widehat{\bm{E}}](\omega,\bm{x})\,,$$

$$G(\omega ,\mathbf{x})=\frac{{\mathrm{e}}^{\mathrm{i}\frac{\omega }{c}|\mathbf{x}|}}{4\pi |\mathbf{x}|},\mathbf{x}\ne 0,\omega \in \text{mathbbm}R$$

$$\check{g}(t)=0\text{ for all }t\le 0.$$

$$\widehat{\bm{E}}^{1}=\widehat{\bm{E}}^{i}+\bm{\mathcal{G}}[\hat{\mbox{\boldmath$\chi$}}\widehat{\bm{E}}^{i}].$$

$$\widehat{\bm{E}}^{2}=\widehat{\bm{E}}^{i}+\bm{\mathcal{G}}[\hat{\mbox{\boldmath$\chi$}}\widehat{\bm{E}}^{1}].$$

$$\widehat{\bm{E}}^{2}=\widehat{\bm{E}}^{i}+\bm{\mathcal{G}}[\hat{\mbox{\boldmath$\chi$}}\widehat{\bm{E}}^{i}]+\bm{\mathcal{G}}\left[\hat{\mbox{\boldmath$\chi$}}\bm{\mathcal{G}}[\hat{\mbox{\boldmath$\chi$}}\widehat{\bm{E}}^{i}]\right],$$

$${\hat{\mathbf{E}}}^2(\omega ,\mathbf{x})$$

$$\displaystyle\phantom{=}+\frac{\omega^{4}}{c^{4}}\int_{\mathbbm{R}^{3}}\int_{\mathbbm{R}^{3}}\bm{G}(\omega,\bm{x}-\bm{y})\hat{\mbox{\boldmath$\chi$}}(\omega,\bm{y})\bm{G}(\omega,\bm{y}-\bm{z})\hat{\mbox{\boldmath$\chi$}}(\omega,\bm{z})\widehat{\bm{E}}^{i}(\omega,\bm{z})\,\mathrm{d}\bm{z}\,\mathrm{d}\bm{y},$$

$$\mathbf{G}(\omega ,\mathbf{x}-\mathbf{y})=G(\omega ,\mathbf{x}-\mathbf{y})\text{mathbbm}1+\frac{c^2}{{\omega }^2}\mathrm{grad}\mathrm{div}(G(\omega ,\mathbf{x}-\mathbf{y})\text{mathbbm}1).$$

$$\widehat{\bm{E}}^{2}(\omega,\rho\bm{\vartheta})=\widehat{\bm{E}}^{i}(\omega,\rho\bm{\vartheta})+\bm{\mathcal{G}}^{\infty}[\hat{\mbox{\boldmath$\chi$}}\widehat{\bm{E}}^{i}](\omega,\rho\bm{\vartheta})+\bm{\mathcal{G}}^{\infty}\left[\hat{\mbox{\boldmath$\chi$}}\bm{\mathcal{G}}[\hat{\mbox{\boldmath$\chi$}}\widehat{\bm{E}}^{i}]\right](\omega,\rho\bm{\vartheta}).$$

$$\bm{\mathcal{G}}^{\infty}[\bm{f}](\omega,\rho\mbox{\boldmath$\vartheta$}):=-\frac{\omega^{2}\mathrm{e}^{\mathrm{i}\tfrac{\omega}{c}\rho}}{4\pi\rho c^{2}}\int_{\mathbbm{R}^{3}}\bm{\vartheta}\times\left(\bm{\vartheta}\times\bm{f}(\omega,\bm{y})\right)\mathrm{e}^{-\mathrm{i}\frac{\omega}{c}\left<\bm{\vartheta},\bm{y}\right>}\,\mathrm{d}\bm{y},$$

$${\mathcal{J}}_l[{\mathbf{E}}^0](t,\mathbf{x})={\mathbf{E}}^{0,ref}(l;t,\mathbf{x}),$$

$${\mathbf{E}}^l(t,\mathbf{x})=\{\begin{array}{ll}{\mathbf{E}}^0(t,\mathbf{x})+{\mathbf{E}}^{0,ref}(l;t,\mathbf{x}),& \text{if }t>0,{\mathbf{x}}_3>l_R,\\ {\mathbf{E}}^0(t,\mathbf{x}),& \text{if }t\le 0,{\mathbf{x}}_3>l_R.\end{array}$$

Full text

1Computational Science Center

University of Vienna

Oskar-Morgenstern-Platz 1

A-1090 Vienna, Austria

2Johann Radon Institute for Computational

and Applied Mathematics (RICAM)

Altenbergerstraße 69

A-4040 Linz, Austria

1. Introduction

Optical Coherence Tomography (OCT) is an imaging technique producing high-resolution images of the inner structure of biological tissues. Standard OCT uses 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 [2, 12] and for a mathematical modeling we refer to [13].

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. [19] and then first applied to produce two-dimensional OCT images [10, 11]. 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 Figure 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 [7, 28]. 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 community [3] or monoclinic in the material science community [27].

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 [21]. 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 [20, 21]. 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 [1, 14, 16, 22]. Recently, the connection between the inverse problem to reconstruct the refractive index and the interior transmission problem has been investigated [4, 5]. For the specific case of an orthotropic medium we refer the reader to the book [3] and to [8, 24] for results concerning the uniqueness and existence of solutions of the inverse problem.

The paper is organized as follows: In Section 2, we derive the integral representation of the scattered field for an orthotropic medium in the accuracy of the second-order Born-approximation in the far-field zone. In Section 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:\mathbbm{R}\to\mathbbm{C}$ be integrable, then the one-dimensional Fourier-transform is defined by

[TABLE]

• •

Let $f:\mathbbm{R}\to\mathbbm{C}$ be integrable, then the one-dimensional inverse Fourier-transform is defined by

[TABLE]

• •

Let $f:\mathbbm{R}^{3}\to\mathbbm{C}$ be integrable, then the three-dimensional Fourier-transform is defined by

[TABLE]

2. 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 ${\bm{E}}:\mathbbm{R}\times\mathbbm{R}^{3}\to\mathbbm{R}^{3}$ and $\bm{H}:\mathbbm{R}\times\mathbbm{R}^{3}\to\mathbbm{R}^{3}$ and the electric displacement $\bm{D}:\mathbbm{R}\times\mathbbm{R}^{3}\to\mathbbm{R}^{3}$ by

[TABLE]

where $c$ is the speed of light. Maxwell’s equations are not sufficient to uniquely determine the fields $\bm{D},{\bm{E}}$ and $\bm{H}$. Therefore additional material parameters have to be specified:

Definition \thedefinition.

• •

An anisotropic medium is called linear dielectric if there exists a function, called the electric susceptibility,

[TABLE]

satisfying

[TABLE]

• •

A linear dielectric medium is called orthotropic [3, 8] if it admits the special symmetric form

[TABLE]

Application of the Fourier transform to Maxwell’s equations (1) and taking into account (2), it follows that the Fourier-transform $\widehat{\bm{E}}$ of ${\bm{E}}$ satisfies the vector Helmholtz equation

[TABLE]

Definition \thedefinition.

We call an electric field ${\bm{E}}^{i}:\mathbbm{R}\times\mathbbm{R}^{3}\to\mathbbm{R}^{3}$ a causal initial field (CIF) with respect to some domain $\Omega\subseteq\mathbbm{R}^{3}$ if

1. i.

its Fourier transform with respect to time solves Maxwell’s equations (1) with a susceptibility \mbox{\boldmath\chi}=0, that is,

[TABLE] 2. ii.

and satisfies $\operatorname{supp}{\bm{E}}^{i}(t,\cdot)\cap\Omega=\emptyset\text{ for every }t\leq 0$.

The second condition means that ${\bm{E}}^{i}$ does not interact with the medium contained in $\Omega$ until the time $t=0.$

Example \theexample.

Let $\Omega\subset\mathbbm{R}^{3}$ be an open set, such that \operatorname{supp}\mbox{\boldmath\chi}(t,\cdot)\subset\Omega for all $t\in\mathbbm{R}.$ Moreover, let $\bm{q}\in\mathbbm{R}^{2}\times\{0\}$ (denoting the polarization vector), $f\in C^{\infty}_{\mathrm{c}}(\mathbbm{R})$ and

[TABLE]

such that

[TABLE]

Then ${\bm{E}}^{0}$ is a CIF.

Proof**:**

To see this note that for arbitrary $\bm{q}\in\mathbbm{R}^{3}$ we get

[TABLE]

And for the particular choice $\bm{q}\in\mathbbm{R}^{2}\times\{0\}$ we even have that $\operatorname{div}{\bm{E}}^{0}=0$. This shows that ${\bm{E}}^{0}$ is a solution of Maxwell’s equation. The second assertion is an immediate consequence of the second assumption. $\square$

Theorem \thetheorem.

Let ${\bm{E}}^{i}$ be a CIF-function and assume that the susceptibility $\chi$ represents a dielectric, orthotropic medium. Then,

1. i.

there exists a solution ${\bm{E}}$ (together with $\bm{H}$) of Maxwell’s equations (1) which satisfies

[TABLE] 2. ii.

For every $\bm{x}\in\mathbbm{R}^{3}$ the function

[TABLE]

can be extended to a square integrable, holomorphic function on the upper half plane

[TABLE] 3. iii.

$\widehat{\bm{E}}$ solves the Lippmann–Schwinger integral equation

[TABLE]

where

[TABLE]

is the fundamental solution of the scalar Helmholtz equation.

The integral operator $\bm{\mathcal{G}}$ is strongly singular and we address its properties in the last section.

Proof**:**

1. ii.

From the initial condition (7) it follows for every solution ${\bm{E}}$ of Maxwell’s equations (1) which fulfills (7) that the inverse Fourier-transform of $g$ satisfies

[TABLE]

Thus, the result is a direct consequence of the Paley–Wiener theorem [23]. 2. i.

The electric field $\widehat{\bm{E}}$ is uniquely defined by (4) together with the assumption that the function $g$ can be for every $x\in\mathbbm{R}^{3}$ extended to a square integrable, holomorphic function on the upper half plane. 3. iii.

The solution of equation (4) can be written as the solution of the integral equation (8), see [3, 25].

$\square$

2.1. Born and Far-field approximation

We assume that the medium is weakly scattering, meaning that \hat{\mbox{\boldmath\chi}} is sufficiently small [6, 7] such that the incident field $\widehat{\bm{E}}^{i}$ is significantly larger than $\widehat{\bm{E}}-\widehat{\bm{E}}^{i}$.

Definition \thedefinition.

The first order Born-approximation of the solution $\widehat{\bm{E}}$ of the Lippmann-Schwinger equation (8) is defined by

[TABLE]

The second order Born-approximation is defined by

[TABLE]

Inserting (9) into (10) gives

[TABLE]

or in coordinate writing

[TABLE]

where now $\bm{G}$ is the Green tensor of Maxwell’s equations [17, 18]

[TABLE]

The physical meaning of the second order Born-approximation is that at a point $\bm{x}$ the total field $\widehat{\bm{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 $\bm{x}=\rho\bm{\vartheta},\rho>0$ and $\bm{\vartheta}\in\mathbb{S}^{2},$ we can replace the above expression by its asymptotic behaviour for $\rho\to\infty,$ uniformly in $\bm{\vartheta},$ see for instance [13, Equation (4.1)], resulting to

[TABLE]

Here we have introduced the operator

[TABLE]

defined for functions $\bm{f}:\mathbbm{R}\times\mathbbm{R}^{3}\to\mathbbm{R}^{3}$.

3. Polarized-sensitive OCT

We describe the standard PS-OCT system in the context of a Michelson interferometer first presented by Hee et al. [19].

The detector array is given by $\mathcal{D}=\mathbbm{R}^{2}\times\{d\}$ with $d>0$ sufficiently large. Moreover, we specify the CIF function to be ${\bm{E}}^{0}$ as defined in Section 2 and we assume that ${\bm{E}}^{0}(t,\bm{x})=0$ for $t\geq 0$ and $\bm{x}\in\mathcal{D}.$

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 [13, 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 $\phi_{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 Figure 2. We formulate this process as a linear operator

[TABLE]

to be specified later. Then, the reference field ${\bm{E}}^{l}$ takes the form

[TABLE]

Sample arm:

The incoming light passes through a QWP (oriented at a different angle $\phi_{2})$ at some time $t<0,$ placed in the plane given by the equation $x_{3}=l_{Q}.$ This process results to a field

[TABLE]

that until $t=0$ does not interact with the medium, see the left picture in Figure 2.

Detector:

The electric field ${\bm{E}}$ which is obtained by illuminating the sample with the incident field ${\bm{E}}^{0,inc}$ is combined with the reference field ${\bm{E}}^{l}.$ We assume here that the backscattered light does not pass through the QWP again. At every point on the detector surface $\mathcal{D}$ we measure the two intensities [13]

[TABLE]

We assume that we do not measure the incident fields at the detector, meaning ${\bm{E}}^{0}(t,\bm{\xi})={\bm{E}}^{0,inc}(t,\bm{\xi})=0$ for $t>0$ and $\bm{\xi}\in\mathcal{D}$ and recalling (16) we obtain $\bm{E}^{l}-\bm{E}^{0}=0$ for $t\leq 0,$ resulting to

[TABLE]

We use Plancherel’s theorem, and since ${\bm{E}}\in\mathbbm{R}^{3}$ it follows that $\widehat{\bm{E}}(-\omega,\cdot)=\overline{\widehat{\bm{E}}}(\omega,\cdot).$ Thus, the above formula can be rewritten as

[TABLE]

3.1. Jones Calculus

Here we describe the operators $\bm{\mathcal{J}}_{l}$ and $\bm{\mathcal{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 \thedefinition.

We define

[TABLE]

where

[TABLE]

is the rotated Jones matrix for a QWP with the fast axis oriented at angle $\phi$ [15].

The above definition summarizes what we described before: In the reference arm, the incoming field passes through the QWP (at angle $\phi_{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 $\phi_{2}.$

We consider the PS-OCT system, presented first in [19] and then considered in [20, 26], where $\phi_{1}=\pi/8$ and $\phi_{2}=\pi/4.$ Then, we obtain

[TABLE]

where $\bm{\eta}=\bm{J}^{2}(\pi/8)\,\bm{q}$ and $\bm{p}=\bm{J}(\pi/4)\,\bm{q}.$ We observe that $\widehat{\bm{E}}^{0,ref}$ is still linearly polarized at angle $\pi/4$ with the linear (horizontal or vertical) initial polarization state and $\widehat{\bm{E}}^{0,inc}$ describes a circularly polarized light.

Now we can define our approximated data. We approximate in (19) the term $\hat{E}_{j}-\hat{E}_{j}^{0,inc}$ by $\hat{E}^{2}_{j}-\hat{E}_{j}^{0,inc}$ and for the term $\hat{E}^{l}_{j}-\hat{E}_{j}^{0}$ we consider (16) and (21).

Definition \thedefinition.

We call

[TABLE]

the second order approximated measurement data of OCT.

4. The inverse problem of recovering the susceptibility

The problem we address here, is to recover \hat{\mbox{\boldmath\chi}} from the knowledge of $\bm{I}^{2}(l,\bm{\xi})$ for $l\in\mathbbm{R},$ $\bm{\xi}\in\mathcal{D}.$ First, we show that the measurements provide us with expressions which depend on \hat{\mbox{\boldmath\chi}} non-linearly.

Proposition \theproposition.

Let ${\bm{E}}^{0}(t,\bm{x})$ be given by the form (6) with $q_{3}=0$ and let the measurement data $I^{2}_{j}$ be given by (22). Then, for every $\omega\in\mathbbm{R}_{+}\setminus\{0\}$ with $\hat{f}(\omega)\neq 0$, the expression

[TABLE]

holds for all $j\in\{1,2\}$, and \mbox{\boldmath\vartheta}\in\mathbb{S}^{2}_{+}:=\{\mbox{\boldmath\mu}\in\mathbb{S}^{2}\mid\mu_{3}>0\}.

Proof**:**

We consider equation (13) where now $\widehat{\bm{E}}^{i}$ is replaced by $\widehat{\bm{E}}^{0,inc}$ for $\omega>0$. Then, we get

[TABLE]

We apply the inverse Fourier transform with respect to $l$ in (22), to obtain

[TABLE]

which for $\tilde{\omega}>0,\,\hat{f}\neq 0$ and $\eta_{j}\neq 0,$ using that $\bm{E}$ and $f$ are real valued, results to

[TABLE]

This identity together with (24), results asymptotically to (23). $\square$

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 \hat{\mbox{\boldmath\chi}}.

Assumption \theassumption.

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.

Then, the expression (3) admits the special form

[TABLE]

where

[TABLE]

for some known $\chi_{0}\in\mathbbm{R}$, a small parameter $\epsilon>0$ and $\psi_{ij}\in C^{\infty}_{\mathrm{c}}(\mathbbm{R}^{3};\mathbbm{C}).$

In the following, we consider this type of media, which are typical for biological tissues, and we assume in addition that the behavior of the homogeneous medium ($\epsilon=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

[TABLE]

for some known functions $M_{j}.$

Proposition \theproposition.

Let the assumptions of Section 4 and the additional Section 4 hold. We define \bm{v}=\tfrac{\omega}{c}(\bm{\vartheta}+\bm{e}_{3}),\,\mbox{\boldmath\vartheta}\in\mathbb{S}^{2}_{+}. Then, the spatial Fourier transform of the matrix-valued function $\bm{\psi}:\mathbbm{R}^{3}\rightarrow\mathbbm{C}^{3\times 3},$ satisfies the equations

[TABLE]

where

[TABLE]

The operators $\bm{\mathcal{K}}$ and $\bm{\mathcal{K}}^{\dagger}$ are defined by

[TABLE]

for functions $\bm{f}:\mathbbm{R}^{3}\rightarrow\mathbbm{C}^{3\times 3},$ with kernels

[TABLE]

Proof**:**

We substitute \hat{\mbox{\boldmath\chi}}, considering Section 4, and (26) in (23) and we equate the first order terms $\bm{\psi}$ and $\bm{M}$ to obtain

[TABLE]

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

[TABLE]

where $m_{j}$ is given by (28). Taking the Fourier transform of $\psi$ with respect to space, we get

[TABLE]

This equation for $\bm{\tilde{m}}(\bm{v}):=\bm{m}(\omega,\bm{\vartheta}),$ using the definitions of the integral operators (29) admits the compact form (27). $\square$

Now, we are in position to formulate the inverse problem: Recover from the expressions

[TABLE]

the matrix-valued function $\bm{\psi}:\Omega\rightarrow\mathbbm{C}^{3\times 3},$ if we assume that we have measurements for every incident polarization.

Let us now specify the polarization vectors $\bm{\eta}$ and $\bm{p}.$ We choose two different incident polarization vectors $\bm{q}^{(1)}=\bm{e}_{1}$ and $\bm{q}^{(2)}=\bm{e}_{2},$ and using the formulas (21) we obtain the vectors

[TABLE]

Remark \theremark:

To find, for instance, the form of the incident wave $\bm{p}^{(1)}\hat{f}(\omega)\,\mathrm{e}^{-\mathrm{i}\tfrac{\omega}{c}x_{3}},$ for $\omega>0,$ in the time domain we have to extend it for negative frequencies and consider its inverse Fourier transform. Then, we have

[TABLE]

If the small spectrum is centered around a frequency $\nu,$ we approximate $\hat{f}(\omega)\simeq\delta(\omega-\nu),$ for $\omega>0,$ to obtain

[TABLE]

We see that $\bm{E}^{(1)}$ describes also a circularly polarized wave with a phase shift.

If we neglect the zeroth third components, we observe that $\bm{\eta}^{(1)},\bm{\eta}^{(2)}\in\mathbbm{R}^{2}$ and $\bm{p}^{(1)},\bm{p}^{(2)}\in\mathbbm{C}^{2}$ form a basis in $\mathbbm{R}^{2}$ and $\mathbbm{C}^{2},$ respectively. The following result shows that measurements for additional polarization vectors $\bm{q}$ do not provide any further information.

Proposition \theproposition.

Let $\bm{\vartheta}\in\mathbb{S}^{2}_{+}$ be fixed and $\bm{q}=\bm{q}^{(1)},\,\bm{q}^{(2)}.$ Then, the equation (27) is equivalent to the system of equations

[TABLE]

where \bm{Y}:=\bm{\tilde{\psi}}(\bm{v})+\mbox{\boldmath\chi}_{0}\bm{\mathcal{K}}[\bm{\tilde{\psi}}](\bm{v})+\bm{\mathcal{K}}^{\dagger}[\bm{\tilde{\psi}}](\bm{v})\mbox{\boldmath\chi}_{0}, $b^{(k)}_{j}:=-\sqrt{2}\tilde{m}^{(k)}_{j},\,k,j=1,2,$ and $\bm{P}_{\bm{\vartheta}}$ denotes the orthogonal projection in direction $\bm{\vartheta}$. The upper index on the data counts for the different incident polarisations.

Proof**:**

The system of equations (27) for $(\bm{q},j)\in\{(\bm{q}^{(1)},1),(\bm{q}^{(1)},2),(\bm{q}^{(2)},1),(\bm{q}^{(2)},2)\}$ is equivalent to the four equations

[TABLE]

Indeed, for arbitrary polarisation $\bm{q}=c_{1}\bm{q}^{(1)}+c_{2}\bm{q}^{(2)}$, $c_{1},c_{2}\in\mathbbm{R}$ the expression $\eta_{j}[\bm{\vartheta}\times(\bm{\vartheta}\times\bm{Y}\bm{p})]_{j}$ can be written as a linear combination of the four expressions $\tilde{m}^{(k)}_{j}$, $k,j=1,2$:

[TABLE]

and similarly

[TABLE]

Decomposing $\bm{Y}\bm{p}=\left<\bm{\vartheta},\bm{Y}\bm{p}\right>\bm{\vartheta}+\bm{P}_{\bm{\vartheta}}\bm{Y}\bm{p}$, where $\bm{P}_{\bm{\vartheta}}\in\mathbbm{R}^{3\times 3}$ denotes the orthogonal projection in direction $\vartheta$, and using that

[TABLE]

the system of equations (34) considering (32) can be written in the form (33). $\square$

For \bm{Y}(\bm{v})=\bm{\tilde{\psi}}(\bm{v})+\mbox{\boldmath\chi}_{0}\bm{\mathcal{K}}[\bm{\tilde{\psi}}](\bm{v})+\bm{\mathcal{K}}^{\dagger}[\bm{\tilde{\psi}}](\bm{v})\mbox{\boldmath\chi}_{0}, where $\bm{v}=\tfrac{\omega}{c}(\bm{\vartheta}+\bm{e}_{3}),$ \mbox{\boldmath\vartheta}\in\mathbb{S}^{2}_{+} , Section 4 shows that the data $\tilde{m}^{(k)}_{j}(\bm{v})$ for $k,j=1,2$ and two different polarisation vectors $\bm{q}=\bm{e}_{1}$ and $\bm{q}=\bm{e}_{2}$ uniquely determine the projections $[\bm{P}_{\bm{\vartheta}}\bm{Y}\bm{p}^{(k)}]_{j}$ for $k,j\in\{1,2\}.$

Moreover, measurements for additional polarisations $\bm{q}$ do not provide any further informations so that at every detector point, corresponding to a direction $\bm{\vartheta}\in\mathbb{S}^{2}_{+}$, only the four elements $[\bm{P}_{\bm{\vartheta}}\bm{Y}\bm{p}^{(k)}]_{j}$, $k,j=1,2$, of the projection influence the measurements.

Remark \theremark:

*In contrast to standard OCT where three polarisation vectors where needed [13, Proposition 11] and to first order Born-approximation where $\bm{Y}=\bm{\tilde{\psi}},$ as we are going to see in the following, the above measurements due to the special form of $\bm{Y}$ allow for reconstructing all the unknowns functions $\psi_{ij}.$ *

Proposition \theproposition.

Let \mbox{\boldmath\vartheta}\in\mathbb{S}^{2}_{\ast}:=\{\mbox{\boldmath\mu}\in\mathbb{S}^{2}\mid\mu_{1}\neq\mu_{2},\,\mu_{3}>0\}. For two given incident polarisation vectors $\bm{q}^{(1)}$ and $\bm{q}^{(2)},$ the system of equations (33) is equivalent to a Fredholm type system of integral equations

[TABLE]

for some compact operator $\bm{\mathcal{C}}:(L^{2}(\Omega))^{3}\rightarrow(L^{2}(\mathbb{S}^{2}))^{3}$ and known right hand side $\bm{b}$ depending on the OCT data. Given the solution of (35), the component $\tilde{\psi}_{33}$ satisfies a Fredholm integral equation of the first kind

[TABLE]

where $\mathcal{C}:L^{2}(\Omega)\rightarrow L^{2}(\mathbb{S}^{2})$ is a compact operator and $b$ depends on the solution of (35).

Proof**:**

In order to reformulate equations (33), first we consider an arbitrary vector $\bm{p}$ and we split the expression $\bm{P}_{\bm{\vartheta}}\bm{Y}\bm{p}$ into the sum

[TABLE]

omitting for simplicity the $\bm{v}$ dependence of the unknown $\bm{\tilde{\psi}}.$

The first term on the right hand side admits the decomposition

[TABLE]

where we observe the independence on $\tilde{\psi}_{33}.$ To analyse the other two terms, we consider (29) and define the operators acting now on the components of the matrix-valued function $\bm{f}:$

[TABLE]

Since we are interested only in the first two components of $\bm{P}_{\bm{\vartheta}}\bm{Y}\bm{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

[TABLE]

where

[TABLE]

The only term where $\tilde{\psi}_{33}$ appears is the last one (as expected), namely

[TABLE]

where

[TABLE]

We can combine now all the above formulas to obtain

[TABLE]

where

[TABLE]

and

[TABLE]

Then, the system of equations (33), considering (32) reads

[TABLE]

We observe that in all equations the coefficient in front of the operator $\bm{\mathcal{M}}$ is the same, which is the only operator applying on the fourth component of $\bm{y}.$ In addition, from (38), we see that $\vartheta_{2}\bm{\mathcal{M}}_{14}=\vartheta_{1}\bm{\mathcal{M}}_{24}.$ Thus, in order to eliminate $y_{4}$ we reformulate the above system as follows: we subtract from equation (39a) the equation (39c), from equation (39b) the equation (39d) and from $\vartheta_{2}\cdot$(39a) the equation $\vartheta_{1}\cdot$(39b), resulting to

[TABLE]

The above system in compact form reads

[TABLE]

where

[TABLE]

and

[TABLE]

We compute the determinant of $\bm{\tilde{\mathcal{I}}}$ which is given by

[TABLE]

Recall that $\bm{\vartheta}\in\mathbb{S}^{2}_{+},$ meaning $\vartheta_{3}>0.$ Then, if in addition we impose that $\vartheta_{1}\neq\vartheta_{2}$ for all $\bm{\vartheta}\in\mathbb{S}^{2}_{+},$ the matrix $\bm{\tilde{\mathcal{I}}}$ is invertible with $\bm{\tilde{\mathcal{I}}}^{-1}=\det(\bm{\tilde{\mathcal{I}}})^{-1}\mbox{adj}(\bm{\tilde{\mathcal{I}}}).$ Then, equation (41) can be written in the form

[TABLE]

which is the Fredholm integral equation of the second kind (35), for $\bm{\mathcal{C}}:=\bm{\tilde{\mathcal{I}}}^{-1}\bm{\mathcal{N}},$ and $\bm{b}:=\bm{\tilde{\mathcal{I}}}^{-1}\bm{\tilde{b}}.$ Once (42) 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:

[TABLE]

for some known function $b,$ depending on $\bm{\tilde{y}}$ and $\bm{\tilde{b}}.$ This is equation (36) for $\mathcal{C}:=\mathcal{M}_{3}.$

To see the compactness of the integral operator $\bm{\mathcal{K}}$, we go back to the definition (29) and we consider the following decomposition:

[TABLE]

The above expression in compact form reads

[TABLE]

for the operators

[TABLE]

The operator $\bm{\mathcal{F}}:L^{2}(\Omega)\rightarrow L^{2}(\mathbb{S}^{2})$ is a modification of the usual far-field operator with smooth kernel thus compact. The operators $\mathcal{G}:L^{2}(\Omega)\rightarrow L^{2}(\Omega)$ and $\bm{\mathcal{G}}_{1}:(L^{2}(\Omega))^{3\times 3}\rightarrow(L^{2}(\Omega))^{3\times 3}$ are also compact due to their weakly singular kernels, see for instance [7, 25], and the operator $\bm{\mathcal{G}}_{0}:(L^{2}(\Omega))^{3\times 3}\rightarrow(L^{2}(\Omega))^{3\times 3}$ is bounded [9]. Thus $\bm{\mathcal{K}}:(L^{2}(\Omega))^{3\times 3}\rightarrow(L^{2}(\mathbb{S}^{2}))^{3\times 3}$ is also compact. The same arguments hold for $\bm{\mathcal{K}}^{\dagger}$ and then we can consider these properties also for the operators acting on the components of the matrix-valued function. $\square$

Remark \theremark:

*Equation (36) reflects the ill-posedness of the inverse problem, due to the compactness of the integral operator. *

5. 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, as a system of Fredholm integral equations (both of first and second kind). Under the assumptions of a non-dispersive, 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.

Acknowledgement

The work of OS has been supported by the Austrian Science Fund (FWF), Project P26687-N25 (Interdisciplinary Coupled Physics Imaging).