Coherent nonlinear optical $\mu $-spectroscopy is a frequently used tool in modern material science as it is sensitive to many different local observables, which comprise, among others, crystal symmetry and vibrational properties. The richness in information, however, may come with challenges in data interpretation, as one has to disentangle the many different effects like multiple reflections, phase jumps at interfaces, or the influence of the Guoy-phase. In order to facilitate interpretation, the work presented here proposes an easy-to-use semi-analytical modeling *Ansatz*, which bases upon known analytical solutions using Gaussian beams. Specifically, we apply this *Ansatz* to compute nonlinear optical responses of (thin film) optical materials. We try to conserve the meaning of intuitive parameters like the Gouy-phase and the nonlinear coherent interaction length. In particular, the concept of coherence length is extended, which is a must when using focal beams. The model is subsequently applied to exemplary cases of second- and third-harmonic generation. We observe a very good agreement with experimental data, and furthermore, despite the constraints and limits of the analytical *Ansatz*, our model performs similarly well as when using more rigorous simulations. However, it outperforms the latter in terms of computational power, requiring more than three orders less computational time and less performant computer systems.

## I. INTRODUCTION

The development of potent and reliable nonlinear optical materials enables a broad field of nonlinear optical (NLO) applications,^{1} such as frequency converters,^{2} sources of squeezed light, as well as heralded single photons,^{3} to name a few examples. The materials need to fulfill the typical demands concerning optical materials, i.e., showing low absorption/scattering, and a reliable quality. These criteria are supplemented by the requirement of large optical non-linearities. Furthermore, it is very important to be able to create functional structures, for example, waveguide-networks, or to combine the material for instance with silicon to form hybrid structures.^{4,5} In typical material systems like lithium niobate (LN), the established technologies of waveguide formation and other functional structures get increasingly replaced by the utilization of thin films (tfs) and other nanosized structures.^{1,6,7}

One of the key tasks in material development and device technology is quality control^{8} and improvement. Here, parametric nonlinear optical processes allow for a non-invasive, fast, and very sensitive probing of different properties. Therefore, they are widely used as a material characterization technique in the form of second- and third-harmonic-microscopy^{9–15} as well as via Coherent Anti-Stokes Raman Scattering^{16,17} (CARS), for example.

When characterizing thin films and using (strong) focusing, the typical textbook description of the nonlinear optical interaction, however, breaks down rapidly or is at least inaccurate. One main issue is, for example, that results are often interpreted based on plane wave models, which do not include focusing. This, for example, leads to a massive overestimation of the coherent interaction length. To counter this, many numerical simulation codes are available^{18} and deliver good results but are often not able to provide results that intuitively allow for disentangling different effects, while also needing significant computing resources.^{11} This limits the ability to scan large parameter spaces. Therefore, this paper tries to clarify the possible applications and limits of the analytical textbook results. Furthermore, we aim to extend these approaches to focused beams and interpret those solutions. Here, we present a modeling toolkit, which enables accurate description for typical characterization setups of (thin film) nonlinear optical materials. At the same time, we aim at conserving and/or extending the established nomenclature and interpretation of typical parameters like the coherent interaction length, the confocal parameter, and the Gouy-phase in order to preserve a simple and intuitive understanding. This approach saves computation power and can help in the development of new methods or rapid analysis of results by allowing quick identification of key parameters in experiments.

This paper is structured into six sections. The sections do not necessarily need to be read in consecutive order as they address different interests such as looking for modeling results for specific processes, discussions of the details of the modeling, or its limits and prospects.

If the reader is interested in a brief overview of established textbook descriptions, on which our semi-analytical model is built, and a discussion on their limitations, then these can be found in Sec. II. This section presents a brief reminder on the textbook description of nonlinear optical interactions in crystals. Here, the plane wave and Gaussian beam approaches are discussed and the relevant parameters are introduced.

Section III contains a short discussion on the general relations and limits between full numerical and semi-analytical solutions in the context of nonlinear interactions in focused beams.

The improved semi-analytical approach is introduced in Sec. IV. This chapter deals with the extension and reformulation of the standard analytical textbook models in order to reliably describe nonlinear optical microscopy in thin films. Crucially, Sec. IV B discusses the role and meaning of the coherent interaction length in focused beams and its differences compared to the well-established theory that bases on plane waves.

Should the reader be mainly interested in the performance of the model vs experimental data, as well as its power for predicting nonlinear optical phenomena in focused beams, then the reader may directly jump to Sec. V. This section contains different applications of the model for bulk crystal second-harmonic generation (SHG), as well as SHG and third-harmonic generation (THG) in thin films.

Section VI closes the paper with a summarizing discussion and outlook.

## II. TEXTBOOK APPROACHES AND THEIR LIMITS

In many lectures on nonlinear optics, the model calculations are introduced on the basis of planar waves, cf. Eq. (1), with a position-dependent amplitude $aq$ as solutions for the nonlinear wave equation

Here, $q$ represents the harmonic order, i.e., $q=2$ belongs to SHG. One then makes use of the slowly varying envelope approximation and solves the approximated version of the nonlinear wave equation, so that the amplitude $a$ of the newly generated wave at a specific position, e.g., where it leaves the crystal, can be computed via Eq. (2). It yields the archetype $sinc(x)$ solution as

Here, the phase mismatch $\Delta k\u2192=\u2211i,ok\u2192i\u2212k\u2192o$, i.e., the difference between wave-vectors of the incoming ($k\u2192i$) and produced waves ($k\u2192o$), is the key parameter one needs to know in order to predict any signal.

This picture may already deliver a good description if its approximations are fulfilled, i.e., there is almost no focusing. In the textbook of Boyd,^{19} the paraxial coupled wave equations, cf. Eq. (3), are solved using a Gaussian *Ansatz* as in Eq. (4) and, therefore, also include weak to medium focusing conditions,

The new key parameters in Eq. (4) are the confocal parameter $b$, as well as the beam waist radius $w0$. This more sophisticated solution, therefore, also justifies the criterion for using plane waves, namely, that the interaction length $L$ should be much smaller than the Rayleigh range $zR=b/2$ or in other words, the Rayleigh range should be very large, which corresponds to weak focusing so that the amplitude is stable and only the exponential of the plane wave (p.w.) *Ansatz* in Eq. (2) determines the size of the amplitude $a$ in Eq. (5),

In Boyd’s textbook,^{19} the resulting integral Eq. (6) is solved for large crystals and focal positions far away from interfaces via the residue theorem

Equation (7) shows that only for positive phase-mismatch, there is a signal when the focus is positioned in the crystal far away from the surfaces. Interestingly, even for full phase-matching, the signal in the bulk vanishes. One has to note here, however, that the solution for the case $\Delta k=0$ is only valid for $q\u22653$, though. For SHG ($q=2$) and perfect phase-matching, the signal does actually not drop to zero. In the phase-matched case, it holds $ei\Delta k=1$. Therefore, one can directly integrate the remaining $1(1+i2z/b)q\u22121$ term; for $q=2$, this leads not to a monomial but a logarithmic dependence, which does not vanish for large arguments. One can also picture the problem by using the analytic expression and trying to decrease $\Delta k$ to vanish. For THG and all higher harmonics, the transition is smooth, but not for SHG, where the result is $0$ or $b\pi $, depending on the order of $\Delta k\u21920$ and $q\u21922$. One can use the direct approach of the integral and successively decrease the symmetric boundaries to zero, whereas the phase term is isolated

For arbitrary boundaries, this difference of large terms is not a very suitable derivation, though. However, one can also argue that the physical electric field is given by the real part, i.e., $2Re(E)=E+c.c$. Therefore, looking at the real and imaginary part of the electric field

The real part is well defined for arbitrary boundaries and the imaginary part vanishes inside large crystals due to point symmetry. That means that SHG does not vanish in bulk if perfectly phase-matched and focused; however, the strength of the signal decreases with larger NA as it depends linearly on the confocal parameter in *b*. To sum up, the analytical approximations in Eq. (6) enable to include focusing into the description of the harmonic generation process. It is also shown that focusing changes the value of the phase mismatch, which is, for example, also observed in the recent experimental work^{11} of the authors. To obtain more realistic predictions, however, one also needs to utilize numerical strategies.

## III. RELATIONS TO FINITE ELEMENTS APPROACH AND AVAILABLE SEMI-ANALYTICAL EXTENSIONS

Since the integrals of Eq. (6) are strictly analytically calculable only for infinitely large crystals, they have to be solved numerically for finite boundaries, using finite numerical elements; one may, however, consider such a common numerically-aided approach which uses intuitive analytical trial expressions, including approximations, as semi-analytical. In contrast, more rigorous numerical toolkits solve the problem using the complete or only partially-approximated nonlinear wave-equation directly via summation over very basic structures.

Such an approach is employed, for example, when using plane waves in the angular spectrum representation, or when propagating an arbitrary polarization distribution via convolution with the point source propagator.^{9,11,18,20,21} Given enough computing power, they have the advantage that they provide more accurate and sometimes more rich results, as they often solve the full 3D-problem without relying on symmetry. The latter enables to make very precise examinations of effects, which are not necessarily included when using specific trial solutions, like focus distortion and depolarization^{9,20–23} for tight focusing conditions, as the trial solutions are mostly tied to a specific parameter regime. In these regimes, however, the semi-analytical approach, albeit being slightly less precise, has two other advantages. First, the needed computing power is much lower, which enables fast estimations and evaluations as well as scans over a larger parameter space of different setup conditions, such as focus depth, crystal thickness or focusing strengths, for example, in order to optimize the setup. Second, apart from final results, i.e., the field distribution or output signal, single phenomenological parameters like the phase-mismatch, Gouy-phase, Rayleigh range, etc., can be associated as beam parameters, and their influence can be understood more easily. This makes these approaches not only important for educational reasons but also when extracting specific observables or improving the experimental method. It is critical, however, to keep track of the specific limits of the approach to make sensible interpretations. Such approaches do exist for many applications.

A similar approach can be found, for example, in the influential work^{24} of Kleinman and Boyd, which discusses in detail the second-harmonic interaction of Gaussian beams and works out strategies to optimize the parametric output. There, the ratio of crystal length $l$ and confocal parameter $\xi =l/b$ is the key parameter, which has to be optimized. The in-depth discussion also uses a semi-analytical Gaussian approach where discussions on limiting cases are possible and provides an intuitive understanding of the situation.

Our work here deals with related topics, as far as parametric interactions of Gaussian beams are concerned. However, the main focus lies in the computation and understanding of parametric processes for material characterization, especially thin films and layered materials, using, e.g., scanning microscopy. Here, varying thicknesses, substrates and also different parametric processes including SHG, THG, or CARS play an important role. The task is not necessarily the output optimization but to understand the influence that the optical setup and sample geometry have on the output signal. In order to be able to make reliable statements on the material properties, this illustrative disentanglement of the interplay of optics and material on the one side and material properties on the other side can be very useful.

## IV. EXTENDING THE PHENOMENOLOGICAL APPROACH

### A. Modeling of finite crystals and thin films

In order to describe nonlinear optical microscopy signals generated in the bulk material or thin films more realistically, i.e., with finite thicknesses, multiple layers, and significantly tight focusing, the *Ansatz* of the modeling presented in Sec. II has to be extended.

In the following, we will use a setup, where a transparent material, with good air/material transmission is bonded to a (reflecting) high index substrate. These conditions can be found, for example, in a lithium niobate slab on silicon, a typical layer structure which is often experimentally examined^{11} and typical in modern (quantum) optics,^{4,6,25} or when investigating 2D materials,^{14} which is depicted in Fig. 1(a). The basic *Ansatz* of the used paraxial solution is taken from Boyd.^{19} Here, an incoming Gaussian fundamental beam is assumed and casted into a compact form as shown in Eq. (18),

Here, $a1$ 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 propagates in the $z$-direction and no angular dependence is present.

The well-known coupled wave equation for nonlinear phenomena, including the slowly varying envelope approximation (SVEA), i.e., $\u22022Aq\u2202z2\u226akq\u2202Aq\u2202z$, will be used in cylindrical coordinates, cf. Eq. (19),

where $Pq$ is the amplitude of the contribution of the nonlinear material polarization, which acts as a source for waves with a frequency of $\omega q=q\omega $.

Here, for $q$th-harmonic generation (QHG),

is the phase mismatch, where $q$ represents the order of the harmonic process, e.g., $q=2$ for SHG and $q=3$ for THG. Note that the refractive index $n=n(\lambda )$ is a quantity, which depends on the wavelength $\lambda $ or frequency $\omega $ of the light. Thus, as in almost every material, the refractive index is not constant, a certain phase mismatch is almost always present. Usually, one also defines the distance $lc,class$ up to which the waves are partially positively interfering and the harmonic signal is growing in the plane wave picture, cf. Eq. (21), where for plane waves, the $\Delta k$ is the only contribution to the signal oscillations, as

This equation can finally be compared to the extended ones of Eqs. (33) and (34). It remains to be seen how the experimental conditions generate deviations from this plane wave case; however, it gives a rough estimation of the values to be expected. In order to solve the differential equation, we use a trial solution which is similar to a Gaussian beam but includes also a z-dependent amplitude, cf. Eq. (22),

By inserting Eq. (22) into the differential equation (19), one obtains in good approximation^{19} a solution for the amplitude in form of an ordinary differential equation.^{19} That can be integrated and yields Eq. (23),

Here, $c=iq\omega 2nc\chi (q)$, where $\chi (q)$ is the $q$th-order nonlinear optical susceptibility and $z$, $z0$ are the start- and end point values of the material (and, therefore, the interaction region) determined relative to the focal position, which is zero in this parametrization. In the work of Boyd,^{19} presented in Sec. II, the integral is solved analytically via extending the borders to infinity in an analysis of thick crystals. It shows the surprising result that even phase-matched THG does not deliver any harmonic signal but needs to have a little positive phase mismatch in order to compensate for the additional Gouy-phase.

In the case of thin films, however, as depicted in Fig. 1(a), such a continuation of the integral is not possible; in fact, even typical scenarios of surface-near scans [cf. Fig. 1(b)] are not described by that approximation. It has to be solved numerically for finite boundaries. Furthermore, especially when samples are mounted on a strongly reflecting substrate like silicon, one has to consider the effects of reflection and transmission as well. Since in a typical case of lithium niobate (LN), the interface between LN and air has a high transmission $T\u22a5(\lambda \u22480.85\cdots 0.9\mu m)\u223c85%$, generally, it is sufficient to take into account only a small number $n$ of reflections; especially for the fundamental beam, interferences of the multiple reflections change the basic properties of our assumptions as they are related to the harmonic in a nonlinear fashion and have to be taken into account in detail; the issue here is that the polarization in the NLO material slab of thickness $d$ would be a large sum of positive direction propagating waves $z\u2192(z+2dn)$ and negative direction propagating waves $z\u2192(2nd\u2212z)$. This looks like the following Eq. (24), modulo the transverse parts:

Here, $cneg,n$ and $cpos,n$ are the amplitude coefficients for propagation in positive direction and negative direction, respectively, which are reflected $n$ times; they consist of the Fresnel coefficients for the corresponding reflections (and damping factors). In our case of a single transparent layer, we see that $cpos,0=1$, corresponding to the first incoming ray with frequency $\omega $, and $cpos,n=[rLN,Si(\omega )rLN,Air(\omega )]n$ corresponding to $n$ round-trips. The coefficients for the negative propagating part are then $cneg,n=rLN,Si(\omega )[rLN,Si(\omega )rLN,Air(\omega )]n$, where the extra reflection coefficient stems from the first reflection of the incoming beam at the substrate.

When searching for the $q$th harmonic, we need the $eiqkz$ and the complex conjugate term. Therefore, we can treat the counter- and co-propagating waves independently and both sums can be simplified: For a high transmission and $L<b$, we can group all elements of a respective sum together in a geometric series, when approximating the denominator in a (negligible) small damping factor $Adamp$; each series is dominated by a constant phase factor of $eik\xd72d$ and the corresponding Fresnel-coefficients then.

Thus, one can calculate the two contributions with two separate integrals for forward and backward propagating fundamental, where the generated harmonics may interfere afterwards. The corresponding thin film interference is included via multiplying the integrand with a geometric series to the power of $q$, which depends on the thickness $d$ of the slab. This dependence can also be implemented for the harmonic light, which shows only the linear interference which can be as well modulated with a geometric series as shown in Eq. (25),

where $kx$ is the fundamental or harmonic k-vector and $cx(\omega )=rLN,Si(\omega )rLN,Air(\omega )$ summarizes the Fresnel-coefficients corresponding to the reflections at the interfaces for the respective frequency $\omega $, i.e., the fundamental or harmonic frequency. In addition, depending on the point of observation, i.e., transmission or reflection geometry of the setup, the integrals obtain further Fresnel-coefficients, e.g., because the signal is detected in reflection geometry, the positive direction contribution gets a global $rLN,Si(\omega q)$, whereas the integral of the negative direction propagating waves naturally shows $rLN,Si(\omega )$, as discussed for $cneg,n$.

In general, further damping coefficients can be introduced as well. The series produces a characteristic pattern of oscillations and for the fundamental wave it holds that the larger $q$, the more significant the oscillation gets, as $q$ sits in the exponent.

One has to note that obviously the beam, when propagating, features an amplitude, which decreases with a polynomial not an exponential dependence. But as long as the transmission out of the material is quite large, the thin film interference can be very roughly approximated via a geometric series using a damping factor as mentioned above, as the damping by the transmission “losses” dominates. In comparison with the high transmission losses, the amplitude change due to the propagation, which evolves with $1(1+i2z/b)q\u22121$, does not result in a significant decrease in crystal thicknesses of order $t\u22642b$. Furthermore, when $L\u226bb$ the additional contribution due to reflected light is actually also getting increasingly weaker as the intensity far away from the focus also generates less signal contribution, so that the *Ansatz* remains viable. For different samples, especially with significantly smaller transmission coefficients between air and material, it might be necessary, however, to use the exact expression for the material polarization, as shown in Eq. (24). This requires a more involved calculation, which is presented in the Appendix and is as well only valid, as long as one can assume that the solution can be approximated by a Gaussian beam as presented in Eq. (22), requiring that we sum beams which do not differ significantly in the beam radius, i.e., due to propagation, otherwise an additional examination of the transverse parts will be necessary.

In fact, the geometric series obviously contributes oscillations much faster than that of the phase mismatch; they are correlated to the thin film interference of the fundamental and harmonic light, and clearly observable in Fourier-transforms^{11,26} of simulation and experiment, but concerning the coherence length, are more of cosmetic nature.

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

For THG, one needs to take into account the signal produced in air, or any other cladding material and, finally, sum over both, so that a more complex interference pattern is expected. This can be described by a sum of integrals, as sketched in Eq. (27),

Note that in such a case of contributions from different layers, the propagation to a common end point needs to be taken into account. Here, the air signal can be reflected, i.e., multiplied with the reflection coefficient or enter the crystal, i.e., multiplied with the same geometric series, the harmonic signal generated in the material encounters.

In this scenario, the exact ratio of the nonlinear optical susceptibility of air $\chi air(3)$ and lithium niobate $\chi LN(3)$ is of critical importance as they define the strength and phase of the respective contribution, as will be shown in Sec. V C.

### B. The phase mismatch and coherent interaction length in Gaussian beams

In fact, the changes in the phase evolution that stem from focusing, should be made quantifiable. As seen in different publications,^{11,19} the focusing conditions influence the observed coherent interaction length and the overall output signal. In some cases,^{19} focusing compensates for the positive phase mismatch, whereas it may also decrease the coherent interaction length, as observed for negative phase mismatch.^{11} In order to get a better understanding of this phenomenon, one can reformulate the integrand of Eq. (23) to the exponential form

Extracting the phase leads to Eq. (31), as

Demanding that the coherence length $lc$ is the length up to which wavelets are added, i.e., $\Delta \Phi (z)=!\pi $, we obtain Eq. (33),

Here, the approximation of Eq. (34) gives a good first estimation and is valid for $|2z/b|<1$, meaning the Rayleigh range is larger than the interaction region. That is the regime, where the $arctan\u2061(2z/b)$ is almost linear. The regime of Eq. (34) and the transition to Eq. (33) can be particularly well prepared and examined, when one is able to control the thickness of the sample, as in examinations of thin films/wedges, as is shown in Secs. V A–V D. Furthermore, one can observe that for high phase mismatches, meaning that $\Delta k\u226b(q\u22121)2/b$, the Gouy-phase contributes only minor corrections, so that the significance of the Gouy-phase is also dependent on the wavelength and the order $q$ of the process.

Note that when contemplating increasing crystal thicknesses or tighter focusing conditions, i.e., $L\u2265b/2$ the approximation is not valid and one better uses the full $arctan\u2061(2z/b)$-term. In the following, we will also use the complete transcendental equation and apply a common root finding algorithm, i.e., Newton’s method, to find an exact solution.

In fact, both formulas imply that on the one hand, a positive phase mismatch $\Delta k>0$, i.e., the missing quasi-momentum, and the Gouy-phase can partially compensate each other, so that the total effective phase mismatch in the relevant focal region can remain rather small, as is also observed in Boyd’s textbook.^{19} On the other hand, for $\Delta k<0$, the total phase mismatch is always increased, so that the coherence length for negative phase mismatch gets reduced for all wavelengths.^{11}

However, the exact phase evolution implies more modifications in comparison to the case without focusing. Therefore, in summary and in order to have a common interpretation of the nomenclature, one has to carefully differentiate the observables, which are commonly associated with the term “coherent interaction length,” or colloquially shortened to “coherence length.” The coherent interaction length $lc$ is the interaction distance up to which the generated wavelets interfere constructively. This corresponds to the distance of a minimum of local NLO intensity and the following maximum. Experimentally, as it is easier to read off, one often uses the distance from minimum to minimum.^{11} Thus, the coherence length, therefore, contains a total phase shift of the wavelets of $\pi $ and for the experimentally determined minimum-to-minimum distance, it is $2\pi $, as visualized in Fig. 5(a). One could call all these different distances in relation to their phase as $l\pi $, $l2\pi $. For the plane wave case, which is often implicitly assumed, as the phase $\varphi (z)=\Delta kz$ changes linearly, it naturally holds that $2\xd7l\pi =l2\pi $, so that both are directly linked. Moreover, **every single oscillation** in the plane wave case is the same; it does not matter if one looks at the first minimum-to-minimum distance or any other. By definition, this behavior also holds true for the linear approximation of Eq. (34). But, as was already shown in Eq. (33), for focused systems, the exact phase evolution is no longer linear. Here, due to the $arctan\u2061(2z/b)$-term, $2\xd7l\pi \u2260l2\pi $; for example, having a negative $\Delta k$, the $\u2212arctan\u2061(2z/b)$ term adds additional phase, so that $\pi $ is reached faster than in the plane wave case. For larger $z$, the $arctan\u2061(2z/b)$ converges to a constant, however, so that the following distances are larger, i.e., $2\xd7l\pi <l2\pi $. This relation then holds for all subsequent oscillations, until the $arctan\u2061(2z/b)$ is almost constant and the oscillation converges to the plane wave case. From the nonlinearity it also follows that $l\Delta \varphi $ also depends on the position of the focus relative to the optical material. When imagining the slab far away from the focus, the extra phase from the $arctan\u2061(2z/b)$ has no effects, but placing it within the Rayleigh range it does, so that even the placement of the focus, e.g., on the interface or at $z=zR$ produces different results.

To clarify which observable is dealt with, we will use $l\Delta \varphi $, e.g., $l2\pi $ for the first minimum-to-minimum distance. Furthermore, the entry point of the beam into the crystal will be given as a reference point, which will in general be the surface at $z0=0$.

## V. CASE STUDIES

In Subsections V A–V D, we study different NLO processes taking place in the model system of thin film LN (TFLN) on silicon and aim at working out on the influence of different parameters on the final signals. This thin film system is schematically depicted in Fig. 1(a), where we have used a wedged z-cut LN sample for variable thickness testing. All following studies will use this system, except the first, which deals with SHG in bulk LN. The values of $nLN(\lambda )$ are calculated via the Sellmeier equation, using the coefficients^{27} as tabulated in the work of Zelmon *et al.*

Furthermore, we will make some simplifications concerning the properties of LN. This concerns its anisotropy in the linear properties, i.e., its birefringence, and the nonlinear properties, i.e., the anisotropy and non-diagonal nature of its nonlinear coefficients. Here, we use the ordinary refractive index only and do neglect birefringence, as we are considering a z-cut sample and do not have a tight focus, but a $NA=0.45$. The latter, combined with good polarization control, also allows us to use a single effective nonlinear optical coefficient $deff$, e.g., $d211$ of the second order nonlinear optical susceptibility in the case of a Z(X,Y)–Z scattering geometry for SHG, as axial focal components are very weak in the paraxial regime.^{21}

### A. SHG in bulk LN

Before promoting LN to a thin film, we want to take a look at a LN bulk crystal once more. As already mentioned, the integral of the amplitude has to be solved numerically, in general. Even for large crystals, in the case of surface-near scans, the setup does not fulfill the condition of Eq. (7) allowing for infinite integration boundaries. These scans, however, are used, for example, for ferroelectric domain mapping.^{9} As the crystal is thick, the upper boundary can be set to infinity, whereas the entering point $z0$ is near zero. The experimental situation is depicted in Fig. 1(b).

The corresponding amplitude can be computed with Eq. (35),

Now, we can in principle observe co- and counter-propagating second-harmonic signal, where the latter has a very large, but positive phase-mismatch $\Delta k=2|k1|+|k2|$. As the crystal is very thick and bounded by an LN/air-interface on both sides, multiple reflections can be neglected due to high interface transmission, scattering, and damping of the signal.

In Fig. 2(a), one can observe the typical significant surface-near SHG.^{9,28} As can be extracted by the bulk integrals of Eq. (7) where the signal vanishes in an infinite bulk crystal, one sees SHG at the surface despite the negative (co-propagation) or very large (counter-propagation) positive mismatch. The counter-propagating signal is more than two orders of magnitudes weaker, though, simply due to the much shorter coherence length. As a side note, this also means that in a backreflection setup, depending on the NA and the actual crystal thickness and, therefore, collection efficiency of the diverging signal, it is possible to observe mainly reflected co-propagating light, as the transmission losses at the backside are roughly “only” $85%$. Here, the co-propagating light-signal would still be more than one order of magnitude larger, provided it can be collected efficiently.

Figure 2(a) also shows a case of a phase matched ($\Delta k=0$) scenario, which may be achieved in birefringent crystals by choosing specific wavelengths and scattering geometries. Since the corresponding signal is much stronger, it is scaled down, to be comparable; one can state that much more signal is generated, but the main peak still lies near the surface. However, it is not exactly centered at the surface. It validates the calculations of Eq. (17) that there remains a (weak) signal from the bulk, but also demonstrates, how strong the influence of the Gouy-phase is, functioning as an effective phase-mismatch here, when compared to a plane wave, where the signal would monotonically grow with crystal thickness.

It is interesting to consider the influence of the NA and the coherence length at this point. The NA, i.e., strength of focusing, on the one hand, is the dominant quantity for the focal length $b$ and for strong focusing this length dominates the effective interaction length and is also responsible for the rate, at which the Gouy-phase changes, as $\varphi Gouy\u221darctan\u2061(2zb)$. On the other hand, the phase mismatch $\Delta k$ determines the length scale on which constructive interference is possible. It is modified, however, by the Gouy-phase, as shown in Eq. (31). Therefore, the final results are determined by the ratio of $\Delta k$ to $2b$.

In this context, the dependence of the exact position of the signal peak needs further consideration. One would naively consider that to get the most signal, the focus of the pump beam should lie inside the material. In a closer examination, one can actually observe, however, that the maximum is mostly symmetrically distributed at the surface, cf. Fig. 2(b). Although the positions of these example cases, motivated by experimental work,^{11} only show tiny deviations, i.e., $\u223c10$ nm, from the surface, when choosing other wavelengths that produce smaller phase-mismatch values or even anomalous dispersion, the deviation can be significant, i.e., $\u223c1\mu m$ and will be easily observable. Therefore, we will also examine scenarios with decreased mismatch in the following, where $\Delta k$ is set to be within the same order of magnitude as the inverse of the focal length, i.e., $1b\u223c\Delta k$.

For $\Delta k\u226b2b$, the signal concentrates at the surface $zfoc\u22480$ and has a width within the same order of magnitude as the focal length, as also observed in experiments.^{9} The latter is due to the fact that with high $\Delta k$, i.e., $lc\u226ab$, the rapid $\Delta k$-oscillations prohibit constructive interference on a larger length-scale. Therefore, whenever the intensity maximum enters the crystal, the strongest signal can be detected. When the rest of the beam enters, the signal decreases again due to corresponding destructive interference.

One may consider that due to the Gouy-phase, the coherence length of counter-propagating light is (slightly) increased in dependence of the NA, but the mismatch for counter-propagating light is in general so large that the effect is almost negligible, i.e., the behavior for co- and counter-propagation is similar in this case.

When $\Delta k$ is significantly smaller, the deviation from the surface is modified by the interplay of the phase mismatch and the Gouy-phase, as both can contribute to the total phase evolution via $\Phi (z)=\Delta k\u2212arctan\u2061(2z/b)$. One can start from two different limiting cases.

When the focusing is significantly stronger or the phase mismatch very small, i.e., $\Delta k\u226a2b$ the Gouy-phase determines the evolution of the phase from $\u2212\pi /2$ to $\pi /2$. It helps to look at the real and imaginary part of the complex integrand, as given in Eq. (36),

As is illustrated in Fig. 3(a), the real part, i.e., the $cos[arctan\u2061(2z/b)]$ part, is symmetric and for $\Delta k=0$ always larger than zero, so that the optimum is situated in the center of the crystal. The $sin[arctan\u2061(2z/b)]$ part is antisymmetric, changes sign at the origin and is the largest when the focus is positioned at the surface. Then, the total optimum compromises both contributions and, as shown in an example calculation in Fig. 3(b) for the phase matched case, is near the surface but slightly inside the material for thick but finite crystals.

Thus, as a side note, the exact position depends on the exact thickness via the cutoff as well, which has to be introduced to regularize the logarithmically divergent integral of the imaginary part. The physical interpretation of this cutoff is nothing else as the crystal thickness; it co-determines the size and ratio of the symmetric and antisymmetric contributions. Obviously, the thinner the crystal, the more the damped oscillatory behavior of the integrand becomes visible in the position of the maximum in dependence of the crystal length, so that even maximum positions slightly above the surface are observable under certain conditions. In order to minimize these oscillatory effects, the numerical analysis of the trend of maximum positions in Fig. 4 is carried out with a realistic bulk crystal thickness of $t=3$ mm, so that the amplitude of the oscillation of the integrand at the upper border of the integral is two to three orders of magnitudes smaller than the maximum.

For scenarios in between, i.e., $2/b\u223c\Delta k$, the relation between $\Delta k$ and $b$ is crucial. When $\Delta k<0$, the phase-mismatch and the Gouy-phase add and the oscillating contributions lead to an overall weaker signal. Furthermore, for increasing phase mismatch, due to the introduced oscillations, cf. Fig. 3(a), the maximum of the integral of the $cos[\Delta kz+arctan\u2061(2z/b)]$ moves from deep into the crystal in the direction of the surface, avoiding the negative contribution left of the position of the first root, which appears on the left of the central maximum. The total maximum, therefore, moves continuously to the surface; it may even move slightly above the surface, when the total negative contribution of the symmetric part is larger, which depends on the cutoff, i.e., the crystal thickness. Then, with a further growing phase mismatch, the maximum position will oscillate around the surface and rapidly converge to the surface position for large $\Delta k$ reaching the case $\Delta k\u226b2b$, where the phase-mismatch dominates, cf. Fig. 2(b).

When $2/b\u223c\Delta k$ but $\Delta k>0$, the phase mismatch and the Gouy-phase can partially compensate. This leads to a very flat evolution of both integrands, as is demonstrated in Fig. 3(a). Although the maximum of the symmetric part moves slowly to the surface, the antisymmetric contribution almost vanishes in the focal region and the total maximum moves closer to the maximum position of the symmetric part. For larger $\Delta k$, the evolution gets sinusoidal again and converges in the direction of the surface again as in Fig. 2(b).

Finally, if in the former medium regime cases, the NA is increased, the dominance of the Gouy-phase grows again. Here, all optimum positions shift increasingly to the surface, cf. Figs. 4(a) and 4(b), as $2b$ increases while also the absolute scale $b$ decreases; the only exception is the $\Delta k>0$ case, where the Gouy-phase is still compensating for the positive phase mismatch.

In conclusion, this means that the position of the surface does not necessarily coincide with the position of strongest intensity, at least not when the phase mismatch is small or of similar magnitude compared to $2b$ and especially when it is positive, i.e., the chosen wavelength range shows anomalous dispersion. This may be relevant when accurate positions in depth need to be determined via SHG microscopy.

### B. SHG in TFLN

Next, we will consider SHG in the thin LN slab on silicon as shown in Fig. 1(a). Here, two components can be expected. We can have co-propagating and counter-propagating second-harmonic light, where the former plays the dominant role. We describe scans along the direction of growing crystal thickness; a probable experimental situation is depicted in Fig. 5(c).

The scans normally are considered with the focus lying at the surface $zin=0$, but can also be described for a variable focus depth. The effects of such differing focal depths are dealt with in the second part of this subsection.

#### 1. Signal evolution for fixed focus position

We now employ the discussed approach building on Eq. (23) with finite boundaries. We also take into account multiple reflections in the slab of thickness $L$ by including the nonlinear polarization generated by subsequent reflections, as shown in Eq. (24), as well as in the modulation of the fundamental and harmonic signal via the geometric series, as shown in Eq. (25), where the series belonging to the fundamental is set to the power of $q$. The resulting approximated integral is given in Eq. (37), whereas the complete integral and its scope of use is discussed in the Appendix,

As a reminder, $l2\pi $ is the distance where the phase of the signal generated in the following spatial slice is again constructive, i.e., the signal grows again. This corresponds to the thickness difference of two minima. In contrast to $lc=l\pi $ it is easier to directly extract $l2\pi $ from experimental data, as the fast oscillations, which increase the uncertainty when determining the position of the extrema, are less pronounced for the minima. In the following, $l2\pi $ will be computed, so that the results can be readily compared to experimental works.^{11} In general, $2lc\u2260l2\pi $ and thus, only if the phase is a linear function of space, then $2l\pi =l2\pi $, which is not given for non-negligible focus conditions in the vicinity of the focus. In Fig. 5(a) one now observes low frequency, but also high frequency oscillations, where the latter belong to the multi-reflections of the fundamental and SHG signal. The low frequency oscillations are due to the phase-mismatch and we see that the coherence length changes with wavelength as it is supposed to do with $n=n(\lambda )$. Considering the experimental values^{11} $l2\pi ,exp(850nm)\u22482.4$ *μ*m and $l2\pi ,exp(950nm)\u22483.55$ *μ*m, we can compare with the simulated data, as summarized in Table I. The plane wave case delivers $l2\pi (850nm)\u22482.73$ *μ*m and $l2\pi (950nm)\u22484.03$ *μ*m which is in the same order of magnitude but over $10%$ to large. Using the approximation of Eq. (34) we obtain $l2\pi (850nm)\u22482.37$ *μ*m and $l2\pi (950nm)\u22483.36$ *μ*m, which is closer to the experimental value and only about $5%$ off. The approximation works better for the lower wavelength, as the condition $L<b$ is better fulfilled. Finally, one can compute the exact solution, using the transcendental Eq. (33) (or extracting it from graph of the simulated data), yielding $l2\pi (850nm)\u22482.42$ *μ*m and $l2\pi (950nm)\u22483.50$ *μ*m, respectively. These values are even closer and lie in the scope of the measurement uncertainty, which was determined to be $\xb1100$ nm in the corresponding paper.^{11} Beware that the experimental publication^{11} calls the observable $2lc$, but it is in fact $l2\pi $, the thickness difference between the first and second minimum. The actual calculation of $2\xd7lc=2\xd7l\pi $ as twice the distance between the first minimum and maximum results in deviations, i.e., $2\xd7lc(850nm)\u22482.39$ *μ*m and $2\xd7lc(950nm)\u22483.40$ *μ*m. Here, the deviation increases with the coherence length, as the rate of phase-change, due to the Gouy-phase decreases with the distance, as the $arctan\u2061(x)$ converges to a constant. The counter-propagating light, on the other hand, does again not play a significant role in this context, because the coherence length is orders of magnitude smaller and the co-propagating light is efficiently reflected from the substrate. Thus, the counter-propagating signal is significantly weaker.

. | l_{2π} (850) (μm)
. | l_{2π} (950) (μm)
. |
---|---|---|

Experiment^{11} | (2.4 ± 0.1) | (3.55 ± 0.1) |

Plane wave | 2.73 | 4.03 |

Approx. Eq. (34) | 2.37 | 3.36 |

Full model | 2.42 | 3.50 |

#### 2. Dependence for varying focus position

The former results used $zf=0$, i.e., the focal plane is situated at the interface air/LN. In fact, for the case at hand, the overall trend of the signal does not change significantly, when changing the focal depth in a certain range. The resulting cross section for SHG in a thin LN slab, which was simulated with the same parameters as before, is shown in Fig. 5(b).

One can observe that the oscillations are continued when focusing into the material, although, in detail, deviations exist. But these are very small in this respect due to $b\u223cL$. Actually, the oscillation period is minimum, when the focus is placed at the center of the film ($zf\u2248t/2$, corresponding to a lower integral boundary of $zin\u2248\u2212t/2$), while slightly growing when focusing closer to either of the interfaces. Indeed, this was also predicted in our previous work by full numerical calculations.^{11} Due to the large computational demand of those calculations, though, this could only be calculated for a few selected cases, while here a full picture can be discussed readily with the semi-analytical model.

### C. THG and the interference of different layers

This section focuses on the generation of third-harmonic light. As practically all materials, including substrates or the air above the sample, show a third order non-linearity, one acquires distinct results that depend significantly more compared to SHG on the setup environment. As is illustrated in Fig. 7(c), for a typical microscopy setup, air (or oil/water) fills the space between lens and sample. As nonlinear optical microscopy operates with high pulse power, the air may as well interact with the light.

Thus, this additional layer has to be taken into account for the total picture of third-harmonic processes, too. Therefore, another integral is added to the computation, featuring the corresponding boundaries, refractive indices, and a different $\chi (3)$. To clearly see the interference effects, we chose the $\chi (3)$ for the air to be of the same order of magnitude, i.e., as large as the one of the LN slab.

#### 1. The role of the *χ*^{(3)}-relations

For THG, the signal of the air above the surface can play a significant role. In air, signal is generated, where on one hand the phase-mismatch is very small due to the low dispersion while on the other hand, part of this generated light enters the LN slab and undergoes the same multireflections as the signal generated in the slab; together they generate a more complex interference pattern governed by the phase relations of the different $\chi (3)$ values. As the nonlinear optical susceptibility can be a complex quantity in general, the result of interference patterns depends on the phase relation $\Delta \Phi =\u2212ilog\u2061(\chi air(3)/\chi LN(3))$ of the contributing $\chi (3)$-elements, cf. Fig. 6(a).

The data underline that the relative phase plays a significant role on the signal development. Although all signal curves show the characteristic low frequency oscillations, which stem from the phase-mismatch, especially the trend for small thicknesses is different. When the relation is a phase factor of $e0=1$, the signal initially drops, as the newly generated third-harmonic signal of LN is seemingly out of phase; for $e\pi =\u22121$ we get a growing signal. The signal drops or grows at first then, because the region, where LN-THG is generated, is growing, whereas its total phase is continuously changing; therefore, the interference pattern is shifted and also develops at small thicknesses at first. When the thickness reaches the size of the confocal parameter, the development becomes more stable, as less signal is generated in the air layer. Note that the figure shows the intensity and not the amplitudes; in the amplitudes, the oscillations of the zero and $\pi $ phase cases are (almost) symmetric. Interestingly, via the additional Gouy-phase, the susceptibilities with the same phase are (almost) destructively interfering at first.

This results demonstrate that the amplitude and phase relation of the nonlinear optical susceptibilities play an important role for the overall signal; this is true as well for the Gouy-phase. It does not only influence the coherence length, as already shown in different works,^{11,19} but also the phase of different interfering signals originating from different layers. Here, our own work^{26} shows in fact that the case of no additional significant phase-shift of the non-linearities is realized in the thin film LN on silicon material system, i.e., the blue signal evolution of Fig. 6(a).

#### 2. Note on the phase evolution

In order to obtain a more intuitive and general picture of the influence of the air signal, one can compare the signal patterns for varying $\chi (3)$ relations, as shown in Fig. 6, to a corresponding plane wave model and try to figure out different phase contributions. In the case of two interfering third-harmonic signals, the air signal has almost no phase mismatch, but presumably a collected Gouy-phase $\varphi 1$, so that one obtains

Where in contrast to the plane wave case in one medium, one finds an additional background $|A|2$ and an interference term $\u221dcos(\Delta kL/2\u2212\varphi 1)$. Obviously, for vanishing $L$ only the background remains. However $\varphi $ steers the further evolution; for non-focused systems, i.e., plane waves, one would expect no additional phase of the air signal and, therefore, a growing signal, similar to SHG with an additional background as shown in Fig. 6(b) in “blue,” although the interference term leads to a phase shift.

For the focused case, an additional phase is added. If the focal plane lies, for example, on the interface, the phase change can be extracted by computing the amplitude (and its phase) for the phase-matched beam in the air layer

Obviously, the plane wave model with $\varphi =\pi /2$ (orange), as in Fig. 6(a), fits perfectly to the case of no phase rotation of the susceptibility of the numerical analysis (blue), as in Fig. 6(b).

One can derive a corresponding relation for arbitrary harmonics and focal positions $z0$, which is given in Eq. (48),

One expects a phase that changes with focal position.

Noteworthy is also the fact that the total phase is the same for all processes for $z0=0$, although the higher harmonic beams pick up more Gouy-phase. This is due to the higher suppression, i.e., lower intensity of the wavelets with larger phase, i.e., wavelets created far from the focal region have low intensity.

The equation does not necessarily hold for SHG. However, we have calculated the corresponding integral for the real part in Eq. (17); doing the same for the imaginary part with the integrand $4z/b1+(2z/b)2$ yields the stem function $ln\u2061(1+(2z/b)2)$, which diverges. The ratio of the imaginary to real part is, therefore, infinite and $arctan\u2061(x\u2192\u221e)=\pi /2$, like for all the other harmonics.

#### 3. Third-harmonic signal for line-scans and cross sections

As expected, $lc$ changes for the different wavelengths, as given by the dispersion of the material, cf. Fig. 7(a). One can also see quite well the interference with a $\varphi 1=\pi /2$ signal due to air, as observed in the experiment.^{26}

Beware that due to the different peak and dip distribution for THG, it is necessary to specify the distance, which is chosen to measure the coherent interaction length, since the rate of phase change and, therefore, the interference pattern is dependent on the thickness and focal position, as was shown in Sec. IV B. In fact $l~2\pi $, which is shown in Fig. 7(a), marks a phase difference of the wavelets of $2\pi $ relative to the first dip, which corresponds to a total phase of $\pi $, i.e., the marked minimum-to-minimum distance corresponds to $l~2\pi =l3\pi \u2212l\pi $. It is in general larger than the preceding maximum-to-maximum distance, i.e., $l~2\pi >l2\pi $; however, $l~2\pi $ can be directly compared with the equivalent minimum-to-minimum values of the experimental work^{26} of Amber *et al.*

A corresponding cross section into the wedge, as depicted in Fig. 7(b), shows a significantly different signal pattern in contrast to SHG.

Since the total phase of the air-signal and the LN-signal change with the focal position, there is also a signal pattern in the depth-direction, which can be observed in Fig. 7(b). This can be seen in Eq. (48), where the focal position as an argument of the $arctan\u2061(x)$ determines the phase of the air signal; equally, the translation of the focal position $z\u2192z\u2212zfoc$ also shifts the integral boundaries of the LN signal. Approximating again with plane waves, one obtains Eq. (51),

Thus, we obtain again a similar oscillation period in depth as for varying thickness, to which in principle both phase evolutions, i.e., of the signal generated in air and of the signal generated inside the material, contribute.

These are governed by the phase mismatch and the Gouy-phase contributions inside and outside the crystal. The z-dependent part of the argument of the cosine, shown in the approximated expression in Eq. (52), i.e.,

allows for approximating the period as $ldepth,2\pi \u22482\pi \Delta k+O(1b)$. The similarity of both periods is visible in Fig. 7(b).

### D. Models vs the experiment

There is a quite rich experimental database for employing the modeling scenarios shown before. First of all, one can find a lot of depth scans done in the context of mapping ferroelectric domain walls, which correspond to the case of a half-infinite crystal. For example, in past works^{9,28} of the authors, one can find depth scans for high numerical apertures of $NA=0.95$ that show the signal peak centered at the surface as calculated in Sec. V A and depicted in Fig. 2.

Then, there are second-harmonic measurements^{11} for wedged thin film lithium niobate on silicon, which agree very well with our simulations, which will be used in the following to compare the two models. One can also compare the depth-evolution as well as the changes of $lc$ with the focusing conditions, cf. Figs. 8(a) and 8(b).

In Fig. 8(a), one observes an increasing match for the experimental $l2\pi $ and the simulated ones, when using the plane wave, first-order-approximation and exact case, in that order. In comparison to the full simulation used in the source-paper,^{11} the effort for the calculation is much lower; that means several seconds to calculate $l2\pi $ for a set of wavelengths, since here only the transcendental Eq. (33) has to be solved, up to minutes for several line-scans on a typical personal computer (Intel Core i7-6500 @2.5Ghz, 8 GB RAM, Intel HD Graphics 520) for the semi-analytical model. In comparison with that, one needs, depending on the targeted precision, at least several minutes to hours for the full numerical calculations for wavelength on a virtual machine running on a server cluster (24 virtual CPU cores @2.4 GHz, 512 GB RAM). A more detailed discussion of the performances of both models is given in a related publication.^{26}

Comparing the gathered data for different NAs, as depicted in Fig. 8(b), one can clearly observe the dependence of the coherence-length, or more specifically $l2\pi $ in this case, on the NA, which is not implemented in the plane wave case.

As a note, the full 3D-numerics^{11} use a focal spot in the center of the layer, whereas the semi-analytical solutions are centered at the surface; in general, the extended models are much better suited as is the plane wave case, although for $NA>0.45$, the exact solution and the first order approximation diverge. Note that the paraxial approximation of the Gauss-beams becomes worse for high NAs, however as $nLN\u22482.2$, the effective NA is reduced, so that $sin[arcsin\u2061(0.8/2.2)]\u22480.36\u22480.8/2.2$ meaning $sin(x)\u2248x$ is valid and $cos[arcsin\u2061(0.8/2.2)]\u22480.93$, so $cos(x)\u22481$ is roughly valid. That means that the paraxial approximation is still sufficiently satisfied, which can be illustrated when comparing the thickness-evolution of the signals in Fig. 9, which use a numerical aperture of $NA=0.8$ and yield a good agreement. Thus, in other words, for an effective NA of $NAeff\u2248NAn\u22480.36$, the approximation is applicable.

One additional observation can be inferred from the properties of the experimental light source, which shows a bandwidth of $\u223c10$ nm. The results of Figs. 8 and 9 as well as of Ref. 26 show that this does not influence the results in a significant way. Here, the bandwidth compression, due to the nonlinearity of the process does on the one side suppress the effects, but more importantly, the quantities, i.e., the oscillation due to the phase-mismatch and multi-interference are then a superposition of different frequency components, which may reduce the contrast, but for small bandwidths and depths not the central values. Furthermore, concerning the fast oscillations, the reflectivity of the resonator/thin film is very small, so that the loss of contrast due to the bandwidth is of second order. For larger depths, the unequal propagation will eventually show more impact.

## VI. CONCLUSION AND OUTLOOK

In this work, we develop and discuss a semi-analytical modeling *Ansatz*, which can be used for the simulation of nonlinear optical processes in stratified, layered, nonlinear media. It is based on an analytical *Ansatz*,^{19} which solves the approximated nonlinear wave equation using Gaussian beams with spatially dependent amplitudes. This *Ansatz* is re-casted to be suitable for finite bulk and (layered) thin film systems. At the same time, we try to conserve the instructive potential of the parameters of the analytical results.

This concerns especially the need for a modified description of the coherent interaction length. It is shown that it is dependent on the focus position and the actual total phase which is acquired, i.e., subsequent oscillations may have different period lengths. A corresponding mathematical formulation of the phase evolution is derived from the differential equation, which follows from the Gaussian *Ansatz*. Thus, it gets obvious that focusing is the reason for the re-interpretation because the Gouy-phase of the fundamental beam introduces a non-linearity in the vicinity of the focus. Therefore, a specific description of this parameter is necessary for comparison, i.e., relative position of the focus, oscillation order, and/or amount of total acquired phase. Only oscillations far away from the focal spot converge to the plane wave case of constant $lc$, as the additional phase due to focusing, the Gouy-phase, tends to be a constant for large distances.

In the second part, we apply the modeling *Ansatz* to several example cases. The *Ansatz* is applied for surface-near bulk scans as well as for thin layers on a reflective substrate.

In the surface-near SHG bulk scans and for large phase-mismatches ($\Delta k\u226b1/b$), we observe a second-harmonic signal, which is concentrated at the interface. Such surface-near SHG was already observed in several experiments.^{9,28} Furthermore, for small phase-mismatches, the peak position of the second-harmonic signal can substantially ($\u223c1$ *μ*m) deviate from the surface.

For thin film materials for SHG and THG, we observe good agreement with rigorous 3D numerics^{18} and experiments.^{11,26} Specifically for THG, the air layer also produces a signal, which interferes with the signal of the nonlinear material such that a specific interference pattern is visible. It is determined by the phase relations of both signals, which depends on the phase between their susceptibilities but also on the difference in the acquired Gouy-phase in air and the material.

Thus, the modeling *Ansatz* can perform very well in comparison to rigorous numerics and delivers an instructive notion of the contributing physical phenomena. Furthermore, the calculations are very fast, i.e., in the order of a minute when using an ordinary personal computer device for a full simulation of several wavelengths and even less for directly calculating specific parameters like the coherent interaction length.

Concerning future research, our approach, or extensions of it, allow for the examination of further nonlinear optical processes besides harmonic generation, provided that they can be mapped to an adapted version of the *Ansatz*. An interesting candidate for examination is the four wave mixing (FWM) process of Coherent Anti-Stokes Raman Scattering (CARS), which constitutes an important tool for material characterization.^{16}

In summary, although we point out the limits of the approach, which is bound to the form and parameter range of the analytical solution and thus cannot, for example, describe focus distortion,^{22,23} it can be a practical complement to perform nonlinear optical analysis. It opens up the possibility to easily identify the different effects, which are at work, yielding an educational and easy-to-use toolkit for the presented experimental examinations, which facilitates to develop the setup and extract specific information from the experiment.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) through Project Nos. 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, the authors would like to acknowledge the excellent support by the Light Microscopy Facility, a Core Facility of the CMCB Technology Platform at TU Dresden, where the SHG/THG analysis was performed.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

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

## DATA AVAILABILITY

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

### APPENDIX: INTEGRAL EQUATION FOR GENERAL THIN FILMS

The integral equation, i.e., (37), given in Sec. V is an approximation of the given scenario, in that sense, as the multi-interference due to the multi-reflections is captured in two geometric series, one for the fundamental and one for the harmonic signal. As discussed in Sec. IV, this approximation is especially justified for a thin film with small reflectivity $R=rLN,SirLN,Air$ and a thickness $t\u226ab$, which guarantees that the change of amplitude due to the polynomial dependence of the denominator is not too significant. On the other hand, if one works with very large $t\u226bb$, the amplitude may have already decreased significantly, such that reflections of the fundamental can be neglected altogether and the reflected parts of the harmonic beam have eventually diverged so much that they do not significantly contribute in back-reflection geometry, which is described by the original textbook integral for finite boundaries in Eq. (35). For some medium regimes, where the thickness and reflectivity are not too large, one may find deviations from Eq. (37), however, and can use the full integral given by Eq. (A1) instead,

Here, $\chi eff$ is the effective optical non-linearity, which is determined by the type of crystal and its orientation relative to the orientation of the incoming polarization of the light beam. Furthermore, $cpos,n=[rLN,Si(\omega )rLN,Air(\omega )]n$ and $cneg,n=rLN,Si(\omega )[rLN,Si(\omega )rLN,Air(\omega )]n$, as well as $cpos,k(\omega q)=rLN,Si(\omega q)[rLN,Air(\omega q)rLN,Si(\omega q)]k$ and $cneg,k(\omega q)=[rLN,Air(\omega q)rLN,Si(\omega q)]k$ determine the corresponding Fresnel-coefficients for the fundamental and harmonic light, respectively. $Nf$ and $Nh$ corresponds to the included number of the infinite many round-trips for the fundamental and harmonic light and can be chosen to be included for the calculation according to the reflectivity of the layer. For $L\u226ab$, the summation of the harmonics and fundamental wavelets tend again to a pure geometric series. For a small reflectivity and $L<b$, only a few round-trips do significantly contribute and one can approximate the ratio of succeeding terms of the sums with $R\xd7ei2kL$ and treat the geometric series as modulating the leading order, as was presented in Eq. (37). Here, the small reflectivity guarantees the fast convergence of the geometric series and that we can approximate the trial solution as a Gaussian beam. In fact, in order to be valid, the beam radii of the sum of the Gaussian beams, must be similar, which is valid as long as the layer is thin, i.e., of the order of the confocal parameter and the reflectivity is low, i.e., contributing round-trips are of the order of the confocal parameter. Thus, the main effect of Eq. (A1) consists of additional phase contributions, which were neglected in the geometric series approach.

## REFERENCES

*et al.*, “