A Quantitative Model for Optical Coherence Tomography

Optical coherence tomography (OCT) is a widely used imaging technique in the micrometer regime, which gained accelerating interest in medical imaging in the last twenty years. In up-to-date OCT literature, certain simplifying assumptions are made for the reconstructions, but for many applications, a more realistic description of the OCT imaging process is of interest. In mathematical models, for example, the incident angle of light onto the sample is usually neglected or a plane wave description for the light–sample interaction in OCT is used, which ignores almost completely the occurring effects within an OCT measurement process. In this article, we make a first step to a quantitative model by considering the measured intensity as a combination of back-scattered Gaussian beams affected by the system. In contrast to the standard plane wave simplification, the presented model includes system relevant parameters, such as the position of the focus and the spot size of the incident laser beam, which allow a precise prediction of the OCT data. The accuracy of the proposed model—after calibration of all necessary system parameters—is illustrated by simulations and validated by a comparison with experimental data obtained from a 1300 nm swept-source OCT system.


Introduction
Optical coherence tomography (OCT) has proved to be a non-invasive, high-precision imaging technique with micrometer resolution.It emerged around 1990 for in-vivo imaging of the human eye [11,13] and gained increasing interest ever since.Nowadays, extensions like angiography [18], polarization sensitive OCT [3] and optical coherence elastography [25] unlocked a wide range of possible applications; for example, blood vessel analysis [19] and cancer margin detection [17], while OCT endoscopes [1] are clearing the way for high resolution imaging of internal organs and their pathologies.Multi-modal imaging techniques [22] often use OCT as morphological guidance.
While many theoretical OCT articles assume a sample geometry that is perfectly perpendicular to the OCT beam [2,9,21], commonly used OCT systems are designed for rough sample surfaces and arbitrary sample inclinations, which yield much less power at the detector.Normal incidence not only oversaturates the detector easily, especially for samples with a high refractive index and a very directed scattering profile, but can also lead to interference between the sample and optical parts inside the setup, e.g. the scan lens.To prevent imaging artifacts, normal incident is therefore usually avoided in OCT.
In addition, most works are based on a plane wave ansatz for describing the sample and reference fields [10,14,16,24].While this approximation is valid for the immediate focus region, it is most often not true for the whole field of view of the setup or even the whole sample area.Common effects like a focus dependent intensity profile inside the sample cannot be described with a plane wave ansatz.While workarounds like multiplying the spectral resolution with a sensitivity factor have been proposed [7], in this work a Gaussian beam, similar to [4], is used as incoming wave for a more precise description of the light beams.
To make a quantitative reconstruction of the optical parameters of the sample possible, we make the simplifying assumption that the sample is not absorbing and can at least locally be described as a layered medium (with layers not necessarily perfectly perpendicular to the incident light), which is a classical assumption in this field.This simplification allows us to analytically calculate the scattered light from the sample, which is then collected by a scan lens and combined with the reference light to produce the interference pattern.We roughly model the effect of the scan lens on the scattered light by discarding plane wave components moving in wrong directions.Together with the layers and the focusing effects introduced by Gaussian beams we derived simulations which have been in much better accordance with the experimental data obtained by a 1300 nm swept-source OCT system compared to a simple plane wave approach.
The paper is structured as follows: In Section 2 we describe the OCT imaging system that was built for this work.The individual experiments performed for investigating the different effects of the system on the data are also introduced.In Section 3, we present the general problem and the governing equations.We give the forms of the sample and reference fields using the near-and far-field representation of back-scattered Gaussian waves.The section ends with the derivation of the formula for the measurement data.In Section 4 we present in two experiments the dependence of the data on the incident angle and the beam focusing and show how to use them as calibration tools to determine otherwise unknown parameter such as the beam radius in the focus.The comparison between experimental and simulated data is presented in the last section.There, we see that our model nicely predicts the behavior of the experimental data with respect to different orientations and positions of the sample relative to the focus.

Experiments
Since there are many different variants of OCT systems around, we briefly describe the system we use for generating the data and the individual experiments.

OCT Setup and Post-Processing.
All measurements where performed with a custom-built fiberbased OCT system with a central frequency of about 1300 nm and 30 nm bandwidth, schematically shown in Figure 1.The core of the setup is the akinetic swept-source from Insight Photonic Solutions, USA, which emits about 60 mW at a repetition rate of up to 500 kHz.This swept-source shows a flat power profile over its whole bandwidth and a high phase stability, making it well-suited for any type of signal analysis.A fiber optics coupler guides 75 % of the laser light into the sample arm and 25 % into the reference arm, where a 4 mm fiber collimator releases it onto a short free-space path, with a moveable mirror at the end, that reflects the light back into the collimator.The custom-built sample arm features a rotatable imaging probe with conjugated scanning and a LSM54-1310 scan lens from Thorlabs, USA, for a flat imaging plane.Circulators are used to guide the reference and sample arm signal to a 50/50 % fiber coupler, where the laser light recombines.A dual-balance-detector (BPD-1, Insight Photonic Solutions, USA), short DBD, records the cross-correlation term and an ATS9360 data acquisition card from Alazar Technologies, Canada is used to digitalize it.
The Insight source supplies a trigger signal, which is used by a field programmable gate array to coordinate galvanometer movement and triggers the data acquisition card.Since the swept-source has a 100 % duty cycle, which the data acquisition card can not keep up with, every second sweep is neglected, to preserve the whole spectrum for imaging.An attenuation wheel in the reference arm free-space beam path is used to control the interference power detected by the DBD and thereby the intensity of the recorded OCT signal.It is used to ensure a high signal, without oversaturation of the detector.The system achieves an axial resolution of 31 µm and a lateral resolution of 24 µm in air as well as an SNR of 105 dB.The OCT control software was written in Labview and data processing was performed in MATLAB.
The recorded spectrum with 700 datapoints can immediately be Fourier transformed into image space, since the wavelength sweep emitted by the laser is already spaced equally in terms of wavenumber.The small bandwidth makes dispersion compensation unnecessary.Background removal is performed via subtraction of the average spectrum of a volume.Zeropadding ensures the small axial pixelsize of 13.7 µm in air and the lateral pixelsize is 9.8 µm in air.The DBD records only the difference between both input signals, thereby removing common-mode noise and centering the signal around zero.Through the digitization process the signal is shifted in height so it can be stored as unsigned integer.This shift is removed again during post-processing through the background subtraction.

Power vs. Angle.
To investigate the influence of the incident angle between sample surface and OCT beam, one of the fibers entering the DBD was connected to an optical power meter (PM100C with S122C, Thorlabs, USA) instead.A mirror was fixed onto a goniometer stage with a Vernier scale (GOH-65A100RUU, OptoSigma, USA) with a kinematic mount.First, the goniometer was aligned to ensure that the OCT beam goes through its Pivot point.Second, the kinematic mount and a vertical stage were adjusted to put the mirror in focus and ensure normal incident of the laser beam on the mirror, using the power meter as guidance.Then the mirror was tilted in 5 arc minute steps back and forth between −1°and 1°and the power was recorded until each angular position was measured 6 times.

Power vs. Focus.
For quantification of the Gaussian behavior of the focus a motorized stage (T-LSM050A, Zaber Technologies, Canada) was used to transport once a mirror and once a microscopy coverglass through the focus of the OCT system, with a roughly fixed tilt of about 2.75°.The coverglass (631-0124, VWR International, USA) has a refractive index of 1.5088 for 1300 nm, which needs to be taken into account during data analysis.At each position of the stage, we use 11 steps for the mirror and 7 for the coverglass, a 3D OCT volume was recorded.These 3D volumes were post-processed according to Section 2.1 and used to determine the exact incident angle and the distance of the center of the sample surface to the position where sample and reference arm would have the same length, called zero delay.

Mathematical Model
Considering the workflow of the used OCT system, described in the previous section, we model the parts shown in Figure 2 separately.
Firstly, in Section 3.1, we model the produced laser illumination.The laser light is split into two beams, one is sent to the sample and the other to the mirror in the reference arm.We model their scattering process in the Sections 3.2.1 and 3.2.5 respectively.
The reflected light in the sample arm is (partially) collected by a scan lens and coupled into a fiber.This aspect is the topic of Section 3.2.4.
After recombination of the scattered light beams, we model the detection via a dual-balancing detector of this superposition in Section 3.3.This in particular is discussed for the measurement related to two experiments explained in Section 2, which in the end are obligatory for the calibration of necessary parameters in the forward simulations.

Gaussian Beam illumination.
The shape of the light produced inside an optical resonator (we ignore at this point the finite size of the resonator and the boundary conditions) of a laser can according to [23], for example, be well described by a Gaussian beam.
We consider a Gaussian beam E : R 3 → C 3 as a monochromatic solution of the electromagnetic wave equation in vacuum which reduces it to Helmholtz equation (usually it is considered as solution of the paraxial approximation of the wave equation which is not necessary here): ( It is characterized by its form in the focal plane {x ∈ R 3 | x 3 = r 0 } for a function f : R 2 → C such that its 2D Fourier transform is compactly supported in D k0 (0) (the open ball with center 0 and radius k 0 ) and a polarization vector p ∈ R 2 × {0}.
R , Section 3.2.5 E (0) , Section 3.1 S , E S Section 3. Theorem 3.1 Let f : R 2 → C be a function such that its two-dimensional Fourier transform f is compactly supported in D k0 (0) and let p ∈ R 2 × {0}.Then for every x ∈ R 3 a solution of the Helmholtz problem (1) is given by with Such a wave describes well the light inside the optical resonator of the laser.Then through one partly transparent mirror of the resonator, we then obtain only the light moving in the negative x 3 −direction of the form Hereby, a reasonable model for the shape of the function f is one which resembles a Gaussian function.This laser light is transported within single-mode fibers through the OCT system, conserving the shape of the Gaussian beam throughout the system.

Backscattered Gaussian Fields. The laser light is split into two waves, E (0)
S for the sample and iE (0) R for the reference arm, as in (5) respectively, by a beam splitter and both remain in the form of a Gaussian beam, with possibly different beam parameters, for example due an optical attenuation wheel inside the reference arm, which causes a difference in light intensities between the beams.

The Sample Field
The beam E (0) S is now directed onto the sample and we say it is of form (5) with f = f S .Then, if the beam is sufficiently focused, meaning that the values of |E (0) S | can be neglected outside a small region, we only need to consider for the scattering process the shape of the sample inside this region.In this subregion, we denote it by Ω, we assume, using the tangent plane approximation, that it can be described by a layered structure.These layers are not necessarily perpendicular to incident beam, but for simplification assumed to parallel to each other.This is modeled by Ω being a finite union of layers: for some unit normal vector ν Ω .Each of these shall be characterized by a constant refractive index Under these conditions, we model the backscattered field E S as solution of Helmholtz equation where n(x) = L j=1 n j χ Ωj (x) + χ R 3 \Ω (x) and appropriate radiation conditions are assumed.The incident field ( 5) is represented as a superposition of plane waves having different wave vectors.Because of the linearity of the equation it is sufficient to solve the problem for every plane wave.The result for these backscattered fields for such a sample is well known in this plane wave case, see [15,8].For the simplest case L = 1, we consider an (arbitrary) plane wave as incident illumination from the top, with amplitude function α : R 2 → C 3 and propagation vector which we consider implicitly as a function of k 1 and k 2 .We obtain the reflected electric field where x Ω denotes an arbitrary point of the top boundary (that is x Ω , ν Ω = a 1 ) of the object, the wave vector and β S the sum of the reflection coefficients β 0, , β 0,⊥ of the differently polarized parts Here, we have decomposed α into its transverse electric and magnetic polarizations, with coefficients λ 1 and λ 2 , respectively.Further, we use Snell's law for the determination of the transmission angle θ t .
Summarizing the scattered (plane) waves for all (k 1 , k 2 ) and then finally results in with β given by (10).

Far Field Method
Since the distance between the scan lens and the sample (which is roughly 6 cm) is much greater than the size of of the sample itself (which is only a few millimeters), we can be tempted to simplify the integrand by using the far-field approximation.Mathematically, this means, that we are approximating (11) by its behavior at some point rs, s ∈ S 2 , as r → ∞ : To compute the dominating term E S,∞ , we apply the method of stationary phase, see Lemma 8.1, which is based on the approximation of the phase function k 0 Ψ, with , by its Taylor series around its critical points.Theorem 3.2 Let E S be a vector field given by (11).Then, its far field approximation takes the form where for c 3 = s 3 − 2 ν Ω , s ν Ω,3 , the vector components k 1 and k 2 are given by and for k 1 , k 2 the reflected vector k r = Φ(k) is given by (9).

Comparing the Near-and the Far-Fields
Comparing the representations ( 11) and ( 12) of the scattered electric field simulations make it obvious that there is a difference in the visible effects provided by these methods.In this work are concerned with the influence of the focus in the scattered field.
In contrast to the far-field representation, the scattered field in the near field regime is heavily depending on the distance between the positions of the layer and of the focus.We neglect the vectorial quantities in (11) for a moment and allow for an amplitude function where the parameter a > 0 is such that the error | f − f χ Dρ 0 (0) |, is negligible.Hereby, D ρ0 (0) is a disk with small radius ρ 0 and center 0.Then, for a single surface (medium) we obtain the scattered field We assume that on the small disk, the reflection coefficient β 0 is approximately constant and due to small deviations of k 1 , k 2 from zero we may approximate the root in the exponents Further, for the sake of simplification, we restrict ν Ω ∈ R 2 × {0}, ν Ω,2 = 0, fix the positions of the focus r 0 and the object x Ω = x Ω,3 , x Ω,3 < 0 below the origin and evaluate at x = 0.This then finally gives for the scattered field where we defined the phase elements ψ 0 , ψ 1 , ψ 2 by Since e(ψ 2 ) = a > 0, we evaluate both integrals in ( 16) and arrive at After complex conjugation in the exponent and taking the absolute value of the field, we find that Considering now (18) for different positions x Ω,3 of the sample, and therefore for varying d, corresponds to different evaluation points x ∞ = re 3 , with r = |x Ω,3 | in the far-field regime.Taking the absolute value of ( 12) we observe that opposed to the near-field representation, the far-field regime is independent of the focus position.Figure 3 provides a comparison between both for different positions of the focus and the surface.Since the far-field approximation does not show the dependence of the electric field on the focus, we stick with the more complicated near-field representation of the scattered light E S in (11).x Ω,3 [m]

The Scan Lens
The backreflected light E S then passes trough the scan lens and is collected by a fiber collimator.Thereby, we loose all plane wave components whose propagation directions are outside a certain angular range of the collimator.We model this by reducing the area of integration in (11) to a set B of those scattered wave directions k r which have an angle to the measurement direction e 3 less than a certain angle of acceptance θ max , that is which finally gives a (scattered) sample field at the scan lens position x D = r D e 3 above the sample.

The Reference Field
Similarly, we model the reference field as solution to the scattering problem (6) with Gaussian incident illumination R and a medium with constant (infinitely) large refractive index.Following the same line that led to (11), we get with f = f R , a field of the form Following the experimental setup the mirror in the reference arm is perpendicular to the incident light, so that the unit normal vector ν M = e 3 , and positioned in the focus of the E (0) R , such that, following Section 3.2.3, the far-field approximation for reference field E R is valid.We thus have a reference field

OCT Measurements.
With the identities of the incident and the backscattered light in hand, we now proceed with the modeling of the detection process inside an OCT system.In order to suppress common-mode noise and enhance the signal-to-noise ratio, a dual-balance-detector is used to record the OCT signal.After recombination of sample and reference arm light, the laser signal is split 50/50% into two different fibers F 1 , F 2 , each entering one of the DBDs optical inputs.The DBD then subtracts one input from the other, thereby removing everything but the cross-correlation term of the interference.Thus, assuming that the sample and the reference fields E (1)

S and iE
(1) R are passing through a perfect splitter, we obtain the forms for the fields in the fibers as We assume, ignoring the travel paths inside the fibers, that these fields are detected at the position x D of the scan lens.These measurements are performed for different wavenumbers k 0 in a scan range [k min , k max ].We therefore indicate explicitly the dependence on k 0 in the measurements: S , E With the identities ( 20) and ( 21) for E (1) S and E (1) R , we obtain where we define the phase function and use B as given in (19).

Calibration of the Forward Model
So far we have investigated both the modeling of the backscattered wave from a (layered) object under Gaussian laser illumination and the measurement process of the OCT system described in Section 2.1.However, simulations based on this explicit model and the following comparison with experimental data presupposes the knowledge of a list of system parameters.Within this list we distinguish between parameters with values known from specifications such as the wavenumber k 0 and parameters which we need to calibrate from the experiment, the beam radius of the Gaussian beam and the angle of acceptance, for example.In order to be capable of extracting these parameters for the simulations, we use two calibration experiments.On the one hand, we consider an experiment which shows the behavior of the backreflected laser power for different surface tilting angles and on the other hand, we consider the influence of varying positions of the object, with respect to the focus, on the measured data.

Angular Dependence of the Measured Power.
We use a mirror as a sample and analyze the influence of the surface angle on the measured intensity of the scattered electric field.Following the measurement process in Sections 2.1 and 2.2, the reference arm is blocked, preventing any light from the reference arm to reach the detector.Furthermore, one of the two fibers, which would normally enter the DBD, is connected to a power meter.The measured data is therefore given as the intensity of the scattered field of this mirror.Hereby, again referring to the experimental setup in Section 2.2, we model the totally reflecting mirror as a sample characterized by an infinitely large refractive index.We parametrize the unit normal vector of the mirror surface by for small values of We describe the measurement process for this experiment in a way, that the scattered light is detected by a single scan lens point, for simplification we say x D = 0, for a selected wavenumber k 0 in the spectrum [k min , k max ].This in the end, yields a measured intensity of the form where E S is given by ( 20) and τ ∈ C accounts for the traveling through the beam splitters.Additionally, we say that the function fS is approximately given as in (14) with a = w 2 0 /4, where w 0 represents the radius of the Gaussian beam at the focus.Following the experimental setup we fix the location r 0 < 0 (below the detector) of the focus and the mirror x Ω,3 and assume that they are equal: r 0 = x Ω,3 .We follow the notation from Section 3.2.3,but approximate this time the exact form of the domain of integration B defined in (19), which is an ellipse, by the rectangular domain with the parameters and assume that this characterization of B still allows for an approximation of directions as in (15).Then, using the definitions of ψ j for j ∈ {0, 1, 2} in (17), we obtain the intensity of the scattered field as a function of θ Ω , w 0 and θ max τ E (1) where si : R → R denotes the unnormalized sinc function given by si(x) = sin(x) x and erfi : C → C is the imaginary error function, defined by erfi(z) = 2 √ π z 0 e ζ 2 dζ.Thus, from measurements M 1 (θ Ω ) as in (26), corresponding to the data provided by our power meter, for different values θ Ω ∈ [θ Ω , θ Ω ], we can extract the beam radius w 0 at the focus and the angle of acceptance θ max as solutions of the minimization problem with the function G given by (27).

Reconstructing Sample Information from an OCT Experiment.
In the previous section the beam radius w 0 at the focus and the angle of acceptance θ max have been found.In order to complete the set of parameters necessary for the reconstruction from a measurement at a single point, we additionally need the normal vector ν Ω of the tangential plane at each layer boundary.
In swept source OCT one in-depth profile of the sample, that is a measurement of the form of M in (22) (in this case centered at x 1 = x 2 = 0), called an A-scan, is acquired during one wavenumber sweep of the laser.To get 3D information, raster scanning in x 1 and x 2 direction over a certain field of view is performed.The data used in the following is considered as a B-scan, a line of A-scans where only x 1 varies at a fixed position x 2 .Since we assume our layer boundaries to be planes with a certain normal vector ν Ω , the surface points fulfill an equation of the form x Ω , ν Ω = c.If we can therefore determine at every raster position x Ω,1 , x Ω,2 the third component x Ω,3 , this determines the normal direction ν Ω .
Since the single A-scans along those lines are performed independently, we treat these A-scan as single measurements.We shift the coordinate system always so that the incident beam is located at x 1 = x 2 = 0 and therefore have x Ω = x Ω,3 e 3 .We recall, that the mirror in the reference arm is modeled as a medium described by an infinitely large refractive index with unit normal vector ν M = e 3 and fixed position at Under these assumptions, we rewrite ( 23) and ( 24) as For fixed mirror position x M , focus r 0 and detector r D , ψ only varies with respect to different depth positions x Ω,3 of the the sample.Thus, if we can determine the function ψ from the measurements M 2 for different A-scans, we also obtain the depth information about the sample.
Under the simplifying assumption that the far-field approximation of the scattered sample field, using s = e 3 in Theorem 3.2, is a reasonable approximation in this case, we rewrite (28) as where the point (k 1 , k 2 ) is defined by (13).Since ψ depends linearly on k 0 , the measurements are then given as a harmonic oscillation with respect to k 0 and with frequency ψ/k 0 .To solve for this frequency, we want to Fourier transform with respect to k 0 , which we define by However, since we only have band-limited data, we will study the function We will show in the following that ψ/k 0 is determined as the argument where the maximum is located, that is, ψ/k 0 = argmax κ I(κ).(The absolute value is used to avoid real-or imaginary parts with higher frequent oscillations in order to stably calculate a maximal point.) We assume that fS and fR in (29) are of exponential form as in ( 14) and define the measurement function with the parameters and interchanging the integral and differentiation in the Fourier transform F ∂ σ 2 e −k 2 0 σ 2 , we find a form for the Fourier integral in (30) as a convolution for k = kmax+kmin 2 and δ = kmax−kmin 2 .To simplify this expression, we introduce the values and observe from Table 1, that σk and σ are of the same order and σ δ is considerably larger compared to both of them, meaning that σk = Qσ, σ σ δ , (33) for some Q ∈ R which is close to one.Writing the functions under the integral (31) in terms of these values gives us with the expression Considering (33), we will expand this around σ = 0.
Lemma 4.1 Let σk, σ δ , σ as in (32) satisfying (33).Further, let f σk,σ be defined as in (34).Then, we have for small values of σ the approximation Thus, by applying Lemma 4.1 to (35), we obtain after changing back to the original system of coordinates Note that the dominant sinc terms are centered symmetrically with respect to the origin.In order to derive an explicit expression for the maximum of (37), we want to assume that Θ 0 is far away from the origin (which can be accomplished experimentally by tuning the position of the sample) then these sinc functions do not influence each other strongly.We shift one of them to the origin by setting κ = κ − Θ 0 and obtain Lemma 4.2 Let F be defined by (38).Then, for Θ 0 → ∞, the function F attains a local maximum at κ = 0.
Shifting back to the original coordinates and using Lemma 4.2 yields that (37) attains a maximum at κ = Θ 0 , that is Θ 0 = argmax κ I(κ), which finally gives a representation of (30) as Thus, from the definition of Θ 0 we can uniquely determine ψ.
We use this information for the reconstruction of the surface angle θ Ω .For two different, but known lateral positions x j Ω,1 , j ∈ {1, 2}, we consider A-scans leading to measurements M j 2 of the form (28), for different depth positions x j Ω,3 , for j ∈ {1, 2}, respectively.By using the above analysis (under the assumption that the far-field approximation is valid), we determine from the Fourier transform of these two the phase contribution ψ in dependence of x 1 Ω,3 and x 2 Ω,3 .Under the assumption that θ Ω is considered small, the subtraction of these two then leads to ). Together with known lateral information and using that the unit normal vector on the surface satisfies ν Ω , (x 2 Ω − x 1 Ω ) = 0, we determine θ Ω as .

Results
Finally, we want to validate our model by comparing the simulations with experimental data.We focus on the two previously addressed experiments (see Sections 4.1 and 4.2, respectively 2.2 and 2.3).This quantitative approach shows the dependence of the data on the surface tilting and the focus position.First we will use the calibration measurement to calculate the beam radius at the focus and the angle of acceptance, then we look (using the just calibrated parameters) at the data from Section 2.3.The main parameters are presented in Table 1.

Power vs. Angle -Experiment.
Following Section 4.1, we first calibrate the beam radius and the angle of acceptance from the experimental data for different values of the surface angle.This procedure is presented in Algorithm 1.The algorithm is based on the approximated form (27) shortening the computation time by evaluating a one-dimensional integral, instead of the two-dimensional integration presented in (20).As discussed in Section 2.2, the sample arm power arriving at one of the DBD entrances was measured M = 6 times at J = 25 angular positions.While the laser power behaves very stable, due to the expected error from the rotational stage and the Gaussian dependence from the angle, some error is observed in the power vs. angle data.Thus, in the following, the data will be plotted with errorbars, representing the standard deviation.
The simulated data M 1 , see (26), is given for θ Ω = θ j , j = 1, ..., J.The match between experimental and simulated data is presented in Figure 4. We remark that both data and simulation follow a Gaussian behavior and attain the maximum at normal incidence, as expected.

Power vs. Focus -Experiment.
Using the calibrated spot size and angle of acceptance from the previous experiment, we compare the simulations with experimental data for multiple B-scans of a mirror and a coverglass as samples of interest.Following Section 4.2, we first determine the surface angle θ Ω and adapt the integration area in (19).The coverglass, which is described by a medium with constant refractive index n 1 = 1.5088 and perfectly parallel surfaces, has a thickness d, which is also determined from the experimental data, see Algorithm 2.
We ignore polarization effects in the following and use the form (20) with β S = 1 for the simulations of the scattered field of the mirror data.However, for the coverglass experiment we extend the form to a

Parameter
Value Unit λ 0 central wavelength 1300 nm w 0 beam radius at focus 14.15 µm θ max angle of acceptance 1.5709 °N A 0.037

µm
Table 1.List of parameters: The central wavelength and the wavenumbers k min and kmax are determined through the specifications of the used swept-source.The beam radius and the angle of acceptance were found through calibration (see Section 5.1).The numerical aperture NA was calculated from the angle of acceptance.The parameter σ is given for a typical tilting angle θ Ω in our experiments.
layer model with two parallel surfaces where we have given the reflection coefficient We assume again that fS can be well approximated by (14) and define B as in (19).
2n1 cos θt ; Algorithm 2: Extraction scheme for the tilting angle θ Ω and the thickness d from the power vs. focus experiment.
Although the dual-balance detection already lowers the noise level in the data, we need to minimize the effects of the residual noise for a meaningful comparison with the simulations.
Thus, for every individual B-scan, we use the mean value of all maximum intensities of it's A-scans, i.e. we consider the map 1 Ñ Ñ n=1 max κ∈R I n (κ), where we use I n = I, as defined in (30), with n ∈ {1, . . ., Ñ } accounting for the number of A-scans in every B-scan and for the different position x n Ω,3 (corresponding to this A-scan), as a reference for our simulations.The errorbars in Figure 5 represent the standard deviation of these variations over each B-scan.The highest value of the experimental data is matched to the highest value of the simulation for comparison.
Figure 5 shows, that the Gaussian behavior of the data for the mirror experiment follows the theory.
In Figure 6, we see the comparison between experimental data and simulations for both boundaries of the coverglass.Unfortunately, the calibrated parameters (w 0 , θ max ) from the previous subsection do not yield optimal results, see Figure 6.Similar to Figure 5 a sufficiently strong decrease of the averaged maximum values away from the focus position can be identified in the simulations for the coverglass experiment as well.At this point, we remark by comparing the experimental data sets, see Figure 7, that the Gaussian curve for the coverglass experiment shows a slightly stretched behavior.This is explained by the measurement of the returning laser light in a diffusive regime originating from the reflection at a slightly rough coverglass boundary surface.
However, updated parameters can be found using an experiment similar to the calibration of w 0 and θ max , see Figure 8.In contrast to the calibration, the power measurement for the coverglass includes information of both boundaries and therefore the measured field intensity is provided as a sum where we consider first order reflections only.Due to additional scattering events inside the coverglass material, the background information E (1) S,1 and therefore neglected for the calibration.A comparison for the coverglass experiment -after updating the system parameters -shows that almost all simulated data points lie inside the estimated range for both boundaries, see Figure 9.

Discussion
We have considered samples with a very distinct scattering profile, but still a slight difference in the angular reflectivity profile could be observed for the mirror and the coverglass.For diffusely scattering samples the proposed description can easily be generalized, especially with the mentioned automatic angular power measurement.An automated measurement would also reduce the error, so that a fast angle scanning procedure could be implemented prior to OCT imaging.
The error in the power vs. focus experiment (see Figures 5 and 6) is caused by power variations inside single B-scans.These are explained by a combination of reasons.The tilt of the sample with respect to the OCT illumination, generates a slight continuous change in distance to the focus for each A-scan inside a B-scan.In addition the scan lens induces a certain curvature of field, resulting in a change of illumination depending on the raster scanning position x 1 .
Although the presented results in this work show suitable correspondence with the provided OCT data, we note that the algorithm for solving the minimal-error-solution problem in extracting the beam radius and the angle of acceptance from the power-angle experiment suffers from the fact that the function in (27) is highly oscillating and it is therefore difficult to find an "optimal" set of parameters.
Nevertheless, for a rather large range of values of these parameters, we get a reasonable match with the experimental data, at least for simple, layered examples.The model should, however, work nicely also for more complicated samples with different geometries (that is, with more and potentially curved layers).We expect that this model can be used as a forward model of the inverse problem of reconstructing the refractive index inside of the layers.

Conclusion
We presented a method to model the image formation in OCT based on a real-life 1300 nm swept-source setup.In contrast to publications based on plane wave models for the OCT system, the proposed model includes the effect of additional system relevant parameters such as the focus and the beam radius of the incident laser light and the angle of acceptance.We also suggested a way how to determine these (not necessarily a priori known) parameters, either from the OCT data itself and from calibration measurements.A comparison between simulation and experiment shows, that the presented model, together with the derived system and sample parameters, produces a quantitatively correct prediction of the OCT data.We therefore expect that this model can serve as a forward model for an inversion algorithm to quantitatively reconstruct from OCT data the optical material properties of the sample, in particular its refractive index.

Appendix
Here, we collect the proofs of the theorems.We start with the one of Theorem 3.1 which describes the Fourier decomposition of a solution of the Helmholtz equation.

Proof (Proof of Theorem 3.1):
To simplify the notation, we shift the coordinate system x 1 = x 1 , x 2 = x 2 , x 3 = r 0 − x 3 such that focal plane is in the origin of the new system.Then, we apply the Fourier transform in the first equation of (1) with respect to the x 1 , x 2 -components, resulting in an ordinary differential equation x 3 ) = 0, for the first two components of the electric field.We know that for j ∈ {1, 2}, is a solution of this problem.Using the Fourier transformed initial data at the plane x 3 = 0, we find that the coefficients α −,j , α +,j are given by So far we have seen that E 1 , E 2 are solutions of the Helmholtz equation, without considering the third component of E. Finally, we use that E is divergence-free to find that Moreover, taking two times the derivative with respect with x 3 , we find  which in the end, using also the second derivatives with respect to x 1 and x 2 , gives Thus, E 3 is also solution of the Helmholtz equation if and only if This is equivalent to the condition Since this condition must hold true for every (k Given the representations of E j , j ∈ {1, 2}, we derive a representation also for the third component of the electric field Finally, we use the original coordinate system and we obtain the desired representations (3) and (4).
Next, we come to the derivation of the far-field representation of the scattered field presented in Theorem 3.2.This is described, for example, in the book [20] and is based on the stationary phase method.We can now apply this stationary phase method, Lemma 8.1, to the integral in (11).
Finally, we come to the derivation of the asymptotic behavior of the intensity of the maxima in the Fourier transform of the OCT signal for small incident angle.
Proof (Proof of Lemma 4.1): Since the analysis of the integral in (36) proceeds along the same lines for g σ δ ,+ and g σ δ ,− , we simply write g σ δ = g σ δ ,± to make the notation easier.
We assume that locally around ζ 0 = 0, we can write g σ δ as its Taylor series Using this in (36), gives • For j = 1, we find in the same way • Similarly, for j = 2, we get which is constant with respect to σ.
• Following the same procedure, the remaining integrals for j ≥ 3 are of the form R for a given pre-factor C(j).
The assumption that σ is small, let us say σ 1, yields that the terms of order σ −2 dominate.Keeping these terms only, results in Then, it is clear that lim Θ0→∞ F (Θ 0 ; κ ) = si(δκ ) 2 , which attains its maximum at κ = 0.

2 . 1 EFigure 2 .
Figure 2. Modeling of the separate parts of the OCT experiment: We start by describing the light beam produced by the laser (the red box) in Section 3.1, we give a representation for the beam in the sample arm which is backscattered from the sample (the orange box) in Section 3.2.1 and is then coupled back into the fiber system via the scan lens (the yellow box) in Section 3.2.4.This is afterwards recombined with the beam from the reference arm (the green box) in Section 3.2.5 and detected by the dual balance detector (the blue box) in Section 3.3.

Figure 3 .
Figure 3.The near-field for different positions of the focus (red, blue, green) vs. the far-field (black) regime for different positions of the surface.At the dotted lines, indicating where the surface and focus position coincide, the intensities of near-and far-field regime are equal.

Figure 4 .
Figure 4. Comparison between the power meter measurements for different angular steps of the mirror (blue dashed curve) as in Section 2.2 and the simulation (red curve) for (27).

Figure 5 .
Figure 5.Comparison of averaged maximum values for all B-scans (of different sample locations) of the experimental data (black with errorbars) and the simulated data points (red) for the mirror experiment.

Figure 6 .
Figure 6.Comparison between the experimental data (black with errorbars) and the simulated (red) data points for calibrated values of w 0 and θmax.Above we see the top boundary surface of the coverglass, below the background surface.

Figure 7 .
Figure 7.Comparison between averaged maximum values for all B-scans of the mirror (red) and the coverglass (blue) experiment.

Figure 8 .
Figure 8. Angular scattering profile of a coverglass: comparing experimental data (blue dashed line) with simulations (red).

Figure 9 .
Figure 9.Comparison between averaged maximum values of the experimental data (black) and simulations (red) for different position through the focus after the recalibration of values w 0 and θmax.Above we see the top surface, below the background surface of the underlying coverglass.