Mathematical Modeling of Soliton-Like Modes at Optical Rectification

We discuss the results of numerical modeling of forming optical-terahertz bullets at the process of optical rectification. Our calculations are based on a generalization of the well-known Yajima - Oikawa system, which describes the nonlinear interaction of short (optical) and long (terahertz) waves. The generalization relates to situations when the optical component is close to a few-cycle pulse. We study the influence of the number of optical pulse oscillations on the formation of an optical-terahertz bullet. We develop original nonlinear conservative pseudo-spectral difference scheme approximating the generalization of the Yajima-Oikawa system. It is realized with the help of FFT algorithm. Mathematical modeling demonstrates scheme efficiency.


Introduction
Studies in the ways of terahertz radiation generation are of increasing interest. Signals of the THz range today find many applications in image processing, security systems, astronomy, biology, medicine and other fields [1,2]. One of the most efficient methods used to generate THz radiation is a method based on the optical rectification mechanism [3,4,5,6,7,8,9]. In this case, the spectrum of the THz signal is broadband, which means that the spectral width of the pulse is comparable with the central frequency of its spectrum [9]. Generated pulse contains about one oscillation period of THz range. Thus, in the optical THz generation method, the signal has the properties of extremely short pulses. Consequently, in a theoretical consideration of the interaction of such pulses the well-known approximation of slowly varying envelope is not further applicable.
First the idea of generating of a terahertz pulse based on the optical rectification mechanism was theoretically proposed in [3]. Experimental studies confirmed this idea [4,5].
Realizing optical rectification, one launch a femtosecond optical pulse to a nonlinear medium, at that, the pulse spectrum contains frequencies at the difference of which THz radiation can be generated. The generation condition can be obtained from the laws of conservation of momentum and energy for elementary acts of scattering. They lead to the Cherenkov condition where v g is the group velocity of the optical pulse, v T is the phase velocity of the THz signal, Θ is the angle between the directions of the propagation of the optical and the THz pulses [3,4,5]. In the collinear mode (Θ = 0) the phase matching condition has the form of the Zakharov -Benney resonance [6,9]: v g (ω) = v T (Ω). Here ω is the optical signal frequancy,Ω is the frequency of the THz component.
Applying the approach of the slowly varying envelope to the optical component and the unidirectional propagation approximation [10] to the terahertz component one have a set of two nonlinear equations [6,7]. A similar system of equations was obtained in [11] for interaction of ion-sound and Langmuir waves in a plasma. From mathematical point of view, this system turns out to be integrable by means of the inverse scattering transform method. It is often called the Yajima-Oikawa system. Yajima-Oikawa system for optical and terahertz components holds at the Zakharov-Benney resonance, when diffraction and nonlinearity in THz component are absent. A soliton solution of this system has been found [6,7]. At that, it has been obtained that in the process of optical rectification the light component experiences a nonlinear "red" shift of the carrier frequency which is proportional to optical radiation intensity. This effect has been demonstrated in experiment [12]. In [13] nonlinearity in THz component and the phase modulation of the optical pulse were taken into account. Integrability of an arising generalized Yajima-Oikawa system has been shown. Optical pulsed radiation with a tilted wavefront has been shown useful to increase (by two-three orders of magnitude) the generation efficiency for media with velocities, which are significantly different in optical and terahertz regions [8,14]. In [15] the optical rectification in a system of resonant anisotropic molecules with a permanent dipole moment was studied taking into account the transverse dynamics of pulses. It has been shown that, if the carrier frequency of a femtosecond optical signal is higher than the eigen frequency of a selected molecular transition by a certain value, the generation occurs in the filamentation regime followed by transformation of these filaments into stable optical terahertz bullets. Authors of [16] have demonstrated that in the generation process of terahertz radiation a filamentation of the optical and terahertz components arises but no stable optical terahertz spatiotemporal solitons are formed. Capturing of optical and terahertz pulses in a waveguide has been considered in [17] and it has been shown that capture occurs under conditions of filamentation with both normal and anomalous group dispersions of the optical signal. In [18] the formation of an optical-terahertz spatiotemporal soliton was studied in a focusing gradient waveguide in the case of few-cycle optical pulses.
Numerical simulation of the effects at optical rectification is based on generalizations of the Yajima-Oikawa system. They have motion integrals, thus, we must develop conservative numerical methods holding difference analogues of the integrals [19,20]. Generalizations of the Yajima-Oikawa system are multidimensional and we should use an effective method saving computational time. In [15,16,17] the numerical simulation of THz pulse generation was performed using Runge-Kutta integration scheme on the propagation coordinate z. Often the splitting technique is applied to such problems [19,21,22]. These two classes of methods are remarcable because of a significant reduction in calculation time. But they are not suitable for the construction of methods preserving several motion integrals [23]. It has been demonstrated that Fast Fourier Transform (FFT) methods are conservative and effective being applied to multidimensional problems of nonlinear optics [24,25,26].
In the present paper we propose a conservative pseudo-spectral difference scheme realized with FFT algorithm. This original scheme is applied to a study of a few-cycle optical pulse generating an optical-terahertz bullet in a waveguide. We show that using a few-cycle pulse launched into a waveguide we can promote the capture of the terahertz component into a parametric bullet.

Model and numerical method
We apply the unidirectional propagation approximation [10] for the optical component and take into account the focusing optical and terahertz waveguides. Generalized system for the optical ψ and THz E T components is as follows χ and χ (2) are the linear and nonlinear optical susceptibilities, respectively, n ω,T (x) are the optical and THz refractive indices, n ω,T are these indices near the waveguide centre. At that, near the waveguide centre the following conditions hold: g ω,T (x) = f ω,T (x) = 0. The limit of a quasi-monochromatic optical pulse in the system (1) corresponds to the conditions In this limit we neglect the second terms in the brackets of the first equation (1) as well. If we also neglect the transversal and waveguide effects (∂ 2 /∂x 2 ≡ g ω ≡ g T ≡ 0), then the system (1) goes into the equations of Yajima-Oykawa. The third, fifth, and sixth terms on the right-hand side of the first equation of the system (1) appear in [13]. The first, second and fourth terms on the right-hand side of the second equation of the system (1) are included in the generalized Yajima-Oykawa system in the same work. The fourth term iq ∂ ∂τ ψ * ∂ψ ∂τ − ψ ∂ψ * ∂τ corresponds to taking into account the phase modulation of the optical pulse when generating terahertz radiation [13]. In addition, in [18] the waveguide and diffraction corrections to the quasi-monochromatic limit are taken into account. These corrections correspond to the term i ω ∂ ∂τ in both brackets of the first equation of the system (1). We perform numerical simulation normalizing the system (1) in the following way. ψ =ψψ 0 , E T =Ē T ψ 0 , ψ 0 is the peak amplitude of the optical component at the entrance of the medium.
Diffraction length of the optical pulse l ω D = nωω c R 2 0 , diffraction length of the THz signal l T D = 2n T cτ 0 R 2 0 . R 0 and τ 0 are the width and duration of the initial optical pulse.
Since pulse shortening may influence greatly all considered processes, it is reasonable to rewrite the above dimensionless parameters through the number N of oscillations under optical pulse carrier. We use the ratio corresponding to the soliton solution [13]: q σ = b a , µ = 2b. As the result, we can show that all processes caused by the optical pulse few-cycle nature depend on N , thus, some coefficients may be expressed as follows: N . Then, the system (1) can be transformed: Only optical pulse is launched into the input of the nonlinear crystal.
Boundary conditions are zero.
In (2-4) L z is the normalized length of nonlinear medium, L τ is the dimensionless time interval during which laser pulse interaction with a medium is analyzed. L x is the normalized length of the transversal domain. Below we omit bars in dimensionless variables' notations. System (2-4) possesses the following integrals of motion: In (5) parameter N 0 is the number of photons in the optical pulse. Each photon gives a part of its energy to the terahertz region. At that, the number of terahertz photons increases and the frequency of optical photons decreases. The latter means the decrease of the optical pulse carrier frequency. These effects may be presented as parametric scattering or optical photon decay into a different optical photon of less frequency and a terahertz signal. Constructing a finite difference scheme for the problem (2-4) we will keep in mind the fulfillment of difference analogs of (5).
Nonlinear symmetric finite difference scheme for Fourier coefficients corresponding to (2-4) is given below. Since the formulae are too cumbersome in what follows we omit the terms of the orders of N −1 and N −2 .
In (7) 0.5 A (z j ) = (A(z j+1 ) + A(z j )/2, λ m = (2πm/L x ) 2 , λ n = (2πn/L τ ) 2 , where A(z j ) is a grid function. Scheme (7) is realized with iterations: −i( Scheme (7) has the approximation of the order of O(h 2 z + h mx x + h mτ τ ), where mx and mτ are the orders of the highest derivatives of the functions ψ, E T with respect to x and τ correspondingly. Scheme holds difference analogs of (5) where integrals are replaced by sums.

Results of numerical modelling
We have performed a series of numerical experiments on the basis of the dimensionless system (2-4) using the numerical scheme (6)(7)(8). Varying the value of third-order dispersion we studied the influence of N on the process of optical rectification. Our goal was to simulate the formation and propagation of a coupled optical-terahertz beam-pulse. The brightest illustration refers to the case of normal dispersion, focusing optical and terahertz waveguides and very short pulse duration. The corresponding coefficients are expressed as follows: N = 5, D k2 = 0.1, D k3 = 0.02, D σ = 0.2, D gω = −1.0, D gT = −0.1, D xω = D xT = 0.1. It is not correct to call this localized state as a soliton but nevertheless it propagates as a whole up to the distance more than ten diffraction lengths, gradually scattering its energy. In Fig. 1, a-c the transversal profiles are given of the optical and terahertz components at different distances up to z = 200. Beam profile is similar to the soliton of the 1D integrable Yajima-Oikawa system. Further the profile is undergoing distortions and the beam-pulse decays. In Fig.1, d the evolution along z-axis is presented of the peak intensities of both optical and terahertz components. We observe beam decay starting from the drastic decay of the terahertz component. After that the optical component loses coupling.

Conclusion
In the present paper we consider optical rectification at a few-cycle input pulse. Our prime interest is to study the capture of the terahertz component into a parametric bullet. Ultra-short pulse duration requires the account of higher order effects, so additional terms are added to the Yajima-Oikawa system. To solve this complicated problem, we develop an original nonlinear conservative pseudo-spectral difference scheme. Mathematical modeling has confirmed scheme efficiency and conservativeness. We demonstrate the formation of a soliton-like bound state propagating over a distance of more than 10 diffraction lengths. Moreover, we show that formation possibility increases when a number of oscillations N decreases. This is explained by the fact that for soliton capture at quadratic nonlinearity and normal dispersion a special medium geometry is required. The simplest way is to use a waveguide. Concerning opticalterahertz interaction it is very important to balance the contributions of diffraction for optical and terahertz components. In the quasi-monochromatic case (N = 100) the influence of a terahertz waveguide is higher than the influence of an optical one and two-component bound state does not form. In the few-cycle case (N=5) the contributions of optical and terahertz diffraction terms are comparable with each other. As a consequence, we observe the formation and further propagation of an optical-terahertz bullet. Work was supported by Russian Science Foundation (Project No. 17-11-01157).