A Sävert

S P D Mangles

M Schnell

E Siminos

J M Cole

M Leier

M Reuter

M B Schwab

M Möller

K Poder

O Jäckel

G G Paulus

C Spielmann

S Skupin

Z Najmudin

M C Kaluza

Direct observation of the injection dynamics of a laser wakefield accelerator using few-femtosecond shadowgraphy

We present few-femtosecond shadowgraphic snapshots taken during the non-linear evolution of the plasma wave in a laser wakefield accelerator with transverse synchronized few-cycle probe pulses. These snapshots can be directly associated with the electron density distribution within the plasma wave and give quantitative information about its size and shape. Our results show that self-injection of electrons into the first plasma wave period is induced by a lengthening of the first plasma period. Three dimensional particle in cell simulations support our observations. Laser-wakefield accelerators (LWFA) operating in the 'bubble'-regime [1] can generate quasimonoenergetic multigigaelectronvolt electron beams [2,3] with femtosecond duration [4,5] and micrometer dimensions [6,7]. These beams are produced by accelerating electrons in laser-driven plasma waves over centimeter distances. They have the potential to be compact alternatives to conventional accelerators [8]. In a LWFA, the short driving laser pulse displaces plasma electrons from the stationary background ions. The generated space charge fields cause the electrons to oscillate and form a plasma wave in the laser's wake. This wave follows the laser at almost c, the speed of light; for low amplitude it has a wavelength of λ p = 2πc ε 0 m e /(n e e 2 ), (1) where n e is the electron density of the plasma. At high amplitude, electrons from the background can be injected into the wake and accelerated, producing monoenergetic electron pulses [9][10][11]. Significant progress has been made regarding achievable peak energy [3], beam stability [12] and the generation of bright X-ray pulses [13][14][15]. Until now, most of our knowledge about the dynamics of the self-injection process has been derived from detailed particle-in-cell (PIC) simulations. These simulations show that self-focusing [16] and pulse compression [17] play a vital role in increasing the laser pulse intensity prior to injection. Furthermore, simulations indicate that self-injection of electrons is associated with a dynamic lengthening of the first plasma wave's period (the 'bubble'). This lengthening can be driven by changes of the electric field structure inside the plasma wave caused by the injected electrons [18]. In contrast, the lengthening may also be due to an intensity amplification of the laser pulse caused by the non-linear evolution of the plasma wave [19,20] or due to a local increase in intensity caused by two colliding pulses [21]. In these latter scenarios, injection is a consequence of the lengthening of the bubble. However, experimental insight into these processes is extremely challenging due to the small spatial and temporal scales of a LWFA.

The plasma wave, a variation in the electron density, has an associated refractive index profile which can be detected using longitudinal [22][23][24] or transverse probes [5]. Longitudinal probes cannot measure the rapid and dynamic evolution of the plasma wave that occurs in nonlinear wakefield accelerators and suffer from the strong refraction caused by the steep refractive index gradients in a plasma wave. However, a small offset angle between pump and probe can measure the evolution, but only in 1D [25]. Furthermore, to properly resolve the wake structure the probe must have a duration (or equivalent bandwidth) less then λ p /c which has not yet been achieved [5,25], meaning that the important details of the wake evolution, e.g. the lengthening of the bubble in relation to the injection process have not yet been resolved. By using a probe pulse shorter than λ p /c to perform high resolution shadowgraphy we are able to show for the first time that, under our experimental conditions, bubble expansion occurs before self-injection starts.

In the present study, the JETI-laser system at the Institut für Optik und Quantenelektronik in Jena, Germany delivered pulses of 750 mJ energy and 35 fs duration (central wavelength λ L = 810 nm). The pulses were focused by an f /13 off-axis parabolic mirror to an elliptical focal spot with dimensions (FWHM) 8.9 µm×12.8 µm containing 27% of the energy, resulting in peak intensities of I L = 6 × 10 18 Wcm -2 (corresponding to a peak normalized vector potential of a 0 = 8.55 × 10 -10 λ L (µm) • I L (Wcm -2 ) ≈ 1.7). A supersonic helium gas jet was used, generating a plasma with electron density n e in the range of (0.5 . . . 2.5) × 10 19 cm -3 . Electrons accelerated during the interaction could be detected using a magnetic spectrometer or a scintillating screen. A small fraction of the laser was split from the main pulse, spectrally broadened in a hollow-core fiber filled with argon to support a transform limited pulse duration of τ FL = 4.4 fs. Using dispersive mirrors and glass wedges to optimize dispersion, probe pulses as short as τ probe = (5.9 ± 0.4) fs were created [26]. These synchronized, few-cycle probe pulses were used to back-light the LWFA perpendicularly to the pump-pulse direction. A high-resolution imaging system produced shadowgraphic images with micrometer resolution on a CCD camera. By varying the delay between pump and probe, different stages of the plasma wave's evolution were recorded on subsequent shots close to the threshold density for self injection. The snapshots shown in Fig. 1 are representative of each stage in the acceleration process. Shots were selected that exhibit similar quasi-monoenergetic electron spectra (Fig. 1(g-l)) and produced a high contrast shadowgram. The latter was affected by jitter in probe duration and pointing fluctuations of the pump which shifts the image out of focus. To reduce modulations induced by the probe pulse's beam profile, the relative intensity modulation was plotted using I norm = (I -I 0 )/I 0 , with I being the pixel value at each individual position and I 0 the value derived from a low order spline fit in the horizontal direction.

Early in the interaction the regions of high and low electron density gradient (dark and light regions in the image) are approximately equal in length, indicating a linear plasma wave, Fig. 1(a). The transverse extent of the plasma wave reduces and the amplitude of the wave increases, Fig. 1(b). Later on, significant increase of curvature of the plasma-wave train and, in particular, in the lengthening of the first plasma period, Fig. 1(c) is apparent. Just ahead of the region where this lengthening starts we observe bright emission from the plasma. This has a broad spectrum (covering at least a range from 600-1000 nm, cf. the 60 nm bandwidth of the pump) and is consistent with 'wavebreaking radiation' [27], which is a direct signature of the onset of self-injection in the experiment. Further propagation enhances the density gradient at the front of the bubble which now appears in the shadowgrams at the beginning of the wave train. After wavebreaking, the wave becomes highly non-linear, as indicated by the reversal in the direction of curvature of the trailing wave periods in the shadowgrams, Fig. 1(d) to 1(f). These features are closely linked to the process of transverse wavebreaking [28].

snapshots reveal that the dynamic process of bubble lengthening is intimately tied to self-injection. We plot in Fig. 2(a) the evolution of the plasma wave's first period. Early in the interaction the length of the first period has already increased as compared to the wavelength for a linear relativistic plasma wave, λ p = 2πc/ω p . The length of the bubble increases up to the point of wavebreaking, cf. Fig 1(c). During a single interaction, this radiation is emitted from a distinct spot on the optical axis with longitudinal position varying slightly in the experiment around (930±67) µm for n e ≈ 1.65×10 19 cm -3 . Beyond this point the shape of the plasma wave varies from shot to shot. This can lead to the formation of a single bubble or to the merging of the first two plasma wave periods due to beam loading resulting in a strong variation of the bubble length after injection (Fig. 2(a)). A linear regression fit to the data (Fig. 2(a)) shows that the bubble starts expanding from a diameter of (9.4±1.0) µm at an expansion velocity of v be = (2.4 ± 1.4) × 10 6 ms -1 . The fact that we observe bubble lengthening before in- jection clearly demonstrates that the initial expansion is not caused by the effect of charge loaded into the wake. This lengthening is therefore most likely caused by intensity amplification of the pulse as it propagates in the plasma wave.

We also measured the length of the second wave period, shown in Fig. 2(b), as a function of density at a fixed position in the plasma (v g t = 1.0 mm). These measurements were made sufficiently far into the gas jet to ensure that it was in the uniform density plateau. At low densities, the length is well matched to Eq. ( 1), but at high densities, λ p is significantly longer. The density at which this transition occurs corresponds to the onset of injected charge, see Fig. 2(c), and to the self-injection threshold predicted in Ref [29]. The second period lengthens due to the relativistic γ-factor of the electrons associated with the large amplitude plasma wave. As n e is increased beyond 1.7 × 10 19 cm -3 , λ p decreases but the rate of expansion increases so that we observed an approximately constant wavelength at this fixed position.

In addition, three dimensional particle in cell (3D-PIC) simulations were performed with the code EPOCH [30]. A laser pulse with τ L = 36 fs duration and λ L = 810 nm was focused to a spot size of 18.8 µm (FWHM), 300 µm into a plasma density profile similar to the experiment (peak density n e = 1.7 × 10 19 cm -3 ). In order to take into account imperfections in the experimental focal spot we set the maximum intensity of the laser pulse to I 0 = 2.5 × 10 18 Wcm -2 , leading to similar energies within the focal spot FWHM in experiment and simulation [31]. The computational domain was a 'sliding window' of size 150 × 70 × 70 µm 3 moving at c. We used 2700 × 525 × 525 cells with two electrons per cell and a stationary ion background. A 6 th order finitedifference-time-domain scheme was employed, together also been fully simulated in 3D with EPOCH. At different time steps during the driver pulse propagation the moving window was stopped and the probe was injected from the side of the box, propagating in the negative y direction, perpendicularly to the pump. The probe had a central wavelength λ probe = 750 nm, a Fourier limited duration of 4.4 fs and a negative linear chirp increasing its duration to 12 fs for a best match to the experimental images. The probe propagated past the wakefield structure, until y ≈ -15 µm. Subsequently, propagation in vacuum was assumed and modeled in Fourier space including the imaging system aperture, sensor sensitivity and image plane position. To adjust for the latter, we propagated the probe pulse backwards and recorded the time-averaged Poynting flux through the object plane, which was at y = 0. The comparison of density maps and PIC generated shadowgrams in Fig. 3 verifies that shadowgrams capture local variations in plasma density and allows a direct interpretation of the experimental shadowgrams. We note that the injected electron bunch appears neither in the simulated nor the experimental shadowgrams. This is caused by the reduction of local plasma frequency due to the γ-factor to the index of refraction. These simulated probe images confirm that few-femtosecond shadowgraphy provides quantitative information about the plasma wave including the plasma wave length, curvature and number of trailing periods. The length of the first plasma period, i.e. the bubble, as taken from the PIC simulations at various positions during the evolution is shown in Fig. 4 together with the maximum amplitude of the pump pulse's Poynting vector and the total injected charge with energy above 20 MeV. The intensity increase due to pump compression and self-focusing [see also Fig. 3(a)] is slow until approximately v g t = 800 µm. Until that point the bubble length grows from 1.2 λ p to 1.5 λ p . After v g t = 800 µm a phase of rapid intensity amplification begins, accompanied by bubble expansion and injection. At this stage, there is no substantial charge injected and so the bubble expansion is not due to beam-loading but due to intensity amplification of the pump and the associated increase of the γfactor of the plasma electrons; λ * p ≈ λ p (1+a 2 0 /2) 1/4 . This intensity amplification manifests itself in an increased visibility of the front of the bubble in the simulation, Fig. 3(b) and experimental (Fig. 1(e-f)) shadowgrams. In the PIC simulation significant charge is only injected into the wake (around v g t = 930 µm) after the length of the bubble has started to increase as also observed previously in simulations [9,19,32]. After v g t = 1000 µm, the intensity decreases, while the length of the bubble keeps increasing due to beam-loading. At approximately v g t = 1300 µm, the simulation shows a merging of the first two periods of the wake. The simulation supports the experimental observation that self-injection in the LWFA is caused by the expansion of the bubble. It predicts an expansion velocity of the bubble length of (4.2 ± 0.4) × 10 6 ms -1 (cf. Fig. 4), which is faster than the experimentally measured value. This could be due to imperfections in the experimental pulse profile leading to less efficient pulse self-compression and self-focusing than in the simulation. Our ability to determine the rate of expansion of the bubble from the experimental shadowgrams could allow the benchmarking and further development of dynamic bubble models in the future.

Applying our approach to visualize the full non-linear evolution of the plasma wave allows the acceleration process to be studied with unprecedented precision. As well as providing greater understanding of acceleration in the bubble-regime, our technique can easily be adopted to more complex acceleration geometries, e.g. staging [12], or for beam-driven acceleration [33][34][35]. Furthermore, increasing the probe pulse's wavelength while keeping a few-cycle duration will increase the sensitivity of our technique to probe low-density plasmas at the same relative resolution. Such low plasma densities are essential for high-energy plasma-acceleration scenarios. As LWFAs are widely expected to become useful sources of ultra short radiation [36], the increased level of understanding of plasma wave evolution and injection that can be gleaned from few-femtosecond shadowgraphy and the improvements in beam properties resulting from this technique are therefore likely to have a large impact on biomedical imaging and ultrafast condensed-matter study.

FIG. 1 .

 1 FIG.1. Left: gas/electron density profile and focus position (dashed line). a-f: Experimental shadowgrams at various positions in the plasma at a background electron density of ne = 1.65 × 10 19 cm -3 . The vertical lines indicate the on-axis position of the plasma wave's peaks as deduced from simulated shadowgrams (cf. Fig.3). g-l: Energy in the electron beam per MeV and spatially resolved in the vertical coordinate corresponding to the above images.

FIG. 2 .

 2 FIG. 2. (a) Length of the 1 st plasma period (bubble) as a function of propagation distance vgt taken from the shadowgrams [error bars represent combination of uncertainty in determining the bubble length (±0.5 µm) in a single shot reduced by √ N where repeat shots are available]. Blue (solid) horizontal line is the expected λp for ne = 1.65 × 10 19 cm -3 . Between vgt = 0.86 and 1.0 mm (grey shaded area) wavebreaking radiation was frequently detected. (b) Wavelength of 2 nd plasma period versus ne at a fixed position vgt = 1.0 mm and corresponding integrated electron bunch charge. Open circles represent measured and averaged data points with the standard error of the mean and red line is according to Eq. (1). (c) Electron beam profiles for different plasma densities averaged over 40-180 shots.

with 5 FIG. 3 .

 53 FIG. 3. (a) Electron density maps (gray) from the simulations at various positions and contours of the envelope of pump magnetic field Bz (red or solid lines) corresponding to 50% and 75% of the peak. (b) Shadowgrams simulated from the images in Fig. 3(a).

W/cm 2 )FIG. 4 .

 24 FIG.4. Bubble length derived from density maps (black squares) and injected charge with energy larger than 20 MeV (red or light gray) from the PIC simulation and maximum instantaneous value Smax of the Poynting vector magnitude (blue or dark grey). Black horizontal line is the expected λp for ne = 1.7 × 10 19 cm -3 .

We thank B. Beleites, W. Ziegler, and F. Ronneberger for running the JETI-laser system. E.S. and S.S. thank V. T. Tikhonchuk for helpful discussions. This study was supported by DFG (grants TR18 A12, B9 and KA 2869/2-1), BMBF (contracts 05K10SJ2 and 03ZIK052), European Regional Development Fund (EFRE), STFC (ST/J002062/1) and EP-SRC (EP/H00601X/1). The collaboration was funded by LASERLAB-EUROPE (grant agreement n • 284464, EC's Seventh Framework Programme). EPOCH was developed under UK EPSRC grants EP/G054940/1, EP/G055165/1 and EP/G056803/1. For the computations, the Max Planck supercomputer at the Garching Computing Center RZG was used.


Illia Thiele
email: illia-thiele@web.de

Stefan Skupin

Rachel Nuter

Modelling arbitrarily shaped and tightly focused laser pulses in electromagnetic codes

Keywords: electromagnetic codes, Maxwell solver, particle-in-cell (PIC) codes, tight focusing, vector beams

Investigation of laser matter interaction with electromagnetic codes requires to implement sources for the electromagnetic fields. A way to do so is to prescribe the fields at the numerical box boundaries in order to achieve the desired fields inside the numerical box. Here we show that the often used paraxial approximation can lead to unexpected field profiles with strong impact on the laser matter interaction results. We propose an efficient numerical algorithm to compute the required laser boundary conditions consistent with the Maxwell's equations for arbitrarily shaped, tightly focused laser pulses.

Introduction

Electromagnetic codes are useful tools to study various problems in microwave engineering, plasma physics, optics and other branches of natural science. Such codes solve Maxwell's equations coupled to constitutive equations describing the matter. In studies of laser matter interaction, external electromagnetic waves (the "laser") have to enter the computational domain in order to interact with the matter. In the case of particle-in-cell (PIC) codes like CALDER [START_REF] Lefebvre | Electron and photon production from relativistic laserplasma interactions[END_REF], PICLS [START_REF] Mishra | Collisional particle-in-cell modeling for energy transport accompanied by atomic processes in dense plasmas[END_REF] or OCEAN [START_REF] Nuter | Prepulse suppression and optimization of backward raman amplification with a chirped pump laser beam[END_REF], it is common practise to prescribe ex-arbitrarily shaped and polarized. Such laser pulses become more and more popular in the context of laser driven radiation and particle sources as well as laser material processing [START_REF] Yew | Second harmonic generation polarization microscopy with tightly focused linearly and radially polarized beams[END_REF][START_REF] Vieira | Nonlinear laser driven donut wakefields for positron and electron acceleration[END_REF][START_REF] Buccheri | Terahertz emission from laser-induced microplasma in ambient air[END_REF][START_REF] Cheng | A review of ultrafast laser materials micromachining[END_REF][START_REF] Hnatovsky | Polarization-dependent ablation of silicon using tightly focused femtosecond laser vortex pulses[END_REF][START_REF] Kammel | Enhancing precision in fs-laser material processing by simultaneous spatial and temporal focusing[END_REF][START_REF] Wang | Creation of a needle of longitudinally polarized light in vacuum using binary optics[END_REF].

The paper is organized as follows. Section 2 details the problem we want to solve. In Sec. 3, the theory of laser propagation in vacuum is reviewed. Section 4 describes in detail our algorithm for the computation of Maxwell consistent LBCs, and in Sec. 5 we present two illustrative examples: A tightly focused Gaussian beam and a longitudinal needle beam. Section 6 summarizes the results and offers perspectives for potential applications.

Schematic presentation of the laser injection

In numerical studies of laser matter interaction, it is common practise to define the laser by its propagation in vacuum, for example, by position and shape of the pulse at focus. In this paper, we choose to prescribe the pulse in a plane P parallel to a boundary of the rectangular numerical box, i.e., typically in the focal plane (see Fig. 1). The laser (red) is passing through the plane P, where the fields 1 E 0 , B 0 are prescribed for all times t. The goal is to calculate the fields E B , B B at the boundary from E 0 , B 0 . As we will see in Sec. 4, choosing P parallel to a boundary allows us to resort to Fast Fourier Transforms (FFTs) in the numerical computation of the LBCs. It is of course possible to prescribe the fields in an arbitrary plane and use the general solution given in the next Section to calculate the LBCs. However, in this case one cannot exploit the advantage of an efficient computation with FFTs (Sec. 4) and will have to evaluate the spatial Fourier integrals directly, for example by performing discrete sums.

Laser field propagation in vacuum

Let E 0 (r ⊥ , t) = E(r ⊥ , z = z 0 , t) and B 0 (r ⊥ , t) = B(r ⊥ , z = z 0 , t) be the electromagnetic fields in the plane P. In the following, we want to compute E, B in the whole space and for all times. We will see that not all components of E 0 , B 0 can be prescribed independently. Moreover, we will comment on how to handle evanescent fields, and finally discuss the paraxial limit.

1 Vectors are typed in bold.

∇ • Ê(r, ω) = 0 ∇ × Ê(r, ω) = iω B(r, ω) (1)
∇ • B(r, ω) = 0 ∇ × B(r, ω) = -iω 1 c 2 Ê(r, ω) . ( 2
)
Here, ω is the frequency variable, c is the vacuum speed of light, and ˆdenotes the Fourier transform with respect to time t. For the definition of the Fourier transforms as used in this paper see Appendix A. The wave equation for the electric field E in frequency space reads (analogue for the magnetic field B)

∆ Ê(r, ω) + ω 2 c 2 Ê(r, ω) = 0 . ( 3
)
Written in spatial Fourier space (with wavevector k as spatial Fourier variables)

Eq. ( 3) would reduce to an algebraic equation and lead to the vacuum dispersion relation k 2 = ω 2 /c 2 . However, we want to describe propagation of E 0 along z.

To this end, we keep the z variable and perform the Fourier transform with respect to the transverse variables r ⊥ only. Transforming Eq. (3) to transversal spatial Fourier space, where k ⊥ = (k x , k y ) T is the transversal wavevector, gives

k 2 z (k ⊥ , ω) Ē(k ⊥ , z, ω) + ∂ 2 z Ē(k ⊥ , z, ω) = 0 , (4)
where

k z (k ⊥ , ω) = ω 2 /c 2 -k 2
x -k 2 y , and ¯denotes the temporal and transverse spatial Fourier domain. The fundamental solutions of Eq. ( 4) are the forward (+) and backward (-) propagating, plane or evanescent waves (analogue for the magnetic field B)

Ē± (k ⊥ , z, ω) = Ē± 0 (k ⊥ , ω)e ±ikz(k ⊥ ,ω)(z-z0) . ( 5
)
It is important to note that E ± 0 , B ± 0 cannot be chosen arbitrarily. In fact, only two out of six vector components (for forward and backward direction, respectively) are independent. For example, we can choose to prescribe E ± 0,⊥ in the plane P. Then, by exploiting Eqs. ( 1) and ( 5), we get

Ē± ⊥ (k ⊥ , z, ω) = Ē± 0,⊥ (k ⊥ , ω)e ±ikz(k ⊥ ,ω)(z-z0) (6) Ē± z (k ⊥ , z, ω) = ∓ k ⊥ • Ē± ⊥ (k ⊥ , z, ω) k z (k ⊥ , ω) (7)
B± (k ⊥ , z, ω) = 1 ωk z (k ⊥ , ω) R ± (k ⊥ , ω) Ē± ⊥ (k ⊥ , z, ω) , (8)
with the matrix

R ± (k ⊥ , ω) =      ∓k x k y ∓ k 2 z (k ⊥ , ω) + k 2 y ± k 2 z (k ⊥ , ω) + k 2 x ±k x k y -k y k z (k ⊥ , ω) k x k z (k ⊥ , ω)      . ( 9
)
Obviously, we are imposing k z = 0, which is implicitly assumed when stating that the laser is passing through the plane P. Thus, the laser must not have any components propagating parallel to P. In complete analogy, one could prescribe the transverse magnetic fields B ± 0,⊥ in the plane P and exploit Eqs. [START_REF] Mishra | Collisional particle-in-cell modeling for energy transport accompanied by atomic processes in dense plasmas[END_REF] to compute B ± and E ± in the whole space. In Appendix B, we give an alternative method for computing of Maxwell consistent laser fields based on the vector potential in the Lorentz gauge. Such description can be advantageous in specific cases, for example radially polarized doughnut beams [START_REF] Salamin | Accurate fields of a radially polarized gaussian laser beam[END_REF], where only one component of the vector potential is sufficient to describe the whole laser.

Evanescent fields and the paraxial limit

For k 2 x + k 2 y > ω 2 /c 2 , k z (k ⊥ , ω
) becomes imaginary and Eq. ( 5) describes evanescent waves, with exponentially growing or decaying amplitude in z direction. In free space propagation, evanescent waves violate energy conservation and are thus unphysical and do not exist. In order to get rid of evanescent waves, the spatial Fourier spectrum of E 0 and B 0 has to be filtered in transverse spatial Fourier space, such that it contains only components with k 2

x + k 2 y < ω 2 /c 2 . This condition is nothing else then ensuring the Abbe diffraction limit [START_REF] Born | Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light[END_REF] for the fields prescribed at z = z 0 , which, for instance, forbids to focus a beam to arbitrary small transverse size.

In contrast, if the spatial Fourier spectrum of E 0 and B 0 is nonzero only for k 2

x + k 2 y ω 2 /c 2 , one can expand k z as a Taylor series and approximate

k z (k ⊥ , ω) ≈ |ω| c - c 2|ω| k 2 x + k 2 y . (10)
Then, Eqs. ( 6)-( 8) simplify as

Ē± ⊥ (k ⊥ , z, ω) ≈ Ē± 0,⊥ (k ⊥ , ω)e ±i[ |ω| c -c 2|ω| (k 2 x +k 2 y )](z-z0) (11) Ē± z (k ⊥ , z, ω) ≈ 0 B± x (k ⊥ , z, ω) ≈ ∓ 1 c Ē± y (k ⊥ , z, ω) (12)
B± z (k ⊥ , z, ω) ≈ 0 B± y (k ⊥ , z, ω) ≈ ± 1 c Ē± x (k ⊥ , z, ω) , (13)
which is well known as the paraxial or Fresnel approximation [START_REF] Goodman | Introduction to Fourier Optics[END_REF].

Implementing the laser boundary conditions

Let us now describe a practical implementation of LBCs based on the solution of Maxwell's equations as derived in the previous Section. In the following, the laser will propagate in forward direction (+) along z, i.e., we inject the laser from the left side of the box (see Fig. 1). We prescribe the electric field E 0,⊥ (r ⊥ , t) in the plane P at z = z 0 , for example a Gaussian profile in t and r ⊥ . Then, we want to calculate the fields E B (r ⊥ , t) and B B (r ⊥ , t) at the boundary z = z B on the numerical grid for all times. Let us consider an equidistant rectangular grid x i , y j , indices i, j running from 1 to N x , N y , respectively, and with spatial resolution δx, δy. We evaluate E 0,⊥ at the grid points x i , y j for equidistant times t n , n is running from 1 to N t , with temporal resolution δt:

E ijn 0,⊥ = E 0,⊥ (x i , y j , t n ) . ( 14
)
The following algorithm computes the electric and magnetic fields E ij B (t) and

B ij B (t) at the boundary z = z B for any given time t ∈ [t 1 -zB-z0 c , t Nt -zB-z0 c ]:
1. Calculate Êijn 0,⊥ via discrete Fourier transforms (DFTs) in time [START_REF] Press | Numerical Recipes in C (2Nd Ed.): The Art of Scientific Computing[END_REF]:

ω n = 2π N t δt - N t 2 + n (15) Êijn 0,⊥ = δt 2π Nt l=1 E ijl 0,⊥ e iω n t l , n = 1, . . . , N t . ( 16
)
2. Calculate Ēijn 0,⊥ via two-dimensional DFTs in transverse space:

k i x = 2π N x δx - N x 2 + i k j y = 2π N y δy - N y 2 + j (17) Ēijn 0,⊥ = δxδy (2π) 2 Nx,Ny l,m=1
Êlmn 0,⊥ e -i(k i x x l +k j y y m ) , i, j = 1, . . . , N x,y . ( 18)

3. Calculate transverse electric field components at the boundary (z = z B ):

k ijn z = (ω n ) 2 c 2 -(k i x ) 2 -(k j y ) 2 (19)
Ēijn B,⊥ =      Ēijn 0,⊥ e ik ijn z (zB-z0) for k ijn z > 0 0 for k ijn z = 0 . ( 20
)
Here, denotes the real part of a complex number. Note that we have set

75 k ijn z ≡ 0 and Ēijn B,⊥ ≡ 0 for indices i, j, n with (k i x ) 2 + (k j y ) 2 ≥ (ω n ) 2 /c 2 ,
in order to suppress evanescent waves (see Sec. 3.2).

Calculate the longitudinal electric field component at

z = z B : E ijn B,z =      - k i x E ijn B,x +k j y E ijn B,y k ijn z for k ijn z > 0 0 for k ijn z = 0 . ( 21
)
5. Calculate the magnetic field at z = z B :

R ijn =      -k i x k j y (k i x ) 2 -(ω n ) 2 c 2 (ω n ) 2 c 2 -(k j y ) 2 k i x k j y -k j y k ijn z k i x k ijn z      (22)
Bijn B =      1 ω n k ijn z R ijn Ēijn B,⊥ for k ijn z > 0 0 for k ijn z = 0 . ( 23
)
6. Calculate Êijn B and Bijn B via two-dimensional inverse DFTs: ]:

Êijn B,⊥ = (2π) 2 N x N y δxδy Nx,Ny l,m=1 Ēlmn B e i(k l x x i +k m y y j ) (24)
Bijn B,⊥ = (2π) 2 N x N y δxδy Nx,Ny l,m=1
E ij B,⊥ (t) = (2π) N t δt Nt n=1
Êijn B e -iω n t (26)

B ij B,⊥ (t) = (2π) N t δt Nt n=1 Bijn B e -iω n t . ( 27
)
The DFTs in steps 1, 2, and 6 can be calculated efficiently by means of FFTs.

There are various FFT libraries available, one of the most popular and efficient implementations is the FFTW [START_REF] Frigo | Fastest Fourier Transform in the West[END_REF]. One has to take into account the particular definitions of spatial and temporal Fourier transform used in this paper (see Appendix A), as well as the conventions of the particular FFT library. For the FFTW [START_REF] Frigo | Fastest Fourier Transform in the West[END_REF], one has to use the forward transform (flag FFTW FORWARD) in step 2, and the backward transform (flag FFTW BACKWARD) in steps 1 and 6.

The Fourier sums in step 7 allow to compute E ij B (t) and B ij B (t) for any given time t by means of discrete Fourier interpolation. In fact, most of the discrete frequencies ω n will have a negligible contribution to the spectrum when we are dealing with not-too-short laser pulses, i.e., a pulse envelope modulated with the centre frequency ω c (see Fig. 2). By taking only the significant summands into account when evaluating the Fourier sums Eqs. ( 26) and ( 27) reduces significantly both memory consumption and execution time.

When using DFTs to approximate continuous Fourier transforms as in the proposed algorithm above, one has to be careful with respect to sampling rates and the inevitable periodic boundary conditions. The initial datum E 0,⊥ has to be well resolved in space and time, and one has to check that the beam fits well in the transverse numerical box for all relevant z (e.g., the beam width may be larger at the boundary z = z B due to diffraction). Finally, one should not forget that Eqs. ( 26) and ( 27) should be evaluated for times t in the interval

[t 1 -zB-z0 c , t Nt -zB-z0 c
] only, otherwise a pulse train will be injected due to periodicity in time.

In a practical implementation, steps 1-6 will be performed by a pre-processor before launching the main simulation. Then, only the relevant (nonzero) contents of the arrays Êijn B and Bijn B (see remark above) will be passed to the main code and step 7 will be calculated at each time step of the main simulation.

Before going on with examples, we want to make a last remark concerning the grid structure of particular Maxwell solvers. For solvers like the "Directional Splitting scheme" [START_REF] Nuter | Maxwell solvers for the simulations of the laser-matter interaction[END_REF], E and B are discretized on the same equidistant grid and the above algorithm can be applied directly. For other solvers, like the "Yee scheme" [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwells equations in isotropic media[END_REF], the fields are described on grids shifted by δ x /2, δ y /2, δ z /2, respectively. In such case, a straight forward work around would be to run the pre-processor several times with transversely shifted grids and/or shifted boundary, in order to compute the desired field components for laser injection.

Examples

Tightly focused Gaussian pulse

Tightly focused pulses are potentially interesting for various kinds of experiments giving the possibility to achieve high intensities at rather low pulse energy or to generate micro-plasmas. Here, we are going to simulate a tightly focused Gaussian pulse and its interaction with an initially neutral gas, that is going to be ionized during the interaction. The electromagnetic fields resulting from LBCs in paraxial approximation Eqs. ( 11)-( 13), as they are often applied in PIC codes, will be compared with LBCs according to the Maxwell consistent approach Eqs. ( 6)- [START_REF] Yew | Second harmonic generation polarization microscopy with tightly focused linearly and radially polarized beams[END_REF]. For sake of computational costs, we restrict ourselves to the two dimensional case, where ∂ y ≡ 0 accounts for translational invariance in transverse y direction. For both cases a linear polarized Gaussian pulse is prescribed in the focal plane z = z 0 by We solve Maxwell's equations numerically using the PIC code OCEAN [START_REF] Nuter | Prepulse suppression and optimization of backward raman amplification with a chirped pump laser beam[END_REF]. In all simulations we consider an argon atmosphere at ambient pressure. The code OCEAN fully accounts for ionization according to the quasistatic ADK theory [START_REF] Nuter | Field ionization model implemented in particle in cell code and applied to laseraccelerated carbon ions[END_REF][START_REF] Ammosov | Tunnel ionization of complex atoms and of atomic ions in an alternating electric field[END_REF][START_REF] Yudin | Nonadiabatic tunnel ionization: Looking inside a laser cycle[END_REF] and uses ionization data from [START_REF] Carlson | Calculated ionization potentials for multiply charged ions[END_REF]. It is thus instructive to inspect the electron plasma generated by the tightly focused laser pulses for paraxial and Maxwell consistent LBCs. The resulting distributions of the

E 0,⊥ (x, t) = E 0 e -x w 0 2 -t t 0 2 cos(ω c t)e x , (28)
(d) (e) (f) (a) (b) (c)

Longitudinal needle beam

In order to demonstrate generality and ease of use of the proposed Maxwell consistent LBCs, let us have a look at a (on the first glance) more complicated example. In [START_REF] Wang | Creation of a needle of longitudinally polarized light in vacuum using binary optics[END_REF], the authors describe the "creation of a needle of longitudinally polarized light" by tight focusing of a radially polarized Bessel-Gaussian beam.

The radial component of the electric field of such beam at focus reads

E 0,⊥ (r, t) = α 0 T (θ) √ cos θ sin(2θ)e -( sin θ sin α ) 2 J 1 2 sin θ sin α J 1 ω c c r sin θ dθ × E 0 cos(ω c t)e r . (29)
Here, the electric field is written in cylindrical coordinates (r, φ, z), and e r is the radial unit vector. The beam profile is given as an integral over the angle θ, where α denotes the acceptance angle of the focusing optic. Following [START_REF] Wang | Creation of a needle of longitudinally polarized light in vacuum using binary optics[END_REF],

we assume a numerical aperture NA = 0.95, corresponding to α ≈ 0.4π. J 1 (x) denotes the corresponding Bessel function. The transmission function T (θ) takes into account a binary-phase optical element, which may further increase the relative longitudinal field strength as well as the length of the needle, however, to the detriment of the optical efficiency. Here, we consider a five-belt optical element [START_REF] Wang | Creation of a needle of longitudinally polarized light in vacuum using binary optics[END_REF] T

(θ) =      1 for 0 ≤ θ < θ 1 , θ 2 ≤ θ < θ 3 , θ 4 ≤ θ < α -1 for θ 1 ≤ θ < θ 2 , θ 3 ≤ θ < θ 4 , (30)
with θ 1 = 0.0275π, θ 2 = 0.121π, θ 3 = 0.19π, and θ 4 = 0.26π. As in the previous example, we consider a laser wavelength of λ c = 0.8 µm.

Figure 5 presents radial and longitudinal electric fields of the longitudinal needle beam from simulations using OCEAN [START_REF] Nuter | Prepulse suppression and optimization of backward raman amplification with a chirped pump laser beam[END_REF] and Maxwell consistent LBCs.

160

In agreement with [START_REF] Wang | Creation of a needle of longitudinally polarized light in vacuum using binary optics[END_REF] we find a longitudinal field amplitude that exceeds the radial one in the focal region along several laser wavelengths (∼ 8λ c ). The maximum longitudinal field amplitude is about 1.6 times larger than the radial one, which achieves its maximum out of focus at z = ±4λ c . This allows the longitudinal field to dominate in the focal plane by a factor of 2.5. 165

Conclusion

Injecting laser pulses into Maxwell solvers requires to prescribe the electromagnetic fields at the boundaries of the numerical box. Often, these fields are calculated by using the paraxial approximation. We have shown that for tightly focused beams this approach does not give the expected results. Instead, Maxwell's equations in vacuum have to be solved rigorously in order to find the proper fields at the boundaries. We proposed an easy to implement algorithm to achieve this goal, which allows to calculate the LBCs from transversal electric or magnetic field components defined in a plane, e.g., the focal plane. The presented algorithm can be parallelized in a straight forward manner and may be used with simulations tools employing domain decomposition.

We successfully employed our approach to simulate a tightly focused Gaussian pulse. An accurate handling of the laser injection turns out to be crucial:

Electron density profiles from ionization of neutral argon atoms due to field ionization are shown to be strongly dependent on the LBCs. Consequently, the LBCs may have significant impact on features like back-reflected radiation or energy deposition in the medium. Furthermore, our algorithm offers a simple way to simulate more complex pulse configurations or even sampled experimental beam profiles. Such "structured light" receives a lot of recent interest from various communities [START_REF]Structured light in structured media: From classical to quantum optics[END_REF]. As an example we demonstrated a longitudinal needle beam, which may be interesting for, among others, laser based material processing or particle acceleration studies. Thus, we believe that our approach will be useful for a larger community working on electromagnetic simulation codes.

By plugging Eq. (B.4) into Eq. (B.1), and using Eq. (B.2) to eliminate φ, electric and magnetic fields can be expressed in terms of the vector potential at z = z 0 : B± (k ⊥ , z, ω) = ik ± (k ⊥ , ω) × Ā± 0 (k ⊥ , ω)e ±ikz(k ⊥ ,ω)(z-z0) (B.5) Ē± (k ⊥ , z, ω) = iω 1 -c 2 ω 2 k ± (ω)k ± (ω) T Ā± 0 (k ⊥ , ω)e ±ikz(k ⊥ ,ω)(z-z0) . (B.6)

In general, the three components of A ± 0 can be chosen independently, however, only two components are necessary to prescribe an arbitrary laser pulse 2 . The use of the vector potential can be nevertheless advantageous, because certain beams, like radially polarized doughnut beams [START_REF] Salamin | Accurate fields of a radially polarized gaussian laser beam[END_REF], can be described by a single (longitudinal) component of the vector potential.

Figure 1 :

 1 Figure 1: Schematic picture of the laser (red) injection problem into the computational domain: Electric and magnetic fields E 0 , B 0 are prescribed in the plane P [here the (x, y)-plane at z = z 0 ]. The fields E B , B B at the boundary (blue) are unknown and have to be calculated.

3 3. 1 .

 31 Propagation of electromagnetic fields and their interdependenciesElectromagnetic fields in vacuum are governed by Maxwell's equations. In frequency or temporal Fourier space they read

Blmn B e i(k l x x i +k m y y j ) . ( 25 ) 0 Figure 2 :

 2502 Figure 2: Sketch of the electric field amplitude for a multi-cycle laser pulse in frequency domain. The spectrum is significantly different from zero only in Nω Nt frequency points.

  with center wavelength 2πc/ω c = λ c = 0.8 µm, pulse duration t 0 = 20 fs, peak intensityI 0 = 0 c|E 0 | 2 /2 = 5 × 10 14 W/cm2 giving E 0 = 61.4 GV/m and beam width w 0 = 0.35 µm. The particular choice of the beam width w 0 implies that non-negligible parts of Ē0,⊥ (k x , ω) are evanescent. These evanescent fields are suppressed in the calculation of E B (r ⊥ , t) and B B (r ⊥ , t) at the boundary z = z B fully compatible with Abbe's diffraction limit (see Sec.3.2). This leads to a 10% larger full-width-at-half-maximum (FWHM) beam width and smaller electric field at focus.

Figure 3 compares

 3 Fig.3(d) shows no side-wings in the beam profile. The maximum transversal electric field amplitude for the paraxial LBCs is significantly lower than that achieved with the Maxwell consistent LBCs. For both LBCs, the longitudinal field amplitude reaches about 30% of the transversal field amplitude, a direct consequence of the steep transversal gradients in the beam profile.

Figure 3 :

 3 Figure 3: Comparison of LBCs in paraxial approximation Eqs. (11)-(13) (a-c) and according to the Maxwell consistent approach Eqs. (6)-(8) (d-e). Snapshots of transversal fields Ex (b,e) and longitudinal fields Ez (c,f) of a tightly focused Gaussian pulse (see text for details) reveal strong distortions in case of the paraxial LBCs. Calculations were performed using the PIC code OCEAN [3], assuming an argon atmosphere at ambient pressure. In (a) and (d) line-outs of the transversal electric field Ex at focus are presented, revealing strong side-wings in the beam profile for the paraxial LBCs. The laser propagates from left to right.

Figure 4 :

 4 Figure 4: Electron densities ne produced by the tightly focused Gaussian laser pulses shown in Fig. 3 (see text for details). The profile produced by paraxial LBCs (a) is even qualitatively different than the one produced by Maxwell consistent LBCs (b). Electron densities are scaled to the initial neutral density n 0 . The laser propagates from left to right.

Figure 5 :

 5 Figure 5: The absolute values of the radial and longitudinal electric fields Er (a) and Ez (b) of a longitudinal needle beam in the focal region. The fields are normalized to the maximum Er,max of the radial field Er in the whole space.

Acknowledgements

Numerical simulations were performed using computing resources at Mésocentre de Calcul Intensif Aquitain (MCIA), Grand Equipement National pour le Calcul Intensif (GENCI, grant no. 2015-056129), and Partnership for Advanced Computing in Europe (PRACE, grant no. 2014112576).

Appendix A. The Fourier transforms

We define the temporal Fourier transform f (r, ω) of a function f (r, t) by

Further one, we define the transverse spatial Fourier transform f (r ⊥ , z, ω) of a function f (r, ω) by

where r ⊥ = (x, y) T and k ⊥ = (k x , k y ) T .

195

Note the difference in the sign of the exponent for temporal and spatial transform, which is common practise in the optical context. In particular when one wants to approximate Fourier integrals by finite sums, and resort to discrete Fourier transformations (DFTs) or even fast Fourier transforms (FFTs) [START_REF] Press | Numerical Recipes in C (2Nd Ed.): The Art of Scientific Computing[END_REF], it is important to keep track of these sign conventions (see Sec. 4).

Appendix B. Generating Maxwell consistent solutions using the vector potential in Lorentz gauge

Introducing electromagnetic potentials A, φ in Lorrentz gauge via

leads to decoupling of φ and the components of A, and one finds (in vacuum) [START_REF] Jackson | Classical electrodynamics[END_REF] k

In analogy to Eq. ( 4), fundamental solutions are the forward (+) and backward (-) propagating, plane or evanescent waves Ā± (k ⊥ , z, ω) = Ā± 0 (k ⊥ , ω)e ±ikz(k ⊥ ,ω)(z-z0) . (B.4)