Thin-film materials from $\mu $m thickness down to single-atomic-layered 2D materials play a central role in many novel electronic and optical applications. Coherent, nonlinear optical (NLO) $\mu $-spectroscopy offers insight into the local thickness, stacking order, symmetry, or electronic and vibrational properties. Thin films and 2D materials are usually supported on multi-layered substrates leading to (multi-)reflections, interference, or phase jumps at interfaces during $\mu $-spectroscopy, which all can make the interpretation of experiments particularly challenging. The disentanglement of the influence parameters can be achieved via rigorous theoretical analysis. In this work, we compare two self-developed modeling approaches, a semi-analytical and a fully vectorial model, to experiments carried out in thin-film geometry for two archetypal NLO processes, second-harmonic and third-harmonic generation. In particular, we demonstrate that thin-film interference and phase matching do heavily influence the signal strength. Furthermore, we work out key differences between three and four photon processes, such as the role of the Gouy-phase shift and the focal position. Last, we can show that a relatively simple semi-analytical model, despite its limitations, is able to accurately describe experiments at a significantly lower computational cost as compared to a full vectorial modeling. This study lays the groundwork for performing quantitative NLO $\mu $-spectroscopy on thin films and 2D materials, as it identifies and quantifies the impact of the corresponding sample and setup parameters on the NLO signal, in order to distinguish them from genuine material properties.

## I. INTRODUCTION

Thin films and nanosized materials with dimensions ranging from a few $\mu $m thicknesses down to atomically thin 2D materials offer various compelling and novel properties for applications in electronics, optics, and many other fields.^{1–6} A key task in material development as well as device fabrication is characterization and quality control.^{7,8}

In this context, parametric nonlinear optical (NLO) interactions allow for the nondestructive and sensitive analysis of many relevant material properties. NLO processes scale, as the name says, nonlinearly as a function of the incident field strength,^{9,10} for example, quadratic or cubic, and are generally very weak. However, they may become very significant, whenever strongly confined light fields, such as in focused and ultra-short pulsed laser beams, are involved. Parametric NLO optical interactions are characterized by the fact that the total energy of in- and outgoing photons match. Here, simple examples are second- or third-harmonic generation (SHG, THG), where two or three photons at a fundamental frequency are annihilated and a single photon at twice or triple the fundamental frequency is created, respectively.^{11–14} More complex processes involving photons at different energies are possible as well, such as sum or difference frequency generation, or four wave mixing.^{15–24} A key feature for all these parametric processes is that they do not require an energy transition, e.g., an electronic or vibrational excitation, to be involved, but rather depend on the symmetry and directionality of the material, only. This makes NLO processes a very elegant tool for probing the symmetry, structure, and directionality of any material with a high sensitivity.^{25}

For example, SHG microscopy is able to visualize changes in symmetry,^{12,16,25–27} such as at ferroelectric domain walls or twin boundaries in transition metal dichalcogenides (TMDs),^{28–32} analyze the direction and symmetry of stacked 2D materials,^{33,34} quantify phase transitions,^{35,36} or unravel localized changes in symmetry, such as non-Ising-type ferroelectric domain walls.^{26,27} The nonlinear scaling of the intensity and the coherent build-up results in a high sensitivity for small changes, e.g., in material thickness.^{11,12,26,27} In the case of 2D materials, THG, for example, is able to determine the number of stacked layers discriminating between a single layer up to two-digits of layers.^{37}

While NLO methods are independent of existing transitions in a material, they can be significantly enhanced in the presence of material resonances, such as electronic or vibrational states. Therefore, measuring the intensity of the NLO process as a function of fundamental wavelength allows for decent spectroscopy. For example, NLO spectroscopy in semiconductors allows us to sensitively probe weak transitions, such as complex excitonic states.^{35,36,38} One of the most widely used methods for parametric NLO spectroscopy is coherent anti-Stokes Raman spectroscopy (CARS), which addresses the challenge of weak intensity of spontaneous Raman scattering by using a resonant four wave mixing process instead. Recently, this method was successfully translated to crystalline materials as well.^{17,18}

The large intensity enhancement of a parametric NLO process compared to their non-coherent counterparts can be explained by their coherent nature, i.e., the fixed phase relation between in- and outgoing light fields. However, as the refractive index at the fundamental and generated frequency, in general, is not the same due to material dispersion, this results in a walk-off between the fundamental and generated waves, and in destructive or constructive interference depending on (1) the material thickness, (2) the focusing conditions, and (3) the refractive indices at the involved wavelengths. The coherent interaction length $lc$, the length up to which the interference is net positive, can range from well below 100 nm for counterpropagating waves up to infinity, for example, for birefringent materials at certain wavelengths or angles. It should be noted for spectroscopic applications, that depending on the analyzed wavelengths, largely different coherent interaction lengths are present within the same spectrum. This can make quantitative analysis challenging.

For thin materials, the situation is particularly complex, because many effects may influence the coherent signal generation in NLO microscopy, including reflections, thin-film interference, phase matching, and dispersion, the nonlinear scaling with interaction length,^{12} as well as the focusing strength.^{39,40} The latter includes the focus placement and, in particular, the *so-called* Gouy-phase, which introduces a phase jump between signals generated above and below the focal plane. Therefore, to accurately interpret quantitative results in NLO $\mu $-spectroscopy on thin films, these effects need to be carefully taken into account, such as by rigorous modeling. While many 2D materials have thicknesses much shorter compared to any co- or counterpropagating coherent interaction lengths, phase(-mis)matching still can play an important role in quantitative analysis, if (spaced) multi-layers are considered or if 2DM are placed on reflective substrates or even nonlinear crystals, such as in recent examples of TMDs on nonlinear, ferroelectric crystals.^{41–43} In previous work,^{11} we have developed a model for SHG microscopy based on a full vectorial description of the focusing process and performed wavelength dependent experiments on wedge-shaped samples.^{11} This allowed us to analyze the influences of resonances and phase-matching in nano-confined media.^{11}

In the present work, first, we have extended this full vectorial modeling and the experiments to THG. Note that due to its nature as a third-order NLO process, THG is present in all materials. Hence, third-order processes may be applied per se to study any material. However, due to this universality, any material above or below the sample of interest, such as air or substrate, must be taken into account as well, since they potentially contribute to a coherent interference with the signal of the sample, resulting in complex responses. This also applies for other third-order processes such as CARS. In this regard, THG is an ideal test-bed to gain fruitful insights into other more complex third-order processes like CARS.

Second, we have applied a recently developed semi-analytical approach and compare this to both the full vectorial numerical simulations and the experiments. While the full vectorial simulation is by definition more accurate because it encompasses various complex effects like focus-induced polarizations^{39,40} or focal distortions when deeply focusing into media,^{44} it comes at the cost of considerable computing resources. In this context, while the semi-analytical approach, can by definition not describe certain effects such as focus-induced polarization changes,^{39} it still allows us to gain broad insights into many experiments at significantly reduced computational costs. Therefore, the comparison will allow us to gain an understanding into the potentials of each modeling approach in SHG, as well as THG.

This paper is structured as follows: in Sec. II, the experimental details, the semi-analytical model, and the full vectorial numerical model will be briefly discussed. Section III will discuss the experimental results, while in Sec. IV, the theoretical models are compared to the experimental findings.

## II. METHODOLOGY

### A. Experiment

The experiments were performed on a Zeiss LSM980 microscope using a femtosecond, tunable laser source (Spectra Physics Insight X3, 690-1300 nm, $<100$ fs pulse duration, up to 5 W). The pump light was focused via an objective lens with a numerical aperture of 0.45 onto the sample, while the generated light was collected in backreflection via the same objective. The SHG or THG light was separated via wavelength appropriate filters and detected with a photomultiplier tube. Imaging is realized via beam scanning through a galvo-mirror system integrated into the microscope. Whenever the scan area is larger than the area determined by the field of view of the objective lens, the sample can be repositioned via a scanning stage. This allows for imaging larger scan areas via accurate stitching.

The sample was fabricated from a piece of z-cut 5$%$ MgO-doped lithium niobate (LNO) (X $\xd7$ Y $\xd7$ Z: 5 $\xd7$ 6 $\xd7$ 0.2 mm$3$), which was bonded to a (100) silicon wafer of 500 $\mu $m thickness. This sample was completely immersed in epoxy and polished at a shallow angle until a LNO wedge is formed. A photograph of the investigated sample is displayed in Fig. 1(a); Subfigures 1(b) and 1(c) show a principle sketch viewed from the top and as a cross section, respectively. Please note that a second piece of silicon was bonded to the right hand side of the LNO crystal to serve as a height reference during polishing. As the silicon substrate below, the LNO has a larger extent in y-direction compared to the LNO, a wedge of epoxy on the silicon substrate is also formed next to the LNO crystal. The LNO wedge angle was measured with a laser scanning microscope and determined to be $\alpha =(0.45\xb10.01)\xb0$. This wedge angle is considered small enough that the lateral displacements of the beam, which could lead to a less of contrast in the interference effects, are much less than 1% of the film thickness and the film interfaces can be considered parallel in the modeling. The surface roughness of the polished interface was determined via laser scanning microscopy to be better than $<25nmRMS$, which is less than a 16th of the involved wavelength and should not influence the measurements, as it was also shown in our previous work.^{11} The lower interface of the LNO should still retain its factory polish with a roughness better than $<1nmRMS$, because it was only glued to the silicon wafer. More details on sample fabrication can be found in our previous publication.^{11}

### B. Paraxial modeling

The ansatz to describe nonlinear optical interaction in focused beams is based on the paraxial solution as presented in Boyd’s textbook^{9} and was adapted and discussed by us for thin films and general NLO interactions in another publication.^{45} Here, we will recall the most important steps. In that framework, an incoming Gaussian fundamental beam is assumed and casted into a convenient form as shown by Eq. (1),

Here, $a$ is a complex amplitude, $b=2zR$ is the confocal parameter, which corresponds to twice the Rayleigh range, $w0$ is the beam waist radius, and $r$ and $z$ are the cylindrical coordinates, where the beam propagation is along the $z$-direction and no angular dependence is present.

The well-known coupled wave equation for nonlinear phenomena, including the slowly varying envelope approximation, will be used here in cylindrical coordinates [cf. Eq. (2)],

Here,

is commonly identified as the phase mismatch, where $q$ represents the order of the harmonic process or the number of involved fundamental photons, respectively, e.g., $q=2$ for SHG and $q=3$ for THG. Note that the refractive index $n=n(\lambda )$ is a quantity that depends on the wavelength of the light. Thus, as in almost every material, the refractive index is not constant, and, therefore, a certain phase mismatch is always present. Usually, one defines the distance up to which the waves are partially positively interfering with the harmonic signal growing, as the coherent interaction length [cf. Eq. (4)]. For plane waves, $\Delta k$ is the only contribution to the signal oscillations,

This changes, however, for focused beams with more general and realistic descriptions. The inclusion of the beam character, e.g., with a Gaussian beam ansatz, introduces more contributions to the walk-off of fundamental and harmonic waves, which are mostly caused by the additional phase evolution in focused beams, the Gouy-phase. Focused light will experience a phase change of $\pi $ when crossing the focal plane as $\Phi Gouy\u221darctan\u2061(zzR)$. For parametric processes, this has important consequences. Due to the exponential nature of the phase, parametric light generated at different positions will acquire additional different phase contributions, resulting in a position-dependent additional phase-mismatch. For many scenarios with phase mismatch, this prominently results in zero SHG and THG intensity from focused beams in bulk (infinitely extending) materials and signals stemming from interfaces^{9,45} as well as vanishing THG intensity, even if there is complete phase-matching from the refractive index point of view.

If we now consider SHG and THG processes, where three or four light fields interact, respectively, i.e., two or three pump fields, and one generated “out” field, one or two Gouy-phase shifts of $exp(\u2212i\Phi Gouy)$ are added to the harmonic, resulting in a total additional shift of $exp[\u2212(q\u22121)\u22c5i\Phi Gouy]$. This may have a significant impact even for thin films, as the phase $\Phi Gouy=\u2212arctan\u2061(zzR)\u2248\u2212zzR+13(zzR)3$ can be seen as a sort of additional phase mismatch, which can be annihilated by the positive phase mismatch of anomalous dispersion.^{46} Furthermore, as the shift is non-linear, the distance between maxima and minima of the SHG/THG signals depends on the location with respect to the focal plane and oscillation order, i.e., the thickness. For normal dispersion, the coherence length will be lowest for a material placed in the focal region, due to the influence of the Gouy-phase.^{11,45}

In Ref. 45, we, therefore, suggest to incorporate the Gouy-phase evolution [cf. Eq. (5)] and, therefore, also the phase reference point $l0$, i.e., the coordinate relative to the focal plane where the interaction begins, into the coherence length computation. This delivers a more exact description as it accounts for all phase-contributions with Eq. (7),

From that, it follows that the distance for a specific phase change can be calculated via the phase difference shown in Eq. (7),

Using the linear approximation of Eq. (6) and considering a situation where entry point and the focus coincide ($l0=0$), we obtain a simple extension to the classical formula for $lc$.

Note that the approximation is only valid for $L\u226ab$. One can well observe how the coherence length decreases in scenarios with negative phase-mismatch^{11,45} and in principle—due to its nonlinear content—also depends on the focal position and the total acquired phase $\Delta \Phi $, i.e., $l(2\pi )\u22602\u22c5l(\pi )$. The latter relation is only fulfilled for plane waves. As a result, instead of calculating $lc=l(\pi )$, we have calculated the oscillation period $l(2\pi )$, as the distance between two minima. As the distance between minima is easier to determine in the data, due to the small amplitude of resonant oscillations present at the point of minima, compared to calculating the distance from a minimum to a maximum, i.e., $lc$, this process introduces the least uncertainty.

In order to solve the differential equation, we use a trial solution,^{9} which is similar to a Gaussian beam, but shows a z-dependent amplitude [cf. Eq. (9)],

Inserting Eq. (9) into the differential equation (2), one obtains in good approximation^{9} a solution for the amplitude in the form of an ordinary differential equation,

Here, $clayer=iq\omega 2nc\chi layer(q)$ depends on the material, specifically its refractive index $n$, and the nonlinear susceptibility $\chi layer(q)$ at order $q$ of the respective layer, whereas $z$, $z0$ are the entry and exit points of the fundamental beam of that layer, so that $d=z\u2212z0$ is the thickness of the layer.

The integral has to be solved numerically for finite boundaries. Furthermore, especially when the sample sits on a strongly reflecting substrate like silicon, one cannot neglect reflection and transmission. In principle, in our case, the interface of LN and air does have a high transmission of $T\u22a5(0.85\mu m)\u223c85%$, and it is sufficient to take into account only a small number of reflections. It can be shown^{45} that one has to compute two integrals, one for forward propagating waves and one for backward propagating waves, which interfere after reflection, where multi-reflections are accounted for by multiplying the integrand with a geometric series to the $q$th power. In an analog way, the multi-reflections of the harmonic light can be implemented with a second geometric series as shown in Eq. (11),

Here, $kx$ represents the modulus of the fundamental or the harmonic k-vector, $nx$ the corresponding refractive index, and $rx$ summarizes the Fresnel-coefficients for the reflections at the interfaces for the specific situation, meaning if it is co- or counterpropagating and fundamental or harmonic, respectively. These series result in oscillations much faster than that of the phase mismatch; they are directly correlated to the fundamental and harmonic wavelength. Therefore, these fast oscillations are easily damped by incoherence of the source, especially due to the spectral width of the pulsed laser, but also by the effects of the optical components such as lenses. Furthermore, with respect to the coherence length, oscillations originating from thin-film interferences are more of cosmetic nature and can be analyzed and filtered out via a Fourier analysis, as we will see in Sec. III.

It should be noted here, that for SHG the calculation is in that sense simpler, as only the lithium niobate slab contributes to the second-harmonic signal due to the vanishing of the $\chi (2)$ in all other materials (air, silicon) at hand. Therefore the parameter $c$ of Eq. (10) is not important and we end up with an according integration of the first incoming and reflected beam, contributing to a signal which is further modulated by the reflections of the harmonics.

For THG, however, one needs to take also into account the signal produced in air and finally sum over both, so that a more complex interference pattern is possible. Here, the exact relation (and phase) of the nonlinear optical susceptibility of air and lithium niobate gets very important^{45} as it defines the strength and phase of the respective contribution, as will be shown in Sec. III.

In general, solving these one-dimensional integrals does not need much computing resources, allowing for fast calculations. For example, simulating a signal evolution for a thin wedged crystal means solving integrals of the type of Eq. (10) for a number of different thicknesses up to a maximum thickness of, e.g., $t\u223c10$ $\mu $m. Taking $\Delta t=5$ nm as the step size and, therefore, $N1=2000$, while doing, e.g., $N2=300$ different total thicknesses which correspond to $\Delta x\u224833$ nm and $N\lambda =6$ different wavelengths, the computation takes approximately 1 min. Should it be necessary to simulate the air-layer, computation time rises significantly, as the field, depending on the parametric process, drops relatively slowly; to acquire good results $z\u226bb$, which may be the case for $100$ $\mu $m. For a $5$ nm step size, this would result in $N1=2000\u2192N1=22000$, i.e., approximately tenfold computation time of $\u2a86$ 10 min. However, as the fastest oscillation is the harmonic $\u223c400$ nm, and the important phase-change of $\Delta k$ is even on a larger scale, a step size of $10\u221220$ nm is viable. In particular, the air layer where only the $arctan(z)$ contributes can be simulated with larger step sizes of $50$ nm without changing the result and brings computation time down to two minutes. In sum, the average computation time is quite short and a full description of thin films spanning ranges of several $\mu $m can be done in a minute or for single line-scans even below that. The calculation of the coherent interaction length alone is even faster, which can be acquired independently from a full signal simulation by solving the transcendental equation (7) with a root finding algorithm, which takes seconds for a set $N\lambda =6$ different wavelengths.

### C. Numerical simulation

To complement both our experiments and the semi-analytical paraxial approach, numerical calculations in a full vectorial approach of the focusing and detection process is performed here. The calculation is executed in terms of the angular spectrum representation, e.g., as found in Novotny and Hecht,^{10} which was modified by Sandkuijl *et al.*^{47} to include reflections from an arbitrary amount of interfaces. The original code used here is based on the works by Sandkuijl *et al.*^{47} and was modified for our purpose.^{11,12,48} The calculation is structured as follows:

The input electric field distribution is calculated in the focal plane. For this work, we always assume a fully filled pupil.

In the second step, at each grid point, the local NLO dipole moment is calculated for SHG or THG, respectively. In a simplification, only a single tensor element ($d22$ or $d222$, respectively) is assumed to be non-zero. As shown in our previous work,

^{11,12}such an assumption is justified, because we assume the fundamental light to be polarized parallel to the*y*axis of the LNO crystal, only, and contributions from off-diagonal elements or from focus induced polarizations are negligible—especially for a low numerical aperture of 0.45 as used in the experiment.In the third step, the NLO polarization at each grid point acts as an emitting dipole, which emits light at the SHG or THG frequency, which then propagates into the far field. Here, again, reflections and interference from an arbitrary number of interfaces is accounted for. The light may interfere in the far field according to the respective phase between the dipoles, which accounts for NLO phase matching.

In the last step, the respective light fields at SHG or THG at the back of the focusing objective is calculated via the reverse process of step (1). The light at the back focal plane is integrated to describe the pointwise detection. The refractive indices for the respective materials are calculated based on the Sellmeier equations from Refs. 49 and 50, respectively. Further details on the simulation, its accuracy and convergence tests, as well as a theoretical discussion of various influence factors, such as small deviations in refractive index or wavelength, can be found in the previous work.

^{11}

In the work here, we have chosen a simulation volume in the lateral direction of $2\xd72\mu m2$ and up to 15 $\mu $m in axial direction to account for the coherent interaction length of the SHG signal at the chosen wavelength. The grid spacing was 16 nm in all spatial dimensions. As shown in previous work, this will yield accurate, convergent results.^{11} In this work, the SHG and THG responses for scans of the wedge-shaped sample, i.e., for different total layer thicknesses, were calculated using increments of 100 nm of total layer thickness. Such a calculation typically takes several minutes of computational time for each increment. Hence, simulating a full scan from a 0 to 15 $\mu $m thickness at a single wavelength may take more than $>10$ h in total calculation time. While the simulation code likely can be further optimized, we can see from this simple example that the semi-analytical model offers significant speed advantages for cases, where the physical accuracy of a full vectorial description is not required. It should be noted that due to the computational limitations in our numerical simulations of the THG process, only a $\chi (3)$ in the LNO layer was accounted for, but not in the air or substrate, in contrast to the semi-analytical model. Accounting for the $\chi (3)$ above the focal plane requires a significantly increased computational and memory demand, as the calculation is performed on a 3D grid.

## III. RESULTS

To analyze the thickness-dependent response of the SHG and THG signal, we imaged along the wedge varying both the sample thickness and the fundamental wavelength from 1250 to 1300 nm of excitation $Z(X,\u22c5)Z\xaf$ . In order to separate the SHG and THG signal from the fundamental, we used wavelength appropriate filters (bandpass filters: 420–480 and 610–670 nm). Figure 2 displays large scale SHG [Fig. 2(a)] and THG [Fig. 2(b)] overview images taken at a 1280 nm fundamental wavelength. In both images [Figs. 2(a) and 2(b)], the wedge thickness gradually increases from left to right, as seen in the sketch at the bottom of each figure. In both images, we can clearly see the characteristic long-period oscillations stemming from the phase matching process, producing wavefronts aligned orthogonally to the thickness increase (i.e., the wedge slope). As we have shown previously, they belong to the co-propagating phase matching process rather than the counterpropagating process and are visible through reflections at the back interface.^{11}

When analyzing these images in detail, several observations can be made. First, we notice that these oscillations are visible over the entire length of the wedge. As the angle (0.45$\xb0$) and approximate length (6 mm) of the wedge is known, we can calculate the thickness of the wedge at the right edge, to measure approximately 50 $\mu $m, where reflections of co-propagating phase-matched light obviously influence the overall SHG as well as THG intensity. This value is more than 2 times larger as compared to the axial optical resolution $\Delta z$ in LN, which is more than 20 $\mu $m for the given wavelength and numerical aperture.

Second, 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 30–40 $\mu $m for both the SHG and THG signal. Normalizing with respect to the inclination angle (0.45$\xb0$) yields thickness differences between neighboring intensity maxima in the range of 200–300 nm. These oscillations are associated with thin-film interference effects of the fundamental and the respective SHG or THG signal as will be discussed below. As it is apparent from these images here, the overall SHG or THG signal is being significantly altered by the thin-film resonance over thickness variations that measure less than 100 nm. Moreover, these oscillations are visible up to thicknesses of at least several ten $\mu $m. This means that, in particular, for quantitative measurements of SHG or THG intensities, reflections, and the co-propagating process need to be carefully evaluated even up to a 50 $\mu $m crystal thickness and likely beyond. Furthermore, the interface reflectivity, thin-film interference, and the coherent interaction length all are governed by the dispersion of the refractive index. Hence, the responses will also be affected when inspecting crystals of a constant thickness while varying the incident wavelength. These findings are excellently supported through our simulations and for SHG have been documented in our previous work in more detail.^{11}

Third, we notice that in the image for SHG [Fig. 2(a)] only the part with the LN wedge produces a significant signal, whereas in the image based on the THG signal [Fig. 2(b)] every part of the image is showing a response. This stems from the fact that every material possesses a third-order non-linear optical susceptibility ($\chi (3)$), including silicon as well as the epoxy resin at the top and bottom of the image. As a result, we observe THG signal from the surrounding epoxy wedge and the bare silicon substrate to the left. For the wedge formed out of epoxy resin, also characteristic oscillations are visible at the bottom of Fig. 2(b), but with a much larger coherent interaction length due to the lower dispersion of the refractive index in the polymer.

To analyze the characteristic oscillation length $l2\pi $ for SHG and THG as a function of wavelength, we performed detailed scans at the thinner left end of each wedge. Scanning at this location also avoids the larger areas of largely uneven glue distribution at the bottom interface of the sample, which can be seen at the center of the sample in both the SHG and THG images. To avoid any ambiguity or error through stitching or rendering, we only used single frames here that at least contained the first full oscillation, as shown by the insets in Figs. 2(a) and 2(b). Here, the focus was placed approximately at a depth of the center of the film at the thin edge of the wedge. As discussed in our previous work,^{11} the first oscillation period is least subjective to variations in focus positioning. If the images are inspected closely, a granular noise can be observed in both the SHG and THG images. As discussed above, the top and bottom interface roughness, respectively, are considered optically flat. We believe this is an effect of small scale variations in the thickness of the glue of much less than 100 nm. As shown in the previous work by simulations,^{11} such small variations only influence the absolute level of the intensity but not the phase matching and interference effects. To reduce the effects of local variations in reflectivity or laser noise, lateral binning was applied in order to obtain averaged intensity profiles over the (yellow) shaded areas in the insets in Fig. 2.

For both cases, the fundamental wavelength was swept from 1250 to 1300 nm in increments of 10 nm, via the microscope control software while the sample was kept at a constant location. The location dependent data were converted to thickness via the known wedge angle (0.45$\xb0$) and is displayed in Figs. 3(a) and 3(b) as a function of thickness for SHG and THG, respectively. What becomes crucial is to exactly define the origin for both cases, i.e., the exact position of zero thickness for each SHG/THG profile as this will influence the measured coherent interaction length. Here, we calculated the origin for both SHG and THG, using the SHG data, by first identifying the location of the first significant SHG intensity ($>0.1$% of the maximum intensity as shown in our previous work^{11}). Based on profilometer scans of the wedge edge, we determine the initial step height of the LN wedge, i.e., the minimum thickness of homogeneous LN still attached to the wedge, to be in the range of approximately 50–100 nm—depending on the exact location. Therefore, we added to each position of significant first SHG intensity an offset of 100 nm thickness and conservatively assume a confidence interval for the converted thickness of $\xb1100$ nm. As before, the assumption of an initial step of 100 nm also fits well with observations from the numerical simulations. It should be noted that a confidence interval of $\xb1100$ nm corresponds to less than 2% (5%) of the total value of $l(2\pi )$ for SHG (THG), which is larger than 7.5 $\mu $m (2 $\mu $m) for SHG (THG). Therefore, any further uncertainty in the measured data, e.g., imposed by the uncertainty in the measured wedge angle of $\xb10.01$, is negligible for wedge thicknesses of $<15$ $\mu $m, which are at most considered for analysis here. This artificial origin, i.e., point of zero thickness, was then applied to the THG data as well, as both measurements were performed exactly at the same location while only switching the filters in the detector path.

For optimal comparability, all datasets are normalized in their magnitude by dividing each data set by its maximum value and displaying it with an artificial offset. Below the experimental data, the data obtained by the semi-analytical paraxial model and the numerical simulation are displayed. Both datasets were normalized and plotted the same way as the experimental data and are displayed in Figs. 3(c)–3(f), respectively.

Overall, for SHG or THG, as well as for experiment and theory, a distinct oscillation period is observed, which increases with larger wavelength as expected from the decreasing dispersion. In detail, for the experimental values for SHG, one can see an increase in the position of the first minimum from $l2\pi \u22487725$ nm to $l2\pi \u22488582$ nm [Fig. 3(a)], whereas the expected values from plane wave theory [using Eq. (4)] are $l2\pi \u22489460$ nm to $l2\pi \u224810540$ nm. For THG, a shorter period is observed [Fig. 3(b)] with the position of the first minimum shifting from $l2\pi \u22481930$ nm to $l2\pi \u22482315$ nm. The expected values from plane wave theory [using Eq. (4)] are $l2\pi \u22482270$ nm to $l2\pi \u22482410$ nm. The values of plane wave description are larger compared to the experimental data as well as the adapted theory, because it does not account for the effects of numerical aperture, focus positioning, and the Gouy-phase.

A very remarkable thing that we observe in the case of THG [Fig. 3(b)] is the presence of a non-zero signal at the origin, i.e., the point of zero thickness of the LN wedge. The signal starts with a value comparable to the value in the next maximum and then decreases until it reaches a minimum. From here on, the usual oscillations of the signal similar to the SHG experiment are observed.

This can be readily explained in terms of the universality of the $\chi (3)$-process, which is present for any material including air above the LNO or substrate.^{45} Due to the coherence of the process, this light can interfere with light generated within the LN wedge being analyzed. For example, if the focal plane lies on the top interface, and, therefore, half the beam interacts with the air layer, the total amount of additionally accumulated phase can be analytically calculated^{45} via computing the amplitude (and its phase) for the phase matched signal in the air to be $\varphi =\pi 2$. We see that the wavefront has not propagated $(q\u22121)\pi 2$, as the total phase comprises of different wavelets originating from different positions, each wavelet contributes different phase, making the total phase smaller than expected. In this case, as we see in Fig. 3(b), the signal first decreases down to a minimum, as the additional phase gained in LNO is negative due to the negative phase-mismatch and, therefore, raises the difference continuously to $\pi $. Subsequently, the signal increases again.

Figures 3(c)–3(f) show the simulated results based on the semi-analytical paraxial model and the numerical simulation, respectively. For the case of SHG [cf. Fig. 3(c)], we see for both models a very good agreement with the experiment. The $l2\pi $ changes similar to the experiment as a function of wavelength from $l2\pi \u22487825$ nm at 1250 nm fundamental wavelength to $l2\pi \u22488698$ nm at 1300 nm. We see a less than 5% difference in the values of the coherence length, with the simulated values being larger. Similarly, in the case of the numerical simulation in Fig. 3(e), we see an increase in the position of the first minimum from $l2\pi \u22487897$ nm to $l2\pi \u22488752$ nm at 1300 nm wavelength. The difference between the values of the coherence length for the experiment and simulation is less than 2%, with the value from simulation being slightly larger than in the experiment, as well. Furthermore, due to the small step size, the paraxial simulation reproduces the short wavelength oscillations due to thin-film interference very well. In contrast, the 100 nm step width in the numerical simulation is not enough to reproduce the oscillations due to thin-film interference.

For the THG case, similar observations are made. Again, both simulations predict the coherent length $l2\pi $ very well, with less than 5% deviations from the experiment in the paraxial case, and less than 7% for the numerical model. Here, the deviations are somewhat larger, which may be explained due to larger uncertainties in accurately determining the minimum, as the short-range thin-film interference effects partly overlap the minima in both experiment and simulations. If one looks closely at the data, there is one key deviation between the paraxial model and the numerical model. For the paraxial model, we assumed—as it is also the case in the experiment—a non-zero nonlinear susceptibility $\chi (3)$ for the air above the sample. In terms of magnitude, we have chosen a value of 1/2 of the susceptibility of the LNO in the semi-analytical calculation. As discussed above, this results in the observed pattern, where the intensity first drops from zero thickness to a minimum and the approximately half period shift of the oscillations. In contrast, in the numerical simulation due to computational limitations, the nonlinear susceptibility of air was not accounted for. Therefore, these data show a very similar behavior to SHG, i.e., first showing an increase in THG signal with the coherent length oscillations following accordingly. However, as one can see, this does not affect the characteristic oscillation length $l2\pi $.

## IV. COMPARISON OF EXPERIMENT AND THEORY

To compare the experimental values with the predictions from the numerical and semi-analytical model, we have determined $l2\pi $ based on all graphs in Fig. 3. As discussed above, $lc$ is not necessarily equal to half the distance between two minima, i.e., $l\pi =1/2\u22c5l2\pi $. However, as the minima are easier to determine, here for simplicity, we have evaluated $l2\pi $ for the following analysis. To be consistent with data analysis throughout the SHG and THG datasets, respectively, we evaluated the oscillation period $l2\pi $ for the SHG data by determining the position of the first minimum, while for the THG datasets, we evaluated the difference between the first and second minimum positions, respectively.

The obtained datasets are plotted in Fig. 4 as a function of fundamental wavelength. Additionally, the graphs contain the predicted analytical behavior for the plane wave case as described by Eq. (4) (continuous black line). As expected, the plane wave description always yields a larger coherent interaction length as compared to the experiment and both simulation approaches. As discussed above, this is an effect of using focused beams as Eq. (4) is valid for plane waves, only. For the case of SHG as seen in Fig. 4(a), we see that datasets of paraxial model and numerical simulation overlap. This could be a result of using the same Sellmeier equation for the refractive indices in the calculations.^{49} This further shows the excellent agreement of the paraxial model and numerical simulations. In the case of THG as seen in Fig. 4(b), the overlap is worse as compared to the SHG results. This is due to the larger error resulting from determining the minima because of the overlapping thin-film interference patterns. Nevertheless, the trends and slopes are very well reproduced.

Moreover, the experiment and both simulation datasets also contain short-period oscillation patterns. These short-period oscillations [cf. Fig. 3] are due to thin-film interference at the fundamental, SHG, or THG wavelength, respectively. 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, SHG, or THG frequency, and $2n\omega $ is its respective refractive index. To investigate these oscillations, we performed a FFT analysis of the experimental datasets for SHG and THG. The results for the 1280 nm fundamental wavelength are displayed in Figs. 5(a) and 5(b), respectively. Both fit exemplarily well to the expected resonance periods for the fundamental and SHG/THG light, respectively.

Figures 5(c) and 5(d) show the extracted thin-film resonance periods as a function of wavelength, together with the analytical calculation based on Eq. (13). In contrast to the analysis of the coherence length, here the use of a focused beam plays no prominent role for the thin-film resonance periods, which we have also seen in previous work.^{11} It should be noted that for THG, the oscillation periods only up to 1280 nm could be analyzed. If one inspects the FFT spectrum in Fig. 5(b) carefully, it can be seen that the amplitude of the THG associated peak at approximately $12\mu \u22121$ is already weak compared to the noise level. For the fundamental wavelengths of 1290 nm and 1300 nm, respectively, the amplitudes were below the noise level and could not be identified with certainty.

The experimental data show the same slope and behavior as the numerical; however, it is systematically shifted to slightly smaller periods for fundamental and SHG/THG wavelengths. A potential explanation might be a slightly different refractive index of the used 5% MgO-doped LN as compared to literature values. Please note that a difference in the absolute $value$ of the refractive index will not affect the calculated coherent length as much, as the latter is dependent on $dispersion$, i.e., differences between the refractive index at the fundamental and SHG/THG wavelength, that are similar even when the absolute refractive index is shifted due to variations in ambient conditions, doping, defects, or fabrication.^{51,52}

## V. CONCLUSION

In this work, we demonstrate that rigorous theoretical analysis allows us to disentangle the parameters that influence the parametric NLO process in focused beams on thin-film materials. For demonstration, we used SHG and THG as two archetypal NLO processes. We compared the experimental finding to two models, a semi-analytical paraxial model and a numerical simulation in the full vectorial approach. As demonstrated before, the numerical simulation in the full vectorial approach is a strong tool for any quantitative analysis, as it encompasses a wide range of physical effects,^{39,40} but it may be limited by computing resources. In contrast, the semi-analytical paraxial model can describe most effects in the model system at considerably lower computational costs and at the same time provides a reasonable insight into the physics behind the NLO process.

In our experiment as well as the two simulation approaches, we demonstrate how thin-film interference, reflection, and phase matching strongly affect the signal strength in SHG or THG microscopy. For both processes, we find that due to the presence of reflections at the back interface, the co-propagating forward signal overpowers the counterpropagating signal, which one may naively expect in a backreflection geometry. The most surprising result is that due to the omnipresence of $\chi (3)$ for any material and the Gouy-phase shift in the focus, THG from thin films first may start to vanish with increasing thickness instead of increasing. This means, for THG experiments on thin films, that for a certain film thickness or wavelength, a much lowered signal may be observed, which has to be taken carefully into account when performing experiments.^{53}

## ACKNOWLEDGMENTS

The authors gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) through projects CRC1415 (ID: 417590517), EN 434/41-1 (TOP-ELEC), INST 269/656-1 FUGG, and FOR5044 (ID: 426703838) as well as the Würzburg-Dresden Cluster of Excellence on “Complexity and Topology in Quantum Matter”—ct.qmat (EXC 2147; ID 39085490). Also, we would like to acknowledge excellent support by the Light Microscopy Facility, a Core Facility of the CMCB Technology Platform at TU Dresden, where the SHG and THG analysis was performed.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Zeeshan H. Amber:** Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). **Kai J. Spychala:** Methodology (supporting); Software (equal); Validation (equal); Writing – original draft (supporting); Writing – review & editing (supporting). **Lukas M. Eng:** Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). **Michael Rüsing:** Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

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

## REFERENCES

*et al.*, “

*2020 Joint Conference of the IEEE International Frequency Control Symposium and International Symposium on Applications of Ferroelectrics (IFCS-ISAF)*(IEEE, 2020), pp. 1–4.

*Solid-State Mid-Infrared Laser Sources*(Springer, 2003), pp. 99–143.

*et al.*, “

*et al.*, “