P González De Alaiza Martínez

G Duchateau

B Chimier

R Nuter

I Thiele

Stefan Skupin

V T Tikhonchuk

Maxwell-consistent, symmetry-and energy-preserving solutions for ultrashort laser pulse propagation beyond the paraxial approximation

Keywords: numbers: 42, 42.25.Bs, 42.55.-f

published or not. The documents may come

Maxwell-consistent, symmetry-and energy-preserving solutions for ultrashort laser pulse propagation beyond the paraxial approximation

I. INTRODUCTION

Spatial and temporal pulse shaping makes the laser a highly versatile tool for a large number of applications such as micromachining and material processing [START_REF] Davis | Writing waveguides in glass with a femtosecond laser[END_REF][START_REF] Gamaly | Laser-matter interaction in the bulk of a transparent solid: Confined microexplosion and void formation[END_REF][START_REF] Patel | Non-paraxial polarization spatio-temporal coupling in ultrafast laser material processing[END_REF], Terahertz generation [START_REF] Hebling | Generation of high-power terahertz pulses by tilted-pulse-front excitation and their application possibilities[END_REF][START_REF] Stepanov | Scaling up the energy of thz pulses created by optical rectification[END_REF], or biological imaging and non-invasive surgeries [START_REF] Chung | Surgical applications of femtosecond lasers[END_REF][START_REF] Zipfel | Nonlinear magic: multiphoton microscopy in the biosciences[END_REF]. Paraxial approximation, which assumes that the light angular spectrum is sufficiently narrow, is widely used to study the propagation of laser beams in weak focusing conditions. However, the applications mentioned above usually require tightly focused ultrashort laser pulses. Modeling the propagation of such laser pulses is a complex problem because the deviation from the principal propagation direction is large and the paraxial approximation is no longer valid.

Electromagnetic codes, such as Particle-In-Cell codes [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF][START_REF] Hockney | Computer Simulation Using Particles[END_REF] or codes based on the Unidirectional Pulse Propagation Equation [START_REF] Kolesik | Nonlinear optical pulse propagation simulation: from maxwell's to unidirectional equations[END_REF][START_REF] Kolesik | Unidirectional optical pulse propagation equation[END_REF], are powerful tools for understanding experiments on laser-matter interaction, where laser field components are commonly known in the far field. In these simulations external electromagnetic waves that enter the computational domain are usually prescribed as paraxial modes on the boundaries, which is not adequate for strongly focused ultrashort laser pulses because the initial distorsion may be increased in the course of propagation, leading eventually to nonphysical fields in the simulation box. Therefore, there is a need to determine analytical solutions of Maxwell equations for tightly focused laser pulses.

Different analytical models, restricted to specific beam shapes or spatial symmetry conditions, have been developed to describe nonparaxial laser beam propagation in several physical contexts, such as perturbative expansions of the wave equation [START_REF] Esarey | Theory and group velocity of ultrashort, tightly focused laser pulses[END_REF][START_REF] Li | Fields of an ultrashort tightly focused laser pulse[END_REF], the angular spectrum method [START_REF] Sepke | Analytical solutions for the electromagnetic fields of tightly focused laser beams of arbitrary pulse length[END_REF], transformation optics [START_REF] Yu | Accessing extreme spatiotemporal localization of high-power laser radiation through transformation optics and scalar wave equations[END_REF] or laserdriven electron acceleration based on Helmholtz equation [START_REF] Marceau | Validity of the paraxial approximation for electron acceleration with radially polarized laser beams[END_REF][START_REF] Martens | Direct electron acceleration with tightly focused tm0,1 beams: boundary conditions and non-paraxial corrections[END_REF]. Lax et al. [START_REF] Lax | From maxwell to paraxial wave optics[END_REF] proposed a simple method which paved the way to introduce the nonparaxial corrections to a given paraxial solution in more general situations. They demonstrated that the paraxial solution is actually the zeroth-order consistent solution to the Maxwell equations, obtained by expanding the wave equation (in their case, for Gaussian linearly-polarized vector potential) using a power series in the beam divergence angle.

The nonparaxial perturbative equations proposed by Lax et al. were subsequently analyzed in more detail, always on the basis of the wave equation applied to the vector potential, by several authors for either Gaussian beams [START_REF] Couture | From gaussian beam to complex-source-point spherical wave[END_REF][START_REF] Davis | Theory of electromagnetic beams[END_REF][START_REF] Salamin | Fields of a gaussian beam beyond the paraxial approximation[END_REF][START_REF] Seshadri | Fundamental electromagnetic gaussian beam beyond the paraxial approximation[END_REF][START_REF] Sheppard | Electromagnetic gaussian beams beyond the paraxial approximation[END_REF] or Hermite-Gaussian and Laguerre-Gaussian beams [START_REF] Agrawal | Free-space wave propagation beyond the paraxial approximation[END_REF][START_REF] Agrawal | Gaussian beam propagation beyond the paraxial approximation[END_REF][START_REF] Laabs | Propagation of hermite-gaussian-beams beyond the paraxial approximation[END_REF][START_REF] Takenaka | Propagation of light beams beyond the paraxial approximation[END_REF]. More recently, Favier et al. took into account spatio-temporal couplings in the wave equation in order to extend Lax perturbative equations to few-cycle pulses [START_REF] Favier | Short pulse laser beam beyond paraxial approximation[END_REF]. In the transversespatial and temporal Fourier domain, they linked the Lax series with a Taylor expansion of the exact solution of the wave equation but their proposed high-order corrections hinged on an arbitrary number of integration constants, which were chosen to match some known nonparaxial solutions.

This paper aims at addressing two problems which remain open despite the advances made in the previous works. The first problem is that all the previous approaches solely dealt with the wave equation (in the cited papers, applied to the vector potential) split into a Lax series, and not with the full set of Maxwell equations when calculating high-order corrections. However, the electric and magnetic field components derived from a vector potential that is a rigorous solution of the wave equation do not necessarily satisfy the whole set of Maxwell equations. In other words, since each component of the electric, magnetic and vector potential fields verifies the scalar wave equation, we expect to obtain a unique solution to the Maxwell equations whatever the component chosen to calculate high-order corrections. The second problem is that, when calculating high-order terms of the Lax series from the solutions at lower orders, spurious homogeneous solutions that are not compatible with Maxwell equations may be added through integration constants. We demonstrate in this paper that removing those spurious homogeneous solutions implies preserving the laser energy through transverse planes. Conservation of energy is a fundamental physical principle that, to the best of our knowledge, had never been considered before in the context of nonparaxial corrections. Indeed, determining these integration constants by making ad hoc assumptions on the nonparaxial corrections at the focal point [START_REF] Couture | From gaussian beam to complex-source-point spherical wave[END_REF][START_REF] Takenaka | Propagation of light beams beyond the paraxial approximation[END_REF] or on the beam axis [START_REF] Salamin | Fields of a gaussian beam beyond the paraxial approximation[END_REF], as done in previous works, may be incompatible with the Maxwell equations [START_REF] Lekner | Nonexistence of exact solutions agreeing with the gaussian beam on the beam axis or in the focal plane[END_REF].

In Sec. II our Lax-series-based analytical solution for all electromagnetic field components is presented. Maxwell equations are satisfied and thus each electromagnetic field component verifies the scalar wave equation. The total energy through transverse planes and the existing symmetry between electric and magnetic fields are preserved. Recursive relations to obtain the terms of our Lax series are given and the resulting solution is succesfully compared with a numerically exact Maxwell solver [START_REF] Thiele | Boundary conditions for arbitrarily shaped and tightly focused laser pulses in electromagnetic codes[END_REF]. In Sec. III we calculate the leading term of the asymptotic limit of our Lax-series-based analytical solution far from the focal plane, for both monochromatic beams and ultrashort laser pulses, which results in paraxial-like expressions. These analytical expressions are a baseline for further developements aiming at obtaining an easy and low-computational-cost means of computing the near fields related to those assumed-known paraxial far fields. Thanks to three-dimensional Maxwellconsistent numerical simulations carried out with the code arctic, based on the Yee scheme [START_REF] Yee | Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media[END_REF], we discuss the adequacy of prescribing ultrashort laser pulses by the leading term of the asymptotic limit at a finite distance from the focal plane. Conclusions and outlooks are drawn in Sec. IV.

II. ANALYTICAL SOLUTIONS OF MAXWELL EQUATIONS

A. Maxwell equations and their properties

Maxwell equations in vacuum read as follows:

∇ ∇ ∇ • E E E = 0, (1)
∇ ∇ ∇ • B B B = 0, (2)
∂ t B B B + ∇ ∇ ∇ × E E E = 0, (3)
∂ t E E E -c 2 ∇ ∇ ∇ × B B B = 0, (4)
where E E E and B B B are the electric and magnetic fields, respectively, and c is the speed of light in vacuum.

Maxwell's equations are highly symmetrical and they place the electric and magnetic fields on equal footing [START_REF] Cameron | Electric-magnetic symmetry and noethers theorem[END_REF]. Indeed, both electric and magnetic fields verify the wave equation:

c -2 ∂ 2 t -∇ ∇ ∇ 2 E E E = 0 0 0, (5)
c -2 ∂ 2 t -∇ ∇ ∇ 2 B B B = 0 0 0. (6)
Note that in this paper we formally present our results in vacuum. For monochromatic or narrow-bandwidth pulses, by replacing c by c/n 0 , where n 0 is a constant refractive index, our results generalize to homogeneous dielectric media. Because our solutions are derived in the Fourier space, it would be straightforward to extend it to shorter pulses with linear dispersion.

B. The wave equation

Throughout this paper we shall work in Cartesian coordinates (x, y, z), where x is the optical propagation axis (also referred to as longitudinal axis) and y and z are the transverse coordinates. The beam focus position is placed at x = 0.

We seek solutions of Maxwell equations that are waves propagating along longitudinal axis according to the following Ansatz:

E E E(x, y, z, t) = E 0   ψ Ex (x, y, z, t) ψ Ey (x, y, z, t) ψ Ez (x, y, z, t)   e i(k0x-ω0t) , (7) B B B(x, y, z, t) = E 0 c   ψ Bx (x, y, z, t) ψ By (x, y, z, t) ψ Bz (x, y, z, t)   e i(k0x-ω0t) , (8)
where ω 0 = 2πc/λ 0 is the angular frequency of the laser field, λ 0 is the wavelength, k 0 = ω 0 /c is the wavenumber, ψ Ex , ψ Ey , ψ Ez , ψ Bx , ψ By and ψ Bz are the spatiotemporal envelopes of E x , E y , E z , B x , B y and B z , respectively, and E 0 represents the electric field amplitude. Note that in this paper we only seek forward-propagating solutions along x axis, as stated by Ansätze ( 7) and ( 8), although Eqs. ( 5) and ( 6) admit, in general, bidirectional solutions. Implicitly, we require that E E E and B B B have no evanescent components. Moreover, because they are complex fields, their temporal bandwidth must be restricted to the positive frequency range. By substituting Eq. ( 7) into Eq. ( 5) and Eq. ( 8) into Eq. ( 6), each of the six spatial envelopes, generically denoted as ψ, verifies the so-called wave equation:

∇ ∇ ∇ 2 ⊥ ψ + 2ik 0 ∂ x + ∂ t c ψ = -∂ 2 x ψ + ∂ 2 t ψ c 2 , (9)
where

∇ ∇ ∇ 2 ⊥ = ∂ 2 y + ∂ 2 z .
It is useful to express Eq. ( 9) in the laser co-moving reference system x = x and t = t -x/c:

∇ ∇ ∇ 2 ⊥ ψ + 2ik 0 1 + i∂ t ω 0 ∂ x ψ = -∂ 2 x ψ. ( 10
)
The paraxial approximation neglects the term on righthand side of Eq. ( 10) by claiming that field variations along x axis are small compared to the wavelength λ 0 (i.e., the wavefront is considered to be almost perpendicular to x axis) and to the transverse variations along y and z axes (i.e., the transverse profile is supposed to remain almost unchanged over a distance of the order of λ 0 ). Considering D 0 the 1/e diameter of the Gaussian solution at beam focus (we assume that the diameter is the same along y and z axis) and x R = πD 2 0 /(4λ 0 ) the associated Rayleigh length, we reformulate Eq. [START_REF] Chung | Surgical applications of femtosecond lasers[END_REF] in the dimensionless coordinates ξ = x /x R , τ = ω 0 t , υ = 2y/D 0 and ζ = 2z/D 0 as follows:

∂ 2 ⊥ ψ + 4i T ∂ ξ ψ = -ε 2 ∂ 2 ξ ψ, (11)
where

∂ 2 ⊥ = ∂ 2 υ + ∂ 2
ζ and the operator T = 1 + i ∂ τ introduces the space-time focusing [START_REF] Bergé | Ultrashort filaments of light in weakly ionized, optically transparent media[END_REF][START_REF] Brabec | Nonlinear optical pulse propagation in the single-cycle regime[END_REF]. Equation [START_REF] Couture | From gaussian beam to complex-source-point spherical wave[END_REF] reveals that the term on right-hand side is actually a small correction of order of ε 2 , where ε = D 0 /(2x R ) is the tangent of the beam divergence angle and is assumed to be small in the paraxial limit. For arbitrary spatial beam shapes, for which the Gaussian 1/e beam diameter D 0 does not apply, one can define ε as the angular spectral width divided by k 0 . Note that monochromatic solutions are given by Eq. ( 11) in the limit T → 1, which means that time variations of the envelopes are negligible.

Equation [START_REF] Couture | From gaussian beam to complex-source-point spherical wave[END_REF] in the transverse-spatial and temporal Fourier domain (see Appendix A) reads:

iκ 2 ⊥ 4 T + ∂ ξ - iε 2 4 T ∂ 2 ξ ψ = 0, ( 12
)
where κ 2 ⊥ = κ 2 y + κ 2 z , κ y = D 0 k y /2, κ z = D 0 k z /2, T = 1 + Ω, and Ω = ω/ω 0 . Restricting the temporal bandwidth of the complex fields E E E and B B B to the positive frequency range implies that Ω 1. The exact forwardpropagating solution of Eq. ( 12), with the boundary condition placed at ξ = 0, reads:

ψ(ξ, κ y , κ z , Ω) = ψ(0, κ y , κ z , Ω) e -2i T ε 2 1-1- ε 2 κ 2 ⊥ 4 T 2 ξ , (13)
which, by abuse of language, will be called general solution of the wave equation all through this paper in spite of its lack of bidirectionality.

Equation [START_REF] Davis | Theory of electromagnetic beams[END_REF] discloses that the exact forwardpropagating solution preserves its complex module in all transverse planes:

| ψ(ξ, κ y , κ z , Ω)| = | ψ(0, κ y , κ z , Ω)|, (14)
whenever εκ ⊥ /(2 T ) ≤ 1 (i.e., for propagating modes).

C. The Lax series approach

A Taylor expansion of Eq. ( 13) in powers of κ ⊥ (around κ ⊥ = 0) and ξ (around ξ = 0), reveals that the general solution of the wave equation depends on powers of ε [START_REF] Favier | Short pulse laser beam beyond paraxial approximation[END_REF]. Motivated from this fact, in order to solve Eq. ( 12) one can express ψ in a series using ε as expansion parameter [START_REF] Lax | From maxwell to paraxial wave optics[END_REF]. Because this perturbative approach is a rearrangement of a Taylor expansion, its convergence is thus guaranteed by Taylor's theorem for any ε if highorder terms are calculated as explained below (i.e., satisfying Maxwell consistency, energy preservation, symmetry between electric and magnetic fields and absence of evanescent modes). For linearly-polarized laser pulses, the transverse components (i.e., ψEy , ψEz , ψBy and ψBz , generically denoted as ψ⊥ ) expand in even powers of ε [START_REF] Davis | Theory of electromagnetic beams[END_REF]:

ψ⊥ (ξ, κ y , κ z , Ω) = ∞ j=0 ε 2j ψ(2j) ⊥ (ξ, κ y , κ z , Ω), (15)
whereas the longitudinal components (i.e., ψEx and ψBx , generically denoted as ψ ) expand in odd powers of ε:

ψ (ξ, κ y , κ z , Ω) = ∞ j=0 ε 2j+1 ψ(2j+1) (ξ, κ y , κ z , Ω), ( 16
)
where the functions ψ(2j) ⊥ and ψ(2j+1) have to be determined.

D. Lax series: Splitting the wave equation

If we subsitute Eqs. ( 15) and ( 16) into Eq. ( 12), the wave equation is split into recursive equations.

The lowest order (j = 0) corresponds to the paraxial equation:

iκ 2 ⊥ 4 T + ∂ ξ ψ(0) ⊥ = 0, ( 17
) iκ 2 ⊥ 4 T + ∂ ξ ψ(1) = 0, (18)
where ψ(0)

⊥ = C (0) 0, ⊥ e -i κ 2 ⊥ 4 T ξ and ψ(1) = C (1)
0, e -i κ 2 ⊥ 4 T ξ are, respectively, their solutions. The coefficients

C (0) 0, ⊥ = C (0) 0, ⊥ (κ y , κ z , Ω) and C (1) 0, = C (1)
0, (κ y , κ z , Ω) do not depending on ξ (see their expressions for Hermite-Gaussian and Laguerre-Gaussian beams in Appendix B).

High-order corrections (j > 0) verify:

iκ 2 ⊥ 4 T + ∂ ξ ψ(2j) ⊥ = i 4 T ∂ 2 ξ ψ(2j-2) ⊥ , (19)
iκ 2 ⊥ 4 T + ∂ ξ ψ(2j+1) = i 4 T ∂ 2 ξ ψ(2j-1) , (20)
with the paraxial differential operator in the lefthand side.

We choose to express the solution to Eqs. ( 19) and ( 20) as the sum of a homogeneous solution Ĥ and a particular solution P :

ψ(2j) ⊥ = Ĥ(2j) ⊥ + P (2j) ⊥ , (21)
ψ(2j+1) = Ĥ(2j+1) + P (2j+1) , (22)
where the homogenous solutions are, respectively:

Ĥ(2j) ⊥ = C (2j) 0, ⊥ e -i κ 2 ⊥ 4 T ξ , (23)
Ĥ(2j+1) = C (2j+1) 0, e -i κ 2 ⊥ 4 T ξ , (24)
where the coefficients

C (2j) 0, ⊥ = C (2j) 0, ⊥ (κ y , κ z , Ω) and C (2j+1) 0, = C (2j+1) 0, (κ y , κ z , Ω) do not depend on ξ. It is important to note that even though Ĥ(2j) ⊥ and Ĥ(2j+1)
formally obey the paraxial equation, they are part of the nonparaxial higher-order corrections.

The particular solutions can be writen as:

P (2j) ⊥ = P (2j) ⊥ (ξ)e -i κ 2 ⊥ 4 T ξ , (25)
P (2j+1) = P (2j+1) (ξ)e -i κ 2 ⊥ 4 T ξ , (26)
where the coefficients P

⊥ (ξ) and P (2j+1) (ξ) do depend on ξ. Since in the neighborhood of the focal plane the form ψ ∼ e -i κ 2 ⊥ 4 T ξ dominates, the particular solutions must vanish in that plane, i.e., P (2j)

⊥ (0) = P (2j+1) (0) = 0.
To evaluate them through a recursive procedure as shown below, they can be constructed as j-order polynomials in ξ:

P (2j) ⊥ (ξ) = j k=1 C (2j) k, ⊥ ξ k , (27)
P (2j+1) (ξ) = j k=1 C (2j+1) k, ξ k , (28)
where the coefficients

C (2j) k, ⊥ and C (2j+1) k,
have to be determined.

From the point of view of Lax recursive equations, homogeneous solutions Ĥ are simply arbitrary integration constants and hence Eqs. ( 19) and ( 20) do not suffice to determine them. These homogenous solutions must be determined from Maxwell equations (see Sec. II E). We demonstrate in this paper that such Maxwell-consistent calculation of the high-order corrections ensures that the overall laser energy through transverse planes is not modified by the Lax series terms of order j > 0 and that the Lax series converges towards an exact solution of the wave equation for all electromagnetic components (see Sec. II F). This is a fundamental difference with respect to previous works, where, for example, in order to determine the high-order corrections, some authors had considered ad hoc assumptions such that they are zero at the beam focal point [START_REF] Couture | From gaussian beam to complex-source-point spherical wave[END_REF][START_REF] Takenaka | Propagation of light beams beyond the paraxial approximation[END_REF], they follow the structure of a spherical wave emanating from such beam focal point [START_REF] Salamin | Fields of a gaussian beam beyond the paraxial approximation[END_REF] or they must match some known nonparaxial solutions [START_REF] Favier | Short pulse laser beam beyond paraxial approximation[END_REF]. Indeed, in the particular solutions proposed by most of these works dealing with Hermite-Gaussian and Laguerre-Gaussian paraxial families, spurious homogeneous solutions are found when a Gram-Schmidt orthogonalization process is applied in the focal plane [START_REF] Gram | Über die entwickelung reeller functionen in reihen mittelst der methode der kleinsten quadrate[END_REF][START_REF] Schmidt | Zur theorie der linearen und nichtlinearen integralgleichungen[END_REF]. These spurious modes make the total power through transverse planes increase exponentially with ε [START_REF] Salamin | Fields of a gaussian beam beyond the paraxial approximation[END_REF], which is not physical.

When substituting Eq. ( 21) into Eq. ( 19), and Eq. ( 22) into Eq. ( 20), the following recursion relations are obtained for the coefficients of the particular solutions for all 1 ≤ k ≤ j and j > 0: where, by notation convention,

C (2j) k, ⊥ = - iκ 4 ⊥ 64 T 3 C (2j-2) k-1, ⊥ k + κ 2 ⊥ 8 T 2 C (2j-2) k, ⊥ + i 4 T (k + 1)C (2j-2) k+1, ⊥ , (29)
C (2j+1) k, = - iκ 4 ⊥ 64 T 3 C (2j-1) k-1, k + κ 2 ⊥ 8 T 2 C (2j-1) k, + i 4 T (k + 1)C (2j-1) k+1, , (30)
C (2j-2) k, ⊥ = C (2j-1) k, = 0 if k = j and C (2j-2) k+1, ⊥ = C (2j-1) k+1, = 0 if k ≥ j -1.
It is important to note that the above recursive relations involve the coefficients C (2j-2) 0, ⊥ and C (2j-1) 0, of the homogeneous solution, which will be determined from the Maxwell equations in the following subsection.

E. Lax series: Splitting Maxwell equations

We split Maxwell equations by substituting the Lax expansions ( 15) and ( 16), together with the Ansätze ( 7) and ( 8), into Eqs. ( 1)-( 4).

The envelopes of all the electromagnetic components at paraxial order (j = 0) must verify simultaneously the following overdetermined system of equations:

T ψ(1) Ex + κ y 2 ψ(0) Ey + κ z 2 ψ(0) Ez = 0, (31)
T ψ(1)

Bx + κ y 2 ψ(0) By + κ z 2 ψ(0) Bz = 0, (32)
T ψ(1) Bx -

κ y 2 ψ(0) Ez + κ z 2 ψ(0) Ey = 0, (33)
ψ(0) By + ψ(0) Ez = 0, (34)
ψ(0) Bz - ψ(0) Ey = 0, (35)
T ψ(1)

Ex + κ y 2 ψ(0) Bz - κ z 2 ψ(0) By = 0, (36)
which has a unique solution whatever two components are prescribed [START_REF] Thiele | Boundary conditions for arbitrarily shaped and tightly focused laser pulses in electromagnetic codes[END_REF]. In this paper, without loss of generality, we choose the paraxial-order field polarized along y axis (note that the solution for any other polarization angle can be obtained by applying a rotation transformation):

ψ(0) Ey = C e -i κ 2 ⊥ 4 T ξ , (37)
ψ(0) Ez = 0, ( 38
)
where C(κ y , κ z , Ω) is a coefficient not depending on ξ.

The rest of the components are then calculated from the system (31)- [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF]:

ψ(1) Ex = - κ y 2 T C e -i κ 2 ⊥ 4 T ξ , (39)
ψ( 1)

Bx = - κ z 2 T C e -i κ 2 ⊥ 4 T ξ , (40)
ψ(0) By = 0, (41)
ψ(0) Bz = C e -i κ 2 ⊥ 4 T ξ . (42)
Similarly, the envelopes at high orders (j > 0) must verify simultaneously the following overdetermined system of recursive equations:

2i T ψ(2j+1) Ex + iκ y ψ(2j) Ey + iκ z ψ(2j) Ez = -∂ ξ ψ(2j-1) Ex , (43)
2i T ψ(2j+1)

Bx + iκ y ψ(2j) By + iκ z ψ(2j) Bz = -∂ ξ ψ(2j-1) Bx , (44)
T ψ(2j+1)

Bx - κ y 2 ψ(2j) Ez + κ z 2 ψ(2j) Ey = 0, (45)
ψ(2j) By + ψ(2j) Ez = κ z 2 T ψ(2j-1) Ex + i 2 T ∂ ξ ψ(2j-2) Ez , (46)
ψ(2j) Bz - ψ(2j) Ey = - κ y 2 T ψ(2j-1) Ex - i 2 T ∂ ξ ψ(2j-2) Ey , (47)
T ψ(2j+1)

Ex + κ y 2 ψ(2j) Bz - κ z 2 ψ(2j) By = 0, (48)
ψ(2j) Bz - ψ(2j) Ey = κ z 2 T ψ(2j-1) Bx + i 2 T ∂ ξ ψ(2j-2) Bz , (49) ψ
(2j) By + ψ(2j) Ez = κ y 2 T ψ(2j-1) Bx + i 2 T ∂ ξ ψ(2j-2) By , (50)
which allows us to calculate the homogeneous parts in Eqs. ( 21) and [START_REF] Kolesik | Unidirectional optical pulse propagation equation[END_REF]. Note that the particular solutions calculated in Sec. II D satisfy all equations (43)-(50). In order to determine a unique homogeneous solution, we have to account for the symmetry existing between the electric and magnetic fields. For forward-propagating linearly-polarized pulses, by observing Eqs. ( 46) and (50) and Eqs. ( 47) and (49), we require that:

Ĥ(2j) Ez = Ĥ(2j) By , (51)
Ĥ(2j) Ey = - Ĥ(2j) Bz , ( 52
)
which indeed is the opposite situation to the paraxial order (compare to Eqs. ( 34) and ( 35)). A posteriori, we will demonstrate in Sec. II F that this symmetry condition prevents high-order corrections from modifying the total laser energy. After some manipulations, taking into account that we have prescribed the transverse electric field as in Eqs. [START_REF] Takenaka | Propagation of light beams beyond the paraxial approximation[END_REF] and [START_REF] Thiele | Boundary conditions for arbitrarily shaped and tightly focused laser pulses in electromagnetic codes[END_REF], we get the following homogeneous solution for orders j > 0:

C (2j+1) 0, Ex = κ 2 ⊥ 16 T 2 C (2j-1) 0, Ex + i 4 T C (2j-1) 1, Ex , (53)
C (2j) 0, Ey = κ 2 z 8 T 2 C (2j-2) 0, Ey - κ 2 ⊥ 16 T 2 C (2j-2) 0, Bz - κ y κ z 8 T 2 C (2j-2) 0, Ez - i 4 T C (2j-2) 1, Bz , (54)
C (2j) 0, Ez = κ 2 y 8 T 2 C (2j-2) 0, Ez + κ 2 ⊥ 16 T 2 C (2j-2) 0, By - κ y κ z 8 T 2 C (2j-2) 0, Ey + i 4 T C (2j-2) 1, By , (55)
C (2j+1) 0, Bx = κ 2 ⊥ 16 T 2 C (2j-1) 0, Bx + i 4 T C (2j-1) 1, Bx , (56)
C (2j+1) 0, By = C (2j) 0, Ez , (57)
C (2j+1) 0, Bz = -C (2j) 0, Ey , (58)
where C

(1)

1, Ex = C (1) 1, Bx = C (0) 1, By = C (0)
1, Bz = 0 by notation convention.

In conclusion, by setting C in Eqs. ( 37) and ( 38) the nonparaxial solution can be calculated in the whole space thanks to the recursive formulae Eqs. ( 29)- [START_REF] Salamin | Fields of a gaussian beam beyond the paraxial approximation[END_REF] and Eqs. ( 53)-(58). By way of example, the correction at order j = 1 reads:

ψ(3) Ex = κ 2 ⊥ 16 T 2 - iκ 4 ⊥ 64 T 3 ξ ψ(1) Ex , (59) ψ
(2) Ey = κ 2 z -κ 2 y 16 T 2 - iκ 4 ⊥ 64 T 3 ξ ψ(0) Ey - κ y κ z 8 T 2 ψ(0) Ez , (60)
ψ(2) Ez = κ 2 y -κ 2 z 16 T 2 - iκ 4 ⊥ 64 T 3 ξ ψ(0) Ez - κ y κ z 8 T 2 ψ(0) Ey , (61)
ψ(3) Bx = κ 2 ⊥ 16 T 2 - iκ 4 ⊥ 64 T 3 ξ ψ(1) Bx , (62)
ψ(2) By = κ 2 y -κ 2 z 16 T 2 + iκ 4 ⊥ 64 T 3 ξ ψ(0) Ez - κ y κ z 8 T 2 ψ(0) Ey , (63)
ψ(2) Bz = - κ 2 z -κ 2 y 16 T 2 - iκ 4 ⊥ 64 T 3 ξ ψ(0) Ey + κ y κ z 8 T 2 ψ(0) Ez . (64)

F. Properties of the solution

The overall laser energy is calculated by integrating the longitudinal component of the Poynting vector (Π x ) over transverse coordinates and time (see Appendix C): where Π x = c 2 ε 0 (E y Bz -E z By ). In terms of the inner product between envelopes defined by Eq. (C5), the normalized total energy 4ω 0 U/(cε 0 E 2 0 D 2 0 ) can be expanded in powers of ε as follows:

U = D 2 0 4 ω 0 +∞ -∞ Π x dυ dζ dτ, (65)
ψ Ey , ψ Bz -ψ Ez , ψ By = ∞ j=0 ε 2j j α=0 ψ (2α) Ey , ψ (2j-2α) Bz - ∞ j=0 ε 2j j α=0 ψ (2α) Ez , ψ (2j-2α) By . ( 66
)
We shall demonstrate that our solution, calculated according to Secs. II D and II E, conserves this total energy. That is to say, the paraxial level sets the energy and highorder corrections do not modify it:

ψ Ey , ψ Bz -ψ Ez , ψ By = ψ (0) Ey , ψ (0) Bz -ψ (0) Ez , ψ (0)
By .

(67) Proving Eq. ( 67) is equivalent to prove that all the ε (2j) -order terms of Eq. (66) are zero, for j ≥ 1. Using Fourier-backtransformed Eqs. (60), (61), ( 63) and (64), we verify that the ε 2 -order term is zero:

ψ (0) Ey , ψ (2)
Bz + ψ (2) Ey , ψ (0) Bz - ψ (0) Ez , ψ (2)
By + ψ (2) Ez , ψ (0) By = 0, (68)
because, in virtue of the properties of the inner product given by Eqs. (C9), (C10) and (C12), the first two terms cancel one another:

ψ (0) Ey , ψ (2)
Bz + ψ (2)
Ey , ψ

(0) Bz = -ψ (0) Ey , ∂ 2 υζ 8 T 2 ψ (0) Ez + ∂ 2 υζ 8 T 2 ψ (0) Ez , ψ (0) Ey =0 + ψ (0) Ey , ∂ 2 ζ -∂ 2 υ 16 T 2 ψ (0) Ey - ∂ 2 ζ -∂ 2 υ 16 T 2 ψ (0) Ey , ψ (0)
Ey =0 + ψ (0)
Ey , -

iξ∂ 4 ⊥ 64 T 3 ψ (0) Ey + - iξ∂ 4 ⊥ 64 T 3 ψ (0) Ey , ψ (0)
Ey =0 = 0, (69)
as well as the last two terms (not shown). We have checked that the higher-order terms of Eq. (66) are also zero, suggesting that Eq. ( 67) is satisfied. Indeed, numerical simulations confirm this energy conservation.

Moreover, as already explained in Sec. II C, the Taylor theorem guarantees us that our j-order solution, for any laser component (E x , E y , E z , B x , B y or B z ), converges, up to an error O(ε 2j+2 ), towards the exact solution of the wave equation, given by Eq. ( 13), in complex module and complex argument. The former is invariant according to Eq. ( 14) and the later, by doing a Taylor series of Eq. ( 13) in κ ⊥ / T , expands in powers of ε as follows:

arg ψ(ξ, κ y , κ z ) -arg ψ(0, κ y , κ z ) = - κ 2 ⊥ ξ 4 T - κ 4 ⊥ ξ 64 T 3 ε 2 + O(ε 4 ). ( 70
)
For instance, we verify the convergence for E y component at order j = 1. From Eqs. ( 15), ( 37), [START_REF] Thiele | Boundary conditions for arbitrarily shaped and tightly focused laser pulses in electromagnetic codes[END_REF] and ( 60), on the one hand we have that:

| ψEy (ξ)| | ψEy (0)| = 1 + - iξκ 4 ⊥ ε 2 64 T 3 1 + (κ 2 z -κ 2 y )ε 2 16 T 2 = 1 + O(ε 4 ), (71)
and on the other hand we have:

arg ψEy (ξ) -arg ψEy (0) = - κ 2 ⊥ 4 T ξ + arctan     - ξκ 4 ⊥ ε 2 64 T 3 1 + (κ 2 z -κ 2 y )ε 2 16 T 2     = - κ 2 ⊥ 4 T ξ - κ 4 ⊥ ξ 64 T 3 ε 2 + O(ε 4 ). ( 72
)
The previous developments have shown the reliability of the proposed analytical solution.

Finally, it is worth mentioning that our solution does converge in the propagating region, given by εκ ⊥ /(2 T ) ≤ 1, but it may diverge in the evanescent one. This comes from the fact that the quotient between a high-order correction and the previous one depends on power of εκ ⊥ /(2 T ), as one observes from the recursive formulae ( 29)-( 30) and ( 53)-(58). This ratio is ≤ 1 in the propagating region and hence the convergence of the series is guaranteed.

G. Example: Monochromatic Hermite-Gaussian and Laguerre-Gaussian beams

We confront our Lax-series-based analytical solution to a numerical algorithm computing Maxwell-consistent solutions [START_REF] Thiele | Boundary conditions for arbitrarily shaped and tightly focused laser pulses in electromagnetic codes[END_REF] (see Appendix D). To do so, monochromatic beams are considered (i.e., T → 1) and the coefficient C(κ y , κ z , Ω) in Eqs. ( 37)-(42) shall refer here to either a Hermite-Gaussian beam (see Eq. (B8)) or a Laguerre-Gaussian beam (see Eq. (B16)). Since the Lax series originates from a Taylor expansion around the beam focus, the best way to proceed is to prescribe our analytical solution in the focal plane, truncated at different orders j, and subsequently measure, for all electromagnetic components, the error between the solution of exact solver ( ψsolver ) and our analytical solution ( ψ) in different transverse planes. We compute errors using the standard Frobenius norm. The local relative error in a transverse plane is quantified as:

e = e(ξ, κ y , κ z ) = k 0 ψsolver -ψ k 2 ⊥ ≤k 2 0 ψsolver 2 dk y dk z , (73)
and the global relative error in the same plane is:

E = E(ξ) = k 2 ⊥ ≤k 2 0 ψsolver -ψ 2 dk y dk z k 2 ⊥ ≤k 2 0 ψsolver 2 dk y dk z . (74)
Figures 1 and2 show the analytical solution built from a (1, 1)-order Hermite-Gaussian and (1, 1)-order Laguerre-Gaussian modes, respectively, in the transverse plane placed at ξ = 1. We take λ 0 = 800 nm and a moderate ε = 0.25 (for which the evanescent power is negligeable). The highest local relative error (see Eq. ( 73)) appear in a ring (i.e., high values of transverse wavenumbers). When increasing the truncation order of the Lax series, this ring becomes narrower and the errors reduce in absolute value (not shown). This confirms numerically the convergence of our Lax-series-based solution seen as a Taylor expansion in κ y and κ z . Figure 3 shows that the global relative error diminishes too in all transverse planes when increasing the truncation order. This also confirms numerically the convergence of our Lax-seriesbased solution seen as a Taylor expansion in ξ.

III. ASYMPTOTIC BEHAVIOR FAR FROM FOCAL PLANE

Let us assume that, following Eqs. ( 37)-( 42), the transverse field components at the paraxial order in the position space are:

ψ (0) Ey = ψ (0) , (75)
ψ (0) Ez = 0. ( 76
)
The paraxial mode ψ (0) expands in the limit ξ → ±∞ as:

ψ (0) = 1 ξ N a 0 + a 1 ξ + a 2 ξ 2 + • • • , ( 77
)
where N > 0 is the leading exponent of the asymptotic limit (in general, N is not necessarily an integer), which implies that a 0 = 0, and all coefficients a j = a j (υ, ζ, τ ) do not depend on ξ. Equation (77) verifies the paraxial equation (B1):

∂ 2 ⊥ a 0 ξ N + ∂ 2 ⊥ a 1 -4i T N a 0 ξ N +1 + ∂ 2 ⊥ a 2 -4i T (N + 1)a 1 ξ N +2 + • • • = 0, (78)
from where we deduce that:

∂ 2(j+1) ⊥ a j = 0, ( 79
)
for all j ≥ 0.

In the scope of this paper, we aim at calculating the asymptotic limit where ξ → ±∞ of our solution, generically denoted as ψ ∞ , at the leading term O(ξ -N ). The particular solutions of high-order corrections for E y and E z components, given by Eq. ( 29), vanish by virtue of Eq. ( 79) at the leading order O(ξ -N ). Hence, only the homogenous solutions of high-order corrections may contribute to the limit ξ → ±∞, given by Eqs. ( 54) and ( 55), at such leading order O(ξ -N ). After some manipulations, the limits for the transverse components are, respectively:

ψ ∞ Ey ψ (0) ∼ 1 + 1 a 0 ∞ j=1 ε 2j A (2j) Ey , ( 80
)
ψ ∞ Ez ψ (0) ∼ 1 a 0 ∞ j=1 ε 2j A (2j) Ez , (81)
where, from Eqs. ( 60) and ( 61) we obtain for j = 1: and, from Eqs. ( 54), ( 55) and ( 79) we obtain the following recursive formulae for j > 1:

A (2) Ey = (∂ 2 υ -∂ 2 ζ ) 16 T 2 a 0 , (82)
A (2) Ez = ∂ 2 υζ 8 T 2 a 0 , (83)
A (2j) Ey = - 1 8 T 2 ∂ 2 ζ A (2j-2) Ey -∂ 2 υζ A (2j-2) Ez , (84)
A (2j) Ez = - 1 8 T 2 ∂ 2 υ A (2j-2) Ez -∂ 2 υζ A (2j-2) Ey . ( 85
)
From Eqs. ( 80)-( 85) we see that the leading terms of the limits of E y and E z where ξ → ±∞ hinge upon the dominant coefficient a 0 in Eq. ( 77). The series in Eqs. ( 80) and (81) must be truncated at order ∼ O(ξ -N ). These limits are first specified below for monochromatic (i.e., T → 1) Hermite-Gaussian (Appendix E 1) and Laguerre-Gaussian (Appendix E 2) families. Then, these limits are calculated with a time envelope coupled to Hermite-Gaussian (Appendix E 3) and Laguerre-Gaussian beams (Appendix E 4).

Our analytical method is able to link the nonparaxial near fields to the far fields. These leading terms of the asymptotic expressions in many cases are simply the paraxial-level term of the Lax series, as shown in Appendices E 1-E 4, and usually mimic quite well experimental conditions, e.g., a broad beam on a focusing mirror.

In order to carry out accurate simulations under highly nonparaxial conditions using Finite-Difference-Time-Domain (FDTD) Maxwell codes, one should compute the fields on the focal plane (where the series converges faster) using our analytical approach and then propagate it to the boundaries using numerically exact Fourier solvers [START_REF] Thiele | Boundary conditions for arbitrarily shaped and tightly focused laser pulses in electromagnetic codes[END_REF]. This would be less expensive than calculating directly the high-order corrections of the Lax series in a transverse plane far from the focal plane with a sufficient precision, as one can see in Fig. 3. Nevertheless, if the boundaries where fields need to be prescribed are very distant from the focal plane (several Rayleigh lengths), which means that the transverse-spatial and temporal windows are very large, spectral algorithms would demand considerable computational ressources. In this situation, prescribing directly the asymptotic limit of the Lax series on boundaries far enough from the focal plane may be computationally interesting.

We verify our analytical results with three-dimensional (3D) Maxwell-consistent numerical simulations performed using the code arctic [START_REF] Thiele | Broadband terahertz emission from two-color femtosecond-laser-induced microplasmas[END_REF]. Maxwell equations are discretized by means of Yee scheme [START_REF] Yee | Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media[END_REF]. The simulation domain is delimited by Bérenger's Perfectly-Matched-Layer (PML) absorbing boundary condition [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Bérenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF]. The laser is injected via E y , E z , B y and B z components in the transverse plane placed right after the PML according to the Total-Field/Scattered-Field technique [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF], as shown in Fig. 4. This boundary is placed several Rayleigh lengths from the beam focus x = x f . The origin of the optical axis (x = 0) is set at the position of the left boundary. The input paraxial-order Gaussian pulse at x = 0 is characterized by its 1/e beam diameter D and numerical aperture (NA). The numerical aperture (0 ≤ NA ≤ 1) of a Gaussian beam is defined as the sine of its divergence angle. Our Lax series expansion parameter ε = (D 0 /2)/x R , that is, the ratio of the 1/e beam radius at focus D 0 /2 and the Rayleigh length x R , represents the tangent of the beam divergence angle. Thus, expressed in terms of ε, the numerical aperture reads:

NA = ε √ 1 + ε 2 . ( 86
)
The beam focal plane situates at x = x f :

x f = λ 0 πε 2 πεD 2λ 0 2 -1. ( 87
)
If injecting the leading term of the asymptotic expansion of our solution at x = 0 (i.e., a simple Gaussian beam) instead of directly injecting the full solution, then the position of the focal plane obtained with the full Maxwell solver may differ from x f .

We simulate a y-polarized Gaussian laser beam at the wavelength λ 0 = 0.8 µm coupled, according to Eq. (E22), with the Gaussian time envelope:

C τ (Ω) = τ p 2 √ π e - τ 2 p Ω 2 4 = F τ e -τ 2 τ 2 p , (88)
where τ p = 16.99 fs is the 1/e duration and thus the Full-Width-at-Half-Maximum (FWHM) duration of the pulse (envelope of intensity) is 20 fs (i.e., 7.49 optical cycles). The 1/e beam diameter at x = 0 is D = 7.31 µm.

We take a numerical aperture of NA = 0.57 in the vacuum, which gives ε = 0.7 corresponding to strong focusing conditions where the nonparaxial regime is completely stablished [START_REF] Salamin | Fields of a gaussian beam beyond the paraxial approximation[END_REF]. The beam focal plane should situate at x f = 5.20 µm according to Eq. (87). Since x R = 0.52 µm, the prescription plane is x f /x R = 10 Rayleigh lengths far from the beam focal plane. An overall input energy of 36 nJ is considered (which corresponds to E 0 = 55.36 GV/m). The resolution chosen in arctic is ∆x = 31.8 nm (25 points per wavelength), ∆y = ∆z = 63.7 nm (13 points per wavelength) and ∆t = 84.9 as (31 points per period). The PML layer is ten cells wide in each direction.

Figure 5 shows the maximum value of |E y | over time in the XY plane (i.e., z = 0). The laser pulse is prescribed according to Eq. (E22). The temporal inverse Fourier transformed is computed from Eq. (A6) using the 64-points Gauss-Legendre quadrature formula in the frequency interval -10/τ p ≤ Ω ≤ 10/τ p . The error of 22% between the position of the beam focal plane given by arctic (x = 4.06 µm) and the theory (x f = 5.20 µm) is due to the fact that only the leading term of the asymptotic solution is taken into consideration. The previous evaluations have been performed within conditions of very tightly focused pulses (ε = 0.7). By decreasing ε to 0.5, the error drops to roughly 10% (not shown), which is acceptable. Therefore, it turns out that prescribing the laser fields at a finite distance implies a contribution of next-to-leading orders of the asymptotic expansion: the smaller ε, the smaller the next-to-leading order contribution.

The maximum values in the focal plane over time of the module of the electric field components are comparable with the module of the corresponding spatial envelopes with T → 1. These latter values, |E 0 ψ|, computed up to order j = 18 using the recursive formulae (37)-( 39) and ( 53)-(55) in the position space, are depicted in Fig. 6(d-f). Note that one must apply a filter to the solution of the Lax series because it does not converge for evanescent modes (see Sec. II F). The peak of E y predicted by our Lax-series-based solution is 47.80 GV/m, which is lower than the peak of the paraxial-order term of the series (E 0 = 55.36 GV/m), as illustrated in Fig. 7, due to the strong focusing conditions. The cuts in the focal plane of the simulation corresponding to Fig. 5 are shown in Fig. 6(a-c). One observes that the results of arctic and our analytical solution qualitatively agree but the amplitudes are ∼20% smaller (the peak of E y is 37.80 GV/m).

IV. CONCLUSION AND OUTLOOKS

Both the wave equation and the paraxial equation possess an infinite number of solutions. In this pa-per, we have demonstrated that from any paraxial solution we can build, in a self-consistent fashion, an exact solution of the wave equation for the six electromagnetic field components, assuming forward-propagating linearly-polarized laser pulses, which, by the first time to best of our knowledge, is consistent with Maxwell equations, conserves the energy transported through transverse planes and preserves the symmetry between the electric and magnetic fields. To do so, we have split, following the procedure by Lax et al. [START_REF] Lax | From maxwell to paraxial wave optics[END_REF] and in the transverse-spatial and temporal Fourier space, both the scalar wave equations applied to each electromagnetic field component and to the Maxwell equations. Highorder corrections have been separated in a homogeneous solution and a particular solution. The particular solution is integrated directly from the wave equation. The homogeneous solution, instead, must be calculated so that the whole set of Maxwell equations is satisfied and the existing symmetry between the electric and magnetic fields is preserved. Only then the total laser energy through transverse planes is conserved. We give simple recursive relations in order to obtain these Maxwellconsistent high-order corrections, which are polynomials on the longitudinal coordinate whose coefficients are paraxial modes related to transverse-spatial and temporal derivatives of the paraxial-order term of the Lax series. The convergence of our Lax-series-based solution towards the exact solution of the wave equation is demonstrated and then numerically illustrated.

Since in experiments fields are usually known far from the focal plane, we have derived the leading term if the asymptotic expansion of the full analytical solution of the Maxwell equations. In the case of a strongly focused 20-fs-FWHM Gaussian laser pulse, numerical simulations confirm the reliability of this asymptotic expression up to an accuracy of 10%. Further developments for next-toleading orders are expected to decrease this error.

Hermite-Gaussian modes

The Hermite-Gaussian modes are a well-known complete family of orthogonal paraxial solutions:

ψ (HG) n,m (F, υ, ζ) = (2F -1) m+n n! m! 2 n+m H n √ 2 F υ √ 2F -1 H m √ 2F ζ √ 2F -1 F e -F ρ 2 , ( B3
)
where n is the order of the Hermite polynomial H n along y axis, m is the order along z axis, ρ 2 = υ 2 + ζ 2 . Hermite polynomials verify:

H n (x) = 2xH n-1 (x) -2(n -1)H n-2 (x), (B4) H n (x) -2xH n (x) + 2nH n (x) = 0, ( B5
)
where accounts for the derivative with respect to the variable x and the first two polynomials are H 0 (x) = 1 and H 1 (x) = 2x. Hermite-Gaussian propagation modes are orthogonal between one another, with the inner product defined by Eq. (C7):

ψ (HG) n,m , ψ (HG) p,q = ψ (HG) p,q , ψ (HG) n,m = π 2 δ p n δ q m , (B6)
where δ p n refers to Kronecker delta function. In the transverse-spatial Fourier space,the (n, m)-order Hermite-Gaussian mode reads:

ψ(HG) n,m = C (HG) n,m e -i κ 2 ⊥ 4 ξ , (B7)
where:

C (HG) n,m = (-i) n+m 4π √ n! m! 2 n+m H n κ y √ 2 H m κ z √ 2 e -κ 2 ⊥ 4 .
(B8) Transverse derivatives of Hermite-Gaussian modes can be expressed as a linear combination of Hermite-Gaussian modes:

iκ y ψ(HG) n,m = - √ n + 1 ψ(HG) n+1,m + √ n ψ(HG) n-1,m , (B9) iκ z ψ(HG) n,m = - √ m + 1 ψ(HG) n,m+1 + √ m ψ(HG) n,m-1 , (B10) where, by notation convention, √ n H n-1 (x) = 0 if n = 0.

Laguerre-Gaussian modes

The Laguerre-Gaussian modes are a well-known complete family of orthogonal paraxial solutions:

ψ (LG) p,l (F, υ, ζ) = (2F -1) p ( √ 2 F ) |l| (p + |l|)! p! (υ + sgn(l) iζ) |l| L |l| p 2ρ 2 F 2 2F -1 F e -F ρ 2 , (B11)
where p ≥ 0 is the radial index and l is the azimuthal index (it can be negative, zero or positive integer) of the generalized Laguerre polynomial L |l| p , and sgn(l) is the sign of l, i.e., sgn(l) = 1 if l ≥ 0 and sgn(l) = -1 if l < 0. Generalized Laguerre polynomials verify:

L |l| p (x) = (2p + |l| -1 -x)L |l| p-1 (x) -(p + |l| -1)L |l| p-2 (x) p , (B12)
where the first two polynomials are

L |l| 0 (x) = 1 and L |l| 1 (x) = 1 + |l| -x.
Laguerre-Gaussian propagation modes constitute an orthogonal set:

ψ (LG) p,l , ψ (LG) q,r = ψ (LG) q,r , ψ (LG) p,l = π 2 δ q p δ r l . (B13)
The Gaussian beam belongs to both Hermite-Gaussian and Laguerre-Gaussian families:

φ (HG) 0,0 = φ (LG) 0,0 . (B14)
In the transverse-spatial Fourier space, the (p, l)-order Laguerre-Gaussian mode reads:

ψ(LG) p,l = C (LG) l,p e -i κ 2 ⊥ 4 ξ , (B15)
where:

C (LG) p,l = (-i) 2p+|l| √ p! 4π 2 |l| (p + |l|)! (κ y + sgn(l) iκ z ) |l| L |l| p κ 2 ⊥ 2 e -κ 2 ⊥ 4 . (B16)
Transverse derivatives of Laguerre-Gaussian modes can be expressed as a linear combination of Laguerre-Gaussian modes:

-κ 2 ⊥ ψ(LG) p,l = -2(2p + |l| + 1) ψ(LG) p,l -2 (p + 1)(p + 1 + |l|) ψ(LG) p+1,l -2 p(p + |l|) ψ(LG) p-1,l , (B17)
where, by notation convention,

√ p L |l| p-1 (x) = 0 if p = 0.
whose integral over the transverse coordinates, calculated by employing Ansätze ( 7) and ( 8), gives the power flux through the transverse planes:

P = cε 0 E 2 0 D 2 0 4 +∞ -∞ (ψ Ey ψBz -ψ Ez ψBy ) dυ dζ. (C3)
Integration of Eq. (C3) over time, assuming that there is a time-dependent envelope, gives the total laser energy, which should be the same through any transverse plane:

U = 1 ω 0 +∞ -∞ P dτ. (C4)
The form of the integral in Eq. (C4) suggests us to define the following inner product of spatial envelopes: We shall adapt the exact Maxwell solver in transversespatial Fourier domain of Ref. [START_REF] Thiele | Boundary conditions for arbitrarily shaped and tightly focused laser pulses in electromagnetic codes[END_REF] to the spatial envelopes given in Ansätze [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF] and [START_REF] Brabec | Nonlinear optical pulse propagation in the single-cycle regime[END_REF]. Only the solver for monochromatic laser beams is presented here. To do so, those Ansätze are substituted into the Maxwell equations and we obtain the following overdetermined system: where k x = k 2 0 -k 2 y -k 2 z is the longitudinal component of the wavevector. We only consider forward-propagating modes (i.e., k x ≥ 0), and hence we require ψ(x, k y , k z ) = 0 if k 2 y + k 2 z > k 2 0 . The system (D1) has a unique solution if we assume that the two transverse components of the electric field, E y and E z , are known: The transverse components of the electric field are prescribed in the transverse plane at x = x 0 and propagated according the following expression: ψEy (x, k y , k z ) = ψEy (x 0 , k y , k z )e -i(k0-kx)(x-x0) , (D6) ψEz (x, k y , k z ) = ψEz (x 0 , k y , k z )e -i(k0-kx)(x-x0) , (D7) which is the exact forward-propagating solution of Eq. [START_REF] Davis | Writing waveguides in glass with a femtosecond laser[END_REF]. totically where ξ → ±∞ like:

           k x k y k z 0 0 0 0 0 0 k x k y k z 0 -k z k y -k 0 0 0 k z 0 -k x 0 -k 0 0 k y -k x 0 0 0 k 0 k 0 0 0 0 -k z k y 0 -k 0 0 -k z 0 k x 0 0 k 0 -k y k x 0                    
ψEx = - k y k x ψEy - k z k x ψEz , (D2) ψBx
ψ (HG) n,m ∼ 1 ξ a 0 + O ξ -1 , (E1) a 0 = - i π √ n! m!
Eq. (E21) is straightforward and thus the paraxial mode in position space reads:

ψ (0) = F -1 τ C τ (Ω) ψ (HG) n,m ( F , υ, ζ) , (E22)
where ψ

(HG) n,m is given by Eq. (B3) and:

F = i i -ξ/ T . ( E23
)
Following Sec. E 1, whenever n and m are not simulneously odd integers the asymptotic limits far from the focal plane are:

ψ ∞ Ey ∼ F -1 τ C τ (Ω) ψ (HG) n,m ( F , υ, ζ) , ( E24
)
ψ ∞ Ez ∼ 0. ( E25
)
When both n and m are odd integers:

ψ ∞ Ey ∼ F -1 τ C τ (Ω) ψ (HG) n,m ( F , υ, ζ) , (E26) ψ ∞ Ez ∼ ε 2 8T 2 υζ , F -1 τ C τ (Ω) ψ (HG) n,m ( F , υ, ζ) . ( E27
)
4. Laguerre-Gaussian laser pulses Analogously to Sec. E 3, we prescribe laser field components based on the following paraxial mode that comes from multiplying a time envelope by a Laguerre-Gaussian mode in the focal plane: is given by Eq. (B16) and F is given by Eq. (E23).

ψ (0) = F -1 τ C τ (Ω) ψ
Following Sec. E 2, the asymptotic limits for each transverse laser components are, respectively: (E30)

ψ ∞ Ey ∼F -1 τ C τ (Ω) ψ (

FIG. 1 :

 1 FIG. 1: Lax-series-based analytical solution ψ constructed from the (1, 1)-order Hermite-Gaussian mode, truncated at order j = 5. We consider λ0 = 800 nm and ε = 0.25. The results are plotted in the transverse plane placed at ξ = 1. We show the spatial envelopes for (a) Ex, (b) Ey, and (c) Ez. The corresponding local relative errors, given by Eq. (73), are shown in (d), (e) and (f), respectively.

FIG. 2 :

 2 FIG. 2: Lax-series-based analytical solution ψ constructed from the (1, 1)-order Laguerre-Gaussian mode, truncated at order j = 5. We consider λ0 = 800 nm and ε = 0.25. The results are plotted in the transverse plane placed at ξ = 1. We show the spatial envelopes for (a) Ex, (b) Ey, and (c) Ez. The corresponding local relative errors, given by Eq. (73), are shown in (d), (e) and (f), respectively.

4 FIG. 3 :

 43 FIG.3: Global relative error [Eq. (74)] between our analytical solution and the exact solution for Ey as a function of the longitudinal coordinate at different truncation orders. The Lax series is built from the (1, 1)-order Hermite-Gaussian mode, taking λ0 = 800 nm and ε = 0.25.

FIG. 4 :

 4 FIG.4: Numerical box in arctic (cut at z = 0). The PML region is colored in cyan. The paremeters of the laser beam are defined at x = 0, namely, the 1/e diameter D and the numerical aperture. The injection plane for the Total-Field/Scattered-Field technique is placed right after the PML (vertical black dashed line). The beam focal plane (vertical blue dashed line) is at x = x f . Solid blue lines illustrate the 1/e beam diameter of the corresponding Gaussian pulse.

FIG. 5 :

 5 FIG.5: Maximum value that |Ey| reaches on the XY plane, for a y-polarized 0.8-µm-wavelength 20-fs-FWHM 36-nJ Gaussian laser beam prescribed at x = 0 with NA = 0.57 and D = 7.31 µm. The horizontal white dotted line represents the optical axis. The length of the PML layer is 0.32 µm along x axis and hence the laser pulse is injected at x = 0.32 µm using the paraxial-order term. The vertical white dashed line indicates the beam focal plane position (x f = 5.20 µm) and the white solid lines depict the profile of the Gaussian pulse (paraxial-order term of the Lax series). The leading term of the asymptotic expansion, given by Eq. (E22), is employed to prescribe the laser pulse.

FIG. 6 :FIG. 7 :

 67 FIG. 6: Comparison of arctic's results (maximum values that the module of the electric field components attain in the focal plane) with our analytical solution. Results corresponding to Fig. 5 in the beam focal plane at x = 4.06 µm: Maximum values of (a) |Ex|, (b) |Ey| and (c) |Ez|. Analytical solution in the focal plane at x f = 5.20 µm, calculated from the homogeneous parts given by Eqs. (37)-(39) and (53)-(55) up to order j = 18, with T = 1 and by filtering the evanescent modes: (d) |E0 ψE x |, (e) |E0 ψE y | and (f) |E0 ψE z |.

2 0= 2 0=

 22 a, b := +∞ -∞ a b dυ dζ dτ, (C5)which gives us the total energy of the laser pulse:4ω 0 U cε 0 E 2 0 D ψ Ey , ψ Bz -ψ Ez , ψ By . (C6)Note that for monochromatic beams (i.e., the envelopes do not depend on time) the inner product is defined as:a, b := +∞ -∞ a b dυ dζ, (C7)and, in this case, ψ Ey , ψ Bz -ψ Ez , ψ By represents the total power flux through transverse planes:4P cε 0 E 2 0 D ψ Ey , ψ Bz -ψ Ez , ψ By . (C8)Following the definition of the inner product, if x is a scalar (i.e., it does not depend on υ and ζ), we have that:xa, b = x a, b , (C9) a, xb = x a, b . (C10)Moreover, it follows from the theory of distributions that odd transverse-coordinate and time derivatives are anticommutative and even transverse-coordinate and time derivatives are commutative. For instance:∂ υ a, a = -a, ∂ υ a , (C11) ∂ 2 υ a, a = a, ∂ 2 υ a ,(C12)provided that a(υ → ±∞) = 0 and ∂ υ a(υ → ±∞) = 0.Appendix D: The exact Maxwell solver in the transverse-spatial Fourier space

- 1 τ

 1 C τ (Ω) ψ (LG) p,l ( F , υ, ζ) × iε 2 l(|l| -1) 8T 2 (υ + i sgn(l)ζ) 2 .

Acknowledgments

This research was supported by the project ASTGV (Amplitude Systèmes Through-Glass Via) from French DGA (Direction Générale de l'Armement) funding, and also by the project ELITAS (ELI Tools for Advanced Simulation) CZ.02.1.01/0.0/0.0/16 013/0001793 from European Regional Development Fund. This work was granted access to the HPC resources of TGCC under the allocation A0030506129 made by GENCI and the allocation 2017174175 made by PRACE. S. Skupin acknowledges support by the Qatar National Research Fund (Grant No. NPRP 8-246-1-060).

Appendix A: Definition of the transverse-spatial and temporal Fourier transform

Using the dimensionless coordinates ξ = x /x R , υ = 2y/D 0 , ζ = 2z/D 0 , κ y = D 0 k y /2, κ z = D 0 k z /2, τ = ω 0 t and Ω = ω/ω 0 we define the transverse-spatial and temporal Fourier transform of ψ, denoted as ψ, as the combination of the transverse-spatial Fourier transform (F ⊥ ) and temporal Fourier transform (F τ ):

where the transverse-spatial Fourier transform is:

and the temporal Fourier transform is:

For monochromatic pulses, the temporal Fourier transform defined in Eq. (A5) reduces to a multiplication by a Dirac delta function δ(Ω) in the temporal Fourier space.

Appendix B: Solutions of the paraxial equations

The paraxial equation is:

where the complex longitudinal variable F = i/(i-ξ) has already been introduced by Salamin [START_REF] Salamin | Fields of a gaussian beam beyond the paraxial approximation[END_REF]. By rewriting Eq. (B1) in the transverse-spatial and temporal Fourier space, we can see that the paraxial solution is of the form:

where

Three families of exact solutions for Eq. (B1) are known when T → 1 (i.e., monochromatic pulses): Hermite-Gaussian modes (often called the free-space eigenmodes), Laguerre-Gaussian modes, and Ince-Gaussian modes [START_REF] Bandres | Ince-gaussian beams[END_REF]. Each of these families constitute a countably infinite set of orthogonal paraxial solutions, and they are complete [START_REF] Agrawal | Free-space wave propagation beyond the paraxial approximation[END_REF].

Appendix C: Laser power and energy transported through a transverse plane and definition of the inner product between spatial envelopes

The Poynting vector is defined as:

where the symbol ¯denotes the complex conjugate. Its longitudinal component is:

Appendix E: Examples of asymptotic expansions (leading term)

Monochromatic Hermite-Gaussian beams

If n and m are both even integers, the (n, m)-order Hermite-Gaussian mode (see Eq. (B3)) behaves asymp-(E8)

Whenever n and m are not simulneously odd integers, by substituting Eqs. (E2), (E4) and (E6) into Eqs. (82) and (83) one deduces that the paraxial-order term dominates far from the focal plane:

When both n and m are odd integers, the substitution of Eq. (E8) into Eqs. ( 82) and (83) yields an extra term ∼ O(ξ -3 ) in the asymptotic limit of E z :

Following Eq. ( 67), it is straightforward to verify that these asymptotic limits contain all the power through the transverse plane of the solution.

Monochromatic Laguerre-Gaussian beams

Laguerre-Gaussian modes (see Eq. (B11)) behave asymptotically where ξ → ±∞ like:

where F = i/(i-ξ), ρe ±iφ = υ±iζ, and, in the cylindrical coordinate system, ρ = υ 2 + ζ 2 represents the radial distance and φ is the azimuth (such that υ = ρ cos φ and ζ = ρ sin φ). Note that L |l| p (0) = 0 for all p ≥ 0 and l. After some manipulations, when substituting Eq. (E14) into Eqs. ( 82) and (83) and Eqs. ( 84) and (85), we have:

Ey ,

which are zero if |l| ≤ 1, and for j > 1:

Therefore, the limits for E y and E z are, respectively:

Hermite-Gaussian laser pulses

In the transverse-spatial and temporal Fourier space, we multiply the (n, m)-order Hermite-Gaussian mode in the focal plane (ξ = 0) by a temporal envelope C τ (Ω), in order to prescribe the transverse fields according to Eqs. (75) and (76) with the following paraxial mode: ψ(0) = C τ (Ω) C (HG) n,m (κ y , κ z ) e is given by Eq. (B8). Since by this choice the temporal and transverse-spatial envelopes are separated in the focal plane, the inverse transverse-spatial Fourier transform of


Henri Delebecque Supelec
email: henri.delebecque@supelec.fr

OPTIM: AN OPEN PLATFORM FOR TEACHING INTERACTIVELY WITH MULTIMEDIA

Keywords:

In this paper, we propose an open framework for teachers and lecturers in science, to help them write their pedagogical documents with both static textual parts, and interactive and multimedia content. This framework will also allow them to factorize the structure common to all this documents into classes, and to free themselves from all other design considerations. Moreover, OPTIM allows authors to build interactive parts of these documents without any programming skills. The author can focus his attention to the pedagogical aspects, leaving the GUI generation to our engine. The computational code the teacher supplied, and that provides numerical values displayed by the GUI can be developed in various programming languages.

INTRODUCTION

We propose a general authoring framework, called OPTIM, which uses HDSML, a meta-language more specifically suited for the definition of structured documents than XML [Cowan 2002]. Moreover, OPTIM will allow the author to define interactive objects into pedagogical documents. This is of great interest in the pedagogical context, where students and learners can visualize the concepts described in an alternate graphical manner, while keeping an additional textual description.

OPTIM Principles

The OPTIM acronym means «Open Platform for Teaching Interactively with Multimedia». OPTIM's objective is to give to content authors a framework that allow them to focus on contents (in all its richness) and ignore presentation's problems, while using their favorite editing tool. OPTIM uses HDML [Delebecque 2003] as structuring language that allows the description of the static part of the document. Authors can add interactive parts into their documents using first the OPTIM's GUI generator, and then its ability to link this GUI with the computational code that supplies the information displayed as curves or dynamic schematics.

To allow authors to focus their attention to content, and ignore presentation problems, a structural model is defined by the editor of the whole pedagogical document. This model defines all the mandatory or optional elements content authors have to put into their documents. HDSML is clearly a meta language, since this model is also used to define a markup language that we will call the content authoring language. This language defines tag names that are very familiar to the content authors, since they are defined by the editor, using the words and terms authors currently use.

We will start with some definitions, which help us while presenting the general principles. For more clarity, the following definitions are done in the context of a large pedagogical document, we will generalize hereafter under the global document term.

Definitions

HDSML can be used by authors following three different roles described in the next paragraph, while defining three kinds of documents, which are:

 The class (the structural model),  The presentation model,  The concrete document (a class instance).  The interactive objects to be inserted (Optim objects)

The class designer, who is also the whole pedagogical document editor, defines the structural model common to all the pedagogical sub-documents. The graphics designer, or the ergonomist, designs the look-and-feel used to display these structure elements on a specific presentation medium, in a presentation model. Last, the content author instantiates the class in concrete documents.

The HDSML engine works in three phases, using the four different kinds of document described above.

HDSML Engine

Meta authoring phase

To ensure the uniqueness of the source document, whatever publishing medium is used, we define the notion of pairs of models. For every meaningful subset of the global document (book, chapter, section…), the document editor can supply a specific class (the structural model) using HDSML, and the graphics designer can provide a specific presentation model for any presentation medium. This is done during the first phase, we call the meta-authoring phase.

The HDSML engine takes all these pairs of models, and then automatically generates a XSLT [Adler et al 2001] style-sheets for each one, building a projection function from the class to the presentation model. Moreover the HDSML engine produces the grammar of the authoring language, which will be used by the content author in the next phase.

Authoring Phase

Then, the content authors will write the concrete documents, using any editing tools they want, including not XML compliant ones. The great advantage of using XML compliant editing tools is that they are able to provide context sensitive help. Obviously, an editing tool able to understand the XML Schema [Brown et al 2001] grammar definition (provided by HDSML during the previous phase) can display the attributes which are valid for the tag currently edited, or the children it will accept.

Pre-publication phase

In this phase, the HDSML engine applies the projection function to the concrete documents that conform to the class structure, generating the documents to be published. These later documents are a mix between items supplied by content authors, and others solely bounded to presentation (i.e. background, logo,…) and provided by the presentation model.

Publication phase

In this last phase, the authors can activate the new functionality OPTIM adds to HDSML, and transform their static documents into interactive ones. Since the concrete documents can reference OPTIM objects, these objects will be activated at publishing time by the student. This activation will produce events that are sent to the GUI defined by the author, and compiled as Java© objects by the OPTIM compiler. Every event is either completely managed by the Java code, or transformed into a request sent to the code that the content author supplied, to achieve the computation required. This code can, for example, compute a new 3D-curve to be displayed, the new state of a data path into a processor, or some animation, sound, etc.

OPTIM versatility 1.3.1 Computational language

OPTIM accepts currently Java computational code. But it will be quickly extended to many other domain specific languages, since its working principle requires only basic exchanges between the GUI and the domain-specific code, as explained in the section devoted to the OPTIM architecture. We plan to interface languages like Mathematica and Matlab.

Concrete documents format

HDSML also accepts as concrete document any file not expressed in one of its authoring language at all. We have designed for a special kind of files (namely Microsoft Word © [Microsoft 2004] ones) a specific front-end. It analyzes the content according to grammatical rules, stored in a specific file, called the front-end grammar file. These rules describe how the bookmarks and styles found in the document have to be interpreted to reveal the document internal structure.

One advantage of XML is the great amount of flexibility its tags offer; however, this flexibility also results in a great variety of semantics given to tags by the various authoring languages based on XML, and a lack of standard in this domain. On the contrary, a general authoring language like DocBook © [Walsh et Muellner 1999] have its drawback, since it defines only very general structural entities, like header, footer, body.

We want to help content authors, not always familiar with programming languages, and even less with XML, to write their concrete documents, by defining an authoring language that will use only their domain specific terms. For example, a mathematician will greatly appreciate to find tags named «theorem», «premises», «conclusion», «proof» in the authoring language he/her uses. This is something that people have to code manually if they use languages like Latex.

HDSML Concepts

Let us start with the main structuring entity that both the class designer and the content authors will use: the element.

The Element

All the authors structure their classes or concrete documents by dividing them in meaningful parts they will put into elements. For HDSML, an element is a XML tree node, such as a header or footer structure description (in the class), its content (in a concrete document) or its appearance (in the presentation model, actually using HTML nodes).

The Class

The class definition, (also called the structural model) shared by all the concrete documents, is the editor's responsibility. The actual structure of the concrete documents can be loosely linked to the model that the defines, since the class designer can describe some elements as optional or mandatory, unique or multiple (like using the Kleene operators).

The class design is done customizing the building blocks that HDSML defines (or instantiating them in a object-oriented meaning), to describe the structure common to the global document.

The Presentation Model

The presentation model should probably be designed by a specific author, since it implies more knowhow in ergonomics and graphics than in the normal content domain. This author can choose the presentation language best suited for his/her presentation medium, as long as it is supported by the HDSML engine (currently all the HTML dialects). The author also can use any authoring tool for this task. This model defines the appearance of the all building blocks defines into the class, and that should or can be displayed on that medium.

The Concrete Document

Content authors can freely edit their concrete documents, using their favorite publishing tool, without the burden of sharing a common document. Despite this complete freedom, the global document editor can always compile all the concrete documents in a single one which references them, and publish it after the projection by the HDSML engine.

Another great advantage of this concrete document set is its uniqueness, even if it has to be published on multiple media, thanks to projection functions generated by the HDSML engine. This uniqueness gives to the content authors the insurance that the information published on every medium will conform to the content they have defined, even if it is very dynamic.

OPTIM Architecture

Every interactive object included into a concrete document should be described in two complementary ways. First, the author have to explain how the graphical sub-objects will answer to end-user stimuli (click, keyboard events,…). On the opposite, the computation code should be able to send refresh orders to the GUI, when new data arrives. These two aspects are managed by a classical MVC architecture evangelized by Java. This very efficient and robust way to separate the Model (in charge of all computational tasks) and the View (that displays the values in the most convenient way for a given context or user) has proved its value since its first definition in the Smalltalk [Goldberg et Robson 1998] object-oriented language.

STRUCTURING ELEMENTS

The class designer builds the class using generic blocks, objects and semantically qualified elements, such as those defining the author, the authoring platform, the document version, or others.

Every element of a class definition inherits the occurrence constraints that HDSML defines. But the class author can modify them, according to the structural freedom he wants to give to his concrete document authors. This way, he will define the occurrence constraints of the content authoring language elements.

The Elements

The Page Element

The page is the XML root node of every HDSML document. A page includes the semantically qualified elements that describe the whole page, and generic ones, that define a common document structure.

The Block Element

The block is the most general element, able to hold other generic ones, like (sub) blocks or objects. The inclusion of semantically qualified elements is not currently supported, since they are theoretically limited to the whole page description.

A block owns a very important «name» attribute, that will give its name to the corresponding tag used by the content author when he writes his concrete document. An optional «type» attribute specifies the kind of block as for example: header, menu, footer, body or form.

The «name» attribute allows the domain centered definition of the content authoring language. This language will be very fine tuned to content authors, since it named the structuring units they will use with their most familiar terms.

Finally, a block can be structured using sub elements taken in the following non exhaustive list:

 Block  Object  Keyword  Import  Optim
The block recursive nature avoids the definition of a specific structuring levels, as defined by some structuring languages like DocBook, while giving it more descriptive power.

The Object Element

The object element is always a leaf of the XML node tree. Its purpose is to describe either non-structured elements (such as picture, sound or movie) or very basically structured elements (such as the «link» or the «field» objects).

The «link» object is a regular association between a text part (the link source) and an URI (the link target). The «field» object is also an association, between a text part (the field name) and an un-typed value (the field value), mostly used in forms. Many other object types exist, and more can be added by class designers.

The Keyword Element

The keyword is a powerful tool for the content author, who can use it to add semantics to potentially every document element.

These fine grain defined keywords are of great utility for information retrieval systems or more generally the Semantic Web Activity [START_REF] Bourda | Henri Delebecque Learning Resources and the Semantic Web Proceedings of Ed-MEDIA 2003[END_REF]] applications. The current full text search tools often return irrelevant documents, since they are unable to use the word definition context for filtering. We all know that the same word can have a lot of meanings, depending on this context: for example, «Washington» can be understood as the last name of a person, a state's name, a city name, …. HDSML allows the definition of keywords in an element context, allowing finer semantically based searches.

The Import Element

The import element allows content authors to factorize document parts (either structured or non structured) in external files, dramatically reducing the global document maintenance cost by removing redundant information. The same feature is of course available to the class designer, to factorize without inheritance, reaching multiple inheritance power.

The Optim Element

The Optim element is the link to interactive objects build using the OPTIM platform. Since the description language of Optim elements is still the authoring language, we keep the content author in a familiar context. Moreover, by allowing the user to supply the computational code in many languages (including those that require specific run-time), we insure the maximal adaptability to most scientific domains.

Example

Let us give an example, with few lines of HDSML defining a class block named «Theorem», followed by an example of one theorem definition in a concrete document. <Block name="Theorem" type="body"> <Object name="Name" type="text" Occurs="1" /> <Object name="Premises" type="text" Occurs="1"/> <Object name="Conclusion" type="text"Occurs="1" /> <Object name="Proof" type="text" Occurs="1" /> <Object name="Keyword" type="text" Occurs="*" /> </Block>

The previous HDSML fragment means that the class designer defines the Theorem generic structure with the following sub-elements:  a mandatory name (which is a text zone),  a mandatory proposition, conclusion and proof,  some optional keywords, as in the following example: <Theorem> <Name>Thales</Name> <Keyword>Geometry</Keyword> <Premises>Let and be two points …</Premises> … </Theorem> We can see that the content authoring language uses domain centered entity names, instead of the generic ones (like block) HDSML defines.

RELATED WORK

HDSML versus Tag Libraries or XSLT

One can think that the powerful capabilities tag libraries like those defined by Struts [Apache group 2000] add to HTML, allowing an easy creation of a coherent set of HTML pages, fulfill our objectives. We think that tag libraries are clearly a good solution, but limited to the HTML production, and probably more oriented towards dynamic HTML. Moreover, they tend to stress on the look-and-feel, when we want to free the authors from such considerations, and stress on semantics and document's structure. CSS share the same drawback.

The XSL sheets offered in the XML tools family, has multiple drawbacks for a content author: he has to fully master XML and XSL. And this mastering has to be regularly exercised for two reasons. First, the author has to write a separate XSL sheet for every presentation medium used, and for every subset of the document. Then, this XSL writing task must be redone every time the user changes the look and feel, and for every presentation medium, and for every subset of the global document.

Moreover, this way of producing web pages doesn't allow the factorization and genericity our tool allows, nor any possibility for handling interrelations between documents. It doesn't provide either any way to manage the document versioning the HDSML engine supports, and which is of great help in managing the production of many contents author.

OPTIM versus current IDEs

One can think that existing Integrated Development Environments, like JBuilder©, or Visual Studio©, are solutions to our problem. We think that these products are really powerful, and maybe too powerful for the rather simple graphical objects authors need. The complexity of the programming task should not be measured in the GUI's design, but rather in the simulation part. And we think that it is very important to leave our content author using his/her favorite programming language in that case. Moreover, IDEs impose often the mastering of many specific programming skills, such as animation techniques, graphical algorithms, which are clearly not trivial, and require a great investment.

E-Learning Development Tools

We can found currently projects, like the one described in [START_REF] Matsuno | Richard Gilbert A Tiered Approach Utilizing Reusable Componentware Methodologies for Multimedia Educational Software Creation and Development Proceedings of[END_REF], that tends to help teachers in designing interactive pedagogical documents. This component-ware based software allows the teacher to create very sophisticated graphical applications without programming. It offers links with databases, and many other services. But it is based on a component-ware architecture, and lacks of a general paradigm, like HDSML. Moreover, this software is currently based on the Delphi© programming language, and seems to have no capability to cooperate with other programming language, for the computational requirements of the simulated objects.

Another project, named SimTool [START_REF] Sassen | Stefan Voß Instructional Design and Implementation of Interactive Learning Tools Proceedings[END_REF], embedded in the more ambitious VORMS project, shares some of our goals. Experiments done using SimTool have proved that students are highly motivated by interactive counterparts to textual presentations. One of the major differences we have found between SimTool and OPTIM is the implementation language, Java, which is, in the case of SimTool, the unique programming language, used for both the GUI construction and the computational tasks. We think that these two complementary parts will be better implemented by different authors. The GUI design involves the mastering of computer animation, graphic algorithms, and the full knowledge of the GUI's object library available. On the opposite, the computational part has to be written by the teacher, which is, in this case, the knowledge reference. Moreover, SimTool is only a tool for building simulation applet, and does not try to help the teacher throughout the whole design of his/her pedagogical document, as HDSML.

Another Hypertext Language?

The author familiar with HTML can have fear to throw away all this heavily acquired knowledge and to learn another hypertext dialect. We have multiple answers.

First, many excellent XML editors exist, for every platform, allowing everybody to find the one that meets its requirements. Such tools can be configured using the HDSML syntax we supplied, allowing a very high learning curve of HDSML, even for beginner editors. Such editing tools tutor users and help them by supplying only the tags allowed in the current context of definition.

The content author can also find all the legal attributes of every element, with a list of their possible values. Moreover, we plan to work on more conversion tools, expanding the way opened with the Word front-end processor already operational [Delebecque 2003]. These tools will be able to extract the structure of documents produced by the most common word processing tools available. Our goal is to simplify authors work and no to ask them to learn a new and strange language.

Is OPTIM universal?

We do not think to solve all the problems that communication between GUI and applicative code can rise, nor pretend to define an universal tool for building interactive applications. We try to describe an architecture that allows an easy separation between GUI and applicative code in an educational context, that have two main characteristics. First the GUI used in a pedagogical context are rather simple, and use a small number of


Hands-on Experience for Teaching Computer A r c h i t e c t u r e H. Delebecque* A b s t r a c t

Teaching computer architecture in our School of Electrical Engineering (Supélec) includes both lectures and practical projects. Regarding the practical projects (each about 50 hours of work) three different pedagogical computers have been defined. The first one is very simple: it has only eight instructions and two addressing modes. It emphasizes the very basic concepts of computer architecture, and all our first-year students must design and build it after completing a course on the fundamentals of combinational and sequential logic. The second pedagogical computer has a more complex architecture and all the senior computer science majors have to design and build either a PDP/11-like or a stack-based microprogrammed processor. The third processor our majors have to design and build is a pipelined RISC processor, including their choice of trade-offs between hardware and software as part of their task. The building part of the projects, which have all to end up with a fully working processor, is possible in the rather short time available by using a pedagogical hardware toolkit which has been developed in our lab for that purpose. This toolkit has been designed to allow students to concentrate on the functional characteristics of digital devices by reducing the practical constraints induced by real circuits to a minimum (fan-in and fan-out, short-circuits, wiring, lay-out, etc.). We are currently developing more advanced add-on devices for the toolkit to allow students to explore further concepts of modern architectures, such as the cooperation of processors in multiprocessor arrays. More information on the contents of the projects, the reactions of the students and the design of the toolkit are given below.

Pedagogical Context

"Computer Architecture" is a two-year general course which is compulsory for all students at Supélec, and includes both lectures and practical sessions. The first year course helps students to assimilate the basic concepts of computer architecture from the hardware and software points of view.

Students become familiar with the usual topics of combinational and sequential logic, and acquire knowledge in high level and assembly programming languages.

The second year course lets students master the fundamental physical entities of computer architecture, namely busses, memory chips and modules, interrupts, and their software counterparts. All students ending their second year at Supélec design and build a pedagogical computer, named OPIP, during a lab session which will be described later.

During their third and last year, all students have to specialize in one of fourteen majors, each in a particular field (e.g. communications systems, power electronics, computer science…). Students who choose to major in computer science will have to master the main fields of modern computer systems, through theoretical lectures and lab work. For that purpose, they have to complete two hardware oriented projects. The first lab-project is devoted to the design and building of a PDP-11-like microprogrammed processor, and the second one is dedicated to a pipelined RISC processor.

Pedagogical aims

We have designed the three lab sessions in Computer Science with the four following pedagogical goals in mind. The first two of them are specific to the Computer Architecture domain, and the two others represent general principles of education in our school of engineering

Students should master the computer's complexity

One of the main purposes of the first lab session is to make students able to design and build a fully working computer within a one-day session. During this session they learn that a computer is neither a "completely closed black box", nor a "magic device", but rather only a complex and comprehensible automaton.

To achieve this goal, we prpose the students to build a very small but realistic computer : the OPIP. It is small since it has only eight instructions, and a 4K-word memory module. However, it remains realistic, because it is based on the same principles as most modern computers.

Students have to build their OPIP using a small set of hardware modules, each devoted to a simple functionality: Arithmetic and Logical Unit, data storage for the memory and register modules, etc. Each module has been designed to remain easily understandable with the theoretical knowledge gained during the lectures in Computer Science. This allows students to keep a clear overall understanding of their computer during the whole building phase, and thus to master easily its behaviour, even in its finest details.

Lab sessions present architectures of increasing c o m p l e x i t y

Computer architecture lectures present computers in increasing order of complexity and the lab sessions follow this progression. Nevertheless, the basic building blocks remain the same, for economical and pedagogical reasons. For this reason, we limit the diversity of the hardware modules, which reduces design and servicing costs and we present our students an homogeneous set of tools, able to synthesize and debug the three kinds of processors with the same paradigm. To achieve these goals, we have focused on modularity and clearness in the design of the hardware toolkit involved in these projects.

Students should use their theoretical knowledge

During lab sessions, we ask students to apply the notions presented in lectures, such as the interpretation of an instruction, the effect of an addressing mode, the purpose and behaviour of the various registers. These lab sessions are complementary to lectures since students learn to solve particular but real problems with the methods taught in theoretical lectures, which in turn helps them to better understand the concepts.

These lab sessions are a good example of "horizontal" projects, where students have to use concepts acquired in "vertical" lectures (e.g in digital electronics, computer architecture, programming languages).

A minimal training time

The time devoted to lab sessions in Computer Science during the two first years is limited to about 50 hours, because the global curriculum has to cover the whole scope of electrical engineering. This time is divided into two modules, one per year. The first one allows students to become familiar with assembly programming languages, whereas the second is the OPIP lab session. To fulfil this requirement, we had to design a toolkit and projects so as to minimize both the training phase, and the time spent by students in auxiliary operations (wiring of modules, loading and downloading the memory, etc).

Lab Material

All these practical sessions use a pedagogical hardware toolkit, called the ESE1000, which has been designed and built at The Computer Science Department. ESE1000 presents modules such as registers, arithmetic and logical unit, or memory in a pedagogical way: inputs and outputs are continuously displayed in hexadecimal, power supply connections are hidden, outputs are protected against short-circuits. Moreover, ESE1000 features modules which are as universal as possible, allowing it to support the three lab sessions mentioned above with only eight kinds of modules.

The Arithmetic and Logic Unit, register and memory modules are common to the three architectures. The ALU module allows sixteen different operations (addition, subtraction, one's or two's complement…). The register module includes a 16-bit wide register (usually devoted to the role of accumulator), and a 16 bit-wide Program Counter. Finally, a 4K-word memory module allows dual access through two independent ports and it has built-in latches for address and incoming data, allowing students to minimize data transfers between the memory and the other modules. There is also a front panel with sockets for Programmable Logic Devices. This panel is able to control the other modules through the equations loaded into PLDs by the students. The OPIP and RISC projects both use a multiphase clock, but the CISC architecture requires a microprogrammed sequencer.

Pedagogical results

The OPIP lab session is proposed to our 360 students, which first design a basic but rather slow processor, and then optimize it to suppress useless phases of instruction interpretation. The OPIP is also given as a one-day vocational training workshop. In both cases, 90% of the students complete the minimal requirement, and 70% the optimized one. Moreover, former pupils asked to mention one of their favourite lab session very often select the OPIP, even if they dislike hardware. They stress the advantage of remaining physically in contact with the computer they are building, and thus understanding everything that happens.

The success rate is lower for the CISC project. We have observed that the length of this project should better be restrained to the only four afternoons available in a week (the remaining time is reserved for lectures and sport). Often the work which is left for the following week is unsatisfactory, because students seem to lose track of what they were working on. Moreover, the programming environment supplied until now leads to lengthy loading and downloading operations which slow down the debugging phase. To overcome that difficulty, we are currently developing a powerful front-end, based on the X-Windows interface, and supported by workstations (Macintoshes or DEC VaxStations)


Hang Zhou
email: hang.zhou@recherch.enac.fr

Rémi Douvenot

Alexandre Chabory
email: alexandre.chabory@recherche.enac.fr

Comparison of N×2D and 3D Split-Step Fourier Methods in Realistic 3D Ducting Conditions

Keywords: 3D propagation, discrete formulation, troposphere propagation, ducts

A 3D split-step Fourier method is presented and tested. It is based on a spectral representation of the discrete wave equation in cylindrical coordinates. The formulations in this method are self-consistent in the discrete domain. This method is applied to simulate the propagation with ducts retrieved after clutter data inversion. In the numerical tests, the difference between N × 2D and 3D models is not significant. Thus, in this case, considering the 3D atmosphere effects is unnecessary.

INTRODUCTION

The split-step Fourier method is widely used to simulate the wave propagation in 2D vertical planes [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF]. To model a 3D scene, a N × 2D model approximation scanning all azimuth angles is classically used. However, N×2D models neglect lateral effects. These latter appear in the case of irregular reliefs and/or non-constant refractive index. Accurate 3D models are required for the wave propagation in complex environments. A 3D split-step Fourier method based on a spectral representation of the discrete wave equation in cylindrical coordinates has been proposed by Zhou et al. [START_REF] Zhou | A 3D Split-Step Fourier Algorithm Based on a Discrete Spectral Representation of the Propagation Equation[END_REF]. This formulation is selfconsistent in the discrete domain. Obtained in cylindrical coordinates, this method is also more suitable for radar propagation than those in Cartesian coordinates. The method can handle a slowly varying inhomogeneous atmosphere using phase screens.

In this paper, this method is applied to a realistic 3D ducts. The ducts over the sea are obtained by inversion of clutter data measurements [START_REF] Douvenot | On the knowledge of radar coverage at sea using real time refractivity from clutter[END_REF]. In the numerical test, both simulations using 3D and N×2D methods are compared.

II. CONFIGURATION

The cylindrical coordinate system (r,θ,z) is used. The sources are assumed inside the cylinder r ≤ r0 and the fields/potentials are known at r = r0. A perfect metallic and infinite planar ground is located at z = 0. The propagation is computed in the region r > r0, z > 0. The fields are represented by TE and TM components with respect to z by means of Hertz potentials. In this paper, only the TE component is considered. The TM component could be simulated in a similar way.

For numerical applications, the computation domain is discretized and of finite size. Thus, the vertical domain is limited to z ∈ ]0,z max [ and the following uniform grid is used:

with ∆r = r max /Nr , ∆z = z max /Nz and ∆θ = 2π/Nθ.

III. 3D SPLIT-STEP FOURIER METHOD

The computation is performed marching on in distance. The potential Ψ on the cylinder at distance r 0 + p r ∆r is expressed in the spectral domain by Fourier transforms. We propagate the potential through a homogeneous medium using the spectral propagator Λ, such that [START_REF] Zhou | A 3D Split-Step Fourier Algorithm Based on a Discrete Spectral Representation of the Propagation Equation[END_REF] with r = r 0 + p r ∆r, where and with

Then the spectrum is expressed in the spatial domain with inverse Fourier transforms. This propagator differs from the one classically used [START_REF] Siegmann | A wide-angle threedimensional parabolic wave equation[END_REF], because it is obtained from a discretized wave equation.

In the spatial domain, we apply a phase screen to account for the refraction index variations. To remove reflections over the top boundary, an apodization is applied with a Hanning window on the top half area in the spatial domain. In the final step, the field is calculated from the potential.

IV. NUMERICAL TEST

A. Refractive index model

To highlight the 3D effects of the atmosphere, a refractive index model in the presence of a duct with azimuthal variations is required. Since we have no measured data of the refractive index with strong azimuthal variations, we choose the ducts retrieved after clutter data inversion by Douvenot et al. [START_REF] Douvenot | On the knowledge of radar coverage at sea using real time refractivity from clutter[END_REF]. In their work, measured clutter data with Spandar radar [START_REF] Ingwersen | Radars for ballistic missile defense research[END_REF] of frequency 2.84GHz are inversed to infer the refractive index along a 90° azimuthal aperture. The retrieved ducts are plotted in Fig. 4. of [START_REF] Douvenot | On the knowledge of radar coverage at sea using real time refractivity from clutter[END_REF]. We focus on the 24° azimuthal sector where the duct has the stronger azimuthal variation. The maximum gradients are 0.91 Munits/m in vertical and 0.98 M-units/deg in azimuth.

B. Configuration of the test

In the simulation, a 2D complex source point [START_REF] Deschamps | Gaussian beam as bundle of complex rays[END_REF] at height 30.78 m is used with the same frequency as Spandar radar. The waist width of the source is chosen as 0.5 m. The center of the source is at azimuth 150°. The simulation parameters are r0 = 2 km, r max = 60 km, ∆r = 0.2 m, z max = 204.8 m, ∆z = 0.2 m. The sectoral propagation method [START_REF] Zhou | A 3D Split-Step Fourier Algorithm Based on a Discrete Spectral Representation of the Propagation Equation[END_REF] is applied in a sector from 138° to 162° , discretized on Nθ = 10000 points. The sea surface is modelled as a perfect electric conductor in the test.

C. Propagation test

The normalized electric field obtained at distance r max = 60 km is shown in Fig. 1. Due to the atmosphere ducts over the sea, we can see a deformation of the final field. The difference between the normalized final fields obtained with the 3D and the N×2D models are shown in Fig. 2. The difference is negligible. The difference of the electric field normalized at each range on the horizontal plane at the height of the source between 3D and N×2D methods is shown in Fig. 3. The difference is also negligible.

The proposed 3D model is theoretically more accurate than the N×2D model. However, the difference here is not significant because the azimuthal gradient of refractivity is minor. N × 2D methods should be preferred, since 3D method is resource intensive. This simulation is performed on a standard desktop computer. The computation time of the 3D model is about 6 hours, and the one of N×2D model is about 1 hour (Nθ = 10000).

V. CONCLUSION

An application of the 3D split-step Fourier method based on a spectral representation of the discretized wave equation in cylindrical coordinates has been presented with realistic ducting conditions. The ducts over a sea retrieved after clutter data inversion have been used. Both simulations using 3D and N × 2D methods have been performed to compare the difference. The 3D model has a better accuracy since lateral effects are considered. However, in ducting conditions, the azimuthal gradient is not significant. Hence, the use of N×2D methods to model this type of scenarios is justified.

Fig. 1 .

 1 Fig. 1. Normalized final electric field (dB) using the proposed 3D method.

Fig. 2 .

 2 Fig. 2. Difference of normalized final electric fields (dB) between 3D and N×2D methods.

Fig. 3 .

 3 Fig. 3. Difference of normalized electric fields (dB) on the horizontal plane at height of the source between 3D and N×2D methods.