Giant fast-response nonlinearities of the refractive index of a vibrational nature have recently been theoretically predicted and experimentally observed in a variety of media in the terahertz (THz) spectral range^[1–7]. This resulted in a particular relevance of the systemic analysis of giant nonlinearity mechanism, the study of the physical processes and wave phenomena in the media with such nonlinearity, and the potential applications of these nonlinear effects of high-speed THz photonics.

Numerous interesting results have already been obtained in the field of nonlinear optics of pulsed THz radiation. Its main features are related to the fact that THz pulses are typically few-cycle ones, as has been observed in many experimental works^[8–11]. An example would be the phenomenon of pulse self-focusing disappearing with just one complete cycle, even at power levels significantly exceeding the critical one^[12]. This is due to the fact that with a few-cycle pulses, the dispersion may prevail over diffraction, and the concept of the critical power of self-focusing loses its physical meaning. However, the strong nonlinearity of the polarization response of materials in the field of THz pulses can lead to the formation of noticeable nonlinear lenses even in very thin layers of the materials with the thicknesses of less than one wavelength of the radiation^[3–5,13]. Such lenses are adequate for measuring nonlinear refractive index coefficients of the media.

Moreover, it appears that for single-cycle THz pulses propagating in a quadratic nonlinear medium, there is no effect of radiation generation at doubled frequencies, typical for quasi-monochromatic pulses in conventional nonlinear optics. Instead, radiation generation occurs at triple frequencies^[14]. A similar result was rather unexpectedly predicted in a theoretical work for cubic nonlinear media^[15]: in the field of a single-cycle wave, radiation generation disappears at triple frequencies and emerges at quadrupled frequencies. This effect was experimentally confirmed in Ref. [16]. In Ref. [16], the authors explained the above phenomenon by the interference of the waves generated in the nonlinear medium at triple frequencies and induced by self-phase modulation. Such interference is unconventional for the pulses in the visible and near-infrared spectral ranges. Even in the generation of spectral supercontinuum in the field of intense femtosecond pulses, their spectrum typically does not extend to the third harmonic of the femtosecond laser radiation. This phenomenon is characteristic of few-cycle pulses.

This study explores the patterns of interference phenomena in the waves generated at triple frequencies and in those caused by self-phase modulation during the propagation of pulsed THz radiation in media with cubic nonlinearity. The investigation considers the influence of the number of oscillations and the temporal shape of the input pulse on the interference. The overlap coefficients are found to depend on the spectra corresponding to self-phase modulation and radiation generation at triple frequencies. The study also examines the attenuation coefficients of the above effects based on the temporal shape of few-cycle pulses, considering both the even and the odd symmetry, including unipolar pulses. It is demonstrated that the interference considered can be either constructive or destructive, and depending on the pulse profile, the effects are either enhanced or attenuated. The theoretical simulation is compared with the experimental results in Refs. [16,17], and it is demonstrated that for both cases, the minimal values are located at triple frequencies and the spectral shapes are in good agreement. The obtained results are relevant for addressing challenges in controlling THz radiation and developing devices operating in the THz spectral range^[6,18–20].

The model of pulsed THz radiation propagation in a transparent dielectric cubic nonlinear medium is based on the following field equation^[21]: ∂E∂z′+1Lp∂E∂t−1Ldisp∂3E∂t3+1LnlE2∂E∂t=1LdiffrΔ⊥∫−∞tEdt′,(1)where z′ is the spatial coordinate along which the radiation propagates, E=E′/E0 is the normalized field of the THz pulse, E′ is the field of the THz pulse, which at the input of the nonlinear medium (at z=0) can be presented in a conventional optical manner^[22], E′(t′,0)=E0sin(ω0t′+φ)exp[−(t′τp)2],(2)where E0 is the amplitude of the THz pulse, t=ω0t′ is the normalized time, t′ is the time, ω0 is the angular frequency, Lp=λ0/(2π) is the characteristic longitudinal pulsed radiation length on which the field changes notably, λ0=2πcn0ω0 is the central wavelength, c is the speed of light in vacuum space, Ldisp=Lpn0/Δndisp is the dispersion length, Δndisp=acω02 is the refractive index variation due to dispersion, n0 and a are the empirical constants describing the dependence of the linear refractive index n on the frequency ω in the following form: n(ω)=n0+caω2, Lnl=Lpn0/(4Δnnl) is the nonlinear length, Δnnl=n2′E02/2=n2I0 is the radiation-induced change in the refractive index of the medium, n2′ and n2 are the coefficients of the nonlinear refractive index of the medium in CGS and SI, I0=cn0ε0E02/2 is the intensity of the THz pulse, Ldiffr=2r02/Lp is the diffraction length, r0 is the beam transverse radius, Δ⊥=∂∂x2+∂∂y2 is the transverse Laplacian, x=x′/r0, y=y′/r0, and x′,y′ are the transverse spatial coordinates.

In the pulse form described in Eq. (2), τp is the pulse duration. Varying the parameter τp changes the number of oscillations within the pulse. Depending on the phase value φ=0 or π/2, Eq. (2) describes either a sinelike or cosinelike function^[10,18,23]. The carrier-envelope phase (CEP) of ultrashort pulses, corresponding to the phase φ in our model, is known to have a significant impact on the nature and efficiency of both linear and nonlinear processes and on the pulse shape^[24,25]. The processes involved include the second- and the third-harmonic generation, pulse synchronization, and other linear and nonlinear processes. Due to the extremely small number of full oscillations in the THz pulse, the concept of an envelope loses its physical implication. However, the term CEP to describe the phase of pulses with extremely low oscillation numbers remains adequate. It should be noted that based on the above, CEP will make a significant contribution to the nonlinear phenomena in this frequency range, which will be discussed below.

In this study, we consider a situation where nonlinearity contributes more significantly to the evolution of the pulse than dispersion and diffraction, i.e., Lnl≪Ldisp,Ldiffr. In this case, the following conditions must be met: Δnnl≫Δndisp/4 and Δnnl/n0≫0.5(Lp/D0)2; D0=2r0 is the transverse beam diameter.

Such a situation is realistic in the THz frequency range. For instance, at an intensity of I0=108W/cm2 in a LiNbO3 crystal for ω0=1THz, r0=10λ0=3.7mm, where n0=5.12, a=4.47·10−36s3/m^[26,27] and n2=5·10−11cm2/W, the estimate holds: Δnnl=2.5·10−3, Δndisp/4=3.3·10−4, and 0.5(Lp/D0)2=3.2·10−5.

Therefore, there are media in which the influence of nonlinear effects becomes more significant than diffraction, making this effect negligible^[28]. In such cases, the dynamics equation, Eq. (1), can be represented as ∂E∂z+∂E∂t+μnl(E2∂E∂t−μdispμnl∂3E∂t3)=0,(3)where z=L0/z′, μnl=Lp/Lnl=4Δnnl/n0, and μdisp=Lp/Ldisp=Δndisp/n0. We mentioned before that influences of dispersion were smaller than those of nonlinearity, but the dispersion can also be significant for the frequency range where the intensity became small.

Thus, for the given example, μnl∼10−3, the above equation can be solved through series expansion based on this small parameter, E(t,z)=E(0)(t,z)+μnl[E(1,nl)(t,z)+μdispμnlE1,disp(t,z)]+…,(4)where E(i) represents the ith term in the field expansion.

The solution to Eq. (3) in Eq. (4) is given by^[28]E(0)(τ,z)=E(0)(τ),(5a)E(1)(τ,z)=−[E(τ)(0)]2∂E(0)(τ)∂τz,(5b)where τ=t−z is the reduced time.

Here we have to mention that the given models of pulse propagation through the nonlinear media must follow Rosanov’s rule about the conservation of the electric pulse’s square^[29–31]. The given Eq. (3) follows this rule.

The normalized pulse field expression at the input of Eq. (2) takes the following form: E(τ,0)=sin(τ+φ)exp[−(τN)2],(6)where N=ω0τp is the number of full oscillations in the pulse^[14]. The increase of N value causes the increase in the oscillation number, which follows the growth of τp value. The spectrum of described pulse takes the form G˜(ω,z=0)=iπN2·[e−iφe−14N2(ω+1)2−eiφe−14N2(ω−1)2].(7)

In Fig. 1, pulse shapes are presented for φ=0 [sine-like function, (a)] and for φ=π/2 [cosine-like function, (b)] at various values of N. The normalized magnitude spectra of the pulses described at the input of the medium for sine-like and cosine-like functions, respectively, are illustrated in Figs. 1(c) and 1(d). For clarity, the illustrations correspond to the central frequency of THz radiation, ω0=1THz.

It should be noted that the field representation given by Eq. (6) at φ=0 and N=1 corresponds, with 97% accuracy (in terms of root-mean-square deviation), to another commonly used representation of a single-cycle pulse as follows^[15,16]: E(τ,0)=τexp(−τ2).(8)

When φ=π/2, the magnitude of the spectrum at zero frequency is different from 0. This is explained by the fact that such pulses are video pulses (or unipolar pulses for small N) or similar to them^[32–35] (quasi-video pulses or quasi-unipolar pulses) for large values of N. The value of the direct current (DC) component at zero frequency is πNexp(−N2/4). An increase in the value of N will evidently lead to a rapid decrease in the DC.

The solution to Eq. (5b) for the input pulse shape in Eq. (6) is given by E(1)(τ,z)=z{3sin(τ+ϕ)−sin[3(τ+ϕ)]}4N2·[2τ−N2cot(τ+ϕ)]exp(−3τ2N2).(9)

Then the spectra of contributions describing the generation of radiation at triple frequencies (TFG) and of the radiation obtained due to the self-phase modulation (SPM) of the initial pulse can be expressed as $$\begin{gathered} Gnl(1)(ω,z)SPM=−zNπω43exp[−N2(ω2+1)12] \\ ⁢[isinφcosh(N2ω6)+cosφsinh(N2ω6)]⁢exp(−iωz), \end{gathered}$$(10)$$\begin{gathered} Gnl(1)(ω,z)TFG=zNπω123exp[−N2(ω2+9)12] \\ ⁢[isin(3φ)cosh(N2ω2)+cos(3φ)sinh(N2ω2)]⁢exp(−iωz). \end{gathered}$$(11)

The spectra of SPM and TFG are shown in Fig. 2. Their different signs for sine-like pulses at the input of the medium [Figs. 2(a) and 2(c)] are due to the fact that they differ in phase by π [Eqs. (10,11)]. In the case of cosine-like pulses [Figs. 2(b) and 2(d)], the SPM and TFG spectra do not have a mutual shift by π, but they are completely imaginary. For clarity, the illustrations correspond to the central frequency of THz radiation, ω0=1THz.

Figure 3 demonstrates the interaction of nonlinear pulse responses. The left column demonstrates the pulses with φ=0, the right is for φ=π/2, the upper row is for N=1, and the lower row is for N=3. It shows that for the sine-like form SPM and TFG destructively interfere, while for cosine-like pulses these components interfere constructively.

It is also worth noting that for both sine-like and cosine-like pulses for small values of N [Figs. 2(a) and 2(b)], there is a significant overlap between the SPM and TFG spectra, noticeably altering the overall spectral contribution. However, for large values of N [Figs. 2(c) and 2(d)], the overlapping region is significantly reduced. To numerically estimate the spectrum overlap, we introduce the overlap coefficient using the expression, s=∫−∞∞Gnl(1)(ω,z)SPMGnl(1)*(ω,z)TFGdω∫−∞∞|Gnl(1)(ω,z)SPM|2dω∫−∞∞|Gnl(1)(ω,z)TFG|2dω,(12)where * means the complex conjugation.

The coefficient introduced in this way is a real quantity in the case considered. This coefficient is similar to the first-order correlation function. The numerator describes the energy of mutual radiation multiplication where the components in the denominator describe the energies of both the SPM and TFG components. A negative value of the coefficient indicates that the TFG and SPM spectra differ in phase by π. Figure 4(a) shows that the overlap coefficients s for φ=0 at N=1, N=3, and N=5 are −0.95, −0.28, and −0.02, respectively, and for φ=π/2 at N=1, N=3, and N=5 are 0.88, 0.27, and 0.02, respectively. Figure 4(b) indicates that for N=1, the spectrum overlap changes from −1 to +1 as φ varies from 0 to π/2, passing through 0 at φ≈0.93rad.

Since for small N values, the spectra of SPM and TFG tend to overlap, it is natural to refer to the interference of the radiation generated at triple frequencies and arising due to SPM. To assess whether these effects attenuate or enhance each other, we will evaluate the coefficient of mutual influence of SPM and TFG, k(N,φ)=∫−∞−∞|Gnl(1)(ω,z)TFG+Gnl(1)(ω,z)SPM|2dω∫−∞−∞|Gnl(1)(ω,z)TFG|2dω+∫−∞−∞|Gnl(1)(ω,z)SPM|2dω.(13)

The physical meaning of Eq. (13) describes the relation of interfered energies of SPM and TFG components to the energies of the sum of separate parts. It demonstrates to us how the interaction of SPM and TFG fields increases or decreases the nonlinear response due to pulse propagation through the medium. With this definition, if the coefficient k is less than 1, the interference of TFG and SPM leads to their mutual attenuation; if it is more than 1, then they interfere constructively, amplifying each other.

In Fig. 5, the dependence of the mutual influence coefficient on the parameter φ is presented for fixed values of N. It can be seen that for N=1, the coefficient k undergoes the most significant changes at the extreme values of the parameter φ: at φ=0, the coefficient is approximately 0.05, while at φ=π/2, its value is 1.7. The value of 0.05 indicates that the effects of SPM and TFG mutually suppress each other significantly. There is almost no transfer of the energy of the main pulse to high-frequency radiation. In the case of φ=π/2, there is some increase in the energy of the radiation caused by SPM and TFG.

Thus, for N=1, CEP has a significant impact on the value of the coefficient k: changing φ from 0 to π/2 results in a 32-fold change in this coefficient. For N=3, the change is 1.8, and for N=5 it equals 1. In other words, the interference of radiations caused by different nonlinear effects is significant for the pulses with the small number of full oscillations, and it is particularly manifested for single-cycle pulses.

To approve the results of theoretical simulations, we present its comparison with the experimental data in Refs. [16,17] in Fig. 6. References [16,17] described THz pulse generation via optical rectification of intense femtosecond radiation in nonlinear crystal of LiNbO3. It is well known that phase quasi-synchronism and therefore THz radiation generation occurs at a distance of the pulse radius^[36], so the THz wave propagates in the remaining part of the crystal and is exposed to nonlinear effects. It was experimentally confirmed in Ref. [16]. To obtain pulse forms in the near-field region, the integration of the field data is completed with E(t,z)=∫−∞tE(t′,z)dt′^[37]. To simulate theoretical results, we used experimental parameters for LiNbO3 and generated THz pulse radiation; for Ref. [17], we set the following parameters: the central frequency of THz radiation is νc=1.05THz, the pulse duration τp is 0.4 ps, pulse intensity is I=107W/cm2, nonlinear refractive index coefficient n2∼10−11cm2/W, and sample thickness is 5 mm. For Ref. [16], the central frequency of THz radiation equals νc=0.55THz, the pulse duration τp is 0.55 ps, pulse intensity is I=108W/cm2, nonlinear refractive index coefficient n2∼10−11cm2/W, and the sample thickness is 5 mm.

Figures 6(a) and 6(c) represent the experimental data from Ref. [17] for φ=11π6. Figures 6(b) and 6(d) demonstrate the experimental data in Ref. [16] for φ=11π36. The following theoretical spectrum curves are obtained with Eqs. (7), (10), and (11). The solid blue curves describe experimental data, and the red dashed curves are for theoretical calculation. As we can see, the theoretical minimum value located on the triple frequencies according to the spectrum maximum is in good agreement with experimental measurements. Similar experimental results can be found in other works as well, for instance, in Ref. [38]. This proves the validity of theoretical model for the interference processes.

This work presents a theoretical analysis of the evolution of the field and the spectra of Gaussian-shaped THz pulses with a small oscillation number under the envelope as they propagate in a medium with cubic nonlinearity.

It has been shown that the changes in the pulse spectrum due to the effects of TFG and SPM overlap significantly when the number of oscillations in the pulse is less than three, and this overlap depends on CEP. For example, the overlap coefficients s for CEP φ=0 with the parameters N=5, N=3, and N=1 are −0.02, −0.28, and −0.95, respectively, and for φ=π/2 with N=5, N=3, and N=1, they are 0.02, 0.27, and 0.88, respectively.

With the mentioned spectrum overlap, the interference between the radiation generated at triple frequencies and that arising due to SPM becomes crucial. Whether these effects attenuate or enhance each other is determined by the CEP. For example, with N=1, the coefficient of mutual interference attenuation of these effects at φ=0 is 0.05, while at φ=π/2, the effects mutually enhance each other by a factor of 1.7. Changing φ from 0 to π/2 results in a 32-fold change in the mutual influence of these effects. For N=3, this change is only 1.8 times. In other words, the interference of radiations caused by different nonlinear effects is significant only for the pulses with the number of complete oscillations of three or less, and it is particularly manifested for single-cycle pulses.

The theoretical simulation was approved with experimental data presented in Refs. [16,17]. The obtained results demonstrate the good agreement of spectral minimum values between experimental data and simulation results.

We assume that the presented results are relevant for addressing the challenges related to the control of pulsed THz radiation and the development of devices operating in the THz spectral range.