## Abstract

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.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

Multiphoton microscopy [1] is widely used for imaging thick biological samples. Beyond its unique performance for in-depth imaging with $\unicode{x00B5}{\rm 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–5], second- and third-harmonic generation (SHG, THG) [6–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–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–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 (${\gt}100\;\unicode{x00B5}{\rm 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–28], and more recently to calculate aberrations induced when light propagates through bone [29] or 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 FDTD-based approach accurately accounts for experimentally observed artifacts, with important consequences for the interpretation of coherent and polarization-resolved images.

## 2. METHODS

#### 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 ($\overrightarrow B = {\mu _0}\overrightarrow H ;\overrightarrow D = {\varepsilon _0}{\varepsilon _r}\overrightarrow E$). The simulations are performed in the time domain, but spectral information can be retrieved using Fourier transforms in a postprocessing step,

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\;\unicode{x00B5}{\rm m}$ and 50 nm bandwidth, and performed calculations over a focal region spanning $18 \times 18 \times 13\;\unicode{x00B5} {{\rm m}^3}$ discretized over 20–40 nm steps with a time resolution of $\approx 0.2\;\rm 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, $({\chi _1^{(1)},\chi _1^{(3)}}),({\chi _2^{(1)},\chi _2^{(3)}})$. In the Lumerical software formalism, a nonlinear polarization term is introduced explicitly as follows: ${P_i}(t) = {\varepsilon _0}\chi _{\textit{ii}}^{(1)}{E_i}(t) + {\varepsilon _0}\chi _{\textit{iii}}^{(2)}E_i^2(t) + {\varepsilon _0}\chi _{\textit{iiii}}^{(3)}E_i^3(t)$ ($i \in \{x,y,z\}$). $\overrightarrow D$ is then defined as $\overrightarrow D (\omega) = {\varepsilon _0}\overrightarrow E (\omega) + \overrightarrow P (\omega)$. 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 $\chi _{\textit{xxxx}}^{(3)} = \chi _{\textit{yyyy}}^{(3)} = \chi _{\textit{zzzz}}^{(3)} = \chi _0^{(3)}$ and only considered excitations polarized along ${\boldsymbol x}$ or ${\boldsymbol 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 \times$ 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 $\approx 1 \!-\! 2\,\rm 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 \times$, 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\;\unicode{x00B5}{\rm 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\;\unicode{x00B5}\rm g/mL$, Sigma-Aldrich). We covered the sample with a $170\;\unicode{x00B5}{\rm m}$ thick borosilicate glass coverslip (Menzel, Thermo Scientific Menzel) to let the fluorescent solution spread and surround the quartz slide. $230\;\unicode{x00B5}{\rm 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 $xy - \rm 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].

## 3. 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\;\unicode{x00B5}{\rm m}$ above the expected focus (i.e., 500 nm away from the interface) and let the beam propagate through a volume of $24 \times 24 \times 12\;\unicode{x00B5} {{\rm m}^3}$. We scanned the position of the interface between $y = - 5\;\unicode{x00B5}{\rm m}$ and $y = + 5\;\unicode{x00B5}{\rm m}$ by steps of 250 nm. We then exported the time-averaged intensity distribution in a $5 \times 5 \times 8\;\unicode{x00B5} {{\rm 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 $xz$ 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 \in [- 6\,\,\unicode{x00B5}{\rm m};6\,\,\unicode{x00B5}{\rm m}]$ on either side of the vertical interface and $z \in [- 4\,\,\unicode{x00B5}{\rm m};12\,\,\unicode{x00B5}{\rm 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 $xz$ image involving two nonlinear materials with ${\chi ^{(3)}}$ values taken to be $1.68 \times {10^{- 5}}$ for water and $2 \times {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 \approx 5\,\,\unicode{x00B5}{\rm 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–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.

One first striking observation is that the experimentally measured P-THG modulation is $\approx 45\%$ at all depths when integrated over a $1\;\unicode{x00B5}{\rm 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 ($\approx 6\,\,\unicode{x00B5}{\rm 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 \pm 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\,\,\unicode{x00B5}{\rm 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 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).

## 4. 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 (${\chi ^{(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].

## Funding

Ligue Contre le Cancer; Agence Nationale de la Recherche (ANR-10-INBS-04, ANR-10-LABX-0073, ANR-15-CE11-0012, ANR-EQPX-0029).

## Acknowledgment

We thank all members of the advanced microscopies group at the Laboratory for Optics and Biosciences and of the zebrafish neurogenetics unit for scientific discussions on multiphoton microscopy. In particular, we thank Marie-Claire Schanne-Klein for proofreading the paper, Emilie Menant for zebrafish husbandry, and Laure Bally-Cuif and Sebastien Bedu for advice on fish mounting.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

The data (experimental and code) are available at [34].

## Supplemental document

See Supplement 1 for supporting content.

## REFERENCES

**1. **W. R. Zipfel, R. M. Williams, and W. W. Webb, “Nonlinear magic: multiphoton
microscopy in the biosciences,” Nat.
Biotechnol. **21**,
1369–1377 (2003). [CrossRef]

**2. **W. Denk, J. H. Strickler, and W. W. Webb, “Two-photon laser scanning
fluorescence microscopy,” Science **248**, 73–76
(1990). [CrossRef]

**3. **P. T. So, C. Y. Dong, B. R. Masters, and K. M. Berland, “Two-photon excitation
fluorescence microscopy,” Annu. Rev. Biomed.
Eng. **2**,
399–429 (2000). [CrossRef]

**4. **N. G. Horton, K. Wang, D. Kobat, C. G. Clark, F. W. Wise, C. B. Schaffer, and C. Xu, “In vivo three-photon
microscopy of subcortical structures within an intact mouse
brain,” Nat. Photonics **7**, 205–209
(2013). [CrossRef]

**5. **K. Guesmi, L. Abdeladim, S. Tozer, P. Mahou, T. Kumamoto, K. Jurkus, P. Rigaud, K. Loulier, N. Dray, P. Georges, M. Hanna, J. Livet, W. Supatto, E. Beaurepaire, and F. Druon, “Dual-color deep-tissue
three-photon microscopy with a multiband infrared
laser,” Light Sci. Appl. **7**, 12 (2018). [CrossRef]

**6. **F. S. Pavone and P. J. Campagnola, *Second Harmonic Generation
Imaging* (CRC Press,
2016).

**7. **Y. Barad, H. Eisenberg, M. Horowitz, and Y. Silberberg, “Nonlinear scanning laser
microscopy by third harmonic generation,”
Appl. Phys. Lett. **70**,
922–924 (1997). [CrossRef]

**8. **M. Müller, J. Squier, K. Wilson, and G. Brakenhoff, “3D microscopy of transparent
objects using third-harmonic generation,” J.
Microsc. **191**,
266–274 (1998). [CrossRef]

**9. **D. Yelin and Y. Silberberg, “Laser scanning
third-harmonic-generation microscopy in biology,”
Opt. Express **5**,
169–175 (1999). [CrossRef]

**10. **D. Débarre, W. Supatto, A.-M. Pena, A. Fabre, T. Tordjmann, L. Combettes, M.-C. Schanne-Klein, and E. Beaurepaire, “Imaging lipid bodies in cells
and tissues using third-harmonic generation
microscopy,” Nat. Methods **3**, 47–53
(2006). [CrossRef]

**11. **J.-X. Cheng and X. S. Xie, “Coherent anti-stokes Raman
scattering microscopy: instrumentation, theory, and
applications,” J. Phys. Chem. B **108**, 827–840
(2004). [CrossRef]

**12. **J.-X. Cheng and X. S. Xie, “Green’s function formulation
for third-harmonic generation microscopy,” J.
Opt. Soc. Am. B **19**,
1604–1610 (2002). [CrossRef]

**13. **J.-X. Cheng, A. Volkmer, and X. S. Xie, “Theoretical and experimental
characterization of coherent anti-stokes Raman scattering
microscopy,” J. Opt. Soc. Am. B **19**, 1363–1375
(2002). [CrossRef]

**14. **N. Olivier, D. Débarre, P. Mahou, and E. Beaurepaire, “Third-harmonic generation
microscopy with Bessel beams: a numerical study,”
Opt. Express **20**,
24886–24902 (2012). [CrossRef]

**15. **I. Gusachenko and M.-C. Schanne-Klein, “Numerical simulation of
polarization-resolved second-harmonic microscopy in birefringent
media,” Phys. Rev. A **88**, 053811 (2013). [CrossRef]

**16. **R. S. Pillai, G. Brakenhoff, and M. Müller, “Analysis of the influence of
spherical aberration from focusing through a dielectric slab in
quantitative nonlinear optical susceptibility measurements using
third-harmonic generation,” Opt.
Express **14**,
260–269 (2006). [CrossRef]

**17. **D. Sandkuijl, A. E. Tuer, D. Tokarz, J. Sipe, and V. Barzda, “Numerical second-and
third-harmonic generation microscopy,” J. Opt.
Soc. Am. B **30**,
382–395 (2013). [CrossRef]

**18. **M. Zimmerley, P. Mahou, D. Débarre, M.-C. Schanne-Klein, and E. Beaurepaire, “Probing ordered lipid
assemblies with polarized third-harmonic-generation
microscopy,” Phys. Rev. X **3**, 011002 (2013). [CrossRef]

**19. **R. Turcotte, D. J. Rutledge, E. Bélanger, D. Dill, W. B. Macklin, and D. C. Côté, “Intravital assessment of
myelin molecular order with polarimetric multiphoton
microscopy,” Sci. Rep. **6**, 31685 (2016). [CrossRef]

**20. **J. Morizet, G. Ducourthial, W. Supatto, A. Boutillon, R. Legouis, M.-C. Schanne-Klein, C. Stringari, and E. Beaurepaire, “High-speed
polarization-resolved third-harmonic microscopy,”
Optica **6**,
385–388 (2019). [CrossRef]

**21. **P. Gasecka, A. Jaouen, F.-Z. Bioud, H. B. De Aguiar, J. Duboisset, P. Ferrand, H. Rigneault, N. K. Balla, F. Debarbieux, and S. Brasselet, “Lipid order degradation in
autoimmune demyelination probed by polarized coherent Raman
microscopy,” Biophys. J. **113**, 1520–1530
(2017). [CrossRef]

**22. **Q. Fang and S. Yan, “Graphics processing
unit-accelerated mesh-based Monte Carlo photon transport
simulations,” J. Biomed. Opt. **24**, 115002 (2019). [CrossRef]

**23. **T. Young-Schultz, S. Brown, L. Lilge, and V. Betz, “FullMonteCUDA: a fast,
flexible, and accurate GPU-accelerated Monte Carlo simulator for light
propagation in turbid media,” Biomed. Opt.
Express **10**,
4711–4726 (2019). [CrossRef]

**24. **X. Cheng, Y. Li, J. Mertz, S. Sakadžić, A. Devor, D. A. Boas, and L. Tian, “Development of a beam
propagation method to simulate the point spread function degradation
in scattering media,” Opt. Lett. **44**, 4989–4992
(2019). [CrossRef]

**25. **M. Weigert, K. Subramanian, S. T. Bundschuh, E. W. Myers, and M. Kreysing, “Biobeam–multiplexed
wave-optical simulations of light-sheet microscopy,”
PLoS Comput. Biol. **14**,
e1006079 (2018). [CrossRef]

**26. **K. Choi, J. W. Chon, M. Gu, and B. Lee, “Characterization of a
subwavelength-scale 3D void structure using the FDTD-based confocal
laser scanning microscopic image mapping technique,”
Opt. Express **15**,
10767–10781 (2007). [CrossRef]

**27. **S. Tanev, J. Pond, P. Paddon, and V. V. Tuchin, “Optical phase contrast
microscope imaging: a FDTD modeling approach,”
Proc. SPIE **6791**,
67910E (2008). [CrossRef]

**28. **P. Török, P. Munro, and E. E. Kriezis, “High numerical aperture
vectorial imaging in coherent optical microscopes,”
Opt. Express **16**,
507–523 (2008). [CrossRef]

**29. **K. F. Tehrani, P. Kner, and L. J. Mortensen, “Characterization of wavefront
errors in mouse cranial bone using second-harmonic
generation,” J. Biomed. Opt. **22**, 036012 (2017). [CrossRef]

**30. **M. Menzel, M. Axer, H. De Raedt, I. Costantini, L. Silvestri, F. S. Pavone, K. Amunts, and K. Michielsen, “Toward a high-resolution
reconstruction of 3D nerve fiber architectures and crossings in the
brain using light scattering measurements and finite-difference
time-domain simulations,” Phys. Rev.
X **10**, 021002
(2020). [CrossRef]

**31. **J. van der Kolk, A. C. Lesina, and L. Ramunno, “Effects of refractive index
mismatch on SRS and CARS microscopy,” Opt.
Express **24**,
25752–25766 (2016). [CrossRef]

**32. **J. N. van der Kolk, S. Bancelin, C. Kioulos, A. C. Lesina, F. Légaré, and L. Ramunno, “Effect of refractive index
mismatch on forward-to-backward ratios in SHG
imaging,” Opt. Lett. **43**, 5082–5085
(2018). [CrossRef]

**33. **C. Liu, Z. Huang, F. Lu, W. Zheng, D. W. Hutmacher, and C. Sheppard, “Near-field effects on coherent
anti-stokes Raman scattering microscopy imaging,”
Opt. Express **15**,
4118–4131 (2007). [CrossRef]

**34. **J. Morizet, G. Sartorello, N. Dray, C. Stringari, E. Beaurepaire, and N. Olivier, “FDTD for NL
microscopy,” Zenodo (2021),
http://doi.org/10.5281/zenodo.4722857.

**35. **A. Taflove and S. C. Hagness, *Computational Electromagnetics: The
Finite-Difference Time-Domain Method* (Artech
House, 1995),
pp. 149–161.

**36. **Í. R. Çapoğlu, J. D. Rogers, A. Taflove, and V. Backman, “The microscope in a computer:
image synthesis from three-dimensional full-vector solutions of
Maxwell’s equations at the nanometer scale,” in
*Progress in Optics*
(Elsevier, 2012),
Vol. 57,
pp. 1–91.

**37. **B. Gallinet, J. Butet, and O. J. Martin, “Numerical methods for
nanophotonics: standard problems and future
challenges,” Laser Photon. Rev. **9**, 577–603
(2015). [CrossRef]

**38. **J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid, J.-Y. Tinevez, D. J. White, V. Hartenstein, K. Eliceiri, P. Tomancak, and A. Cardona, “Fiji: an open-source platform
for biological-image analysis,” Nat.
Methods **9**,
676–682
(2012). [CrossRef]

**39. **L. A. Royer, M. Weigert, U. Günther, N. Maghelli, F. Jug, I. F. Sbalzarini, and E. W. Myers, “ClearVolume: open-source live
3D visualization for light-sheet microscopy,”
Nat. Methods **12**,
480–481
(2015). [CrossRef]

**40. **N. Olivier and E. Beaurepaire, “Third-harmonic generation
microscopy with focus-engineered beams: a numerical
study,” Opt. Express **16**, 14703–14715
(2008). [CrossRef]

**41. **N. Dray, S. Bedu, N. Vuillemin, A. Alunni, M. Coolen, M. Krecsmarik, W. Supatto, E. Beaurepaire, and L. Bally-Cuif, “Large-scale live imaging of
adult neural stem cells in their endogenous niche,”
Development **142**,
3592–3600
(2015). [CrossRef]

**42. **R. M. White, A. Sessa, C. Burke, T. Bowman, J. LeBlanc, C. Ceol, C. Bourque, M. Dovey, W. Goessling, C. E. Burns, and L. I. Zon, “Transparent adult zebrafish as
a tool for in vivo transplantation analysis,”
Cell Stem Cell **2**,
183–189
(2008). [CrossRef]

**43. **D. Débarre and E. Beaurepaire, “Quantitative characterization
of biological liquids for third-harmonic generation
microscopy,” Biophys. J. **92**, 603–612
(2007). [CrossRef]

**44. **M. J. Farrar, F. W. Wise, J. R. Fetcho, and C. B. Schaffer, “In vivo imaging of myelin in
the vertebrate central nervous system using third harmonic generation
microscopy,” Biophys. J. **100**, 1362–1371
(2011). [CrossRef]

**45. **R. Genthial, E. Beaurepaire, M.-C. Schanne-Klein, F. Peyrin, D. Farlay, C. Olivier, Y. Bala, G. Boivin, J.-C. Vial, D. Débarre, and A. Gourrier, “Label-free imaging of bone
multiscale porosity and interfaces using third-harmonic generation
microscopy,” Sci. Rep. **7**, 3149 (2017). [CrossRef]