Validity of the unidirectional propagation model: application to laser-driven terahertz emission

The Unidirectional Pulse Propagation Equation (uppe) is often used to compute the forward component of ultrashort light pulses in nonlinear materials. While its accuracy has frequently been reported for many applications in ultrafast optics, its validity can be questionable for ionizing pulses, in particular in the frequency domain below the electron plasma frequency ω pe . Inaccuracy of the uppe model in this frequency range may be detrimental to, e.g., a correct description of laser-driven terahertz emission. Here we demonstrate, analytically and numerically, that over long enough propagation paths, the one-dimensional solutions of both uppe and the full wave equation match in the entire spectrum, including the spectral range ω ≤ ω pe . Our findings confirm the reliability of the uppe solutions for a wide class of propagation problems.

which involves the plasma current density J z t , ( ) and the macroscopic polarization P z t , ( ). The latter quantity includes the nonlinear polarization, P E NL 0 3 3  c = ( ) , with 3 c ( ) representing the third-order susceptibility assumed to be instantaneous. 0  and 0 m are the electric permittivity and the magnetic permeability in vacuum, respectively. For simplicity, we shall neglect linear gas dispersion in the forthcoming analysis. At moderate laser intensities and on short time scales, ion motion can be neglected. In the non-relativistic limit the electron current density reads as where W E (| |) is a field-dependent photoionization rate and N a is the initial density of neutral atoms. Generally, solving full Maxwell models accounting for multiple ionization processes as well as optical and plasma nonlinear effects is computationally expensive when simulating long propagation ranges as, e.g., in the for an integrable function f and k c w = . Equation (5) is used in ultrafast nonlinear optics to describe light pulses with ultrabroad spectra [5] and in higher dimensions it can describe pulse diffraction associated with transverse wave numbers k ⊥ through the formal change k k k 2 2 2  -^ [6].
The solution to the partial differential equation (1) is commonly advanced in time from a couple of initial/ boundary conditions on E and its derivative. In contrast, equation (5) only requires the incident field value E z t 0, = ( )and it is solved along the longitudinal direction z 0 > . Equation (5) may be refound from an easy factorization of equation (1), expressed in Fourier domain as [7][8][9] k k E Q i i . 7 Assuming weak enough nonlinearities, the field Ẽ is decomposed as the sum of (quasi) linear forward (+) and backward (−) modes, E E E e e kz kz i i = + + --˜˜. Substituting the previous quantities into equation (7), applying the approximation on the backward propagator: k E kE i e 2i e z kz kz i i ¶ + » + + ( )˜˜, and conjecturing that both Ẽ and Q contain essentially the part corresponding to E + and Q + , we then arrive to equation (5). In this simple derivation, the amplitude E +  is supposed to be slow compared to e kz i . This assumption is not necessary if more involved projection techniques are employed and if the PDE (1) is treated as an inhomogeneous second-order ODE [3,4]. A pending issue about the solutions of equations (5) and (7) is thus their equivalence in the whole spectrum under the formal substitution k k i 2i z ¶ +  . This equivalence is implicitly assumed a priori when one uses a UPPE model. Proving this statement a posteriori, by computing explicit solutions of these two models, remains to be done.
Over the past decade, the increasing potential of terahertz (THz) radiation has stimulated intensive efforts to develop efficient sources based on the ionization of gases by ultrashort laser pulses [10,11]. In order to properly calibrate dedicated experiments, nonlinear propagation codes have to be particularly accurate in assessing the lowfrequency part of the pulse spectrum. Different experiments and plasma conformations (e.g., micro-plasmas) have been recently examined, for which relevant THz emissions in frequency ranges below the electron plasma frequency were reported [12][13][14]. So far, however, the validity of the UPPE approach has received no numerical confirmation in this specific spectral range, which is often accessed in laser-driven THz experiments [15].
To highlight more this problematics, let us consider for a moment the elementary situation where the electron density is constant, as met in a preformed homogeneous plasma. Usually, solving equation (1) needs both initial and boundary conditions on E z t , ( ), the evolution of which is monitored by the dielectric function that characterizes the mediumʼs response to the electric field. In such a plasma with negligible collisions, the dielectric function expresses as 1 )is the electron plasma frequency and ω is any Fourier frequency of the electromagnetic wave. When a monochromatic wave with frequency ω is incident on a vacuum-plasma boundary, only the fraction with real wave number k c 1 T pe [16]. Frequency components with pe  w w are evanescent and vanish over the plasma skin depth. So, for Fourier components pe  w w , the incident pulse is mainly reflected and we expect this spectral range to damp out in the transmitted pulse along the propagation path, consistently. In contrast, for the UPPE model, plasma opacity is not present. This can be better appreciated when looking at the linear modes of the two equations for N e = const., namely, for equation (1), and for equation (5). Obviously these two linear solutions do not match, in particular in the domain pe  w w . The WE solution (8) vanishes exponentially as z exceeds about one plasma skin length c pe pe d w º [17]. Accordingly, the UPPE solution (9) can in principle hold in the frequency range pe w w  only. In current applications, however, in which pe w varies through the time variations of the electron density N t e ( ), UPPE is supposed to be valid in the entire frequency range.
Whereas the solving for equation (1) can be performed exactly for a constant plasma frequency [16], the same problem becomes tricky when Q is a nonlinear function of E. This situation applies to a number of important issues in nonlinear optics and plasma physics [18,19]. Here, we report that, surprisingly, both solutions to equations (1) and (5) become generically identical when the plasma is self-generated in time by the laser pulse. In this paper we indeed show that the THz pulse spectra and fields, when they are described either by a full Maxwell-fluid model-encompassing both backward and forward propagations-or by its UPPE approximation-modeling only the forward wave-match as the pulse propagates over distances larger than the plasma skin depth and its dynamics is mainly driven by the nonlinearities.
The paper is organized as follows: we derive in section 2 exact analytical solutions for the two WE and UPPE 1D models describing the secondary fields produced inside a collisionless plasma. In the framework of a perturbative approach, WE and UPPE solutions are shown to converge for long enough propagation distances. In section 3 we display numerical evidence by means of our 1D UPPE and MAXFLU codes that over distances z pe d > , both WE forward and UPPE solutions indeed match in the THz spectral range and therefore provide similar THz fields. This statement is established for two-color femtosecond pulses and few-cycle single-color pulses as well, at intensities privileging Kerr or plasma responses. Section 4 concludes our investigation. Appendix describes the numerical implementation of the UPPE and MAXFLU codes.

The analytical solution
where I 0 is the mean pump intensity, L t is the full-width-at-half-maximum (FWHM) duration of the fundamental pulse (the second harmonic has half pump duration), r denotes the relative intensity ratio of the second harmonic and f is the relative phase between the two harmonics. One chooses a fundamental wavelength of 1 0 l = μm related to the carrier wave frequency c 2 2 0 0 0 w pn p l º = , and an intensity ratio of r 0.1 = (two-color pulses) or 0 (single-color pulses). The interaction gas is argon at ambient pressure.The ionization rate W E (| |) is given by the well-known quasi-static tunnel (QST) rate where the constants α and β depend on the ionization potential of the gas [20,21]. More complex rates describing multiple ionization events could be implemented as well [22]. These will be employed when investigating the highest intensity levels.
Because of the complex nonlinearities involved, we apply the same perturbation approach as in [23]. We decompose E into an unperturbed laser field, E L , and an induced perturbation (secondary field), E d , namely, so that the nonlinearity Q reduces to an inhomogeneous source term only evaluated on the laser field. Using equation (1) one obtains Here, the perturbation E d contains the THz wave field, which can be extracted by means of a suitable frequency filtering. The electron plasma frequency pe w is a function of time through its dependency on N z t , e ( ). In the tunnel ionization regime, the electron density increases steplike and produces THz radiation by coupling with the high-frequency field E z ct L -( ) [6]. In the framework of UPPE (equation (5)), the perturbation instead obeys Equations (13) and (15) respectively. It is important to notice that E L is only a function of ξ and thereby N e and pe 2 w , computed on the same laser field, only depend on ξ too. Our initial conditions for the WE read yielding no constraint on the spatial derivatives. The pulse is sited in the half-plane z ct < ( 0 x < for causality reasons) and enters the plasma region at time t=0. To satisfy the requirements (18), the input pulse (z=0) is positioned in the plasma domain at and (e)). Such boundary conditions offer flexibility to treat both UPPE and WE solutions in the same analytic framework. However, they differ from the inputs used in the numerical scheme of the UPPE1D and MAXFLU1D codes discussed in appendix. Therefore, we shall have to propagate over long enough distances c 2 L  t , ensuring the full development of the nonlinearities, in order to perform quantitative comparisons between the perturbed fields E d and their numerical counterparts. These boundary conditions also differ from those used in [23], where the progression of the laser pulse towards the plasma was described and the source term Q was tuned to zero for z 0 < in this reference.
For the WE, we now apply the Laplace transform f p f s s e d ps 0 After some algebra we find straightforwardly where J 1 is the Bessel function of the first kind and ). In the original frame variables, this solution expresses as and its spectrum is obtained after taking Fourier transform in time t.
Applying the same treatment to equation (17), we get From a first glance, equations (22) and (24) only differ by their linear kernel As the UPPE/WE pulses advance in the z t , ( ) plane, the integration variable x¢ can cover all values from 0 to ct 0 min max x ¢ = -, t max being fixed by the boundary of the time window (t 3.3 max = ps here). Both UPPE and WE solutions converge towards each other at coordinates satisfying z z ct This requirement can indeed be understood from applying simple Taylor expansions in the argument z ct z z ct in the UPPE solution in the same limit. These behaviors show that the UPPE solution is missing a dependency in the z ct + ( ) propagation variable associated with a backward component. Applying the second inequality z z ct - | |in the previous approximations leads to the convergence of the two solutions.
This property is reflected by figure 1(a), which plots the argument of the Bessel function for a two-color 50 fs Gaussian pump pulse with 150 TW cm −2 intensity. Here, the pulse region extends over 100 x¢ < | | μm and we can observe the convergence of the two Bessel arguments when the distance z is increased from 100 to 500 μm. Concerning their potential discrepancies, one can observe that, at large times ( 100 x¢ > | | μm), the plasma frequency pe w is constant, so that G z ct z ct , ) . The oscillations in the WE/UPPE solutions are then dictated by those of the Bessel function J 1 . At times corresponding to 600 x¢ »μm and z=100 μm, these oscillations relax on the plasma period as J X t t c sin  Figure 1(b) thus displays evidence of a minimum frequency smaller than pe w in the UPPE solution, whose value increases with z until reaching the electron plasma frequency.
To get more insight into the convergence dynamics, we may also consider a much simpler situation by assuming a constant inhomogeneity Q. With the help of equations (22) and (24), this yields the ratio which can be useful to understand the differences induced by the linear propagators. Figure 1(c) shows equation (27) as a function of z and t. The lower-right part z ct > ( ), in which our solutions are not defined, is set to zero for causality reason. Both WE and UPPE solutions converge as long as z z ct - | |, which includes the laser region. In contrast, the same solutions depart from each other in the sharper limit ct z  . Let us now imagine that, in the vicinity of the pulse head, the source term Q has a certain finite extent, schematically delimited by the two grey solid lines in figure 1(c). The WE and UPPE solutions converge near the laser head, where they are dominated by the nonlinearities computed on the laser profile. In the opposite domain, the solutions are mainly driven by their linear propagators that behave differently over large times. However, the larger the propagated distance z, the broader the convergence domain, which spans a cone in the z t , ( ) plane (see blue area in figure 1(c)). This is confirmed by figures 1(d) and (e) that detail the field amplitudes in the z t , ( ) plane computed from the complete expressions (22) and (24). One can see that the UPPE field contours differ from the WE contours in the spatio-temporal domain ct z  . In particular, a hyperbolic distribution occurs, associated with the longer periods of figure 1(b) and with the fact that UPPE does not admit pulse components varying with z+ct. Nevertheless, the solutions achieve the same dominant component near the nonlinearity region, where they mutually converge and whose area grows with z (see figures 1(d) and (e) where convergence is reached inside the cones delimited by black dashed lines). Figures 2(a) and (b) show some examples of analytical UPPE/WE spectra and fields for a two-color Gaussian pulse with FWHM duration of 50 fs, 150 TW cm −2 overall intensity. The nonlinearities consist of plasma generation alone. Here and in the following the THz fields shown as insets are computed from an inverse Fourier transform of E d˜in the frequency window 2 90  n w p º THz. The plotted propagation distances are z=100 μm and z=1 mm. One can observe that (i) the spectral region 2 pe pe n n w p < º becomes depleted as z increases, (ii) the minimum frequency marking the UPPE spectrum, min n , increases in turn, and (iii) in the rear part of the pulse (beyond the laser head) the UPPE linear mode consistently develops longer periods (see inset of figure 2(a)). At smaller intensity (I 50 0 = TW cm −2 ) and with zero phase angle ( 0 f = ), the Kerr response is expected to be a key player [24]; so we now include it. A typical Kerr-plasma spectrum at weak intensity is illustrated in figure 2(c), where the low plasma frequency, 0.53 pe n = THz, highlights the lesser contribution of photoionization and a parabolic spectral shape characterizes the Kerr signature in the band of higher THz frequencies ( 10 20 n » -THz) [11]. For a mean pump intensity I 1 0 = PW cm −2 , in contrast, several electrons can be extracted from their atom. In this high intensity regime, the ionization of the successive electron shells of argon is described by the multiple-ionization model built in [22] and based on a field-dependent cycle-averaged rate computed from the Perelomov-Popov-Terent'ev (PPT) theory [25]. At short propagation distances, unlike the WE solution, the UPPE spectrum is not peaked at the plasma frequency 55 pe n » THz, but it develops spectral oscillations in the region pe n n < , as expected ( figure 2(d)). At larger distances, discrepancies are amplified (not shown), because E d starts to break the underlying hypothesis of our perturbative approach, E E L d  , as the perturbation itself produces optical frequencies , 2 , ... 0 0 w w ( ) through the nonlinearities, e.g., the photoionization. Since the optical pump pulse is not depleted along propagation, our formalism cannot assure a proper conservation of the electromagnetic energy. As a validity criterion we consider that our analytical solution stops to hold whenever the spectral intensity of E d at optical frequencies becomes comparable to 75% of the laser spectral intensity. Such limitations are of course absent in the results of the full numerical simulations, as can be inferred from, e.g., figure 6.

Numerical results
Our theoretical expectations are tested by running the MAXFLU1D and UPPE1D codes, whose respective numerical schemes are detailed in appendix. For both codes the input condition at z=0 is the two-color Gaussian pulse defined above by equation (10). The pump intensity is alternatively set to 50, 150 and 1000 TW cm −2 in order to investigate various ionization degrees. The phase angle f is valued as 0 f = to enhance the Kerr effect in the 50 TW cm −2 case and 2 p otherwise [24]. For Gaussian pulses with moderate laser intensity, from 50 to 150 TW cm −2 , and undergoing single ionization, we employ the QST rate (11). Consistently with section 2, when dealing with 1 PW cm −2 pulses, multiple ionization will be described from the multi-ion model of [22] employing the field-dependent PPT ionization rate.
To start with, only the Kerr response is accounted for n c 3 4 1 10 cm W and we first ignore plasma generation (N 0 e = ) and collisions. Figures 3(a) and (b) show the spectra of the THz fields produced in argon by a 50 fs two-color pulse with 1 μm fundamental pump wavelength at increasing propagation distances, when using the UPPE and the MAXFLU codes. Although the WE and UPPE spectra may not perfectly match over the shortest propagation distances, e.g., z=10 μm, an excellent agreement is found at further distances z 50  μm for both intensities. These simulations show that, in the absence of plasma generation, both WE and UPPE solutions match in the whole spectral domain over relatively short distances 10 » μm. This property is independent of the pulse duration, which has been counterchecked by another simulation using longer pump duration, 300 L t = fs (see figure 3(c)). Here, the two UPPE and MAXFLU spectra match again over distances less than the pulse length ∼90 μm. Minor early discrepancies are linked to small differences in the initialization of the numerical codes. The convergence speed between the WE and UPPE spectra driven by a Kerr response alone thus does not depend on the distance propagated over the pulse length c L t . This behavior is rather logical, as the Kerr nonlinearity is just treated as a perturbation in the source term Q and does not impact the frequency range of the linear modes in equations (8) and (9). Next, in figure 4, only plasma generation is taken into account, similarly to figures 2(a) and (b). So, the Kerr response and collisions are set equal to zero. The selected intensity level is I 150 0 = TW cm −2 . When the backward-propagation operator is dropped out, the fundamental linear modes beating at the electron plasma frequency pe w are lost and no plasma opacity is allowed, which results in the development of oscillatory components in the frequency range pe n n < of the UPPE spectrum. This is consistent with the linear mode of equation (9) that admits non-zero frequency components in the range 2 pe pe   w w w . In contrast, the MAXFLU spectrum is dominated by plasma current oscillations, which prevail as long as the propagation distances remain of a few plasma skin depths (here, 3.3 m pe d m = ), as evidenced by figure 4(b). Over longer distances, however, both UPPE and MAXFLU spectra merge in the range pe n n > , as photocurrents become the dominating source in the THz generation process. Out off the laser head, long oscillations over longer times proceed from the Bessel function discussed above. Besides the good agreement between our analytical solutions shown in figures 2(a) and (b) and the numerical solutions of figures 4(c) and (d), we can observe that: • At large times the field oscillations are slower for the UPPE solution than for the WE solution (see insets). The oscillation frequency increases as the optical path is augmented. • Whereas, as expected, the spectral region pe n n < is flat in the forward WE solution from z=3 μm, the UPPE spectrum develops oscillations from a minimum frequency min n that increases with the propagation distance z.
• Convergence is almost reached at z=1 mm. The equality min pe n n = (associated to the plasma response of the input field) is met at z=1 cm and spectra match for all frequencies (not shown). At 150 TW cm −2 intensity, one reports a comparable matching rate between the two spectra and fields, being even sped up by the damping of oscillations at low frequencies 10 < THz and the decrease of the current density in time. Indeed the collision term damps the free electron current in equation (3) and thus both UPPE and MAXFLU solutions are also damped to zero over long times ( 190  fs) beyond the laser head. The green curve shows the backward spectrum collected at z 10 =μm from the vacuum-plasma interface in the MAXFLU simulation. This spectrum occupies the region pe  n n , as expected [5,23], since it is emitted by plasma current oscillations over the plasma skin depth. This backward spectrum remains unchanged over propagation in vacuum.
Similar properties of matching can be refound between the solutions of the two models for pulse configurations favoring either a weaker plasma response (thus a more efficient Kerr effect) at smaller intensities or a stronger plasma response achieved at higher intensities. Figures 5(c) and (d) display the evolution of the same two-color pulse having an input intensity of 50 TW cm −2 . The pulse is undergoing an effective Kerr response combined with plasma generation in argon. The corresponding plasma frequency is very weak, 0.53 pe n = THz, which is related to a long plasma skin depth 90 pe d  μm. Even at rather weak spectral amplitudes, the WE and UPPE spectra and fields approach to each other over distances exceeding far this depth, at least from 1 mm, which confirms the important role of the plasma skin depth in the matching process. The spectral shape follows its analytical counterpart plotted in figure 2(c) for z=1 mm. The numerical UPPE/WE spectra merge from z=5 mm until perfectly overlapping at z=1 cm ( figure 5(d)).
In the opposite range of pulse intensities, I 1 0 = PW cm −2 , the peak plasma density increases and a much shorter plasma skin depth-0.75 pe d » μm for 65 pe n » THz as imposed by the incident pulse-should lead to a quicker merging between the UPPE and WE solutions. Matching is indeed achieved at about z=50 μm, i.e., over a few tens of pe d (see figures 6(a) and (b)). The optical field distortions induced by self-steepening and plasma generation are plotted as inset in figure 6(b). They also show a good agreement between the WE and UPPE models at distances 10 > μm. The spectral distributions and THz field amplitudes (inset of figure 6(a)) still reasonably agree with those computed at z = 10 μm from our analytical solution ( figure 2(d)). Besides, for the same high intensity level, it has been recently shown that photocurrents could be the main player for single-color pulses, provided that the pulse duration be short enough, i.e., few-cycle [13]. For this purpose, figures 6(c)-(e) show two spectra computed from a 1 μm Gaussian pump pulse with 8 fs FWHM duration at 1 PW cm −2 intensity and evolving from z=8 μm ( 1 pe d » μm). Again the UPPE and WE spectra nicely approach to each other in the range pe n n > from z 100 » μm (not shown) and they overlap in the whole spectral range at z 1 mm » . In order to illustrate the influence of the electron dynamics that can vary through laser distortions, figures 7(a) and (b) finally show the electron density corresponding to figures 5(a), (b) and 6(c). Between the  two plotted distances, one reports for the 8 fs pulse with 1 PW cm −2 peak intensity a decrease by half the maximum field value resulting in one decade decrease in the electron density. For comparison, the 50 fs pulse with 150 TW cm −2 peak intensity keeps comparable density levels. The matching distance of the UPPE/ WE solutions for those two configurations then become similar: at z=1 mm, the 8 fs pulse with 1 PW cm −2 peak intensity has almost the same charge level as the 50 fs pulse with 150 TW cm −2 peak intensity, and has thus a comparable plasma skin length ( 3.8 pe d = versus 5.4 μm) at this distance. The two skin depths only differ from each other by a factor 2 and, therefore, the two fields display comparable convergence speed between the two models. This justifies that the number of skin depths needed for matching the two solutions is not universal, as it also depends on the changes in the electron density induced by the distortions of the laser field along propagation. Similar behaviors could be reported from longer (300 fs) pulses (not shown).

Conclusion
In summary, we have derived one-dimensional analytical solutions describing both bidirectional (WE) and unidirectional (UPPE) propagating light pulses by means of a perturbative approach. Structural differences between the WE and UPPE solutions have been explored thanks to those analytical solutions, especially the shape of the THz spectra over mm-range propagation distances. The convergence between both solutions in the z t , ( ) plane has been examined from direct numerical computations integrating Maxwell-fluid equations and the unidirectional pulse propagation model. Even if discrepancies in the linear propagators occur at large times beyond the laser head (ct z  ), the UPPE solution matches its WE counterpart in the z t , ( ) region where the nonlinearity is effective. The extent of the convergence region increases with the propagation distance. Numerical simulations covering a wide range of pulse configurations confirm that, over a propagated distance larger than some plasma skin depths, the UPPE and Maxwell-fluid solutions superimpose to one another. As a result, WE and UPPE spectra match remarkably well over all the spectrum, including the range pe  w w . To conclude, we demonstrated that, in a one-dimensional geometry, the UPPE model, which only governs the forward pulse component, is able to provide similar spectra to a bidirectional Maxwell-fluid model over distances where Kerr nonlinearities as well as photocurrents drive THz pulse generation. Further studies should aim at testing this property in full 3D propagation geometries.
contributions. Since we are interested in THz generation, one should use simultaneously a fine spectral resolution and a fine time step in order to correctly describe the low frequency spectrum below pe n and the two-color laser pulse components including its higher harmonics generated along propagation. The time window of our simulations is, therefore, set to 3.33 ps corresponding to a frequency step of for the UPPE simulations. The highest resolutions used in the MAXFLU code have been employed when it was necessary to decrease the background noise in the lowest parts of the pulse spectrum (e.g., for a Kerr response alone).
Let us finally notice that, so far, we have neglected linear dispersion P E L 0 1  c = *