Theoretical aspects of intense field harmonic generation

The authors present theoretical studies of high-order harmonic generation in a rare-gas medium. The experimental results obtained at Saclay with a 1064 nm Nd-YAG laser in the 1013 W cm-2 intensity range are summarized. The harmonic emission strengths, first decrease rather steeply for the first orders, then form a long plateau up to the 21st harmonic in xenon, or up to the 33rd harmonic in argon, before decreasing again rather abruptly. The theoretical description of these experiments consists first in the calculation of the photoemission spectra emitted by a single atom. The spectra are obtained by numerically integrating a time dependent Schrodinger equation for the laser-excited rare-gas atom. Second, one must account for collective effects in the medium, described by Maxwell's equations. A theoretical framework for describing the generation and propagation of harmonics in strong laser fields is developed. An numerical solution of the propagation equations for the harmonic fields in xenon at 1064 nm provides results which agree well with experimental data.


Introduction
' An intense laser focused into a 10 Torr rare-gas vapour leads to the generation of very high-order odd harmonics of the pump field. Experiments performed at the University nf!!!inok e! Chicagn (Xtphersnn e! n! !Os?, Xosman e! a! 1988) shnwed !he generz!ion of up to the 17th harmonic of a 248 nm KrF laser in a neon vapour. The 33rd harmonic in argon and the 21st harmonic in xenon have been observed at Saclay using a 1064 nm Nd-YAG laser (Ferray er al 1988, Li er al 1989. This was recently extended to the 25th harmonic of a KrF laser by Sarukura er al (1991) and to !he 53rd harmonic of a 1 ps 1053 nm Nd-glass laser by L'Huillier ef a/ (1991a), both produced in a neon vapour. These experiments performed with short pulses and a! high laser intensities result in very short wavelength coherent radiation, 10 nm (125 eV) for the 25th harmonic of a 248 nm laser, 20 nm (62 eV) for the 53rd harmonic of a 1053 nm laser. Efficiencies are of the order of 10-l'-lO-q. These frequency conversion processes might lead to useful sources of short pulse, short wavelength coherent radiation. They seem to be much more promising than previous experimental investigations performed with longer p!se ir?frared !Ism !Me!rhkov e! a! 1977, Groseva e! a! Wi!denauer !?U! or ultraviolet pump fields (Reintjes et nl 1978(Reintjes et nl , 1981, Bokor er Q / 1983) which were limited by various effects such as the ionization of the medium, the absorption of the generated radiation in optically thick media and/or the breakdown of phase-matching conditions. The harmonic intensities in the strong-field regime exhibit a characteristic distribution. After the expected rapid decrease for the first orders, there is a long plateau, which ends up with a rather sharp cutoff. This behaviour is surprising for the following reasons. If generation of, for example, the 25th harmonic becomes as probable as fifth-order harmonic generation, as is the case for a Nd-YAG laser at an intensity of 3 X lOI3 W cm-' in a 10 Torr argon vapour (Li et a/ 1989), it means that the weak-field picture (see e.g. Gontier and Trahin 1982, Gao and Starace 1989, Potvliege and Shakeshaft 1989a. Pan et a/ 1989, which has been successfully applied for describing (non-resonant) non-linear optical phenomena before, does not apply anymore. One has to go beyond lowest-order perturbation theory for providing a correct description of these non-linear optical phenomena where the field is no longer a weak perturbation of the atomic medium. Time-dependent approaches involving the numerical solution of the Schrodinger equation (Kulander and Shore 1989, DeVries 1990, LaGattuta 1990 and Floquet calculations (Potvliege and Shakeshaft 1989b) have been successfully employed. Classical methods (Bandarage e t a / 1990, Chu et a/ 1990), one-dimensional approximations (Eberly et a/ 1989a, b, c) and many other model calculations (Shore and Knight 1987, Biedenharn et a/ 1989, Becker et a/ 1990, Sundaram and Milonni 1990) also provide some insight into the physics involved. The single-atom photoemission spectra from all these calculations generally reproduce qualitatively the experimental data. This is another interesting and surprising aspect of the problem. Phase-matching conditions are known to play an important role in the overall response of the medium and to be very sensitive to various parameters such as the atomic density, the focusing geometry, the laser frequency and the process order. On the basis of weak-field calculations of harmonic generation, one would expect phase matching to severely degrade with the order Kulander 1989, L'Huillier et a/ 1990). In contrast, the comparison between experimental data and single-atom calculations seems to indicate that propagation effects either play no role or affect all the harmonics in the same way.
Thus, the interpretation of the experimental results is not an easy task because it involves both the single-atom response to a strong laser interaction and the many-atom response, the capability of the medium to ensure proper phase matching between the (non-perturbative) non-linear polarization induced by the incident field and the propagating harmonic radiation. The purpose of the present paper is to discuss this double aspect (microscopic and macroscopic) of harmonic generation processes. The conversion efficiency depends both on the quantities which govern the interaction of a laser with a single atom, the atomic system, the laser wavelength and intensity, and also on the parameters which affect the propagation of the fields such as the atomic density or the interaction geometry. Our emphasis will be on this latter part (propagation), which has received little attention up to now, whereas in contrast, a lot has been done concerning harmonic generation by a single atom (see the references mentioned above).
Studies of harmonic generation in gases up to 1987 have been extremely well documented in reviews and textbooks (Hanna er a/ 1979, Reintjes 1984, Shen 1984, Arkhipkin and Popov 1987, Delone and Krainov 1988. The experimental and theoretical results obtained more recently have been described by L'Huillier et a / (1991b) and we refer the reader to this work for a comprehensive review. In the present paper, we chose to focus our attention on fewer results which we analyse as thoroughly as possible.
We present the theoretical method used for describing these harmonic conversion processes. It consists first in the calculation of the single-atom photoemission spectrum performed by numerically integrating a time-dependent Schrodinger equation Shore 1989, 1990). Then, we solve the paraxial propagation equation using as a source the non-linear polarization induced by the radiating atomic dipoles.
Our method is applied to the interpretation of one particular result, harmonicgeneration in a 15 Torr xenon vapour irradiated by a 1064 nm Nd-YAG laser at intensities of the order of 10'' Wcm-2 (L'Huillier et a/ 1991~). The paper is organized as follows. In section 2, we present in parallel experimental results and the single-atom calculations. Although the emphasis of this paper is essentially theoretical, we thought it useful to include a summary of the experimental method and of the main results. We describe in particular the results obtained in xenon (Li et a / 1989, LomprC et a / 1990) which provide the basis for our theoretical analysis. We present the corresponding single-atom emission spectra. This is the first part of our theoretical analysis, which will allow us to separate the contribution of the single-atom response and that of phase matching, an essential step in the understanding of the experimental results. We also give some details about the method used for integrating the time-dependent Schrodinger equation Shore 1989,1990). In section 3, we develop a formalism for describing the generation and propagation of harmonics in a medium exposed to a strong laser field. We show how this can be applied to the analysis of the experiments performed in xenon at 1064 nm. Our final results agree well with the experimental observations, reproducing in particular the plateau behaviour. We try to unravel the role of phase matching in strong-field harmonic generation.
2. Experimental results and single atom photoemission spectra 2.1. Experimental siudies 2.1.1. Method. A harmonic generation experiment consists of focusing an intense laser radiation into a rather dense rare-gas medium (a few Torr) and then analysing along the propagation axis the vuv light emitted during the interaction. A schematic picture of the experiment principle is shown in figure 1. Most of the experiments reported here have been done by using a mode-locked Nd-YAG laser (40 ps pulse width, 1064 nm wavelength), with a maximum energy of 1 GW at a IO Hz repetition rate. The linearly polarized laser is focused into the interaction chamber. A useful quantity which characterizes the focus is the confocal parameter b = 2?iw,3A, wo denoting the beam radius and A the wavelength. For a Gaussian beam, b is equal to twice the distance on the propagation axis over which the beam section increases by a factor of two.
A pulsed gas jet (Kung 1983, Bokor el a / 1983) produces a well collimated atomic beam with a 1 5 Torr pressure at 0.5 mm below the nozzle of the jet (Lompri er a1 1988). The atomic density distribution in the jet has been measured to be very close to a Lorentzian distribution with a 1 mm full width at half maximum ( L ) . An advantage of gas jets compared with differential pumping systems is that low background pressures can be kept in the interaction chamber and in the detection chamber (see figure 1) with reasonable pumping equipment. Moreover, the interaction region is small (typically 1 mm long) and it can be exposed to an intense laser pulse in an almost collimated beam geometry. This can be very advantageous for phase matching in a positively dispersive medium (L'Huillier et a1 1988) and for avoiding significant absorption of the emitted radiation. The vuv detection system, which covers a broad spectral range from 10 to 350 nm, is described in figure 1. A grazing incidence holographic grating separates the different components of the light emitted along the propagation axis. The light is then detected by photomultipliers [ A > 120nm) or electron multipliers ( A < 120 nm). This monochromator has a good detection efficiency owing to the lack of an entrance slit ( s e e figure 1). The number of photons produced at each frequency can be estimated to within one order of magnitude by using the spectral efficiency of the photon converter, the grating efficiency and the absolute electron multiplier gain.

Results.
Typical spectra consist of series of harmonic peaks superposed on a broad background. Figure 2 shows, for example, the spectrum between 70 nm and 110 nm obtained in Xe at 1064nm, 3 x 10'' W cm-* with a 75 mm focal lens. Although the monochromator did not have a sufficient resolution for separating all the lines, most of the features could be attributed to discrete transitions in Xe, Xe', Xe2+. This fluorescence background is probably due to dielectronic recombination processes, following the formation of a plasma by multiphoton ionization. These processes occur on a much longer time scale than harmonic generation and at an intensity high enough for ionizing-at least partially-the non-linear medium. This light emission is incoherent and probably isotropic. Whether there is an underlying continuous background is an open question: these processes should be studied with better resolution and at another detection angle (e.g. perpendicularly to the laser axis).
In figure 3, the number of photons produced in a 15 Torr xenon vapour is plotted as a function of the harmonic order for several laser intensities between 5 x IO" W cm-' and 3 x 1013 W cm-2 (Lompr6 et a/ 1990). This result has been obtained with a 200 mm focal length (the confocal parameter b is estimated to be 4 mm). Only odd harmonics are observed, which is to be expected for harmonic generation in an isotropic gaseous monic signal decreases with the order. As the intensity increases, a plateau followed by a rather abrupt cutofi appears. Its length increases with the laser power, up to the intensity at which the medium becomes ionized with a probability close to unity, above 1.3 x 10" W cm-'. Above this intensity, the signal increases much less rapidly, the distribution becomes smoother but the maximum observable harmonic order remains constant (equal to 21 for this pulse length). The vertical scale gives an order of magnitude estimate of the number of photons produced at each laser shot. This means a power efficiency of 10-8-10-9 for the plateau harmonics at the highest laser intensity. The brightness is estimated to be 10'' photons/s 8, mrad'.
Another way of looking at these results is to plot the number of photons as a function of the laser intensity. Figure 4 shows, for example, the behaviour of the 15th harmonic. All the intensity dependences of the harmonics present a common feature: the number of photons increases first rapidly, then saturates when the medium gets ionized. There are two reasons why the ionization of the gas limits harmonic generation (Miyazaki andKashiwagi 1978, Reintjes 1984, L'Huillier el a/ 1990). The main medium responsible for harmonic generation (the neutral atoms) gets depleted when the medium becomes ionized. Harmonics are still produced in the periphery of the interaction volume or at the beginning of the laser pulse. Ions could also generate harmonics, but their response is expected to be less efficient at these intensities. The second effec! that might limit harmonic conversion efficiencies is the breaking of phase matching owing to the presence of free electrons in the medium. These free electrons have a nonnegligible effect on the refractive index (at the fundamental frequency and at the medium, with inversion symmetry. . A! the !owest in!cnsity ( 5 U 10'2 w cm?), the har-Harmonic order harmonic frequencies). They induce a large positive phase mismatch between the generated beam and its driving polarization, which can reduce the conversion efficiency (see section 3.3).
Beiow saturation, the harmonics vary as some power law of the laser intensity (Lompr6 et a / 1990). The lowest order harmonics ( q s 9) vary as predicted by lowest order perturbation theory. The 11th and 13th harmonics exhibit a more complex intensity dependence, which might be due to the influence of discrete resonances and which is much less rapid than the I" or perturbative power laws. The highest harmonics vary all in the same way, approximately as the 12th power of the laser intensity. iihis is consistent with the piateau observed in the harmonic intensity distributions (figure 3). Indeed, the fact that there is a plateau means that all the harmonics must have approximately the same power law. Because the intensity dependences deviate from a I4 power law, where I is the laser intensity and q the process order, the atomic response cannot be described within the weak-field limit. The measured power law reflects, however, both the single-atom contribution and the coiieciive response of medium (phase maichingj which may also be power depcndeni. Similar results have been obtained with the other heavy rare gases (krypton and argon). The distributions obtained in Xe, Kr, Ar at 3 x l O I 3 W cm-* are shown in figure 5 (Li er n/ 1989). The conversion efficiency decreases from Xe to AI, which is not surprising, since xenon is more polarizable than lighter rare gases. However, the maximum order that can be observed increases from 21 in Xe, 29 in Kr to 33 in Ar. inc J~I U riarriiunic ( J L nm, JO e v , 1s L~L -bnuriesi wavcicrrgrri I~U I~U U I I p i u u u c~u W L L I I the Nd-YAG laser system (limited, however, to about 20 mJ in 40 ps). Atoms with higher ionization energies have in general a lower conversion efficiency but they can produce more harmonics. Moreover, they can experience a higher laser intensity without being ionized. Shortening the laser pulse length also increases the intensity at which an atom ionizes and therefore might lead to the production of higher order harmonics. Recently, the 53rd harmonic of a 1 ps, 1 p m Nd-glass laser was reported in neon at an intensity of 5 x 10" W cm-2 (L'Huillier et a / 1991a).
Other experimental studies involve the variation of the parameters that influence the macroscopic aspect of the interaction, the atomic density and the focusing conditions. Let us briefly summarize the main conclusions. First, the harmonic signal is found to vary as the square power of the atomic pressure (from 1 to 25 Torr) independent of the gas, the harmonic order and the laser intensity (Li et a / 1989). This is the signature of a coherent process. Indeed, a n incoherent sum of single-atom dipole radiation would give a linear dependence of the signal on the atomic density. The square dependence shows that the measured light results from a coherent sum of the radiating dipoles. (In the first case, one sums the intensities; in the second case, one sums the amplitudes and then squares.) The second study is of the dependence of the signal on the interaction geometry (i.e. on the volume within which the harmonics are created). This helps in understanding phase matching, which is expected to be strongly affected by focusing (Bjorklund 1975, L'Huillier el a/ 1990). As will he shown in section 3, the measured number of photons is proportional to b'lFq(b)12, where b is the confocal parameter, proportional to the focal section S and where the quantity denoted IF,(b)/*, the square of the phase matching function, reflects the propagation of the fields throughout the medium. The b' factor has the following origin: b' arises from the coherence of the process, the fact that the harmonic signal vanes as the square of the number of atoms involved. Since the medium is limited in the propagation direction by the length of the gas jet, this must be understood as the number of atoms in the transverse direction, in the focal plane, so that the number of photons then varies as S2 or b'. The detected signal is a number of photons, proportional to the focal section of the harmonic field, which is itself proportional to the focal section of the incident field S. Hence the additional factor of b. Experiments have been performed with different geometries (Lompri er al 1990) varying from confocal focusing ( b = L ) to a plane-wave situation ( b >> L). In figure 6, we compare the harmonic intensity distributions obtained with b = 1, 4 and 6 mm and L= 1 mm at a 1.3 x IO" W cm-* laser intensity. The signal has been divided by the global factor b'. The difference between the three results, which then reflects the variation of the phase matching factor lF,(b)l2 with b, lies within the experimental error bar. With b = 1 mm, harmonics higher than the 15th could not be detected. That is simply due to a harmonic/background ratio which was barely above one in this case. This result shows that phase matching (lF,(b)l2) does not depend much on the focusing geometry. This contradicts predictions derived from the weak-field limit (L'Huillier et a1 1990; see section 3) and emphasizes the need for a general, nonperturbative, description of harmonic generation processes. Another important practical conclusion is that using a loosely focused geometry considerably enhances the conversion efficiency (as 6').

Harmonic generation by a single atom
The electrons in an atom oscillate in response to a strong laser field. This acceleration of the charge density causes the emission of radiation at odd multiples of the driving frequency. If the field is not too strong, the rate of emission can be calculated using standard perturbative methods. Following the pioneering work of Manakov et a1 (l975), Manakov and Ovsyannikov (1980) and Gontier and Trahin (1982), very accurate high-order non-linear susceptibilities have been obtained by Potvliege and Shakeshaft (1989a). Pan ef nl (1989,1990) and Gao and Starace (1989) for hydrogenic systems.
The calculations of perturbative emission rates for multielectron systems are substantially more difficult, with only the lower order susceptibilities of some of the rare gases having been reported (see e.g. Sitz and Yaris 1968, Manakov et a/ 1975, Manakov and Ovsyannikov 1980, Bishop et al 1988. The strong-field emission from atomic systems has been successfully investigated using two different approaches. First, the expansion of the electronic wavefunction in a Floquet basis has been used by Potvliege and Shakeshaft (1989b) to obtain nonperturbative harmonic emission rates by hydrogen for driving frequencies of 1064 and 532 nm. This method is valid for intensities up to the point where the ionization rate becomes comparable with the laser frequency. For intensities beyond this regime, the assumption of a constant intensity pump is invalid and the atomic response becomes very sensitive to the pulse shape.
The second approach for modelling photoemission from laser excited atoms is the direct solution of the time-dependent Schrodinger equation. This has been accomplished for a number of one-dimensional model systems by Eberiy et al (1989a, b, c) and by Sacks and Szoke (1991); for several model systems by Becker er a Snndaram andMilonni (1990); for a classical electron in the field of a proton by Bandarage et al (1990), Chu er al (1990) and for realistic, three-dimensional atoms by Shore (1989, 1990), DeVries (1990) and LaGattuta (1990). and then discretizing along the radial axis with 5 = (j-0.5)Ar. L is adjusted as required to achieve convergence. Using a threepoint formula for the second derivative, and defining g{ = r,+j, we can rewrite (2.1) as where Ho, the atomic Hamiltonian, couples radial values j to j , jzt 1 and is diagonal in I, while HI, the interaction term, couples angular momenta I to I + 1 and is diagonal on j . In treating many-electrons atoms, we have used effective potentials (Kulander and Rescigno 1991)  where I is the unit matrix and T = A t / 2 . The interaction term is evaluated at the midpoint of the time step. Since the matrices in the brackets are all tridiagonal, the multiplication and inversion can be achieved with vector operations. The propagator is accurate to second order in the time step and is approximately unitary. The computational effort in solving (2.4) is linear in the number of grid points. We can propagate more than 1 . 2~ lo6 spacetime points per second on a Cray Y/MP machine. The prompt emission by an atom in an intense field is due to the oscillation of the electron charge density in the vicinity of the nucleus. Electrons at large radial distances cannot emit high energy photons because momentum cannot be conserved. Therefore, we determine the time-dependent wavefunction in a limited volume near the nucleus by removing the flux which reaches the edges of our grid with a mask function which forces the amplitude smoothly to zero at the boundary. This mask function is applied after each integration step.
The calculation proceeds as follows. We choose a pulse shape,f(t) in (2.11, which rises as the square of a sine function over five optical cycles, then is unity for the next 15-30 cycles. During the first part of the constant intensity interval the transient excitations which occurred during the ramp decay by ionization. Then during the latter part of this interval, we determine the photoemission spectrum and rates. These we obtain by Fourier transforming the time-dependent induced dipole, at odd multiples of the driving frequency above a broad background. The width of the harmonic peaks is determined by the lesser of the two intrinsic time scales in the problem, the ionization time or the pulse length. The background is quite sensitive to the integration parameters (Krause er al 1991) and has no constant phase relationship to the driving field. Therefore, no phase matched, coherent background signal can be expected. The harmonics do have a well defined but intensity-dependent phase delay which weakly affects the phase matching of these fields. We have performed calculations for a fine grid of intensities between 5~1 0 '~ and 5 x IO" W c t f 2 . In figure 8, the harmonic emission strengths for several of these intensities are displayed. This figure shows the emergence of a plateau which increases in strength and extent with increasing pump intensity. The harmonics initially rise very quickly with intensity as predicted by perturbation theory, then rise substantially more slowly with numerous oscillations due to intermediate resonance effects. The intensity dependences of the 3rd, 9th, 17th and 23rd harmonics are shown in figure 9. Only the third harmonic behaves approximately as I 3 over this intensity range. The approximate power laws for the higher harmonics are much lower than their orders, being approximately I 6 for those in the plateau. They are also lower than the experimental power laws.

harmonic order
Similar studies of other rare-gas atoms qualitatively reproduce the observed trends in the harmonic conversion efficiencies. The highest conversion efficiency, for a given pump intensity, was found in xenon which is the most polarizable of the atoms studied. The broadest plateau in the spectrum was found in the system with the highest ionization potential. In figure 10 we show a comparison of the single-atom spectra for Ar, Kr and Xe at 3 x 10" W c K 2 . Allowing for the difference in ionization rates for these three atoms, the shapes of the harmonic spectra in the figure agree very well with the measurements shown in figure 5 [ Li et a1 1989). In order to compare the yields directly, however, the effects of phase matching during the propagation of the harmonic fields through the medium must he determined.

Harmonic generation by an assembly of atoms exposed to an intense field
We first give a general theoretical framework for the description of harmonic generation in a gaseous medium. Then, we analyse the phase matching conditions for the experiments reported in section 2.1 within the weak-field limit, using traditional non-linear optics arguments. In section 3.3, we give the results of a more complete calculation for the propagation of the harmonic fields, which uses as a source the response of an atom exposed to a strong field (see section 2.2). These two calculations lead to opposite conclusions. We find that phase matching of the high harmonics is considerably enhanced in a strong-field regime as compared with the weak-field limit. As will be explained in section 3.4, this is mainly an amplirude effect, related to the harmonics power law, which is much lower in a strong-field situation than in the perturbative case.

Theoretical framework
We start from the general wave equation describing the propagation of an electromagnetic field %' ( r, 1 ) in an isotropic dielectric medium characterized by a n electronic polarization 9 ( r , t ) .
wnen rhe incident fieid is iineariy poiarized, this is a scaiar equation, in the direction of the laser field. 8(r, 1 ) and P(r, 1 ) can be expanded as The polarization of the medium is the distribution of dipole moments induced by the total electric field (and not simply by the laser field). It is the sum of the linear polarization P:=Nx'(-qw; q w ) gq where x ' ( -q w ; q w ) is the field-free dipole po!irizibi!ity ind the rmr?-!inezr po!irizitioa " P: ', which itse!f contzins z number of terms involving the fundamental and the harmonic fields. Among those terms, we can distinguish the non-linear polarization induced by the fundamental field, which is the driving term for harmonic generation, denoted by 9:. There are higher order terms to the linear polarization at frequency qw, which lead to intensity-dependent corrections to the refractive index (Zych and Young 1978, Mahon and  For example, the fifth harmonic can be created from the interaction of the third harmonic field and the fundamental (absorption of one photon at frequency 3w, plus absorption of two photons at frequency 0). The ninth harmonic can be created from third harmonic generation of the third harmonic field, or from the seventh harmonic, plus absorption of two laser photons, etc. The number of possible mixing processes increases very rapidly with the order. In the high-order harmonic generation experiments in gaseous media discussed in the present paper, however, the harmonic conversion efficiency remains weak, so that these indirect processes can be neglected. The study of the pressure dependence of the harmonic signal also indicates that such processes which would exhibit a much higher pressure dependence than N 2 cannot be very important. In the same way, we neglect depletion of the pump field, which has to be taken into account for high conversion efficiencies (Tomov and Richardson 1976, h e l l el a/ 1976, Kilda1 and Brueck 1980). These approximations allow us to decouple the equations describing propagation of the harmonics. Introducing the intensity-dependent dipole polarizability x l ( q w , /gal2), for the fundamental, and for the harmonic fields.
Equation (3.4) describes the propagation of an intense laser beam in a medium, including non-linear effects such as self-focusing (Grishkowsky 1970, Akhmanov et nl 1972, Marburger 1975, Shen 1984. We neglect these effects, which are weak for long incident wavelengths and at relatively low atomic densities, as well as the (defocusing) effects which may arise from the presence of free electrons. %' ] is then simply the laser field propagating into a medium characterized by the refractive index n, (wavevector k, =w/c(l+21rNx'(-w, a))= n,w/c).
Equation (3.5) describes harmonic generation in a medium exposed to a strong laser field, for weak harmonic conversion efficiencies. However, its solution with non-perturhative polarizations and refractive indices remains a formidable problem, and we shall make an additional approximation which will allow us to use a simple procedure for solving numerically the propagation equation. We shall neglect higher. order corrections to the refractive index at frequency qo as we d o for the fundamental, (in the resonance region) and above. However, as will he shown below, the effect of dispersion (i.e. the fact that the refractive index varies with the wavelength, whether it is intensity dependent or not) remains small in the experiments performed at Saclay, at a relatively low pressure. Equation The homogeneous part of this differential equation is linear and we can therefore use an integral representation (Kleinman 1962 , Bloembergen 1965, Lago et a / 1987). The outgoing Green function associated with the homogeneous part of (3.6) is e'*qR/4.rrR, with R=lr-r'l. (3.6) can be written as (Jackson 1975) (From now on, we drop the d index in Pd,.) We assume that the field is emitted close to the propagation axis (paraxial approximation). We are interested in the harmonic field outside the medium, in the far Jield. Let r = ( x , y , z) be a point in the medium and r ' = ( x ' , y ' , z') the observation point, far from the sources. The far-field approximation implies that ( z ' -z ( >>lx'-xl, ( y ' -y l so that R = 2'-z + [(x'-x )~+ (y'-y)2]/2(z'-z ) . Equation  where Ak, = k, -qk, is the phase mismatch between the polarization field and the harmonic field (Akq << k,, qk,). Absorption and (z-dependent) atomic density distrihutions (Rettner er al 1984, Lago er a / 1987) can he accounted for simply by replacing Ak,z in the exponential in (3.9) by the expression t m (3.10) where ~~( 2 ) = Im(Ak,(z)) is the absorption coefficient at frequency qw, equal to X(z)uq/2, where U,, is the photoionization cross section and N ( z ) the atomic density distribution.
The solution of the propagation equation thus reduces to a straightforward threedimensional integral over the non-linear medium. The last problem that needs to he understood in order to calculate propagation of the harmonics, particularly in a strong-field situation, is how to relate the Fourier components d, (=d(qw) of (2.6)) of the dipole moment d ( t ) calculated for a real field g( f ) = ' 8 cos wf, to the propagating non-linear polarization Pq = P,(r) eiqkhz resulting from the interaction of a focused beam with the medium. Consider, for example, an incident Gaussian beam (this is what we shall use in our calculation, though this argument is much more general). The polarization PJr, z) can therefore be obtained from the time-dependent dipole momentcalculatedforafieldstrengthIE,(r, z)l andevaluatedatatimer'= f+q(r, z ) / w . The expression for the envelope PJr, z ) is then Pq(r, z)=2#(2)dq(r,z)exp t a n -' ( 2~/ b ) -~ 2k r 2 z ) ] (3.14) b2+4z2 where dJr, z) denotes the qth harmonic component of the time-dependent dipole moment evaluated for a field strength lE,(r, z ) l . In the weak-field limit, dq(r, z ) is related to the qth-order non-linear susceptibility ,yq (Reintjes 1984) by (3.15) In this case, the integral in (3.9) can be performed almost completely analytically. In a more general situation, where dq(r, z) cannot be expressed simply as a function of the incident field, (3.9) needs to be performed numerically. Note that by making use of the revolution symmetry of the problem, it actually reduces to a two-dimensional integral where Jo denotes the zero-order Bessel function. To conclude this section, we give the expression for the total number of photons emitted, Nq, which is the quantity measured experimentally. Assuming that the laser pulse width is long compared with the light period, N,, is obtained by integrating the harmonic intensity profile both spatially and temporally N,, =-4h9w r'lE@(r', z', t')I2 dr'dt'. (3.17)

3.2.
Phase marching in the weak-field limit In the perturbative limit, the integral in (3.9), with the polarization field Pq(r, z) given by (3.14) can be performed analytically. The 9th harmonic field E, is equal to (Ward and New 1969, Miles and Harris 1973, Bjorklund 1975)

Eo(r', z', 1 ' ) = -i7ik,bXqN~E,4Z'-'GcP(r', z')%?:(t')F,,(b, Ak), (3.18)
Eo is the peak laser intensity. GZm denotes a Gaussian beam envelope for a field oscillating at frequency qw. with wavevector k,, and with a confocal parameter equal to b (3.11). The generated harmonic field E, is Gaussian with the same confocal parameter and beam waist location as the fundamental field. Its focal section (equal to .irb/2kq = bh/49) decreases with increasing order. 9: is the 9th power of the incident pulse distribution. The pulse width of the 9th harmonic field is equal to T/& for a shaped (Gaussian) pulse, T denoting the laser pulse width. Phase matching is described by the factor F,(b,Ak), defined by a one-dimensional integral over the non-linear medium as where T,, 94,(1) dt. Note that the single-atom response at the peak intensity ld,,12=lxq(Eo/2)412 and the effect of propagation in the medium, described by the dimensionless factor IF,(b, Ak)l*, can be factored out and are therefore independent, We have studied the behaviour of F,(b,Ak) as a function of the process order q, the focusing geometry (characterized by the confocal parameter b ) and also the phase mismatch Ak. The gas density distribution used in all the calculations presented in this paper is described by a truncated Lorentzian in the z-direction, p ( z ) = 1/(1+4z2/L2) for I z l s L and p ( z ) = O for ( z ( > L , with a width at half maximum of L= 1 mm. Figure I l ( a ) shows the variation of IF,,(b,Ak)I2 with the phase mismatch Ak (assumed to he real) for the seventh harmonic and confocal parameters b = 4 mm, b = l m m , b=0.4mm. Figure l l ( b ) shows lF,,(b,Ak)('for b = l m m and 9 = 3 , q = 7 , q = 13. Although varying b or 9 leads of course to different results, the general trend is similar: as the process increases in order or as the geometry becomes more focused, the maximum of the phase-matching function shifts towards a higher negative phase 2s !ong 2s the we.k-fie!d approxima!ion rem2ins va!id.  Phase matching is optimized when the phase in (3.19) remains small. The optimum (negative) phase mismatch, which is close to 2( 1 -q)/b, increases with increasing order and decreasing confocal parameter (see figure 11). Indeed, the phase variation across the focus of the non-linear polarization induced by focusing becomes more rapid, leading to a larger phase lag between the generated field and its driving polarization. It must then be compensated by a greater phase mismatch Ak.
The amplitude factor in (3.19) affects the shape of 1FJ2, particularly in the wings region. In the limit of an infinite medium and small coherence length, the harmonic fields generated before the focus are exactly cancelled by those created after the focus yielding no net harmonic generation. For a finite medium, exact cancellations occur only for some periodic values of the phase mismatch, which leads to oscillations in the wings of the phase matching function. The average level of these oscillations depends on how rapidly the amplitude term in (3.19) gets damped away from the focus. Therefore, it decreases as the geometry becomes more focused or, equivalently, as the process order increases. One might say that for the same confocal parameter of the incident beam, the interaction geometry becomes more and more of the tight focus type as the order of the harmonic increases. The phase mismatch Ak is proportional to the difference between the dynamic polarizabilities at the harmonic frequency and at the fundamental frequency: Ak = 2rrqwX[,y'(-qw; q w ) -, y ' ( -w ; w ) ] / c . We show in table 1 the values of the phase mismatch at 15 Torr (5.3 x IO1' atoms/cm') in Xe for the harmonics emitted at multiples of the YAG frequency (LHuillier et a / 1990). By comparing these values with the horizontal scale in figure 11, one sees that perfect phase matching is not realized in the high-order harmonic generation experiments. In a first approximation, one can simply take the intersection of the functions plotted in figure l l ( b ) with the line Ak =O. IFq(b,Ak)12=lFq(b,0)12. Since this is near the positive wing of IF,(b, Ak)l', phase matching in the experimental conditions is essentially determined by the variation of the amplitude of the polarization field throughout the medium, and not SO much by the variation of its phase. (As will be shown below, this conclusion will also be true in a strong-field situation.) In the weak-field limit the polarization field varies as the qth power of the incident field so that the volume in which the harmonics are generated becomes smaller as the process order increases. Phase matching is predicted to decrease with increasing order, leading to the expectation that any plateau in the single-atom response would be destroyed by propagation in a focused geometry.
As a further illustration, we show in figure 12 IF,(b, Ak)I2 as a function of the process order q for the phase mismatches Ak indicated in table 1 corresponding to the experimental conditions in xenon. Two geometries have been investigated: b = 1 mm (circles) and b = 4 mm (squares). The full curves indicate the weak-field limit calculations, the broken curves with symbols are strong-field calculations, which will be discussed later. The broken curves without symbols indicate the results obtained within the weak-field limit by neglecting dispersion (Ak = 0). The comparison between (F,(b, Ak)I2 and IF,(b, 0)l' for the two geometries shows that dispersion does not play an important role here, owing to the rather low pressure used in the experiments. In contrast, the effect of focusing leads to a substantial decrease of lF,(b, Ak)I2 with q, reflecting the fact that the polarization amplitude is more rapidly damped away from the focus. The importance of this effect increases as the laser is more strongly focused, as shown by comparing IFJb = 1, Ak)12 and IF,(b =4, Ak)I2. Note in particular the twelve orders of magnitude difference between the 3rd and 25th harmonics in the tight focusing case. This weak-field analysis of the phase-matching conditions in the high-order harmonic generation experiments leads to the following conclusions: (i) phase matching decreases with increasing order, thus destroying any plateau observed in single-atom emission spectra; (ii) IF,(b, Ak)l' varies with b and this variation depends strongly on the process order q. Consequently, the harmonic signal (b3lF,(b)l2) should not have a simple b3 scaling. These conclusions obviously contradict the experimental ObseNations.

Phase matching in a strong-jeld regime
We now present our calculations of harmonic generation in a strong laser field. As explained in section 2.2, the single atom response, d,, is obtained from the wavefunction generated by numerically integrating a time-dependent Schrodinger equation Shore 1989, 1990). We have concentrated our effort on the xenon atom, which has been experimentally investigated in great detail. In order to describe the non-linear polarization throughout the medium, i.e. for a distribution of intensities, we have calculated the qth component of the time-dependent dipole moment, d,, over a fine intensity grid, between 0.5 and 5 x 10'' W cm-' (assuming lowest-order perturbation theory to be valid below 5 x IO'' W cm-').
The macroscopic parameters of the interaction are chosen to mimic the experimental conditions described in section 2.1. The incident laser beam is assumed to be Gaussian and we consider two cases with confocal parameters b = 4 mm or b = 1 mm (equation 3.11)). The laser pulse shape is taken to be Gaussian with a 36ps width at half maximum. The phase mismatch (Ak) is assumed to he intensity independent and is given in table 1. This is of course not valid for a 40 ps pulse, but it gives an upper estimate of the free electrons, influence. As already pointed out in section 3.1, we also neglect the defocusing effect that these free electrons might have on the propagating fields, considering only the variation of the phase mismatch.
The number of photons emitted at a given harmonic frequency is obtained by calculating the generated field Eq(r', 1') (with 1 ' -1 = (2'-z ) / c ) and then by integrating temporally and spatially the intensity profile lEq(r', t')I2. The results of our calculations in xenon for several laser intensities ranging from 5 x w cm-' to 3 x 10" w cm-2 and for a b = 4 mm laser confocal parameter are shown in figure 13. They can be compared with the experimental results in absolute value (figure 3 ) and also with single atom emission spectra (figure 8). At the highest laser intensity 3 x 10" W cm-2, where ionization becomes significant, we show results obtained with and without including the effect of free electrons (plotted respectively with stars and open circles). Although the additional positive phase mismatch induced by these free electrons is large, especially for the high harmonics, the resulting decrease in efficiency remains weak, at most a factor of two, and sometimes non-existent. The theoretical curves are generally higher than the experimental results particularly for the highest intensity. They are, however, mostly within the experimental error bar, estimated to be one order of magnitude. The general behaviour of the harmonic spectrum as a function of the laser intensity is well reproduced. The effect of phase matching is seen to be most important at the lowest intensity, i.e. in the weak-field limit. High harmonics are not as well phase matched as low-order ones. At the lowest intensity, the cutoff occurs at a much lower order (9 instead of 15) in the many-atom response than in the single-atom response (compare figures 8 and 13). However, at laser intensities above 9 x 10" W c K 2 , the length of the plateau is quite similar to that of the single-atom response (see figure 8). as if all the harmonics were equally phase matched. From the lowest intensity to about 10" W crK2 (below the onset of saturation), the number of photons rises more rapidly with the laser intensity than Id,/*: the effects of imperfect phase matching are lessening with increasing laser intensity.
The dotted curve in figure 4, which was presented in section 2, shows the result of our calculation for the 15th harmonic. There is a disagreement of a factor of two between the calculated and experimental saturation intensities and of a factor of seven for the harmonic signal. However, the relative agreement between both results is remarkable (see the full curve). Our calculation can reproduce fairly well the behaviour of the saturation, which is mainly due to depletion, and also the high power law ( -1 " ) helow the onset of saturation. This power law is higher than the one obtained in the single atom response, again indicating a substantial improvement of phase matching conditions as the intensity increases and as one departs from the weak-field limit.
Similar results as those presented in figure 13 have been obtained with a tighter focusing geometry ( b = 1 mm). The calculations show the same b' scaling that was observed in the experiments (see section 2.4). The shape of the harmonic distribution at high intensity is relatively insensitive to the geometry, in contrast to the weak-field limit predictions discussed previously.
Some information can also be gained by studying the spatial and temporal profiles of the generated fields, which have not been measured experimentally so far. Figure  14(a) presents the spatial profile in the far field of the 15th harmonic in a strong field situation at 3 x 10'' W cm-2. Here, for simplicity, we exclude the effect of free electrons and we take a 'snapshot' of the profile at the maximum of the pulse at time f = 0. The results are shown for two focusing geometries, collimated ( b = 4 m m , broken curve) and confocal ( b = 1 mm, chain curve). All the profiles have been normalized to unity at the maximum. The horizontal scale has been chosen so that the perturhative Gaussian profiles corresponding to the two geometries, which are shown by the full curve in the figure, be superposed. At high laser intensity, for b = 4 mm, the far-field harmonic profile becomes narrower, which means that the harmonic beam is defocused compared with the perturbative limit. Rings appear, particularly in the tight focus geometry. However, the integrated signal does not depend strongly on how the energy is spatially distributed. In the perturbative limit, the high harmonics become more focused with increasing order (focal section varying in I j q ) . By contrast, in the strong-field regime, their focal section remains approximately constant for the harmonics of the plateau region.
We have also studied the temporal profiles of the harmonics, obtained by integrating the spatial distributions [ j lEq(r', t')12r' dr']. These profiles, which are quite insensitive to the geometry, are generally larger than in the weak-field limit. They are of the order of 15 ps for the plateau harmonics for a 36 ps incident laser pulse. Figure 14(b) shows the temporal intensity profile of the 15th harmonic at b = 1 mm and 3 x lo" W cm-'. The asymmetry of the profile arises from the depletion of the population of neutral atoms during the pulse, which reduces the conversion efficiency for positive times (at the end of the pulse). Finally, we found it useful for a better understanding of these results to separate the role of propagation from the single atom contribution in the calculations, as we did in the perturbative case. We define IF,(b, Ak)I2 by analogy with equation (3.19) as (3.22) where Id$ denotes the dipole moment evaluated at the maximum intensity of the pulse (at best focus). Note that, in this definition, IF,(b,Ak)l' not only characterizes the coherence length of the process but also reflects the spatial (and temporal) profile of the emitted harmonics in comparison with the weak-field limit. Here we assume a r Iarb. units1 square laser pulse which is short enough for ionization to be negligible (so that the separation of the single-atom response and propagation remains meaningful). The broken curves at the top in figure 12 (see section 3.2) show the phase matching factor lF,(b, Ak)I2 at 3 x lo" W cm-'. Although in the weak-field limit [FJb, Ak)lL decreases rapidly with increasing order, and is strongly dependent on the geometry, in contrast, at 3 x IO1' W cm-', IF,(b, Ak)I2 remains between 1 and independent of the order. Moreover, it is relatively insensitive to the geometry. Our calculations yield results that are indeed in good agreement with the experimental conclusions and in strong disagreement with the weak-field predictions.

Znterpretation
We now try to explain why going to a strong-field regime for the laser-atom interaction significantly improves the phase matching of the high harmonics. Compared with a perturbative picture, the dipole moment dq obtained from the time-dependent calculation (or from any non-perturbative method) varies generally less rapidly with the laser intensity and similarly for most of the harmonics (which gives rise to the plateau behaviour). It shows structures and resonances (see figure 9). Finally, it has an intensity-dependent phase lag relative to the fundamental. What influences phase matching is the first property, namely the fact that the dependence of the dipole moment on the electric field is (on average) much weaker than that predicted within lowest-order perturbation theory. In order to illustrate this idea and to get some insight into the reason why the high harmonics are phase matched in the same way, we shall now consider a simple model: we assume that the dipole moment d,, varies as the pth power of the incident laser field, p denoting an effective order of non-linearity, lower than the harmonic order (d,(r, z ) = ( q ] E , ( r , z)l"). As in the perturbative limit, the integral in (3.9) can he performed partly analytically, though the expressions are more cumbersome. The harmonic field can be written as (L'Huillier et al 1991h) We have introduced the notations b'= pblq b"= qb/p 4k IF,(b, Pk)lz =$ 1!Fq(r', z')I2r' dr'. (3.25) In the loose focusing limit (b >> L), a ( z , z') reduces to I +2iz'/b". The first term in the integrand of (3.24) is the Gaussian distribution Gx:(r', z'), which can be taken out of the z-dependent integral. The generated harmonic field Eq is therefore Gaussian with a confocal parameter b" larger than that predicted in the weak-field limit. The focal section of the generated beam, equal to nb"/2kq = bU4p depends on the effective order p, rather than on the non-linear order q. The rest of the integral is quite close to (3.19), apart from the amplitude term, which varies approximately as (1+ 422/b2)"-"'2, instead of ( 1 + 4~~/ b~) (~-' " '~. When p" q, this term varies much less rapidly throughout the medium than in the perturbative limit. As discussed previously, phase matching around A k = 0 is determined by the variation of the polarization amplitude throughout the medium. It depends on the effective order of non-linearity (p) of the dipole moment d,,. Therefore it remains relatively constant for all the harmonics which have the same intensity dependence.
In a more general case, the generated harmonic field is not Gaussian. The situation is quite similar to difference-frequency processes (Bjorklund 1975). The term exp[-k,r"/b"a(z, z')] leads to a strong distortion o f t h e spatial profile which develops rings. From (3.24), one sees that phase matching can be favoured for some r' values, which may be different from zero (off axis). The number of these rings increases as the geometry becomes more focused. Typically, for b =4mm, there is a central spot and an external ring, whereas for b = 1 mm, up to four or five rings may appear (see figure 1 4 ( a ) ) . Although phase matching in this case cannot be reduced to a onedimensional problem along the propagation axis, the conclusion reached in the loose focusing case remains also valid: since phase matching depends on the variation of the amplitude of the polarization field throughout the medium which does not depend on the process order, it remains constant for all the harmonics in the plateau. Figure   15 compares IFJb, Ak)12 for the Sth, 13th and 21st harmonics in the weak-field limit (full curve) and for the 13th and 21st harmonics in the model situation where ldq12 is assumed to vary with the laser field strength as / E $ . At A k = O , IFq(b,0)l2 remains relatively constant for the model harmonics, whereas it decreases rapidly with q in the perturbative limit. This simple model, which can be carried out analytically for the most part, reproduces fairly well the main aspects of the numerical calculation using as a source the non-perturbative atomic dipole moment. It shows that the main difference between the strong-and weak-field regimes is the rate at which the magnitude of the polarization field varies within the medium. This has an important effect on the conversion efficiency because, in the experimental conditions, phase matching does not depend much on the phases of the interfering fields, but it does depend on the variation of the amplitudes of the fields throughout the medium. In a strong-field regime, the harmonics get defocused compared with the perturbative limit and, moreover, the phase matching is found to be relatively constant with increasing order.

Conclusion
As was pointed out in the introduction, these harmonic generation processes are difficult to grasp because of the interplay between the microscopic response (emission of harmonics by a single atom) and the collective one (propagation, phase matching of the propagating fields). In conclusion, we would like to summarize the role played by the many parameters of the problem, related either to the pump field or to the non-linear medium, influencing either the atomic response or propagation, or both. These must be viewed as a few tentative guidelines, most of them requiring a more thorough investigation, rather than definite conclusions.
Laser intensity. As the incident pump intensity increases, the plateau increases rapidly both in length and in amplitude, because these processes are highly non-linear. This remains true up to the intensity at which the atoms are depleted by ionization.
Laserpulse length. Shortening the pulse length increases the intensity at which the atoms are ionized. Consequently, the medium can produce more harmonics because it can experience a higher intensity. Note that the pulse length of generated radiation follows that of the pump field, being smaller by about a factor of two.
Laser wauelength. This aspect has not been discussed in this paper where we have focused on the results at 1 pm. On the basis of experiments at 248 nm (McPherson el a1 1987) and calculations of single atom spectra performed at 532 nm and 1064 nm (Potvliege andShakeshaft 1989b, Krause et al 1991), one might expect the harmonic conversion efficiency to increase for shorter incident wavelengths and the plateau to remain approximately of the same length in photon energy.
Laser confocal parameter. This provides a simple way of increasing the harmonic generation efficiency, which was found to vary as b', in the range b = 1-6 mm and L = 1 mm. Note that in the plane-wave limit the b-scaling should become less rapid, being approximately linear.
Atom. Atoms with higher ionization energies can produce more harmonics but with a reduced efficiency. A still open question is whether ions can efficiently produce harmonics.
Atomic density. The signal varies as N 2 independent of the order. This is true between 1 and 25Torr at 1 pm. Saturation effects (Rosman er a1 1988) are to be expected for shorter incident wavelengths and higher gas pressures.
Finally, let us emphasize the main conclusion reached in this paper. These highorder harmonic generation processes are twice the signature of a strong-field interaction with a non-linear medium, which cannot be described using weak-field approximations. The emission of radiation by a single atom has a non-perturbative behaviour: it exhibits the same plateau that was observed in the experiments; the harmonics do not vary as Iq, where q is the process order. Moreover, the plateau is conserved in the many-atom response, a fact that cannot be reconciled with a weak-field approach to the problem. It comes from the variation of the non-linear polarization with the laser intensity which is lower than in the weak-field limit and also similar for all the harmonics in the plateau. Extending harmonic generation into a strong-field regime has two unexpected and fortunate consequences: a plateau forms in the single-atom response, which means that generation of, say, the 17th harmonic is as probable as the generation of the 5th harmonic; and all the generated harmonics become equally and efficiently phase matched.
pari under the auspices of the Department of Energy at the Lawrence Livermore National Laboratory under contract number W-7405-Eng-48.
Nore added in proox The laser focal section obtained with the / = 200 mm lens has been recently measured and has been found ro be ahout a factor of two smaller than the one estimated previously. T h i s does not change any o f the conclusions of the paper but leads to a better agreement between the calculated saturation intensity and the experimental one (increased by a factor of two) and also hetween the calculated and measured numbers o f harmonic photons.