Quantifying the Coherent Interaction Length of Second-Harmonic Microscopy in Lithium Niobate Confined Nanostructures

Nonlinear optical (NLO) microscopy techniques such as second-harmonic generation (SHG), third-harmonic generation (THG), or coherent Anti-stokes Raman spectroscopy (CARS), see widespread use for both imaging and the spectroscopic analysis from cell biology to material sciences [1, 2, 3, 4, 5]. All these methods rely on the coherent NLO interaction of electromagnetic waves at a single or multiple fundamental frequencies, in order to generate a NLO signal at a different frequency. The NLO-generation efficiencies strongly depend on the underlying, local material properties. Conversely, they equally provide a label-free contrast mechanism for an in-depth analysis of distinct material properties using NLO imaging. For example, CARS spectroscopy is sensitive to molecular vibrations and phonons in solids, thus fingerprinting the local chemical composition or local structural changes [6, 7, 8]. Bulk SHG, in contrast, works only in non-centrosymmetric materials (e.g. piezoelectrics and ferroelectrics) and strongly depends on the material symmetry and incident light polarization [9, 10, 11, 12, 13]. SHG polarimetry thus witnesses a widespread use to analyze the orientation and crystallographic structure of appropriate (nano)systems. In this context, SHG might be even applied to study phase transitions, to visualize ferroelectric (FE) domains and domain walls (DWs) therein [14, 15, 16, 17, 18], or to analyze novel topologies[9, 10].

NLO effects per se are coherent optical processes. In the absence of dispersion, we witness a monotonic and nonlinear increase of the generated NLO signal as a function of interaction length. However, due to dispersion, each of the interacting waves propagates with different velocities, inherently resulting in an increasing phase difference between the fundamental and NLO signal beams. This limits the overall efficiency of any NLO process and is commonly referred to as the phase-matching problem. In this context, the coherent interaction length $l_{c}$ specifies exactly the length up to which a NLO process is still constructive [19]. Notably, the overall NLO signal amplitude will oscillate with twice the length $l_{c}$.

SHG is the archetypical NLO process where only two waves are involved, the fundamental (pump) beam at frequency $\omega_{inc}=\omega$ and the SHG wave $\omega_{SHG}=2\omega$. Phase matching of these two waves in $k$-space requires

with $\Delta k$ the phase mismatch. Here, two main processes may be differentiated (see Fig. 1), (i) co-propagation with $+k_{2\omega}$ being collinear, and (ii) counter-propagation with $-k_{2\omega}$ aligning anti-collinearly with respect to the fundamental wave vector $k$, respectively. Using $k=2{\pi}n_{\omega}/{\lambda_{\omega}}$, with $n_{\omega}$ the refractive index at the given $\lambda_{\omega}$, and ${\Delta}k={\pi}/l_{c}$, we obtain for $l_{c}$:

Note the positive $/$ negative sign in $k_{2\omega}$ in Eq. 2, standing for the counter- and co-propagating SH wave, respectively. To underline the importance of such a differentiation between (i) and (ii), a quick calculus for the lithium niobate (LN) sample used here, pumped at 800 nm (using a Ti:Sapphire laser) reveals the following two interaction lengths for co- and counter propagating waves: $l_{c,co}\approx 1280$ nm and $l_{c,counter}\approx 44$ nm. Since SHG intensities scale quadratically with interaction length, co-propagating SH light intensities are approximately three orders of magnitude larger as compared to the counter-propagating SH light intensities (assuming a sample size of exactly $l_{c}$). Similar arguments can be made for other NLO processes, such as THG or CARS.

Experimentally, co- and counter-propagating light in general is easily separated in SHG microscopy, due to their opposite $k$-vector behavior. Nevertheless, when mounting the sample under inspection on a substrate support (as is usually necessary when investigating thinner films or 2D materials)[3, 11], a significant amount of the co-propagating light can be reflected off the sample/substrate interface, and thus scatters into backward direction. Due to the potentially huge intensity difference between the co- and counter-propagating waves even a few percent of reflectivity can amount to a significant additional signal strength in backward-direction. Even more importantly, reflective interfaces may also resonantly enhance both the fundamental and/or the SH beams, hence additionally complicating experimental interpretation.

Note that the above concerns might have negligible impact when applying standard NLO microscopy for instance to cell biology investigations. Firstly, active layers might be much smaller in size than the coherent interaction length $l_{c}$, while secondly, only few and weakly reflective interfaces are present that might cause back-reflection or interface enhancement. In modern-type material science, however, a SHG-microscopy investigation of nanostructures and samples with confining structures is of utmost interest, hence automatically needing to address the above issues. Therefore, effects of reflection, resonant enhancement, and/or different phase matching scenarios, all need to be carefully taken into account when aiming at the quantitative analysis of the measured SHG contrast.

In the work presented here, we systematically investigate these effects by comparing our experimental findings from SHG microscopy to our rigorous numerical simulations. Both x- and z-cut LN wedges mounted on a Si and epoxy substrate, respectively, are analyzed, allowing us to elegantly change the sample thickness $t$ from tens of nanometers up to the bulk threshold, through simply scanning the incident linearly polarized photon beam across the wedges’ slope (see Fig. 1). Our study paves the way for the systematic understanding of both reflections and the role of phase matching when applying NLO microscopy to (ultra)thin films and nanostructures, i.e. ferroelectric and ferroic thin films, or even 2D materials. Due to the archetypical nature of the SHG process, these findings can readily be translated to any other NLO process.

We use single-crystalline 5-$\%$MgO-doped congruent LN for our SHG studies here, prepared into wedges as x- and z-cut thin films and mounted onto an epoxy and Si substrate, respectively (for details see Supplementary Information I.B.). LN is one of the most common NLO material and sees widespread use in electronics as surface-acoustic wave filters, high-frequency electro-optical modulators, or in quantum optics. Moreover, the current boost in implementing LN into CMOS-compatible architectures can undoubtably be attributed to thin-film lithium niobate (TFLN) that provides many benefits as compared to the standard bulk LN platform [20, 21, 22, 23, 24, 25, 26, 27]. LN in the form as prepared here into wedges, hence, provides an ideal test system for quantifying the coherent interaction length $l_{c}$ and therefore analyzing the phase-matching process as needed in many NLO applications.

SHG scattering is described through the induced SHG polarization $\vec{P}$, with $I_{SH}\propto\vec{P}^{2}$ the SH intensity. $\vec{P}$ is readily calculated through

with $E_{x}$, $E_{y}$, and $E_{z}$ denoting the electric field components at the fundamental wavelength along the crystal axes x,y,z, respectively, and $\chi^{(2)}$ being the second order NLO susceptibility tensor given for LN as:

Due to symmetry, only three components $d_{ij}$ are independent in LN, i.e. $d_{21}=d_{22}=d_{16}=3$ pm/V, $d_{31}=d_{32}=d_{24}=d_{15}=6$ pm/V, and $d_{33}=36$ pm/V, which can be assumed approximately constant over the relevant wavelength range from 850 to 950 nm [28, 29, 30]. Moreover, LN is an uniaxial birefringent material, and hence features two independent refractive indices, the ordinary and extraordinary refractive indices $n_{o}$ and $n_{e}$ for light polarized in the x-y plane and along the crystallographic z-axes, respectively [31]. Combining these two refractive indices with the SHG susceptibility tensor $d_{ij}$, we identify three distinct scenarios as listed in Tab. 1 and schematically depicted in Fig. 2, in order to investigate $l_{c}$ by SHG microscopy along the x- and z-cut LN wedges:

• i)

  Case I: When using z-cut TFLN, the linearly-polarized incident beam propagates along the z-axis [see Fig. 2(a)], with hence the backreflected SHG light stemming from the SHG coefficients $d_{21}$ or $d_{22}$ within the x-y-plane. Since both parameters yield similar results (see Supplementary Information Section II.), we consider here the case of y-polarized light, only. Due to strong dispersion, $l_{c}$ will be relatively small.

• ii)

  Case II: For x-cut TFLN with the incident x-aligned beam being polarized along the extraordinary z-axes [see Fig. 2(b)], the SHG response is exclusively generated through the largest tensor element $d_{33}$ and polarized along the extraordinary axes. We thus expect a slightly larger coherent interaction length as compared to case I.

• iii)

  Case III: Here, the incident beam hits the x-cut TFLN under ordinary polarization (y-direction) while SHG light leaves the crystal z-polarized (extraordinary polarization), as mediated via the $d_{32}$ element [see Fig. 2(c)]. Due to LN’s strong birefringence, an almost 3$\times$ larger $l_{c}$ results, that converges to infinity for pump wavelengths close to 1080 nm for 5%MgO:LN [32]. This is the well-known case of birefringence-phase matching in NLO, and hence widely exploited experimentally when fabricating highly-efficient NLO devices.

To study the coherence length for the three cases listed in Tab. 1, wedge-samples have been fabricated in collaboration with SGS Fresenius GmbH as shown in Fig. 1 from z- and x-cut LN, respectively. The z-cut sample was bonded to a (001)-silicon wafer, providing an enhanced stability that allowed wedge polishing to a smaller angle. For polishing, both samples were completely immersed in epoxy and polished until interference patterns became visible at the thinner end of the wedges. indicating sub-wavelength thicknesses. The thickness and angle of the wedges were determined via Confocal laser scanning microscopy at a 405 nm wavelength (Olympus LEXT OLS4000). The inclination angle $\alpha$ of the z- and x-cut sample was determined over the first 2 mm to $\alpha_{z}=0.45$° $\pm 0.01$° and $\alpha_{x}=1.97$° $\pm 0.01$°, respectively, while the wedge surface roughness evaluated from optical measurement was less than $<25$ nm_RMS. Details on the fabrication process and angle extraction can be found in the Supplemental Information sections I.B and I.C.

All SHG experiments were performed on a commercial microscope (Zeiss LSM980) equipped with a tunable Ti:Sapphire laser light source (see Supplemental Information section I.A). In our experiments, we employed two different objective lenses with numerical apertures NA = 0.45 and NA = 0.8. The incident pump light is always linearly polarized while the back-reflected SHG light contains all polarizations (unpolarized). Note, that each of the above discussed three cases is dominated by one major SHG process as denoted in Tab. 1, making polarization analysis obsolete. Moreover, this is also confirmed with our numerical calculations.

Numerical calculations based on a code previously developed by Sandkujil et al. were carried out to complement our experimental findings [33, 34, 11, 35, 36]. These simulations feature a full 3-dimensional vectorial model for both, the focusing and collection processes, as well as a matrix formalism in order to describe influences of (multiple) reflections and/or (self) interferences occurring for focal beams. The model hence even accounts for larger k-vector spreads when using high-NA objectives; for example, this allows to address and quantify the impact of mixed polarization components $d_{24}$ or $d_{15}$ on the SHG efficiency by solely using linearly-polarized incident light [35, 37]. We thus performed all calculations with a linear polarization input and using the full $\chi^{(2)}$-tensor. Details on our numerical simulations can be found in the Supplement (section I.D).

To analyze the thickness-dependent coherent interaction length $l_{c}$ of the SHG signal we carefully investigated the above three cases I, II, and III. Figure 3(a), (c), and (e) display large-scale (stitched) SHG overview images taken at a 850 nm fundamental wavelength. As the complete sample cannot be focused at the same time for each tile a z-stack was recorded and projected onto the 2D plan. In all these images, the wedge thickness $t$ linearly increases from left to right, as denoted by the sketch at the bottom of the figure. Note that Fig. 3(a) and (c) were recorded from the z-cut wedge with NA = 0.45 and NA = 0.80, respectively, while Fig. 3(e) shows the SHG response across the x-cut wedge at a NA = 0.80. In general, all images show characteristic oscillations in the SHG response, producing wavefronts that are aligned orthogonal to the thickness increase (i.e. the wedge slope). Notably, the periodicity measured for the long-wavelength oscillations is much shorter for the x-cut wedge ($\approx 70~{}\mu$m) as compared to the z-cut wedge ($\approx 250~{}\mu$m), a fact that solely stems from the roughly 4$\times$-larger wedge angle in case II, and not associated to different $l_{c}$-values.

When analyzing these images in more detail, several more observations can be made. Firstly, we notice for the x-cut wedge in Fig. 3(e) that oscillations are visible over more than a 2 mm width, with the contrast vanishing in the last third of the image to the right, only. A rough calculation allows to estimate the crystal thickness at that spot to measure approximately 70 $\mu$m, where reflections of co-propagating phase-matched light obviously influence the overall SH intensity. This value is more than 10$\times$ larger as compared to the axial optical resolution $\Delta z$ in LN, with the latter given to:

Secondly, apart from the long-period oscillation associated with the co-propagating phase matching process, short-range oscillations are visible. Most noticeably, they have wavelengths of approximately 20-25 $\mu$m and 5 $\mu$m for z- and x-cut LN wedges, respectively. Normalizing with respect to the inclination angles $\alpha_{z}$ and $\alpha_{x}$ yields thickness differences between neighboring intensity maxima in the range of 160-190 nm. These oscillations are associated with thin film interference effects of the fundamental and/or the SH signal, and also prominently appear in the numerical calculations as discussed below in further detail.
As is apparent from these cross sections here, the thin film resonances significantly alter the overall SH signal intensity by a factor of 2 or more, over thickness variations that measure just a few ten nm. Even more dramatic, these oscillations are visible up to thicknesses of at least several ten $\mu$m. This means that in particular for quantitative measurements of SHG intensities, reflections and the co-propagating process need to be carefully evaluated even up to a 50 $\mu$m crystal thickness and beyond. Moreover, interface reflectivity, thin film interference, and the coherent interaction length $l_{c}$ all are governed by the dispersion of the refractive index. Hence, they will also be affected when inspecting crystals of a constant thickness while varying the incident wavelength. These findings are excellently supported through our numerical simulations and documented in the Supplement Information section. III.D.

The third striking feature appears when analyzing the number of interference fringes, i.e. 12 and 15 oscillations, as counted in case I for the two different NAs used here, over the same crystal width [see Figs. 3(a) and (c)]. When strictly interpreting every oscillation period as twice the coherent interaction length $l_{c}$ as governed by Eq. 2, the one and only sound conclusion is that $l_{c}$ must depend on the NA of the objective lens in use; a larger NA thus must result in a shorter $l_{c}$, a fact that is absolutely not reflected by Eq. 2. Nevertheless, this - at first glance counter-intuitive - result can be readily explained when taking the nature of a focusing objective into account:
Eqs. 1 and 2 describe the phase matching condition for planar waves only, while focal beams always show a severe k-vector spread. In a simple picture this means that rays stemming from the outer rim within the focused beam, hit the sample to inspect under non-orthogonal conditions and hence take a longer effective trajectory when traveling through the crystal, as compared to central rays within the beam. Different parts of the focused beam, therefore, will be phasematched at different depths of the layer [11]. A strongly focusing objective (large NA) thus will observe an overall perceived shorter coherent interaction length, as parts of the beam will travel at a particularly larger angle through the material.

We have simulated this scenario for different NAs and display them in section III.C in the Supplemental Information. Note that a NA of $\approx$ 0.1 already very well reproduces the plane-wave condition as relevant to Eq. 2. As a conclusion, we thus recommend here to introduce the effective coherent interaction length $l_{c,eff}$ when dealing with focused beams, rather than solely using $l_{c}$ as being valid strictly speaking for planar waves, only.

To analyze $l_{c,eff}$ for the three different cases as a function of wavelength and/or NA, we performed detailed scans at the thinner end of each wedge. To avoid any unambiguity or error through stitching or rendering, we used here single frames only that at least contained the first full oscillation, as shown by the insets in Fig. 3(a), (c), and (e). Here, the focus was placed approximately at the center of the film at the thin edge of the wedge. As discussed in the supplement, the first oscillation period is least subjective to variations in focus positioning (See Supplemental Sec. IIIE). Moreover, lateral binning was applied in order to obtain the averaged intensity profiles, which is shown as the bold line in all these insets and displayed as cross sections in Figs. 3(b), (d), and (e).
For every case, the fundamental wavelength was swept from 850 nm to 950 nm via the microscope control software while the sample was kept at a constant location. Examples for as-obtained raw SHG data at the 850 nm fundamental wavelength, are shown in Figs. 3(b), (d), and (f). For comparison reasons, all such SH intensity data were normalized with respect to the inclination angles $\alpha_{z}$ and $\alpha_{x}$, and then displayed in Fig. 4 as a function of thickness rather than lateral width. What becomes crucial, however, is to exactly define the origin in all these graphs, i.e. the exact position of zero thickness for each SHG profile. To do so, we first identified the location of the first significant SHG intensity ($>0.1$ % of the maximum intensity). Then, profilometer scans across the wedge edge allowed us to quantify the residual step thickness at the thinner end of every LN wedge to approximately 50-100 nm. We therefore added to each position of significant first SHG intensity an offset of a 100 nm thickness, and conservatively assumed a confidence interval for the converted thicknesses of $\pm 100$ nm. Note, that the assumption of an initial step of 100 nm also fits well to our observations from numerical simulations, where the first local SHG maximum and hence the first significant SHG intensity as defined above, is typically observed for thicknesses of approximately 90 nm. This corresponds to the first maximum of the thin film interference of SHG light at 425 nm. A proof of that is seen from the excellent overlap of the measured and simulated data sets as illustrated later in Fig. 5. It should be noted that the confidence interval of $\pm 100$ nm corresponds to less than 5% of the total value of $2l_{c}$, which is always larger than 2.5 $\mu$m for all three cases I, II, and III. Therefore, any further errors in the measured data, e.g. imposed by the uncertainty in the wedge angles of $\pm 0.01$°, are negligible for wedge thicknesses of $<20$ $\mu$m, which are at most considered for detailed quantitative analysis here.

Figure 4 shows the obtained data sets as converted to the respective material thickness for all three cases. Each subfigure contains in total eleven datasets measured at fundamental wavelength from 850 - 950 nm in 10-nm increments. For optimal comparability, all datasets are normalized in their magnitude and displayed as a waterfall plot. Overall, for each case a distinct oscillation period is observed, that increases with increasing wavelength as expected from the decreasing dispersion.
In detail for Case I, one can see an increase of the position of the first minimum from $2l_{c,eff}\approx 2400$ nm to $2l_{c,eff}\approx 3550$ nm for the lower NA [Fig. 4(a)], while a slightly shorter period is observed for the larger NA [Fig. 4(b)] with the position of the first minimum shifting from $2l_{c,eff}\approx 2280$ nm to $2l_{c,eff}\approx 3460$ nm. In contrast, the corresponding periods for Case II [Fig. 4(c)] are somewhat lower as compared to Case I (from $2l_{c,eff}\approx 2220$ nm to $2l_{c,eff}\approx 3450$ nm for increasing wavelengths), although expected to be larger (see Tab. 1). As discussed later and also clearly underlined by the numerical calculations, the influence of focusing on the effective coherent interaction length is less pronounced in Case II as compared to Case I, which is likely due to the slightly lower value of the overall refractive index of the extraordinary index resulting in less refraction, and hence a longer path length for the outer rays in a focused beam.
Furthermore, we notice that the short-scale oscillation amplitudes related to the thin film interferences, are much lower for the x-cut sample (Cases II and III) as compared to the z-cut Case I. This is mostly due to the significantly different substrate reflectivity: The silicon substrate used to prepare the z-cut sample has a refractive index of $n=3.63$ to $n=5.02$ at 850 nm to 425 nm wavelengths, respectively, while the epoxy features a moderate refractive index of $n\approx 1.52$ without any pronounced dispersion over the relevant wavelength range. This results in much larger reflections and hence thin film interference effects for the z-cut wedge used in Case I.
Figure 4(d) shows the data related to Case III. Compared to the other cases a more complex shape with more than one distinct minimum is observed. The dominating minimum, for example, observed at $2l_{c,eff}\approx 6400$ nm for 850 nm shows a much longer coherent interaction length compared to the cases discussed so far, and belongs to Case III as described in Tab. 1. The position of this minimum shifts to higher wavelengths and can be identified up to $\lambda$ = 930 nm, where a length of $2l_{c,eff}\approx 13200$ nm is observed. Above that the contrast becomes too weak, while other processes may overlap with this effect. Apart from this minimum, a second set of oscillations is visible with minima observed, for example, at approximately 2.5 $\mu$m and 5 $\mu$m at 850 nm. They belong to a Case I process described by $P_{y}=d_{22}E_{y}^{2}$, and appear for this case here due to the unpolarized detection. The displayed data for Case III show that even when using unpolarized light, very complex SHG responses may be expected that need to be taken into account in every SHG experiment.

For quantitative comparison of the experimental data we have performed thorough numerical simulations. Figure 5 shows such an example comparing the numerically-calculated profiles to the experimental data of Case I with $NA$ = 0.8 at selected wavelengths. All datasets were normalized for maximum intensity and plotted with an artificial offset for better visibility. To allow for a detailed comparison of the overall shapes the SHG intensity was calculated in thickness increments of only 8 nm. Overall, the simulation reproduces the general shape, as well as the position, period, and magnitude of the short-scale and large-scale oscillations very well. These short-range oscillations belong to thin film interferences of either the fundamental and/or SHG beam as discussed in detail below. It should be noted that the experiment and simulated peaks show some slight deviations. For example at 950 nm, the double peaks visible around the overall maximum show a reversed intensity ratio in the experiment as compared to simulation. We believe that this results from the different reflectivity observed in experiment and simulation at the respective wavelength, due to the additional layer of glue present in the experiment, which is well backed up by numerical simulations as shown in the Supplemental Information section III.G. Furthermore, slight differences between the calculated refractive indices given via the Sellmeier equations [31, 38] and the real refractive indices of the sample, might play an additional role. Moreover, our numerical calculations were always carried out with monochromatic beams, while the spectral width of pulse width effectively measures approximately 10 nm. However, we expect this to have a negligible impact on the simulated coherent length, because a larger wavelength will result in a longer coherent interaction length, while the opposite will be observed for the shorter wavelength of the pulse. Therefore, these effects are expected to cancel and only the contribution of the central wavelength will remain. In conclusion, this result demonstrates that the simulation is able not only to reproduce the large-scale oscillations originating from the coherent interaction length, but also various other features.

As explained above for the experimental data it was considered that the thinnest, homogeneous thickness of the LN wedge is in the range of 50 to 100 nm based on mechanical profilometer scans. Therefore, to all experimental datasets the origin (0 nm thickness) was shifted by a constant thickness of 100 nm from the point of first significant SHG intensity, and this was accounted for in a respective thickness error bar ($\pm 100$ nm). This shift is also seen in Fig. 5, where the experimental data are plotted from 100 nm onward only. This previously artificial offset yields a very good agreement with the numerical data, adding another independent argument to the plausibility of the minimum wedge thickness in the range of 50 nm to 100 nm, as well as the plausibility of the confidence interval.

In order to theoretically evaluate the coherent interaction length for all 3 cases, we have performed calculations as the one displayed in Fig. 5. To save time, thickness increments of 100 nm (Cases I and II) and 200 nm (Case III) were used, however, using the same FEM grid for all calculations. This allows to analyze a much larger parameter space, as well as to test the influence of various different parameters as discussed in the Supplemental Information section III. Increments of 100 nm or 200 nm were deemed enough, as they allow to describe the coherent interaction length with a similar accuracy as in the experiment, as shown by the error bars in Fig. 6.

To compare the experimental and numerical data, we have determined the minima for all datasets. This data is plotted in Fig. 6 as a function of fundamental wavelength. Additionally, the graphs contain the predicted analytical behavior for plane wave incidence as describe by Eq. 2 (continuous blue line). First, we immediately notice that the analytical result always yields a larger coherent interaction length as compared to the experiment and simulation. As discussed above, this is an effect of using focused beams, while Eq. 2 is valid for plane waves, only. Indeed, the numerical data also display a shorter effective coherent length when compared to the analytical data, and fit well to our experimental results. In particular, the simulation is even able to reproduce the effect of a decreasing $l_{c,eff}$ with increasing $NA$, as displayed in Fig. 6(a) for Case I. While the confidence intervals for the two $NA$ datasets in Fig. 6 do overlap, the experimentally deduced $l_{c,eff}$ for $NA$ = 0.8 is always larger as compared to simulation. Note, though that the effective coherent interaction length for $NA$ = 0.8 is only slightly smaller as compared to the $NA$ = 0.45 case, for all wavelengths.

To exclude a systematic error, e.g. having chosen a wrong offset or different locations, we have repeated the Case I experiments for both NAs while rotating the sample by 90°. Then the detected SHG intensity is governed by $P_{x}=-d_{21}E_{x}^{2}$, however, leading to the same coherent interaction length and process as in Case I. As shown in the Supplemental Information section II, the determined $l_{c,eff}$ for both NAs and both orientations are the same, despite being independent measurements. Therefore, the slight deviations encountered here might stem from a sub-optimal illumination of the entry pupil of the NA = 0.8 objective. Please note, that our simulation always considers an ideal planar wave back illumination of the objective. However, a lower illumination factor yields an effectively reduced numerical aperture, as shown in the Supplement. Hence, the observed effective coherent interaction length will be affected.

Nevertheless, the trends and slopes are very well reproduced. The excellent agreement between the behavior predicted by numerical calculation and the experiment confirms that indeed the co-propagating signal is detected in back-reflection geometry. It further shows that the numerical simulation is able to accurately reproduce many of the effects, and presents an important tool to improve the understanding of NLO microscopy experiments.

As mentioned above, the short-period oscillations present in most datasets in Figs. 4 and 5, are due to thin film resonance effects of the fundamental and SHG beams. A resonant thickness $d$ or multiples $m$ of it, can be described by the well known relation

where $\lambda_{\omega}$ is the wavelength at the fundamental or SHG frequency and $n_{\omega}$ its respective refractive index. To investigate the origin of these oscillations, we have performed a FFT analysis of the datasets for Case I [Fig. 4(a)], which shows the strongest resonances, and the respective high resolution numerical data (Fig. 5). The results for the 850 nm fundamental wavelength are displayed in Fig. 7(a) and (b), respectively. The FFT of both datasets, numerical and experimental, show two prominent peaks at spatial frequencies of approximately 5.5 $\mu$m^-1 and 11 $\mu$m^-1 belonging to periods of about 180 nm and 90 nm, respectively. Both fit exemplary well to the expected resonance periods for the fundamental and SHG light, respectively.

Figure 7(c) shows the extracted thin film resonance periods as a function of wavelength, together with the analytical calculation based on Eq. 6. In contrast to the analysis of the effective coherent length, here the numerical aperture plays no prominent role for the thin film resonances period. Moreover, the numerical and analytical calculations align well and yield similar results. This is expected as the oscillations (a) are aligned parallel to the interface planes, as seen in the focal field plots in the supplement, and (b) the refractive indices for both, the analytical and numerical calculations, are based on the same refractive index [31]. The experimental data shows the same slope and behavior as do the numerical and analytical data, however, is systematically shifted to slightly smaller periods for both, the SHG and fundamental wavelength. A potential explanation is a slightly different refractive index of the used 5% MgO-doped LN as compared to literature values. Please note, a difference in the magnitude of the refractive index will not influence the calculated coherent length, as this is dependent on dispersion, i.e. differences between the refractive index at the fundamental and SHG wavelength, which are similar even when the absolute refractive index is shifted due to variations in ambient conditions, doping, defects or fabrication [39, 40]. Nevertheless, the excellent agreement of the FFT data (a),(b) and the slopes in (c) demonstrate the capability of the simulation to also describe resonance effects.

In this work we have investigated the impact of confinement and reflective interfaces on the SHG signal in standard back reflection microscopy. In such a setup it is commonly assumed that the detected SHG intensity exclusively reflects the counter-propagating phase-matched SHG signal, generally stemming from a tiny effective interaction volume only, due to the short coherent interaction length $l_{c}$. In addition, that phase mismatch $\Delta k$, hence, must always be positive independent of wavelength, which has severe implications on the expected SHG signatures as shown by Kaneshiro et al. [32].

However, whenever reflective interfaces of a few percent reflectivity and/or confining (nano)structures are involved, which thin films are per se, also the co-propagating phase-matched signal is reflect back into the detector, and superposes to the counter-propagating signal. Even more dramatically, $l_{c}$ for the co-propagating signal reaches significantly larger values as compared to the counter-propagating case. Additionally for confined structures, significant enhancement of the SHG response can be observed due to thin-film interference effects of either the fundamental or SHG beam.

As our experiments and numerical simulations clearly show, these confinement effects are not just limited to the film thickness on the order of the wavelength, as it is often assumed for linear microscopy, but play a major role even up to macroscopic thicknesses of a few ten micrometers. Therefore, from the viewpoint of quantitative SHG microscopy, a crystal may only be considered bulk whenever its thickness measures at least several times the axial resolution, or when spurious reflections are significantly suppressed, e.g. by index matching oil. Our results are particularly relevant when obtaining quantitative SHG values, not only of samples with varying thicknesses, but also when tuning over a broader wavelength range, solely due to dispersive and polarization effects. This can be of particular relevance when examining for instance ferroelectric thin films with thicknesses ranging from a few nm up to several $\mu$m, where SHG microscopy and polarimetry is widely used to analyze local bulk and domain wall symmetries [10, 41, 26, 42].

In this work, we have considered the SHG effect as an archetypal nonlinear optical (NLO) effect. While phase-matching length and wavelengths will depend on the specific NLO effect, it can be expected that similar results will influence the NLO response in other processes as well, e.g. CARS or THG. To judge on the specific influence and magnitude of confinement and reflections in those cases, further theoretical and experimental studies are needed.

In conclusion, our work demonstrates that reflections can play a crucial part in NLO experiments and should not be neglected. This improved understanding of NLO microscopy, hence, will lead to a much more detailed experimental interpretation and - as already suggested for SHG on TFLN [11, 33] - potential super-resolution applications of NLO microscopy.

Additional details on the experimental setup and the sample preparation, additional experiments and a numerical discussion on the influence of various parameters can be found in the Supplemental file.

The data that support the findings of this study are available from the corresponding author upon reasonable request.