## Abstract

Precise control of the temperature rise is a prerequisite for proper photothermal therapy. In retinal laser therapy, the heat deposition is primarily governed by the melanin concentration, which can significantly vary across the retina and from patient to patient. In this work, we present a method for determining the optical and thermal properties of layered materials, directly applicable to the retina, using low-energy laser heating and phase-resolved optical coherence tomography (pOCT). The method is demonstrated on a polymer-based tissue phantom heated with a laser pulse focused onto an absorbing layer buried below the phantom’s surface. Using a line-scan spectral-domain pOCT, optical path length changes induced by the thermal expansion were extracted from sequential B-scans. The material properties were then determined by matching the optical path length changes to a thermo-mechanical model developed for fast computation. This method determined the absorption coefficient with a precision of 2.5% and the temperature rise with a precision of about 0.2°C from a single laser exposure, while the peak did not exceed 8°C during 1 ms pulse, which is well within the tissue safety range and significantly more precise than other methods.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

## 1. Introduction

Controlling the temperature rise in tissue during laser therapy is fundamental to many medical procedures, from thermal ablation of tumors to dermatological and retinal treatments. Retinal photocoagulation, for instance, has been a standard procedure for decades to treat a number of retinal pathologies, such as proliferative diabetic retinopathy [1], macular edema [2], central serous chorioretinopathy [3], and others [4]. Extent of heating and the induced therapeutic effect depend on the temperature course (magnitude and duration), which, in turn, depends on the laser parameters (wavelength, power, duration, spot size, etc.) and the tissue properties (optical absorption and scattering, heat capacity, and conductivity). Recent studies have suggested that the therapeutic window, typically quantified by the Arrhenius integral Ω, of the non-damaging hyperthermia is narrower (0.1 < Ω < 1) than in conventional photocoagulation [5–7], and thus such treatments require higher precision in the temperature control.

In the retina, light absorption is primarily governed by the pigmentation in the retinal pigment epithelium (RPE) and choroid, associated with an optical absorption that is proportional to the melanin concentration [8]. Measuring the optical absorption in the retina (and in most tissues) *in vivo* is, however, particularly challenging since there is no direct optical access to the transmitted light, and it must therefore be done in reflection. Furthermore, the melanin concentration can vary locally within the retina and from patient to patient by a factor of two to four [9]. These variations, if not accounted for, can result in significant temperature differences and thus undesired lesions or sub-therapeutic exposures. To ensure the laser exposure be within the non-damaging therapeutic range (0.1 < Ω < 1), one can estimate using previously published finite-element thermal models of the retina [8,10] that the absorption coefficient must be determined with accuracy no worse than ±20% (see the Non-damaging treatment range section in the Supplemental Document and corresponding Dataset 1 [11]).

To assess the retinal absorption in reflection geometry, one can measure the tissue response to laser heating. A few approaches have been developed along these lines. One of them, an optoacoustic method, is based on the detection of the pressure waves at the cornea, generated by the expansion of melanosomes during repeated nanosecond laser exposures [12–15]. During the thermal and acoustic confinement, the detected pressure is proportional to the absorbed energy density, related to the absorption coefficient, and the Grüneisen parameter [16]. Precision of the temperature measurement in this approach has been limited to a few degrees. Further, the signal is sensitive to the position of the piezo sensor on the cornea and hence the measurements require a careful calibration to convert the measured pressure waves into temperature. Another approach is based on detecting the changes in back-scattered light due to microbubbles formed by vaporization, which lead to the RPE cell death [17]. However, due to associated tissue damage, this technique is not suitable for non-damaging therapies. Optical coherence tomography (OCT) has also being explored as a non-contact, high spatial (3D) resolution imaging method for *in-situ* monitoring of the tissue deformation [10,18–21]. In particular, phase-sensitive OCT [22] allows tracking of the nanometer-scale optical path length changes, and has been proposed for precise temperature control during photothermal therapy [18]. Yet, to date, experimental investigations have not directly and quantitatively related the phase signal to thermo-mechanical changes in tissues, from which the temperature could be derived. Instead, temperature profiles have been separately estimated via the optoacoustic method mentioned above [18] or analytically, using an *a priori* value for the absorption coefficient [19,20]. The calculated temperature profiles were then qualitatively compared with the phase profiles recorded via OCT, without any direct validation. Recent advances in pOCT demonstrated capabilities of detecting nanometer-scale deformations with sub-ms temporal resolution [23,24], which should enable rapid and precise optical thermometry.

In this study, we present a methodology to precisely determine the optical (absorption coefficient) and thermal (heat conductivity and coefficient of thermal expansion) properties of layered materials using line-field phase-resolved optical coherence tomography (pOCT). In this approach, we monitored the vertical displacements and the changes in the optical path lengths (ΔOPLs) of reflective and scattering layers with nanometer precision in a non-damaging regime. We developed an analytical model to the coupled thermo-elastic problem and solved it in the Hankel-Laplace transformed domain. This approach significantly increases the computation speed compared to finite element simulations, and therefore allows faster fitting of the material properties to match the spatial and temporal aspects of the experimental ΔOPLs.

## 2. Materials and methods

#### 2.1 Optical setup

The optical layout (Fig. 1) followed the design by Pandiyan et al. [25] for high-speed line-scan OCT imaging. We here summarize the main components of the design; further details on the optical system can be found in Ref. [25]. The illumination is provided by a fiber-coupled super-luminescent diode (M-T-850-HP-I, Superlum, 9-mW power) with a center wavelength of 840 nm and a bandwidth of 130 nm, providing an axial resolution of 2.4 μm in air. After collimation with a reflective collimator (RC) to a diameter of 1.7 mm (FWHM of the Gaussian intensity profile), the beam is shaped using a cylindrical lens (CL1) to yield a line-illumination at the tissue phantom, which is held around its periphery and whose front and back surfaces are in contact with air.

After CL1, the beam is split between the reference and sample arms with a 30:70 (R:T) beam splitter. In the sample arm, the beam, after going through magnification telescopes, is focused along the y direction to yield a line field of 0.8-mm in length (FWHM of the Gaussian intensity profile) in the x direction for an approximated focal spot size of 10 μm (beam waist diameter) and a Rayleigh range of 95 μm. A galvo-scanner (8310K Series Galvanometer, Cambridge Technology) is used to steer the beam laterally (*y* direction) at the sample surface and a deformable mirror (DM69-15, ALPAO) is used for fine focus adjustment and for the future adaptive optics implementation. The entrance pupil plane (P), the galvo-scanner, and the deformable mirror are optically conjugated using afocal telescopes (L1 + L2, L3 + L4).

The detection path, shown in two views in Fig. 1(b),(c) follows an anamorphic imaging configuration by using two cylindrical lenses to allow independent control of magnification in the *x* and *y* dimensions, and thus independent control of the spatial and spectral resolutions. The OCT spectrometer is constructed using a 600 line/mm grating (Wasatch Photonics) and a high-speed camera (Phantom v641, 10-μm pixel size) for a resolution of 512${\times} $ 768 pixels (spatial ${\times} $ spectral). B-scans, with a static illumination field, are recorded for 500 ms with a frame rate of 10 kHz and 98-µs exposure times for a full imaging field of view of 2.3 mm laterally and 1.0 mm in depth, one camera pixel corresponding to 4.5 µm laterally and 2.7 µm in depth. Flat glass windows are also added to the reference arm (not shown in schematic) for coarse chromatic dispersion compensation between the reference and sample arms.

For laser heating, a blue diode beam (450-nm wavelength, 1-ms duration, 260-mW peak power) is coupled into a 200-µm multimode fiber. The fiber output is imaged onto the absorbing layer with a 400-µm diameter spot of uniform intensity. An iris is introduced at an intermediate image plane to sharpen the beam image edge and adjust the spot size. The beam is coupled into the OCT path using a long pass dichroic mirror and manually steered so that the heating beam diameter overlaps with the line OCT illumination. The pulse duration and the delay for the OCT camera acquisition are controlled by a pulse/delay generator.

#### 2.2 Tissue phantom fabrication and materials properties

A five-layer tissue phantom, emulating a simplified retina, was fabricated by spin coating several layers of poly(dimethysiloxane) (PDMS) to different thicknesses on top of a glass substrate. PDMS is a common tissue simulant material and has been used for retina phantom assemblies [26]. PDMS layers were fabricated using two-part, 10:1 mixing ratio, Sylgard 184 kits from Dow Corning, on top of a 1-inch borosilicate glass substrate (N°2 coverslip, ChemGlass) (Layer 5). The spin-coating recipes are given in Table 1. Layers were assembled in reverse layer label order, Layer 1, the top surface, being the last layer coated (see Fig. 1(a) inset). After spin-coating, samples were degassed in a vacuum desiccator for 30 minutes. All layers were then cured for 48 hours at room temperature before the next coating step. Curing at room temperature was chosen so not to degrade the dye, which visibly deteriorated at temperatures above 40°C. Since the elastic modulus is dependent on the curing temperature [27], the same curing temperature was applied to all layers to simplify the mechanical model. To make the top, scattering PDMS layer (Layer 1), titanium particles (2.5wt%, 21-nm nominal diameter, Sigma Aldrich) were added to the base PDMS, previously diluted in hexane (3:1 weight ratio 3:1) to reduce the viscosity and hence facilitate the nano-powder mixing. The hexane was left to fully evaporate during stirring and regular sonication of the mixture for 24 hours prior to mixing with the curing agent and spin-coating. Layer 3, which aims to simulate the thin scattering and absorbing RPE layer, was fabricated using a PDMS/dye/titanium oxide mixture. First, a saturated solution of hexane with Oil Red O dye, strong absorber at 450 nm, was prepared. The solution was then mixed with the PDMS base (1:1 weight ratio) to obtain a low viscosity mixture. Lastly, the titanium oxide nano-powder was added (2.5 wt% with respect to PDMS base) to the mixture, which was left under stirring and regular sonication for 24 hours. The mixture was then spin-coated at high speed.

To measure the absorption coefficient of the dye-doped PDMS layer, a similar mixture but without the titanium particles (to avoid scattering), was prepared and spin-coated onto a glass substrate following the same recipe as for Layer 3. Light absorption in this layer, measured using the same 450-nm diode laser as in the heating experiments, was 8%. To calculate the absorption coefficient *μ _{a}*, the thickness of that layer was measured by a profilometer. For simplicity the mechanical and thermal properties of all PDMS layers were considered identical. All material properties used in the modeling and references are listed in Table 2.

#### 2.3 OCT signal processing

OCT image reconstruction included the standard steps of background subtraction (of the sample and reference arms), wavenumber resampling, and time referencing to cancel the local arbitrary phase [25]. The remaining dispersion, not compensated by the glass windows, is measured by placing a mirror instead of the tissue phantom. It is then numerically compensated during the OCT reconstruction [31]. We measured a system sensitivity of 75 dB by placing a mirror at the image plane in the sample arm. The sensitivity is comparable to the system developed by Pandiyan et al. [25].

In polar coordinates (*r*,*z*), the spatially static B-scan OCT signal can be described as a complex number:

The phase difference $\mathrm{\Delta }\varphi $ between two pixel coordinates (*r,z*)_{1} and (*r,z*)_{2} as a function of time is calculated by taking the argument of the product of the complex and complex conjugate numbers of the OCT signals at these locations:

The OCT intensity image of the top 100 μm of the phantom is shown in Fig. 2(a). The signal-to-noise ratio (SNR), which is directly related to phase sensitivity [32], is about 20–30 dB for the scattering layers of the tissue phantom. Because the top surface is reflective, while the other layers are scattering, it shows the highest SNR (maximum value ∼30 dB) and therefore would yield the lowest phase noise. The change in phase, Δ*φ*, can be converted to change in single-pass optical path length, ΔOPL, by:

*λ*is the central wavelength of the OCT illumination. The

_{OCT}*z*-coordinate of the top surface of the tissue phantom as a function of radius was found by tracking the peak OCT intensity averaged over 100 B-scan acquisitions, thus defining the plane A coordinates. Plane B and C were defined by translating plane A coordinates by 25 and 60 μm in depth, respectively.

#### 2.4 Solutions to the axisymmetric thermo-mechanical problem

The transient response of the tissue phantom to laser heating can be solved using a thermo-mechanical model for an elastic medium in cylindrical coordinates [33]. The problem can be described using the governing partial differential equations of elasticity, the generalized thermo-elastic Hooke’s law including thermal stress contributions, and the heat diffusion equation [34]. Details regarding the problem formulation can be found in the Supplemental Document. We modeled the laser heating as an external heat source flux ${q_e}$ in the absorbing layer (layer 3), using the Beer-Lambert law, where the heat flux is:

*P*, the incident laser power on layer 3,

_{0}*a*is the radius of the top-hat intensity profile, ${\phi _{laser}}$ is the spatio-temporal profile of the laser beam described below and ${\mu _a}$ is the optical absorption coefficient.

_{laser}*P*was calculated by taking into account the transmission through layer 1, which is scattering. The transmission of layer 1 was measured with the same diode laser used for heating and a power meter by comparing the transmitted power through an assembly comprising layers 1 to 5 (the tissue model) and a second assembly comprising layers 2 to 5 only to isolate the contribution from layer 1. Transmission of layer 1 is 50.0 ± 0.5%.

_{0}The beam spatio-temporal profile was modeled as a square pulse in time and a square top-hat profile in space convolved with a Gaussian kernel to simulate blur and defocus [35]:

*H*is the Heaviside function,

*t*is the pulse duration (

_{0}*t*= 1 ms), and ${\sigma _{defocus}}$ is the Gaussian kernel width. For simplicity, we converted the volumetric heat source to a surface heat source of the following form:

_{0}*l*is the thickness of layer 3. Given the diffusivity of PDMS, the volumetric heat source and the surface heat source become equivalent after the characteristic time of diffusion over the thickness of layer 3. This characteristic time is about 2 ms so fitting of the ΔOPL data after 2 ms using a surface heat source is justified and material parameters can be accurately fitted with this simplification. It is important to keep in mind that the calculated temperature at the absorbing layer differs whether a volumetric or surface heat source is used. Using a volumetric or a surface heat source does not, however, affect temperature precision, which depends on the uncertainty of the fitted values.

To solve for the displacements and temperature fields, which are necessary to calculate the changes in optical paths (ΔOPL), we used the state space approach applied to layered media, introduced by Bufler [36] and Bahar [37]. This approach consists in reducing the partial differential equations of coupled thermo-elasticity using integral transforms to a state space equation, taking advantage of the planar layered structure. The relationship between states vectors can then be described by a transfer matrix [37,38] or a stiffness matrix [39] obtained by solving the state space equation. This method has been extensively used in macro-scale multi-layer problems in geophysics [40,41], with applications to seismology and nuclear waste management, but never at the microscale for biomaterials. The thermo-elastic problem can be analytically solved, in the transformed domain, for planar layers in a semi-infinite half space, which makes it particularly appropriate for the retina, as long as we assume a heat source radius being relatively small compared to the retina radius of curvature. In comparison, finite element methods (FEM) treat semi-infinite space by significantly extending the computational domain, therefore increasing the computation time. Moreover, for an axisymmetric configuration, inversing the solution from the transformed domain requires integration over two dimensions (the radial and temporal dimensions), whereas FEM perform integrations over three dimensions (the radial, axial, and temporal dimensions), which increases the computation time. The increased computation efficiency of the state space approach has been discussed in Refs. [42,43] suggesting a >10 × acceleration. Validation of the method and efficiency considerations can be found in the Supplemental Document (Stiffness matrix method (SMM) validation section and corresponding Dataset 2 [44]). Rapid computing of the solution is critical for fitting convergence, fitting quality evaluation via bootstrapping, and precision assessment via Monte-Carlo simulations, all of which involve calculations of many solutions. We further note that modeling a volumetric heat source is possible using the stiffness matrix method [34] but at the expense of longer computation time, which is acceptable after material properties are fitted.

We followed the formulation by Ai et al. [34] to derive the state vector ${[{{{\tilde{u}}_r}^1,{{\tilde{u}}_z}^0,\; {{\tilde{\theta }}^0}} ]^T}$, where ${\tilde{u}_r}^1$ and ${\tilde{u}_z}^0$ are the 1^{st} and 0^{th}-order Hankel-Laplace (HL) transforms of the normal displacement components in the *r* and *z* directions, and ${\tilde{\theta }^0}$ is the 0^{th}-order HL transform of the temperature change. Numerical integration was then performed to find ${u_r}$, ${u_z}$, and $\theta $. The derivation of the state vector and the inverse numerical transform procedures are described in the Supplemental Document (Transfer and stiffness matrix derivation section).

#### 2.5 Arrhenius integral for determination of the therapeutic window

Thermal damage in retinal cells during transient hyperthermia has been approximated following the Arrhenius model, which assumes a first-order kinetics in denaturation of the critical molecular component, the ‘weak link’. Its concentration *D(t)* decreases at the rate [8]:

*A*is the rate constant,

*E**is the activation energy,

*R*is the gas constant (8.314 J.K

^{-1}.mol

^{-1}) and

*T*is the temperature as a function of time. The total decrease of

*D*relative to its initial value

*D*over the course of hyperthermia with duration

_{0}*t*can be expressed by the Arrhenius integral $\mathrm{\Omega }$:

_{h}Values for *A* (1.6×10^{55} s^{-1}) and *E** (340 kJ/mol) have been calibrated and characterized in previous studies so that $\mathrm{\Omega } = 1$ corresponds to the RPE damage threshold [8]. No cellular response was observed for $\mathrm{\Omega } < 0.1$ thus the non-damaging therapeutic window was defined as $0.1 < \mathrm{\Omega } < 1$ [7]. The Arrhenius integral can be directly calculated using the temperature course at any chosen point in the tissue.

## 3. Results

#### 3.1 High-speed imaging of laser heating

The phase difference can be measured between any arbitrary pair of points. However, because we are interested in measuring the absolute vertical displacements of the layers, it is prudent to use a reference point whose phase remains largely unchanged following laser heating. Here, we have chosen a point at the surface of the sample (plane A, *z *= 0) far enough (*r _{ref }*= 505 μm) laterally to be least affected by the thermal expansion while maintaining a high SNR (Fig. 2(a)). Selection of such a reference point is possible since, with the line-field OCT, we have simultaneous access to the phase at all radial positions in a common-path configuration (in the sample arm only) [45]. In contrast, point-scan OCT only allows common-path phase referencing along the z axis (axial referencing) at a single radial position.

The complex values of the OCT signal at the three planes were averaged over 3 pixels in depth around the plane *z*-coordinates (± 1 pixel) before computing the phase values $\varphi ({r,{z_A},t} )$, $\varphi ({r,{z_B},t} )$, $\varphi ({r,{z_C},t} )$. The phase reference $\varphi ({{r_{ref}},{z_A},t} )$ was obtained by averaging the complex values of the OCT signal at the plane A *z*-coordinate and a radius offset of ${r_{ref}}$ = 505 µm (± 1 pixel in *r* and *z*). The phase at ${r_{ref}}$ = -505 µm was also calculated to retrieve and cancel the surface tilt angle in the *r*-*z* plane.

The optical path difference was tracked as a function of time for the three planes of interest: the free surface (plane A, *z _{A}* = 0), an intermediate plane in layer 1 (plane B,

*z*

_{B}_{ }= 25 μm), and at the absorbing layer (plane C,

*z*= 60 μm). For those planes, the ΔOPLs can be expressed as:

_{C}Figure 2(b) shows the $\overline {\textrm{AA}} $, $\overline {\textrm{AB}} $, and $\overline {\textrm{AC}} $ at the center of the heating beam (*r* = 0). The experimental OPL noise, 5.1 nm at *r* = 0, *z* = 0, is close to the theoretical value, 4.0 nm for a phase difference calculated between points with SNR values of 30 dB and 22 dB (lateral offset point) [32].

Assuming for simplicity of interpretation that the reference phase is unchanged, $\overline {\textrm{AA}} $ is directly opposite to the top surface displacement. Upon laser heating, the temperature of layer 3 increases and heat diffuses to adjacent layers. Heating leads to thermal expansion, so the layers deform and become thicker. It causes the top layer to move up (toward negative *z*) and $\overline {\textrm{AA}} $ to decrease by the vertical displacement of the plane A (${u_z}({r,\; {z_A} = 0} )$). Meanwhile, as the temperature increases by *θ*, the index of refraction of PDMS decreases by $\Delta n$, with $\Delta n = {\alpha _{TO}}\theta $, where ${\alpha _{TO}}$, the thermo-optic coefficient of PDMS, is negative. Therefore, as plane B moves up, the optical path $\overline {\textrm{AB}} $ decreases by the additive effects of a shorter physical distance and a lower index of refraction. Since plane A also moves up and air is replaced by the polymer over the displacement thickness, the optical path $\overline {\textrm{AB}} $ decrease is slightly reduced by this effect. Same goes for $\overline {\textrm{AC}} $. Because the thermal expansion (${\alpha _{TE}})\; $ and thermo-optic (${\alpha _{TO}}$) coefficients are intimately related, their combined effects resulted in $\overline {\textrm{AA}} $, $\overline {\textrm{AB}} $, and $\overline {\textrm{AC}} $ being of similar amplitude (Fig. 2(b)). The simplified expressions of $\overline {\textrm{AA}} $, $\overline {\textrm{AB}} $, and $\overline {\textrm{AC}} $ are:

*u*is normal displacement components in the

_{z}*z*direction, ${n_{RT}}$ is the index of refraction of PDMS at room temperature, and $\overline {\Delta n} $ is the average change in refractive index induced by heating between the planes of interest. The exact expression and the derivation of these ΔOPLs can be found in the Supplemental Document (Eqs. S2 and S7). For comparison, the change in optical path length that would be measured between the planes A and C using axial referencing, i.e. the phase difference between the planes at the same radial coordinate, would be:

*u*and $\theta $, used to calculate $\Delta n$, which allow exact calculations of the ΔOPLs.

_{z}#### 3.2 Retrieval of the material properties

We assumed that the properties of the glass substrate are well known and, given the separation between the top of the assembly and the glass, that the contribution from the glass to the response was very small. We considered for simplicity that all PDMS layers have the same isotropic thermal and mechanical properties, including the layers containing titanium oxide particles, given their low weight fraction (2.5 wt%). The density *ρ*, the Poisson’s ratio *ν*, and the specific heat capacity *c _{p}* are all well characterized for PDMS [28,29] and are considered fixed parameters. Similarly, we assumed the index of refraction

*n*at room temperature (23°C) and the thermo-optic coefficient

*α*[30] to be known. On the other hand, the coefficient of thermal expansion

_{TO}*α*strongly depends on the curing temperature, which leads to uncertainty in its value [28], and a range of coefficient of thermal conductivity

_{TE}*κ*can be found in the literature. We therefore considered

*α*and

_{TE}*κ*as fitting parameters. Although the Young’s modulus

*E*depends on the curing temperature [27], it has little effect on the state vector (see Supplemental Document

*,*Elastic modulus effects section). Lastly, the absorption coefficient ${\mu _a}$ was fitted through the deposited heat magnitude ${q_0}$ with ${q_0} \propto 1 - {e^{ - {\mu _a}l}}$, where

*l*is the absorbing layer thickness. We further assumed that these properties are temperature independent, given the small temperature range of the study (<10°C). The Gaussian blur to the ideal disk (top-hat) intensity profile (Eq. (5)) [35] added two geometrical fitting parameters to the model (width of the blur

*σ*and top-hat profile radius

_{defocus}*a*). This brought the total number of fitting parameters to five:

_{laser}*α*,

_{TE}*q*,

_{0}*κ, σ*.

_{defocus}, and a_{laser}Fitting of the parameters was performed in two steps. First, these properties were fitted by minimizing the sum of the squares of the residuals between the model and the experimental data $\overline {\textrm{AA}} $, which showed the highest SNR, using the exact ΔOPL expression (equation S2). The laser geometrical parameters, *σ _{defocus} and a_{laser}*, were mainly fitted through the spatial shape of the response, and the coefficient of thermal conductivity was mainly fitted through the temporal decay of the response. Importantly, however, $\overline {\textrm{AA}} $ scales linearly with both

*α*and

_{TE}*q*so fitting of $\overline {\textrm{AA}} $ did not allow decorrelation between

_{0}*α*and

_{TE}*q*, therefore the product ${\alpha _{TE}} \times {q_0}$ was fitted through the magnitude of the response. Second, to find

_{0}*α*and

_{TE}*q*, we used the experimental data $\overline {\textrm{AC}} $, which is only proportional to

_{0}*q*. The fitted parameters found in the first step were used to solve for $\overline {\textrm{AC}} $, and then amplitude of $\overline {\textrm{AC}} $ was fitted to determine the absolute values of

_{0}*α*and

_{TE}*q*, as described in the Supplemental Document (Material property fitting procedure section)

_{0}*.*The comparison between the fitted model and the data is shown in Fig. 3 for $\overline {\textrm{AA}} $ and $\overline {\textrm{AC}} $. The data are well reproduced in both the spatial (radial) and temporal dimensions, thus suggesting high fidelity of the model. The uniformity of the residuals (see Supplemental Document) over the entire data set also attests for the accuracy of the model and the absence of systematic errors. The resulting fitted values are listed in Table 3. We note a good agreement between the fitted and initial values for

*α*(286 vs 310 K

_{TE}^{-1}),

*κ*(0.175 vs 0.17 W.m

^{-1}.K

^{-1}), and

*μ*(20.3 mm

_{a}^{-1}vs 19.9 mm

^{-1}). The uncertainties in parameter determination were obtained by fitting bootstrapped data, resampled over a total of 200 iterations, and are listed in Table 3. The obtained uncertainties are an order of magnitude lower than for the starting values of

*α*and

_{TE}*κ*, hence demonstrating the high precision of the parameter determination using this method. We also note that the obtained precision on

*μ*(2.5%) is much higher than the required accuracy for therapeutic titration (20%). We finally verified the low correlation between the parameters (see Parameter correlations section in Supplemental Document) obtained by repeated fitting of the bootstrapped data.

_{a}#### 3.4 Temperature estimation and precision

Since we are ultimately interested in the temperature, the best-fit parameters can be used to compute the temperature increase *θ* at any point in space and time. For non-damaging retinal therapy, the temperature rise at the RPE, represented here by the absorbing PDMS layer, is of utmost importance. Figure 4(a) shows the temporal course of the temperature at that layer. Temperature increased during the 1 ms laser pulse, and then decreased slowly due to heat diffusion. Figure 4(b) shows two snapshots of the temperature distribution at the end of the heating pulse (*t* = 1 ms), when heat is relatively confined, and after 4 ms, when heat has diffused radially and axially. To estimate the temperature precision of the calculations, we ran Monte-Carlo simulations [46] using the fitted parameters with their uncertainty over 1000 iterations to compute the peak temperature (*t *= 1 ms) at the center of the beam at the absorbing layer (*r* = 0, *z _{C}* = 60 μm), and found a temperature uncertainty of ±0.25°C. These simulations at a single point in time and space took about 20 seconds. For comparison, running the same number of iterations using the COMSOL model would take approximately 6 hours.

## 4. Discussion

Determination of the temporal and spatial distribution of the temperature rise in tissue is critical for evaluating the safe therapeutic window and the extent of potential thermal damage. Both are typically assessed using the Arrhenius equation, which describes the first order kinetics of the thermal denaturation. The therapeutic window for non-damaging hyperthermia corresponds to the Arrhenius integral value between 0.1 and 1, where values above 1 correspond to permanent damage. With the 2.5% precision of the absorption coefficient determination in our approach, the uncertainty in $\mathrm{\Omega }$ is about 0.17, which is well within the therapeutic window. Uncertainty of 0.50 in Arrhenius integral, approaching the boundary of the therapeutic window corresponds to 2.5 times lower precision in determination of the absorption coefficient (i.e., 6.25%). At the current laser settings, this would translate to a temperature rise uncertainty of about 0.65°C. Precision of this method exceeds the other methodologies for temperature-controlled retinal treatment and is more than sufficient for titration the non-damaging thermal therapy. In the current example, the absorption coefficient determination and the temperature calculations are performed way below the damage threshold (Ω < 10^{−3}), which enables repeated measurements (multiple heating pulses) for further data averaging.

We are working on adapting this approach to *in-vivo* measurements and hope to achieve a similar level of the temperature precision despite the potential difficulties associated with tissue movements and reduced OCT signal strength, which directly affects the phase noise. Here, with the tissue phantom, we only considered the phase information from two planes (A and C), as they enable sufficient precision in determination of the material parameters. However, more data points from additional planes (e.g., plane B) are available for fitting, which would increase the precision and allow additional fitting parameters. On one hand, retina offers many scattering planes from which phase can be extracted, while on the other hand, we expect a larger number of variables for the retina. It remains to be seen whether the increased number of experimental data points provided by the additional layers will compensate for the increased number of unknown parameters.

Retinal thermal expansion coefficient is often assumed to be equal to that of water [19], but a recent study suggested otherwise and estimated it to be four times larger [18]. We can also expect this coefficient to be anisotropic (albeit transverse isotropic) for the retina, given its anisotropic cellular structure. The same might apply to the coefficient of thermal conductivity and, for the same reason, layers could exhibit different Poisson’s ratios. Conversely, the index of refraction of the retina is relatively uniform within 1% across the layers and has been reasonably determined [47]. Questions remain, however, regarding the thermo-optic coefficient (${\alpha _{TO}}$), which might more strongly depend on the solid content [48,49]. It also has been shown that the elastic modulus varies throughout the retina [21], and is two-to-three orders of magnitude lower than for PDMS (kPa instead of MPa). We therefore expect it to play a more significant role than in the tissue phantom. In addition, the present model assumes a uniform laser irradiance and a uniform absorption coefficient within the laser spot. Additional contributions from metabolic heat generation and convection due to blood perfusion can be implemented in the model, although the metabolic heat is generally considered negligible [50]. The characteristic time constant of the perfusion heat transfer is on the order of seconds, and therefore it is more relevant for seconds-long laser exposures [51] rather than ms-short pulsed therapies. For *in-vivo* measurements, we can ignore potential cellular-scale variations (${\sim} $ 30 μm) in the absorption coefficient and laser-intensity variations (${\sim} $ 10 μm) in the heat source profile. since the variations would be homogenized within 1 ms of the laser pulse, as heat diffuses by about 30 μm over that duration.

## 5. Conclusion

We present a methodology for determination of the optical and thermal parameters in multi-layered materials, applicable to retinal photothermal therapies. We demonstrated that, by matching the optical path length changes, recorded with line-field phase-resolved optical coherence tomography, across the beam width and along the full course of heating and cooling, one can precisely identify the material properties. Based on model fittings, we estimated the temperature with a precision of about 0.2°C for a heating amplitude not exceeding 8°C during 1 ms—well within the safety range of the retina [5]. *In-situ* measurements of the retinal heating in every patient should allow precise titration of the energy deposition for photocoagulation and non-damaging laser therapy. We finally note that the present methodology is also generally applicable to any planar multi-layered assembly, such as skin, for example, where the method could complement the current temperature monitoring techniques, typically based on thermal cameras [52,53]. Beyond biomaterials and tissues, the method could be extended, for example, to coating and surface metrology [54] or artwork studies [55].

## Funding

National Institutes of Health (P30 EY 001730, R01 EY 029710, U01 EY 025501, U01 EY 032055); Air Force Office of Scientific Research (FA9550-20-1-0186).

## Acknowledgements

We thank Fabio Feroldi, Austin Roorda, and Junya Hattori for discussions regarding the optical design. DV thanks Ludwig Galambos for profilometry measurements, Zhi Yong Ai for insights on the transfer and stiffness matrix method, and Mohajeet Bhuckory for discussions on the translation to *in-vivo* measurements.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

All study data are included in the article and/or in Supplement 1. The thermo-mechanical COMSOL models for a rodent retina and the tissue phantom are available online as Dataset 1 Refs. [11], Dataset 2, respectively, Refs. [44].

## Supplemental document

See Supplement 1 for supporting content.

## References

**1. **C. M. Moorman and A. M. P. Hamilton, “Clinical applications of the MicroPulse diode laser,” Eye **13**(2), 145–150 (1999). [CrossRef]

**2. **C. M. Lee and R. J. Olk, “Modified grid laser photocoagulation for diffuse diabetic macular edema,” Ophthalmology **98**(10), 1594–1602 (1991). [CrossRef]

**3. **E. Burumcek, A. Mudun, S. Karacorlu, and M. O. Arslan, “Laser photocoagulation for persistent central serous retinopathy: Results of long-term follow-up,” Ophthalmology **104**(4), 616–622 (1997). [CrossRef]

**4. **J. Chhablani, Y. J. Roh, A. I. Jobling, E. L. Fletcher, J. J. Lek, P. Bansal, R. Guymer, and J. K. Luttrull, “Restorative retinal laser therapy: Present state and future directions,” Surv. Ophthalmol. **63**(3), 307–328 (2018). [CrossRef]

**5. **C. Sramek, M. Mackanos, R. Spitler, L. S. Leung, H. Nomoto, C. H. Contag, and D. Palanker, “Non-damaging retinal phototherapy: dynamic range of heat shock protein expression,” Investig. Ophthalmol. Vis. Sci. **52**(3), 1780–1787 (2011). [CrossRef]

**6. **D. Lavinsky and D. Palanker, “Nondamaging photothermal therapy for the retina: Initial clinical experience with chronic central serous retinopathy,” Retina **35**(2), 213–222 (2015). [CrossRef]

**7. **D. Lavinsky, J. Wang, P. Huie, R. Dalal, S. J. Lee, D. Y. Lee, and D. Palanker, “Nondamaging retinal laser therapy: Rationale and applications to the macula,” Investig. Ophthalmol. Vis. Sci. **57**(6), 2488–2500 (2016). [CrossRef]

**8. **C. Sramek, Y. Paulus, H. Nomoto, P. Huie, J. Brown, and D. Palanker, “Dynamics of retinal photocoagulation and rupture,” J. Biomed. Opt. **14**(3), 034007 (2009). [CrossRef]

**9. **S. Y. Schmidt and R. D. Peisch, “Melanin concentration in normal human retinal pigment epithelium: regional variation and age-related reduction,” Investig. Ophthalmol. Vis. Sci. **27**(7), 1063–1067 (1986).

**10. **G. Goetz, T. Ling, T. Gupta, S. Kang, J. Wang, P. D. Gregory, B. Hyle Park, and D. Palanker, “Interferometric mapping of material properties using thermal perturbation,” Proc. Natl. Acad. Sci. U. S. A. **115**(11), E2499–E2508 (2018). [CrossRef]

**11. **D. Veysset, T. Ling, Y. Zhuo, V. P. Pandiya, R. Sabesan, and D. Palanker, Interferometric imaging of thermal expansion for temperature control in retinal laser therapy: dataset 1, figshare (2021), https://doi.org/10.6084/m9.figshare.17036288.

**12. **R. Brinkmann, Stefan Koinzer, K. Schlott, L. Ptaszynski, M. Bever, A. Baade, S. Luft, Y. Miura, J. Roider, and R. Birngruber, “Real-time temperature determination during retinal photocoagulation on patients,” J. Biomed. Opt. **17**(6), 061219 (2012). [CrossRef]

**13. **A. Baade, C. von der Burchard, M. Lawin, S. Koinzer, B. Schmarbeck, K. Schlott, Y. Miura, J. Roider, R. Birngruber, and R. Brinkmann, “Power-controlled temperature guided retinal laser therapy,” J. Biomed. Opt. **22**(11), 1 (2017). [CrossRef]

**14. **G. Schüle, G. Hüttmann, C. Framme, J. Roider, and R. Brinkmann, “Noninvasive optoacoustic temperature determination at the fundus of the eye during laser irradiation,” J. Biomed. Opt. **9**(1), 173 (2004). [CrossRef]

**15. **K. Schlott, S. Koinzer, L. Ptaszynski, M. Bever, A. Baade, J. Roider, R. Birngruber, and R. Brinkmann, “Automatic temperature controlled retinal photocoagulation,” J. Biomed. Opt. **17**(6), 061223 (2012). [CrossRef]

**16. **V. E. Gusev and A. A. Karabutov, * Laser Optoacoustics* (American Institute of Physics, 1993).

**17. **E. Seifert, J. Tode, A. Pielen, D. Theisen-Kunde, C. Framme, J. Roider, Y. Miura, R. Birngruber, and R. Brinkmann, “Selective retina therapy: toward an optically controlled automatic dosing,” J. Biomed. Opt. **23**(11), 1 (2018). [CrossRef]

**18. **H. H. Müller, L. Ptaszynski, K. Schlott, C. Debbeler, M. Bever, S. Koinzer, R. Birngruber, R. Brinkmann, and G. Hüttmann, “Imaging thermal expansion and retinal tissue changes during photocoagulation by high speed OCT,” Biomed. Opt. Express **3**(5), 1025 (2012). [CrossRef]

**19. **S. Makita and Y. Yasuno, “In vivo photothermal optical coherence tomography for non-invasive imaging of endogenous absorption agents,” Biomed. Opt. Express **6**(5), 1707 (2015). [CrossRef]

**20. **V. Y. Zaitsev, A. L. Matveyev, L. A. Matveev, G. V. Gelikonov, A. I. Omelchenko, O. I. Baum, S. E. Avetisov, A. V. Bolshunov, V. I. Siplivy, D. V. Shabanov, A. Vitkin, and E. N. Sobol, “Optical coherence elastography for strain dynamics measurements in laser correction of cornea shape,” J. Biophotonics **10**(11), 1450–1463 (2017). [CrossRef]

**21. **Y. Qu, Y. He, Y. Zhang, T. Ma, J. Zhu, Y. Miao, C. Dai, M. Humayun, Q. Zhou, and Z. Chen, “Quantified elasticity mapping of retinal layers using synchronized acoustic radiation force optical coherence elastography,” Biomed. Opt. Express **9**(9), 4054 (2018). [CrossRef]

**22. **Y. Zhao, Z. Chen, C. Saxer, S. Xiang, J. F. de Boer, and J. S. Nelson, “Phase-resolved optical coherence tomography and optical Doppler tomography for imaging blood flow in human skin with fast scanning speed and high velocity sensitivity,” Opt. Lett. **25**(2), 114 (2000). [CrossRef]

**23. **V. P. Pandiyan, A. Maloney-Bertelli, J. A. Kuchenbecker, K. C. Boyle, T. Ling, Z. C. Chen, B. H. Park, A. Roorda, D. Palanker, and R. Sabesan, “The optoretinogram reveals the primary steps of phototransduction in the living human eye,” Sci. Adv. **6**(37), eabc1124 (2020). [CrossRef]

**24. **K. C. C. Boyle, Z. C. C. Chen, T. Ling, V. P. P. Pandiyan, J. Kuchenbecker, R. Sabesan, and D. Palanker, “Mechanisms of light-induced deformations in photoreceptors,” Biophys. J. **119**(8), 1481–1488 (2020). [CrossRef]

**25. **V. P. Pandiyan, X. Jiang, A. M. Bertelli, J. A. Kuchenbecker, R. Sabesan, A. Maloney-Bertelli, J. A. Kuchenbecker, U. Sharma, and R. Sabesan, “High-speed adaptive optics line-scan OCT for cellular-resolution optoretinography,” Biomed. Opt. Express **11**(9), 5274 (2020). [CrossRef]

**26. **J. Baxi, W. Calhoun, Y. J. Sepah, D. X. Hammer, I. Ilev, T. Joshua Pfefer, Q. D. Nguyen, and A. Agrawal, “Retina-simulating phantom for optical coherence tomography,” J. Biomed. Opt. **19**(2), 021106 (2014). [CrossRef]

**27. **I. D. Johnston, D. K. McCluskey, C. K. L. Tan, and M. C. Tracey, “Mechanical characterization of bulk Sylgard 184 for microfluidics and microengineering,” J. Micromech. Microeng. **24**(3), 035017 (2014). [CrossRef]

**28. **A. Müller, M. C. Wapler, and U. Wallrabe, “A quick and accurate method to determine the Poisson’s ratio and the coefficient of thermal expansion of PDMS,” Soft Matter **15**(4), 779–784 (2019). [CrossRef]

**29. **G. Zhang, Y. Sun, B. Qian, H. Gao, and D. Zuo, “Experimental study on mechanical performance of polydimethylsiloxane (PDMS) at various temperatures,” Polym. Test. **90**(2020), 106670 (2020). [CrossRef]

**30. **H. Gao, H. Hu, Y. Zhao, J. Li, M. Lei, and Y. Zhang, “Highly-sensitive optical fiber temperature sensors based on PDMS/silica hybrid fiber structures,” Sensors Actuators A Phys. **284**, 22–27 (2018). [CrossRef]

**31. **A. F. Fercher, C. K. Hitzenberger, M. Sticker, R. Zawadzki, B. Karamata, and T. Lasser, “Dispersion compensation for optical coherence tomography depth-scan signals by a numerical technique,” Opt. Commun. **204**(1-6), 67–74 (2002). [CrossRef]

**32. **B. Hyle Park, M. C. Pierce, B. Cense, S. Yun, M. Mujat, G. J. Tearney, B. E. Bouma, and J. F. de Boer, “Real-time fiber-based multi-functional spectral-domain optical coherence tomography at 13 µm,” Opt. Express **13**(11), 3931 (2005). [CrossRef]

**33. **M. A. A. Biot, “Thermoelasticity and irreversible thermodynamics,” J. Appl. Phys. **27**(3), 240–253 (1956). [CrossRef]

**34. **Z. Y. Ai, Z. Zhao, and L. J. Wang, “Thermo-mechanical coupling response of a layered isotropic medium around a cylindrical heat source,” Comput. Geotech. **83**, 159–167 (2017). [CrossRef]

**35. **S. Zhuo and T. Sim, “Defocus map estimation from a single image,” Pattern Recognit. **44**(9), 1852–1858 (2011). [CrossRef]

**36. **H. Bufler, “Theory of elasticity of a multilayered medium,” J. Elast. **1**(2), 125–143 (1971). [CrossRef]

**37. **L. Y. Bahar, “Transfer Matrix Approach to Layered Systems,” J. Eng. Mech. Div. **98**(5), 1159–1172 (1972). [CrossRef]

**38. **Z. Q. I. Yue and J. H. Yin, “Backward transfer-matrix method for elastic analysis of layered solids with imperfect bonding,” J. Elast. **50**(2), 109–128 (1998). [CrossRef]

**39. **L. Wang and S. I. Rokhlin, “Stable reformulation of transfer matrix method for wave propagation in layered anisotropic media,” Ultrasonics **39**(6), 413–424 (2001). [CrossRef]

**40. **J. Wang and S. Fang, “The state vector solution of axisymmetric Biot’s consolidation problems for multilayered poroelastic media,” Mech. Res. Commun. **28**(6), 671–677 (2001). [CrossRef]

**41. **Z. Y. Ai, Z. Y. Cheng, and J. Han, “State space solution to three-dimensional consolidation of multi-layered soils,” Int. J. Eng. Sci. **46**(5), 486–498 (2008). [CrossRef]

**42. **Z. Y. Ai, Z. K. Ye, and J. J. Yang, “Thermo-mechanical behaviour of multi-layered media based on the Lord-Shulman model,” Comput. Geotech. **129**, 103897 (2021). [CrossRef]

**43. **Z. Y. Ai and Y. C. Cheng, “Extended precise integration method for consolidation of transversely isotropic poroelastic layered media,” Comput. Math. with Appl. **68**(12), 1806–1818 (2014). [CrossRef]

**44. **D. Veysset, T. Ling, Y. Zhuo, V. P. Pandiya, R. Sabesan, and D. Palanker, Interferometric imaging of thermal expansion for temperature control in retinal laser therapy: dataset 2, figshare (2021), https://doi.org/10.6084/m9.figshare.17036279.

**45. **Z. Yaqoob, W. Choi, S. Oh, N. Lue, Y. Park, C. Fang-Yen, R. R. Dasari, K. Badizadegan, and M. S. Feld, “Improved phase sensitivity in spectral domain phase microscopy using line-field illumination and self phase-referencing,” Opt. Express **17**(13), 10681 (2009). [CrossRef]

**46. **P. Doubilet, C. B. Begg, M. C. Weinstein, P. Braun, and B. J. Mcneil, “Probabilistic sensitivity analysis using monte carlo simulation: a practical approach,” Med. Decis. Mak. **5**(2), 157–177 (1985). [CrossRef]

**47. **E. Chen, “Refractive indices of the rat retinal layers,” Ophthalmic Res. **25**(1), 65–68 (1993). [CrossRef]

**48. **I. Y. Yanina, E. N. Lazareva, and V. V. Tuchin, “Refractive index of adipose tissue and lipid droplet measured in wide spectral and temperature ranges,” Appl. Opt. **57**(17), 4839 (2018). [CrossRef]

**49. **Y. L. Jin, J. Y. Chen, L. Xu, and P. N. Wang, “Refractive index measurement for biomaterial samples by total internal reflection,” Phys. Med. Biol. **51**(20), N371–N379 (2006). [CrossRef]

**50. **L. T. D. Truong, P. J. Lesniewski, and B. Wedding, “Heat transfer simulation in laser irradiated retinal tissues,” Biomed. Phys. Eng. Express **8**(1), 015027 (2022). [CrossRef]

**51. **J. Kandulla, H. Elsner, R. Birngruber, and R. Brinkmann, “Noninvasive optoacoustic online retinal temperature determination during continuous-wave laser irradiation,” J. Biomed. Opt. **11**(4), 041111 (2006). [CrossRef]

**52. **J. Kosir, D. Vella, and M. Jezersek, “Non-contact monitoring of the depth temperature profile for medical laser scanning technologies,” Sci. Rep. **10**(1), 20242 (2020). [CrossRef]

**53. **J. Kosir, D. Vella, M. Lukac, and M. Jezersek, “Towards personalized and versatile monitoring of temperature fields within heterogeneous tissues during laser therapies,” Biomed. Opt. Express **12**(7), 4530 (2021). [CrossRef]

**54. **D. Stifter, “Beyond biomedicine: a review of alternative applications and developments for optical coherence tomography,” J. Micromech. Microeng. **88**(3), 337–357 (2007). [CrossRef]

**55. **P. Targowski, M. Góra, and M. Wojtkowski, “Optical coherence tomography for artwork diagnostics,” Laser Chem. **2006**, 1–11 (2006). [CrossRef]