Coherent nonlinear optical μ-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.

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 control8 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-microscopy9–15 as well as via Coherent Anti-Stokes Raman Scattering16,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 available18 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.

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

E(z,t)=aq(z)×(ei(kqzωqt)+c.c.).
(1)

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

aeiΔkzdz.
(2)

Here, the phase mismatch Δk=i,okiko, i.e., the difference between wave-vectors of the incoming (ki) and produced waves (ko), 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,

2ikqAqz+T2Aq=ωq2ε0c2Pqe(iΔkz),
(3)
A(z,r)=aq(z)1+i2zbe(qr2w02(1+i2z/b))ei(kqzωqt)+c.c.
(4)

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),

adzeiΔkz(1+i2zb)q1p.w.Lb1.
(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

I(Δk,z0,z)=eiΔkz(1+i2zb)q1dz
(6)
={0,Δk0,b22π(q2)!(bΔk2)q2ebΔk/2,Δk>0.
(7)

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 Δk=0 is only valid for q3, 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Δk=1. Therefore, one can directly integrate the remaining 1(1+i2z/b)q1 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 Δk to vanish. For THG and all higher harmonics, the transition is smooth, but not for SHG, where the result is 0 or bπ, depending on the order of Δk0 and q2. One can use the direct approach of the integral and successively decrease the symmetric boundaries to zero, whereas the phase term is isolated

11+i2z/bdz=[ln(1+2iz/b)2i/b]
(8)
=limz1z2[ln(1+(2z/b)2ei×arctan(2z/b))2i/b]z2z1
(9)
=limz1z2[ln(1+(2z1/b)2ei×arctan(2z1/b)1+(2z2/b)2ei×arctan(2z2/b))b2i]
(10)
=z2=z1limz1b2[arctan(2z1b)arctan(2z1b)]=πb2.
(11)

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

Re(E)12dz11+i2z/b+112iz/b
(12)
=12dz1i2z/b+1+i2z/b1+(2z/b)2
(13)
=12dz21+(2z/b)2
(14)
=12dxb1+x2
(15)
=12[b×arctan(x)]=πb/2,
(16)
Im(E)12dz4z/b1+(2z/b)2=0.
(17)

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 work11 of the authors. To obtain more realistic predictions, however, one also needs to utilize numerical strategies.

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 depolarization9,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 work24 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 ξ=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.

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 examined11 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),

Af(r,z)=a11+i2z/b×exp(r2w02(1+i2z/b)).
(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.

FIG. 1.

Illustration of (a) the z-cut LN wedge sample used in the case study to examine the thickness dependence of the nonlinear signal for thin films. The angle α is exaggerated to improve visibility. Subfigure (b) visualizes the geometry of a depth scan into a semi-infinite crystal, which is used as the basic geometry to record the data displayed in Figs. 2 and 4. Note that zf=0 corresponds to a focus position at the top interface and positive zf-values correspond to a focus position zf inside the crystal. From the viewpoint of the Gaussian beam of Eq. (18), which is centered at the origin, the integral boundaries (or the integration variable) have to be adapted for the actual calculation, e.g., zin=zf as the adapted lower boundary.

FIG. 1.

Illustration of (a) the z-cut LN wedge sample used in the case study to examine the thickness dependence of the nonlinear signal for thin films. The angle α is exaggerated to improve visibility. Subfigure (b) visualizes the geometry of a depth scan into a semi-infinite crystal, which is used as the basic geometry to record the data displayed in Figs. 2 and 4. Note that zf=0 corresponds to a focus position at the top interface and positive zf-values correspond to a focus position zf inside the crystal. From the viewpoint of the Gaussian beam of Eq. (18), which is centered at the origin, the integral boundaries (or the integration variable) have to be adapted for the actual calculation, e.g., zin=zf as the adapted lower boundary.

Close modal

The well-known coupled wave equation for nonlinear phenomena, including the slowly varying envelope approximation (SVEA), i.e., 2Aqz2kqAqz, will be used in cylindrical coordinates, cf. Eq. (19),

2ikqAqz+T2Aq=ωq2ε0c2Pqe(iΔkz),
(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 ωq=qω.

Here, for qth-harmonic generation (QHG),

Δk=qkωnωkqωnqω
(20)

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(λ) is a quantity, which depends on the wavelength λ or frequency ω 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 Δk is the only contribution to the signal oscillations, as

lc,class=πΔk=λω2q1nωnqωwith(q=2,SHGq=3,THG).
(21)

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),

Aq(r,z)=aq(z)1+i2z/b×exp(qr2w02(1+i2z/b)).
(22)

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

aq(z;z0)=cz0za1q×eiΔkz(1+i2z/b)q1dz.
(23)

Here, c=iqω2ncχ(q), where χ(q) is the qth-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(λ0.850.9μm)85%, 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(z+2dn) and negative direction propagating waves z(2ndz). This looks like the following Eq. (24), modulo the transverse parts:

PNL[cpos,0eikz1+i2z/b+n=1(cneg,neik(2ndz)1+i2[2ndz]/b+cpos,neik(z+2nd)1+i2[z+2nd]/b)+c.c.]q.
(24)

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 ω, and cpos,n=[rLN,Si(ω)rLN,Air(ω)]n corresponding to n round-trips. The coefficients for the negative propagating part are then cneg,n=rLN,Si(ω)[rLN,Si(ω)rLN,Air(ω)]n, where the extra reflection coefficient stems from the first reflection of the incoming beam at the substrate.

When searching for the qth 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×2d 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),

0qn=11qforq<1
(25)
q=cx×ei2kxd,
(26)

where kx is the fundamental or harmonic k-vector and cx(ω)=rLN,Si(ω)rLN,Air(ω) summarizes the Fresnel-coefficients corresponding to the reflections at the interfaces for the respective frequency ω, 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(ωq), whereas the integral of the negative direction propagating waves naturally shows rLN,Si(ω), 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)q1, does not result in a significant decrease in crystal thicknesses of order t2b. Furthermore, when Lb 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-transforms11,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 χ(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),

cz0zmaxdzcairz0=zLNdz+cLNzLNzmaxdz.
(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 χair(3) and lithium niobate χ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.

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

a(z)=eiΔkz(1+i2z/b)q1dz
(28)
eiΔkz(1+i2z/b)q1=eiΔkz[1i2z/b1+4z2/b2]q1
(29)
=eiΔkzei×(q1)×arctan(2z/b)(1+4z2/b2)q1.
(30)

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

ϕ(z)=Δkz+(q1)×arctan(2z/b)GouyPhase
(31)
Δkz(q1)×2z/b.
(32)

Demanding that the coherence length lc is the length up to which wavelets are added, i.e., ΔΦ(z)=!π, we obtain Eq. (33),

ΔΦ=π=ϕ(x1)ϕ(x0),wherelc=x1x0,
(33)
withx0=0,Lblc|πΔk(q1)×2/b|.
(34)

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(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 AV D. Furthermore, one can observe that for high phase mismatches, meaning that Δk(q1)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., Lb/2 the approximation is not valid and one better uses the full arctan(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 Δ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 Δ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 π and for the experimentally determined minimum-to-minimum distance, it is 2π, as visualized in Fig. 5(a). One could call all these different distances in relation to their phase as lπ, l2π. For the plane wave case, which is often implicitly assumed, as the phase ϕ(z)=Δkz changes linearly, it naturally holds that 2×lπ=l2π, 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(2z/b)-term, 2×lπl2π; for example, having a negative Δk, the arctan(2z/b) term adds additional phase, so that π is reached faster than in the plane wave case. For larger z, the arctan(2z/b) converges to a constant, however, so that the following distances are larger, i.e., 2×lπ<l2π. This relation then holds for all subsequent oscillations, until the arctan(2z/b) is almost constant and the oscillation converges to the plane wave case. From the nonlinearity it also follows that lΔϕ 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(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Δϕ, e.g., l2π 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.

In Subsections V AV 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(λ) are calculated via the Sellmeier equation, using the coefficients27 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 

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),

aq=cLNz0a1q×eiΔkz(1+i2z/b)q1dz.
(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 Δ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.

FIG. 2.

Results for surface-near second-harmonic depth scans [cf. Fig. 1(b)] for co- and counter-propagating signal in a semi-infinite LN crystal. The fundamental wavelength of 850 nm and a NA of 0.45 are employed for this simulation, alongside an actual crystal thickness of t=3 mm. One observes the experimentally well-known signal peak at the surface, typical to large phase mismatch. A hypothetical phase-matched case is also given for discussion; (a) shows an overview down to 50μm into the crystal, whereas (b) shows a close-up of the interface region. Here, with this moderate focusing, one observes that the phase-matched case shows the peak to be well inside the material, while the other cases have their peak approximately at the surface.

FIG. 2.

Results for surface-near second-harmonic depth scans [cf. Fig. 1(b)] for co- and counter-propagating signal in a semi-infinite LN crystal. The fundamental wavelength of 850 nm and a NA of 0.45 are employed for this simulation, alongside an actual crystal thickness of t=3 mm. One observes the experimentally well-known signal peak at the surface, typical to large phase mismatch. A hypothetical phase-matched case is also given for discussion; (a) shows an overview down to 50μm into the crystal, whereas (b) shows a close-up of the interface region. Here, with this moderate focusing, one observes that the phase-matched case shows the peak to be well inside the material, while the other cases have their peak approximately at the surface.

Close modal

Figure 2(a) also shows a case of a phase matched (Δ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 ϕGouyarctan(2zb). On the other hand, the phase mismatch Δ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 Δ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., 10 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., 1μm and will be easily observable. Therefore, we will also examine scenarios with decreased mismatch in the following, where Δk is set to be within the same order of magnitude as the inverse of the focal length, i.e., 1bΔk.

For Δk2b, the signal concentrates at the surface zfoc0 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 Δk, i.e., lcb, the rapid Δ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 Δ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 Φ(z)=Δkarctan(2z/b). One can start from two different limiting cases.

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

aqz0ei(Δkxarctan(2xb))1+(2xb)2dx=z0cos(Δkxarctan(2xb))1+(2xb)2+isin(Δkxarctan(2xb))1+(2xb)2dx.
(36)

As is illustrated in Fig. 3(a), the real part, i.e., the cos[arctan(2z/b)] part, is symmetric and for Δk=0 always larger than zero, so that the optimum is situated in the center of the crystal. The sin[arctan(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.

FIG. 3.

Example simulations with confocal parameter b=5 and crystal thickness of t=0.3 mm. Subfigure (a) shows the real and imaginary part of the integrand of Eq. (36) for different values of phase mismatch Δk. In subfigure (b), you can see an example of the absolute squares of the computed integrals for Δk=0, where the contributions of real and imaginary part and the total signal are illustrated separately.

FIG. 3.

Example simulations with confocal parameter b=5 and crystal thickness of t=0.3 mm. Subfigure (a) shows the real and imaginary part of the integrand of Eq. (36) for different values of phase mismatch Δk. In subfigure (b), you can see an example of the absolute squares of the computed integrals for Δk=0, where the contributions of real and imaginary part and the total signal are illustrated separately.

Close modal

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.

FIG. 4.

Surface-near SHG depth scans [cf. Fig. 1(b)] for different values of Δk and different NAs in a semi-infinite LN crystal of t=3 mm thickness. The fundamental wavelength of 850 nm and a NA of 0.45 is employed for the simulation and the original Δk values for co- (Δkco) and counter-propagation (Δkcount) at 850 nm are changed to Δk, to illustrate possible different signal behavior. One observes that (a) for Δk=0, the maximum moves from inside to the surface with increasing NA; (b) for small negative mismatch (Δk0.04Δkco), the maximum is slightly inside or outside the surface, to evade the additional mismatch due to the Gouy-phase, and moves to the surface for increased focusing; (c) for small positive Δk (Δk0.01Δkcount), the maximum also lies inside the crystal and moves inward for higher NAs as the Gouy-phase compensates the mismatch until its gets too large and can be approximated with case (a).

FIG. 4.

Surface-near SHG depth scans [cf. Fig. 1(b)] for different values of Δk and different NAs in a semi-infinite LN crystal of t=3 mm thickness. The fundamental wavelength of 850 nm and a NA of 0.45 is employed for the simulation and the original Δk values for co- (Δkco) and counter-propagation (Δkcount) at 850 nm are changed to Δk, to illustrate possible different signal behavior. One observes that (a) for Δk=0, the maximum moves from inside to the surface with increasing NA; (b) for small negative mismatch (Δk0.04Δkco), the maximum is slightly inside or outside the surface, to evade the additional mismatch due to the Gouy-phase, and moves to the surface for increased focusing; (c) for small positive Δk (Δk0.01Δkcount), the maximum also lies inside the crystal and moves inward for higher NAs as the Gouy-phase compensates the mismatch until its gets too large and can be approximated with case (a).

Close modal
FIG. 5.

SHG signal of (a) a line-scan for different LN thicknesses with zf=0. The wavelength is varied from 850 to 950 nm and the NA is 0.45. One observes that the coherence length is influenced similar to experiments.11 In (b) the corresponding cross section of the wedge is shown for λ=850 nm. In the depicted depths the evolution changes only marginally. In (c) the measurement setup is sketched, i.e., for (a) a line-scan on the surface was carried out.

FIG. 5.

SHG signal of (a) a line-scan for different LN thicknesses with zf=0. The wavelength is varied from 850 to 950 nm and the NA is 0.45. One observes that the coherence length is influenced similar to experiments.11 In (b) the corresponding cross section of the wedge is shown for λ=850 nm. In the depicted depths the evolution changes only marginally. In (c) the measurement setup is sketched, i.e., for (a) a line-scan on the surface was carried out.

Close modal

For scenarios in between, i.e., 2/bΔk, the relation between Δk and b is crucial. When Δ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[Δkz+arctan(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 Δk reaching the case Δk2b, where the phase-mismatch dominates, cf. Fig. 2(b).

When 2/bΔk but Δ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 Δ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 Δ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.

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,

S(L,zf=0;Δk)tAir,LN(ω)tLN,Air(ωq)×[0LrLN,Si(ωq)a1qeiΔkz(1+i2z/b)q1dz+0LrLN,Siq(ω)a1qeiΔk(z+L)(1+i2(z+L)/b)q1dz]×(11cx(ω)ei2kL)q(11cx(ωq)ei2kqL).
(37)

One observes in Fig. 5(a), a similar behavior as in the experiment.11 

As a reminder, l2π 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π it is easier to directly extract l2π 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π will be computed, so that the results can be readily compared to experimental works.11 In general, 2lcl2π and thus, only if the phase is a linear function of space, then 2lπ=l2π, 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(λ). Considering the experimental values11l2π,exp(850nm)2.4μm and l2π,exp(950nm)3.55μm, we can compare with the simulated data, as summarized in Table I. The plane wave case delivers l2π(850nm)2.73μm and l2π(950nm)4.03μm which is in the same order of magnitude but over 10% to large. Using the approximation of Eq. (34) we obtain l2π(850nm)2.37μm and l2π(950nm)3.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π(850nm)2.42μm and l2π(950nm)3.50μm, respectively. These values are even closer and lie in the scope of the measurement uncertainty, which was determined to be ±100 nm in the corresponding paper.11 Beware that the experimental publication11 calls the observable 2lc, but it is in fact l2π, the thickness difference between the first and second minimum. The actual calculation of 2×lc=2×lπ as twice the distance between the first minimum and maximum results in deviations, i.e., 2×lc(850nm)2.39μm and 2×lc(950nm)3.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(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.

TABLE I.

Comparison between the experimental data of l2π for two selected wavelengths in the LN thin film and corresponding computed data, using the simulation data, as well as the approximation of Eq. (34) and the plane wave approximation.

l2π (850) (μm)l2π (950) (μm)
Experiment11  (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 
l2π (850) (μm)l2π (950) (μm)
Experiment11  (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 

A detailed comparison of the coherent interaction lengths for different wavelengths and an NA=0.45 as a benchmark case is presented in Fig. 8(a), where it is compared to experimental data and data sets of more rigorous numerical simulations.11 They are discussed in Sec. V D.

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 bL. Actually, the oscillation period is minimum, when the focus is placed at the center of the film (zft/2, corresponding to a lower integral boundary of zint/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.

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 χ(3). To clearly see the interference effects, we chose the χ(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 χ(3) values. As the nonlinear optical susceptibility can be a complex quantity in general, the result of interference patterns depends on the phase relation ΔΦ=ilog(χair(3)/χLN(3)) of the contributing χ(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π=1 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 π phase cases are (almost) symmetric. Interestingly, via the additional Gouy-phase, the susceptibilities with the same phase are (almost) destructively interfering at first.

FIG. 6.

THG signal strength (a) in dependence of different phases between the χ(3) of air and LN calculated with the Gaussian Ansatz. The scan corresponds to a zin=0 line-scan in Fig. 7(c). Here, the amplitude of the susceptibility of air is equal to that of LN. The importance of the phase for the signal evolution is well observable; (b) shows an analog calculation for the layer system in the plane wave model of Eq. (40). To obtain the (experimentally observed)26 results of the Gaussian beam simulations with no phase-shift, the air-signal needs to have a phase of π/2, i.e., the blue curve of (a) corresponds to the orange curve of (b).

FIG. 6.

THG signal strength (a) in dependence of different phases between the χ(3) of air and LN calculated with the Gaussian Ansatz. The scan corresponds to a zin=0 line-scan in Fig. 7(c). Here, the amplitude of the susceptibility of air is equal to that of LN. The importance of the phase for the signal evolution is well observable; (b) shows an analog calculation for the layer system in the plane wave model of Eq. (40). To obtain the (experimentally observed)26 results of the Gaussian beam simulations with no phase-shift, the air-signal needs to have a phase of π/2, i.e., the blue curve of (a) corresponds to the orange curve of (b).

Close modal
FIG. 7.

(a) Simulation of THG signal for increasing material thickness for different fundamental wavelengths, corresponding to a line-scan on the surface (zin=0) as depicted in (c). One observes an evolution of the coherence length which agrees well with experimental data.26 The corresponding cross-section scan of the wedge (b) illustrates the dependence of the interference pattern on the focal depth. In (c), the measurement setup is sketched, i.e., for (a) a line-scan on the surface was carried out.

FIG. 7.

(a) Simulation of THG signal for increasing material thickness for different fundamental wavelengths, corresponding to a line-scan on the surface (zin=0) as depicted in (c). One observes an evolution of the coherence length which agrees well with experimental data.26 The corresponding cross-section scan of the wedge (b) illustrates the dependence of the interference pattern on the focal depth. In (c), the measurement setup is sketched, i.e., for (a) a line-scan on the surface was carried out.

Close modal
FIG. 8.

(a) Comparison of the evolution of l2π for different calculation schemes with experimental data.11 Then, in (b), the evolution of l2π for different NAs and calculation schemes are shown in relation to experimental data. For both, the full vectorial model11 is used besides the plane wave (P.W.) and the (semi-)analytical calculations of this paper, which can be compared to experimental values.11 Note that the legend shows how the corresponding model calculates lc=lπ. Here, we determined the confocal parameter b using tan(β)=limzw(z)/z=2w0/b and the relation b2=πw02nLNλ, so that with NA=sin(α)=nLNsin(β)=nLNsin(arctan(2w0/b))nLN×2w0/b, i.e., b=2λπnLN(tan(arcsin(NA/nLN))22λnLNπNA2.

FIG. 8.

(a) Comparison of the evolution of l2π for different calculation schemes with experimental data.11 Then, in (b), the evolution of l2π for different NAs and calculation schemes are shown in relation to experimental data. For both, the full vectorial model11 is used besides the plane wave (P.W.) and the (semi-)analytical calculations of this paper, which can be compared to experimental values.11 Note that the legend shows how the corresponding model calculates lc=lπ. Here, we determined the confocal parameter b using tan(β)=limzw(z)/z=2w0/b and the relation b2=πw02nLNλ, so that with NA=sin(α)=nLNsin(β)=nLNsin(arctan(2w0/b))nLN×2w0/b, i.e., b=2λπnLN(tan(arcsin(NA/nLN))22λnLNπNA2.

Close modal

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 work26 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 χ(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 ϕ1, so that one obtains

I=|A×eiϕ1+B×0LeiΔkzdz|2,withA,BR+
(38)
=|Aeiϕ1+B×eiΔkL/22Δk×sin(ΔkL/2)|2=|A|2+|B2Δksin(ΔkL/2)|2
(39)
+AB4Δkcos(ΔkL/2ϕ1)sin(ΔkL/2).
(40)

Where in contrast to the plane wave case in one medium, one finds an additional background |A|2 and an interference term cos(ΔkL/2ϕ1). Obviously, for vanishing L only the background remains. However ϕ 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

aair01(1+i2z/b)2dz
(41)
=[b2i11+i2z/b]0
(42)
=b2i=b2eiπ2
(43)
ϕ1=π2.
(44)

Obviously, the plane wave model with ϕ=π/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),

aair,q(z0)z01(1+i2z/b)q1dz
(45)
=[b2i(q2)1(1+i2z/b)q2]z0
(46)
=b2(q2)eiπ2i(q2)arctan(2z0/b)1+(2z0/b)2q2
(47)
ϕ1=π2(q2)arctan(2z0/b).
(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(1+(2z/b)2), which diverges. The ratio of the imaginary to real part is, therefore, infinite and arctan(x)=π/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 ϕ1=π/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π, which is shown in Fig. 7(a), marks a phase difference of the wavelets of 2π relative to the first dip, which corresponds to a total phase of π, i.e., the marked minimum-to-minimum distance corresponds to l~2π=l3πlπ. It is in general larger than the preceding maximum-to-maximum distance, i.e., l~2π>l2π; however, l~2π can be directly compared with the equivalent minimum-to-minimum values of the experimental work26 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(x) determines the phase of the air signal; equally, the translation of the focal position zzzfoc also shifts the integral boundaries of the LN signal. Approximating again with plane waves, one obtains Eq. (51),

I=|A×eiϕ1+B×0LeiΔk(zzfoc)dz|2,withA,BR+
(49)
=|Aeiϕ1+B×eiΔkL/2iΔkzfoc2Δksin(ΔkL/2)|2=|A|2+|B2Δksin(ΔkL/2)|2
(50)
+AB4Δkcos(ΔkL/2Δkzfocϕ1)sin(ΔkL/2).
(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.,

ϕdepth(zfoc)=Δkzfoc+ϕ1(b,zfoc)
(52)

allows for approximating the period as ldepth,2π2πΔk+O(1b). The similarity of both periods is visible in Fig. 7(b).

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 works9,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 measurements11 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π 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π 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π in this case, on the NA, which is not implemented in the plane wave case.

As a note, the full 3D-numerics11 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 nLN2.2, the effective NA is reduced, so that sin[arcsin(0.8/2.2)]0.360.8/2.2 meaning sin(x)x is valid and cos[arcsin(0.8/2.2)]0.93, so cos(x)1 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 NAeffNAn0.36, the approximation is applicable.

FIG. 9.

Second-harmonic signal data of line-scans, cf. Fig. 7(c), for different wavelengths using a NA of 0.8 together with the simulated thickness evolution using the full11 and the paraxial model.

FIG. 9.

Second-harmonic signal data of line-scans, cf. Fig. 7(c), for different wavelengths using a NA of 0.8 together with the simulated thickness evolution using the full11 and the paraxial model.

Close modal

One additional observation can be inferred from the properties of the experimental light source, which shows a bandwidth of 10 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.

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 (Δk1/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 (1μm) deviate from the surface.

For thin film materials for SHG and THG, we observe good agreement with rigorous 3D numerics18 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.

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.

The authors have no conflicts to disclose.

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).

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

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 tb, 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 tb, 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,

S(L,zf=0;Δk)χeff×tAir,LN(ω)tLN,Air(ωq)×0L(1+i2z/b)×eiΔkz[n=0Nfcpos,neik2nL1+i2[z+2nL]/b]q×j=0Nhcpos,j(ωq)eikq2jL1+i2[+2jL]/b+0L(1+i2(z+L)/b)×eiΔk(z+L)[n=0Nfcneg,neik2nL1+i2[(2n+1)L+z]/b]q×j=0Nhcneg,j(ωq)eikq2jL1+i2[+2jL]/b.
(A1)

Here, χ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(ω)rLN,Air(ω)]n and cneg,n=rLN,Si(ω)[rLN,Si(ω)rLN,Air(ω)]n, as well as cpos,k(ωq)=rLN,Si(ωq)[rLN,Air(ωq)rLN,Si(ωq)]k and cneg,k(ωq)=[rLN,Air(ωq)rLN,Si(ω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 Lb, 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×ei2kL 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.

1.
S.
Yuan
,
C.
Hu
,
A.
Pan
,
Y.
Ding
,
X.
Wang
,
Z.
Qu
,
J.
Wei
,
Y.
Liu
,
C.
Zeng
, and
J.
Xia
, “
Photonic devices based on thin-film lithium niobate on insulator
,”
J. Semicond.
42
,
041304
(
2021
).
2.
M.
Allgaier
,
V.
Ansari
,
L.
Sansoni
,
C.
Eigner
,
V.
Quiring
,
R.
Ricken
,
G.
Harder
,
B.
Brecht
, and
C.
Silberhorn
, “
Highly efficient frequency conversion with bandwidth compression of quantum light
,”
Nat. Commun.
8
,
14288
(
2017
).
3.
D.
Bonneau
,
M.
Lobino
,
P.
Jiang
,
C. M.
Natarajan
,
M. G.
Tanner
,
R. H.
Hadfield
,
S. N.
Dorenbos
,
V.
Zwiller
,
M. G.
Thompson
, and
J. L.
O’Brien
, “
Fast path and polarization manipulation of telecom wavelength single photons in lithium niobate waveguide devices
,”
Phys. Rev. Lett.
108
,
053601
(
2012
).
4.
X.
Wang
,
P. O.
Weigel
,
J.
Zhao
,
M.
Ruesing
, and
S.
Mookherjea
, “
Achieving beyond-100-GHz large-signal modulation bandwidth in hybrid silicon photonics Mach Zehnder modulators using thin film lithium niobate
,”
APL Photonics
4
,
096101
(
2019
).
5.
M.
He
,
M.
Xu
,
Y.
Ren
,
J.
Jian
,
Z.
Ruan
,
Y.
Xu
,
S.
Gao
,
S.
Sun
,
X.
Wen
,
L.
Zhou
,
L.
Liu
,
C.
Guo
,
H.
Chen
,
S.
Yu
,
L.
Liu
, and
X.
Cai
, “
High-performance hybrid silicon and lithium niobate Mach-Zehnder modulators for 100 Gbit rs1 and beyond
,”
Nat. Photonics
13
,
359
364
(
2019
).
6.
G.
Chen
,
K.
Chen
,
R.
Gan
,
Z.
Ruan
,
Z.
Wang
,
P.
Huang
,
C.
Lu
,
A. P. T.
Lau
,
D.
Dai
,
C.
Guo
, and
L.
Liu
, “
High performance thin-film lithium niobate modulator on a silicon substrate using periodic capacitively loaded traveling-wave electrode
,”
APL Photonics
7
,
026103
(
2022
).
7.
P. V.
Tien
, “
Light waves in thin films and integrated optics
,”
Appl. Opt.
10
,
2395
2413
(
1971
).
8.
H.
Pulker
, “
Characterization of optical thin films
,”
Appl. Opt.
18
,
1969
1977
(
1979
).
9.
K. J.
Spychala
,
P.
Mackwitz
,
M.
Rüsing
,
A.
Widhalm
,
G.
Berth
,
C.
Silberhorn
, and
A.
Zrenner
, “
Nonlinear focal mapping of ferroelectric domain walls in LiNbO3: Analysis of the SHG microscopy contrast mechanism
,”
J. Appl. Phys.
128
,
234102
(
2020
).
10.
S.
Cherifi-Hertel
,
H.
Bulou
,
R.
Hertel
et al., “
Non-ising and chiral ferroelectric domain walls revealed by nonlinear optical microscopy
,”
Nat. Commun.
8
,
15768
(
2017
).
11.
Z. H.
Amber
,
B.
Kirbus
,
L. M.
Eng
, and
M.
Rüsing
, “
Quantifying the coherent interaction length of second-harmonic microscopy in lithium niobate confined nanostructures
,”
J. Appl. Phys.
130
,
133102
(
2021
).
12.
J. A.
Squier
,
M.
Müller
,
G.
Brakenhoff
, and
K. R.
Wilson
, “
Third harmonic generation microscopy
,”
Opt. Express
3
,
315
324
(
1998
).
13.
Y.
Shen
, “
Surface properties probed by second-harmonic and sum-frequency generation
,”
Nature
337
,
519
525
(
1989
).
14.
L.
Zhou
,
H.
Fu
,
F.
Lv
,
C.
Wang
,
H.
Gao
,
D.
Li
,
L.
Deng
, and
W.
Xiong
, “
Nonlinear optical characterization of 2D materials
,”
Nanomaterials
10
,
2263
(
2020
).
15.
P. A.
Hegarty
,
H.
Beccard
,
L. M.
Eng
, and
M.
Rüsing
, “
Turn all the lights off: Bright- and dark-field second-harmonic microscopy to select contrast mechanisms for ferroelectric domain walls
,”
J. Appl. Phys.
131
,
244102
(
2022
).
16.
S.
Reitzig
,
F.
Hempel
,
J.
Ratzenberger
,
P. A.
Hegarty
,
Z. H.
Amber
,
R.
Buschbeck
,
M.
Rüsing
, and
L. M.
Eng
, “
High-speed hyperspectral imaging of ferroelectric domain walls using broadband coherent anti-Stokes Raman scattering
,”
Appl. Phys. Lett.
120
,
162901
(
2022
).
17.
F.
Hempel
,
S.
Reitzig
,
M.
Rüsing
, and
L. M.
Eng
, “
Broadband coherent anti-Stokes Raman scattering for crystalline materials
,”
Phys. Rev. B
104
,
224308
(
2021
).
18.
D.
Sandkuijl
,
A. E.
Tuer
,
D.
Tokarz
,
J. E.
Sipe
, and
V.
Barzda
, “
Numerical second- and third-harmonic generation microscopy
,”
J. Opt. Soc. Am. B
30
,
382
(
2013
).
19.
R. W.
Boyd
,
Nonlinear Optics
(
Academic Press, Inc.
,
San Diego
,
1992
).
20.
L.
Novotny
and
B.
Hecht
,
Principles of Nano-Optics
(
Cambridge University Press
,
2012
).
21.
K. J.
Spychala
,
P.
Mackwitz
,
A.
Widhalm
,
G.
Berth
, and
A.
Zrenner
, “
Spatially resolved light field analysis of the second-harmonic signal of χ(2)-materials in the tight focusing regime
,”
J. Appl. Phys.
127
,
023103
(
2020
).
22.
S.
Hell
,
G.
Reiner
,
C.
Cremer
, and
E.
Stelzer
, “
Aberrations in confocal fluorescence microscopy induced by mismatches in refractive index
,”
J. Microsc.
169
,
391
405
(
1993
).
23.
M. J.
Nasse
and
J. C.
Woehl
, “
Realistic modeling of the illumination point spread function in confocal scanning optical microscopy
,”
J. Opt. Soc. Am. A
27
,
295
302
(
2010
).
24.
G. D.
Boyd
and
D. A.
Kleinman
, “
Parametric interaction of focused Gaussian light beams
,”
J. Appl. Phys.
39
,
3597
3639
(
1968
).
25.
S.
Saravi
,
T.
Pertsch
, and
F.
Setzpfandt
, “
Lithium niobate on insulator: An emerging platform for integrated quantum photonics
,”
Adv. Opt. Mater.
9
,
2100789
(
2021
).
26.
Z. H.
Amber
,
K. J.
Spychala
,
M.
Eng
, and
M.
Rüsing
, “
Nonlinear optical interactions in focused beams and nanosized structures
,”
J. Appl. Phys.
132
,
213102
(
2022
).
27.
D. E.
Zelmon
,
D. L.
Small
, and
D.
Jundt
, “
Infrared corrected Sellmeier coefficients for congruently grown lithium niobate and 5 mol.% magnesium oxide-doped lithium niobate
,”
J. Opt. Soc. Am. B
14
,
3319
3322
(
1997
).
28.
K. J.
Spychala
,
G.
Berth
,
A.
Widhalm
,
M.
Rüsing
,
L.
Wang
,
S.
Sanna
, and
A.
Zrenner
, “
Impact of carbon-ion implantation on the nonlinear optical susceptibility of LiNbO3
,”
Opt. Express
25
,
21444
21453
(
2017
).