Highly efficient XUV generation via high-order frequency mixing

The efficient generation of the coherent XUV light via frequency conversion of intense laser drivers is a problem of both fundamental and technological importance. Increasing the intensity of the generated high harmonics by raising the intensity of the driving field works only up to a point: at high intensities, rapid ionisation of the medium limits the conversion efficiency. Considering the combined effect of the phase-matching and of the blue shift of the driving field during its propagation in a rapidly ionising medium, we show that the latter can be the dominant limiting mechanism. We introduce a new spatial scale, the blue-shift length, which sets the upper bound for the quadratic intensity growth of the generated harmonics. Moreover, we show that this seemingly fundamental restriction can be overcome by using an additional generating weak mid-IR field. For specific combinations of frequencies of the generating fields, the corresponding high-order frequency-mixing process does not suffer from the blue shift of the drivers and phase mismatch, and thus its efficiency grows quadratically with propagation distance. Our results thus open a new route for highly efficient generation of coherent XUV light.


Introduction
The importance of the macroscopic aspects of high-order harmonic generation (HHG) has been understood from the very beginning of the study of this process. In particular, some key phase-matching properties specific to the intense-field domain were investigated already in early papers [1][2][3][4]. The essential feature of the HHG phase-matching problem, which makes it more complicated than the well-studied case of low-order harmonic generation in crystals [5], is the temporal variation of the refractive index of the medium. Namely, the HHG by an intense laser field is tightly connected with the rapid medium ionisation by this field, and the free electrons appearing due to the ionisation give a substantial-and usually the main-contribution to the dispersion. This leads to a blue shift of the fundamental and to the temporal variation of the HHG phase mismatch.
In this paper we study the phase matching of HHG and HFM analytically and numerically. In the analytical study we consider both the blue shift of the fundamental and the HHG transient phase matching. Having the same origin-the temporal variation of the refractive index-these two phenomena turn out to influence the HHG process in a way that can be characterised with a single spatial scale, which we introduce as the blue-shift length. We demonstrate that the quadratic growth of the macroscopic HHG signal is limited by the blue-shift length under typical HHG experimental conditions. Moreover, we show that this limitation for HHG can be overcome in HFM using an additional weak low-frequency field with a specific frequency. In this process, for some HFM spectral components, the blue shift as well as the phase mismatch vanish.
The generality of our approach brings us to the conclusion that the macroscopic efficiency of HFM is not limited by the fundamental phase-matching aspects that define the HHG efficiency. This makes phase-matched HFM a prospective way to intense coherent XUV generation.
We compare these analytical results with numerical simulations based on the direct integration of the propagation equation for the fields in a 1D geometry. The nonlinear polarisation is calculated numerically by solving the 3D time-dependent Schrödinger equation (TDSE) for a single-active electron atom in the field at each propagation step. This approach allows us to capture all features of 1D propagation (at least within the single-active electron approximation). The calculations are done using IR, visible and near-UV fundamental fields. The numerical results confirm our analytical findings concerning both the HHG and HFM spectral evolution with the propagation, as well as the total XUV energy generated due to these two processes.

Theory: interplay between phase matching and blue shift
Already in the early studies of HHG [1] the following length scales characterising macroscopic aspects of the process have been introduced: (i) the coherence length where the phase mismatch for the generation of the qth harmonic is Δk = k q − qk 0 , k q and k 0 are the harmonic and fundamental wavevectors, respectively; and (ii) the absorption length where N at is the gas density and σ is the ionisation cross section for the harmonic field. We proceed with an accurate description of HHG for conditions where absorption is negligible (or the absorption length L abs is long), which is a valid assumption for most experimental conditions and sufficiently high harmonics, and where the coherence length L coh can significantly vary within the pulse due to the ionisation of the medium during the generation.

High-order harmonic generation
We start from the wave equation for the field E(x, t) linearly polarised along the z-axis and propagating along the x-axis through the medium with the nonlinear polarisation P(x, t). The wave equation (3) can be simplified [5] to the reduced first-order wave equation by applying the slowly-varying amplitude approximation under the conditions of modest medium density and laser intensity, which are discussed in detail in [29,31]. If we define the Fourier transform for some function G(x, t) as and slowly-varying amplitudesG(x, ω) of the spectral components G(x, ω) as the reduced propagation equation can be written in the form In this work we solve the reduced propagation equation (6) including the influence of the transient ionisation of the medium. Neglecting the modification of the fundamental pulse envelope F 0 during the propagation in the medium, we can write the fundamental field as where ζ = t − x/c. Here we take into account the temporal variation of the refractive index of the fundamental and assume it to be linear in time: wheren = n(0) is the refractive index at the centre of the pulse (t = 0) and a = ∂n(0)/∂t. According to this, the coherence length (1) corresponding to the centre of the pulse is defined as L coh = πc/(qω 0 |Δn|), where Δn = 1 −n (for XUV the real part of the refractive index differs negligibly from unity). For the refractive index (8), we can present the nonlinear polarisation due to the qth-order harmonic of the field (7) in the form whereF(t) = F(t) exp(iqω 0 t) and F(t) is the microscopic response. Note that from equation (9) one can see that The Fourier transform (4) of the polarisation (9) is where F (ω ) is the spectrum of the microscopic response F(t). Substituting the polarisation (10) into equation (6), taking into account equation (5), and integrating it over x from 0 to L with a change of variable to ω = ω − qω 0 ax/c, we find This qth-order response is non-zero near the harmonic frequency, i.e. when ω ≈ qω 0 . However, the harmonic line is broadened due to the combined effect of the phase matching and the blue shift of the fundamental.
Assuming the Gaussian lineshape and integrating equation (11), we derive the field amplitudeẼ q (L, ω) in the form where (14) defines the harmonic lineshape, γ = (ω − qω 0 )/Δω is detuning from the exact harmonic frequency, and κ = Δωτ is a constant defined by the pulse shape. For the Gaussian pulse (12) where the constant h characterises the chirp of the harmonic pulse. Here we use h = 0 and the case of h = 0 is discussed in detail in appendix A. The harmonic lineshape can vary, but the dependence of the harmonic energy on the propagation distance is independent of the chirp and can be considered as for the non-chirped pulse. The harmonic pulse chirp originates in the dependence of the microscopic response phase on the fundamental intensity. The coefficient in this dependence is inversely proportional to the cube of the fundamental frequency [32], therefore the chirp of the harmonic generated by visible or UV fundamental is much less than the one for the IR fundamental.
The key length scale for the field amplitude (13) propagation is the 'blue-shift length', where n i and n f are the fundamental refractive index at the beginning and end of the generation, respectively. We take into account that for the Gaussian lineshape |n f − n i | = aτ , where τ is the pulse duration, or more precisely, it is the time interval within which the ionisation takes place. Results for other pulse lineshapes are presented in appendix B.
We study the full energy of the qth harmonic as a function of the propagation distance, as shown in figure 1(b). There are three distinct regimes of the harmonic energy's behaviour with propagation distance, corresponding to the different behaviours of the refractive index illustrated in figure 1(a). The first situation [red lines in figures 1(a) and (b)] is when the detuning from the phase matching is much more than its variation during the generation, that corresponds to L bs L coh . This case is well studied because it takes place for low-order harmonic generation in crystals [5]. For HHG this condition is realised when the geometrical or atomic dispersion dominates; however, note that it is not very typical for HHG, except for the cases where a tightly-focused laser beam is used, leading to high geometrical dispersion. In this case, the harmonic intensity initially grows quadratically, saturates at L = L coh , and then decreases.
The second case [green lines in figures 1(a) and (b)] is usual for HHG by an IR field when, already under low ionisation degree (at the beginning of the generation), the dispersion of the free electrons compensates the one of the neutral atoms. For further generation the detuning from the phase matching is defined mainly by the free electrons, and Δn = |n f − n i |/2 or L bs /L coh = 1/2.
Finally, for HHG using high-pressure gases [33], or visible/UV fundamental [34], the free-electron dispersion can compensate the atomic one near the centre of the pulse [blue lines in figures 1(a) and (b)]. Thus, for L bs < L coh , the HHG signal grows quadratically up to L bs . For longer propagation distances the signal continues to grow, but it does so more slowly due to the transient compensation of the phase mismatch.
In order to define the character of this growth, we consider the harmonic spectrum (13) shown in figure 1(c). One can see that initially the spectrum is centred near ω = qω 0 . During further propagation, the blue shift of the fundamental δω 0 = ω 0 aL/c (18) results in the blue shift qδω 0 of the harmonic polarisation response; therefore, the new XUV components are generated at the blue side of the line, and the red side does not change. After some propagation, this shift exceeds the harmonic response linewidth, and this XUV does not add coherently to the XUV generated initially, i.e. the total signal does not increase quadratically.
Using the definitions of a and κ, we find that for L L bs the blue shift of the harmonic response is much larger than the harmonic response linewidth: qδω 0 /Δω = πL/(κL bs ) 1. For these conditions most of the XUV energy is confined within the spectral region [see also figure 1(c)]. Within the spectral region (19) the integration limits in equation (11) can be changed to ±∞, and then the integral can be taken analytically for arbitrary lineshape F (ω ): Outside of the spectral region (19) (or, more precisely, for ω < qω 0 − Δω and for ω > qω 0 + qδω 0 + Δω) the integrand in equation (11) vanishes, so the intensity (20) vanishes as well.
Taking into account that within the spectral range (19) the XUV intensity (20) is approximately constant [given by equation (20)] and vanishes outside of this spectral range, we find that the energy (17) emitted due to qth-order process is W q ≈ 4π 2 a 2 |F(−Δn/a)| 2 qδω 0 . Hence, for the condition L L bs assuming that contributions (20) of different q do not spectrally overlap and using equation (18), we can finally rewrite the total energy as Here the presence of the term F(−Δn/a) has a clear physical interpretation: this is the amplitude of the response at the time instant t = −Δn/a, when the phase matching is achieved. The term 1/a appears because the faster the ionisation is, the shorter the temporal window of the phase matching is, and thus the lower the intensity of the generated field becomes. One can conclude thus that for L > L bs the blue shift of the fundamental plays a crucial role in the harmonic line formation. Its peak intensity does not grow during further propagation and its width grows linearly with the propagation distance, leading to the linear growth of the total harmonic energy.

High-order frequency mixing
As a possible solution to overcome this limitation in the energy growth, we suggest the XUV generation by an intense fundamental field together with an added, much weaker, low-frequency field. This results in an HFM process, in which q photons from the fundamental (ω 0 ) field and m photons from the low-frequency (ω 1 ) one are converted into single photon with the frequency Note that q and m are integers, and |q| + |m| has to be odd. In [9,11,17,35] it was shown (in [9,11,35] for the case of a weak second field, and in [17] for arbitrary intensities and processes' orders) that the phase-matching problem for HFM can be considered exactly in the same way as it was done earlier for low-order processes: the detuning from the phase matching is where k l (l = 0, 1, 2) is the wavevector at the frequency ω l , indices 0, 1, 2 correspond to fundamental, low-frequency and HFM fields, respectively. Below we consider co-directed k 0 and k 1 . If (i) plasma and/or capillary contribution to the dispersion dominates and (ii) the frequencies of the generating fields ω 0,1 and of the generated one ω 2 sufficiently exceed the plasma frequency ω pl ω l , then the refractive index for a frequency ω l is n l = 1 − ω 2 pl /(2ω 2 l ) and the detuning from the phase matching for the process can be written as As far as the last term in (23) is negligible for high q, one can see that for given numbers q and negative m there is a specific frequency for which the detuning (23) vanishes irrespective of the plasma density. Thus the HFM can be phase-matched during the whole laser pulse even when the ionisation degree changes significantly during the pulse. This phase matching was experimentally demonstrated in [8] for ω 0 = 2ω 1 , q = 6 and m = −3.
Let us consider the blue-shift contribution to the HFM process. Denote the ionisation-induced blue shift of the fundamental by δω 0 . Since the blue shift δω 1 of the low-frequency field is inversely proportional to its frequency (24), it is therefore given by δω 1 = δω 0 q/|m|. Hence, the polarisation response at the frequency ω 2 under condition (24) is not blue-shifted, because the frequency shifts of the two generating waves compensate each other: Let us also comment briefly on the blue shifts for the two field components. The blue shift of the fundamental δω 0 and of the low-frequency field δω 1 are defined by the temporal variation of the refractive indexes at these frequencies: δω 0,1 ∝ ω 0,1 ∂n 0,1 /∂t, and ∂n 0,1 /∂t ∝ (1/ω 2 0,1 )∂(ω 2 pl )/∂t. Here we again assume that the plasma dispersion dominates. Thus, the blue shift of a field is inversely proportional to its frequency, so for the low-frequency field given by equation (24) the blue-shift is δω 1 = δω 0 q/|m|.
Taking into account also the neutral gas dispersion contribution to Δk q,m , one would find non-zero phase mismatch. Denoting this contribution as Δk at q,m , we find the coherence length (1) for this process as L coh(q,m) = π/|Δk at q,m |. However, for the conditions considered here the coherence length is typically much longer than the absorption length (2). Therefore, HFM efficiency in our conditions is limited by the XUV absorption.
Thus, from this analytical consideration, we conclude that the HFM macroscopic signal can keep the quadratic growth up to much longer propagation distances than the HHG signal, because HFM can be phase-matched for varying ionisation degrees and the HFM line is not broadened by the plasma-induced blue shift of the generating fields.

Numerical method
To verify our analytical study, we perform numerical calculations, based on the solution of the reduced propagation equation (6) for the external field along with the 3D TDSE for the quantum-mechanical atomic response on this field, in a self-consistent fashion as suggested in [5] and applied in [36,37].
The spectrum of the nonlinear polarisation response P(x, t) included into the reduced wave equation (6) isP(x, ω) = −N at d z (x, ω), where N at is the initial atomic density and the dipole moment d z in a spectral representation is calculated via the electron acceleration or the total force f z (x, ω) acting on the electron as d z (x, ω) = −f z (x, ω)/ω 2 . Therefore, the polarisation response can be rewritten as The choice to consider the polarisation within the total force f z gives the opportunity to include naturally the correct contribution of the free electrons due to the atomic potential and also to take into account the depletion of the electronic wavepacket rigorously. We solve numerically the TDSE for the model atom in the single-active electron approximation in the external field E(x, t). The numerical method is described in [38]. A soft Coulomb potential [39] providing the ionisation energy equal to the one of Ar is used. The size of the numerical box for TDSE integration for the fundamental wavelengths 600 nm and 800 nm is chosen so that only the short electronic trajectories contribute to the microscopic response, as was done in [40]; when suppressing the long trajectory contribution we take into account the fact that experimentally this contribution produces a very divergent XUV beam and thus does not contribute much to the observed XUV signal on the beam axis. For the fundamental wavelength of 330 nm the contribution of all trajectories is taken into account. The initial field is a sum of the fundamental field with amplitude E 0 and frequency ω 0 , and the much weaker low-frequency field with amplitude E 1 and frequency ω 1 . Sine-squared envelopesã 0,1 (t) are used for both components of the field. Then we calculate the expectation value of the time-dependent total force f z (x, t) as and apply the Fourier transform to it [f z (x, ω)]; from there we obtain the nonlinear polarisation response (25). We then substitute the polarisation response into the reduced wave equation (6), which allows us to find the total field E(Δx, t) after passing through a thin layer of matter of thickness Δx. This field is used again to get the nonlinear polarisation response by repeating this procedure in an iterative way for every new thin layer of medium.

The XUV spectra generated due to HHG and HFM
We consider a near-UV fundamental with central wavelength 330 nm, and a weak low-frequency field with central wavelength 1980 nm [ω 1 is given by equation (24) with m = −1 and q = 6, hence ω 1 = 1 6 ω 0 ]. The peak intensity of the fundamental is 2 × 10 14 W cm −2 and the low-frequency field is 100 times less intense for all our calculations, unless mentioned otherwise. The pulse duration (defined by intensity FWHM) is 100 fs. We use a medium density 3 × 10 18 cm −3 throughout. As it propagates through the medium, the fundamental pulse ionises it and this leads to the blue shift of incident fields [see figure 2(a)]. Note that the low-frequency component is too weak to ionise the medium. The free-electronic density is below the critical density even for the low-frequency field, so the latter can propagate in the ionised gas. However, its blue shift is much more pronounced than for the fundamental, as seen in figure 2(a).
The spectra of the XUV generated due to the HFM with frequency 6ω 0 − 1 6 ω 0 and the spectrum of the 7th harmonic are shown in figure 2(b) [shown also as a colourmap in figure 3(a)]. Here we see that the linewidth of the harmonic grows with the propagation distance, and the peak intensity does not change much after approximately 0.26-0.52 mm, in agreement with L bs = 0.2 mm, found from (16) for these conditions. Moreover, the modification of the high harmonic lines with the propagation length agrees very well with the theoretical prediction presented in figure 1(c). In contrast, the line generated due to HFM remains narrow (despite the blue shift of both generating fields), and its width does not change much as a function of the propagation distance, as a result of the blue-shift compensation. Figure 4 shows the generated XUV spectrum [shown also as a colourmap in figure 3(b)] in a frequency range near the 17th harmonic for a 600 nm fundamental field and a low-frequency field with frequency ω 1 = 2 17 ω 0 (λ 1 = 5.1 μm), using a pulse duration of 90 fs. The lower fundamental frequency allows us to consider processes of higher order, and a larger number of them. Different peaks of the spectrum correspond to HFM processes with different q and m and also to HHG processes for m = 0. The difference in the intensity of the spectral components originates from the difference in the microscopic response as well as from the difference in the phase matching of the generation.
It is also important to notice that in figures 3(b) and 4 the HFM enhancement is restricted to only a few peaks. This is caused by the relatively low value of q in use here, and if the targeted harmonic order q is chosen at higher values, the enhancement bandwidth increases significantly, see also [41].

Efficiency as a function of the propagation length
The propagation behaviour of two spectral components, the HHG signal and the HFM signal, under three different conditions, is shown in figure 5. One set of results is obtained for the same parameters as in figure 2 above, for 330 nm fundamental wavelength and 1980 nm low-frequency field wavelength. Two more sets are for the same fundamental as in figure 4 (600 nm wavelength) and for low-frequency fields of two different frequencies: ω 1 = 2 17 ω 0 (λ 1 = 5.1 μm) and ω 1 = 2 13 ω 0 (λ 1 = 3.9 μm). One can see that high-harmonic energy as a function of the propagation length tends to saturate after the propagation distance near L bs , shown in figure 5 with arrows. A higher fundamental frequency and lower harmonic order (blue lines) lead to longer L bs . Moreover, the spectral components with the combined frequencies ω 2 keep growing almost quadratically for the whole calculated propagation distance, demonstrating perfect phase matching for the HFM process. The efficiency of the HFM for even longer propagation distances should be limited by the XUV absorption, as well as by nonlinear modification of the low-frequency generating pulse.

Comparison of analytical and numerical results
Here we discuss the effect of absorption (which is included in our numerical study) in comparison to the analytical theory (which neglects it). Figure 6 shows in linear scale the comparison of the harmonic energy calculated numerically and analytically for two different values of the L bs /L coh ratio: for graph (a) we use an IR pump, where the electronic dispersion dominates, while for graph (b) we use near-UV pump, where the electronic and atomic dispersion are comparable and compensate each other at a certain time instant. We see that for graph (a) the harmonic energy growth becomes slower as the propagation distance exceeds the blue-shift length, and that such sudden change is absent for conditions of graph (b), in both cases in agreement with the theoretical prediction. Note that for the propagation distances exceeding 0.3 mm the XUV absorption reduces the harmonic energy (the XUV wavelength is about 47 nm in both cases), therefore the calculated energy is less than the one predicted by our theory, which omits the absorption.   . XUV spectra generated by the field with the initial fundamental wavelength 600 nm, the initial frequency of the low-frequency field ω 1 = 2 17 ω 0 at the two propagation distances. Numbers of quanta from the fundamental field q and from the low-frequency one m are shown for the several peaks generated via frequency mixing.

Effect of the generating pulse intensity and duration
We now turn to the effect of the fundamental peak intensity and pulse durations on the HHG and HFM processes. Figures 7(a) and (b) presents the energy of the spectral peaks generated due to HHG and HFM using 330 nm fundamental and 1980 nm weak field (ω 1 = 1 6 ω 0 ) for the peak intensities 2 × 10 14 W cm −2 [graph (a)] and 3 × 10 14 W cm −2 [graph (b)]. Namely, for conditions of figure 7(a) we find L bs = 0.32 mm for the short pulse and L bs = 0.20 mm for the long one. For conditions of figure 7(b) due to high fundamental intensity used in this calculation the ionisation saturates even for the short pulse, so both for the short and long pulses we have L bs = 0.16 mm. Again, these distances agree with the behaviour of the harmonic energy shown in the figure.
In figures 7(a) and (b) one can notice that HFM spectral component grows almost quadratically for the propagation distance where HHG saturates. The efficiency of the HFM for even longer propagation distances should be limited by the XUV absorption. Note, however, that deep ionisation of the medium, Figure 5. Energies of harmonics (dashed) and HFM components (solid) as functions of the propagation distance. The initial fundamental field wavelength is 600 nm, the initial frequency of the low-frequency field is ω 1 = 2 17 ω 0 (green lines) and ω 1 = 2 13 ω 0 (red lines). For blue lines the initial field is the same as in figure 2(a). The blue-shift lengths are shown with arrows of corresponding colours.  Energies of harmonics (dotted) and spectral components with the combined frequencies satisfying equation (24) (solid) as functions of the propagation distance. The initial central fundamental wavelength is 330 nm, the initial intensity is (a) 2 × 10 14 W cm −2 and (b) 3 × 10 14 W cm −2 , the long pulse duration is 100 fs and short pulse is 50 fs long, whether for 'combined pulses' the fundamental pulse is 50 fs and the low-frequency one is 100 fs long, the initial frequency of the low-frequency field is ω 1 = 1 6 ω 0 . (c) The same as (b) but for ω 1 = 2 13 ω 0 .
from approximately 50% to 99% depending on the pulse duration and intensity, lowers the absorption due to the reduced number of neutral atoms. In figure 7(c) we show the results obtained in the similar conditions as considered before [see figure 7(b)], but for nonlinear optical processes of higher order. One can see that for long propagation distances the HFM signal both for the long and short fundamental pulses stops growing and demonstrates a complicated oscillating behaviour. The origin of this is the strong modification, both in amplitude and in phase, of the low-frequency pulse due to ionisation of the medium by the fundamental. This modification includes not only a blue shift, but also pulse stretching, walk-off and other effects. The natural solution to reduce the influence of such effects is to use a longer low-frequency pulse. Figure 7(c) shows that the application of combined pulses indeed leads to an improvement of the HFM process.

Conclusions and perspectives
In this paper we study macroscopic aspects of HHG and HFM analytically and numerically. The latter study is based on the integration of the 1D propagation equation with the nonlinear polarisation found by numerical 3D TDSE solution; we calculate the XUV spectra generated at different propagation distances using fundamental wavelegth in IR (800 nm), visible (600 nm) and near UV (330 nm) range, different intensities and pulse durations. We find that the blue shift of the fundamental is one of the main factors limiting the quadratic growth of the macroscopic HHG signal as a function of the propagation distance. As HHG is always accompanied by the ionisation of the medium, the blue shift of the fundamental is inevitably present during the propagation. It results in the shift of the central frequency of the harmonic polarisation response, which equals to the fundamental's blue shift multiplied by the harmonic order. After some propagation this shift exceeds the linewidth of the harmonic response. Further propagation leads to XUV generation in the spectral range above this linewidth, so this XUV does not add coherently to the XUV generated initially, and the total signal does not increase quadratically any more.
We introduce a new spatial scale characterising the macroscopic HHG process, namely, the blue-shift length L bs = πc qω 0 |n f −n i | . After propagation over a length L bs , the frequency shift of the harmonic polarisation response is equal to the initial harmonic response linewidth. We show that the quadratic growth of the HHG signal is limited by the shortest of the following distances: the blue-shift length, the coherence length and the absorption length. Moreover, for weakly absorbed XUV the limitation by the blue-shift length is very typical. For longer propagation distances the HHG signal grows linearly, even when the atomic dispersion is transiently compensated by the electronic one.
We also demonstrate that for HFM, at certain frequency ratios between the two generating waves, both the phase mismatch and the blue shift vanish for arbitrary ionisation degree. This leads to the generation of narrow spectral lines whose intensity keeps growing quadratically above distances where the HHG signal saturates. The HFM signal is limited by the XUV absorption, as well as by nonlinear modification of the low-frequency generating pulse. We show that the influence of the latter factor can be reduced using longer low-frequency pulses.
In this way, HFM using an intense laser field and a weaker low-frequency field is a prospective way of intense coherent XUV generation. Its experimental realisation becomes very feasible with the current progress of mid-IR laser sources, as well as secondary far-IR and THz sources: indeed, the related phenomenon of high-order parametric generation has already been observed recently [41]. Moreover, since HFM components of different orders can be phase-locked, our work also paves the way to the use of HFM for the generation of attosecond pulses.
The phase of the harmonic microscopic response ϕ q depends on the fundamental intensity I, and for every quantum path contribution this dependence is assumed to be close to linear, as ϕ q = −αI [42,43]. This behaviour leads to the approximately linear harmonic chirp, which is described by term h in equation (15). Figure 8 presents the comparison of the spectra calculated for non-chirped pulse [shown in figure 1(c)] and for the chirped one. The latter value corresponds to maximal fundamental intensity of 2 × 10 14 W cm −2 and α = 1.2 × 10 −14 cm 2 W −1 . Such α value is typical for the short quantum path for 800 nm fundamental.
The harmonic line modification due to the chirp of the harmonic response can be understood as following. At the rising edge of the laser pulse, the generation is better phase-matched because L bs /L coh = 0.4 corresponds to the transient phase matching, achieved in the beginning of the pulse. The harmonic blue shift is added with the fundamental blue shift, providing the spectral maximum in the blue side of the harmonic line, see figure 8. For the generation at the falling edge, the red shift of the harmonic response leads to emission at the low-frequency part of the harmonic line. This emission is weak because of the weak phase matching at the falling edge of the pulse. Thus, even within the approximation of the constant fundamental blue shift our calculations reproduce the 'foot-shape' of the harmonic line, observed in numerous experiments. As we mentioned above, the α coefficient decreases with the fundamental  Finally, the dependence ϕ q (I) does not influence at all the harmonic energy as a function of the propagation distance presented in figure 1(b).

Appendix B. Pulse shape of the harmonic response
The microscopic harmonic response F(t) appearing in equation (9) is temporally confined within time interval when the fundamental intensity is high enough so that the cut-off energy I p + 3U p exceeds the photon energy for the considered harmonic. Therefore, this temporal interval is usually shorter than the ionisation duration τ , at least for the harmonics in the higher part of the plateau. For κ (15) the duration of the harmonic response (full duration at the half level of the intensity) is τ/2. However, even in this case the wings of the Gaussian pulse F(t) are outside of the time interval where the assumption of the linear variation of the refractive index (8) is applicable, and this might be a potential source of serious mistakes in our analytical consideration. To check this we consider another dependence Using equation (B2), we calculate numerically the XUV spectrum (11) and the full harmonic energy (17). We find that results are very similar to ones found using Gaussian spectrum (12), except the case L bs /L coh = 0.5. Figure 9 presents the comparison of the full energies calculated using a Gaussian pulse [same as figure 1(b)] and the temporally-confined pulse given by equation (B1). One can see that for L bs /L coh = 0.4 for both pulses the energy goes on growing for L > L bs although with slightly different slopes. For L bs /L coh = 0.5 the energy increases for the Gaussian pulse case but saturates for the temporally-confined pulse case. This can be understood from equation (21): the slope is defined by the factor F(−Δn/a), and taking into account that F(−Δn/a) = F(−τ L bs /L coh ), one can see that for L bs /L coh = 0.5 this factor is zero in case of the temporally-confined pulse (B1) but non-zero for the Gaussian one. However, being defined by the wing of the Gaussian, this slope is mild. For L bs /L coh > 0.5 the slope is even milder and for L bs /L coh > 0.77 the signal oscillates in both cases. In practice, the difference between the temporally-confined and the Gaussian pulses is not crucial, and thus a Gaussian pulse is a reliable approximation in our approach.