Third harmonic imaging contrast from tubular structures in the presence of index discontinuity

Accurate interpretation of third harmonic generation (THG) microscopy images in terms of sample optical properties and microstructure is generally hampered by the presence of excitation field distortions resulting from sample heterogeneity. Numerical methods that account for these artifacts need to be established. In this work, we experimentally and numerically analyze the THG contrast obtained from stretched hollow glass pipettes embedded in different liquids. We also characterize the nonlinear optical properties of 2,2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$'$$\end{document}′-thiodiethanol (TDE), a water-soluble index-matching medium. We find that index discontinuity not only changes the level and modulation amplitude of polarization-resolved THG signals, but can even change the polarization direction producing maximum THG near interfaces. We then show that a finite-difference time-domain (FDTD) modeling strategy can accurately account for contrast observed in optically heterogeneous samples, whereas reference Fourier-based numerical approaches are accurate only in the absence of index mismatch. This work opens perspectives for interpreting THG microscopy images of tubular objects and other geometries.

www.nature.com/scientificreports/ by analyzing the polarized THG (PTHG) contrast near a simple vertical index-mismatched interface 24 . As a further step towards understanding the THG contrast observed on optically complex structures such as myelinated axons 9,18,25 or biological cylindrical structures 26 , in the present work we analyze the model geometry of an in-plane tubular object made of isotropic material. Following 27,28 , we consider the case of a stretched hollow glass pipette in the presence of internal and external index liquids. We use experiments and numerical simulations to investigate how the THG and PTHG responses are governed by the combined effect of index mismatch and cylindrical geometry. Experimentally, we stretch borosilicate glass capillaries to a diameter about twenty microns, and use 2,2 ′ -thiodiethanol (TDE) as a water-soluble index-matching medium. TDE has been used in fluorescence microscopy studies 29 as a tunable-index medium to reduce spherical aberrations due to its high index (1.52) and solubility in water. We estimate the previously unknown third-order nonlinear susceptibility χ (3) of TDE by performing z-scan experiments, and we use this value to implement simulations based on the FDTD 24 and ASR/Green 15 formalisms. We find that the index mismatch has a strong influence on the PTHG contrast of vertical interfaces and on the relative THG visibility of horizontal and vertical interfaces. We show that a FDTD modeling approach can be used to accurately describe THG from cylindrical structures in the presence of index mismatch. This analysis also allows us to resolve some ambiguities and discrepancies between previous studies analyzing THG from cylinders and vertical interfaces 15,20,24,27 .

Material and methods
Pipettes preparation. We used borosilicate capillaries (Phymep, model GC100-10, ∅ ext = 1.0 mm, ∅ int = 0.58 mm) which are classically used to perform microinjections in biological samples such as zebrafish embryos. We pulled the capillaries with a micropipette puller equipped with a platinum heating coil (P-1000, Sutter Instrument). By combining high heating temperatures and low stretching forces 30 , we obtained thin elongated pipette tips with an inner diameter of approximately 12 µm . After removal of the obstructed section at the very end of the pipettes, we filled them with an index liquid ( n = 1.33 or 1.47) containing a fluorescent marker (Rhodamine B, Sigma). The tips of the stretched pipettes were then sealed by melting the glass with a heated resistance. Finally, we placed the pipettes horizontally in Petri dishes with bottom glass plates and immersed them in nonfluorescent index liquids ( n = 1.33 or 1.47). The index liquids used in this study were distilled water ( n = 1.33 ) and TDE diluted in water ( n = 1.47 ). The refractive indices of these solutions were controlled with a refractometer at = 589 nm. The fluorescent marker inside the pipettes was used to check for proper sealing of the pipette tips after immersion.
Pipettes imaging and data processing. Multiphoton imaging was performed on a custom-built upright microscope as described in 24 . Excitation was provided by a femtosecond laser source (80 MHz, 1100 nm, 100 fs, Insight X3, SpectraPhysics, USA) and a water immersion objective (25× , 1.05NA, Olympus, Japan). We imaged the pipette tips sequentially with THG and 2PEF contrasts. THG detection was performed in transmission with a high NA condenser and a UV filter (Semrock FF01-377/50). 2PEF was detected in the epi direction with a red filter (Semrock FF01-590/20). The signal level was kept in a range that avoided saturation of the photon-counting detection chain, i.e., less than one photon detected every four laser pulses. We performed two types of acquisitions for each combination of internal and external liquids. First, we acquired z-series of xy images, where z is the optical axis and the incident polarization is parallel to the pipette axis y. The xyz series were resliced to obtain cross-sectional (xz) profiles, and averaged along y. In addition, we acquired polarization-resolved xy series in the equatorial plane of the pipettes. Lateral and axial sampling were typically 0.23µm and 1.0µm . Polarimetric images were used to extract PTHG profiles by averaging 5 consecutive cross-sectional profiles along the pipette axis. The PTHG modulation was extracted from the polarization series using a unidirectional FFT analysis 31 . χ (3) measurements. We estimated the third-order nonlinear susceptibility χ (3) of TDE using the classical z-scan method [32][33][34] . In this approach, horizontal interfaces between materials of known (water/glass) and unknown (e.g. TDE/glass) nonlinear susceptibilities are successively z-scanned across the focus, and the measured amplitudes of the THG maxima are used to infer the χ (3) properties of the unknown medium. Samples were prepared by sandwiching a drop of liquid (water, TDE, or diluted TDE) between a microscope slide and a borosilicate glass coverslip (N-BK7, Schott). We then recorded THG profiles by scanning the samples along the z-direction from top to bottom with a step of �z = 0.5 µm . The focal volume successively scanned immersion water, the glass coverslip, and the liquid sample. THG z-profiles therefore displayed two peaks: a first peak with a maximum intensity I water/glass was observed at the water/glass interface, and a second peak with a maximum intensity I glass/liquid was observed at the glass/liquid interface. The ratio between these two peak intensities I glass/liquid /I water/glass referred as I norm glass/liquid was used to determine χ where I norm glass/air was estimated from a z-scan from a sample where air was substituted to an index liquid solution. χ (3) glass was assessed by applying the previous formula to intensities retrieved from a glass-water z-scan experiment and using the value for χ (3) water established in the literature (See "Results" section). Simulations using FDTD model. We performed FDTD simulations [35][36][37] of THG from tubular structures immersed in water and TDE by elaborating on the approach described in 24 . For the simulations, we neglected (1) www.nature.com/scientificreports/ dispersion by keeping the linear indices constant and we used the χ (3) values listed in Table 1, based on reference values for glass and water, and our experimental characterization of TDE properties. We worked with a commercial FDTD implementation (Lumerical Device suite, Ansys Inc, Canonsburg, PA) to facilitate future comparisons. In this implementation, the electric and magnetic fields are calculated on every point of a 3D grid at successive times by solving discretized Maxwell equations (Eq. 2) for non-magnetic materials The simulations are performed in the time domain, and spectral information can be retrieved using Fourier transforms in a post-processing step.
We used FDTD calculations to estimate THG as a focused beam is laterally scanned in the equatorial plane of cylinder with a 12µ m inner and 20µ m outer diameter. 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 24 × 24 × 16 µm 3 discretized over 20 − 40 nm steps with a time resolution of ≈ 0.5 fs . In the Lumerical software formalism, a nonlinear polarization term is introduced explicitly as follows (Eq. 3): One limitation of the current version of the software we used is that nonlinear materials can exhibit only diagonal non-zero tensor elements. In order to simulate isotropic materials, we therefore set χ 0 and only considered excitations polarized along x or y . Since storing the entire simulated fields would require several tens of gigabytes per simulation, we instead 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 ). Figure 1.a summarizes the simulation geometry and parameters. FDTD simulations were performed on a Dell Precision 7920 (2 × Intel 6140 CPUs, 384 GB RAM, DDR4 2666 MHz) running Lumerical version 2019b, 2020a, or 2021R2.4 with a typical computation time of ≈ 2 hours per condition. Table 1. n and χ (3) values in the simulations. (ii) the spatially-dependent induced nonlinear polarization is calculated in the focal region; (iii) this polarization is used as a source term and in turn propagated to the detector plane (see 15,21  www.nature.com/scientificreports/ Simulations using ASR/green model. For comparison, we performed ASR/Green's function calculations of THG from the same geometries using the method described in previous studies 15,[20][21][22]34 . This simulation strategy has been used to simulate nonlinear microscopy in the absence of aberrations. It relies on an angular spectrum representation (ASR) to calculate the excitation field distribution near the focus, assuming the same linear index everywhere 38,39, and on Green's functions 39,40 to propagate the nonlinear response from the focal region to the detector plane (again, assuming constant linear indices). Figure 1b summarizes the simulation geometry and parameters.

Results
THG from interfaces as a function of orientation and index mismatch. We first investigated the THG visibility of internal, external, horizontal and vertical interfaces as a function of index mismatch. We used TDE as a water-soluble index medium and diluted it (26%) to approximately match the refractive index of borosilicate glass (1.47 at 589 nm). We recorded xz sections of pipettes in four different conditions where the internal and external liquids were either index-mismatched (water) or index-matched (diluted TDE) with the borosilicate pipettes. Rhodamine was diluted in the internal liquid in order to verify that the pipette tips were efficiently sealed, and to easily identify the internal part of the pipettes. The experimental geometry and typical images are shown in Fig. 2. An obvious phenomenon is that the visibility of the different interfaces depends on both their orientation and the index mismatch. More specifically, we observed that: • Glass-TDE interfaces efficiently produce THG despite index matching. This suggests that TDE and glass have different nonlinear susceptibility χ (3) . In such case, standard models of THG microscopy near interfaces 15 predict that the observed signal scales as (χ • Glass-TDE horizontal interfaces are more visible than glass-water horizontal interfaces. This suggests that TDE and water also have different χ (3) . www.nature.com/scientificreports/ nonlinear susceptibilities of the two media: THG ∝ �(χ (3) ) 2 . These measurements also suggest that TDE could be used as a contrast agent in THG imaging.
Modelling THG contrast in the equatorial plane. Using this measured value of the third-order nonlinear susceptibility of TDE, we implemented FDTD simulations of the THG contrast as the excitation focus is laterally scanned in the equatorial plane of a 20 µm diameter glass pipette immersed either in TDE or in water. We also performed corresponding ASR/Green model-based calculations of the same geometry. Figure 4 presents experimental and simulated profiles obtained in the index-mismatched (glass-water) and index-matched (glass-TDE) cases. The index-mismatched case is of particular interest, because focus distortions can produce non-trivial THG contrast, in particular near interfaces 24 . One first observation is that internal vertical interfaces experimentally produce 3-to-5 times more THG than external interfaces (Fig. 4a). This phenomenon is likely due to the fact that the scanned beam is more aberrated when focused on the outer edge of the pipette in its equatorial plane, resulting in reduced excitation intensity. Interestingly, this THG profile is reproduced by the FDTDbased calculation (Fig. 4b). This confirms that FDTD is a valid formalism for reproducing the propagation of the excitation field in an index-mismatched cylindrical geometry. In contrast, ASR/Green calculations predict that inner and outer THG peaks should have similar amplitudes (Fig. 4c), in contradiction with the experimental observation. We also note that the shape of the peaks predicted by ASR/Green calculations is distorted when using realistic values for χ (3) . This is a numerical artifact due to the fact that the signal scales as the difference between the squared χ (3) ( �χ (3) ) while the background signal due to the finite integration volume scales as the average value of the squared χ (3) , so realistic values add an interference term between these two contributions. This artifact is classically avoided in published studies 15,20 by numerically increasing the χ (3) contrast and setting χ (3) = 0 or 1. As shown in Fig. 4d, this numerical trick removes the artifacts. However, this is at the cost of using unrealistic χ (3) values, which limits the quantitative nature of the numerical predictions. In the index-matched case of a pipette immersed and filled with TDE ( Fig. 4e-g), inner and outer THG peaks have similar amplitudes, and all simulation strategies produce accurate predictions.
PTHG contrast in the presence of index discontinuity. We then explored experimentally and numerically the polarization-resolved THG contrast in the equatorial plane of the pipettes in the two previous cases (Fig. 5): we recorded THG profiles with an excitation polarization in the directions perpendicular (x direction) and orthogonal (y direction) to the pipette axis, in index-mismatched (water) and index-matched (TDE) environments. The experimental observations are somehow intriguing: in line with previous reports 24 , indexmismatched vertical interfaces produced a strong PTHG modulation ( > 50% ), with maximum THG observed when the polarization is orthogonal to the interface (Fig. 5a,a2); in contrast, we observed a greatly reduced PTHG modulation ( < 8% ) near index-matched interfaces (Fig. 5d,g). This phenomenon is necessarily caused by polarization-dependent aberrations, and should of course be taken into account in the context of PTHG microscopy. Fortunately, this modulation dependence on index mismatch is reproduced by the FDTD calculations (Fig. 5b,b2,e,h). In contrast, the phenomenon is not reproduced by ASR/Green calculations (Fig. 5c,c2,f,i). www.nature.com/scientificreports/

PTHG modulation inversion on vertical interfaces.
Looking more carefully at the index-matched case, we note that both simulation models predict a limited dependence of THG on polarization (modulation amplitude ≤ 15 %) with maximum signal obtained for a polarization parallel to the interface 15 . In the following, we will speak of positive (resp. negative) PTHG modulation in cases where maximum THG is obtained for an incident polarization orthogonal (resp. parallel) to the interface. In experiments on index-matched pipettes, we did experimentally confirm a greatly reduced PTHG modulation near vertical interfaces, but we did not unambiguously measure the small negative modulation predicted by the simulations. This discrepancy may be due to non-perfect index matching, since we measured the optical indices at a different wavelength than the one use for imaging. Since PTHG modulation near a vertical interface appears to strongly depend on the optical contrast between the two materials, we investigated this effect more precisely. We revisited experimentally and numerically the case of a simple vertical interface 24 , this time with various values of the index mismatch. We measured PTHG modulation at the interface between glass and TDE solutions of increasing concentrations, with indices ranging from 1.33 (pure water) to 1.52 (pure TDE). For increasing concentrations of TDE, we measured (Fig. 6) PTHG modulation values ranging from largely positive (+45%) to null and even to increasingly negative (-25%) values. This previously undescribed phenomenon is well reproduced by the FDTD simulations for all experimental cases (Fig. 6). In the well-studied case of perfect index matching, a small PTHG negative modulation of approximately 10-15% is predicted by the classical ASR model 15 , which is due to the presence of the interface that breaks the circular symmetry of the problem. In the case of an index-mismatched interface, the physical interpretation is more complicated, as additional mechanisms alter THG efficiency. Most importantly, propagation along an index-mismatched vertical interface can severely distort the excitation field distribution at the sub-wavelength scale, and such aberrations are polarization-dependent. As an illustration, Fig. 6 shows polarization-resolved views of the distorted excitation field distribution in the xy and xz planes, for both matched ( n 1 = n 2 = 1.458 ) and mismatched ( n 1 = 1.458, n 2 = 1.33 ) cases. Aberrations introduce a polarization-dependent field asymmetry near the focal region, and therefore a polarization-dependent source of THG signal. THG efficiency is then governed by the interplay between this altered field distribution and the spatial distribution of nonlinear susceptibility χ (3) of the sample. In the most general case, PTHG modulation near a vertical interface is possibly a function of linear index mismatch, χ (3) distribution, as well as materials anisotropy 41 near focus. As indicated by our measurements and simulations, this complexity can in principle be investigated with advanced numerical methods. Finally, double lateral THG peaks were occasionally observed near index-mismatched vertical interfaces (Figs. 2b, 5a2), and were properly predicted by FDTD simulations (Fig. 5b2), as discussed in a previous study 24 .
Relative visibility of horizontal and vertical interfaces. Finally, we analyzed the relative visibility of horizontal (H) and vertical (V) interfaces, summarized in Fig. 7. Experiments (Fig. 2c) showed that the visibility www.nature.com/scientificreports/ ratio between H-interfaces and V-interfaces is low ( H/V = 0.1 ± 0.01 ) in index-mismatched cases (glass/water interfaces), and gets close to 1 ( H/V ≈ 0.5 − 1 ) in index-matched cases (glass/TDE). We ran FDTD and ASR/ Green calculations of THG by horizontal and vertical interfaces in a tubular geometry Fig. 7, averaged for x− and y− polarizations. We found that FDTD calculations properly reproduced the index-mismatch (glass/water) experiment ( H/V = 0.1 ± 0.02 ), in contrast with the ASR/Green simulations ( H/V ≈ 1 ± 0.2 ). In the indexmatched case (pipettes immersed in TDE), both models produced reasonable predictions ( H/V ≈ 0.7 − 9 ). The higher relative visibility of vertical index-mismatched interfaces could be partly due to the fact that this geometry results in a severe distortion of the focused field distribution (see e.g. Fig. 6e, in turn creating an effective symmetry breaking resulting in more efficient THG (i.e. less efficient Gouy shift-induced destructive interference). In a sense, THG microscopy may be viewed as a aberration detector in which signal can be enhanced by focus distortions, at least up to a certain amount of aberration. Our data indicate that this phenomenon should be tractable with FDTD calculations.

Discussion
In this work, we have shown that an FDTD modeling strategy can accurately account for the THG contrast observed in optically heterogeneous structures, while the faster classical Green function-based simulation methods are reliable only in the absence of index mismatch. This analysis opens perspectives for the analysis of the interplay between sample geometry 42 , excitation NA 19 , beam profile 20,21 , and field distortion induced by index heterogeneity 23,24 in THG microscopy -and ultimately to develop more quantitative interpretations of THG images, with the possibility of extracting submicron structural parameters. Our data show that the polarized THG contrast near interfaces is particularly sensitive to index mismatch. An interesting observation is that  , f). Bottom row, PTHG modulation amplitudes measured from the experimental data (g), FDTD (h) and ASR/Green simulations (i). Index mismatch results experimentally in a strong PTHG modulation (50-60% with more THG for an incident polarization perpendicular to the interface) and in lateral THG peak doubling near vertical interfaces (a2, g). These phenomena are predicted by FDTD calculations (b2, h) and not by the ASR/Green model (c2, i).
Both models perform well in the case of an index-matched samples. www.nature.com/scientificreports/  www.nature.com/scientificreports/ index-mismatched vertical glass-water interfaces, often encountered in microscopy experiments, produce a strong PTHG modulation ( > 50 %) with maximum THG observed when the polarization is orthogonal to the interface, whereas index-matched interfaces produce a weak inverted modulation with maximum THG when the polarization is parallel to the interface, as originally predicted by Green-based simulations 15 . This phenomenon is caused by polarization-dependent field-distributions, and can lead to misinterpretations if the images are not analyzed with an appropriate model. This is a particularly important observation in the context of PTHG microscopy. Indeed, it should be emphasized that these strong PTHG modulations are demonstrated here in the case of isotropic materials, and they can be stronger than the contribution due to materials anisotropy (linear and nonlinear). The interpretation of some PTHG 41,43 and possibly PCARS [44][45][46] experiments in optically heterogeneous samples should therefore probably be reconsidered in light of these findings. Fortunately, our data also show that these subtle effects can be analyzed with numerical FDTD approaches, even in the case of a complex sample geometry such as a double cylinder. These result pave the way for accurate modeling of PTHG from arbitrary objects. In particular, the numerical approach presented here could be extended to account for optical anisotropy, which would be appropriate for the analysis of THG from organized materials such as ordered lipid assemblies 41,47 , corneal stroma 48,49 , myelin 9,18,25 , or biocrystals 43,50 . www.nature.com/scientificreports/