Modeling nonlinear microscopy near index-mismatched interfaces

Nonlinear microscopy is widely used to characterize thick, optically heterogeneous biological samples. While quantitative image analysis requires accurately describing the contrast mechanisms at play, the majority of established numerical models neglect the influence of field distortion caused by sample heterogeneity near focus. In this work, we show experimentally and numerically that finite-difference time-domain (FDTD) methods are applicable to model focused fields interactions in the presence of heterogeneities, typical of nonlinear microscopy. We analyze the ubiquitous geometry of a vertical interface between index-mismatched media (water, glass, and lipids) and consider the cases of two-photon-excited fluorescence (2PEF), third-harmonic generation (THG) and polarized THG contrasts. We show that FDTD simulations can accurately reproduce experimental images obtained on model samples and in live adult zebrafish, in contrast with previous models neglecting field distortions caused by index mismatch at the micrometer scale. Accounting for these effects appears to be particularly critical when interpreting coherent and polarization-resolved microscopy data.


INTRODUCTION
Multiphoton microscopy [1] is widely used for imaging thick biological samples. Beyond its unique performance for in-depth imaging with µm 3 resolution, one of its distinctive advantages is the variety of nonlinear processes that can be used as contrast mechanisms to provide complementary information on biological samples: two-and three-photon-excited fluorescence (2PEF, 3PEF) [2][3][4][5], second-and third-harmonic generation (SHG, THG) [6][7][8][9][10], or coherent Raman anti-Stokes scattering (CARS) [11]. While for fluorescence imaging based on 2PEF (respectively, 3PEF) the signal is simply described as the squared (resp. cubed) intensity of the excitation beam convoluted with the concentration of fluorophores, the contrast mechanisms are less intuitive in the case of coherent nonlinear imaging relying on harmonic generation or coherent Raman processes. In coherent nonlinear microscopy, the signal and scattering directions result from the interplay between the excitation field distribution and the sample microstructure, and quantitative image interpretation therefore requires modeling. A simulation strategy based on the angular spectrum representation (ASR) to calculate the excitation field distribution near focus and on Green's functions to propagate the nonlinear response from the focal region to the detector plane has gained popularity and has been extensively used for that purpose [12][13][14][15]. However, this approach and most other established models neglect the influence of sample refractive index heterogeneity near the focus, concentrating instead exclusively on the distribution of nonlinear indices. One exception is the specific geometry of interfaces normal to the direction of propagation where semianalytic solutions have been derived [16,17]. While near-focus linear propagation effects have usually been neglected for the sake of simplicity, it can be intuited that they have important consequences on multiphoton signals, and even more on their polarization-resolved versions [15,[18][19][20][21]. This is particularly expected near vertical interfaces between two dielectric media (e.g. water/lipids), which happen to be present in most microscopy images of biological samples.
Recent developments in numerical methods have opened interesting perspectives for investigating phenomena such as beam propagation at large depths (>100 µm) inside tissues. In particular, approaches such as graphics processing unit (GPU)-accelerated Monte Carlo [22,23] or hybrid approaches [24,25] have been successful in accurately modeling beam propagation for light-sheet fluorescence microscopy, albeit with assumptions that are not appropriate for the tightly focused beams used in point-scanning nonlinear microscopy due to the loss of phase information.
Along a complementary direction, the finite-diffference time-domain (FDTD) method has been used to simulate several microscopy techniques, including wide-field, confocal, and phase-contrast microscopy [26][27][28], and more recently to calculate aberrations induced when light propagates through bone [29] or 2334-2536/21/070944-08 Journal © 2021 Optical Society of America brain tissue [30]. Finally, there have been recent demonstrations of using FDTD to evaluate theoretically the consequences of index mismatch on CARS and SHG microscopy signal level [31,32], and near-field effects close to a gold nanoparticle [33]. These studies, however, did not include experimental validations and relied on custom FDTD implementations that hindered their reproduction and extension by other investigators.
In this work, we revisit experimentally and numerically the commonly encountered geometry of a vertical interface between index-mismatched media imaged using incoherent (2PEF), coherent (THG) and polarization-resolved coherent (P-THG) multiphoton microscopy. We show that the ASR/Green model fails to reproduce experimental observations because it neglects field distortion near-focus. In contrast, we show that an FDTDbased approach accurately accounts for experimentally observed artifacts, with important consequences for the interpretation of coherent and polarization-resolved images.

A. FDTD Calculations
We outline in the following section the models used for FDTD and ASR/Green calculations. All scripts and codes used to generate the figures are available on Zenodo [34].
The FDTD family of methods is now well documented in several books [35,36] and reviews [37], which describe the different implementations with their advantages and their limits. We chose to work with a standard implementation designed for photonic multiphysics calculations (Lumerical Device suite, Ansys Inc, Canonsburg, PA), to facilitate further comparisons. In short, the electric and magnetic fields are calculated on every point of a 3D grid at successive times by solving discretized Maxwell equations [Eq. (1)] for nonmagnetic materials ( . The simulations are performed in the time domain, but spectral information can be retrieved using Fourier transforms in a postprocessing step, (1) We performed three different FDTD calculations in order to estimate (i) the field distribution when focused on a vertical glass/water interface; (ii) the resulting 2PEF signal assuming the water to be fluorescent; and (iii) the resulting THG. We considered an incoming focused (NA = 1.0) Gaussian beam with a central wavelength of 1.2 µm and 50 nm bandwidth, and performed calculations over a focal region spanning 18 × 18 × 13 µm 3 discretized over 20-40 nm steps with a time resolution of ≈ 0.2 fs. Since storing the entire simulated fields would require several tens of gigabytes, we instead saved the relevant data using different types of monitors to analyze the results: • For linear (focal fields) simulations, we solved the equations for two simple nondispersive materials (n 1 , n 2 ) and we saved the time-integrated electric field components (E x , E y , E z ) for a 3D volume encompassing the expected focal point.
• For 2PEF simulations, we performed similar calculations with the same type of materials, but we used a time monitor to compute the integral of the squared intensity as a function of time.
• For THG simulations, we considered two nonlinear materials, ( . In the Lumerical software formalism, a nonlinear polarization term is introduced explicitly as follows: One limitation of the current version of the software we used is that nonlinear materials can exhibit only diagonal nonzero tensor elements. In order to simulate isotropic materials, we therefore set 0 and only considered excitations polarized along x or y. We then considered a 2D array of detectors in a slab of linear material located at the end of the simulation domain, and computed the integrated intensity around 1/3 of the incoming wavelength (385-415 nm). See Fig. S1 in Supplement 1 for the schematics of the FDTD simulation domain for THG at interfaces. FDTD simulations were run on a Dell Precision 7920 (2× Intel 6140 CPUs, 384 GB RAM, DDR4 2666 MHz) and on a Dell Precision 5820 (Intel Xeon W2175 CPU, 128 GB RAM, DDR4 2666 MHz) running Lumerical 2019b or 2020a with typical computing time of ≈ 1−2 h per condition. Visualization 1 was made using Fiji [38] and ClearVolume [39].

B. ASR/Green Model Calculations
For reference and comparison, we performed ASR/Green-based calculations of THG from the same geometries using the method described in previous studies [12,14,40]. We used the same formalism to compute the excitation intensity for 2PEF simulations and assumed isotropic emission in that case.
The MATLAB code used to generate the figures is available on Zenodo [34].

C. Microscopy
Experiments were performed on an upright lab-built multiphoton microscope equipped with a femtosecond laser source (80 MHz, 1100 nm, 100 fs, Insight X3, Spectra Physics, U.S.), galvanometric scanners (GSI Lumonics, U.S.) and a water immersion objective (25×, 1.05NA, Olympus, Japan). THG detection was performed in transmission using a high NA condenser (Olympus) and a UV filter (Semrock FF01-377/50). 2PEF was detected in the epi direction with appropriate filters (Semrock FF01-590/20). Detection was performed using photon-counting detectors (SensTech, UK), lab-designed counting electronics, and LabVIEW-written software (National Instruments, U.S.). We used an excitation power of 20-100 mW at the sample surface, depending on the experiments. The signal level was kept in a range avoiding saturation of the detection chain, i.e., less than one detected photon every four laser pulses. The incident polarization was linear, and its direction in the focal plane was controlled using a rotating waveplate located just before the objective, as described in [18].
The experimental data used to generate the figures are available on Zenodo [34].

D. Model Samples
In order to obtain a well-defined vertical interface, we cut one edge of a quartz slide (thickness 200 µm, SPI Supplies, U.S.) using a diamond knife. We deposited the coverslip on a glass slide, adjacent to a drop of rhodamine B solution (2.5 µg/mL, Sigma-Aldrich). We covered the sample with a 170 µm thick borosilicate glass coverslip (Menzel, Thermo Scientific Menzel) to let the fluorescent solution spread and surround the quartz slide. 230 µm thick spacers (Bytac surface protector, Saint-Gobain, France) were inserted between the bottom glass slide and the borosilicate coverslip to ensure its horizontality and minimize aberrations.
For lipid/water experiments, oil droplets (olive oil, Puget, France) were embedded in an agarose gel to prevent motion. The agarose solution was prepared by dissolving 1.0 g of low melting point agarose powder (Invitrogen 9012-36-6 ) into 50 mL of distilled water under heating, mixed with the oil to create an emulsion, and a small volume was deposited between a glass slide and a coverslip for imaging.

E. Zebrafish Imaging
Live adult zebrafish imaging was performed using the protocol described in [41]. 12-month-old casper (depigmented) zebrafish [42] were used in the study. Anesthesia was initiated by soaking the fish for 90 s in water containing 0.02% MS222 (Sigma-Aldrich). Fish were then transferred into a water solution of 0.005% (v/v) MS222 and 0.005% (v/v) isoflurane to maintain the anesthesia during imaging, mounted in a plastic dish and held between pieces of sponge. We first recorded in a mosaic manner a large-scale THG map from the skin and skull above the telencephalon and in the dorsal region to locate lipid droplets, and then recorded P-THG images of such droplets. Signals were epidetected using the excitation objective. Overall, one imaging session lasted 30-45 min. All animal experiments were carried out in accordance to the official regulatory standards of the department of Paris (agreement number A914772 to N.D.) and conformed to French and European ethical and animal welfare directives (project authorization from the Ministère de l'Enseignement Supérieur, de la Recherche et de l'Innovation to N.D.). Before imaging the zebrafish, a starch granule was imaged as a polarization reference (see Fig. S2 in Supplement 1).

F. Image Analysis
THG images recorded with different polarization directions were converted into a x y − P stack, with an optional denoising/dedrifting step. Then each pixel was fitted with a cosine-square function to extract the average intensity, the amplitude, and the phase (polarization angle yielding maximum signal) of the modulation. The goodness-of-fit parameter was also saved, and together with the average intensity and amplitude, was used to generate a binary mask, which was then applied to the phase and amplitude maps for Visualization 1. (See Fig. S2 in Supplement 1 for an overview of the image analysis workflow). The FIJI macro used for the image analysis is also on Zenodo [34]. Fig. 1. Principles of the ASR/Green and FDTD strategies for nonlinear microscopy simulations, and excitation intensity distributions calculated around a near-focus vertical interface. (a) Schematic representation of the ASR/Green three-step strategy (see [12,14]); (b) schematic representation of an FDTD simulation: the incident field is propagated until the beginning of the FDTD volume (assuming no distortion), and this initial distribution as well as the induced response are then propagated numerically through a discretized volume for successive time steps; (c) near-focus intensity distribution calculated with the ASR formalism assuming an unaberrated incident wavefront and neglecting index heterogeneity near the focus; (d) near-focus intensity distribution calculated using FDTD when propagating a similar beam 7 µm along a vertical interface located near the focus. The FDTD calculation predicts strong local distortions. Scale bar = 1 Scale bar = 1 µm. See also Visualization 1.

RESULTS
We first used FDTD calculations to evaluate the field distortions experienced by a beam near a vertical interface. As depicted in Fig. 1, we start by considering a tightly focused Gaussian beam (NA = 1) propagating along a water-quartz interface. We chose these materials because their linear (n = 1.33 for water, n = 1.45 for quartz) and nonlinear properties are well described [43], and similar index mismatches are found in biological tissues. We calculated the intensity distribution as the beam was focused 7 µm below the sample surface and laterally scanned across the interface. We initiated the FDTD calculations in a water domain starting 7.5 µm above the expected focus (i.e., 500 nm away from the interface) and let the beam propagate through a volume of 24 × 24 × 12 µm 3 . We scanned the position of the interface between y = −5 µm and y = +5 µm by steps of 250 nm. We then exported the time-averaged intensity distribution in a 5 × 5 × 8 µm 3 region around the theoretical focal spot; the results are displayed in Fig. 1(d) and Visualization 1. Equivalent calculations performed using the ASR model and a flat incident wavefront on the objective are shown in Fig. 1(c) for comparison. FDTD calculations highlight the dramatic distortions experienced by the focal field distribution as the beam is scanned across the vertical interface, transforming from a diffraction-limited spot in water to a double-peaked distribution close to the interface, and finally, to an aberrant focus characterized by defocus and spherical aberration inside the quartz domain (Visualization 1). Such distortions occurring on micrometer scales in the sample are not accounted for by simple propagation models [ Fig. 1(a)] and are largely neglected in nonlinear microscopy studies. One can, however, anticipate that they have important consequences on image contrast, as we will now confirm.
To provide a ground truth to test our simulations, we recorded 2PEF and THG 3D images of a vertical interface between quartz and a rhodamine-water solution (see Fig. 2). In x z reprojections, both 2PEF [ Fig. 2(b)] and THG [ Fig. 2(c)] images exhibit nontrivial artifacts in the vicinity of the interface, reminiscent of the beam distortions shown in Fig. 1. We then investigated whether FDTD methods can account for these nonlinear responses, and first simulated a 2PEF process, assuming a homogeneous concentration of fluorophores in water and no fluorescence in glass. We computed the squared intensity in the fluorescent medium, and integrated these values over the simulation volume to estimate the 2PEF signal. We scanned a region of x ∈ [−6 µm; 6 µm] on either side of the vertical interface and z ∈ [−4 µm; 12 µm] in depth around the horizontal interface [see Fig. 2(b) for the case of a linear incident polarization and Fig.S3 for the case of a circular polarization] to obtain a theoretical image. As a reference, we also performed numerical simulations based on the ASR/Green formalism, neglecting index heterogeneities near-focus. Figure 2(b) shows both experimental and numerical data using the same look-up table normalized to the maximum image intensities. The FDTD simulation is remarkably similar to the experimental data, and accurately reproduces nontrivial artifacts such as the axial decrease in intensity near the interface and the lateral shift of the fluorescence maximum with respect to the interface. In contrast, these experimentally observed artifacts are not reproduced when using the ASR/Green formalism [ Fig. 2(b), right].
We then adapted the simulation to the case of coherent nonlinear processes such as THG [ Fig. 2(c)] and SHG (Fig. S4 of Supplement 1). We performed a FDTD calculation of a THG x z image involving two nonlinear materials with χ (3) values taken to be 1.68 × 10 −5 for water and 2 × 10 −5 for quartz in order to keep the ratio between the two materials consistent with experimentally reported nonlinear susceptibilities [43]. We monitored the spectrum in a square array of positions at the end of the simulation volume and placed the numerical detectors within a linear material larger than the wavelength to avoid near-field effects (see Fig. S1 of Supplement 1 for detailed schematics of the simulation). Figure 2(c) shows experimental THG images of the interface along with the corresponding ASR/Green and FDTD simulations. While the ASR/Green function simulation predicts a simple single-peaked signal at both horizontal and vertical interfaces, the FDTD simulation accurately reproduces the experimentally observed lateral shift of the THG maximum outside the quartz material, as well as the emergence of a second lateral peak starting at a depth z ≈ 5 µm along the vertical interface.
We then investigated the relevance of FDTD strategies in the context of polarization-resolved nonlinear microscopy. The polarization properties of SHG, THG, or CARS signals are increasingly used to probe material properties at the micrometer scale [15,[19][20][21]. This approach, however, requires appropriate models to accurately interpret the measurements. We therefore analyzed the influence of the incoming polarization on THG signals arising from a vertical interface. As a reminder and reference, past theoretical studies using ASR/Green formalism have predicted that the THG signal from a vertical interface between isotropic media should exhibit 10−15% intensity variations when the incident linear polarization is rotated in the xy plane [12]. We revisited our quartz-water interface experiment by recording polarization-dependent THG images at successive depths along the interface. Following [20], we performed a simplified analysis of the P-stacks to extract maps of the P-THG modulation amplitude ("P -THG modulation" defined as (I max − I min )/I max ) and of the polarization angle resulting in maximum THG ("P-THG angle"). Results are shown in Fig. 3(a), along with corresponding FDTD and ASR/Green calculations. Since our current FDTD implementation enables only diagonal nonlinear tensor elements, we simulated the nonlinear response of the interface with an incident polarization successively along the x axis and the y axis, from which we estimated the modulation amplitude and the angle of maximum signal, assuming the maximum signal is along one of the two axes, in accordance with the experimental data. Fig. 3. P-THG response at the interface between two isotropic media. (a) Experimental geometry: a rotating waveplate installed before the objective is used to control the direction of the incident linear polarization in the x y plane. (b) x z projections of P-THG images from a vertical water-quartz interface. From top to bottom: THG intensity averaged over all incident polarizations; in-plane polarization angle providing strongest THG; P-THG modulation defined as the THG intensity variation when the incident polarization is turned between x and y . From left to right: experimental data; zoom on the region corresponding to the simulations; FDTD calculation; ASR/Green calculation. Scale bars =1 µm.
One first striking observation is that the experimentally measured P-THG modulation is ≈ 45% at all depths when integrated over a 1 µm distance. This value is much larger than predicted by the ASR/Green model, but accurately estimated by FDTD calculations. Moreover, while the ASR/Green model predicts that the maximum THG should be observed when the incident polarization is parallel to the vertical interface [see [12] and Fig. 3(a)], the experimentally observed THG is in fact maximized when the incident polarization is perpendicular to the interface, as properly predicted using the FDTD approach. This behavior most likely results from the index mismatch at both fundamental and harmonic wavelengths, resulting in modified contrast mechanisms. This strong polarization response at the interface between isotropic media should of course be properly modeled when interpreting polarization data in terms of sample optical properties.
We then analyzed a model sample more closely mimicking a biological sample. For that purpose, we recorded P-THG modulation and angle maps from oil droplets (≈ 6 µm radius) immobilized in a low concentration agarose gel (Fig. 4). The linear and nonlinear properties of water and olive oil are known [43], which enables quantitative comparison with the simulations. We consistently observed a maximum modulation in the 55 − 60% range (57 ± 6%, n = 5), with maximum THG signal recorded when the polarization is orthogonal to the droplet surface. Both phenomena are accurately reproduced by FDTD calculations and largely missed by ASR/Green calculations [ Fig. 4(b), last column].
Finally, to confirm the general relevance of our modeling approach for biological imaging, we recorded epidetected P-THG images of lipid droplets in anesthetized adult zebrafish. Figures 5  and S5 show representative images of 10 − 15 µm droplets located 50-100 µm under the skin in the dorsal region and in the head. These experiments were reminiscent of the measurements on model oil droplets: lipid droplets in vivo produced a P-THG modulation ranging between 30% and 60% and a radial angular profile with maximum THG observed for a polarization orthogonal to the droplet surface. Of note, the P-THG modulation observed on lipid droplets in vivo exhibited more heterogeneity than the model case of oil droplets in water considered in the previous figure. This may be attributed to various factors, including the different composition and optical properties of zebrafish lipid Fig. 4. P-THG from immobilized oil droplets in agarose. (a) Geometry of the experiment; (b) left: raw experimental THG from x -polarized and y -polarized excitation; right: experimental data and corresponding simulations using FDTD and ASR/Green models. The intensity maps shown (top row) correspond to the average THG intensity as the incident polarization is rotated in the x y plane, the amplitude maps (middle row) correspond to the fraction of the THG signal that depends on the polarization direction, and the angle maps (bottom row) correspond to the polarization direction yielding maximum THG signal. Scale bar = 5 µm; (c) profiles of THG intensity along the dotted lines shown in (b) for incident polarizations along the x and y directions and corresponding FDTD calculations showing quantitative agreement. droplets, their nonperfectly spherical shape and possible heterogeneity at the submicrometer scale, and the lower signal-to-noise ratio in live images where illumination was kept low to minimize phototoxicity. Nethertheless, the P-THG modulation amplitudes and angles observed in vivo are again consistent with FDTD calculations and cannot be reproduced with the ASR/Green model, assuming that these droplets are composed of an isotropic lipid-like medium surrounded by water (Fig. 5, left).

DISCUSSION
In summary, we have shown that a simulation strategy based on FDTD in the focal region can accurately describe the nonlinear microscopy response of a vertical interface between nonlinear materials with different linear properties. We chose a generic geometry and well-characterized materials (water, lipids, and glass with n ranging from 1.33 to 1.47) so that we could confirm the relevance of the approach for quantitative analyses of microscopy data. We have found that commonly used simulation strategies fail at reproducing artifacts experimentally observed on vertical interfaces imaged with 2PEF, THG, and P-THG contrasts. We estimate that this may be a major limitation of ASR/Green model when interpreting polarimetric microscopy data. Indeed, vertical interfaces between aqueous and lipid structures are present in a majority of cell and tissue microscopy images.
We point out that the geometries and optical index variations considered here are encountered in many biological samples at spatial scales similar to the ones considered here. Indeed, lipidic structures such as vesicles and myelin are an important source of THG contrast [10,44] in many organisms and typically exhibit sizes ranging between 0.1 and 20 µm, while connective tissues and bones [45] also exhibit large index mismatches (up to n = 1.55) at similar scales. Fortunately, we find that FDTD methods appear as a well-suited complementary numerical framework to investigate complex contrasts in this context, with the potential to enable us to extract more information from the images. For example, we anticipate that FDTD simulations will provide a way to better understand the contributions of both linear (index mismatch, birefringence) and nonlinear (χ (3) tensor symmetries) sources of contrast in polarization-resolved imaging of index-mismatched objects [18]. More generally, our findings have important implications for the analysis of fluorescence images recorded in the presence of local index mismatch and, most likely, for the exploitation of nonlinear techniques involving several excitation beams such as coherent Raman microscopy [31].