Olivier Roussel
email: rousselo@coria.fr

Kai Schneider

Adaptive multiresolution computations applied to detonations

Keywords: multiresolution, adaptivity, wavelet, finite volume, detonation

A space-time adaptive method is presented for the reactive Euler equations describing chemically reacting gas flow where a two species model is used for the chemistry. The governing equations are discretized with a finite volume method and dynamic space adaptivity is introduced using multiresolution analysis. A time splitting method of Strang is applied to be able to consider stiff problems while keeping the method explicit. For time adaptivity an improved Runge-Kutta-Fehlberg scheme is used. Applications deal with detonation problems in one and two space dimensions. A comparison of the adaptive scheme with reference computations on a regular grid allow to assess the accuracy and the computational efficiency, in terms of CPU time and memory requirements.

Introduction

Real world industrial or environmental problems, e.g., management of industrial risks, typically involve physical and chemical phenomena having a multitude of dynamically active spatial and temporal scales. Their direct numerical modelling thus leads to prohibitive computational cost. Introducing adaptivity can be understood in the sense that the computational effort is concentrated at locations and time instants where it is necessary to ensure a given numerical accuracy, while efforts may be significantly reduced elsewhere.

Adaptive methods are in many cases more competitive than schemes on regular fine grids, in particular for solutions of nonlinear PDEs exhibiting a non-uniformly distributed regularity of the solution. Reliable error estimators of the solution are essential ingredients of fully adaptive schemes. They are based for example on Richardson ideas of extrapolation, adjoint problems or gradient based approaches. For evolutionary problems, a major task is the time evolution of the grid and its reliable prediction for the next time step. However, to become efficient, adaptive methods require a significant effort on implementing data structures, which are typically based on graded trees, hash-tables or multi-domains. Moreover, the computational cost per cell is significantly increased with respect to uniform discretizations. Hence, an adaptive method is only efficient when the data compression is large enough to compensate the additional computational cost per cell. Fortunately, for problems exhibiting local discontinuities or steep gradients, adaptive computations are faster than fine grid computations.

Adaptive discretization methods for solving nonlinear PDEs have a long tradition and can be tracked back to the late seventies [START_REF] Brandt | Multilevel adaptive solutions to boundary value problems[END_REF]. Adaptive finite element methods have a long history, in particular for elliptic problems. For chemically reactive flow with detailled chemical reaction in three dimensions [START_REF] Braack | Stabilized finite elements for 3D reactive flows[END_REF][START_REF] Braack | Solving multidimensional reactive flow problems with adaptive finite elements[END_REF] proposed stabilized finite elements with adaptive mesh refinement. The equations are treated fully coupled with a Newton solver and the solution of large linear non-symmetric, indefinite systems becomes necessary for which a parallel multigrid solver is used. Moving grid techniques have been applied successfully to combustion problems [START_REF] Hyman | An adaptive moving mesh method with static rezoning for partial differential equations[END_REF]. A posteriori error estimators have also been studied for a long time to improve the grid, since the early work of Babuska and Rheinboldt [START_REF] Babuska | A posteriori error estimates for the finite element method[END_REF]. However, adjoint problems have to be solved which are linear although the original PDE can be nonlinear [START_REF] Becker | An optimal control approach to a-posteriori error estimation[END_REF]. Fully adaptive finite element discretizations of reaction-diffusion problems encountered in electrocardiology have been proposed in [START_REF] Lang | Two-dimensional fully adaptive solutions of reaction-diffusion equations[END_REF][START_REF] Franzone | Adaptivity in space and time for reaction-diffusion systems in electrocardiology[END_REF]. For time adaptivity a stepsize control with linearly implicit time integrators is used. In space a multilevel finite element method is combined with a posteriori local error estimators.

The main challenge is to estimate and control the error of adaptive schemes with respect to the exact solution, or at least with respect to the same numerical scheme on an underlying uniform grid. Self adaptive methods are preferred as they automatically adjust to the solution. The block-structured adaptive mesh refinement technique (AMR or SAMR) for hyperbolic partial differential equations has been pioneered by Berger and Oliger [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF]. While the first approach utilized rotated refinement grids that required complicated conservative interpolation operations, AMR denotes today especially the simplified variant of Berger and Collela [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF] that allows only refinement patches aligned to the coarse grid mesh. The striking efficiency of this algorithm, in particular for 3D instationary supersonic gas dynamics problems, has been demonstrated by Berger et al. in [4].

Recently, multiresolution (MR) techniques have become popular for hyperbolic conservation laws, going back to the seminal work of Harten [START_REF] Harten | Multiresolution algorithms for the numerical solution of hyperbolic conservation laws[END_REF] in the context of finite volume schemes and cell-average MR analysis. Starting point is a finite volume scheme for hyperbolic conservation laws on a regular grid. Subsequently a discrete multiresolution analysis is used to avoid expensive flux computations in smooth regions, first without reducing memory requirements, e.g. for 1D hyperbolic conservation laws (Harten [30]), 1D con-servation laws with viscosity (Bihari [START_REF] Bihari | Multiresolution schemes for conservation laws with viscosity[END_REF]), 2D hyperbolic conservations laws (Bihari and Harten [START_REF] Bihari | Multiresolution schemes for the numerical solution of 2-D conservation laws I[END_REF]), 2D compressible Euler equations (Chiavassa and Donat [START_REF] Chiavassa | Point value multiscale algorithms for 2D compressible flow[END_REF]), 2D hyperbolic conservation laws with curvilinear patches (Dahmen et al [START_REF] Dahmen | Multiresolution schemes for conservation laws[END_REF]) and unstructured meshes (Abgrall and Harten [START_REF] Abgrall | Multiresolution representation in unstructured meshes[END_REF], Cohen et al [START_REF] Cohen | Multiresolution finite volume schemes on triangles[END_REF]). A fully adaptive version, still in the context of 1D and 2D hyperbolic conservation laws, has been developed to reduce also memory requirements (Gottschlich-Müller and Müller [START_REF] Gottschlich-Müller | Adaptive finite volume schemes for conservation laws based on local multiresolution techniques[END_REF], Kaibara and Gomes [START_REF] Kaibara | A fully adaptive multiresolution scheme for shock computations[END_REF], Cohen et al [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF]). This algorithm has been extended to the 3D case and to parabolic PDEs (Roussel et al [START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF], Roussel and Schneider [START_REF] Roussel | An adaptive multiresolution method for combustion problems: application to flame ball-vortex interaction[END_REF]), and more recently to self-adaptive global and local time-steppings (Müller and Stiriba [START_REF] Müller | Fully adaptive multiscale schemes for conservationlaws employing locally varying time stepping[END_REF], Domingues et al [START_REF] Domingues | An adaptive multiresolution scheme with local time stepping for evolutionary PDEs[END_REF][START_REF] Domingues | Space-time adaptive multiresolution methods for hyperbolic conservation laws: Applications to compressible Euler equations[END_REF][START_REF] Domingues | An adaptive multiresolution method for parabolic PDEs with time-step control[END_REF]). Therewith the solution is represented and computed on a dynamically evolving automatically adapted space-time grid. Different strategies have been proposed to evaluate the flux without requiring a full knowledge of fine grid cell-average values. Applications to shock waves in compressible flows addressing the issue of shock resolution have been presented in [START_REF] Bentaleb | Adaptive multiresolution methods for the simulation of shocks/shear layer interaction in confined flows[END_REF], and extensions to the Navier-Stokes equations in the weakly compressible regime can be found in [START_REF] Roussel | Coherent Vortex Simulation of weakly compressible turbulent mixing layers using adaptive multiresolution methods[END_REF]. Adaptive MR methods with operator splitting have been proposed for multiscale reaction fronts with stiff source terms in [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multi-resolution and dedicated high order implicit/explicit time integrators[END_REF][START_REF] Duarte | Adaptive numerical methods in time and space for the simulation of multiscale reaction fronts[END_REF]. The numerical analysis of the above higher order operator splitting techniques has been performed in [START_REF] Descombes | Adaptive time splitting method for multi-scale evolutionary PDEs[END_REF][START_REF] Duarte | Parareal operator splitting techniques for multi-scale reaction waves : numerical analysis and strategies[END_REF] and it was shown that the splitting time step can be even larger than the times scales involved in the PDEs. Adaptive MR computations using the above method with complex chemistry and including detailed transport can be found in [START_REF] Descombes | Analysis of operator splitting in the non-asymptotic regime for nonlinear reaction diffusion equations. Application to the dynamics of premixed flames[END_REF], further applications involving various stiffness levels have been presented in [START_REF] Duarte | Time-space adaptive numerical methods for the simulation of combustion fronts[END_REF] [START_REF] Duarte | A new numerical strategy with space-time adaptivity and error control for multi-scale streamer discharge simulations[END_REF].

The MR approach has also been used in other contexts. For instance, the Sparse Point Representation (SPR) method was the first fully adaptive MR scheme, introduced by Hölmstrom [START_REF] Holmström | Solving hyperbolic PDEs using interpolating wavelets[END_REF] in the context of finite differences and point-value MR analysis, leading to both CPU time and memory reduction. In the SPR method, the wavelet coefficients are used as regularity indicators to create locally refined grids, on which the numerical solution is represented and the finite difference discretization of the operators is performed. Applications of the SPR method have been published in [START_REF] Domingues | Adaptive wavelet representation and differenciation on block-structured grids[END_REF][START_REF] Pinho | Interpolating wavelets and adaptive finite difference schemes solving Maxwell's equations: gridding effect[END_REF]. Discontinuous Galerkin methods, they have been applied to hyperbolic conservation laws in [START_REF] Calle | Wavelets and adaptive grids for the discontinuous Galerkin method[END_REF] using Haar wavelet indicators to decide where to refine or to coarsen the meshes. These publications reveal that the multiresolution concept has been applied by several groups with success to different stiff problems. For comprehensive literature about the subject, we refer to the books of Cohen [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF] and Müller [START_REF] Müller | Adaptive multiscale schemes for conservation laws[END_REF].

The objective of the paper is the extension of the adaptive multiresolution method [START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF][START_REF] Domingues | Space-time adaptive multiresolution methods for hyperbolic conservation laws: Applications to compressible Euler equations[END_REF] to the numerical simulation of detonation waves. As model we use here a one-step chemical reaction involving two chemical species only. Since the chemical source term is stiff, a time splitting has to be made between the convective and source terms and each term needs to be computed with a different time step and a different time integration method. An error-controled time step based on a Runge-Kutta-Fehlberg approach is used for the source term integration. The application concerns Chapman-Jouguet detonations in one dimension with different stiffness values and instabilities of detonations waves due to an interaction with a pocket of partially burnt gases in two dimensions. These detonation problems motivate the use of the Euler equations, although extensions of the method to compute viscous flows using the Navier-Stokes equations is straightforward, as proposed in [START_REF] Roussel | Coherent Vortex Simulation of weakly compressible turbulent mixing layers using adaptive multiresolution methods[END_REF].

The outline of the paper is the following: first, we present the set of reactive Euler equations for a simplified detonation model. Then we describe the finite volume discretization. A Strang splitting technique is utilized to account for temporal scales in the source term that do not influence the hydrodynamics. We also briefly summarize the multiresolution strategy. Finally, we show the numerical results for the test problems in one and two space dimensions and we set the conclusions together with some perspectives for future work.

Governing equations

For modelling the combustion process, we use the reactive Euler equations, as described in [START_REF] Courant | Supersonic flow and shock waves[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]. The simplest description of a chemically reacting gas flow assumes that the gas mixture is made only of two chemical species, the burnt gas, denoted with subscript b and the unburnt gas, denoted with subscript u. The unburnt gas is converted to burnt gas via a single irreversible reaction. We represent the mixture state by a single scalar variable Z corresponding the mass fraction of the unburnt gas. We also assume gases in the mixture to be ideal polytropic gases with equal specific heat ratio γ and specific gas constant r. The system of equations in two dimensions, which has been non-dimensionalized in a suitable way, may be written as

∂Q ∂t + ∂F ∂x + ∂G ∂y = S, (1)
where Q = (ρ, ρv x , ρv y , ρe, ρZ) T and

F =         ρv x ρv 2 x + p ρv x v y (ρe + p)v x ρv x Z         , G =         ρv y ρv x v y ρv 2 y + p (ρe + p)v y ρv y Z         , S =         0 0 0 0 -K(T )ρZ         (2)
Here ρ denotes the mixture density, V = (v x , v y ) T the mixture velocity, e the mixture total energy per unit of mass, p the pressure, T the temperature and k the chemical reaction rate. The two equations of state completing the model are

p = ρrT (3)
and

e = p ρ(γ -1) + V 2 2 + Q 0 Z (4)
where Q 0 denotes the amount of heat per unit of mass released in the chemical reaction. The reaction rate k(T ) of the irreversible chemical reaction is expressed in Arrhenius form as

k(T ) = A exp - T A T (5)
where the pre-exponential coefficient A and the activation temperature T A are empirical constants. When the reaction source term is stiff, however, the reaction rate may be simplified by adopting the so-called ignition temperature kinetic model, i.e.

k(T ) = 1 τ if T ≥ T i 0 if T < T i (6)
where T i denotes the ignition temperature and τ the characteristic time of the chemical reaction, which determines the stiffness of the problem. This formulation has been chosen in the applications presented in the numerical results section. However, the numerical method is not limited to this simplified model and its extension to the Arrhenius law or even more complex chemical reactions is possible.

Numerical method

In this section, we first describe the classical Strang splitting and the space discretization of the convective terms. Subsequently, the time integration is discussed and, first, for the convective terms, a time step depending on the CFL is chosen, then, for the source term, we split the convective time step into error-controlled time steps using an explicit Dormand-Prince method. A different choice has been proposed by Duarte et al in [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multi-resolution and dedicated high order implicit/explicit time integrators[END_REF], based on implicit and explicit Runge-Kutta methods and a posteriori error estimators. We also briefly recall the multiresolution method, previously published in [START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF][START_REF] Domingues | An adaptive multiresolution scheme with local time stepping for evolutionary PDEs[END_REF].

Strang splitting

We denote by C(Q) the operator of the convective terms. Equation (1) becomes

∂Q ∂t = C(Q) + S(Q). (7)
Discretizing explicitly with first order in time, we get

Q n+1 = Q n + ∆t [C(Q n ) + S(Q n )] (8)
where n denotes the time instant and ∆t the convective time step.

The splitting relies on the separation of the convective and source terms operators. To get a first-order discretization, it writes

Q = Q n + ∆t S(Q n ) (9) Q n+1 = Q + ∆t C(Q ) (10)
i.e.

Q n+1 = C o1 ∆t S o1 ∆t Q n ( 11
)
where S o1 ∆t denotes the first-order accurate source term operator with time step ∆t and C o1 ∆t the first-order accurate convection operator with time step ∆t.

To get second-order accuracy, which corresponds to the classical Strang, the following procedure can be applied

Q n+1 = S o2 ∆t/2 C o2 ∆t S o2 ∆t/2 Q n ( 12
)
where S o2 ∆t/2 denotes the second-order accurate source term operator with time step ∆t/2 and C o2 ∆t the second-order accurate convection operator with time step ∆t. We note that the splitting time step in the above method is fixed. Techniques to introduce adaptive splitting time steps based on a posteriori error estimators have been introduced in [START_REF] Descombes | Adaptive time splitting method for multi-scale evolutionary PDEs[END_REF][START_REF] Duarte | Adaptive numerical methods in time and space for the simulation of multiscale reaction fronts[END_REF] allowing automatic error control.

Space discretization of the convective terms

Discretization in space is made using a finite volume method, to ensure conservative flux computations. Convective terms are discretized using the AUSM+ scheme [START_REF] Liou | A sequel to AUSM: AUSM+[END_REF]. In this procedure, pressure terms are computed separately. The Euler flux F writes

F =         ρv x ρv 2 x + p ρv x v y (ρe + p)v x ρv x Z         = M c         ρ ρv x ρv y ρh ρZ         +         0 p 0 0 0         (13)
where M denotes the Mach number, c the speed of sound and h the enthalpy per unit of mass. Denoting by Φ the purely convective term and by Π the pressure term, we discretize the Euler flux F in space following

F i+ 1 2 = M i+ 1 2 c i+ 1 2 Φ i+ 1 2 + Π i+ 1 2 ( 14
)
where the indices refer to nodal values.

The interface speed of sound is c i+ 1 2 = √ c i c i+1 and the interface convective term is

Φ i+ 1 2 = Φ i if M i+ 1 2 ≥ 0 Φ i+1 otherwise (15)
The terms M i+ 1 2 and p i+ 1 2 follow

M i+ 1 2 = M + i + M - i+1 ( 16
)
p i+ 1 2 = P + i p i + P - i+1 p i+1 (17)
where

M ± i = 1 2 (M i ± |M i |) if |M i | ≥ 1 ± 1 2 (M i ± 1) 2 ± 1 8 (M 2 i -1) 2 otherwise ( 18
)
and

P ± i = 1 2 (1 ± sign(M i )) if |M i | ≥ 1 1 4 (M i ± 1) 2 (2 ∓ M i ) ± 3 16 M i (M 2 i -1) 2 otherwise (19)
Third-order accuracy of the spatial discretization far from discontinuities is obtained using a MUSCL interpolation [START_REF] Van Leer | Towards the ultimate conservative difference scheme V: a second order sequel to Godunov's method[END_REF], together with a Koren slope limiter [START_REF] Koren | A robust upwind discretisation method for advection, diffusion and source terms[END_REF]. Denoting by q one of the conservative quantities (i.e. density, momentum, energy, partial mass of the unburnt gas), the corrected value via a MUSCL interpolation is

q i = q i + 1 6 φ (r i ) (q i -q i-1 ) + 1 3 φ 1 r i (q i+1 -q i ) (20)
and

q i+1 = q i+1 - 1 3 φ (r i+1 ) (q i+1 -q i ) - 1 6 φ 1 r i+1 (q i+2 -q i+1 ) (21)
where

r i = q i+1 -q i q i -q i-1 (22)
and φ(r) = max 0, min 2r,

1 + 2r 3 , 2 (23)
Nevertheless, since we use a second-order accurate Strang splitting in time, the global accuracy of the scheme remains second-order.

Time integration of the convective terms

Time integration of the convective term is made using a classical third-order TVD Runge-Kutta scheme, i.e.

Q = Q n + ∆t C (Q n ) (24) Q = 1 4 [3 Q n + Q + ∆t C (Q )] (25)
Q n+1 = 1 3 [Q n + 2 Q + 2 ∆t C (Q )] (26)
The corresponding Runge-Kutta tableau is given in Table 1.

0 1 1 1 4 1 4 1 4 0 1 6 1 6 2 3
Table 1: Butcher tableau corresponding to the compact TVD third-order Runge-Kutta method.

The convective time step ∆t is chosen to satisfy the Courant-Friedrichs-Levy (CFL) condition. In the computations, CF L = 0.5 is used.

Time integration of the stiff source terms

Due to the stiffness of the chemical source terms in case of detonation computations, a high order time integration is chosen to ensure the numerical stability. As the chemical time scale is much smaller than the convective one we apply a much smaller time step ∆t d than for the convection terms. However, this is only required in a small area of the computational domain, i.e. in the reaction zone. In the rest of the computational domain, the source term is almost equal to zero.

For this reason, we introduce an adaptive local time step ∆t d = ∆t 2N , N being the number of substeps required by the source term computation. In order to adapt this local time step with time, a fourth-fifth order embedded Runge-Kutta method is used. The computational error between the fourth-order and the fifth-order method enables to decide wether the local time step needs to be increased or decreased, as proposed by Fehlberg [START_REF] Fehlberg | New high-order Runge-Kutta formulas with step size control for systems of first order and second order differential equations[END_REF].

Q n+1 o4 = Q n + ∆t d s i=1 b i,o4 k i ( 27
)
Q n+1 o5 = Q n + ∆t d s i=1 b i,o5 k i (28) (29)
where

k i = S   t n + c i ∆t d , Q n + ∆t d s i=1 s j=1 a i,j k j   (30)
In this article, an embedded explicit Dormand-Prince method [START_REF] Dormand | A family of embedded Runge-Kutta formulae[END_REF] was chosen. The coefficients are computed in order to minimize the error of the fifth-order solution. This is the main difference with the Fehlberg method, which was constructed so that the fourth-order solution has a small error. For this reason, the Dormand-Prince method is more suitable when the higher-order solution is used to continue the time integration. The coefficients are given in Table 2. The local relative error of the quantity q is denoted by e and the accepted tolerance by ε. The optimal time step σ ∆t d can be determined by multiplying the scalar σ times the current time step ∆t d . The scalar σ is given by

σ = ε ∆t d 2 e ∆t max 1 4 = ε ∆t d e ∆t 1 4 (31)
In practice, however, the new time step ∆t d is determined in function of σ, following

∆t d =      2 ∆t d if σ > 2 ∆t d 2 if σ < 1 ∆t d otherwise (32)

Multiresolution method

The principle in the multiresolution (MR) setting is to represent a set of function cell averages as values on a coarser grid plus a series of differences at different levels of nested grids. The differences contain the information of the function when going from a coarse to a finer grid.

A tree data structure is an efficient way to store the reduced MR dataas it allows to reduce the memory with respect to a finite volume (FV) scheme on the finest level. This representation could also increase the speed-up during the time evolution because it reduces the time-searching of the elements.

In the following we consider a hierarchy of regular grids in 2D, Ω , 0 ≤ ≤ L. The root cell is Ω 0,0,0 = Ω and corresponds to a rectangle with side lengths h x and h y . The different node cells at a level > 0 forming Ω are denoted by Ω ,i,j where (i, j) ∈ Λ . The ensemble of indices of the existing node cells on the level is Λ . Note that Ω ,i,j are rectangles with side lengths h x, = 2 -h x and h y, = 2 -h y . In the tree terminology, the refinement of a parent node cell Ω ,i,j at level produces four children nodes Ω +1,2i,2j , Ω +1,2i,2j+1 , Ω +1,2i+1,2j and Ω +1,2i+1,2j+1 at level + 1, as illustrated in Figure 1. The cell-average value of the quantity u on the cell Ω ,i,j is given by ū ,i,j = 1 |Ω ,i,j | Ω ,i,j u(x, y) dx dy. and correspondingly the ensemble of the existing cell-average values at level by Ū = (ū ,i,j ) (i,j)∈Λ . The projection (or restriction) operator

P +1→l : Ū +1 → Ū .
estimates the cell-averages of a level from the ones of the level + 1. The parent cellaverage is the weighted average of the children cell-averages

ū ,i,j = 1 4 (ū +1,2i,2j + ū +1,2i,2j+1 + ū +1,2i+1,2j + ū +1,2i+1,2j+1 )
and thus the projection operator is exact and unique. To predict the cell-averages of a level + 1 from the ones of the level , we use a prediction (or interpolation) operator

P →l+1 : Ū → U +1 .
This operator yields an approximation U +1 of Ū +1 at the level + 1. In this paper, we use third order interpolation given by a tensor product approach [START_REF] Bihari | Multiresolution schemes for the numerical solution of 2-D conservation laws I[END_REF]. For n, p ∈ {0, 1}, we define

u +1,2i+n,2j+p = ū ,i,j + 1 8 (-1) n (ū ,i+1,j -ū ,i-1,j ) + 1 8 (-1) p (ū ,i,j+1 -ū ,i,j-1 ) (33)
+ 1 64 (-1) np [(ū ,i+1,j+1 -ū ,i+1,j-1 ) -(ū ,i-1,j+1 -ū ,i-1,j-1 )] .
First, this prediction is local, since it is made from the cell average ū ,i,j and the eight nearest uncles ū ,i±1,j±1 . Second, it is consistent with the projection, i.e. P +1→ • P → +1 = Id. The difference between the exact and the predicted values at three children cells yields the wavelet (or detail) coefficients. The sum of the four details in the children cells is equal to zero [START_REF] Bihari | Multiresolution schemes for the numerical solution of 2-D conservation laws I[END_REF]

d +1,2i,2j+1 = ū +1,2i,2j+1 -u +1,2i,2j+1 d +1,2i+1,2j = ū +1,2i+1,2j -u +1,2i+1,2j (34)
d +1,2i+1,2j+1 = ū +1,2i+1,2j+1 -u +1,2i+1,2j+1 .
This implies that the knowledge of the cell-average values on the children Ū +1 is equivalent to the knowledge of the cell-average values on the parents Ū and the wavelet coefficients D +1 = ( d +1,2i,2j+1 , d +1,2i+1,2j , d +1,2i+1,2j+1 ) (i,j)∈Λ . The so-called multiresolution transform on the cell-average values is obtained by repeating this operation recursively on L levels [START_REF] Harten | Multiresolution algorithms for the numerical solution of hyperbolic conservation laws[END_REF],

ŪL ←→ ( DL , DL-1 , . . . , D1 , Ū0 ).

In conclusion, knowing the cell-average values of all the leaves ŪL is equivalent to knowing the cell-average value of the root Ū0 and the details of all the other nodes of the tree structure.

In the MR scheme, instead of using the representation on the full uniform grid Ω L , the numerical solution Ū n M R = Ū n L,M R is formed by cell averages on an adaptive sparse grid Γ n = Γ n L . Grid adaptivity in the MR scheme is related with an incomplete tree structure, where cell refinement may be interrupted at intermediate scale levels. This means that Γ n is formed by leaf cells Ω ,i,j , 0 ≤ ≤ L, (i, j) ∈ L(Λ ), which are cells without children.

Here L(Λ ) denotes the ensemble of indices for the existing leaf cells of the level .

Three basic steps are undertaken to evolve the solution from

Ū n M R to Ū n+1 M R , Refinement: Ū n+ M R ← R Ū n M R Evolution: Ǔ n+1 M R ← E M R Ū n+ M R Coarsening: Ū n+1 M R ← T( ) Ǔ n+1 M R
The refinement operator R is a precautionary measure to account for possible translation or creation of finer scales in the solution between two subsequent time steps. As the regions of smoothness or irregularities of the solution may change with time, the grid Γ n may not be convenient anymore at the next time step t n+1 . Hence, before doing the time evolution, the representation of the solution should be extended onto a grid Γ n+ , which is expected to be a refinement of Γ n , and to contain Γ n+1 . Then, the time evolution operator

E M R = E M R (∆t) is applied. The subscript MR in E M R
means that only the cell-averages on the leaves of the computational grid Γ n+ are evolved in time, and that an adaptive flux computation

F M R ( Ū n+ M R
) is adopted at interfaces of cells of different scale levels. Finally, a thresholding operation T( ) (coarsening) is applied in order to unrefine those cells in Γ n+ that are unnecessary for an accurate representation of Ū n+1 M R . The choice of the threshold value is motivated by an error analysis equilibrating the perturbation and discretization errors. For details we refer to [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF][START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF].

To compress data in an adaptive tree structure, while still being able to navigate through it, gradedness is required. For instance, for a given node in the dynamic tree structure we assume that:

• its parent and eight nearest uncles are in the tree (if not, we create them as nodes);

• for flux computations, if Ω ,i,j is a leaf, its four nearest cousins Ω ,i±2,j and Ω ,i,j±2 in each direction are in the tree (if not, we create them as virtual leaves);

• if a child is created, all its brothers are also created;

For more details of these procedures, we refer to [START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF].

In the tree structure, the thresholding operator T( ) is defined by removing leaves where details are smaller than a prescribed tolerance , while preserving the gradedness property, and the refinement operation R adds one more level as security zone, in order to forecast the evolution of the solution in the tree representation at the next time step. These two operations are performed by the following procedure.

We denote by Λ the ensemble of indices of the existing tree nodes in Γ n+ , by L(Λ) the restriction of Λ onto the leaves, and by Λ the restriction of Λ to a level , 0 ≤ < L. For the whole tree, from the leaves to the root:

• Compute the details on the nodes d ,i,j , (i, j) ∈ Λ -1 , by the multiresolution transform;

• Define the deletable cells, if the details on the corresponding nodes and their brothers are smaller than the prescribed tolerance.

For the whole tree, from the leaves to the root:

• If a node and its children nodes are deletable, and the children nodes are simple leaves (without virtual children), then delete their children.

• If the node and its parents are not deletable, and it is not at the maximum level, then create the children for this node.

To illustrate the adaptive flux computation, we consider the leaf Ω +1,2i+1,2j , sharing an interface with another leaf Ω ,i+1,j at a lower scale level, as illustrated in Figure 2. For the calculation of the outgoing numerical flux on the right interface, we use the cell width in the x direction h x, +1 as step size. The required right neighboring stencils are obtained from the cousins Ω +1,2i+2,2j and Ω +1,2i+3,2j+1 , which are virtual cells. For conservation, the ingoing flux on the leaf Ω ,i+1,j is set equal to the sum of the outgoing fluxes on the neighbour leaves of level + 1. For more details on the implementation of this procedure we refer to [START_REF] Roussel | A conservative fully adaptive multiresolution algorithm for parabolic PDEs[END_REF].

Non-stiff case

In this first part, we consider a detonation with a time coefficient τ = 10 -1 . As observed in Figure 3, the detonation front propagates from left to right with constant maximal values of density, pressure and velocity, which correspond to the von Neumann state. The slopes for the density, pressure, and velocity, are moderate on the left side of the shock, i.e. in the burnt gas zone, but become sharp on the right side, at the interface with the unburnt gases. For all the displayed quantities, the curves fit with the reference computation, which has been performed on a fine grid with L = 14 scales, i.e. with 16384 grid points. Table 3 assembles results for computations performed with different maximum number of scales for both the MR and the FV method. The relative error is computed on the von Neumann value of the density. The numerical value is averaged in time to damp its oscillations and is compared to the theoretical one. We observe for the MR computations a CPU time compression rate growing with increasing number of scales, while the relative error remains comparable with the one of the same computation on the regular fine grid. As expected, the error is reduced when increasing the number of scales L.

Method

Influence of the stiffness

In this part, we perform computations of stiff problems and set τ first to 10 -2 , then to 10 -3 . In order to resolve correctly the detonation front, 13 scales are required for the case τ = 10 -3 . The results are shown in Figures 4 and5. We observe a much sharper and thinner reaction zone, reduced to a few computational points in the case τ = 10 -3 . We observe a good agreement with the reference computations, which were computed on a regular fine grid with L = 15 scales. Table 4 gives the results of the computations for τ = 10 -2 and τ = 10 -3 . As expected, a high compression rate is reached when 13 scales are used. But as it is seen in Figure 5, the error on the von Neumann density remains quite large, even if the detonation front is well tracked.

Two-dimensional detonation

Now we consider a two-dimensional problem. The interaction between a detonation wave and a pocket of unburnt gas is simulated. The computational domain is [0, 4] × [-1, 1], the initial planar detonation front is located at x = 1.7 and the initial spot of unburnt gas is centered in (x 0 , y 0 ) = (2, 0). The radius of the circular pocket is r 0 = 0.1. The parameters of the detonation front are the same as in the non-stiff case, i.e. we choose τ = 10 -1 . The dimensionless final time is t = 1.3.

In Figure 6, we observe the destabilization of the planar front by the pocket of unburnt gas. We remark circular structures due to the expansion of the pocket in the three directions and their interaction with the detonation wave going from left to right. Recirculations are observed in the center of the pocket, which are advected by the flow. We also observe reflexions of the waves on the free-slip boundaries in y = -1 and y = 1.

On the right side of Figure 6, we observe that the adaptive mesh follows well the structures of the different waves that interact in the flow. A zoom in the center of the mesh (Figure 7) enables to see how well the mesh adapts to all the structures of the flow. The computation, performed on L = 10 scales i.e. a maximum of 2 2L = 1, 048, 576 points, requires 2 h 59 min of CPU time on a 16 Core PC. The same computation on the regular fine grid lasts 14 h 35 min. Hence the CPU time compression is around 20 %.

Figure 8 shows a cut of the density at y = 0. We observe an excellent agreement between the MR and FV computations with L = 10 scales, and a good agreement of both curves with the FV reference computation obtained with L = 11 scales, which validates the grid convergence of the computation. The total number of grid points in the FV computations is 2 2L .

Conclusion and perspectives

We presented an extension of the adaptive multiresolution method for the reactive Euler equations, able to deal with fast chemical reactions. A finite volume scheme with second order shock capturing schemes is used for space discretization. Classical Strang splitting is applied to deal with the stiffness of the physical problem in time. Space and time adaptivity are then introduced using multiresolution analysis and Runge-Kutta-Fehlberg schemes, respectively. The implementation uses dynamic memory allocation with tree data Applications to detonation problems in one and two space dimensions validated the algorithm and illustrated the efficiency of the adaption strategy. The adaptive computations yield both speed-up of CPU time and memory compression with respect to uniform grid computations while the precision is automatically controlled using suitable thresholding techniques.

As perspective we plan the extension of the method to three space dimensions and also to include multi-species chemical reactions and viscous effects. The parallelization of the adaptive method using tree data structures to handle the adaptive grid remains a challenging task for the near future.

Table 2 :

 2 Butcher tableau corresponding to the Dormand-Prince method. The first row of b coefficients gives the fourth-order accurate solution, and the second row yields order five.

Figure 1 :

 1 Figure 1: Dyadic grid refinement in 2D.

Figure 2 :

 2 Figure 2: Adaptive numerical flux computation in 2D.

Figure 3 :

 3 Figure 3: Chapman-Jouguet detonation front: density (top, left), pressure (top, right), temperature (middle, left), partial mass of the limiting reactant (middle, right), velocity (bottom, left) and adaptive mesh (bottom, right) at t = 2.2, L = 11 scales, = 2.5 • 10 -3 , τ = 10 -1 .

Figure 4 :

 4 Figure 4: Chapman-Jouguet detonation front: density (top, left), pressure (top, right), temperature (middle, left), partial mass of the limiting reactant (middle, right), velocity (bottom, left) and adaptive mesh (bottom, right) at t = 2.2, L = 12 scales, = 1.25 • 10 -3 , τ = 10 -2 .

Figure 5 :

 5 Figure 5: Chapman-Jouguet detonation front: density (top, left), pressure (top, right), temperature (middle, left), partial mass of the limiting reactant (middle, right), velocity (bottom, left) and adaptive mesh (bottom, right) at t = 2.2, L = 13 scales, = 6.25 • 10 -4 , τ = 10 -3 .

Figure 6 :

 6 Figure 6: Time evolution of the interaction between a detonation wave an a pocket of unburnt gas computed with the adaptive MR method. Numerical Schlieren of the density gradient (left) and corresponding meshes (right) at t = 0 (top), t = 0.65 (middle) and t = 1.3 (bottom).

Figure 7 :

 7 Figure 7: Adaptive MR computation of the interaction between a detonation wave an a pocket of unburnt gas: Zoom on the adaptive mesh at t = 1.3.

Figure 8 :

 8 Figure 8: Interaction between a detonation wave an a pocket of unburnt gas. Onedimensional cuts of density for y = 0 at t = 1.3. Comparison of the adaptive MR and the FV computations for L = 10 with the FV reference computation for L = 11. The total number of grid points in the FV computations is 2 2L .

Table 3 :

 3 Chapman-Jouguet detonation front: CPU time, CPU compression and relative error on the von Neumann density for different scales and for the FV and MR methods, τ = 10 -1 .

  -3 .

	Method	τ L	CPU time CPU comp.	rel. error
	MR	10 -2 12 1.25 • 10 -3	43 s	20,08 % 1.336 • 10 -1
	FV	10 -2 12	3 min 37 s	100,00 % 1.342 • 10 -1
	MR	10 -3 13 6.25 • 10 -4	2 min 03 s	7,56 % 3.057 • 10 -1
	FV	10 -3 13	27 min 13 s	100,00 % 3.057 • 10 -1

Table 4 :

 4 Chapman-Jouguet detonation front: CPU time, CPU compression and relative error on the von Neumann density for different values of the stiffness coefficient and for the FV and MR methods.

Acknowledgements

This article is dedicated to Professor Henning Bockhorn on the occasion of his 70th birthday, thanking him cordially for the great time we had in his group, first in Kaiserslautern and then in Karlsruhe. We also thank Margarete Domingues, Ralf Deiterding and Sonia Gomes for constructive discussions on the topic and fruitful interactions. KS and OR thankfully acknowledge financial support from the ANR projects SiCoMHD (ANR-Blanc 2011-045) and MAPIE (ANR-13-MONU-0002).

Numerical results

In all numerical simulations, we consider an initially planar Chapman-Jouguet detonation moving with constant unitary speed through the unburnt gas to the right of the domain. Setting γ = 1.4, r = 1, Q 0 = 1 and T i = 0.22 and assigning the following burnt gas values

where

the corresponding von Neumann state past the shock wave is [START_REF] Tosatto | Numerical solution of under-resolved detonations[END_REF][START_REF] Helzel | A modified fracitonal step method for the accurate approximation of detonation waves[END_REF]

and the unburnt gas values are

The resulting temperature of the unburnt gas T u = 0.215995 is only slightly lower than the chosen ignition temperature T i . The initial condition for an initial front at x = x 0 is

where a is set to 1 τ . Here q stands for ρ, v x , p, and Z, while v y = 0 everywhere. This way the initial condition is close to the classical Zeldovich-von Neumann-Dring (ZND) solution of the reaction Euler equations.

One-dimensional detonation

First we consider a one-dimensional setting. The computational domain is Ω = [-3, 1]. In this case, the location of the initial front is set to x 0 = -2. In order to follow the detonation wave and perform longer computations, a change of variables is made for v x . We set v x = v x -δ and omit the superscript everywhere. Since the initial transition happens in the segment [-3, -1], we only focus on the domain [-1, 1]. The dimensionless final time is t = 2.2.


Léonid Arantchouk

G Point

Yohann Brelet

Jean Larour

Jérôme Carbonnel

Yves-Bernard André

A Mysyrowicz

Aurélien Houard

Y.-B André

Compact 180-kV Marx generator triggered in atmospheric air by femtosecond laser filaments

come   L'archive ouverte pluridisciplinaire

Proposed in 1924 by Erwin Otto Marx, the Marx generator is now commonly used in high-energy physics (electron and ion beam acceleration, pinch plasma, microwave generation, etc) and for pulse power applications (industry, medical, environment, agriculture, etc). Marx generator (MG) consists in a set of parallel stages including capacitor banks and switches. The capacitors are first charged in parallel to a voltage V 0 , and are then connected in series through closing switches. In the traditional operation mode, commutation of the lower stage induces an overvoltage and self breakdown on the following stages. In this way, all individual switches close, producing an output voltage theoretically given by the product of the charging voltage V 0 by the number of stages N: V MG = NV 0 . In practice, due to the intrinsic and stray capacitance of the switches, the voltage multiplication across the switch gap is lower than the theoretical value. This results in rather complex construction solutions and reduction of the generator efficiency.

The switch performance is related to its design and triggering conditions. The reader is directed to ref. [1][2][3][START_REF] Calico | Pulsed Power Conf. IEEE, Digest of Techn. Papers[END_REF][START_REF] Dragt | Pulsed Power Plasm. Science, PPPS-2001[END_REF][START_REF] Goerz | Digest Technical Papers[END_REF][START_REF] Neuber | Pulsed Power Conf. IEEE[END_REF] to find some practical realizations of Marx generators with high-level characteristics. Several techniques were developed to improve the spark gap (SG) performances: placing the gap in an appropriate gas mixture at high pressure (up to 20 bars), inserting Marx modules in metallic housing to increase the stray capacitance or using a line-of-sight geometry to obtain a stable Marx erection due to ultra-violet photo-ionization.

In several publications a jitter of 1 ns or less has been achieved by means of rather complicated methods. In ref. [START_REF] Kusama | [END_REF] the Marx stages are operated very close to self-break regime (95%) in dry air at 4 bars. A KrF excimer laser beam passes through a metal mesh installed in the switch electrode holes. Laser interaction with the mesh produces photoelectrons, which in turn trigger a discharge in the two adjacent switches simultaneously. In ref. 9 the spark gaps breakdown was considerably accelerated by connecting the individual stages with coupling capacitors, while a needle auxiliary corona discharge was used to reduce the breakdown temporal jitter. A subns jitter was reported in [START_REF] Kellog | Pulsed Power Conf. IEEE[END_REF] with a mini-Marx using a UV laser with a gas mixture containing SF 6 at 3 bars in the gap switches. A similar jitter was also reported in ref. [START_REF] Mayes | Pulsed Power Plasm. Science, PPPS-2001[END_REF] .

Triggering the gap switch using femtosecond laser filament is based on a different concept. Laser filamentation is a nonlinear phenomenon appearing during propagation of intense femtosecond laser pulses provided their peak power exceeds a critical threshold P cr of the order of 5 GW in air (at 800 nm). A dynamic competition between nonlinear self-focusing and multiphoton ionization leads to a contraction of the laser beam into thin light channels of 100 µm [START_REF] Couairon | [END_REF] . This filamentary laser pulse leaves in its wake a long plasma string with an electron density of 10 16 cm -3 that has the ability to trigger and guide electric discharges over several meters [13][14][15][16] . In ref. 17 we presented a gap switch triggered with femtosecond laser filaments. A sub-nanosecond jitter was obtained in atmospheric air even with a 1-cm gap when the switch was charged at 70% of the self breakdown limit. According to ref. [START_REF] Mesyats | Pulsed Power[END_REF] , the characteristic time for discharge establishment in a 10 mm spark gap in atmospheric air is about T SG ~ 4 ns. If several SGs are triggered with the same laser beam and if their jitter << T SG , one should obtain then an ideal Marx generator with simultaneous switching of all stages. This result has been the motivation to develop a MG triggered by a single filamentary laser beam.

In this experiment the ENSTAmobile laser with a maximum energy per pulse of 300 mJ was used to create the plasma filaments in air 17,[START_REF] Brelet | [END_REF] . The laser pulse duration was fixed to 700 fs corresponding to a peak power about 85 P cr . In this regime, the laser beam focused by a 5 m focal lens produced a tight bundle of ~50 plasma filaments extending over 1.5 m before and 0.5 m after the geometrical focus of the beam. The diameter of this plasma column was measured to vary between 2 and 5 mm. A two meter long plasma channel allows triggering practically simultaneously all the gap switches installed in line-of-sight geometry on the laser propagation axis.

A single stage of our MG (Fig. 1) is composed of two parts. The first one presents six 2 nF "door knob" ceramic capacitors in parallel connected by two metal discs. The second part is a spark gap made by two 35-mm diameter cylindrical electrodes with 5 mm radius curved border. The laser beam enters the spark gap through a 3-4.5 mm diameter hole drilled axially in the electrodes. The inter-electrode gap D is fixed by a Plexiglas ring which serves at the same time as a spacer between the Marx's stages. We worked with inter-electrode spacing D =11.25 mm and 12 mm corresponding respectively to self-breakdown voltages V cr of 34 kV and 36 kV. A DC power supply charged the stage capacitors to a voltage V 0 = 30 kV, corresponding to a self break spacing D = 9.7 mm. In practice, to avoid breakdown hazard and reduce the pre-fire safety factor it is necessary to work as far as possible from the self-breakdown limit. The assembled stages were successively slided over six nylon allthread rods and were tightened up with nuts. design allows arranging the gap switch holes in an optical lineof-sight with a good precision so that the laser beam passes through all the stages.

Our Marx generator is composed of up to 6 stages, with a diameter of 250 mm and a total length of 55 cm. It generates up to 180 kV output voltage. The authors of ref. 20 have proposed to classify MG into three categories according to their overall parameters. In accordance with such parameters our MG belongs to the ultra compact designs, but higher output voltage let us relate the generator to the compact category.

To characterize the Marx generator we studied two parameters: the runtime, which is the delay between the laser triggering and the breakdown of the generator, and the jitter, defined as the standard deviation of the runtime. The statistics presented here were based on 10-15 shots. For high voltage measurements we used a voltage divider made up from HVR resistors (HVR International). The discharge current was monitored with a current viewing resistor (CVR) in series with the load, and the laser firing time was monitored with a fast photodiode.

The output voltage signal of the 6-stages Marx with V 0 = 27 kV presents an exponential decay τ = RC ≈ 5 µs when the generator discharges on a high resistance load R 1oad = 2.7 kΩ (Fig. 2(a)). The front edge of the voltage pulse is shown in Fig. 2(b). To avoid influence of electromagnetic noises at the discharge beginning, we define the time T Marx corresponding to the 50-percent point of the voltage maximum as a reference for the Marx discharge beginning. The point T 0 = 0 corresponds to the laser beam arrival. Runtime delay corresponds then to T Marx -T 0 . We measured a voltage amplitude corresponding with 5% precision to the value NV 0 .

Runtime and jitter for the 6-stages Marx are presented in Fig. 3 as a function of the normalized breakdown voltage V 0 /V cr . The runtime decreases from 60 ns to 20 ns with increase of V 0 /V cr from 60% to 80%. The jitter is less or about 1 ns for V 0 /V cr ≥ 65%. We also measured the Marx performances when the laser was shot from the output side of the Marx. A runtime about 3 ns longer was observed in this case, while the jitter remained unchanged. Several configuration of spark gap holes were tested, with a constant diameter for all stages or with decreasing diameter ranging from 4.5 to 3 mm. Finally the best characteristics were obtained with the hole configuration Ø 3 mm, as shown in Fig. 3.

To confirm the effect of the filament in the triggering process we performed a test while blocking the laser beam at the third stage of the MG. We observed in this case an increase of the runtime from 60 to 350 ns with Marx charged up to 60% self-breakdown voltage.

Finally we determined the minimum laser energy required to maintain a subnanosecond jitter in the generator. With a charging voltage V 0 = 27 kV we were able to keep the same runtime of 20-25 ns with a ns jitter by decreasing the laser energy down to 30 mJ, for a pulse duration of 700 fs. This corresponds to a peak power of 43 GW, which is close to single filamentation regime. Similar characteristics were also obtained with laser energy of 6 mJ with a pulse duration of 100 fs (60 GW).

The laser induced filament beam was used not only to trigger the Marx, but also to produce a guided discharge (GD) at the output of the generator. For this purpose a rather important output load resistance was chosen, to ensure the voltage keeps a high value for a reasonably long time. The advantage of this setup is that the laser is automatically synchronized with the high voltage pulse. To produce a long guided discharge a second electrode E 2 was installed at a distance L d from the Marx output electrode E 1 as shown in Fig. 1(b). The laser beam was aligned so that the filaments touched tangentially both electrodes. The second electrode was loaded on a resistance R d ~ 300 Ω.

Typical signals of the discharge current and the output MG voltage are depicted in Fig. 4. Up to the moment of the long discharge beginning, the MG discharges through the high load R load = 2.7 kΩ. After breakdown occurs between E 1 and E 2 , a current peak appears over the quasi-pedestal created by the previous high resistance discharge. Since the resistance of the long discharge is relatively small, this current peak amplitude is defined by the resistance R d = 300 Ω and by the Marx voltage when the discharge occurs. In our experiments the long discharge maximal current was measured to be ~450 A.

We studied the time delay between Marx voltage pulse and appearance of the guided discharge for two configurations. In the first one, the laser passed through the Marx first and then through the region of the guided discharge (MGD configuration). In the second geometry the laser passed between the two long discharge electrodes first and then entered into the MG (GDM configuration). As it was mentioned above, the Marx characteristics are weakly sensitive to the laser direction.

The measured delays between the Marx output voltage and the guided discharge are shown in Fig. 5 as a function of the reduced discharge field E disch for both configurations. Measurements were performed with discharge lengths ranging from 12 cm to 21 cm, showing that the delay mostly depends on E disch rather than on the discharge length L d . The best triggering characteristics were achieved when the laser enters the generator from its output side (GDM configuration) as shown in Fig. 5(b). In this configuration delay and jitter are almost two times smaller than for the other configuration. As a reference, a delay Marx-Discharge smaller than 300 ns with a jitter less than 60 ns is measured when the discharge field is higher than 10 kV/cm. This higher guiding efficiency in GDM geometry can be explained by the higher laser energy participating in the preparation of the guided discharge, since about 50% of the laser energy is blocked by the sparks gaps holes of the Marx.

Finally one should mention that the discharge triggering characteristics are sensitive to output Marx voltage polarity. The maximal distance L d = 21 cm was only observed with a positive polarity.

In conclusion, we have presented a design of Marx generator triggered by means of laser filaments produced in air by a TW femtosecond laser source. High voltage pulses of 180 kV with a runtime ≤ 20 ns and a temporal jitter ≤ 1 ns could be produced with a charging voltage higher than 70% of self-breakdown. We determined that the Marx performance in terms of jitter and delay start to degrade when the laser input is lower than 50 GW. The same bundle of laser filaments was also used to produce guided discharges up to 21 cm long. The delay between the discharge and the Marx pulse is mostly dependent on the ratio E disch = V Marx /L d . The best discharge characteristics were observed using a geometry where the laser beam enters from the output generator side and with a positive voltage polarity.

FIGURESFIG. 1 .

 1 FIGURES

a) Author to whom correspondence should be addressed. Electronic mail: leonid.arantchouk@lpp.polytechnique.fr.

ACKNOWLEDGMENTS

This work was supported by the French Direction Générale de l'Armement (grant EPAT n°2007 95091 and grant DGA-Ecole Polytechnique 2009). The authors would like to thank J.-L. Charles for his technical help.


Guillaume Point

Yi Liu
email: yi.liu@ensta-paristech.fr

Yohann Brelet

Sergey Mitryukovskiy

Pengji Ding

Aurélien Houard

André Mysyrowicz
email: andre.mysyrowicz@ensta-paristech.fr

Lasing of ambient air with microjoule pulse energy pumped by a multi terawatt IR femtosecond laser

We report on the lasing action of atmospheric air pumped by an 800 nm femtosecond laser pulse with peak power up to 4 TW. Lasing emission at 428 nm increases rapidly over a small range of pump laser power, followed by saturation above ~ 1.5 TW. The maximum lasing pulse energy is measured to be 2.6 J corresponding to an emission power in the MW range, while a maximum conversion efficiency of

is measured at moderate pump pulse energy. The optical gain inside the filament plasma is estimated to be excess of 0.7/cm. The lasing emission shows a doughnut profile, reflecting the spatial distribution of the pump-generated white-light continuum that acts as a seed for the lasing. We attribute the pronounced saturation to the defocusing of the seed in the plasma amplifying region and to the saturation of the seed intensity.

'

OCIS Codes: (320.7120) Ultrafast phenomena (140.4130) Molecular gas lasers Stimulated emission from air plasma columns has attracted wide attention because of its potential application for remote sensing [1][2][3][4][5][6][7][8][9][10][11][12]. Backward Amplified Spontaneous Emissions (ASE) originated from both nitrogen and oxygen molecules have been obtained under different laser excitation conditions [1][2][3]. In 2011, J. Yao and co-workers reported several strong stimulated emission lines in the forward direction in ambient air or pure nitrogen [7]. These radiation lines at 391, 428, and 471 nm can be activated either by a wavelength tunable midinfrared laser (1.29 m -4 m) or an 800 nm femtosecond pulse together with an additional seed pulse at a proper wavelength [7][8][9]. Forward single pass lasing of these emission lines has been attributed to the following mechanism: rapid population inversion between the   u B 2 and   g X 2 states of the  2 N induced by the pump laser, followed by seeded amplification. The seed is either the 3 rd and/or 5 th harmonic signal of the mid-infrared pump or an external seed pulse around 400 nm in the case of the 800 nm pump.

Very recently, three groups have independently reported self-seeded forward lasing action from  2 N in air or pure nitrogen pumped by just one single pulse at 800 nm [10][11][12]. Seeding in this case is provided by the generation of 2 nd harmonic and supercontinuum white-light inside the gas [10][11][12]. Discovery of the self-seeded air lasing effect pumped by the widely available 800 nm femtosecond laser simplifies the experimental setup and offers the possibility of scaling up to Terawatt (TW) class pump lasers. In previous experiments so far, however, the IR pump power was limited to a few tens of GW [7][8][9][10][11][12]. Consequently, the energy of the lasing emission was too weak to be measured.

In this letter, we study self-seeded lasing from laser filaments in ambient air driven by 800 nm femtosecond pulses with a peak power of 4 TW. In the experiment, we used the Enstamobile laser system which delivers 50 fs laser pulses at a repetition rate of 10 Hz with pulse energy up to 300 mJ. The beam exhibits a super-Gaussian profile with a diameter of 40 mm (FMHM). The output pulses were focused by a convex lens of focal length f = 100, 150 and 200 cm. Best results were obtained with the f = 100 cm lens. After the plasma, two pieces of glass filters (BG 40) were installed in the beam to remove the strong infrared driving pulse and the accompanying white light with wavelength longer than 650 nm. For detection, a convex lens of f = 10 cm was used to collect the radiation into a fibre spectrometer (HR 4000, Ocean Optics). The energy of the lasing emission was measured with a sensitive laser pulse energy meter (model: OPHIR, NOVA, PE9-C). We also measured the spatial profile of the lasing radiation with a charge-coupled device (CCD) camera.

In Fig. 1, we present the measured spectrum of the forward UV radiation as a function of incident laser pulse energy. At ambient pressure, the only line showing significant gain is at 428 nm. It appears as a narrow-band emission peak superimposed upon the tail of a broad continuum extending to larger wavelength for pump pulse powers larger than 150 GW (7.5 mJ). Upon increase of the incident laser power, the radiation at 428 nm increases progressively and shows saturation beyond P = 1.5 TW, as shown in Fig. 2. The measured width (FWHM) of this radiation was 1.2 nm, limited by the resolution of the spectrometer. Note that the real spectral width is on the order of ~ 0.1 nm [13]. This radiation line has been identified as due to the transition between the second excited state(

  0 2     u B ) of the  2 N and its ground state (   1 ' 2     g X
), where and'denote the molecular vibrational quantum number of the corresponding states [10][11][12]. In all our experiment, lasing radiation at 391 nm was rarely observed and was very weak compared to that around 428 nm. . This gain value is probably seriously underestimated here due to the limited spectral resolution of the spectrometer, which results in an order of magnitude decrease of the observed peak emission intensity at 428 nm. The gain per unit length for a high pressure gas laser reads (1) where  

N N  ' is the population inversion per unit length, v and  Z v  
are the frequency of the laser emission and the line width with Z the molecule collision frequency,  is the radiative lifetime [14].

Taking s Z 9 10 52 . 6  
at standard temperature and pressure and  ≈ 100 ns [15], we calculated   N N  ' to be 1.8 × 10 16 /cm. The population inversion density per volume can then be estimated to be 4 × 10 17 /cm 3 , by considering a transverse plasma size of ~ 1.2 mm which will be discussed below. We measured the energy of the lasing radiation at 428nm with the laser energy meter equipped with an interference filter at 428 nm (10 nm bandwidth). The net energy of the lasing radiation was calibrated with careful consideration of the spectrum weight of the lasing emission in the total transmitted radiation, the reflection loss of the exit fused silica windows and the collecting lens, the transmission property of the BG 40 color filter and of the interference filter. A calibrated energy of 2.6 J was deduced for the 428 nm lasing emission with Pin = 2.4 TW. The maximum conversion efficiency was found to be We measured the spatial profile of the 428 nm lasing radiation with a CCD camera equipped with the same in terference filter. The results ar e presented in Fig. 3(a)-3(d). The radiation shows a doughnut-shaped pattern, with a slight increase of its opening angle for higher incident pulse energies. In Fig. 3(a), the divergence angle of the emission was estimated to be 18.6 mrad. Such a small angle confirmed again the stimulated radiation nature of this emission. In a previous work, we tentatively attributed this doughnut-shaped radiation to the conical shape of the seeding white-light pulses [10].

Here we examine this hypothesis in more details. First, we monitored the spatial distribution of the broadband continuum at a wavelength close to the seed at 428 nm. The profile of the white-light component around 471 nm (10 nm bandwidth), shown in Fig. 3(e), displays a similar doughnut shape. Next, we checked the profile of the entire white-light emission on a screen installed 3 meters away from the geometric focus of the f = 1 m focusing lens. The result presented in Fig. 3(f), shows a darker central area surrounded by visible radiation with different colors. Based on the above two observations, we conclude that the spatial distribution of the white-light seeding pulse is at the origin of the doughnut shaped 428 nm single pass lasing. Such conical emission has been widely observed during filamentation in air and solid and several mechanisms have been proposed for its explanation [16]. Delete Four key factors for the single pass lasing from filaments can be identified, namely the incident seed pulse energy, the density of population inversion, the volume of the plasma amplifier, and the spatial overlap between the seeding pulse and the amplifying plasma media. The role of the first three factors are expected to be similar to that in a traditional laser amplifier [17], while the last factor is particular for this filaments amplifier due to the non uniform transverse profile of the laser intensity , the plasma nature of the amplifier, and the formation of multiple filaments [18]. In the following, we examine the evolution of these factors as a function of the incident IR pulse energy. We first checked the volume of the filaments plasma amplifier. The images of the filaments taken by a CCD camera are presented in Fig. 4 as a function of incident laser pulse energy. The length and width of the filaments are seen to increase at higher pulse energy. When the pulse energy exceeds 84 mJ, clear evidence for multiple filaments formation at the initial stage of filamentation is observed (Fig. 4 (e)-4(g)). We numerically integrated each rows (lines) of the fluorescence plasma track and defined the plasma channel length (width) by 1/e 2 of the maximum value.

        N N c g
The results are presented in Fig. 5 (a). The plasma volume is estimated as

  L d V   2 2



, where d and L are the length and width of the plasma respectively. We show the results in Fig. 5(b), where a monotonous increase is observed at higher pump pulse energies. We next discuss the density of population inversion. The mechanism for population inversion is still a matter of controversy [7][8][9][10][11][12]. As a result, a direct evaluation of population inversion inside the plasma channel is impossible at this stage. Nevertheless, it is reasonable to assume that the population inversion is related to the pump laser intensity, based on the observation that the population inversion is achieved within the driven pulse duration [8,9]. It is known that the laser intensity inside a single filament remains almost unchanged, or "clamped" for a given focusing geometry above a critical power threshold, due to the interplay between the Kerr self-focusing effect and the plasma defocusing effect [19]. In the case of many filaments compressed to a bundle, the peak intensity can increase above the clamping value, as we will report elsewhere. Therefore, we expect that the density of population inversion is either roughly constant or even increases in our experiment. Since the amplifying medium increases while its population inversion density increases or stays constant, the saturation must have another origin. We now discuss the dependence of the seed pulse energy with pump energy. Lasing action occurs at the spectrum position of the seed pulse, which prohibits a net measurement of the seed pulse energy as a function of pump laser energy in the absence of gain. Therefore, we analyzed the evolution of the seed pulse en ergy at anoth er c losely spac ed wavelength, at 471 nm with a bandwidth of 10 nm.

The results are presented in Fig. 2. Upon increase of the pump pulse energy, the white light spectrum around 471 nm increases and progressively saturates, but the saturation is less pronounced than for the lasing signal. Therefore, the strong saturation of the 428 nm lasing radiation for Pin > 1.5 TW is not due solely to the saturation of the seeding pulse. Finally, we discuss the degree of spatial overlap between the seeding pulse and active amplifying region. Since the amplifying region is a plasma, it acts as a defocusing medium. Under the joint action of self-focusing and plasma defocusing, the fundamental laser pulse and the accompanying white-light supercontinuum can form complex transverse distribution inside the filaments [18], which deteriorates the spatial overlap between the seeding pulse and amplifying plasma. As presented in Fig. 3(e) and 3(f), the non uniform character of the seeding pulse transverse profile is apparent in the far field distribution. In addition, the formation of multiple filaments can exacerbate the transverse non-uniformity of the plasma amplifier. We therefore speculate that the plasma defocusing effect also plays an important role for the saturation of the lasing yield when Pin becomes larger than 1.5 TW.

To conclude, the above experimental observations show that upon increase of pump pulse energy, the density of population inversion keeps unchanged or increases, while the plasma amplifier volume increases considerably. Consequently, the gain saturation cannot be explained by complete depletion of the population inversion. To the contrary, our analysis indicates that there is room for an increase in lasing efficiency in air with TW lasers. On the other hand, the white-light seed pulse shows pronounced saturation and most probably experiences defocusing during its propagation in the amplifying plasma medium. We believe that these last two factors are at the origin of the lasing yield saturation.

In summary, we observed lasing action of atmospheric air pumped by 800 nm laser pulses with peak power up to 4 TW. The radiation yield first increases rapidly but saturates for pulse power larger than 1.5 TW. The maximum lasing radiation energy was 2.6 J and the maximum conversion efficiency was . We attribute the lasing energy saturation to the sub-linear increase of the seed pulse with increasing pump power and the effect of plasma defocusing on the seed propagation.

Fig. 1 .

 1 Fig. 1. (a), Spectrum of the forward emission after filamentation of laser pulses in ambient air. The focal length was 100 cm. (b), Zoomed spectrum around 428nm in the case of Ein = 65 mJ. The dotted line shows the estimated contribution of the supercontinuum white-light. We estimated the optical gain inside the filaments plasma as   L I I g seed out ln 

Fig. 2 .

 2 Fig. 2. Dependence of the 428 nm lasing signal and the 471nm white-light spectrum component energy as a function of the incident laser pulse energy.

Fig. 3 .

 3 Fig. 3. (a)-(d), Spatial profiles of the 428 nm lasing radiation for different incident laser energies. The incident laser pulse energy was 18, 30, 103, and 205 mJ, respectively. (e), Profile of the white-light component around 471 nm. (f), Image of the total white-light conical emission on a screen situated 3 m away from the focal length. For (e) and (f), the incident laser energy was 80 mJ.

Fig. 4 .

 4 Fig. 4. Images of the filaments with different incident pulse energies, as indicated in each panel. The arrow in (a) denotes the direction of laser propagation.

Fig. 5 .

 5 Fig. 5. (a), Length and width of the plasma channel as a function of the incident laser energy. (b), Calculated filaments volume.


Sergey I Mitryukovskiy

Yi Liu

Bernard Prade

Aurélien Houard

André Mysyrowicz
email: andre.mysyrowicz@ensta-paristech.fr

Effect of an external electric field on the coherent terahertz emission from multiple filaments in air

We describe the terahertz emission from multiple filaments in air in the presence of an external electric field. A strong enhancement of the radiated terahertz energy is obtained by the combined effect of terahertz field interference and presence of a static electric field.

Introduction

Terahertz (THz) radiation from plasma in air induced by laser pulses was first demonstrated by Hamster and co-workers about 20 years ago [1]. Since then, the THz radiation from a long plasma channel, a filament, generated by an intense femtosecond laser has attracted much attention [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. For a filament formed by a single femtosecond laser pulse, the THz radiation is produced by the longitudinal plasma oscillations in the wake of the moving ionization front. These plasma oscillations are excited by the ponderomotive force of the laser pulse and are damped by electron collisions on a picosecond scale. The radiation is emitted by a dipole like structure moving at the speed of light; it is therefore a Cherenkov-like emission [3]. A laser filament as a THz source has several unique advantages comparing to the traditional techniques for THz generation including photoconductive antenna and optical rectification with electro-optical crystal. Being a femtosecond filament a THz source can be positioned far away from the laser system, thereby avoiding the absorption of the THz field by water vapor during propagation in ambient air. At the same time, high energy laser pulses can be employed for THz generation without damage to the emitter, which is impossible with the traditional techniques. Unfortunately, the conversion efficiency of this method is very low. It has been measured to be on the order of 10 -9 [4]. However, several methods have been demonstrated to enhance the THz field. An enhancement of a three orders of magnitude has been reported with application of a static or optical electric field [4][5][6][7][8][9][10][11].

In our previous letter [12], we have shown that it is possible to add coherently the THz emission from multiple filaments, leading to an increase of the radiated THz energy scaling like N², where N is the number of filaments. In this letter, we consider experimentally and theoretically the effect of a static electric field on the THz emission from multiple filaments.

Theoretical model and calculations

We first consider the effect theoretically. Because the longitudinal electric field merely increases the THz intensity without changing the radiation pattern [8] it is not necessary to consider it. We therefore focus our attention on the case of transverse electric field, where the action of the applied field leads to new emission patterns from filaments.

Our calculations are based on the Cherenkov model presented in [7]. The filament propagates along the z axis (see Fig. 1(a) and Fig. 2). The THz field distribution in the far field from a single filament in the presence of a transverse electric field in the x direction is described by:

                     . cos 1 1 sin ) ( 4 1 , cos 1 1 cos cos ) ( 4 1 , 0 ) cos 1 ( 0 0 ) cos 1 ( 0 0                  ikL ikR ikL ikR R e R e j E e R e j E E (1)
here  L, k represent the frequency of the considered THz component, the length of the filament, and the THz wave number, respectively. R is the distance from the filament to the detector, θ is the polar angle and φ is the azimuthal angle (spherical coordinates) (Fig. 2). The spectral current amplitude at frequency  is given by

         e pe pe S i E j    2 2 2 0 , (2)
where is the electron plasma frequency, ν e is the electron collision frequency.

In the presence of an applied transverse field, the THz radiation pattern gradually evolves from a conical pattern to a single lobe forward oriented emission for E > 3 kV/cm [7]. In what follows, we always consider applied fields with E > 3 kV/cm. The 3 dimension (3D) emission diagram for the single filament is presented in Fig. 1(a).

The Poynting vector of the THz emission scales like

      2 2 2 2 2 2 0 0 2 cos 1 2 sin ) cos 1 ( sin cos cos 32 1 2 1               j c L E H H E S                , (3)
therefore the intensity of THz emission depends on external electric field like 2 S E . In [7] the quadratic dependence of the energy of the THz radiation on external electric field has been verified experimentally.

One can extend the model to the case of N identical filaments. The interference of the THz fields will give:

2 1 2 1 ) , , ( ) , , (       N i i i i i N i i i i i total i i e R E e R E S         , ( 4
)
Here i i   
is the accumulated phase difference of the THz fields due to the temporal delay  i between the laser pulses and R i is the distance from the i th filament emitter to the detector.

The normalized calculated 3D distribution of the THz field from two parallel filaments is presented in Fig. 1(b). An enhancement of the THz intensity by a factor 4 is achieved comparing with the case of the single filament in the presence of the same applied field.

Experimental results

We now compare our calculations with experimental measurements. We study experimentally how the THz emission depends on the variations of temporal delay  p and spatial separation d between two parallel filaments in the presence of a transverse field. In our experiments, a femtosecond laser system delivering 40 fs pulses with energy up to 15 mJ at 800 nm was employed. The experimental set-up is presented schematically in Fig. 2. We used a Mach-Zehnder interferometry scheme to obtain two parallel propagating pulses. The spatial separation and temporal delay of the two pulses could be controlled by a steering mirror of the interferometer and an integrated mechanic delay line. The two pulses with equal energy of 1.1 mJ were focused with two convex focusing lenses of f = 100 cm to form two 10 mm long parallel, laterally separated filaments in ambient air. The forward THz radiation was detected by a heterodyne detector working at 0.1 THz [13]. The distance between the filaments and the detector was around 30 cm. In our experiment the detector measures only the component in the plane of the two filaments (XOZ plane). The static electric field was applied to the ionized region transversally by placing two electrodes around the filaments. The electrodes were square copper plates with dimensions of 5 cm × 5 cm. We have studied two cases: when the direction of the electric field is parallel (configuration A) or perpendicular (configuration B) to the plane defined by the filaments. The applied field was 4 kV/cm. We first studied the interference of two synthesized ( p = 0) THz fields for different spatial separations in configuration A. The results are presented in Fig. 3. Upon increase of the spatial separation d, the THz radiation pattern becomes narrower, but remains symmetric about the propagation direction of the laser pulses. Note that there is no decrease of the THz intensity within the range studied here. This differs from the results obtained in the absence of an electric field [12], where the THz intensity starts to decrease with the increase of the spatial separation between filaments. We attribute it to the different size of the intersection areas of two THz radiation patterns from each individual filament with and without external electric field. Without external electric field the THz emission from a filament has a narrow conical distribution, thus the overlap of the THz radiation patterns from two filaments is a sensitive function of separation between them. With the transverse electric field, the broader forward lobe allows larger filament separations while interference of the THz emissions is still present.

Next we studied the interference of two synthesized THz fields for different temporal delays with a fixed spatial separation d = 2 mm between two filaments, always in configuration A.

The results are presented in Fig. 4. Upon increase of the delay  p , the total radiation pattern becomes asymmetric, alternatively directed to either side of the filament propagation direction. The change in direction reflects the partial destructive interference between the THz fields due to the phase shift induced by the delay. At the same time, the total THz radiation decreases. Finally ( p = 5.5 ps), the total THz radiation pattern comes to near complete cancelation of the THz field, keeping only two small lobes in the non-intersection area.

Eventually, for  p = 11 ps, the THz radiation recovers the value without time delay. Note that the strongest enhancement of THz radiation is always achieved with no temporal delay. By changing the temporal delay one can redirect the THz radiation along another direction that is not collinear with the laser beams. However, in most cases, the peak intensity decreases.

We have also studied the case of configuration B, when the direction of the transverse external electric field is perpendicular to the plane of filaments and found results similar to the ones presented in Fig. 3 and Fig. 4.

The comparison between the experiments and calculations is excellent, as can be seen in Fig. 3. We have considered theoretically the THz emission at a frequency of 0.1 THz, because it corresponds to the spectral component measured by our detector. However the THz emission from filaments peaks at higher frequencies. In Fig. 5 we show the temporal waveform of a THz pulse and the corresponding spectrum obtained from a 5 mm long filament in air in the presence of an external field of 8 kV/cm. The THz waveform was obtained by the standard time domain THz spectroscopy system based on a 1 mm thick ZnTe crystal [14]. The THz power spectrum reaches a maximum at 0.5 THz. It is therefore instructive to calculate the emission pattern at other frequencies. In Fig. 6(b) the result for 0.5 THz is presented for two parallel filaments separated by 500 µm. Such spatial separation was chosen to keep the ratio d/λ on the same order as for 0.1 THz (d is the spatial separation between filaments, λ is the wavelength of the THz emission). This parameter is important for the constructive interference because the THz radiation pattern becomes narrower with higher THz frequency and thus the intersection area of two THz emission from two filaments decreasing if the separation between the filaments is fixed. For  p = 0 ps, the peak THz intensity is 4 times larger than for the single filament with the same external field (see Fig. 6(a)) showing again complete constructive interference when there is no time delay between the laser pulses. The behavior with time delay between the laser pulses is similar to the case presented in the Fig. 4.

It is periodic, but the period is different. The period corresponds to 1/ d , where  d is the considered THz frequency.

Extrapolation to the case of N filaments

We have calculated the THz radiation distribution further for N = 4, 8, 9, 12 and 16 parallel filaments organized in a square grid in the presence of an external transverse electric field. THz emissions from individual filaments interfere constructively with each other if there is no time delay between the filamentary pulses. The strong enhancement of THz emission from an array of parallel laser filaments in the presence of a transverse static electric field is achieved for all THz frequencies at the same time because the condition of zero time delay needed for the constructive interference applies to every frequency. The total THz energy is channeled toward the direction of the laser. The corresponding THz peak intensity scales like N 2 , as shown in Fig. 7. The 3D distribution of the THz radiation for 16 parallel filaments is presented in Fig. 6(d) as an example. Its peak intensity is found to be ~ 250 times larger than that of a single filament in the presence of the same external field and more than 10 5 higher than that of a single filament without electric field. In principle, with a larger number of filaments, the THz peak intensity may be increased further.

Note that varying of the temporal delay  p leads to significant changes in the radiation diagram. For the example given in Fig. 6(c) for two parallel filaments with temporal delay  p = 1 ps the direction of the radiation at 0.5 THz is completely different from that of the laser beams, forming two lobes on the two sides from the filament direction.

Conclusion

In summary, we have discussed the THz emission from femtosecond laser filaments in air under the combination of two effects: the application of static electric field and the organization of filaments into an array. An enhancement of THz radiation is found, scaling like E² in the presence of a static transverse electric field and scaling like N 2 in the presence of N laser filaments organized in an array. With no temporal delay between the laser pulses forming the filament array, the enhancement is the same at every THz frequency, leading to strong enhancement of a broadband THz radiation in the forward direction. There is also a possibility to direct the THz emission along preferential directions by playing with the time delay between the laser pulses. Dependence of peak THz intensity on the number of filaments N, normalized to the THz intensity of a single filament subject to the same electric field E s > 3 kV/cm. Red squarestheoretical upper limit (quadratic dependence). Black circlesmaximum value obtained in calculations.

Figure captions : Fig. 1 .

 :1 Figure captions:

Fig. 2 .

 2 Fig. 2.Experimental set up. Inset: spatial orientation of the two parallel filaments and direction of the transverse electric field used for calculations and in experiment. The filaments are separated laterally by a variable distance d, and  p is the adjustable time delay between the two filamentary pulses.

Fig. 3 .

 3 Fig. 3. Emission diagrams of the THz radiation at 0.1 THz obtained from two filaments as a function of separation d. The measurements were performed in the plane of the filaments. Top row: calculations. Bottom row: experimental results. All THz intensities are normalized to the peak THz intensity of a single filament.

Fig. 4 .

 4 Fig. 4. Emission diagrams evolution of the THz radiation in the plane of filaments, as a function of the temporal delay between the two laser pulses. The two filaments are separated by d = 2 mm. Top row: calculations. Bottom row: experimental results. THz intensities are normalized to the peak THz intensity of a single filament.

Fig. 5 .

 5 Fig. 5. (a) THz waveform measured by a time domain THz spectroscopy system. (b) Corresponding spectrum of the THz field below 5 THz.

Fig. 6 .

 6 Fig. 6. Calculated THz radiation diagram for the frequency of 0.5 THz from: (a) single filament, (b) two parallel filaments separated by 500 µm without time delay, (c) two parallel filaments separated by 500 µm with a time delay  p = 1 ps, (d) a square matrix of 16 parallel filaments separated by 500 µm without time delay.

Fig. 7 .

 7 Fig. 7.

Figure 1 .

 1 Figure 1.

Figure 2 .

 2 Figure 2.

Figure 3 .

 3 Figure 3.

Figure 4 .

 4 Figure 4.

Figure 6 .

 6 Figure 6.

Figure 7 .

 7 Figure 7.

Acknowledgements

The research leading to these results has received funding from ANR (grant n° ANR-2010-JCJC-0401-01) and from LASERLAB-EUROPE (grant agreement n° 284464, EC's Seventh Framework Programme).