High Efficiency Third-Harmonic Generation in a Medium with Quadratic Susceptibility Due to Cubic-like Nonlinearity Caused by Cascaded Second-Harmonic Generation

Abstract

Third-harmonic generation (THG) is of interest for its various applications. We propose using the cascaded second-harmonic generation (SHG) to implement the frequency conversion process, which is similar to that occurring in a medium with cubic susceptibility. Physically, the process is based on consecutive generation of the second-harmonic and the sum frequency in the same crystal with quadratic susceptibility at large phase mismatching between the fundamental wave and the second-harmonic wave. In this case, at phase matching between the fundamental wave and the third-harmonic wave, THG occurs with high efficiency. To demonstrate such a possibility theoretically, we apply the multi-scale method to a set of Schrödinger equations, describing a three-wave interaction with the frequencies $\omega ,2\omega$ and $3\omega$ in a medium with quadratic susceptibility, to derive modified equations describing the frequency tripling process. These equations are solved without using the fundamental wave energy non-depletion approximation. A THG efficiency equal to $94.5\%$ is predicted theoretically. The analytical solution is confirmed by computer simulation results. We study how various factors, such as the incident pulse intensity, phase mismatching between interacting waves, group velocity mismatching of the pulses, and second-order dispersion of the wave packets influence the THG process.

1. Introduction

Third-harmonic generation (THG) attracts permanent attention from researchers because it is one of the effective ways to obtain laser radiation with a frequency in the ultraviolet range. Recently, the frequency tripling process was applied in the National Ignition Facility laser system, in which the first ever controlled laser fusion with a positive output energy balance was achieved [1]. There are well-known schemes for reaching the high-efficiency THG. Among them, we first emphasize a scheme using two crystals with quadratic susceptibility ${\chi }^{\left(2\right)}$. The first crystal is used for phase-matched second-harmonic generation (SHG), and the second one is used for phase-matched sum frequency wave generation (SFG). Let us note that such a scheme of THG is often called the cascaded generation of the third harmonic (TH). Using this scheme, a conversion efficiency of $80\%$ has been reached in physical experiments [2,3]. With this aim, two Type II potassium dihydrogen phosphate (KDP) crystals were used, each 70 mm in diameter and 12 mm thick, polished to a flatness of $\lambda /4$ at 0.5 μm, and with a wedge of 7 arcsec. It should be noted that the two-crystal scheme required a certain ratio between the fundamental wave (FW) and the second-harmonic (SH) intensities in the input of the second crystal which is equal to 1:2. Another obstacle is the influence of group velocity mismatching between FW and SH, which can strongly decrease the frequency conversion efficiency of a short pulse. A possible solution was proposed in [4], where a doubled frequency was generated using a pulse with a duration of 1 ps at the fundamental frequency (FF). Then, the SH pulse was compressed until a duration of 200 fs (its power density reached 40 GW/cm${}^2$) and was mixed with the pulse at the FF. The frequency conversion efficiency reached about 40%. Another experimental setup was proposed in [5], where a KDP crystal was used to delay the extraordinary FW with a pulse duration of 0.6–1 ps and the incident power density was varied in the 6–75 GW/cm${}^2$ range.

THG has also been widely studied in a single crystal with quadratic susceptibility. A review of different approaches can be found in [6]. Nevertheless, many papers devoted to this problem have been published during passed years, and we briefly discuss them here.

One of the most popular ways to carry out frequency tripling for optical radiation is via frequency conversion using the quasi-phase matching (QPM) technique, in which both the SHG and SFG processes are quasi-phase-matched (the description of this technique can be found, for example, in [7]). In such a structure, various materials are used to reach a high efficiency in frequency conversion. It should be emphasized that lithium niobate (LiNbO${}_3$) is the most popular crystal. So, a 2D periodically poled LiNbO${}_3$ crystal was used for frequency tripling of a 118 fs pulse at a 1568 nm wavelength with a maximum power density of 3.2 GW/cm${}^2$ [8]. An $8\%$ conversion efficiency was reached. The same frequency conversion efficiency was also achieved in a 16 mm chirped periodically poled LiNbO${}_3$ crystal [9] for a 3.5 ns pulse with average power changing between 5 and 20 mJ. THG in a 12 mm short-range-ordered crystal was also obtained in [10] for a 16 ps pulse with 1.73 μJ energy and 12% frequency conversion efficiency.

Another popular nonlinear crystal used for QPM is the lithium tantalate (LiTaO₃) crystal. A frequency conversion efficiency of about 27% was reached in a 15 mm quasi-periodic optical super-lattice for an 1444 nm wavelength and with an average power of 20 mW [11]. An 18 mm periodically poled LiTaO₃ crystal was used for frequency tripling of a 30 ns pulse with a 1342 nm wavelength to reach the frequency conversion efficiency of about 19%. Applying the QPM realized using the electro-optic effect in a quasi-periodic optical super-lattice LiTaO₃ resulted in an 18% frequency conversion efficiency for a 90 ns pulse with an average power of 1.72 W [12]. To date, the highest conversion efficiency (40%) for this approach was obtained in a 10 mm quasi-periodically poled potassium titanyl phosphate (KTiOPO₄ or KTP) crystal for a 15 ps pulse with a 1.5 GW/cm² power density [13]. A similar crystal was also used for frequency tripling of a femtosecond (200 fs) pulse with an average power of 2 W [14]. A conversion efficiency of about 8% was reached. This frequency conversion tool was also studied theoretically. The optimal ratio between the nonlinear coupling coefficients and refraction indexes that defines SHG and SFG, allowing to reach 100% conversion efficiency, was discussed theoretically in [15]. This ratio was further used in the numerical simulation of frequency conversion in a bulk medium for a beam with a waist radius of 15 μm propagating in a 2 cm super-lattice [16], and a frequency conversion efficiency of about 96% was achieved. The QPM configurations of the various types of waves’ interactions were studied in [17]. The disadvantages of this approach are more difficult crystal fabrication and higher cost in comparison with ordinary nonlinear crystals [18]. Furthermore, QPM is harder to implement for short wavelengths, such as those of visible or ultraviolet light, and this technique gives a lower conversion efficiency [19]. Another disadvantage appears during ultraviolet light generation: crystal heating caused by optical energy absorption becomes essential if the light power exceeds a certain value (typically, several watts) [20,21]. Evidently, this restricts the frequency conversion efficiency.

One more approach for reaching high-efficiency THG (about 40%) in a single crystal was demonstrated in [18]. An FW pulse with a duration of 40 ps and a beam diameter 10 mm entered the KDP crystal (20 mm · 20 mm · 10 mm), in which the SH was generated. The FW energy varied from 0.5 to 4.5 mJ. Then, a quarter-wave plate was used to change the SH wave polarization: this was necessary to achieve phase matching at the SFG. After that, the FW and SH waves, passed through the crystal, were reflected by the mirror and propagated in the same crystal, in which THG occurred. A beam with a sufficiently big diameter was used in this experiment to avoid the influence of its diffraction. This approach was also implemented for ${\mathrm{Gd}}_x{\mathrm{Y}}_{1-x}\mathrm{COB}$ [22], ammonium-dihydrogen-phosphate (ADP) [23], and $\beta$-barium borate (BBO) [24] crystals with similar results.

In our opinion, there is another way to achieve high-efficiency THG in a single crystal with quadratic susceptibility. THG occurs under the phase-matched condition between FW and a wave at tripled frequency; meanwhile, both the SHG and SFG processes occur at large phase mismatching. To achieve this, a birefrigent crystal can be used. This means that the cascaded ${\chi }^{\left(2\right)}$ process (SHG at large phase mismatching) is used together with SFG (the second process that is named also as cascaded one as we mentioned above). Various aspects of such an interaction of waves and its influence on THG and the measurement of third-order susceptibility were previously discussed in the literature [25,26,27,28,29,30,31,32,33,34] (we discuss these papers below).

Let us note that the cascaded ${\chi }^{\left(2\right)}$ nonlinear phenomenon remains of interest for many years due to its various practical applications. This phenomenon was theoretically predicted in 1974 by A.P. Sukhorukov and Yu.N. Karamzin [35]. Beginning in 1991, there has been continued interest in this phenomenon. It allows us to obtain quadratic solitons: the ${\chi }^{\left(2\right)}$ temporal soliton [36] and the spatio-temporal soliton [37]. In [38], a tilted pulse was used to decrease the negative influence of group velocity mismatching between interacting waves during temporal soliton formation due to SHG in a quadratic medium. The appearance of a spatio-temporal soliton was experimentally demonstrated in [39].

A comprehensive review of cascading ${\chi }^{\left(2\right)}$ solitons and their practical applications can be found in [40]. Recently, a special type of solitons—Akhmediev breathers—in a medium with quadratic susceptibility accounting for the dispersion of a pulse or optical beam diffraction was studied in [41] and [42], respectively.

One of the cascading SHG applications is soliton compression. By choosing an appropriate sign of phase mismatching, pulse compression at the FF can be reached. Soliton compression at a wavelength of 1030 nm at its propagating in a LiNbO${}_3$ crystal with a length equal to 10 cm has been discussed in [43]. Few-cycle pulse’s formation at the wavelength of 1250 nm via soliton compression of the pulse with a duration of 100 fs and a power density of 30 GW/cm${}^2$ to the pulse with a duration of 6 fs using a defocusing quadratic medium was presented in [44]. In [45], such a pulse compression was observed in a BBO crystal with length of about 3 cm. Limits to pulse compression in cascaded quadratic soliton compressors were discussed theoretically in [46] using the nonlocal temporal response of a medium. Nonlinear phase shift [47], induced by quadratic cascaded processes, can also be used for compensation of self-phase modulation, caused by Kerr nonlinearity, as was shown numerically and experimentally in [48]. Spectral broadening, induced by intracavity self-phase modulation, in a Nd:YLF regenerative amplifier was overcome [49] with correct choice of phase mismatching between FW and SH waves.

Cascaded SHG may lead to self-focusing or defocusing of the beam at the FF, as was shown experimentally in [50], and then theoretically in [51]. In the axial-symmetric case, the beam intensity may increase about 70 times [52]. Another important application of the cascading SHG is the suppression of random fluctuations of both duration and maximum power density of the sub-pulses generated by laser operating in the so-called free-generation mode [53,54,55,56]. Only about 2–3% of the incident energy is lost when suppressing those fluctuations because of the frequency conversion to SH.

Thus, the phenomena inherent to the Kerr effect in a medium with quadratic susceptibility at the cascaded SHG (for brevity, we will name it as induced cubic nonlinearity or susceptibility) was observed in many papers. Herewith, those authors who studied various manifestations of the self-focusing effect accounted only for the phase self-modulation of the beam or pulse at the FF. As mentioned above, THG can be provided via the cascaded processes (SHG and then SFG) at the phase matching between FW and TH. However, in this case, there is a large phase-mismatching between FW and SH as well as between SH and TH. Note that this scheme of the frequency conversion refers to cascading processes similar to the cascading SHG discussed above.

One of the first papers accounting for a contribution of the cascaded THG at frequency tripling due to cubic susceptibility, was the paper [25], in which a measurement of the cubic susceptibility was made. This contribution was estimated in the framework of the energy non-depletion approximation for waves with FF and SH and without accounting for phase self- and cross-modulation of interacting waves. A few studies should be mentioned that refer to cascaded THG in the BBO crystal. A BBO crystal with a length of 7 mm was used for the frequency tripling of the pulse with a duration of 5 ps and a power density varying between 1 and 50 GW/cm${}^2$ [26]. Absolute value of phase mismatching between SH and FW varies between 187 and 600 mm${}^{-1}$. The reached frequency conversion efficiency was about 1%. For a longer pulse (45 ps) with a smaller incident power density (1.5 GW/cm${}^2$), the frequency conversion efficiency was only 0.007% [27]. Using a femtosecond pulse with a duration of 300 fs and a power density changing between 25 and 250 GW/cm${}^2$ allowed to reach the frequency conversion efficiency about 6% in a BBO crystal with a length of 1 mm, while phase mismatching equals −240 or 330 mm${}^{-1}$ [28,29]. Evidently, at such a pulse intensity, the cubic susceptibility also played a significant role in the frequency tripling. As shown in [30,31] for a KTP crystal with a length of 1.1 mm, the cubic susceptibility provides major contribution to the conversion efficiency, which achieved 2.9%, for a pulse of 22 ps with an incident power density of 45 GW/cm${}^2$ at the phase mismatching between SH and FW as low as −127 mm${}^{-1}$. An 11% conversion efficiency was reached in a 3 mm bismuth triborate (BiB₃O₆ or BiBO) crystal for the pulse with a duration of 80 fs and an incident power density of 700 GW/cm${}^2$ with an SHG phase mismatching of several hundred mm${}^{-1}$ [32].

From the obtained results [25,26,27,28,29,30,31,32], it follows that the frequency conversion efficiency is smaller than those reached using other setups. However, this situation can be changed. First of all, in all mentioned papers, the crystal length was too small (a few millimeters) for practical implementation of cascading processes occurring at large phase mismatching between FW and SH. Secondly, the theoretical studies made in those papers, as well as in [33], were provided in a framework of the pump-wave energy non-depletion approximation. Evidently, this approximation can be applied only at weak frequency conversion. Moreover, another condition often used in these papers was $|A_1|>>|A_2|>>|A_3|$, where $A_j,j=1,2,3$ are the complex amplitudes of FW, SH and TH, respectively. A main feature of this condition concludes in lower intensity of the TH in comparison with those for other waves. Although this condition was not used in [33], the analytical derived formulas were very complicated ones and their use is difficult. The numerical simulation was made in [34] for the crystal length equal to 1 mm that is not sufficient for observing frequency tripling due to the physical mechanism proposed above. Nevertheless, from our point of view, one of the main results of these theoretical studies was the demonstration of the appearance of $A_1^3$ in the solution, describing the amplitude of TH, together with ${\left({\chi }^{\left(2\right)}\right)}^2$ and phase mismatching ${\Delta }_{21}k$ between SH and FW upon using the pump-wave energy non-depletion approximation. However, other terms inherent to the cubic nonlinear response of a medium leading to self- and cross-modulation are not taken into account.

Thus, one can conclude that the authors of these papers did not account for all processes that occur at the cascaded THG under the condition of large phase-mismatching between SH and FW. Moreover, the used approach (the pump-wave energy non-depletion approximation) and length of the crystal did not allow us to perform deep analysis of the possibilities of the frequency conversion physical mechanism discussed below. The current study removes this gap because we present the study of THG in the crystal with quadratic susceptibility via the cascading SHG and SFG processes under the condition of the of large phase-mismatching between SH and FW without using the above mentioned restrictions. Based on multi-scale method (see, for example, [57]), we demonstrate that the cascaded SHG in a crystal with the quadratic susceptibility leads to a response of a medium similar to the cubic susceptibility allowing to observe THG. It should be emphasized that this method was applied in [51] in which a self-modulation of the pulse with FF was demonstrated under cascaded SHG. As is shown theoretically below, SHG efficiency is very low for the conditions under consideration—about 2–3%, or less than is comparable with losses using a mirror [2,3] without special features allowing to decrease its reflection on both wavelengths. The residual pump energy can be converted into the TH. Furthermore, in contrast to [3], the changing polarization of the electric field of the SH is not required because only the phase matching between TH and FW should be fulfilled. Based on developed analytical solution of the derived set of equations, we also show that to reach the maximum efficiency of the THG, the phase matching between TH and FW is not required, and weak phase mismatching between these waves is optimal.

Our paper is arranged as follows. In Section 2, we write the set of Schrödinger equations which describes three-wave interaction (FW, SH and TH) in a medium with quadratic susceptibility. The parameters of the crystals used for optical frequency tripling are also presented in this part of the paper. In Section 3, we derive the modified problem (Paragraph 3.1) using the multi-scale method presented in the Appendix A. Another, simpler method is presented in this part of the paper. In Paragraph 3.2, the analytical solution of the modified problem is derived. Here, the influence of phase mismatching between interacting waves on THG efficiency is examined theoretically in long pulse duration approximation without using the pump-wave energy non-depletion approximation. In Section 4, the computer simulation results are shown to support the developed analytical solution of the modified problem as well as the conclusions made on its base. This study is provided in Paragraph 4.1. In Section 4.2, a comparison between our approach and the fundamental pulse’s energy non-depletion approximation for description of the frequency conversion of short pulses is made.

In Section 5, the computer simulation results for an interaction of the shorter pulses are presented mainly for their propagation in a medium with the normal dispersion (Paragraph 5.1). Nevertheless, we also present some results for a medium with the anomalous dispersion at the FF (Paragraph 5.2). It is important to note that the highly effective THG is observed for a wide range of the problem parameters in both cases of the dispersion.

In Paragraph 5.3, Paragraph 5.5 and Paragraph 5.7, the influence of the phase matching between both SH and FW and TH and FW on the THG conversion efficiency is studied. In Paragraph 5.4, the developing of the modulation instability (MI) at negative phase mismatching between SH and FW is shown. It strongly decreases the frequency conversion efficiency. In Paragraph 5.6, it is shown that the nonzero phase mismatching between TH and FW can be preferable in some cases to reach higher THG efficiency. An influence of the incident FW intensity on the THG efficiency is also accounted for.

In Section 6, an influence of the group velocity mismatching (GVM) on the THG efficiency is examined. It is shown that the GVM between FW and SH does not remarkably influence the THG efficiency (Paragraph 6.1), despite the changing shapes of the FW and TH pulses. In contrast, a large GVM between FW and TH can dramatically decrease the conversion efficiency (Paragraph 6.2). Nevertheless, there are some possibilities which can lead to decreasing of the GVM action. One of them concludes in using of the pulses with a picosecond duration. Another one refers to introducing certain shifts of the FW carrier frequency (Paragraph 6.3).

In Section 7, we formulate the main results of the paper, discuss them, and compare them with the results achieved in various schemes of frequency tripling. Section 8 contains the paper’s conclusion.

2. Problem Statement

To clarify our scheme of the frequency tripling, let us consider the following frequency conversion processes, shown schematically in Figure 1 for a crystal with quadratic (${\chi }^{\left(2\right)}$) and cubic (${\chi }^{\left(3\right)}$) susceptibilities. A symbol $\omega$ denotes the FW frequency. Thus, we consider three processes producing a wave with tripled frequency.

The first process is a generation of both an optical wave with doubled frequency ($\omega +\omega =2\omega ,{\chi }^{\left(2\right)}$) and THG based on the SFG ($\omega +2\omega =3\omega ,{\chi }^{\left(2\right)}$) due to quadratic nonlinear response of a crystal. The second one is the directly tripled frequency wave generation ${\chi }^{\left(3\right)}$ due to cubic nonlinearity of a crystal. In comparison with the frequency doubling process, this process requires more intensive incident pulse. The third (new) process of frequency tripling is based on frequency doubling at large mismatching, and this process can occur at low intensities of the incident wave in comparison to those occurring at cubic nonlinearity of a crystal. This process is denoted as $3\omega =\omega +\omega +\omega$ with symbol ${\chi }_{cas}^{\left(3\right)}$. In Figure 1, we also show another frequency conversion process ($3\omega =2\omega +2\omega -\omega ,{\chi }^{\left(3\right)}$), which occurs in a cubic medium. Each of these frequency conversion processes is accompanied by own phase matching and own value of the crystal susceptibility at the corresponding frequencies. Therefore, the susceptibilities ${\chi }^{\left(3\right)}$ and ${\chi }_{cas}^{\left(3\right)}$ are different. As a rule, only one of various conversion processes can occur at phase matching: in Figure 1 this process is THG (${\chi }^{\left(3\right)}$). It should be emphasized again that the cascading susceptibility ${\chi }_{cas}^{\left(3\right)}$, including the ratio between ${\left({\chi }^{\left(2\right)}\right)}^2$ and phase mismatching between waves at the SHG, was demonstrated in several papers [25,26,27,28,29,30,32]. However, we emphasize again that these studies were made in the framework of the pump non-depletion approximation, and many factors essentially influencing the THG due to induced cubic effects were not taken into account.

We suppose that the angle range for achieving the phase matching of the corresponding generation process is small enough to avoid its overlapping with a phase matching angle range of other processes. In this case, it is possible to differ the process of THG, which occurs due to the SFG process (this process occurs at the corresponding phase matching) and due to the cascading SHG which occurs at phase matching between FW and TH and large phase mismatching for others processes of the frequency conversion.

We suppose also that the FW possesses an ordinary polarization and the TH has an extraordinary one. Therefore, the phase matching condition of the process $\omega \left(o\right)+\omega \left(o\right)+\omega \left(o\right)=3\omega \left(e\right)$ can be satisfied in nonlinear crystals discussed below. As for the SH wave, we believe that both cascaded processes $\omega \left(o\right)+\omega \left(o\right)=2\omega \left(o\right),\omega \left(o\right)+2\omega \left(o\right)=3\omega \left(e\right)$ and $\omega \left(o\right)+\omega \left(o\right)=2\omega \left(e\right),\omega \left(o\right)+2\omega \left(e\right)=3\omega \left(e\right)$ are possible, so we do not precise the polarization of the SH wave.

Starting from Maxwell’s equations, one can derive the following set of dimensionless Schrödinger equations with respect to slowly varying envelopes of the wave packets:

$$\begin{array}{c}\frac{\partial A_1}{\partial z}+i{D}_1\frac{{\partial }^2{A}_1}{\partial t^2}+i\left({\gamma }_{12}{A}_1^{\ast }{A}_2{e}^{-i{\Delta }_{21}kz}+{\gamma }_{23}{A}_2^{\ast }{A}_3{e}^{-i({\Delta }_{31}k-{\Delta }_{21}k)z}\right)=0,\hfill \\ \frac{\partial A_2}{\partial z}+{\nu }_{21}\frac{\partial A_2}{\partial t}+i{D}_2\frac{{\partial }^2{A}_2}{\partial t^2}+i\left({\gamma }_{11}{A}_1^2{e}^{i{\Delta }_{21}kz}+2{\gamma }_{13}{A}_1^{\ast }{A}_3{e}^{-i({\Delta }_{31}k-{\Delta }_{21}k)z}\right)=0,\hfill \\ \frac{\partial A_3}{\partial z}+{\nu }_{31}\frac{\partial A_3}{\partial t}+i{D}_3\frac{{\partial }^2{A}_3}{\partial t^2}+3i{\gamma }_{21}{A}_1{A}_2{e}^{i({\Delta }_{31}k-{\Delta }_{21}k)z}=0,0<z\le L_z,0<t<L_t\hfill \end{array}$$

(1)

which describes the optical frequency tripling in a medium with combined nonlinear response in the framework of the plane wave approximation (without taking into account a laser beam diffraction and absorption of the crystal).

Here, $A_1,A_2$ and $A_3$ are the complex amplitudes of waves with the FF ($\omega$), the doubled frequency ($2\omega$), and the tripled frequency ($3\omega$), respectively, normalized at $A_{01}=\sqrt{I_{01}}$ (a square root of the incident pulse maximum intensity at the FF). Parameter ${\gamma }_{jk}$ is the interacting waves’ coupling coefficient due to the quadratic nonlinear response at the corresponding frequencies. In general case, this coefficient may be different for various generation processes. However, for simplicity, we choose the same value (${\gamma }_{jk}=\gamma$) for all these processes involved in the set of equations (1) Parameter ${\Delta }_{21}k=k_2-2k_1$ characterizes the phase mismatching on the length $Z_n=1$ mm, in which units a longitudinal coordinate z is measured. $k_1$ and $k_2$ denote dimensionless wave numbers for waves with the FF and doubled frequency. Parameter ${\Delta }_{31}k=k_3-3k_1$ characterizes the phase mismatching on $Z_n$ at THG. The symbol $k_3$ denotes dimensionless TH wave number. Coordinate t is a time, measured in units of the pulse duration ${\tau }_P$ at the FF. Parameters ${\nu }_{21}$ and ${\nu }_{31}$ characterize a dimensionless GVM between SH and FW, and TH and FW, respectively. Parameters $D_j,j=1,2,3$ characterize dimensionless second-order dispersion (SOD) of the corresponding wave packet. Parameters $L_t$ and $L_z$ define the time interval, during which the laser pulse interaction and the maximum distance of the pulse propagation are analysed, respectively.

The initial distributions of the complex amplitudes are the following:

$$\begin{array}{cc}\hfill & A_1(0,t)=A_{10}\left(t\right),A_2(0,t)=0,A_3(0,t)=0,t\in [0,L_t],\hfill \end{array}$$

(2)

and corresponds to a finite distribution of the laser pulses.

Let us discuss the terms in the set of equations (1). The term $\gamma A_1^2{e}^{i{\Delta }_{21}kz}$ (${\chi }^{\left(2\right)}(2\omega ;\omega ,\omega$)) in the second equation corresponds to the generation of a wave with doubled frequency for the process ($\omega +\omega \to 2\omega$). The term $\gamma A_1^{\ast }{A}_2{e}^{-i{\Delta }_{21}kz}$ in the first equation is responsible for the opposite process ($2\omega \to \omega +\omega$). A process of the THG via the SFG is described by the term $\gamma A_1{A}_2{e}^{i({\Delta }_{31}k-{\Delta }_{21}k)z}$ (${\chi }^{\left(2\right)}(3\omega ;\omega ,2\omega$)) in the third equation. In turn, the term $\gamma A_2^{\ast }{A}_3{e}^{-i({\Delta }_{31}k-{\Delta }_{21}k)z}$ in the first equation and the term $\gamma A_1^{\ast }{A}_3{e}^{-i({\Delta }_{31}k-{\Delta }_{21}k)z}$ in the second equation correspond to the reverse processes ($3\omega \to 2\omega +\omega$).

Let us note that the introduced dimensionless parameters and functions are expressed through the physical ones, denoted by the containing line above them, with accordance to a rule:

$$\begin{array}{c}{A}_j=\frac{{\overline{A}}_j}{\sqrt{I_{01}}},D_j=-\frac{1}{2}\frac{{\partial }^2\overline{k}}{\partial {\overline{\omega }}^2}{|}_{{\overline{\omega }}_j}\frac{Z_n}{{\tau }_P^2},{\omega }_j={\overline{\omega }}_j{\tau }_P,j=1,2,3,t=\frac{\overline{t}}{{\tau }_P},z=\frac{\overline{z}}{Z_n},\hfill \\ {\nu }_{j1}=\left(\frac{\partial \overline{k}}{\partial \overline{\omega }}{|}_{j\overline{\omega }}-\frac{\partial \overline{k}}{\partial \overline{\omega }}{|}_{\overline{\omega }}\right)\frac{Z_n}{{\tau }_P},{\Delta }_{j1}k={\Delta }_{j1}\overline{k}{Z}_n,j=2,3\hfill \\ \gamma =\frac{2\pi {\chi }^{\left(2\right)}\overline{k}\sqrt{I_{01}}}{n^2}{Z}_n,\hfill \end{array}$$

(3)

where n is a refractive index at FF. Parameters ${\overline{k}}_j$ and ${\Delta }_{j1}\overline{k}$ are dimensional wave numbers and phase mismatching, respectively. ${\tau }_P$ denotes the duration of the incident FW pulse, and $I_{01}$ denotes its maximum intensity. $Z_n$ is a normalization length equal to 1 mm.

For computer simulation, we choose the following values of the physical parameters: ${\tau }_P$ is changed from 100 fs to 10 ps (meaning the pulses’ dispersion actions can be neglected) and the wavelength corresponding to wave with FF is equal to ${\lambda }_1=1064$ nm. The angle $\theta$ of the phase matching for the THG and the second-order susceptibility of selected crystals dependent on the angles $\theta$ and $\phi$ as well as the incident pulse intensity at which the coupling coefficient $\gamma$ is equal to 1 at the phase matching angle $\theta$, are shown in

Table 1. Phase matching angle $\theta$ of the process $\omega \left(o\right)+\omega \left(o\right)+\omega \left(o\right)=3\omega \left(e\right)$, dimensional phase mismatching between FW and SH. The second-order susceptibility of some crystals and the incident pulse intensity $I_0$, corresponding to value $\gamma =1$, in dependence of the laser pulse propagation direction with respect to the crystal axis: $(\theta ,\phi )$ are spherical coordinates. Phase matching angles are provided for FW length $\lambda =1064$ nm.

Crystal	THG Matching Angle $\mathit{\theta }(°)$	${\mathbf{\Delta }}_{21}\overline{\mathit{k}}$ (mm${}^{-1}$)	${\mathit{\chi }}^{\left(2\right)}$ (pm/V)	${\mathit{I}}_{0(\mathit{\gamma }=1)}$ (GW/cm${}^2$)
KDP [62]	64.97	147.403	$0.76sin\theta sin\left(2\phi \right)$	6.06
ADP [62]	66.58	163.755	$0.76sin\theta sin\left(2\phi \right)$	6.177
BBO [62]	37.48	140.934	$4.32(sin\theta -0.07cos\theta sin(3\phi \left)\right)$	0.695
KBBF [63]	31.16	87.371	$0.98cos\theta cos\left(3\phi \right)$	3.78
RBBF [61]	33.77	86.507	$0.9cos\theta cos\left(3\phi \right)$	4.84

. The angles are measured between the crystal axis in the uni-axial crystal and wave vector of the FW. References to the papers which contain the quadratic susceptibility coefficient are placed after the names of the crystals. Phase matching angles are calculated from the Sellmeier equations, which can be found in [58,59,60,61], respectively.

In Table 1, the dimensional phase mismatching between SH and FW is shown. Those values for the potassium beryllium fluoroborate (KBe₂BO₃F₂ or KBBF) and rubidium fluoroberyllium borate (RbBe₂BO₃F₂ or RBBF) crystals differ insignificantly. We use these values for computer simulation. The computer simulation results obtained for the other three crystals are similar, but they develop at longer distance in comparison with those required using the KBBF and RBBF crystals. It is important to note that at chosen intensities, the cubic nonlinearity only has a weak influence (the corresponding dimensionless coefficient equals a few hundredths), allowing us to neglect it.

The set of Equations (1) possess some invariants. The first invariant (the energy preservation law):

$$\begin{array}{cc}\hfill & I_1=\underset{0}{\overset{L_t}{\int }}\left(|A_1{|}^2+|A_2{|}^2+{|A_3|}^2\right)dt\hfill \end{array}$$

(4)

The second invariant:

$$I_2=\underset{0}{\overset{L_t}{\int }}\left(\sum _{j=1}^3{p}_j{A}_j\frac{\partial A_j^{\ast }}{\partial t}\right)dt,p_1=6,p_2=3,p_3=2$$

(5)

The third invariant (Hamiltonian):

$$\begin{array}{c}{I}_3=\underset{0}{\overset{L_t}{\int }}(-\sum _{j=1}^3{p}_j{D}_j{\left|\frac{\partial A_j}{\partial t}\right|}^2-p_2{\nu }_{21}Im\left(A_2^{\ast }\frac{\partial A_2}{\partial t}\right)-p_3{\nu }_{31}Im\left(A_3^{\ast }\frac{\partial A_3}{\partial t}\right)\\ +6\gamma Re\left(2A_1{A}_2{A}_3^{\ast }{e}^{i({\Delta }_{31}k-{\Delta }_{21}k)z}+A_1^2{A}_2^{\ast }{e}^{i{\Delta }_{21}k·z}\right)+3{\Delta }_{21}k|A_2{|}^2+2{\Delta }_{31}k{\left|A_3\right|}^2)dt,\end{array}$$

(6)

Let us notice that the integration limits change on $\{-\infty ,\infty \}$, respectively, if the laser pulse interaction is analysed in the unbounded time domain.

3. Induced Cubic Nonlinear Processes Based on Cascading SHG

3.1. Derivation of the Modified Equations

For simplicity of deriving modified equations describing THG, let us make a substitution of the TH complex amplitude as follows:

$$A_3={\overline{A}}_3{e}^{i{\Delta }_{31}kz}$$

(for brevity, the bar is omitted below). The Equations’ set (1) is transformed to the form

$$\begin{array}{c}\frac{\partial A_1}{\partial z}+i{D}_1\frac{{\partial }^2{A}_1}{\partial t^2}+i\gamma \left(A_1^{\ast }{A}_2{e}^{-i{\Delta }_{21}kz}+A_2^{\ast }{A}_3{e}^{i{\Delta }_{21}kz}\right)=0,\hfill \\ \frac{\partial A_2}{\partial z}+{\nu }_{21}\frac{\partial A_2}{\partial t}+i{D}_2\frac{{\partial }^2{A}_2}{\partial t^2}+i\gamma \left(A_1^2{e}^{i{\Delta }_{21}kz}+2A_1^{\ast }{A}_3{e}^{i{\Delta }_{21}kz}\right)=0,\hfill \\ \frac{\partial A_3}{\partial z}+{\nu }_{31}\frac{\partial A_3}{\partial t}+i{D}_3\frac{{\partial }^2{A}_3}{\partial t^2}+3i\gamma A_1{A}_2{e}^{-i{\Delta }_{21}kz}+i{\Delta }_{31}k{A}_3=0,0<z\le L_z,0<t<L_t.\hfill \end{array}$$

(7)

As we mentioned above, the three-wave interaction in a medium with quadratic susceptibility at cascading SHG leads to a response of a medium similar to the response of a medium with the cubic susceptibility. It is easy to show by using the multi-scale method [57]. Because the derivation is rather complicated, we place it in Appendix A. Here, we present the final result, in which the complex amplitudes $A_j,j=1,2,3$ are written as:

$$\begin{array}{cc}\hfill & A_1=U,A_2=\frac{1}{{\Delta }_{21}k}(-\gamma (U^2+2U^{\ast }W)e^{i{\Delta }_{21}kz}+v_1),A_3=W{e}^{i{\Delta }_{31}kz},\hfill \end{array}$$

(8)

accounting for the terms of the order of $O(1/{\Delta }_{21}k)$. The functions $U,W,v_1$ are governed by the set of nonlinear equations

$$\begin{array}{c}\frac{\partial U}{\partial z}+i{D}_1\frac{{\partial }^{2}U}{\partial t^2}-i\frac{{\gamma }^2}{{\Delta }_{21}k}{\left(\right|U|}^{2}U+3U^{\ast 2}W+{2U\left|W\right|}^2)=0,\hfill \\ \frac{\partial W}{\partial z}+{\nu }_{31}\frac{\partial W}{\partial t}+i{D}_3\frac{{\partial }^{2}W}{\partial t^2}-3i\frac{{\gamma }^2}{{\Delta }_{21}k}(U^3+{2\left|U\right|}^{2}W)+i{\Delta }_{31}kW=0,\hfill \\ \frac{\partial v_1}{\partial t}+{\nu }_{21}\frac{\partial v_1}{\partial t}+i{D}_2\frac{{\partial }^2{v}_1}{\partial t^2}=0\hfill \end{array}$$

(9)

with the initial conditions

$${U|}_{z=0}=A_{10}\left(t\right),W{{|}_{z=0}=0,v_1|}_{z=0}=\gamma A_{10}^2\left(t\right).$$

(10)

and zero-valued boundary conditions, because all pulses have a finite initial distribution and we consider the bounded domain in the z-coordinate. One can see that a term proportional to the third power of the FW complex amplitude is present in the second equation of the problem (9), and this term is responsible for THG. We call the below problems (9) and (10) the modified problem.

It should be emphasized that the first two equations in (9) can be obtained in much simpler way. For this aim, we propose to use the following representation of the complex amplitude $A_2$:

$$A_2=(A_2^{\left(0\right)}+A_2^{\left(1\right)}+A_2^{\left(2\right)}+\dots )e^{i{\Delta }_{21}kz}.$$

Substituting this series into the second equation of the set (7), we obtain the following equation in the first order of approximation on parameter $(1/{\Delta }_{21}k)$:

$$\frac{\partial A_2^{\left(0\right)}}{\partial z}+{\nu }_{21}\frac{\partial A_2^{\left(0\right)}}{\partial t}+i{D}_2\frac{\partial A_2^{\left(0\right)}}{{\partial }^2{t}^2}+i{\Delta }_{21}k{A}_2^{\left(0\right)}+i\gamma (A_1^2+2A_1^{\ast }{A}_3)=0,$$

that can be rewritten as

$$A_2^{\left(0\right)}+\frac{\gamma }{{\Delta }_{21}k}(A_1^2+2A_1^{\ast }{A}_3)-\frac{i}{{\Delta }_{21}k}\frac{{\partial }^2{A}_2^{\left(0\right)}}{\partial z}+\frac{{\nu }_{21}}{{\Delta }_{21}k}\frac{\partial A_2^{\left(0\right)}}{\partial t}+i\frac{D_2}{{\Delta }_{21}k}\frac{\partial A_2^{\left(0\right)}}{\partial t^2}=0.$$

The last three terms contain both complex amplitude $A_2^{\left(0\right)}$ and $1/{\Delta }_{21}k$. Therefore, they are much smaller than the two first terms, and they can be neglected if we account for the terms of order $1/{\Delta }_{21}k$ only. Thus, we obtain the following equality:

$$A_2^{\left(0\right)}=-\frac{\gamma }{{\Delta }_{21}k}(A_1^2+2A_1^{\ast }{A}_3).$$

(11)

Evidently, the terms in the right part describe phase gratings, induced at doubled frequency, where $A_1^2$ will lead to THG and FW self-modulation. The term $2A_1^{\ast }{A}_3$ will lead to cross-modulation of the waves and to reverse energy conversion from TH to FW. Substituting representation (11) into the first equation and the third equation of set (7)

$$\begin{array}{cc}& \frac{\partial A_1}{\partial z}+i{D}_1\frac{{\partial }^2{A}_1}{\partial t^2}-i\frac{{\gamma }^2}{{\Delta }_{21}k}\left(\right|A_1{|}^2{A}_1+3A_1^{\ast 2}{A}_3+2A_1|A_3{|}^2)=0,\hfill \\ & \frac{\partial A_3}{\partial z}+{\nu }_{31}\frac{\partial A_3}{\partial t}+i{D}_3\frac{{\partial }^2{A}_3}{\partial t^2}-3i\frac{{\gamma }^2}{{\Delta }_{21}k}(A_1^3+2|A_1{|}^2{A}_3)+i{\Delta }_{31}k{A}_3=0,\hfill \end{array}$$

coincides with the first two equations of set (9).

The equations in (9) possess the energy conservation law:

$$I_1=\underset{0}{\overset{L_t}{\int }}\left({\left|U\right|}^2+{\left|W\right|}^2\right)dt=const.$$

(12)

The values of the energy’s invariant (12) and (4) coincide with each other: they are equal to $\underset{0}{\overset{L_t}{\int }}{\left|A_{10}\right|}^{2}dt$. However, because the invariant for the set of the modified equations does not contain the SH wave, the total energy of two others waves will be greater than those computed using the original problem.

The Hamiltonian of the Equations’ set (9) is written in the form:

$$\begin{array}{cc}\hfill I_3& =\underset{0}{\overset{L_t}{\int }}(2{\nu }_{31}Im\left(W^{\ast }\frac{\partial W}{\partial t}\right)-6D_1{\left|\frac{\partial U}{\partial t}\right|}^2-2D_3{\left|\frac{\partial W}{\partial t}\right|}^2-\hfill \\ & \frac{3{\gamma }^2}{{\Delta }_{21}k}\left(4Re\left(U^3{W}^{\ast }\right)+{\left|U\right|}^4+{4\left|U\right|}^2{\left|W\right|}^2\right)+2{\Delta }_{31}k{\left|W\right|}^2)dt=const.\hfill \end{array}$$

(13)

It should be mentioned that the physical mechanism of a response, which inherent cubic nonlinearity, at the laser pulse propagation in a medium with the quadratic susceptibility under the cascading SHG is discussed in [51]. The authors stress that the self-lensing term ${\left|A\right|}^{2}A$, which describes the nonlinear phase shift, appeared because of sequences of cascaded photon processes, each involving up-conversion to $2\omega$ followed by back-conversion to $\omega$. This two-step process leads to an effective cubic susceptibility ${\chi }^{\left(2\right)}(\omega ;2\omega ,2\omega ){\chi }^{\left(2\right)}(2\omega ;\omega ,\omega )\to {\chi }_{eff}^{\left(3\right)}(\omega ;\omega ,2\omega ,\omega )$ that is characteristic of the effective four-photon interaction $\omega =\omega -\omega +\omega$. However, in that paper as well as in most others papers, the authors consider only the first two equations of the set (7) with respect to FW and SH.

We believe that the physical mechanism of this phenomenon can also be explained on the base of the analysis of the optical phase gratings induced by the interacting waves. So, if the phase mismatching ${\Delta }_{21}k$ is large, then the SH complex amplitude is proportional to $A_1^2$ (see above). This means that the response of a medium contains the phase grating at the doubled frequency. The optical wave at FF scatters on the grating and the optical wave at the tripled frequency appears. If the phase matching between the TH and FW occurs then the frequency conversion will be efficient. However, this mechanism of the frequency conversion has not been discussed.

Mathematically, we can obtain the term responsible for THG in the equation if we analyse all three equations in (7) and use the term $A_1^2$ instead of a small amplitude $A_2$ in the third equation. This is rough mathematical reasoning. Proper mathematical expression is derived on the basis of the multi-scale method (see Appendix A) or using our approach proposed above, which results in (11).

3.2. Exact Solution of the Modified Problem

In the framework of the long pulse duration approximation, the functions $U,W$ depend only on z coordinate, and we can rewrite the problems (9) and (10) in the form:

$$\begin{array}{cc}& \frac{dU}{dz}-i\frac{{\gamma }^2}{{\Delta }_{21}k}{\left(\right|U|}^{2}U+3U^{\ast 2}W+{2U\left|W\right|}^2)=0,\hfill \\ & \frac{dW}{dz}-3i\frac{{\gamma }^2}{{\Delta }_{21}k}(U^3+{2\left|U\right|}^{2}W)+{\Delta }_{31}kW=0,0<z\le L_z\hfill \\ & U\left(0\right)=1,W\left(0\right)=0.\hfill \end{array}$$

(14)

The conservation laws (12) and (13) are written as follows:

$$\begin{array}{cc}& I_1={\left|U\right|}^2+{\left|W\right|}^2,\hfill \\ & I_3=\frac{{\gamma }^2}{{\Delta }_{21}k}\left(-12Re\left(U^3{W}^{\ast }\right)-{3\left|U\right|}^4-{12\left|U\right|}^2{\left|W\right|}^2\right)+2{\Delta }_{31}k{\left|W\right|}^2.\hfill \end{array}$$

(15)

The energy conservation law of the problem (14) equals unity ($I_1=1$) due to the initial conditions.

To write the intensity evolution for the interacting waves, we use a well-known representation for the complex amplitudes:

$$U\left(z\right)=a_1{e}^{i{\phi }_1\left(z\right)},W\left(z\right)=a_3{e}^{i{\phi }_3\left(z\right)},$$

(16)

where $a_j,{\phi }_j$ are the real functions describing amplitudes and phases of the waves. In the new notations, the equations in (14) are re-written as follows:

$$\begin{array}{cc}& \frac{d{a}_1}{dz}=-\frac{3{\gamma }^2}{{\Delta }_{21}k}{a}_1^2{a}_{3}sin\phi ,\hfill \\ & \frac{d{a}_3}{dz}=\frac{3{\gamma }^2}{{\Delta }_{21}k}{a}_1^{3}sin\phi ,\hfill \\ & \frac{d\phi }{dz}-\frac{3{\gamma }^2}{{\Delta }_{21}k}\left(a_1\left(\frac{a_1^2}{a_3}-3a_3\right)cos\phi +a_1^2-2a_3^2\right)+{\Delta }_{31}k=0,\hfill \end{array}$$

(17)

where $\phi ={\phi }_3-3{\phi }_1$ is a phase difference between TH and FW. The invariants (15) are transformed to the following form:

$$\begin{array}{cc}& I_1=a_1^2+a_3^2=1,\hfill \\ & I_3=\frac{{\gamma }^2}{{\Delta }_{21}k}\left(-12a_1^3{a}_{3}cos\phi -3a_1^4-12a_1^2{a}_3^2\right)+{\Delta }_{31}k{a}_3^2.\hfill \end{array}$$

To express $cos\phi$, we modify the third invariant:

$${\tilde{I}}_3=I_3+3\frac{{\gamma }^2}{{\Delta }_{21}k}{I}_1^2=-12\frac{{\gamma }^2}{{\Delta }_{21}k}{a}_1^3{a}_{3}cos\phi -6\frac{{\gamma }^2}{{\Delta }_{21}k}{a}_1^2{a}_3^2+3\frac{{\gamma }^2}{{\Delta }_{21}k}{a}_3^4+2{\Delta }_{31}k{a}_3^2=0,$$

which equals zero because the TH incident intensity is equal to zero: $a_3\left(0\right)=0$. Hence, we can express $cos\phi$:

$$cos\phi =\frac{a_3(-6a_1^2+3a_3^2+2q)}{12a_1^3},$$

(18)

where the parameter q is defined as:

$$q=\frac{{\Delta }_{21}k{\Delta }_{31}k}{{\gamma }^2},$$

which characterizes relation between phase incursions because of phase mismatching between SH and FW, and TH and FW, respectively, and because of the nonlinear coupling coefficient $\gamma$. Obviously, the inequality

$$|cos\phi |\le 1$$

must be valid. This leads to restriction on the intensities changing. Then, we find $sin\phi$ from equality (18) and substitute it into the second equation of (17), which is preliminarily multiplied by $a_1$. As a result, we write the following equation with respect to the TH intensity:

$$\frac{d{a}_3^2}{dz}=\pm 6\frac{{\gamma }^2}{{\Delta }_{21}k}{a}_1^3{a}_3\sqrt{1-\frac{a_3^2{(-6a_1^2+3a_3^2+2q)}^2}{144a_1^6}}.$$

Using the new variable $p_3=a_3^2$ and expressing $a_1^2$ from the first invariant $a_1^2=1-p_3$, we finally obtain

$$\begin{array}{cc}& \frac{d{p}_3}{dz}=\pm 7.5\frac{{\gamma }^2}{{\Delta }_{21}k}\sqrt{p_{3}f\left(p_3\right)},\hfill \\ & f\left(p_3\right)=\frac{16}{25}({(1-p_3)}^3-\frac{1}{144}{p}_3{(9p_3+2q-6)}^2).\hfill \end{array}$$

(19)

where the sign in the right side of the equation corresponds to different processes—the plus sign corresponds to the conversion of the FW energy to TH energy and the minus sign to reverse energy conversion from the wave with tripled frequency to the FW. So, since TH incident intensity equals zero, we choose the plus sign.

The intensity evolution is defined by the parameter q because it defines the number of real roots of the equation

$$f\left(p_3\right)=0,$$

which can be found, for example, by using the Sturm theorem. We denote the roots as $P_{31},P_{32},P_{33}$, where $P_{31}$ is always real. If two other roots are also real, then there is an inequality $P_{31}\le P_{32}\le P_{33}$ between them. Obviously, this assumption does not restrict our consideration. It should be also emphasized that a solution of the Equation (19) is expressed through the elliptical functions or, in certain cases, through elementary functions. Below, we briefly discuss the problem’s solution corresponding to various values of the parameter q.

Moreover, we develop also the TH phase evolution. Because of the complexity of the obtained formulas, we do not write them here. However, using the formulas with respect to the phase and the intensity, it is possible to approximate the solution of problems (1) and (2) with good accuracy. We show it below using the computer simulation results.

• Single frequency conversion mode: domain $q<-1.5,q>3+\frac{9}{\sqrt[3]{2}}-9\sqrt[3]{2}$; two complex roots $P_{32},P_{33}$.

Only root $P_{31}$ is positive. As a result, the TH intensity $p_3\left(z\right)$ changes between zero and $P_{31}$ in accordance to the following formula:

$$p_3\left(z\right)={\left|A_3\left(z\right)\right|}^2=\frac{c{P}_{31}(1-cn(7.5\frac{{\gamma }^2}{{\Delta }_{21}k}\kappa z,m))}{(d-c)cn(7.5\frac{{\gamma }^2}{{\Delta }_{21}k}\kappa z,m)+(c+d)},$$

(20)

where

$$\begin{array}{cc}& c=\sqrt{r^2+s^2},d=\sqrt{{(p_{31}-r)}^2+s^2},\kappa =\sqrt{cd},m=\frac{{\kappa }^2+(r+P_{31})r-s^2}{2{\kappa }^2},\hfill \\ & r=Re\left(P_{32}\right),s=Im\left(P_{32}\right).\hfill \end{array}$$

This domain of the parameter q contains an important case of THG corresponding to the phase matching between FW and TH (${\Delta }_{31}k=0$). Under this condition, there is only a unique branch of the TH intensity changing between 0 and $P_{31}=0.8$. Two other roots are defined by a solution of the equation $5p_3^2-8p_3+4=0$, and they are the complex ones ($P_{32,33}=0.8\pm 0.4i$). As a result, the TH intensity evolution is described by:

$$p_3\left(z\right)=\frac{4(1-cn(\frac{3\sqrt[4]{5}{\gamma }^2}{{\Delta }_{21}k}z;0.5(1-\frac{1}{\sqrt{5}})))}{(5+\sqrt{5})-(5-\sqrt{5})cn(\frac{3\sqrt[4]{5}{\gamma }^2}{{\Delta }_{21}k}z;0.5(1-\frac{1}{\sqrt{5}}))}.$$

(21)

As follows from (21), the coupling coefficient $\gamma$ and the phase mismatching between FW and SH (${\Delta }_{21}k$) only influence the intensity oscillation period. Let us note that in all other modes (at ${\Delta }_{31}k\ne 0$) of THG, there is more complicated dependencies of the TH maximum intensity and period on ${\Delta }_{21}k$. Therefore, we do not discuss this below.

• Bistable frequency conversion mode: domain $-1.5<q<3+\frac{9}{\sqrt[3]{2}}-9\sqrt[3]{2}$; three real roots.

In this case, all three roots $P_{3j}$ are positive. Therefore, there are two regimes of the wave interaction: the low-efficiency mode and the high-efficiency mode. In the low-efficiency regime, the TH intensity changes between 0 and $P_{31}$ as follows:

$$\begin{array}{cc}& p_3\left(z\right)=\frac{P_{31}{P}_{33}{sn}^2(3.75\frac{{\gamma }^2}{{\Delta }_{21}k}\sqrt{P_{32}(P_{33}-P_{31})}z,m)}{P_{31}{sn}^2(3.75\frac{{\gamma }^2}{{\Delta }_{21}k}\sqrt{\left(P_{32}\right)(P_{33}-P_{31})}z,m)+(P_{33}-P_{31})},\hfill \\ & m=\sqrt{\frac{P_{31}(P_{33}-P_{32})}{P_{32}(P_{33}-P_{31})}}.\hfill \end{array}$$

(22)

In the high-efficiency regime, the TH intensity changes between values of the roots $P_{32}$ and $P_{33}$ according to the rule:

$$\begin{array}{c}{p}_3\left(z\right)=\frac{P_{32}{P}_{33}}{(P_{33}-P_{32}){sn}^2(3.75\frac{{\gamma }^2}{{\Delta }_{21}k}\sqrt{P_{32}(P_{33}-P_{31})}z,m)+P_{32}}\hfill \end{array}$$

(23)

with the same coefficient m as in (22).

• Single frequency conversion mode and mode of preserving TH maximum intensity ($P_3=1$): the first special case $q=-1.5$; $0<P_{31}<P_{32}=P_{33}=1$.

For this parameter value q, the TH intensity $p_3\left(z\right)$ changes again between 0 and $P_{31}=0.64$. However, the solution is expressed through the cosine function:

$$p_3\left(z\right)=\frac{P_{31}(1-cos(6\frac{{\gamma }^2}{{\Delta }_{21}k}\sqrt{1-P_{31}}z))}{-P_{31}cos(6\frac{{\gamma }^2}{{\Delta }_{21}k}\sqrt{1-P_{31}}z)+2-P_{31}}=\frac{8(1-cos(4.5\frac{{\gamma }^2}{{\Delta }_{21}k}z))}{-8cos(4.5\frac{{\gamma }^2}{{\Delta }_{21}k}z)+17}.$$

(24)

Another mode, which corresponds to the case $q=-1.5$, is the mode of the unchangeable intensity: $p_3\left(z\right)\equiv 1$ (correspondingly, FW intensity equals zero). To reach this mode, it is necessary to switch the THG mode from this maximum value $P_{31}=0.64$ of the frequency conversion efficiency to its value $P_{32}=P_{33}=1$. In practice, it may require strong conditions.

• Self-switching of generation modes from low-efficiency to high-efficiency: the second special case $q=3+\frac{9}{\sqrt[3]{2}}-9\sqrt[3]{2}$; $0<P_{31}=P_{32}<P_{33}<1$.

In the case $q\approx -1.2$, there are two regimes of the frequency conversion. In both of them (low-efficiency and high-efficiency), the TH intensity $p_3$ evolution is governed by:

$$\begin{array}{cc}\hfill & z\left(p_3\right)=\frac{1}{7.5\frac{{\gamma }^2}{{\Delta }_{21}k}\sqrt{P_{31}(P_{33}-P_{31})}}ln\frac{{\left(\sqrt{(P_{33}-P_{31})p_3}+\sqrt{P_{31}(P_{33}-p_3)}\right)}^2}{|p_3-P_{31}|P_{33}},\hfill \end{array}$$

(25)

where the modulus is expanded with a positive sign at the highly effective regime $P_{31}<p_3<P_{33}$ and with a negative sign at the low-effective regime $0<p_3<P_{31}$.

From (25) together with computations of the equation $f\left(p_3\right)=0$ roots, it is possible to compute the THG maximum efficiency: its value is more than $95\%$ and it is achieved at $q=3+\frac{9}{\sqrt[3]{2}}-9\sqrt[3]{2}$ (Figure 2).

Therefore, the phase matching between TH and FW is not an optimal condition for THG because at ${\Delta }_{31}k=0$ ($q=0$) only $80\%$ of the FW energy converts to the TH energy. One can see also in Figure 2 the abrupt decreasing TH maximum intensity near the optimal negative value of the parameter q. With increasing positive phase mismatching ${\Delta }_{21}k$, the negative value of the phase mismatching ${\Delta }_{31}k$ tends toward zero. Let us note that the efficiency of THG is the same if we change the signs of phase mismatches to opposite ones because the value of parameter q does not change at the same time (for example, ${\Delta }_{21}k=20,{\Delta }_{31}k=-1$ is changed to ${\Delta }_{21}k=-20,{\Delta }_{31}k=1$).

It should be emphasized that the THG based on the proposed mechanism differs from that occurring in a medium with a cubic nonlinear response by absence of the self-modulation process at the tripled frequency, which usually decreases frequency conversion efficiency (see the set of equations (9)). The cross-modulation affect, arising from FW, is three times less in comparison with THG based on the cubic nonlinear response. Moreover, a certain choice of ${\Delta }_{31}k$ can decrease the negative influence of pulse self-action that one can see in Figure 2.

4. Influence of Pulse Duration of FW on Frequency Conversion Due to Cascading Processes

The results obtained in the framework of the modified problem are supported below (Figure 3) by the computer simulation provided using the original problem with the parameters $\gamma =2$ and ${\Delta }_{21}k=87$. The phase mismatching between waves with FF and tripled frequency (${\Delta }_{31}k$) is varied for achieving the corresponding value of the parameter q to realize the THG modes described above.

First, we present computer simulation results using long pulse duration approximation to demonstrate the applicability of provided analysis. Then, we demonstrate an achievement of high efficiency’s frequency conversion for the shorter pulses of picoseconds or even femtoseconds.

4.1. Long Pulse Duration Approximation

In Figure 3a, the conversion efficiency evolution is shown if the phase matching between TH and FW (${\Delta }_{31}k=0$) occurs. Its maximum value ($80\%$) is achieved in the section $z=15$. Both lines coincide with each other with high accuracy. However, at phase mismatching ${\Delta }_{31}k=-0.0549878$ corresponding to the parameter $q=3+\frac{9}{\sqrt[3]{2}}-9\sqrt[3]{2}$, the coincidence between both lines occurs only until the pulse propagation distance equaling $z=12.5$ and the maximum intensity $|A_3{|}^2=0.7$ is achieved at computing based on the original problem (Figure 3b). It occurs because this phase mismatching corresponds to value dividing two modes of the frequency up-conversion. Therefore, for adequate description of the tripling frequency generation, it is necessary either to take into account the next order of the series corresponding to $O\left({\Delta }_{21}{k}^{-2}\right)$ or to increase computer simulation accuracy. Its insufficient accuracy leads to the switching between the up-conversion modes because of computation round-off and its accumulation. For confirmation of this hypothesis, we compute additional frequency conversion at phase mismatching which differs a little bit from those. So, at slight changing phase mismatching (${\Delta }_{31}k=-0.053$), the THG efficiency increases up to $94.5\%$ (Figure 3c) that is essentially greater than in the previous case. The difference in solutions appears at much longer distances of the wave interaction. At further increasing phase mismatching ${\Delta }_{31}k=-0.05$ (Figure 3d), the frequency conversion efficiency decreases insufficiently ($93\%$), but the coincidence between the two solutions improves essentially.

The insets in Figure 3a,c also confirm the partial compensation of the phase nonlinear incursion by the phase mismatching ${\Delta }_{31}k$. The energy conversion from FW to TH occurs if the phase difference $\phi$ is negative. In Figure 3a, this takes place if the distance of laser pulse interaction is less than 15 dimensionless units. Meanwhile, that distance is about $z=45$ in Figure 3c. Both mentioned sections of the z-coordinate corresponds to the phase difference equal to $\phi =0$.

To highlight an advantage of our approach to the solution of the problem, we also show in Figure 3 the original problem’s solution in the pump energy non-depletion approximation. One can see that in all cases, the solution of the modified problem approximates the solution of the original problem much better than the pump energy non-depletion solution. In fact, the last one does not even show key features of the problem solution, and is close to two other solutions only when the TH intensity is smaller than $0.2$.

4.2. Pulses with Picosecond or Femtosecond Duration

In the case under consideration, the SOD of the pulses and especially GVM restricts the frequency’s conversion efficiency. Nevertheless, a significant part of the FW energy still converts to the wave with tripled frequency. For illustration of the advantage of the derived solution of the problem, we show also the frequency conversion’s efficiency computed using the FW energy non-depletion approximation.

The incident pulse at the FF possesses the Gaussian distribution of its shape without any chirp:

$$A_{10}\left(t\right)=A_{01}{e}^{-{(t-0.5L_t)}^2},$$

(26)

where $A_{01}$ is equal to unity due to chosen dimensionless parameters.

Computer simulation results are presented in Figure 4. We show FW with red lines, SH with blue lines, and TH with green lines. The computer simulation results of the original problem are shown with solid lines, the computer simulation results of the modified problem with dots, and results of using the wave energy non-depletion approximation with dashed lines. The incident FW pulse distribution is shown by the black line with triangles.

The problem parameters are equal to $\gamma =2,{\Delta }_{21}k=87,{\Delta }_{31}k=-0.05$. The SOD and GVM coefficients equal ${\nu }_{21}=0.07,{\nu }_{31}=0.21$ and $D_1=-0.0000032,D_2=-0.000029,$$D_3=-0.000043$, respectively, corresponding to a 1 ps pulse duration of the FW with a 1064 nm wavelength at its propagation in the KBBF crystal. We show the pulses’ amplitudes to show better differences between the problem solutions obtained in various approximations.

In Figure 4a, one can see the that coincidence between the solutions of the modified and original problems is perfect. Because of the GVM of presence, the TH pulse moves from the FW one, which results in an asymmetrical shape of the TH pulse. The non-depletion energy approximation for the pump pulse fails to give the adequate description of the wave interaction predicting the super Gaussian pulse shape with much more energy. The only thing that resembles the computer simulation results of the problem is the TH maximum intensity, which is a little bit higher than it should be. As shown in Figure 4d, by using the pulse’s energy non-depletion approximation at the wave with FF, the TH energy exceeds initial FW energy near the section $z=20$, which is completely wrong. In fact, the THG efficiency is about $51\%$, reaching this value in the section $z=40$, and then it remains stable.

If the pulse duration at the FF is ten times less (100 fs), then the SOD and GVM coefficients increase 100 and 10 times, respectively. In this case, GVM strictly decreases the frequency conversion efficiency: only $11\%$ of the FW energy transfers to TH pulse (Figure 4e) with a small maximum intensity (Figure 4e) and long duration in comparison with the initial FW pulse (Figure 4b). Due to the lower conversion efficiency, the pump wave energy non-depletion approximation works better; however, it is still worse than our approach.

It can be seen in Figure 4e, that the dashed line (energy non-depletion approximation) is close to the solid line (numerical simulation of the original problem), but the dotted line (our approach) coincides fully with the solid one.

The conversion efficiency may be increased by the growing incident intensity of the FW. For example, if the parameter $\gamma$ is increased up to 4 (corresponds to increasing incident intensity in four times), then the conversion efficiency reaches $38\%$ (Figure 4f). Its change along the z-coordinate is shown in Figure 4d. The frequency conversion efficiency grows until the section $z=10$ and then it remains almost constant. The TH pulse shape in Figure 4c looks similar to the one shown in Figure 4a, but its maximum intensity is two times less. The pump wave energy non-depletion approximation demonstrates bad results as follows from Figure 4c. In turn, our approach, based on multi-scale method, shows good coincidence with the problem solution: the dotted lines in Figure 4c,f coincide with the solid one.

It should be noted that cubic nonlinear response of a medium may play an important role if the incident intensity of the FW becomes high enough. However, for brevity, we do not study this influence in the current paper.

5. THG Efficiency in Dependence of SOD Coefficients and Phase Mismatching between Interacting Waves

Below, we discuss how each of the factors (phase mismatching, incident intensity of FW, SODs coefficients) influences the THG efficiency due to the cascading SHG under the condition of the pulses’ group velocity matching or under long enough pulses’ durations to avoid GVM’s influence: ${\nu }_{21}={\nu }_{31}=0$. Computer simulation is provided at ${\Delta }_{31}k=0$ using the sets of equations in (1) and (9), respectively, and the corresponding simulation results are compared. We pay attention to reaching the high-efficiency THG with a smooth pulse shape and smooth spectrum of a wave at the tripled frequency.

It should be emphasized that below, in Paragraph 5.1 and Paragraph 5.2 and Paragraph 5.3, and also in Figures 11 and 14, the solution of the problems (1) and (2) is shown as a solid line while the solution of the modified problems (9) and (10) is shown as a dotted line. Numbers $1,2,3$ in the figures denote FW, SH, and TH, respectively. In some cases, we also demonstrate the applicability of the analytical solution derived from the long pulse duration approximation and this solution is shown as a dashed line. The initial spectrum distributions (as well as the FW initial pulse shape) in some figures are shown via solid lines with triangles.

5.1. Normal Dispersion of a Medium at FW Wavelength

First of all, we consider the frequency tripling of the pulse with wavelength $\lambda =800$ nm and propagating in a BBO crystal with the normal dispersion, which corresponds in our notation to negative values of the dimensionless SOD coefficients. In Figure 5a–c, the pulse shapes are shown for the frequency up-conversion at a relatively small phase mismatching between SH and FW (${\Delta }_{21}k=20$) and at the incident pulse duration of about 100 fs. It should be stressed that a positive sign of phase mismatching ${\Delta }_{21}k$ corresponds to the FW pulse compression under its propagation in a medium with the normal SOD if the THG is not essential.

We see a good coincidence between the FW and TH pulse shapes, computed using the original problem solution and the modified one, until the section $z=13$. Despite a medium possessing only the quadratic nonlinearity, the TH intensity is high: it achieves $0.75$ dimensionless units at the pulse center. This growth is accompanied by the spectral broadening of the FW and TH pulses in comparison with the incident FW pulse spectrum. However, the FW pulse spectrum has two local minimums in this section of a medium.

After the distance equals 10 units (Figure 5e), the FW intensity begins its growth. This happens due to the sub-pulse formation at the FF. Indeed, in the section $z=7$ (Figure 5a), the FW pulse shape is smooth, and it becomes close to plate distribution near the pulse center. After this section of a medium, the intensity grows in this part of the pulse (Figure 5b), and its increase becomes very remarkable (Figure 5c).

The conversion efficiency achieves $51\%$ and, at the same distance, the TH spectra (Figure 5d) has a smooth shape without any local minima. The analytical solution (21) perfectly describes the wave interaction until the section $z=7.5$. Then, we see that the FW pulse self-focusing and TH efficiency, predicted by the theory, are slightly higher than the corresponding values obtained from computer simulation results. Nevertheless, the set of equations (9) correctly describes a self-focusing of the pulses, as well as the presence of the effective cubic nonlinearity under wave interaction occurring at the cascading SHG.

If the sign of the phase mismatching between FW and SH is negative (${\Delta }_{21}k=-20$), then the energy conversion from the FW to TH also occurs with high efficiency (Figure 6). Nevertheless, there are new essential features which accompany this process.

Therefore, this sign of the phase mismatching corresponds to FW pulse decompression due to the induced cubic nonlinear response in a medium with normal dispersion. Therefore, the pulse duration of the FW and TH is greater than those with positive phase mismatching (compare Figure 5 and Figure 6). However, despite the FW pulse decompression, the bandwidth of the TH pulse spectrum is greater than the FW incident pulse spectrum. This means that the essential spectral broadening of the TH pulse occurs and this pulse can be additionally compressed. Moreover, the energy conversion reaches the value of about $60\%$ at increased crystal length in comparison with the previous case: the crystal length must be about 3 cm. Other features of the THG process are similar to the generation process upon positive phase mismatching.

5.2. Anomalous Dispersion of a Medium at FW Wavelength

Let us briefly discuss THG in a medium with the anomalous dispersion occurring at the frequency of FW. The SOD coefficient $D_1$ is positive. We consider a propagation of the incident pulse with duration about 100 fs and with wavelength about 1064 nm in KDP crystal (see, for instance, [64]). If ${\Delta }_{21}k>0$ that corresponds to decompression of the FW pulse and compression of the pulse at tripled frequency, then the efficiency of the FW energy conversion to the TH achieves about $60\%$ in the section $z=15$ (Figure 7). It is accompanied with good TH pulse shape and its broadened spectrum in comparison with the incident pulse spectrum at the FF. The FW pulse shape has three intensity maxima and possesses smaller intensity at the pulse center in comparison with its value occurring at the pulse propagation in a medium with the normal dispersion at the FF. Nevertheless, we can see in Figure 7 that the FW pulse spectrum is more broadened in comparison with the incident pulse spectrum.

If the phase mismatching between FW and SH is negative (${\Delta }_{21}k<0$), then the pulse compression at FF occurs, and decompression of the SH and TH pulses is observed. This leads to highly effective THG (Figure 8): more than $80\%$ of the FW energy converts to the TH wave. However, the TH intensity is less than its value at a positive phase mismatching ${\Delta }_{21}k$, and the TH spectrum is smooth and broadened.

5.3. THG at Large Positive Phase Mismatching between SH and FW

If the phase mismatching ${\Delta }_{21}k$ increases until ${\Delta }_{21}k=250$, then the required distance at which the frequency tripling high efficiency is achieved also increases (Figure 9).

Moreover, the maximum TH intensity is not so high as in previous Figures and does not exceed $5\%$ from the incident FW maximum intensity. On the other hand, the TH pulse duration is approximately ten times greater in comparison to the incident FW pulse duration. The energy conversion efficiency is greater than $50\%$ in the section $z=500$ (about 50 cm), which is a very long distance. For shorter crystal length ($z=100$ or 200), the energy conversion efficiency is about $25\%$ and $40\%$, respectively. The pulse shape of TH is not a Gaussian one: it is close to the flat distribution near the pulse center and to the hyperbolic function near the pulse front and near its travelling part, respectively (Figure 9a).

It should be emphasized that using an incident pulse with a longer duration (for example, approximately 1 ps), one can achieve an approximate $62\%$ frequency conversion efficiency (Figure 10) for the crystal with length equal to 20 cm. At the same time, the maximum intensity of the TH pulse becomes high (regarding unity, compare Figure 9c and Figure 10c) and, in contrast to the previous case, the TH pulse spectrum bandwidth is greater than the bandwidth of the incident FW pulse spectrum (Figure 10b). Therefore, the TH pulse may be compressed further. It is also essential to emphasize that the analytical solution describes the frequency tripling process with high accuracy. This statement is especially valid for the pulse spectra, for which the differences between all three distributions are negligible (see Figure 10).

5.4. Developing of FW MI at Negative ${\Delta }_{21}k$

The THG efficiency can be higher at a negative phase mismatching between SH and FW than its corresponding value achieving at a positive one. However, the frequency conversion efficiency may be restricted by MI developing. For example, this occurs if the pulse duration is equal to 1 ps.

The MI develops firstly at the wave with FF (Figure 11). So, despite the pulse being smooth in the section $z=68$ (Figure 11a), the weak oscillations are already observed in the section $z=74$ (Figure 11b) near the pulse center. Then, the intensity oscillations enhance and a certain periodic structure appears in the section $z=82$ (Figure 11b) with about 40 sub-pulses with the dimensionless duration $0.03$ (approximately $30$ fs). Nevertheless, the TH pulse shape remains smooth at it propagates along certain distance, during which the THG efficiency grows continuously despite the development of the oscillating process. Furthermore, the situation changes dramatically in the section $z=100$ (Figure 11d), for example, the TH maximum intensity becomes weak. In contrast, the SH intensity begins to grow and we observe two sub-pulses which leave the area of the wave interaction. Their velocities of motion are constant, and their phase distributions are linear. Moreover, the sub-pulses propagate without any changes to their shape nor maximum intensities as one can see in Figure 12b. Thus, we see the soliton propagation of a wave at doubled frequency.

To prove that the MI causes these oscillations, we present the computer simulation results in

Table 2. Developing of the MI in dependence on incident FW pulse duration and its maximum intensity.

$\mathit{\gamma }$	${\mathbf{\Delta }}_{21}\mathit{k}$	${\mathit{D}}_1$	${\mathit{D}}_2$	${\mathit{D}}_3$	MI	${\mathit{z}}_{\mathit{MI}}$	$|{\mathit{A}}_3({\mathit{z}}_{\mathit{MI}},0.5{\mathit{L}}_{\mathit{t}}){|}^2$	${\mathit{\eta }}_{\mathit{max}}$
1	−250	−0.000032	−0.000083	−0.000199	Yes	74	0.4225	0.34
1	−250	−0.00032	−0.00083	−0.00199	Yes	160	0.2894	0.51
1	−250	−0.00064	−0.00166	−0.00398	No	–	0.1667	0.5
$\sqrt{2}$	−250	−0.00064	−0.00166	−0.00398	Yes	116	0.231	0.62

. As is known, MI’s development depends on the SOD values. So, if the influence of the SOD increases and the other parameters remain the same, then MI does not appear. With this goal, we increase the SOD coefficients ten times. However, the MI is still observed, but the required distance for its appearing increases up to $z=160$. As a consequence, the THG efficiency becomes higher: about $51\%$. After that, we choose SOD coefficients 20 times greater, and the MI does not appear at all. However, because the SOD influence enhances in comparison with the previous case, then the THG efficiency is a bit lower ($50\%$). Increasing $\gamma$ by the factor $\sqrt{2}$ (which means that the effective cubic nonlinearity is two times greater) results in the MI appearing again. Due to enhancing the nonlinear coefficient, the THG efficiency grows ($62\%$) at the decreased length of a medium (see the fourth row in Table 2). Thus, if the phase mismatching between FW and SH is negative, one must carefully choose the incident pulse intensity to avoid the MI developing and to achieve high conversion efficiency.

It should be emphasized briefly that the physical reason for MI appearing differs from the corresponding one occurring in a medium with the quadratic nonlinear response [65,66], because in our case it is observed in a medium in which all three pulses experience decompression. This means that a medium possesses the normal dispersion at all frequencies. However, this condition does not restrict MI’s appearance, and it may be observed at various combinations of SOD types.

5.5. Dependence of the Crystal Length on ${\Delta }_{21}k$ for Reaching Maximum THG Efficiency

To clarify the dependence of the energy conversion efficiency on the phase mismatching between FW and SH, we provide computer simulations of the problem for different values of ${\Delta }_{21}k$, changing from $-200$ to 200, at a propagation of the pulse with its incident duration being equal to $1ps$. Our aim is the definition of the highest conversion efficiency and the corresponding crystal length required for this. The obtained results, computed at a phase matching between TH and FW, are shown in Figure 13.

One can see that at appropriate choice of the parameters, the highest efficiency can be achieved if the phase mismatching between FW and SH is close to zero, but does not equal to zero: for example, ${\Delta }_{21}k=-4$ (Figure 14).

The maximum TH intensity achieves a value close to unity, and the TH pulse shape and its spectrum are smooth, as are as the spectra of other pulses. More than $80\%$ of the incident FW pulse energy transfers to the TH even in the section $z=6.5$. However, taking into account the phase mismatching between TH and FW (${\Delta }_{31}k=0$), such a situation may be hardly realised (or it is impossible at all). So, if we could achieve phase matching between FW and SH as well as between FW and TH simultaneously, we could obtain the comb SH and TH generation. However, their sum efficiency will be less than $80\%$.

We see in Figure 13 that for large positive values of phase mismatching between SH and FW, the TH maximum conversion efficiency remains at a high level despite the changing of ${\Delta }_{21}k$ in a wide range. However, the crystal length for its achievement also increases because SH is absent; therefore, this mode corresponds to cascading SHG. For a small negative value of the phase mismatching ${\Delta }_{21}k$, the maximum efficiency of the FW energy conversion to TH is about 3–4% higher than for the positive one. However, this requires increasing of the crystal length. Moreover, at large negative values of ${\Delta }_{21}k$, the efficiency conversion is restricted by an appearance of the MI, as described above.

If the incident pulse duration is decreased until $100fs$, then the SOD role becomes more remarkable (Figure 15), and in this case, a longer distance is required for the achievement of the maximum THG energy at a positive phase mismatching ${\Delta }_{21}k$. On the other hand, at a negative ${\Delta }_{21}k$, the MI appears in rather small propagation distance. Nevertheless, it is also possible to achieve high efficiency of the energy conversion. The verification of these statements is shown in Figure 16.

Indeed, the MI develops in the section $z=50$ at ${\Delta }_{21}k=-10$ (Figure 16); meanwhile, there is not any MI at ${\Delta }_{21}k=-12$ until the section $z=75$ (Figure 16b). Essentially, the high-efficiency THG is achieved before the section of the MI develops. At ${\Delta }_{21}k=-14,$ the MI development starts only from the section $z=95$ (Figure 16c). However, the maximum THG efficiency decreases because of the decreasing of the effective cubic nonlinear response and increasing of SOD influence.

Another feature of interest, which is presented in Figure 15, is very long distance, required for the achievement of the high-frequency conversion efficiency if ${\Delta }_{21}k<-40$. In fact, this distance is longer than $z=500$ ($50$ cm), until which the computation was provided. One can see two modes of the wave interactions in Figure 16d–f. In the first one, the fast energy conversion from FW to TH occurs, but in the second one, the energy converts between waves very slowly. Moreover, a part of TH energy converts to FW after the achievement of the maximum THG efficiency at ${\Delta }_{21}k=-20$ (Figure 16d) and $-40$ (Figure 16e), and then the quasi-stable mode of the wave propagation occurs. If the phase mismatching between FW and SH is decreased until ${\Delta }_{21}k=-80$ (Figure 16f), this quasi-stable mode appears before the maximum THG efficiency achievement due to domination of the SOD’s influence over influence of the effective cubic nonlinearity on the pulse propagation. In the section $z=500$, the THG efficiency is higher than its value reached in the section $z=100$. However, their difference is tiny.

5.6. Remark about the Frequency Up-Conversion Efficiency in Dependence of Incident Pulse Intensity

The corresponding results are presented in

Table 3. Influence of $\gamma$ on the frequency tripling process at ${\Delta }_{31}k=0,D_1=-0.0032,$$D_2=-0.0083,D_3=-0.0199$.

$\mathit{\gamma }$	${\mathbf{\Delta }}_{21}\mathit{k}$	${\mathit{\eta }}_{\mathit{max}}$	${\mathit{z}}_{\mathit{max}}$	${\mathit{\tau }}_3$	${\mathit{\tau }}_3^{\mathit{\omega }}$	TH Pulse Shape
0.5	20	0.65	80	1.671	1.557	Smooth
0.5	−20	0.48	100	10	3.519	Not smooth
1	20	0.51	13	0.804	3.15	Smooth
1	−20	0.67	30	4.925	5.924	Rather smooth
2	20	0.61	4.2	0.7176	3.649	Smooth
2	−20	0.83	7	1.772	4.280	Smooth

. At low-incident FW intensity (and, consequently, small-value $\gamma$), it is necessary to increase the length of the medium (fourth column of the table) to achieve the high-frequency conversion efficiency. At $\gamma =0.5$ and negative ${\Delta }_{21}k$, the conversion efficiency is less than achieved efficiency at positive ${\Delta }_{21}k$, and it is lower than $50\%$. This happens due to a big pulse decompression. However, upon the increasing of the incident intensity ($\gamma =1$), the situation changes completely due to action of the self- and cross-modulation. So, the THG efficiency decreases down to $51\%$ at positive ${\Delta }_{21}k$. In contrast, at negative ${\Delta }_{21}k$, the enhancement of the incident pulse intensity results in high energy conversion efficiency ($67\%$).

For a doubled value of the parameter $\gamma$ and at the positive phase mismatching between SH and FW (sixth row), the energy conversion efficiency already achieves about $60\%$ in the section $z=4.2$ (in comparison with $z=13$ at $\gamma =1$). Nevertheless, we observe a more or less remarkable generation of a wave at doubled frequency at opposite signs of the phase mismatching (${\Delta }_{21}k=-20$, seventh row), and more than $80\%$ of the incident pulse energy is converted to TH. However, this requires a longer crystal ($z=7$) than at ${\Delta }_{21}k=20$. The growing incident intensity at ${\Delta }_{21}k<0$ may also result in MI’s appearance, which restricts the TH intensity.

5.7. Crucial Influence of Phase Mismatching between TH and FW on Frequency Conversion Efficiency

Below, we consider the influence of weak phase mismatching between FW and TH (${\Delta }_{31}k=\pm 0.01,\pm 0.1,\pm 0.2,\pm 0.4$) under a rather large phase mismatching between FW and SH (${\Delta }_{21}k=\pm 20$) on THG’s efficiency. Computer simulation results, obtained for the incident pulse duration $100$ fs in both signs of SOD, are shown in

Table 4. Influence of ${\Delta }_{31}k$ on the frequency tripling process at $\gamma =1,{\Delta }_{21}k=20$.

Normal SOD at FF						Anomalous SOD at FF
${\Delta }_{31}k$	${\eta }_{max}$	$z_{max}$	${\tau }_3$	${\tau }_3^{\omega }$	TH pulse shape	${\Delta }_{31}k$	${\eta }_{max}$	$z_{max}$	${\tau }_3$	${\tau }_3^{\omega }$	TH pulse shape
0	0.51	13	0.804	3.15	Smooth	0	0.58	15	0.79	3.675	Smooth
0.01	0.49	13	0.802	3.104	Smooth	0.01	0.54	16	0.827	3.578	Smooth
−0.01	0.54	14	0.787	3.105	Smooth	−0.01	0.58	17	0.796	3.661	Smooth
0.1	0.48	11	0.849	3	Smooth	0.1	0.43	12	0.873	3.52	Smooth
−0.1	0.45	14	0.83	3.798	Smooth	−0.1	0.45	15	0.612	4.219	Smooth
0.2	0.31	9	0.888	3.072	Smooth	0.2	0.33	9	0.847	3.523	Smooth
−0.2	0.27	10	0.732	4.112	Smooth	−0.2	0.26	9	0.626	4.067	Smooth
0.4	0.18	6	0.842	3.354	Smooth	0.4	0.19	6.5	0.877	3.253	Smooth
−0.4	0.14	5.5	0.67	4.24	Smooth	−0.4	0.13	5.4	0.672	4.084	Smooth

.

Its influence on the THG efficiency is weak enough for the positive phase mismatching between SH and FW if the absolute value of ${\Delta }_{31}k$ is less than $0.01$. However, when choosing the corresponding sign of the phase mismatching ${\Delta }_{31}k$, one can even increase the THG efficiency by $3\%$ (${\Delta }_{31}k=-0.01$).

With the growth of $|{\Delta }_{31}k|$ until 0.1, the energy conversion efficiency decreases, and this decreasing depends on the ${\Delta }_{31}k$ sign. Nevertheless, even in this case, about $45\%$ of the incident energy is converted to the TH energy. A sign of the SOD at FF does not remarkably change the conversion efficiency. However, the TH pulse spectrum width is greater if a medium possesses anomalous dispersion at the FF.

Comparing Table 4 and

Table 5. Influence of ${\Delta }_{31}k$ on the frequency tripling process at $\gamma =1,{\Delta }_{21}k=-20$.

Normal SOD at FF						Anomalous SOD at FF
${\Delta }_{31}k$	${\eta }_{max}$	$z_{max}$	${\tau }_3$	${\tau }_3^{\omega }$	TH pulse shape	${\Delta }_{31}k$	${\eta }_{max}$	$z_{max}$	${\tau }_3$	${\tau }_3^{\omega }$	TH pulse shape
0	0.67	30	4.925	5.924	Rather smooth	0	0.83	30	2.186	4.764	Not smooth
0.01	0.58	25	4.249	6.045	Smooth	0.01	0.86	40	2.551	5.423	Not smooth
−0.01	0.81	40	6.361	5.578	Rather smooth	−0.01	0.79	25	1.984	4.541	Rather smooth
0.1	0.32	11	1.866	4.08	Smooth	0.1	0.35	11	1.324	3.897	Smooth
−0.1	0.62	18	3.248	7.287	Rather smooth	−0.1	0.55	14	1.48	4.356	Smooth
0.2	0.2	7.5	1.422	3.803	Smooth	0.2	0.22	8	1.181	3.756	Smooth
−0.2	0.44	12	2.148	7.142	Smooth	−0.2	0.38	10	1.202	4.337	Smooth
0.4	0.11	5	1.13	3.698	Smooth	0.4	0.12	5	0.917	3.829	Smooth
−0.4	0.23	7.5	1.354	5.914	Smooth	−0.4	0.21	7	0.963	4.42	Smooth

, one can see the more pronounced influence of the phase mismatching between FW and TH on the THG efficiency if the phase mismatching between FW and SH is negative: ${\Delta }_{21}k<0$ (Table 5). In this case, one can differentiate two modes occurring for small $|{\Delta }_{31}k|$ ($0.01$ or smaller) and rather large $|{\Delta }_{31}k|$ values: ($0.1$ and larger). In the first case, a product of $D_1$ and ${\Delta }_{31}k$ must be positive to increase the THG efficiency. As follows from Table 5, at normal dispersion of a medium, the maximum conversion efficiency may be increased by $14\%$ by the appropriate ${\Delta }_{31}k$ in comparison with its corresponding value achieved at ${\Delta }_{31}k=0$. For the anomalous dispersion of a medium at FF, this conversion efficiency enhancement is not so big; however, it can be observed.

At larger values of the phase mismatching between FW and TH, the sign of ${\Delta }_{31}k$ strongly influences the THG efficiency regardless of normal or anomalous dispersion of the medium at FF. In both cases, the negative value of ${\Delta }_{31}k$ is more preferable than the positive one. For example, the maximum THG efficiencies are equal to $0.44$ and $0.38$ at ${\Delta }_{31}k=-0.2$, respectively, which are bigger than $0.32$ and $0.35$, reaching at ${\Delta }_{31}k=0.1$ (the seventh row of Table 5).

If phase mismatching between the FW and SH is much larger (for example, ${\Delta }_{21}k=127$), then acceptable values of the phase mismatching between the FW and TH ${\Delta }_{31}k$ become less. We see in

Table 6. Influence of ${\Delta }_{31}k$ on the frequency tripling process at $\gamma =1,{\Delta }_{21}k=127$.

Normal SOD at FF						Anomalous SOD at FF
${\Delta }_{31}k$	${\eta }_{max}$	$z_{max}$	${\tau }_3$	${\tau }_3^{\omega }$	TH pulse shape	${\Delta }_{31}k$	${\eta }_{max}$	$z_{max}$	${\tau }_3$	${\tau }_3^{\omega }$	TH pulse shape
0	0.69	200	4.291	1.007	Smooth	0	0.45	90	1.484	2.093	Smooth
0.01	0.51	105	2.162	1.306	Smooth	0.01	0.38	70	1.327	2.257	Smooth
−0.01	0.61	160	7.171	1.508	Smooth	−0.01	0.52	110	1.828	2.14	Smooth
0.1	0.08	28	1.169	2.251	Smooth	0.1	0.08	26	1.002	2.959	Smooth
−0.1	0.1	32	3.31	5.29	Smooth	−0.1	0.1	29	0.864	5.176	Smooth

that the difference between THG efficiency at ${\Delta }_{31}k=\pm 0.01$ and ${\Delta }_{31}k=\pm 0.1$ is very large. So, at ${\Delta }_{31}k=\pm 0.01$, the THG efficiency remains rather high. Moreover, it is even bigger at ${\Delta }_{31}k=-0.01$ than at ${\Delta }_{31}k=0$ if a medium possesses the anomalous dispersion at the FF.

Nevertheless, the phase mismatching ${\Delta }_{31}k$ must not be higher because already at ${\Delta }_{31}k=\pm 0.1$ the energy conversion efficiency becomes low: only about $10\%$ of the incident pulse energy is converted to the TH.

Thus, a very important conclusion for practice follows from Table 4, Table 5 and Table 6: one can achieve essential spectral broadening of the TH pulse by choosing the phase mismatching between TH and FW.

6. Nonequivalent Influence of Pulses’ GVM on the THG Efficiency under Cascading SHG

Here, we discuss how GVM between the interacting waves influences the TH conversion efficiency. We show that the GVM between FW and SH does not remarkably influence THG’s efficiency. However, it can change the shapes of interacting pulses. In contrast, the GVM between TH and FW can strongly decrease the generation process efficiency. Consequently, it is important to find a way to compensate (or, at least, decrease) the GVM influence in this case. One of the possible methods is proposed below.

6.1. GVM between the SH and FW (${\nu }_{21}\ne 0,{\nu }_{31}=0$): Weak Influence on the Frequency Conversion Efficiency

Let us suppose first that the FW and TH have the same group velocity, which means that the dimensionless parameter ${\nu }_{31}$ is equal to zero: ${\nu }_{31}=0$. We explore three different cases: SH has the same group velocity as the FW (${\nu }_{21}=0$). The next case is the following: the SH group velocity does not much differ from the first one; therefore, all three pulses propagate together (${\nu }_{21}=0.01$) along the propagation distance under consideration. The last case corresponds to large GVM between SH and FW (${\nu }_{21}=1$). Consequently, the pulse at doubled frequency leaves a time domain occupied by two other waves at the pulse propagation in a linear medium.

The computer simulation results computed for the parameters $\gamma =1,{\Delta }_{21}k=20,$${\Delta }_{31}k=0,D_1=-0.000032,D_2=-0.000083,D_3=-0.000199$ which correspond to the incident pulse duration about $1$ ps are presented in Figure 17. In other cases, they are proportionally changed (GVM parameters vary as a linear function of ${\tau }_P^{-1}$, SOD coefficients vary as ${\tau }_P^{-2}$).

One can see that there is no visible difference between the energy conversion efficiency occurring at ${\nu }_{21}=0$ and ${\nu }_{21}=0.01$ because the pulse center shift is small. If the GVM is large (${\nu }_{21}=1$) for the remarkable shift of the pulse centers on chosen propagation distance, then it slightly disturbs the pulse shapes at FF and tripled frequency. On the other hand, its influence on the conversion efficiency is positive because we do not observe any conversion efficiency oscillations in Figure 17i which are present in Figure 17g,h. At the same time, the maximum conversion efficiency remains the same (about $67\%$). It should be underlined that the spectra of FW and TH pulses are broadened in comparison with the incident pulse spectrum. This means that it is possible to compress the pulse with tripled frequency.

If the incident pulse duration is decreased until ${\tau }_P=100$ fs, for which the SOD coefficients are equal to $D_1=-0.0032,D_2=-0.0083,D_3=-0.0199$, then the GVM influence on the shape pulses is a little bit stronger (Figure 18).

As in the previous case, there is no visible GVM influence if GVM is weak. For large GVM (${\nu }_{21}=1$, right column), the pulse deformation can be seen for the right sub-pulse of FW. The THG maximum energy conversion is lower in Figure 18i than in Figure 18g,h, and there are not any energy conversion efficiency oscillations as in Figure 17i. However, the TH pulse spectrum is narrower in comparison with the previous case. With further a further decrease in the pulse duration (${\tau }_P=50$ fs), the interaction of waves remains the same. However, due to greater SOD influence, pulse shapes are nearly symmetric, even in the case ${\nu }_{21}=1$.

In the case of anomalous dispersion of a medium for the FF (for example, $D_1=0.0006,D_2=-0.0035,D_3=-0.0062$), the pulse shape deformations are nearly the same, as one can see in Figure 19.

In Figure 19a,b the TH pulse shape is smooth and even symmetric in contrast with one shown in Figure 17c, where the TH pulse shape deformation is seen, especially on the trailing part of the pulse. Nevertheless, from the point of maximum energy conversion efficiency, the difference between values shown in Figure 19d–f is negligible. Moreover, this efficiency is higher in comparison when its value reached that in Figure 18i. The TH spectrum bandwidth (Figure 19f) is greater than its corresponding value as shown in the previous Figures.

Thus, a large GVM between FW and SH wave can (and does) influence the pulse shapes, but in all cases, its influence on the energy conversion efficiency is weak. If there is an anomalous dispersion of a medium at the FF, then the TH pulse shape is symmetric, and its spectrum is broader than in the case of normal dispersion of a medium.

6.2. GVM between the TH and FW (${\nu }_{21}=0,{\nu }_{31}\ne 0$): Strong Influence on the Frequency Conversion Efficiency

In contrast to the previous case, the GVM between FW and TH ${\nu }_{31}$ can dramatically decrease the frequency tripling efficiency. Let us introduce the computer simulation results. The pulse shapes and the conversion efficiency at ${\Delta }_{31}k=0$ for the pulse incident duration 100 fs,$1$ ps$,10$ ps (${\nu }_{31}$ equals $1.02,0.102,0.0102$, respectively) of the FW are shown in Figure 20. Since we showed above that the GVM between FW and SH influences weakly on the frequency tripling process, we can choose ${\nu }_{21}=0$ without any restrictions. Below, the parameter $\gamma$ equals unity: $\gamma =1$.

Firstly, we present the results for a rather small phase mismatching ${\Delta }_{21}k$ ($|{\Delta }_{21}k|=20$ in Figure 20, Figure 21, Figure 22 and Figure 23) because, in this case, the distance required for the achievement of the maximum conversion efficiency is also quite small. Then, some results are discussed for higher phase mismatching between FW and SH (${\Delta }_{21}k=127$), which leads to more dramatic influence of the GVM. The parameters at time derivatives are equal to ${\nu }_{31}=1.02,$$D_1=-0.0032,D_2=-0.0083,D_3=-0.0199$ for the incident pulse duration 100 fs.

So, we see in Figure 20 that at ${\Delta }_{21}k=20$, the TH intensity is low at large GVMs between FW and TH in comparison to the incident pulse intensity (Figure 20a). The TH pulse duration is many times greater than the incident FW pulse duration. This means that we obtain at least the pulse with a picosecond duration if the incident pulse with a duration of 100 fs falls on the nonlinear medium. The THG energy grows with the propagation distance. However, the maximum intensity remains at the same level (about $0.025$).

In two other cases, the TH intensity and THG efficiency are high. In fact, at ${\tau }_P=1$ ps (Figure 20b,e,h), these characteristics of the wave interaction are even better than those achieving at ${\tau }_P=10$ ps (Figure 20c,f,i). The possible reason for this is the following: due to rather large GVM between FW and TH, the TH wave is generated and leaves a time domain of wave interaction. Therefore, there is no energy reverse conversion in contrast to Figure 20i, where the TH conversion efficiency oscillations are observed. On the other hand, the maximum energy conversion efficiency is reached at a smaller z-coordinate if the pulse duration is greater than in the previous case.

To clarify the GVM influence on the frequency tripling process, we provide a computer simulation for the coefficients ${\nu }_{31}$ ($1.02,0.102,0.0102$, respectively), but for fixed values of SOD coefficients, $D_1=-0.0032,D_2=-0.0083,D_3=-0.0199$ in all three cases. The computer simulation results are shown in Figure 21 (${\Delta }_{21}k=20$).

The GVM influences both the pulse shapes and the energy conversion efficiency. We see that in Figure 21b, the TH pulse shape is less smooth than in Figure 20b and the energy, converted to TH (Figure 21h), is about $5\%$ lower than the energy shown in Figure 20h. In fact, the lines shown in Figure 21c,f,i are very close to the ones in Figure 5, because the same distance of the maximum THG efficiency achievement ($z=13$), as well as its value (about $51\%$), occurs. A single difference is observed in the pulses shapes (Figure 21c). Thus, the GVM remarkably influences the pulse spectra: in Figure 21e the TH spectrum possesses local minimum in contrast to Figure 20e.

If the phase mismatching between SH and FW is negative (${\Delta }_{21}k=-20$), then the situation does not differ from (Figure 21) for ${\nu }_{31}=1.02$ and ${\nu }_{31}=0.0102$ (see Figure 22).

However, at ${\nu }_{31}=0.102$ the growing of the energy conversion efficiency finishes near the section $z=40$ due to appearing of MI. If the laser pulse with its incident duration ${\tau }_p=100$ fs propagates in a medium with anomalous dispersion at FF ($D_1=0.0006,$$D_2=-0.0035,D_3=-0.0062$, then at ${\Delta }_{21}k=20$ the THG efficiency (Figure 23) increases very slowly with growing of the propagation distance beginning from the certain section of a medium (Figure 23g). At the other two durations of the incident pulse under consideration, the frequency conversion process occurs similar to the case of normal dispersion of a medium at FF.

At larger phase mismatching between FW and SH (${\Delta }_{21}k=127$, Figure 24), if all three waves propagate in a medium with normal dispersion, the THG efficiency under large GVM is only about $2\%$ (Figure 24g, ${\tau }_P=100fs,{\nu }_{31}=1.02$). This happens because the GVM length is much shorter than the nonlinear coupling length and, therefore, the generation of a wave at tripled frequency is very weak. If the incident pulse duration is increased up to $1ps$ (${\nu }_{31}=0.102$), then the generated TH pulse duration becomes very long. Meanwhile, its maximum intensity is small, the THG efficiency is high enough and reaches about $33\%$. For longer incident pulse duration (10 ps, the right column in Figure 24), the dimensionless GVM coefficient is not very large (${\nu }_{31}=0.0102$) and the THG efficiency becomes higher (about $70\%$).

6.3. Compensation of GVM Action by Choosing Phase Mismatching and Shifting of FW Carrier Frequency

We see in the previous paragraph that the GVM between FW and TH may be the main physical mechanism restricting the THG efficiency. Therefore, we analyse briefly another (together with duration increasing) possible way to minimise (or remove at all) an influence of this GVM. Below we show that it is possible to implement by choosing the incident pulse phase distribution. Let us note that we do not consider the incident chirped pulse because it is of interest for another study. Our current interest is the influence of the frequency shift on the energy conversion efficiency in a medium with a quadratic nonlinear response. Previously, we have shown in [67] how it is possible to increase the energy conversion efficiency of the SHG process in a medium with the combined nonlinear response by choosing the appropriate incident frequency shift and the phase mismatching between FW and SH. Here, we provide a similar derivation for the THG process under the cascading SHG.

The incident complex amplitude distribution is chosen as follows:

$$A_1(0,t)=e^{-{(t-0.5L_t)}^2}{e}^{i\phi \left(t\right)},A_2(0,t)=A_3(0,t)=0.$$

(27)

We transform the equations in (1) by using the following substitution of the complex amplitudes:

$$A_j(z,t)={\tilde{A}}_j(z,t)e^{ij\phi \left(t\right)},j=1,2,3.$$

(28)

In fact, one can consider a more general case in which the frequency shifts differ from $j\phi \left(t\right)$, but this is outside of our current consideration.

It is easy to write a new set of equations:

$$\begin{array}{cc}\hfill & \frac{\partial {\tilde{A}}_1}{\partial z}-2D_1{\phi }^{\prime }\left(t\right)\frac{\partial {\tilde{A}}_1}{\partial t}+i{D}_1\frac{{\partial }^2{\tilde{A}}_1}{\partial t^2}+i\gamma \left({\tilde{A}}_1^{\ast }{\tilde{A}}_2{e}^{-i{\Delta }_{21}kz}+{\tilde{A}}_2^{\ast }{\tilde{A}}_3{e}^{-i({\Delta }_{31}k-{\Delta }_{21}k)z}\right)\hfill \\ \hfill & +i(-D_1{\left({\phi }^{\prime }\left(t\right)\right)}^2+i{D}_1{\phi }^{″}\left(t\right)){\tilde{A}}_1=0,\hfill \\ \hfill & \frac{\partial {\tilde{A}}_2}{\partial z}+({\nu }_{21}-4D_2{\phi }^{\prime }\left(t\right))\frac{\partial {\tilde{A}}_2}{\partial t}+i{D}_2\frac{{\partial }^2{\tilde{A}}_2}{\partial t^2}+i\gamma \left({\tilde{A}}_1^2{e}^{i{\Delta }_{21}kz}+2{\tilde{A}}_1^{\ast }{\tilde{A}}_3{e}^{-i({\Delta }_{31}k-{\Delta }_{21}k)z}\right)\hfill \\ \hfill & +i(-4D_2{\left({\phi }^{\prime }\left(t\right)\right)}^2+2i{D}_2{\phi }^{″}\left(t\right)+2{\nu }_{21}{\phi }^{\prime }\left(t\right)){\tilde{A}}_2=0,\hfill \\ \hfill & \frac{\partial {\tilde{A}}_3}{\partial z}+({\nu }_{31}-6D_3{\phi }^{\prime }\left(t\right))\frac{\partial {\tilde{A}}_3}{\partial t}+i{D}_3\frac{{\partial }^2{\tilde{A}}_3}{\partial t^2}+3i\gamma {\tilde{A}}_1{\tilde{A}}_2{e}^{i({\Delta }_{31}k-{\Delta }_{21}k)z}\hfill \\ \hfill & +i(-9D_3{\left({\phi }^{\prime }\left(t\right)\right)}^2+3i{D}_3{\phi }^{″}\left(t\right)+3{\nu }_{31}{\phi }^{\prime }\left(t\right)){\tilde{A}}_3=0,0<z\le L_z,0<t<L_t.\hfill \end{array}$$

(29)

Comparing the coefficients at the first derivative with respect to time in the first equation and the third one, we can choose the condition of synchronous motion of the pulses at FF and tripled frequency:

$$\phi \left(t\right)=\frac{{\nu }_{31}}{2(3D_3-D_1)}t.$$

(30)

Of course, inequality $D_1\ne 3D_3$ must be valid. Otherwise, it is impossible to compensate the GVM influence by choosing of frequency shift in the form (28).

If the equality (30) is valid, then the set of equations in (29) takes the form

$$\begin{array}{cc}\hfill & \frac{\partial {\tilde{A}}_1}{\partial z}-\frac{D_1{\nu }_{31}}{3D_3-D_1}\frac{\partial {\tilde{A}}_1}{\partial t}+i{D}_1\frac{{\partial }^2{\tilde{A}}_1}{\partial t^2}+i\gamma \left({\tilde{A}}_1^{\ast }{\tilde{A}}_2{e}^{-i{\Delta }_{21}kz}+{\tilde{A}}_2^{\ast }{\tilde{A}}_3{e}^{-i({\Delta }_{31}k-{\Delta }_{21}k)z}\right)\hfill \\ \hfill & -i\frac{D_1{\nu }_{31}^2}{4{(3D_3-D_1)}^2}{\tilde{A}}_1=0,\hfill \\ \hfill & \frac{\partial {\tilde{A}}_2}{\partial z}+({\nu }_{21}-\frac{2D_2{\nu }_{31}}{3D_3-D_1})\frac{\partial {\tilde{A}}_2}{\partial t}+i{D}_2\frac{{\partial }^2{\tilde{A}}_2}{\partial t^2}+i\gamma \left({\tilde{A}}_1^2{e}^{i{\Delta }_{21}kz}+2{\tilde{A}}_1^{\ast }{\tilde{A}}_3{e}^{-i({\Delta }_{31}k-{\Delta }_{21}k)z}\right)\hfill \\ \hfill & -i\left(-D_2\frac{{\nu }_{31}^2}{{(3D_3-D_1)}^2}+\frac{{\nu }_{21}{\nu }_{31}}{3D_3-D_1}\right){\tilde{A}}_2=0,\hfill \\ \hfill & \frac{\partial {\tilde{A}}_3}{\partial z}-\frac{D_1{\nu }_{31}}{3D_3-D_1}\frac{\partial {\tilde{A}}_3}{\partial t}+i{D}_3\frac{{\partial }^2{\tilde{A}}_3}{\partial t^2}+3i\gamma {\tilde{A}}_1{\tilde{A}}_2{e}^{i({\Delta }_{31}k-{\Delta }_{21}k)z}\hfill \\ \hfill & +3i\frac{{\nu }_{31}^2(3D_3-2D_1)}{4{(3D_3-D_1)}^2}{\tilde{A}}_3=0,0<z\le L_z,0<t<L_t,\hfill \end{array}$$

(31)

which means a presence of the additional phase mismatching between the interacting waves. Further, we make another substitution

$$\begin{array}{cc}& {\overline{A}}_1={\tilde{A}}_1{e}^{-i{D}_1\frac{{\nu }_{31}^2}{4{(3D_3-D_1)}^2}z},\hfill \\ & {\overline{A}}_2={\tilde{A}}_2{e}^{i\left(-D_2\frac{{\nu }_{31}^2}{{(3D_3-D_1)}^2}+\frac{{\nu }_{21}{\nu }_{31}}{3D_3-D_1}\right)z},\hfill \\ & {\overline{A}}_3={\tilde{A}}_3{e}^{3i\frac{{\nu }_{31}^2(3D_3-2D_1)}{4{(3D_3-D_1)}^2}z},\hfill \end{array}$$

and for new functions, we write the following set of equations:

$$\begin{array}{cc}\hfill & \frac{\partial {\overline{A}}_1}{\partial z}-\frac{2D_1{\nu }_{31}}{6D_3-2D_1}\frac{\partial {\overline{A}}_1}{\partial t}+i{D}_1\frac{{\partial }^2{\overline{A}}_1}{\partial t^2}+i\gamma \left({\overline{A}}_1^{\ast }{\overline{A}}_2{e}^{-i{\overline{\Delta }}_{21}kz}+{\overline{A}}_2^{\ast }{\overline{A}}_3{e}^{-i({\overline{\Delta }}_{31}k-{\overline{\Delta }}_{21}k)z}\right)=0,\hfill \\ \hfill & \frac{\partial {\overline{A}}_2}{\partial z}+({\nu }_{21}-\frac{4D_2{\nu }_{31}}{6D_3-2D_1})\frac{\partial {\overline{A}}_2}{\partial t}+i{D}_2\frac{{\partial }^2{\overline{A}}_2}{\partial t^2}+i\gamma \left({\overline{A}}_1^2{e}^{i{\overline{\Delta }}_{21}kz}+2{\overline{A}}_1^{\ast }{\overline{A}}_3{e}^{-i({\overline{\Delta }}_{31}k-{\overline{\Delta }}_{21}k)z}\right)=0,\hfill \\ \hfill & \frac{\partial {\overline{A}}_3}{\partial z}-\frac{2D_1{\nu }_{31}}{6D_3-2D_1}\frac{\partial {\overline{A}}_3}{\partial t}+i{D}_3\frac{{\partial }^2{\overline{A}}_3}{\partial t^2}+3i\gamma {\overline{A}}_1{\overline{A}}_2{e}^{i({\overline{\Delta }}_{31}k-{\overline{\Delta }}_{21}k)z}=0,0<z\le L_z,0<t<L_t,\hfill \end{array}$$

(32)

where

$$\begin{array}{cc}& {\overline{\Delta }}_{21}k={\Delta }_{21}k+D_1\frac{{\nu }_{31}^2}{2{(3D_3-D_1)}^2}-D_2\frac{{\nu }_{31}^2}{{(3D_3-D_1)}^2}+\frac{{\nu }_{21}{\nu }_{31}}{3D_3-D_1},\hfill \\ & {\overline{\Delta }}_{31}k={\Delta }_{31}k+\frac{3{\nu }_{31}^2}{4(3D_3-D_1)}.\hfill \end{array}$$

So, if we choose the phase mismatching between TH and FW in appropriate way:

$${\Delta }_{31}k=-\frac{3{\nu }_{31}^2}{4(3D_3-D_1)},$$

(33)

then the possible highest efficiency of the frequency tripling process is achieved. In turn, the phase mismatching as well as GVM between SH and FW also change according to the rules given by:

$$\begin{array}{cc}\hfill & {\overline{\nu }}_{21}={\nu }_{21}-\frac{4D_2{\nu }_{31}}{6D_3-2D_1},\hfill \\ \hfill & {\overline{\Delta }}_{21}k={\Delta }_{21}k+D_1\frac{{\nu }_{31}^2}{2{(3D_3-D_1)}^2}-D_2\frac{{\nu }_{31}^2}{{(3D_3-D_1)}^2}+\frac{{\nu }_{21}{\nu }_{31}}{3D_3-D_1}.\hfill \end{array}$$

(34)

We see from (33) that the phase mismatching ${\Delta }_{31}k$, which is necessary to choose for the GVM compensation, does not depend on the incident pulse duration and is determined only by the difference between group velocities of the waves, and their SOD coefficients. However, it depends on the normalization length $Z_n$ of the z-coordinate. The product of GVM and the length of normalization characterize the time shift between the centers of two pulses on this propagation distance.

Now, let us estimate the parameters for GVM compensation. For definiteness, we choose the incident pulse duration being equal to 100 fs for KDP crystal. So, SOD parameters are $D_1=0.0006,D_2=-0.0035,D_3=-0.0062$, and the GVM parameters are ${\nu }_{21}=0.67,{\nu }_{31}=1.02$. The phase mismatching between SH and FW equals 127. According to the provided analysis, the initial phase modulation must be chosen as follows:

$$\phi \left(t\right)=\frac{{\nu }_{31}}{6D_3-2D_1}t=-26.5625t.$$

To compensate the phase mismatching induced by GVM, we choose the phase mismatching between TH and FW as follows:

$${\Delta }_{31}k=-\frac{3{\nu }_{31}^2}{4(3D_3-D_1)}=40.640625.$$

For practice, it is important to compare this frequency shifting with FW carrier frequency. Let us write the complex amplitude in physical variables:

$${\overline{A}}_1(\overline{z},\overline{t})=\sqrt{I_{01}}{e}^{-\frac{{(\overline{t}-0.5{\overline{L}}_t)}^2}{{\tau }_P^2}}{e}^{-i(2\pi (\nu +\Delta \nu )\overline{t}-\overline{k}\overline{z})},$$

(35)

where ${\overline{A}}_1,\overline{z},\overline{t},{\tau }_P,\overline{k}$ are introduced in (3), and $\nu$ is the FW carrier frequency. For a 1064 nm wavelength, it is approximately equal to 300 THz. $\Delta \nu$ is the frequency shift. The dimensionless value of $2\pi \nu \overline{t}$ is defined as follows:

$$2\pi \nu \overline{t}=2\pi \nu {\tau }_{P}t\approx 186t.$$

which means that the frequency shift $\Delta \nu$ is about $14.2\%$ of the FW carrier frequency $\nu$, and may be reached in some cases, as we believe. It is also should be stressed that the ratio $\Delta \nu$:$\nu$ does not depend on the pulse duration ${\tau }_P$:

$$\frac{\Delta \nu }{\nu }=-\frac{\left(\frac{\partial \overline{k}}{\partial \overline{\omega }}{|}_{3\overline{\omega }}-\frac{\partial \overline{k}}{\partial \overline{\omega }}{|}_{\overline{\omega }}\right)}{4\pi \left(3\frac{{\partial }^2\overline{k}}{\partial {\overline{\omega }}^2}{|}_{3\overline{\omega }}-\frac{{\partial }^2\overline{k}}{\partial {\overline{\omega }}^2}{|}_{\overline{\omega }}\right)\nu }.$$

To verify the theoretical analysis, we compare the computer simulation results obtained by using the equation set (1) with the corresponding results obtained by using the equation set (32). For this aim, we also have to find the value of the modified GVM coefficient between SH and FW:

$${\nu }_{21}-\frac{4D_2{\nu }_{31}}{6D_3-2D_1}=0.298125.$$

Therefore, the value of the modified phase mismatching between SH and FW is chosen as follows:

$${\overline{\Delta }}_{21}k={\Delta }_{21}k-4D_2\frac{{\nu }_{31}^2}{{(6D_3-2D_1)}^2}+\frac{2{\nu }_{21}{\nu }_{31}}{6D_3-2D_1}+2D_1\frac{{\nu }_{31}^2}{{(6D_3-2D_1)}^2}=102.13085975.$$

The corresponding results are shown in Figure 25.

First of all, let us note that two solutions of the problem (1) coincide with each other. This means that our theoretical consideration is correct.

Secondly (and most importantly), the THG efficiency is high and reaches a value of about $48\%$, which is comparable with the efficiency obtained in Section 5 (for example, Figure 3).

Of course, in the general case, the necessary frequency shift of the FW may not always be reached in practice. To our knowledge, the wavelength shift of the FW may be about $20$ nm without losses of the conversion efficiency, if it is necessary to compensate the GVM influence. Therefore, we propose the following way to obtain the high-efficiency THG in one crystal. First of all, one can introduce the frequency shift. Then, by changing of an angle between the crystal axis and the pulse propagation direction, one may choose the optimal phase mismatching condition. If the obtained result is not enough, then one can also increase the pulse duration in such a way that the GVM influence becomes weaker. The TH pulse spectrum is wider than the incident one for the FW in many cases, as described above. So, finally, it is possible to compress the TH pulse and obtain a highly effective TH pulse with a short duration.

On the other hand, in some cases, the GVM’s influence on the THG process is positive, as we show above, and therefore, there is no necessity for its full compensation. Moreover, sometimes a nonzero value of the phase mismatching ${\Delta }_{31}k$, chosen by Formula (33), leads to an increase of the maximum THG efficiency, as can be seen in Figure 26.

7. Discussion

Let us discuss the advantages of the presented scheme of frequency tripling in comparison with other schemes.

First of all, a single nonlinear crystal is required for the effective frequency tripling at a phase matching of the process $\omega \left(o\right)+\omega \left(o\right)+\omega \left(o\right)=3\omega \left(e\right)$. Evidently, this setup is simpler than the setup with two nonlinear crystals and the setup using the QPM technique. In accordance with the computation for the crystals presented in Table 1, the phase-matched THG can be achieved for the wavelength of FW which is greater than 903 nm in the ADP crystal, 890 nm in the KDP crystal, 651 nm in the BBO crystal, and 552 nm in the RBBF crystal. Thus, the generation of optical harmonics with a wavelength lower than 200 nm is possible. Figure 2 and Table 4, Table 5 and Table 6 show that the effective THG can be achieved even if the phase matching condition is not achieved.

Secondly, because one requires phase matching for the THG and big phase mismatching for the SHG, the highly intensive SH pulse is not required. In all computer simulation results, the SH intensity does not exceed 2–3% of the incident FW intensity. Importantly, the GVM between FW and SH does not remarkably influence the THG process (Figure 17, Figure 18 and Figure 19).

Thirdly, the high-efficiency THG can be observed even for rather low intensities (about $1$ GW/cm${}^2$ and low) of the incident wave at FF. Hence, pulses with picosecond or even nanosecond durations can be used to minimize the influence of the GVM and the SOD on the frequency tripling. Importantly, the frequency tripling accompanies the TH spectrum broadening (Figure 6, Figure 7, Figure 8 and Figure 17), which is absent for the well-known scheme of cascaded THG.

The sign of effective cubic nonlinearity (see the set of equations (9)) is defined by the phase mismatching ${\Delta }_{21}k$, which can be useful in practice. In particular, it allows self-focusing of the pulses even in a medium with a normal SOD.

The THG efficiency equal to $94.5\%$ for the long pulses (Figure 3) and more than $80\%$ for the short pulses (Figure 8, Figure 13, Figure 15 and Figure 16, Table 3, Table 4 and Table 5), reached theoretically, is greater than the corresponding value for the scheme based on the QPM technique and the scheme with changing of the SH polarization, and is comparable with the scheme using two crystals.

It should be emphasized that the length of a crystal in the scheme under consideration should be large—at least 1 cm, as the theoretical estimation showed. That is why in the above mentioned papers which account for cascading SHG, the researchers could not observe a high-efficiency frequency tripling when using the frequency conversion process proposed here because the length of their crystals did not exceed several millimeters.

The predicted value of $94.5$ conversion efficiency was computed without accounting for the absorption of the crystal. The absorption coefficient may change from quite a high value (0.1 cm${}^{-1}$) for the ADP crystal and (0.03 cm${}^{-1}$) for the KDP crystal to a very low value (<0.001 cm${}^{-1}$) for the BBO crystal for the optical pulse with a wavelength ${\lambda }_1=$ 1064 nm [68]. At ${\lambda }_3=$ 353 nm, the absorption coefficient of the KDP crystal is lower at 0.003 cm${}^{-1}$, and of the KBBF crystal is 0.1 cm${}^{-1}$ at the optical pulse wavelength 330 nm [69]. Therefore, the absorption can be neglected in some cases.

It is important to notice that the proposed scheme for optical frequency conversion can be used for other frequency conversion problems, such as fifth-harmonic generation [70] or frequency down-conversion [71].

8. Conclusions

We proposed and described a new physical scheme of THG based on cascading SHG, which leads to observing of the nonlinear processes inherent to a cubic nonlinear response. Physically, it is based on the consecutive generation of SH and a wave with sum a frequency in the same crystal with quadratic susceptibility under large phase mismatching between the FW and the SH wave. This leads to a possibility of the high-efficiency frequency tripling of optical radiation with the incident power density being enough for the quadratic nonlinear response of a medium to occur. This value can be essentially less than that requiring for the THG in a medium with cubic susceptibility (according to Table 1). The crystal length at which the efficient THG occurs depends on the phase mismatching between the SH and FW. It varies from 2 cm to 4 cm for a pulse with a long duration (Figure 3) and from 1 cm to 20 cm for a pulse with a short duration (Figure 13 and Figure 15). The required crystal length decreases with increasing power density of the incident pulse at FF. The pulse self-modulation at the tripled frequency is absent. Its cross-modulation occurs and, consequently, the TH spectrum is broadened. Therefore, this pulse can further be compressed. These two features distinguish the THG under consideration from cascaded THG through the SHG and SFG as well as from the THG in a medium with cubic susceptibility.

The THG conversion efficiency can reach $94\%$ (Figure 3) at the optimal parameters of wave interaction. The high-efficiency THG is reached for both signs of the phase mismatching between FW and SH. The THG can demonstrate a bistable mode under certain conditions.

The frequency tripling process via cascading SHG was investigated theoretically. On the basis of the multi-scale method, the modified set of equations was derived. Using its integrals of motion, an analytical solution was developed in the framework of the long-pulse duration approximation. This solution approximates the solution of the Schrödinger equations (1) with high accuracy during the pulse propagation distance, which is many times greater than the distance occurring using the well-known wave energy non-depletion approximation (Figure 3 and Figure 4). The modified equations, as well as the analytical solution, were confirmed by the computer simulation results obtained on the basis of the equations in (1).

We studied the influence of other parameters on the efficiency of the frequency tripling, SOD, and the phase mismatching between TH, FW, and GVM between both SH and FW as well as TH and FW. The computer simulation results of frequency tripling accounting for the SOD of the pulses showed the possibility of reaching high THG efficiency for both positive (Figure 5) and negative (Figure 6) phase mismatching between SH and FW. At its negative sign, the MI may appear, and as a result of this, a dramatic decreasing of the THG’s efficiency occurs (Figure 11). However, at the beginning stage of the MI, a regular sequence of very short sub-pulses with self-similar shapes may appear (Figure 12).

The high THG efficiency can be achieved for both types of a medium: anomalous (Figure 8) and normal (Figure 6) dispersion at the FF.

The GVM between SH and FW pulses only influences the pulse shapes without strong changing of the THG efficiency (Figure 17, Figure 18 and Figure 19). In contrast, the GVM between TH and FW remarkably influences the frequency conversion process (Figure 20, Figure 21, Figure 22, Figure 23 and Figure 24). To compensate for this negative influence, we proposed and discussed some physical schemes for three-wave interaction (Section 6.3). Moreover, since only quadratic nonlinearity is required for THG through the cascading SHG, it is possible to use pulses with durations of picoseconds or nanoseconds. It allows us to decrease a negative influence of the GVM and SOD. Another way is using so-called “tilted pulses” [72,73].