Generation of High Order Harmonics in Heisenberg-Euler Electrodynamics

High order harmonic generation by extremely intense, interacting, electromagnetic waves in the quantum vacuum is investigated within the framework of the Heisenberg-Euler formalism. Two intersecting plane waves of finite duration are considered in the case of general polarizations. Detailed finite expressions are obtained for the case where only the first Poincar\'e invariant does not vanish. Yields of high harmonics in this case are most effective.

Present and forthcoming high power laser developments [1][2][3] open a vast area in the study of nonlinear physics phenomena related to the electromagnetic field interaction with matter and vacuum [4][5][6][7][8][9][10][11][12]. Vacuum polarization effects leading to high order harmonics generation in vacuum in the course of the collision of the extreme intensity laser pulses is one of the more important directions in these investigations [13][14][15][16]. This problem has attracted substantial attention [17][18][19][20][21][22][23][24][25][26][27][28][29] because it sheds light on the dynamical properties of Quantum Electrodynamics (QED) vacuum in strong electromagnetic fields. The generation of high order harmonics plays an important role in the steepening of a nonlinear electromagnetic wave, in the intersecting of strong laser beams [30] and in the formation of the relativistic electromagnetic solitons in the QED vacuum [31,32].
In spite of a number of publications [19,[22][23][24][25][26][27] devoted to the detailed theoretical analysis of the harmonic generation, the theory of this process is far from complete. There are several theoretical questions that need to be more thoroughly clarified. They concern the role of the various small parameters governing the problem that should be ordered and considered simultaneously, taking into account a specific form of the Heisenberg-Euler Lagrangian [33][34][35][36] describing the QED vacuum polarization effects.
The present paper is devoted to the theoretical consideration and clarification of the high order harmonics generation mechanism within the framework of the Heisenberg-Euler electrodynamics. We consider here the case where the high order harmonics are generated by two crossing plane electromagnetic waves in vacuum under the conditions where only one of the two Poincaré invariants of the electromagnetic field do not vanish. We assume that the finite length electromagnetic pulses are infinite in the two transverse directions. Applying a Lorentz transformation we can find a boosted frame of reference where the crossing electromagnetic pulses ap-pear in the form of two counter-propagating waves (e.g. see [38]). Thus the problem becomes one-dimensional, non-stationary and is described by two independent variables. We consider the case of electromagnetic waves with wavelengths much larger than the electron Compton scattering length and with amplitudes substantially lower than the Schwinger field, which corresponds to the conditions of validity of the Heisenberg-Euler approach [33][34][35][36]. In this context we recall that the problem of extending the Heisenberg-Euler electrodynamics has attracted much interest from the quantum field theory community. One of the ways along which such a generalization can be performed is by taking into account the two-loop correction to the ground state energy of the virtual electronpositron sea in an almost constant external electromagnetic field, see for example [37] and references therein.
Here we prefer to consider only the lowest non-vanishing order in terms of powers of fine-structure constant α, assuming the amplitudes of in-coming waves in terms of the Schwinger field as given. This approach allows us to restrict our investigation to the original Heisenberg-Euler Lagrangian and to simplify considerably our results in the framework of the Heisenberg-Euler electrodynamics approach.
The paper is organized as follows. In Section II we present the nonlinear wave equations within the framework of the Heisenberg-Euler electrodynamics which are used throughout the paper. Section III describes the electromagnetic configuration and the Dirac light-cone coordinates convenient for analysing the problem under consideration. Section IV is devoted to the derivation of a convenient form of nonlinear wave equation. The formulation of the scattering problem and the introduction of a perturbation theory which is used for obtaining final results are in Section V. In Section VI we write the expressions giving the intensity of the high order harmonics. Section VII is devoted to a general case corresponding to different polarizations. The discussion and the summary of the results obtained are presented in Section VIII.

II. NONLINEAR ELECTRODYNAMICS EQUATIONS DESCRIBING THE QUANTUM VACUUM
The analysis of the electromagnetic wave interaction in the QED vacuum is based on the Heisenberg-Euler Lagrangian density, L HE [33][34][35][36]. The sum of the Lagrangians, describes the electromagnetic field in the long-wavelength limit. Here is the Lagrangian of classical electrodynamics with the electromagnetic field tensor F µν defined in terms of the 4-vector potential A µ as [39] Here and in the following we use dimensionless variables with = c = 1, m e and e are the electron mass and elementary electric charge, and α = e 2 /(4π) ≈ 1/137 is the fine structure constant. Thus the electromagnetic fields are measured in units of m 2 e /e , i.e. are normalized on the QED critical field E S = m 2 e c 3 /e = 1.32 × 10 18 V /m. In the Heisenberg-Euler theory, the radiative corrections are described by the L HE term on the right hand side of Eq.(1). It can be written as [36] The invariant fields e and b are expressed in terms of the Poincaré invariants respectively. The dual tensorF µν is defined bỹ where ε µνρσ is the Levi-Civita symbol in four dimensions (ε 0123 = −ε 0123 = −1). In the 3D notation the Poincaré invariants are As explained in Ref. [36], the Heisenberg-Euler Lagrangian in the form given by Eq.(4) should be used for obtaining an asymptotic series over the field invariant e and b assuming that they are small. In this limit, the function M(e, b) in Eq.(4) can be expanded for small arguments as: Here Γ(x) is the Euler Gamma function [40]. The properties of the expansion of M(e, b) in series of powers of the fields e and b are discussed in Ref. [41]. The expression (9) yields for the Lagrangian L HE in the weak field approximation (see also Ref. [43]) with the constant κ = m 4 /8π 2 . In the Lagrangian given by Eq.(10) the first, second and third terms on the right hand side correspond to four-, six-and eight-photon interactions, respectively.

III. THE ELECTROMAGNETIC FIELD CONFIGURATION IN THE DIRAC LIGHT-CONE COORDINATES
As explained in the introduction we consider two counter propagating waves along the x-axis. The vector potential corresponding to the electromagnetic field can be presented in the form Below we use Dirac's light cone coordinates x − and x + defined as (see e.g. Ref. [45]) The Lorentz transform of the light-cone variables under a boost along x with velocity β is given by with a prime denoting variables in the boosted frame of reference and η equal to The following shorthand stands for the partial derivatives: Then the derivatives are transformed as: respectively. For simplicity we consider in this section interacting electromagnetic waves of the same polarization. The general polarization case is considered in Section VII. In the case of the same polarizations, the invariant G vanishes identically. The coordinate system can be chosen so that A y = 0. We will use a notation a(x, t) for A z (x, t).
Calculating the derivatives of a(x, t) with respect to the coordinates x + and x − , u = ∂ − a and w = ∂ + a, (18) we find the relationships between u and w and the electric, e z = −∂ t a, and magnetic, b y = −∂ x a, field components. They read The field variables u and w are not independent of each other. Due to commutation of partials, ∂ −+ a = ∂ +− a, the field variables u and w obey the equation The Lorentz transformation of the fields u and w shows that the w and u are transformed as the − and + components of a contravariant 4-vector, respectively. The field product uw = (b 2 y − e 2 z )/2, is Lorentz invariant and is equal to the first Poincaré invariant F.

IV. DESCRIPTION OF THE WAVE INTERACTION IN TERMS OF CANONICAL MOMENTUM CONSERVATION EQUATIONS
In the case of interacting electromagnetic pulses with the same polarization, the Lagrangian given by Eq.(1) can be written in terms of the field variables u and w, defined by Eq.(18), as is the normalized Lagrangian, where it is assumed that uw 1 and The function Q(ζ) can be represented in the form of the series proportional to the Bernoulli numbers [40], B n . These coefficients can be obtained from the general expansion of M(e, b) in series of powers of e, and b, (see Ref. [41]) setting e = 0 and b = √ 2ζ, where ζ is the argument of the function Q(ζ) [see Eqs. (26) and (27)]. It represents actually the field invariant F in Eq. (8). Several leading order coefficients b m in the expansion (27) are presented in Tab. I. Using the asymptotic dependence of the Bernoulli numbers at m → ∞ (see Ref. [40] (Sec. 24.11)) we obtain the asymptotic expression for the coefficients We see that the coefficients of the series for Q(ζ) in Eq.
Thus the condition m m 0 is very weak allowing the use of a substantially large number of the terms in Eq. (27). This problem was considered from different points of view also in Refs. [41,42]. The Euler-Lagrange equation corresponding to the Lagrangian (24) and the equation (20) gives the field equations. They can be written as a current conservation equation We can define the canonical momentā where Q (ζ) denotes differentiation of Q(ζ) with respect to its single argument. Then Eq. (32) takes the form Using the above obtained relationships we can find the explicit form of the field equation. They can be written as We can use Eq. (37) to exclude either ∂ − w or ∂ + u from Eq. (36). As a result we obtain a set of two equivalent symmetric field equations: This form of the field equations is convenient for applying a perturbation approach because in the region where w = 0 Eq. (38) takes the form whereas in the region where u = 0 Eq. (39) becomes The momenta introduced above can be represented in the same regions as respectively. It is also useful to introduce the energy momentum tensorT where a λ = ∂ λ a and δ µν is the Kronecker delta. The energy momentum tensor is useful to express the conservation of energy and momentum of the field In terms of the variables u and w the tensorT µ λ components can be written as We see that the diagonal components,T + + =T − − , are equal to each other due the tensor symmetry T tx = T xt , see also Ref. [44].

V. PERTURBATION THEORY TO OBTAIN THE SOLUTIONS OF THE SCATTERING PROBLEM
We formulate the scattering problem for the field equations (36)(37) or (38)(39).
In the case of non-interacting waves, the functions u and w do not change since they are equal to u = u(x − ) and w = w(x + ). Let us consider the case where the functions u = u(x − ) and w = w(x + ) describe the two pulses of finite length. The length is L + for the pulse propagating in the x + direction and it equals L − in the case of the pulse propagating in the in the x − direction. We assume that the amplitude of the u(x − ) function is approximately equal to u 0 at |x − | L − , and it becomes exponentially small as compared with u 0 outside this region. We assume also that there are no spatial structures in the u(x − ) pulse with scale-length much less than L − . We set analogous assumptions for the w(x + ) function with the obvious substitution: − → +. As was noticed above, the pulse amplitudes are substantially small, u 0 w 0 1. We will consider the case where the nonlinear effects, leading to harmonics generation of different orders, can be treated perturbatively. Thus, we should exclude too small and too large values of the ratio L − /L + . This point will be considered below. Here we present only a simplified condition for the applicability of the perturbation approach used below: This strong inequality is Lorentz-invariant, that means that all its terms have the same dependence on the boost parameter η given by Eq. (15). This criterion will be made more precise below when we will further specify the parameters of the pulses.
Let us now take into account the interaction between the waves u and w. It takes place only in the following 'rectangular' domain, Ω, on the (x + , x − )-plane, where both the following inequalities are valid: In this 'rectangular' region Ω, both amplitudes u and w depend on both the independent coordinates: x + and x − . Outside this region, when |x − | L − or |x + | L + , the waves do not interact with each other, and w depends only on the coordinate x + , whereas u depends only on the coordinate x − . We may distinguish incoming and outgoing waves in the latter combined region: Outside all these regions, including the Ω region, either the field u or the field w, or both are exponentially small or vanish. The scattering problem can be formulated now as follows: to find the outgoing waves w out (x + ) and u out (x − ), knowing the incoming waves w in (x + ) and u in (x − ).

A. Perturbative solution of 1st order
Equations (38) and (39) can be rewritten in the following integral form that takes into account the initial conditions (52) and (53): In view of Eqs. (54) and (55) we can formally send x ± to ±∞ in Eqs. (56) and (57) and obtain the following scattering relationships: For sufficiently small uw we are able to treat Eqs. (58) and (59) perturbatively, and to substitute into the integrands u in and w in instead of u and w, respectively. Then we obtain: (61) For this perturbation procedure to be be relevant, the second terms in the right hand sides of these equations must be much less than the first terms. If we take for this estimation only leading term in Q(ζ → 0) ∝ ζ 2 , then we obtain the strong inequality (50).
The leading order of Q(ζ) ∝ ζ 2 , caused by 4-photon interaction, leads to the following delays and of the u out and w out pulses and with respect to u in and w in , respectively. This corresponds to the phase shift between the interacting waves discussed in Refs. [44,46]. The strong inequality (50) ensures that the delays are much shorter than L − and L + , respectively. These delays do not correspond to harmonic generation, although they determine the applicability of our approach. The next, subleading orders of Q(ζ → 0) cause generation of the harmonics.
Our approach uses the following set of small parameters: α, u 0 w 0 , αw 2 0 L + /L − , and αu 2 0 L − /L + . Here we address the question whether the next order corrections to Eqs. (60) and (61), as they would appear after the next iterations of Eqs. (56) and (57), or the subleading terms in Q are more important for harmonic generation.
To answer this question we use the structure ofL, given by Eqs. (24), (26) and (27), and Eqs. (56) and (57) as the equations of our perturbation theory starting from zero order: u = u in and w = w in . The combined result of the iterations of Eqs. (56) and (57) and of retaining the higher order terms in the expansion of the function Q(u, w) can be written symbolically as whereÛ m andŴ m are linear integro-differential operators that do not contain any small or large parameters. They act implicitly on all arguments x ± to the right of them separately. Here we use the expansion (27) of Q(ζ).
Due to the conservation of the kinematic momentum, the harmonic with the number n of the wave u will be generated by the terms proportional to u n , u n+2 , etc.
The largest term of this type corresponds to the term in Eq. (64) with ( , m) = (1, n − 1). All other terms contributing to the same n-harmonic will have higher orders of α or/and higher orders of uw. Thus the leading order contribution to the amplitude of the n-harmonic comes from the term of the form w(αuw)u n−1 w n−1 in agreement with Eq. (60). The same conclusion can be drawn for the w-wave.
We may conclude that the higher order terms in the expansion of the function Q(uw) determine the leading order contribution to harmonic generation, and hence the expressions (60) and (61) are sufficient to calculate the intensities of the harmonics in the leading order relatively to the small parameters.

A. General considerations
We consider here two plane waves of the same polarization, when the incident spectrum corresponding to the wave u in is narrow enough. In this case we may choose u in in the form: where U 0 is a smooth function of its argument with maximum valueŪ 0 , and decaying at least exponentially outside the region |x − | L − . The carrier frequency of the wave is ω − . We assume that ω − L − is significantly larger than one. We assume that U 0 is almost constant at the spatial-time scales of the order of 1/ω − . The spectral relative half width of the wave is about Instantaneous intensity of this wave (averaged over its carrier period) is equal to Here the asterisk ' * ' denotes the complex conjugation. This expression is in accordance to the expressions for the energy momentum tensor (46)-(49) in terms of u and w. We will skip analogous comments below. It is obvious that the spectral width of the n-harmonic, ∆ω (n) − grows with n. We calculate ∆ω (n) − below in this section. We are able to separate the harmonics with adjacent numbers, when ∆ω (n) − is less than, for example, half spectral distance between the harmonics, ω − : We will not be interested here in the detailed form of the harmonic spectral line for separated harmonics, but only in the total intensity in the harmonic of order n, integrated over its spectral line.
As we explain above, the leading term of the nharmonic wave field, u out (n) (x − ), is determined by the following expression, obtained from Eq. (58): where and Only the first term in Eq.(72) contributes to u out (n) (x − ), the other terms contribute to the n − 2-, n − 4-and lower harmonics. Hence we obtain: Then the leading term of the 'instantaneous' intensity I out n (x − ) of the n-harmonic averaged over its carrier period can be expressed as: We recall that this intensity is integrated over the spectral line of the n-harmonic. The normalized intensity becomes This expression is invariant relative to boosting. We may rewrite the latter expression by order of magnitude: where a coefficient depending only on n has been dropped. The last two terms in Eq.(76) are small parameters ensuring the validity of our approach. It breaks completely even on a qualitative level, when αw 3 0 u 0 L + ω − becomes approximately equal to unity. In this case, secular effects in evolution of the wave shape become very important. The theory presented in Refs. [30,31] treated the case when the latter parameter is approximately equal to or much higher than unity. Equation (74) is written for a general profile of w in (x + ) that does not contain any additional spatial parameter besides L + . However, when this profile has the form similar to Eq. (66) we find that with ω + L + 1, then (w in ) n+1 becomes considerably less than W n+1 0 for even values of n. It means that in this case harmonics with even numbers will be significantly suppressed with respect to the estimation (76).
The spectrum of the n-harmonic is determined by the multiplier U 0 (x − ) n exp(inω − x − ) in the expression (73). If we approximate U 0 (x − ) in the region close to its maximum as whereŪ 0 is a constant equal to an amplitude of the u inwave, then for sufficiently high n the spectrum of the n-harmonic will be determined by the following x − dependent factor: Its Fourier transform is proportional to Hence the spectral width is In this case, the condition (69) for harmonic spectral separation is fulfilled, when The latter inequality is non-restrictive. It allows us to skip consideration of the opposite case, when n is considerably is above n max .

B. Formulas for the n-harmonic generation in physical units
We use in this section the presentations (66) and (77) for the in-waves, with the profiles given by We introduce a relationship between the intensities of two incoming waves, averaged over carrier periods, and the parameters |U 0 | and |W 0 | as follows: where I u and I w are maximal intensities of the u and w in-waves, respectively. The results described by Eqs. (74) and (75) contain the dimensionless parameter (w in ) n+1 that depends on the w-pulse profile. We present the parameters below in the particular case where W 0 is a real number: We do not consider here phase envelope effects that are important for a short enough w-pulse. We set actually a certain phase to have maximal yields for the even harmonics of the u-wave. The coefficient g n (ω + L + ) in Eq. (87) depends only on the harmonic number and on the parameter ω + L + describing the in-pulse of the wwave. The plots of g n versus ω + L + are shown in Figs. 1 and 2 separately for odd and even n. For even n the coefficient g n becomes equal to zero at ω + L + → ∞, whereas for odd n they tend in this limit to constant values: The ratio of the energy of the n-harmonic, E out n , in the outgoing u-wave to the initial incoming energy of this wave is equal to where f n = n 2 (n + 1) 2 b 2 Tab. II presents the values of f n and g n (∞) for several n.
The dimensionless and Lorentz invariant factor in the brackets in Eq. (89) is one of the small parameters of the theory. It is used instead of the simplified form αw 2 0 L + /L − introduced above to provide a more accurate description in the limit ω − L − 1. The limit ω + L + → 0 in Eq. (89) corresponds to the case when the w-wave becomes actually a constant amplitude cross electromagnetic field. This fact is in accordance with the results of Refs. [30,31], where the generation of the second harmonic due to 6-photon interaction plays a main role in the steepening of nonlinear waves on the constant background of a counter-propagating cross electromagnetic field.

A. Definitions and governing equations
We use here the same frame of reference that was introduced in Section III. Taking into account the general expression (11) we may define the fields where, and below in this section, i, j, . . . span the set {y, z}. The electric and magnetic field components can be expressed as Here e ij is the skew-symmetric 2-matrix with e zy = 1.
The field invariants f and g can be written as The Lagrangian for this problem can be written in a way similar to the case given by Eq. (24). It reads and AlthoughL and Q are now functions of four field components they can be presented in the form of functions of two arguments which are two field invariants, f and g. We note that the LagrangianL and Q, introduced in Section IV, that describe the interaction of two electromagnetic waves of the same polarization are functions of the two field components u and w, which were combined in one field invariant uw. The field invariant g is a pseudo-scalar, hence L should be an even function of it. As a result all coefficients b m, n in the sum given by Eq. (97) with odd n vanish. It is possible to obtain explicit analytic expressions for the coefficients b m, n . However the final expressions are too cumbersome for their presentation in this paper. We will present below their particular values, appearing in the problem for perpendicular polarizations.
In a way similar to that used in Section IV, the normalized Lagrangian (96-97) gives the equations of the Heisenberg-Euler electrodynamics They can be written as As a result, using Eq. (96), we obtain (100) To get a closed system of the field equations in terms of u i and w i (instead of A i ) we may use the following direct consequence of Eq. (91) The following system of the field equations, that is equivalent to Eqs. (100-101), is more suitable for the perturbation theory in terms of α: and Integrating Eqs. (102) and (103) along the lines x ∓ = const respectively, we obtain the following exact relationships between in and out states: and

B. Perturbative approach
Here we consider the case of small |f | and |g|. We substitute in the expressions on the right hand sides of Eqs. (104) and (105) the unperturbed u and w i.e. u in (x − ) and w in (x + ). Assuming strong inequalities similar to those given by Eq. (50), we obtain for the out-wave The subscript in at the bracket means that f and g entering it are calculated with u in i and w in i . Here the u in i waves do not depend on the x + coordinate. Similar expressions can be written for w out i − w in i waves.

C. Harmonic generation by two counter-propagating waves with perpendicular polarizations
To demonstrate how polarization can affect harmonic generation we consider here the case where the incoming waves have only u z and w y components, i.e. where u in y = w in z = 0. In this case Eq. (106) becomes: Note that f = 0 and g = u z w y exactly in this case.
Hence we need to know only the coefficients b 0, n in the expansion (97). They are equal to:  We see that expression (107) is quite similar to what we have in Section V A, see for comparison Eq. (60). The differences are only in the subscripts y and z and in the appearance of this different coefficients in the present expansion of Q gg in comparison with the coefficients in the expansion of Q f f (in our present terminology). These correspondences allow us to use our previous expressions (75) and (89) to obtain new expressions for the intensities of the high order harmonics for the present combinations of polarizations.
In this way we obtain the following results. The normalized intensity of the n-harmonic (integrated over the separate spectral line of the harmonic) becomes (n ≥ 3): Here the subscripts y and z indicate the direction of the electric fields of the u-and w-waves respectively. The subscripts 1 and n indicate the number of the harmonic in incoming and outgoing u z -waves. The parameters ω − and U z, 0 (x − ) describe the incoming u z -wave in accordance with definition (66). We imply the obvious insertion of the subscript z into that expression. The parameters L + and (w in y ) n+1 , describing the incoming w y -wave are defined by Eq. (71) with the obvious insertion of the y subscript. Here, L + is the finite length of the w y -wave in the x + direction.
For the total relative energy of the n-harmonic in the outgoing u z -wave we now have E out z, n E in z, 1 =f n g 2 z, n αω − L + |W y, 0 | 2 2 Ū z, 0Wy, 0 2(n−1) .
(111) wheref n = n 2 (n + 1) 2 b 2 0, n+1 andŪ z, 0 ,W y, 0 , and g z, n (ω + L + ) are defined analogously to Eqs. (83), (84) and (87). Plots for the functions g z, n (ω + L + ) are presented in Fig. 1, whereas Tab. IV presents the values off n for several n. The limits g z, n (∞) are presented in Tab. II. We may recall here that only the odd harmonics is generated for such a combination of polarizations at least in this approximation, whereas the suppression of the even harmonics for parallel polarization was only approximate even in our present order. We see also that the numerical coefficients f n are considerably larger than the numerical coefficientf n with the same n: compare Tabs. II and IV. It means that when the polarizations of the two waves are parallel the high harmonic generation, is much more effective than in the case of perpendicular polarizations.

VIII. DISCUSSION AND CONCLUSIONS
We analyse the high order harmonics generation during the interaction in vacuum of two strong plane electromagnetic waves. The harmonics generation occurs inside the region of the waves intersection. The results obtained are presented in the frame of reference where the waves are counter-propagating.
The theory presented in our paper takes into account the nonlinear polarization of the QED vacuum caused by the existence of a sea of virtual electron-positron pairs. The theory is developed within the framework of the Heisenberg-Euler electrodynamics implying relatively low amplitude and long-wavelength interacting waves.
The regime of high order harmonics generation is described by using the perturbation theory giving the expressions for the high order harmonics to leading order in terms of the small parameters. We show that the leading contributions to the number n harmonic are obtained by expanding the Heisenberg-Euler Lagrangian in power series of the electromagnetic field amplitudes and not by iterations of the nonlinear contribution to the Heisenberg-Euler Lagrangian. In terms of the QED theory we can say that the amplitude of the n-th harmonic is determined by the 2(n + 1)-photon Feynman diagram with one electron loop, but not by the combination of several diagrams with fewer photon lines attached to the electron loops.
Although this statement is based on the consideration of two interacting plane waves we may expect that it is valid in a general case. To see this, one can evaluate the impact of nonlinear currents entering the electromagnetic field equations of the Heisenberg-Euler electrodynamics, using only unperturbed electric E in (x, t) and magnetic B in (x, t) fields. As a result, the harmonic intensities can be calculated as the electromagnetic radiation emitted by these nonlinear currents. The unperturbed E in (x, t) and B in (x, t) fields can be obtained by solving the linear Maxwell equations. The emitted electromagnetic fields, E out (x, t) and B out (x, t), can also be found by solving the linear Maxwell equation with an external source produced by the nonlinear currents. This approach gives the leading order intensity of all harmonics for any polarization of the in-waves (see also Refs. [22,24,26,47]).
From the expressions obtained above it follows that to generate one quantum of the 5-th number harmonic in the collision of two counter-propagating laser beams focused in the 1µm focus spot a laser intensity approximately equal to 5 × 10 26 W/cm 2 is needed. Here we assume a 1µm wavelength laser beam with 30 fs duration. To have the same photon amount for the 3-rd harmonic an intensity approximately equal to 3.3 × 10 25 W/cm 2 is required. This estimate corresponds to interacting waves with parallel polarizations. The interaction of perpendicular polarized pulses is significantly less efficient in generating high order harmonics.