M D Croitoru

A I Buzdin

Soliton formation in the FFLO phase

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

If the superconducting state in high magnetic fields is paramagnetically limited, the phase transition from a homogeneous superconducting state to a normal state may occur either directly, at the Pauli limiting field H P = 0 / √ 2 with 0 superconducting gap, or via an interposed wedged spatially inhomogeneous state, characterized by the simultaneous coexistence of normal and superconducting phases and known as the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state [START_REF] Zwicknagl | BCS: 50 Years[END_REF]. The latter is predicted for clean spin-singlet superconductors in the vicinity of the upper critical magnetic field H c2 , beyond the field set by the Pauli paramagnetic limit as a result of the competition between pairing correlations, favoring antiparallel spin alignment, and the Zeeman effect, favoring parallel spin alignment along the field, when the Pauli pair breaking dominates the orbital pair breaking [START_REF] Buzdin | [END_REF], with the limiting field H orb c2 = 0 /2πξ 2 . In this state the order parameter varies periodically in space and Cooper pairs acquire a nonzero momentum. Another requirement for the FFLO state is that the superconductor should be in the clean limit, l ξ 0 = v F /π 0 , where l is the mean free path of the quasiparticles and ξ 0 is the superconducting coherence length at T = 0 and H = 0 [3][4][5]. Furthermore, the highly anisotropic Fermi surface [6,7] favors the FFLO phase formation. This is an inherent property of layered conductors that exhibit a highly anisotropic structure and hence have features of systems with reduced dimensionality [8].

There has been much attention focused on discovering the FFLO state in various superconductors under a high magnetic field since its prediction in 1964 [9,10]. Nevertheless, the conclusive evidence for a nonhomogeneous state has been provided only recently. Several signatures indicating the experimental realization of the FFLO state in organic superconductors were observed for the in-plane external magnetic field when the flux penetrates between the layers in the form of Josephson vortices, thus limiting orbital suppression: (i) the anomaly in the thermal conductivity [11], and finer details in the phase diagram obtained by the tunnel diode oscillator rf penetration depth measurements and pulsed field techniques [12] for the clean organic sample λ-(BETS) 2 GaCl 4 ; (ii) the calorimetric and magnetic torque evidence for the appearance of an additional first-order phase transition line within the superconducting phase in the in-plane high-field regime for the organic sample κ-(BEDT-TTF) 2 Cu(NCS) 2 [13,14]; (iii) evidence for phase transition within the superconducting phase obtained from the local electronic spin polarization in a 13 C NMR spectroscopy experiment [15], NMR detection of spin-polarized quasiparticles [16] forming the Andreev bound states spatially localized in the nodes of the order parameter [6,17,18], and phase transitions in the (H -T ) phase diagram that are consistent with the FFLO phase obtained from rf penetration depth measurements [19] for an organic sample κ-(BEDT-TTF) 2 Cu(NCS) 2 ; (iv) NMR relaxation rate evidence for an additional phase transition line [20] and the clear upturn beyond the Pauli limit in the magnetic field and angulardependent high-resolution specific-heat measurements for the organic material β -(ET) 2 SF 5 CH 2 CF 2 SO 3 [21]; (v) magnetic torque evidence for the tricritical point between the FFLO, homogeneous superconducting, and paramagnetic metallic phases in the 2D magnetic-field-induced organic superconductor λ-(BETS) 2 FeCl 4 [22]; and (vi) measurements of the temperature and angular dependencies of H c2 in pnictide superconductor LiFeAs [23] and KFe 2 As 2 [24] probably suggesting that these systems may have realized the FFLO phase as well.

Similarly, quasi-one-dimensional (1D) superconducting compounds have been studied extensively for exhibiting the inhomogeneous state. The field-amplitude and field-angle dependence of the superconducting transition temperature T c (H ) of the organic superconductor (TMTSF) 2 ClO 4 in a magnetic field applied along the conduction planes have been reported [25]. The authors observed an upturn of the curve of the upper critical field at low temperatures. (An enhancement of almost two times over H P 27 kOe is observed, H c2 50 kOe [25].) Moreover, an unusual in-plane anisotropy of H C2 in the high-field regime was observed. Both observations were interpreted as an evidence of the FFLO state formation [26].

The above overview of the recent experimental developments in hunting inhomogeneous states shows that so far mainly thermodynamic and microscopic due to spin relaxation (NMR) evidences for the existence of an FFLO state have been provided. However, the direct evidence of FFLO states has not been observed yet in experiments. In Ref. [26] we have demonstrated that the temperature evolution of the in-plane critical magnetic field anisotropy as well as of its fine structure may give important information about the FFLO state and unambiguously prove its existence. Our subsequent analysis [27] of the properties of the organic superconductors revealed that in some of them the resonance between the period of the FFLO modulation and the period of the interlayer coupling modulation induced by the external field is possible, which results in anomalous cusps in the temperature and angular dependencies of the in-plane critical field. Furthermore, at the resonances the interplay between the orbital and the paramagnetic effects may result in a commensurability effect between the FFLO and magnetic wave vectors [22,28]. These effects are directly based on the main feature of the inhomogeneous state, the spatial modulations of the order parameter.

In this work we keep studying the effects of orbital contribution to the FFLO modulation, leading to unambiguous proof of the FFLO phase. We study the orbital effect in quasi-1D superconductors under the the influence of the external magnetic field and demonstrate that in the FFLO phase there are two competing length scales: the period of the magnetic field potential and the characteristic length associated with the FFLO modulation [29]. The interplay between these two lengths results in locking phenomena, when the modulation wave vector remains equal to the period (or matches m periods) of the potential through some range of the magnetic field values. To study effects of competing periodicities in the system we use the Frenkel-Kontorova (FK) model [30] and show that the solution in the proximity of the locking transition can be described by cos [Qx + ϕ(x)], with ϕ(x) a phase soliton lattice. The surprising ability of this model to describe many physically important phenomena, such as the dynamics of absorbed layers of atoms on crystal surfaces [START_REF] Floria | [END_REF], charge-density waves [32], ferro-or antiferromagnetics [START_REF] Braun | The Frenkel-Kontorova Model, Concept, Method and Application[END_REF], Josephson junctions chains [START_REF] Lucci | [END_REF], transport properties of vortices [35], and and superionic conductors [36], has attracted much attention in condensed matter physics and nonlinear dynamics. With this work we further extend this series of described phenomena.

II. THEORETICAL MODEL

We consider a quasi-1D conductor with the following electron spectrum:

E p = p 2 x /2m x + 2t y cos(p y d y ) + 2t z cos(p z d z ), (1)
where d y and d z are the interchain distances along the y and z axis, respectively. We assume that the couplings between chains are small, i.e., t z < T c0 and t y < T c0 , but sufficiently large to suppress the charge-and spin-densitywave transitions, to stabilize the superconducting long-range order, and to make the mean-field treatment justified, T 2 c0 /E F t z , T 2 c0 /E F t y [37]. Here T c0 is the critical temperature of the system at H = 0.

In quasi-1D superconductors the orbital effect is extremely weak for the magnetic field applied along the chains [38].

Indeed, near T c0 the parallel (along the x axis) upper critical field is

H x c2 = φ 0 2πξ y ξ z T c0 -T T c0 ,where ξ y = d y t y T c0 and ξ z = d z t z
T c0 , while the perpendicular upper critical field (i.e., along the y axis) is

H y c2 = φ 0 2πξ x ξ z T c0 -T T c0 , with ξ x ∼ v F T c0 . We see that H x c2 H y c2
and the orbital effect for the parallel component of the magnetic field is weakened by the factor

t y T c0 d y ξ x 1.
Therefore, in our analysis of the orbital effect it is enough to take into account only the perpendicular component of the magnetic field because it provides the dominant contribution to the orbital effect. For better clarity we consider the magnetic field to be aligned in the xy plane, making the angle θ with the x axis. The perpendicular component of the magnetic field H y = H sin θ and the corresponding vector potential may be chosen as A z = -xH y = -xH sin θ , A y = 0, A x = 0. The vector potential varies slowly at the interchain distances, because the diamagnetic screening currents are very weak and the magnetic induction is practically uniform and equal to the external field, H . With the help of the Eilenberger equations we obtain the following linearized Ginzburg-Landau equation (for details, see the Appendix):

α (x) =β∂ 2 x (x) -δ∂ 4 x (x) -48t 4 ν sin 4 (Qx) (x) -1 2 t 2 μ[1 -7 cos(2Qx)]Q 2 (x) + 2t 2 μ[3 sin(2Qx)Q∂ x (x) + 3 sin 2 (Qx)Q 2 (x)] -4t 2 γ [1 -cos(2Qx)] (x). (2)
Here Q = -πd z H y /φ 0 = -πd z H sin θ/φ 0 , with

φ 0 = π c/e, t ≡ t z , α ≡ ln T c T c0 -[K 1 (T c ) -K 0 1 (T c )] = T c -T cP AT c
, and α < 0 in the uniform superconducting state. Equation (2) contains the magnetic-field-induced potential, which is periodic in real space with the period λ H = 2φ 0 /d z H y . The corresponding functional reads as

F sn = 1 L x dx α| (x)| 2 + β|∂ x (x)| 2 + δ ∂ 2 x (x) 2 + 1 2 t 2 μ[1 -7 cos(2Qx)]Q 2 | (x)| 2 + 6t 2 μ sin 2 (Qx)|∂ x (x)| 2 + 4t 2 γ [1 -cos(2Qx)]| (x)| 2 , ( 3
)
where the coefficients are given in the Appendix [39]. The functional (3) is obtained for t < v F Q and in this limit the term 48t 4 ν sin 4 (Qx) (x) can be simply neglected.

Without the orbital contribution in 1D, the solution of the Ginzburg-Landau equation is an elliptic sine function (see Appendix) [40], which in the proximity of the line of the transition takes the simple form 0 cos(qx), where q is the absolute value of the modulation wave vector. The linearized Ginzburg-Landau functional (3) and Eq. (2) in the paramagnetic limit (when neglecting the orbital contribution) reduce to

F 1D sn = 1 L x dx α| (x)| 2 + β|∂ x (x)| 2 + δ ∂ 2 x (x) 2 (4)
and

α (x) = β∂ 2 x (x) -δ∂ 4 x (x). (5)
The expression for the modulation wave vector is readily obtained from the solution of the linearized equation in the paramagnetic limit, (x) = 0 cos(qx), q 2 = -β/2δ = 2K 3 (T cP )/K 5 (T cP )v 2 F > 0. This wave vector increases from 0 at the tricritical point to a very large value at T -→ 0.

The orbital effect weakly modifies this solution except for the resonance region, q ≈ Q, where it may strongly influence the FFLO structure [28]. To investigate the influence of the orbital effect on the solution we apply the theory developed by Dzyaloshinskii [41][42][START_REF] Yu | Neutron Diffraction in Long-Periodic Structures[END_REF]. Hence we look for the solution of the type 0 cos [Qx + ϕ(x)], where ϕ(x) is a slowly varying function ϕ (x)

Q. This implies that the orbital effect weakly modifies the exact one-dimensional solution, which is justified for (T c -T )/T c0 (t/T c0 ) 2 . We insert this solution into the functional (3) and perform the average over the period

Q -1 as ... = Q 2π π/Q -π/Q [...]dx, (6)
and taking into account that ϕ(x) Qx, we obtain

F sn Q = 2 0 L x dx α 2 + 1 2 β[Q + ϕ ] 2 + 1 2 δ[Q + ϕ ] 4 + t 2 γ [2 -cos 2ϕ(x)] + 1 8 t 2 μQ 2 [2 -7 cos 2ϕ(x)] + 6 8 t 2 μ[2 -cos 2ϕ(x)][Q + ϕ ] 2 . ( 7
)
In the adopted approximation the functional F sn depends only on the space-varying phase, ϕ(x). Since we are interested in the behavior of the system in the vicinity of the resonance q ≈ Q, we can write the functional in a simpler form (normalized by -πK 5 (T c ) 2 0 , note that K 5 (T c ) < 0)

F ϕ sn Q = 2 0 L x dx -α 2 + 1 8 v 4 F Q 2 (ϕ -δq) 2 + 5 8 t 2 v 2 F Q 2 - 11 16 t 2 v 2 F Q 2 cos[2ϕ(x)], (8)
where δq ≡ Qq is the detuning, or the relative misfit between two periodicities in the system. The phase ϕ(x) is the shift of the FFLO modulation relative to the minima in the potential ∼ cos [2ϕ(x)]. The state ϕ(x) = 0 is the commensurate phase between Q and q. Equation ( 8) is a continuum limit approximation of the Frenkel-Kontorova Hamiltonian introduced by Frank and Van der Merwe [START_REF] Frank | [END_REF].

III. FFLO SOLITON LATTICE

Taking into account the periodicity of the function ϕ(x) we obtain (for the purpose of clarity the expression is normalized by

-πK 5 (T c )v 2 F Q 2 2 0 ) F ϕ sn Q,ϕ = 1 K(κ 2 ) K(κ 2 ) 0 dx v 2 F 8 ϕ ,2 (x) - v 2 F 4 δqϕ (x) + 11t 2 16 cos[2ϕ(x)]. (9)
Therefore, the ground state which minimizes F ϕ sn Q,ϕ , Eq. ( 9), is given by the solutions of the exactly integrable sine-Gordon equation for the phase ϕ(x)

2 ∂ 2 ϕ ∂x 2 + v sin[2ϕ(x)] = 0, ( 10
)
where v ≡ 11t 2 /v 2 F is the effective interlayer coupling parameter. In the absence of this coupling, v = 0, Eq. ( 10) has a solution, ϕ(x) = Cx, that describes unperturbed one-harmonic incommensurate phase, 0 cos (Q + C)x. On the other hand, Eq. ( 10) has a trivial solution, ϕ = 0, that corresponds to a commensurate structure 0 cos Qx. For a finite v Eq. ( 10) should describe some inhomogeneous distribution of the phase of the order parameter. The first integral of this equation is given by

∂ϕ ∂x 2 = v κ 2 {1 -κ 2 sin 2 [ϕ(x)]}, (11)
where κ 2 ≡ 2v/(ε + v) and ε is the constant of integration. Then the exact solution of Eq. ( 11), expressed in terms of the Jacobi amplitude, is

ϕ(x) = am √ v κ x,κ 2 , ( 12
)
where κ is the constant of integration and must be found from the energy minimum of the system. The solution is the regularly spaced solitons, a soliton lattice. The soliton lattice is a compromise between the term cos [2ϕ(x)], which favors ϕ(x) = const, and the derivative part that favors ϕ(x) = Cx.

For κ → 1 it reduces to

ϕ(x) = 2 tan -1 tanh √ v 2κ x , ( 13
)
that describes a domain wall, which separates two commensurate regions, one with phase ϕ = -π/2 and other with ϕ = π/2. Within this wall (soliton) a π change of the phase occurs. The length of the soliton can be small, meaning the fast change of the phase. Substituting the solution ϕ(x) for the order parameter into Eq. ( 8), the normalized energy of the state with the spatial distribution of the order parameter given by Eq. ( 12) is expressed in terms of complete elliptic integrals E an K as

F ϕ sn = δq 4 π 2 √ v κK(κ 2 ) - v 16 1 - 2 κ 2 + 4 E(κ 2 ) κ 2 K(κ 2 ) . ( 14
)
Here it is necessary to take into account the fact that the function K 5 (T c ) < 0 for the region of interest. Minimization of the energy with respect to κ leads to the following equation for κ:

κ E(κ 2 ) = v v c ≡ ζ, ( 15
)
where v c = π 2 δq 2 /4. The numerical solution of this equation is given in Fig. 1. The parameter κ varies in the range 0 κ 1 as long as v varies in the range 0 v v c . In this plot we use the notation ζ = √ v/v c . The variation of the parameter κ results in drastic changes in the behavior of the phase ϕ(x). Figure 2 displays the phase for several values of the parameter ζ ∈ {0.2,0.6,0.8,0.9,0.99,0.999 999}. As ζ is increased a plateau section appears within the period of the soliton lattice L = 2K(κ 2 )κ/ √ v or δqL = (2/π ) 2 E(κ 2 )K(κ 2 ). With increasing ζ , not only the period of the soliton lattice increases but also the relative fraction of the plateau width. On this plateau the phase remains almost constant, but it strongly changes beyond the plateau at the end of the periods. Thus the system can be represented as a periodic structure of domains of the commensurate/locked phase separated by solitons. We can give to Eq. ( 14) another form taking into account the distance L between the domain walls, or since q(κ 2 ) = π 0 /κK(κ 2 ) and h(κ 2 ) = 2E(κ 2 ) 0 /π κ is the parametric definition of the modulation wave vector and the critical field, respectively [45],

F ϕ sn = π 4 δq L - v 16 1 - 2 κ 2 + 4 κ 2 E(κ 2 ) K(κ 2 ) , ( 16
F ϕ sn = δq 8 √ vq(κ 2 ) - v 16 1 - 2 κ 2 + 2q(κ 2 )h(κ 2 ) . ( 17
)
The first term is proportional to the soliton density (inversely proportional to the distance between the solitons) and the second one describes the interaction with the field-induced potential.

One can see that as δq -→ 0 the parameter ζ increases its value. In the limit ζ -→ 1 the period of the soliton lattice diverges, lim

ζ -→1 L = 2κ √ v ln 4 √ 1 -κ 2 , ( 18
)
and the modulation wave vector is negligible. For ζ 1 there is no real valued solution of Eq. ( 15). Consequently, there exists a limiting value of δq c ≡ 2 √ 11t/πv F beyond which the phase is always locked, and the commensurate phase remains stable as δq decreases further. That means that the modulation vector coincides with the modulation induced by the field. Therefore at δq c we have a lock-in of the FFLO modulation to the field-induced potential, i.e., an incommensurate/commensurate-type phase transition. That corresponds to L → ∞ and the FFLO modulation is described by 0 cos Qx. So far we have discussed the quasi-1D case. However, the theory developed is directly applicable to the quasi-2D materials as well. The orientation of the FFLO modulation vector is arbitrary in the pure Pauli-limited s-wave superconductor in the case of a Fermi surface with an elliptical cross section, which can be mapped by scaling transformation to the isotropic case [46]. However, the direction of the FFLO wave vector in real crystals is determined by the crystal field effects [47,48] and/or the pairing symmetry [49]. Therefore, when the in-plane magnetic field is perpendicular to the FFLO wave vector q we reduce the quasi-2D case to the above-described quasi-1D case, Eq. ( 3), but with the renormalized coefficients:

β ≡ π 8 K 3 (T )v 2 F , δ ≡ -3π 128 K 5 (T )v 4 F , μ ≡ -π 4 K 5 (T )v 2 F , γ ≡ π 4 K 3 (T )
. The coefficients γ , ν have the same values. The resonance conditions could be realized, for example, in (TMTSF) 2 ClO 4 with in-plane perpendicular to q H b c2 (0) = 3 T, d z = 1.31 nm, ξ a 0 = 45 nm that result in λ H = 2π/Q = 1050 nm and λ 0 FFLO = 2π/q = π v F / 0 = π 2 ξ a 0 = 444 nm, i.e., λ 0 FFLO < λ 0 H at T = 0. Therefore it guarantees that while approaching a tricritical point the resonance λ FFLO (T ) = λ H (T ) will occur. In the case of the quasi-2D magneticfield-induced organic superconductor λ-(BETS) 2 FeCl 4 the perpendicular H c2 (0) = 4 T, d z = 1.85 nm, ξ 0 = 6.5 nm, and λ 0 FFLO = 65 nm. The superconductivity appears at the in-plane field H 18 T and reaches the maximum critical temperature, T C = 4 K, at H ∼ 32 T. This means that λ H decreases from λ H = 556 nm at T = 0 to 312 nm at T C . This again guarantees the fulfill of the resonance condition in the FFLO phase.

In conclusion: Starting from the Eilenberger equations for a quasi-1D superconductor, which takes into account the orbital effects of an external field, in the proximity to the resonance between the period of the FFLO modulation and the period of the interlayer coupling modulation induced by the external field, we have reduced them to the continuum limit approximation of the Frenkel-Kontorova model equation. This equation is obtained by minimization of the Ginzburg-Landau functional (analogous to the extended Lawrence-Doniach theory for layered superconductors [27]) for an inhomogeneous distribution of the order parameter. As a main result we show that, besides the anomalous cusps in the temperature and angular dependencies of the in-plane critical field, when the FFLO modulation wave vector is in the resonance with the magnetic wave vector, predicted in our previous work, the following occurs: (i) The order parameter is defined by 0 cos [Qx + ϕ(x)], with ϕ(x) describing regularly spaced phase solitons. (ii) The proximity to the resonance condition is characterized by the field-induced phase transition from the incommensurate to the commensurate structure (C-I) of the spatial modulation of the order parameter with respect to the magnetic wave vector structure. The periodic potential tends to "lock-in" the FFLO modulation into a commensurate configuration. (iii) In the incommensurate phase near the C-I transition the periodicity of the phase solitons or phase-domain walls depends continuously on the temperature and external field strength. This soliton lattice is a compromise between two periodicities of the system. The soliton represents a "phase slip" of the order parameter modulation. (iv) The distance between the phase slips logarithmically diverges at critical detuning or misfit, δq c .

It would be promising to study anomalies of the superconducting properties accompanying the field-induced passage of the modulation wave vector of the order parameter through the commensurate value. The formation of the phase soliton lattice could be observed in the NMR spectroscopy measurements similar to those performed in Ref. [15]. The NMR signal, which is proportional to the local electronic spin susceptibility, may be very different in the phase domain wall as compared to the commensurate phase. The appearance of the domain walls in the FFLO phase may enhance its pinning and therefore could be revealed in transport measurements. as

1 λ (x,k ⊥ ) = 2πT Re ω>0 f ω (x,k ⊥ ,p ⊥ ) , ( A5
)
where λ is the pairing constant and the brackets denote averaging over p y and p z and v F x , ... ≡ N (1) π dy with N (1) ≡ m/2πp x d y d z . The temperature unit is so that the Boltzmann constant k B = 1. The maximum critical temperature corresponds to k y = 0, and we have freedom in choosing k z , for example, we can make k z = 0. For the purpose of simplicity, choosing the external magnetic field in the x-y plane, finally we obtain

≡ v F ∂ ∂x + 4it z sin(p z d z ) sin(χ A ), ( A7
)
where χ A ≡ Q x x-k z 2 d. In sequel we make two basic assumptions. First, it is required that the system is near the phase transition that allows expansion of Green functions in powers of /ω n . However, the assumption used is that not only the order parameter is small but also the spatial gradients of the order parameter are small. If we assume that the system resides near the phase transition, then g ω (x,k ⊥ ,p ⊥ ) sgnω n , and this condition gives rise to

n + 1 2 f ω (x,k ⊥ ,p ⊥ ) = , ( A8
)
where n = ω niμ B H . It is convenient to write the procedure as the following iterative scheme:

f (k+1) ω (x,k ⊥ ,p ⊥ ) = n g (k) ω - 1 2 n f (k) ω , (A9) f †(k+1) ω (x,k ⊥ ,p ⊥ ) = * n g (k) ω + 1 2 n f †(k) ω . (A10)
The zero-order approximation (the absence of the superconducting state) results in g

(0) ω (x,k ⊥ ,p ⊥ ) = 1 and f (0) ω (x,k ⊥ ,p ⊥ ) = 0.
Then for the first-order approximation the iterative scheme produces f (1) ω (x,k ⊥ ,p ⊥ ) = / n , f † (1) ω (x,k ⊥ ,p ⊥ ) = * / n . For the iterative scheme we need the iterative expression for the g (k) ω function. If we write it, for example, for g (5) ω , as g (5) ω = g (0) ω + λ g (1) ω + λ 2 ω g (2) ω + λ 3 ω g (2) ω + λ 4 g (4) ω + λ 5 g (5) ω and f (5) (4) ω + λ 5 f (5) ω , then making use of the normalization condition,

ω = f (0) ω + λ f (1) ω + λ 2 f (2) ω + λ 3 f (3) ω + λ 4 f
g 2 ω + f ω f † ω = 1, (A11)
we obtain the following system of coupled expressions for the kth iteration,

g (0) ω (x,k ⊥ ,p ⊥ ) = 1, (A12) k i=0 g (i) ω g (k-i) ω = - k i=0 f (k-i) ω f †(i) ω . (A13)
The second iteration of the iteration scheme leads to

f (2) ω (x,k ⊥ ,p ⊥ ) = (x,k ⊥ ) n - 1 2 2 n (x,k ⊥ ), (A14) f †(2) ω (x,k ⊥ ,p ⊥ ) = * (x,k ⊥ ) n + 1 2 2 n * (x,k ⊥ ). (A15) Therefore here f (2) ω = - (x,k ⊥ )/2 2 n , f †(2) ω = * (x,k ⊥ )/2 2 n .
From the normalization condition we can obtain also

g (1),2 ω + 2 g (0) ω g (2) ω = -f (2) ω f †(0) ω -f (1) ω f †(1) ω -f (0) ω f †(2) ω , (A16) 2 g (0) ω g (2) ω = -f (1) ω f †(1) ω , (A17) which gives g (2) ω = -f (1) ω f †(1) ω /2 = -f (1) ω f †(1)
ω /2 and g (1) ω = 0. Therefore the normal Green function g (2) ω (x,k ⊥ ,p ⊥ ) acquires the following form:

g (2) ω (x,k ⊥ ,p ⊥ ) = g (0) ω + g (1) ω + g (2) ω = 1 - f (1) ω f †(1) ω 2 = 1 - | (x,k ⊥ )| 2 2 2 n .
Following similar steps, then performing averaging over p y and p z and Fermi surface v F x , we obtain on the fifth iteration for f (5) ω ≡ f (5) ω (x,k ⊥ ,p ⊥ ),

f (5) ω = (x,k ⊥ ) n - (x)| (x,k ⊥ )| 2 2 3 n + 1 4 3 n 2 (x,k ⊥ ) + 3 (x,k ⊥ )| (x,k ⊥ )| 4 8 5 n + 1 16 5 n 4 (x,k ⊥ ) - 4 8 5 n | (x,k ⊥ )| 2 2 (x,k ⊥ ) - 1 8 5 n 2 (x,k ⊥ ) 2 * (x,k ⊥ ) - 3 8 5 n * (x,k ⊥ ) (x,k ⊥ ) (x,k ⊥ ) - 2 8 5 n (x,k ⊥ ) (x,k ⊥ ) * (x,k ⊥ ).
Substituting the found averages into the self-consistency relation given by (5) (x,k ⊥ )/λ = πT n f (5) ω (x,k ⊥ ,p ⊥ ) and using the standard regularization rule (A20)

1 λ = ln T T c0 + 2πT ω n >0 1 ω n , ( A18
As shown in Ref. [39], a solution of this equation is the Jacobi elliptic sine function (x) = v F ξ (k) ksn( x ξ (k) ,k), expressed in terms of the modulus k, which is determined upon minimizing the free energy at the fixed external field parameter h. The linearized version of this equation, α (x) = β∂ 2

x (x) -δ∂ 4 x (x) -48t 4 ν sin 4 (χ A ) (x)

- t 2 2 μ[1 -7 cos(2χ A )][∂ x χ A ] 2 (x) + 2t 2 μ[3 sin(2χ A )∂ x χ A ∂ x (x) + 3 sin 2 (χ A )∂ 2 x (x)] -8t 2 γ sin 2 (χ A ) (x), (A21)
is used in the main part of the article.

FIG. 1 .

 1 FIG. 1. Minimization parameter as a function of the anisotropy parameter.

FIG. 2 .

 2 FIG. 2. The phase shift of the FFLO modulation for several values of the anisotropy parameter ζ . The straight line corresponds to an unperturbed incommensurate structure.

Figure 3

 3 illustrates the behavior of the spatially varied order parameter in the vicinity of the IC-C phase transition when ζ → 1. One can see the domain walls (phase solitons) in the spatial distribution of the order parameter.

FIG. 3 .

 3 FIG.3. The spatial distribution of the order parameter in the vicinity of the IC-C phase transition for Q = 0.7 (upper panel) and Q = 1.7 (lower panel).

4 - 4 *

 44 )we obtain the following extended version of the Ginzburg-Landau equation (for clarity we dropped down k ⊥ dependence), if by one understands its three-dimensional version,0 = (x) π K 1 (T ) -K 0 1 (T )ln introduced the notations [39] K m (T ) = 2T ω>0 (ω niμ B H ) -m and K 0 m (T ) = 2T ω>0 ω -m n .This result is similar to the one obtained by Houzet and Mineev[5]. Now we recall the original notations of Eq. (A7) and after the averaging procedure we can write (neglecting the k ⊥ dependence)α (x) = β∂ 2 x (x) -δ∂ 4 x (x) -2γ (x)| (x)| 2 -3ν (x)| (x)| 4 + μ| (x)| 2 ∂ 2 x (x) ∂ x (x)∂ x * (x) -t 2 μ 2 × [1 -7 cos(2χ A )][∂ x χ A ] 2 (x) -80t 2 ν sin 2 (χ A ) × (x)| (x)| 2 + 2t 2 μ 3 sin(2χ A )∂ x χ A ∂ x (x) + sin(2χ A )∂ 2 x χ A (x) + 3 sin 2 (χ A )∂ 2 x (x) -8t 2 γ sin 2 (χ A ) (x) -48t 4 ν sin 4 (χ A ) (x), where α ≡ ln T T c0π[K 1 (T ) -K 0 1 (T )], β ≡ π 4 K 3 (T ) v 2 F x , δ ≡ -π 16 K 5 (T ) v 4 F x , γ ≡ π 4 K 3 (T ), μ ≡ -π 2 K 5 (T ) v 2 F x , ν ≡ -π 8 K 5 (T ).Here, the derivative signs mean the following:∂ 2 (x) = v 2 F x ∂ 2 x (x), ∂ 4 (x) = v 4 F x ∂ 4 x (x), ∂χ A ∂ (x) = v 2 F x ∂ x χ A ∂ x (x),and [∂χA ] 2 = v 2 F x [∂ x χ A ] 2 . Ifone neglects the interchain hopping, then the previous equation reduces to 0 = α (x) -β∂ 2 x (x) + δ∂ 4 x (x) + 2γ (x)| (x)| 2 + 3ν (x)| (x)| μ| (x)| 2 ∂ 2 x (x) -(x)[∂ x (x)] 2 -μ 2 (x)∂ x (x)∂ x * (x).

ACKNOWLEDGMENTS

We acknowledge fruitful discussion with L. Bulaevskii, S. Uji, and J. Wosnitza. This work was supported by the French ANR SUPERTRONICS: Grant No. ANR-15-CE24-0008-02.

APPENDIX: DERIVATION OF THE EXTENDED NONLINEAR LAWRENCE-DONIACH EQUATION

We start from the Eilenberger equations for a quasi-1D superconductor [START_REF] Kopnin | Theory of Nonequilibrium Superconductivity[END_REF]:

where

obtained for an arbitrary direction of the external magnetic field with the chosen gauge (see Sec.

II) resulting in

where α accounts for the field direction, with amplitude H , in the y-z plane from the positive y axis, while θ is the angle accounting the field direction from the positive x axis. This choice leads to Q y = πd y H z /φ 0 = πd y H cos α sin θ/φ 0 , Q x = -πd z H y /φ 0 = -πd z H sin α sin θ/φ 0 , with φ 0 = π c/e. Here k ⊥ = (k y ,k z ), p ⊥ = (p y ,p z ), v F = v F x i is the Fermi velocity along the x axis, n ≡ ω n + 1/2τihsgn(ω n ) with h = μ B H is the Zeeman energy, and we have assumed that the vector potential varies slowly at the inter-chain distances. (This assumption means that we neglect the diamagnetic screening currents and take the magnetic field as uniform and given by the external field H .) These equations depend on k ⊥ , which takes into account the dependence of the order parameter on the center-of-mass coordinate (z + z )/2 and (y + y )/2. The order parameter is defined self-consistently


Axel Modave
email: axel.modave@ensta-paristech.fr

Jonathan Lambrechts

Christophe Geuzaine

Perfectly matched

come

Introduction

Nowadays, the numerical resolution of wave-like problems set on infinite or very large domains remains a challenging task. When using classical schemes based on finite difference, finite volume or finite element methods, a common strategy consists in computing the numerical solution only on a truncated domain, and using an adequate treatment at the artificial boundary to preserve the original solution. This treatment is supposed to simulate the outward propagation of signals and perturbations of all kinds generated inside the truncated domain, even if they are not a priori known. For this purpose, a lot of artificial boundary conditions, artificial layers and alternative techniques have been developed, studied and used for decades (see e.g. the review papers [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF] Colonius | Modeling artificial boundary conditions for compressible flow[END_REF][START_REF] Givoli | Numerical Methods for Problems in Infinite Domains[END_REF][START_REF] Givoli | Computational absorbing boundaries[END_REF][START_REF] Hagstrom | Radiation boundary conditions for the numerical simulation of waves[END_REF][START_REF] Hagstrom | Radiation boundary conditions for Maxwell's equations: a review of accurate time-domain formulations[END_REF][START_REF] Tsynkov | Numerical solution of problems on unbounded domains. A review[END_REF] and references therein). Among them, the high-order absorbing boundary conditions [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF][START_REF] Givoli | High-order local non-reflecting boundary conditions: a review[END_REF][START_REF] Hagstrom | Radiation boundary conditions for Maxwell's equations: a review of accurate time-domain formulations[END_REF][START_REF] Hagstrom | The double absorbing boundary method[END_REF][START_REF] Schmidt | Non-conforming Galerkin finite element methods for local absorbing boundary conditions of higher order[END_REF] and the perfectly matched layers (PMLs) [START_REF] Appelö | Perfectly matched layers for hyperbolic systems: General formulation, well-posedness, and stability[END_REF][START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Bérenger | Perfectly Matched Layer (PML) for Computational Electromagnetics[END_REF][START_REF] Bermúdez | Perfectly matched layers for time-harmonic second order elliptic problems[END_REF][START_REF] Hastings | Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation[END_REF][START_REF] Hu | On absorbing boundary conditions for linearized euler equations by a perfectly matched layer[END_REF][START_REF] Hu | Development of PML absorbing boundary conditions for computational aeroacoustics: A progress review[END_REF][START_REF] Katz | Validation and extension to three dimensions of the Bérenger PML absorbing boundary condition for FD-TD meshes[END_REF] provide treatments of the artificial boundary with arbitrarily-high accuracy.

The PML method has been introduced by Bérenger in the 90s [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] for transient electromagnetic problems, and has been quickly applied to other wave-like problems. With this method, the truncated domain is extended with a layer, where the governing equations are modified in such a way that outgoing waves are damped. In addition, at the interface between the truncated domain and the layer, the outgoing waves are transmitted without any reflection, whatever the angle of incidence. The combination of both these properties made the success of the method.

The key ingredient of the PML method is the set of governing equations defined inside the layer. With Bérenger's strategy, valid for squared and cuboidal truncated domains, the equations are built by splitting the original equations written in Cartesian coordinates, and using specific dissipations terms [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. In the 1 time-harmonic context, this strategy corresponds to a stretch of the Cartesian coordinates in the complex plane [START_REF] Chew | A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates[END_REF][START_REF] Rappaport | Interpreting and improving the PML absorbing boundary condition using anisotropic lossy mapping of space[END_REF], and can be interpreted as a change of the metric [START_REF] Lassas | Complex Riemannian metric and absorbing boundary condition[END_REF][START_REF] Teixeira | Differential forms, metrics, and the reflectionless absorption of electromagnetic waves[END_REF]. The PML equations then involve complex metric coefficients. In alternative formulations, the original time-harmonic equations are used in the PML without any change, but complex anisotropic material parameters are defined [START_REF] Gedney | An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices[END_REF][START_REF] Sacks | A perfectly matched anisotropic absorber for use as an absorbing boundary condition[END_REF][START_REF] Zhao | GT-PML: Generalized theory of perfectly matched layers and its application to the reflectionless truncation of finite-difference time-domain grids[END_REF]. These PMLs, frequently called uniaxial PMLs, are interpreted as anisotropic material absorbers. Time-dependent formulations can finally be obtained by taking the inverse Fourier transform in time of time-harmonic PML equations, and using convolution products or additional fields. All the above mentioned strategies have been developed in the electromagnetics community in the 90s. During the last 20 years, they have been used and adapted to design PML formulations for increasingly complex problems in various physical contexts, such as aeroacoustics [START_REF] Hu | On absorbing boundary conditions for linearized euler equations by a perfectly matched layer[END_REF][START_REF] Hu | Development of PML absorbing boundary conditions for computational aeroacoustics: A progress review[END_REF][START_REF] Hu | Absorbing boundary conditions for nonlinear Euler and Navier-Stokes equations based on the perfectly matched layer technique[END_REF][START_REF] Nataf | A new approach to perfectly matched layers for the linearized Euler system[END_REF], geophysical fluid dynamics [START_REF] Lavelle | A pretty good sponge: Dealing with open boundaries in limited-area ocean models[END_REF][START_REF] Navon | A perfectly matched layer approach to the linearized shallow water equations models[END_REF], elastodynamics [START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF][START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF][START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Chew | Perfectly matched layers for elastodynamics: A new absorbing boundary condition[END_REF][START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Meza-Fajardo | A nonconvolutional, split-field, perfectly matched layer for wave propagation in isotropic and anisotropicelastic media: Stability analysis[END_REF] and quantum mechanics [START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equations[END_REF].

Most of the PML formulations are written in Cartesian coordinates and only deal with truncated domains that have straight artificial boundaries (e.g. squared or cuboidal). However, some problems are naturally written in other coordinate systems and, since the choice of the truncated domains is a priori arbitrary, it could be advantageous to take domains with non-classical shapes. This motivated the derivation of PML systems in alternative coordinate systems, firstly with cylindrical and spherical coordinates [START_REF] Chew | Complex coordinate stretching as a generalized absorbing boundary condition[END_REF][START_REF] Collino | The perfectly matched layer in curvilinear coordinates[END_REF][START_REF] Petropoulos | Reflectionless sponge layers as absorbing boundary conditions for the numerical solution of Maxwell equations in rectangular, cylindrical, and spherical coordinates[END_REF][START_REF] Teixeira | PML-FDTD in cylindrical and spherical grids[END_REF][START_REF] Teixeira | Systematic derivation of anisotropic PML absorbing media in cylindrical and spherical coordinates[END_REF]. Some PML versions dealing with generally-shaped convex domains have been proposed in time-harmonic contexts by Teixeira and Chew [START_REF] Teixeira | Analytical derivation of a conformal perfectly matched absorber for electromagnetic waves[END_REF][START_REF] Teixeira | Complex space approach to perfectly matched layers: a review and some new developments[END_REF], Lassas and Somersalo [START_REF] Lassas | Analysis of the PML equations in general convex geometry[END_REF][START_REF] Lassas | Complex Riemannian metric and absorbing boundary condition[END_REF], Zschiedrich et al [START_REF] Zschiedrich | A new finite element realization of the perfectly matched layer method for Helmholtz scattering problems on polygonal domains in two dimensions[END_REF] and Matuszyk and Demkowicz [START_REF] Matuszyk | Parametric finite elements, exact sequences and perfectly matched layers[END_REF]. Strategies have been presented for time-dependent simulations with non-Cartesian finite difference schemes [START_REF] Roden | Efficient implementation of the uniaxial-based PML media in three-dimensional nonorthogonal coordinates with the use of the FDTD technique[END_REF][START_REF] Teixeira | Conformal PML-FDTD schemes for electromagnetic field simulations: A dynamic stability study[END_REF], mixed finite element schemes [START_REF] Donderici | Conformal perfectly matched layer for the mixed finite element time-domain method[END_REF] and discontinuous Galerkin schemes [START_REF] Alvarez | A leap-frog discontinuous galerkin time-domain method for hirf assessment[END_REF][START_REF] Dosopoulos | Interior penalty discontinuous Galerkin finite element method for the timedependent first order Maxwell's equations[END_REF] for Maxwell's equations. Alternative approaches with layers have also been proposed by Guddati et al [START_REF] Guddati | Continued fraction absorbing boundary conditions for convex polygonal domains[END_REF] and Demaldent and Imperiale [START_REF] Demaldent | Perfectly matched transmission problem with absorbing layers: Application to anisotropic acoustics in convex polygonal domains[END_REF] for polygonal domains.

In this paper, we present two PML formulations for transient acoustic problems defined on convex truncated domains with regular curved boundary. Following strategies used in the electromagnetics community, the formulations are obtained using, respectively, a complex stretch of coordinates and complex material properties in the time-harmonic version of the pressure-velocity system. We derive the tensor form of the final time-dependent equations and we provide the explicit definition of the coefficient tensors, which facilitates implementation with classical discretization methods. Finite element implementations based on a discontinuous Galerkin method are then proposed for both formulations and tested by means of three-dimensional numerical simulations. This paper is organized as follows. In section 2, PML formulations based on the pressure-velocity system are derived in both time-harmonic and time-dependent contexts for generally-shaped convex domains. Section 3 is dedicated to numerical simulations in the time domain. After describing numerical schemes based on a discontinuous Galerkin method, time-dependent formulations are tested and compared by means of a reference three-dimensional benchmark. An illustration of application is finally proposed.

Design of PML formulations

In this section, we derive two families of time-harmonic and time-dependent PML formulations for the acoustic wave system

       ∂p ∂t + ρc 2 ∇ • u = 0, ∂u ∂t + 1 ρ ∇p = 0, (1)
where p(t, x) is the pressure, u(t, x) is the velocity, ρ is the reference density and c is the propagation speed of the medium. The complete original problem consists in finding the fields p(t, x) and u(t, x) that are governed by system (1) for t > 0 and x ∈ R d , with initial conditions given for both fields at t = 0. The spatial dimension d is equal to 2 or 3. For the modified problem, the fields are governed by system (1) only inside the truncated domain Ω ⊂ R d , which is surrounded with the PML Σ (see e.g. Fig. 1a).

The time-harmonic PML formulations are obtained from the original equations by stretching a spatial coordinate in the complex plan, which introduces a directional damping of waves. In order to deal with generally-shaped truncated domains, this stretch is performed in a specific local coordinate system. With this strategy, first used by Teixeira and Chew [START_REF] Teixeira | Analytical derivation of a conformal perfectly matched absorber for electromagnetic waves[END_REF], the truncated domain must be convex with a regularity condition on its boundary, and the PML thickness is constant. The time-dependent formulations are then obtained by defining supplementary fields and applying an inverse Fourier transform in time.

The local coordinate system is described in section 2.1. The time-harmonic and time-dependent PML formulations are derived in sections 2.2 and 2.3, respectively.

Coordinate system associated with the domain boundary

We consider a local coordinate system defined in a layer Σ surrounding the domain Ω with a constant thickness δ. The system is based on lines (in two dimensions) or surfaces (in three dimensions) parallel and perpendicular to the interface Υ (Fig. 1a), which is assumed to be regular enough. This system has been used to derive PMLs [START_REF] Lassas | Analysis of the PML equations in general convex geometry[END_REF][START_REF] Lassas | Complex Riemannian metric and absorbing boundary condition[END_REF][START_REF] Teixeira | Analytical derivation of a conformal perfectly matched absorber for electromagnetic waves[END_REF][START_REF] Teixeira | Complex space approach to perfectly matched layers: a review and some new developments[END_REF] and absorbing boundary conditions [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF].

In two dimensions, the interface Υ is a curve and the coordinate system is denoted (ξ 1 , ξ 2 ). Since Ω is convex, each point x(ξ 1 , ξ 2 ) of the layer Σ has a unique closest point on Υ. We define the coordinate ξ 1 as the distance between the two points, while the coordinate ξ 2 is given by a local parametrization of Υ. We consider a coordinate patch of Υ defined as p : V ⊂ R → Υ. Each point of p(V) ⊂ Υ then is given by p(ξ 2 ), where ξ 2 ∈ V. The coordinate patch is chosen in such a way that

dn dξ 2 = κt, with t = dp dξ 2 ,
where n(ξ 2 ) is the unit outward normal, t(ξ 2 ) is a unit tangent vector and κ(ξ 2 ) is the curvature of Γ at p(ξ 2 ). The coordinates ξ 1 and ξ 2 form an orthogonal curvilinear system, and the set of vectors (e 1 , e 2 ) = (n, t) constitutes an orthonormal frame. For each point of the layer Σ, we can then write

x(ξ 1 , ξ 2 ) = p(ξ 2 ) + ξ 1 e 1 (ξ 2 ),
which is illustrated in Fig. 1b.

The three-dimensional coordinate system (ξ 1 , ξ 2 , ξ 3 ) is adapted from the two-dimensional version. For each point of the layer Σ, the coordinate ξ 1 is the distance with the closest point on the surface Υ, and the coordinates ξ 2 and ξ 3 are provided by a local parametrization of Υ. We consider a coordinate patch p : V ⊂ R 2 → Υ. Each point of p(V) ⊂ Υ then is given by p(ξ 2 , ξ 3 ), where (ξ 2 , ξ 3 ) ∈ V. There exists a coordinate patch that gives (see [START_REF] Antoine | Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape[END_REF][START_REF] Do | Differential geometry of curves and surfaces[END_REF])

dn dξ i = κ i t i , with t i = dp dξ i , for i = 2, 3, (2)
where n(ξ 2 , ξ 3 ) is the unit outward normal, t 2 (ξ 2 , ξ 3 ) and t 3 (ξ 2 , ξ 3 ) are two unit tangent vectors in the principal directions and κ 2 (ξ 2 , ξ 3 ) and κ 3 (ξ 2 , ξ 3 ) are the principal curvatures of the surface Γ at p 2 (ξ 2 , ξ 3 ). The coordinates ξ 1 , ξ 2 and ξ 3 form an orthogonal curvilinear coordinate system, and the set of vectors (e 1 , e 2 , e 3 ) = (n, t 2 , t 3 ) constitutes an orthonormal frame. For each point of the layer Σ, we can then write

x(ξ 1 , ξ 2 , ξ 3 ) = p(ξ 2 , ξ 3 ) + ξ 1 e 1 (ξ 2 , ξ 3 ). (3)
In orthogonal curvilinear coordinates, system (1) can be written

         ∂p ∂t + ρc 2 1 ∏ k h k ( ∑ i ∂ ∂ξ i (( ∏ k̸ =i h k ) u i ) )
= 0,

∂u i ∂t + 1 ρ 1 h i ∂p ∂ξ i = 0, for i = 1, . . . , d (4)
where u i denotes a component of u in the coordinate system, and h i is the scale factor associated with the coordinate ξ i , which is defined by

h i = ∂x ∂ξ i . ( 5
)
The lower and upper bounds of summation and product symbols are 1 and d, respectively. For the sake of clarity, they are not written. Using the definition (5) together with equations ( 2) and (3) gives the scale factors for the three-dimensional coordinate system described in this section,

h 1 = 1, h 2 = 1 + κ 2 (ξ 2 , ξ 3 ) ξ 1 , h 3 = 1 + κ 3 (ξ 2 , ξ 3 ) ξ 1 .

Coordinate stretch and time-harmonic PML systems

Two PML systems are derived for the time-harmonic acoustic wave system

         -ıω p + ρc 2 1 ∏ k h k ( ∑ i ∂ ∂ξ i (( ∏ k̸ =i h k ) ûi ) ) = 0, -ıωû i + 1 ρ 1 h i ∂ p ∂ξ i = 0, for i = 1, 2, 3, (6)
where ω is the angular frequency and the hat ˆdenotes the Fourier transform in time. Both will be used to derive different time-dependent PML formulations in section 2.3.

In a time-harmonic context, a classical way to derive PML systems from the original system consists in stretching one coordinate in the complex plane, where the coordinate corresponds to the direction where waves must be damped in the layer. In our case, this corresponds to replacing the real coordinate ξ 1 ∈ [0, δ] with the complex one ξ1 ∈ U, where U is a curve in the complex plane (Fig. 2). We consider the parametrization of the curve

ξ1 (ξ 1 ) = ξ 1 - 1 ıω ∫ ξ1 0 σ(ξ ′ 1 ) dξ ′ 1 , with ξ 1 ∈ [0, δ], (7)
where σ(ξ 1 ) is the so-called absorption function, which is positive. The effect of the complex stretching with this specific parametrization can be interpreted by considering the plane wave solution. For the original system (6), the plane wave solution reads

e ı(k•x-ωt) ,
where the wave number k is related to the angular frequency through the dispersion relation ω = c ∥k∥.

Replacing ξ 1 with ξ1 in this solution gives with the damping factor

e ı(k•x-ωt) e -γ(ξ1) , ℜe ℑm 0 δ • • (a) Original medium ℜe ℑm U 0 δ • • (b) PML medium
γ(ξ 1 ) = 1 c k • e 1 ∥k∥ ∫ ξ1 0 σ(ξ ′ 1 ) dξ 1 .
Except grazing waves (i.e. k • e 1 = 0), all the plane waves supported by the original system are damped in their direction of propagation. For outgoing plane waves (i.e. k • e 1 > 0), the damping factor is positive and increases when ξ 1 increases. By contrast, for ingoing plane waves (i.e. k • e 1 < 0), the damping factor is negative and increases when ξ 1 decreases.

When replacing ξ 1 with ξ1 in the time-harmonic system (6), the partial derivative with respect to

ξ 1 becomes ∂ ∂ξ 1 → ∂ ∂ ξ1 = 1 1 -σ/(ıω) ∂ ∂ξ 1 ,
and the scale factors h 2 and h 3 , which depend on ξ 1 , become

h 2 → h2 = 1 + κ 2 ξ1 , h 3 → h3 = 1 + κ 3 ξ1 .
The other partial derivatives do not change, nor does the scale factor h 1 = 1. Nevertheless, it is convenient to introduce the complex scale factor defined as h1 = 1 -σ/(ıω). Indeed, the time-harmonic PML system can then be written

         -ıω p + ρc 2 1 ∏ k hk ( ∑ i ∂ ∂ξ i (( ∏ k̸ =i hk ) ûi ) ) = 0, -ıωû i + 1 ρ 1 hi ∂ p ∂ξ i = 0, for i = 1, 2, 3. (8)
An alternative PML system is obtained by defining the new unknowns

û⋆ i = ∏ k̸ =i hk ∏ k̸ =i h k ûi , for i = 1, 2, 3. System (8) then becomes            -ıω p + ρc 2 ∏ k h k ∏ k hk 1 ∏ k h k ( ∑ i ∂ ∂ξ i (( ∏ k̸ =i h k ) û⋆ i ) ) = 0, -ıωû ⋆ i + 1 ρ h i hi Π k̸ =i hk Π k̸ =i h k 1 h i ∂ p ∂ξ i = 0, for i = 1, 2, 3. (9)
Both time-harmonic PML systems ( 8) and ( 9) have interpretations that are well-known in the electromagnetics community. For system (8), the complex coordinate stretch corresponds to a change of the metric of the space. In an orthogonal coordinate system, the original scale factors h i are indeed simply replaced with the complex ones hi (see e.g. [START_REF] Lassas | Complex Riemannian metric and absorbing boundary condition[END_REF][START_REF] Teixeira | Differential forms, metrics, and the reflectionless absorption of electromagnetic waves[END_REF][START_REF] Teixeira | Complex space approach to perfectly matched layers: a review and some new developments[END_REF]). The tensor form of system [START_REF] Bécache | On the analysis of Bérenger's perfectly matched layers for Maxwell's equations[END_REF] reads

   -ıω p + ρc 2 α∇ • û⋆ = 0, -ıω û⋆ + 1 ρ M∇p = 0,
where the scalar α and the second-order tensor M are defined by

α = ∏ k h k ∏ k hk and M = ∑ i h i hi ∏ k̸ =i hk ∏ k̸ =i h k (e i ⊗ e i ) .
This system is identical to the original one (system (1)), except that the (real) material parameters ρc 2 and 1/ρ have been replaced with the complex parameter ρc 2 α and the anisotropic complex tensor M/ρ. The time-harmonic PML can then simply be interpreted as an anisotropic absorber with specific complex material properties [START_REF] Gedney | An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices[END_REF][START_REF] Sacks | A perfectly matched anisotropic absorber for use as an absorbing boundary condition[END_REF][START_REF] Teixeira | Analytical derivation of a conformal perfectly matched absorber for electromagnetic waves[END_REF]. The same material properties have been obtained for the Helmholtz equation by Matuszyk and Demkowicz [START_REF] Matuszyk | Parametric finite elements, exact sequences and perfectly matched layers[END_REF] in a more general framework.

Time-dependent PML systems

Two time-dependent PML systems are obtained from the time-harmonic ones ( 8) and ( 9) by using an inverse Fourier transform in time. This transform can be performed thanks to the introduction of additional fields and equations. For the sake of clarity, we define the real functions

σ 1 = σ(ξ 1 ), σ 2 = κ2 (ξ 1 , ξ 2 , ξ 3 ) σ(ξ 1 ), (10)
σ 3 = κ3 (ξ 1 , ξ 2 , ξ 3 ) σ(ξ 1 ), ( 11
) with κi (ξ 1 , ξ 2 , ξ 3 ) = κ i 1 + κ i ξ 1 , with i = 2, 3, σ(ξ 1 ) = ∫ ξ1 0 σ(ξ ′ 1 ) dξ ′ 1 .
These functions allow us to rewrite the complex scale factors as

hi = ( 1 - σ i ıω ) h i . ( 12
)
Let us firstly consider system [START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF]. By using equation ( 12), the second term of the first equation of system (8) successively becomes

ρc 2 1 ∏ k hk ∑ i ∂ ∂ξ i (( ∏ k̸ =i hk ) ûi ) = ρc 2 1 ∏ k hk ∑ i [ ( ∏ k̸ =i ( 1 -σ k ıω ) ) ∂ ∂ξ i (( ∏ k̸ =i h k ) ûi ) + (( ∏ k̸ =i h k ) ûi ) ∂ ∂ξ i ( ∏ k̸ =i ( 1 -σ k ıω ) ) ] = ∑ i ( 1 -σi ıω ) -1 [ ρc 2 1 ∏ k h k ∂ ∂ξ i (( ∏ k̸ =i h k ) ûi ) ] + ∑ i [ ρc 2 ( ∏ k ( 1 -σ k ıω ) -1 ) ûi h i ∂ ∂ξ i ( ∏ k̸ =i ( 1 -σ k ıω ) ) ] = ∑ i ( 1 -σi σi-ıω ) [ ρc 2 1 ∏ k h k ∂ ∂ξ i (( ∏ k̸ =i h k ) ûi ) ] + ∑ i   ρc 2 (-ıω) σ i -ıω ∑ k̸ =i [ 1 σ k -ıω ∂σ k ∂ξ i ] ûi h i   .
Using the definitions of σ 1 , σ 2 , σ 3 and h 1 , the second term can be simplified as

∑ i   ρc 2 (-ıω) σ i -ıω ∑ k̸ =i [ 1 σ k -ıω ∂σ k ∂ξ i ] ûi h i   = ρc 2 û1 (-ıω) σ 1 -ıω [ 1 σ 2 -ıω ∂σ 2 ∂ξ 1 + 1 σ 3 -ıω ∂σ 3 ∂ξ 1 ] + ρc 2 (-ıω) (σ 2 -ıω)(σ 3 -ıω) [ û2 h 2 ∂σ 3 ∂ξ 2 + û3 h 3 ∂σ 2 ∂ξ 3 ] = ρc 2 û1 (-ıω) σ 1 -ıω [ σ 1 -σ 2 σ 2 -ıω κ2 + σ 1 -σ 3 σ 3 -ıω κ3 ] + ρc 2 (-ıω) (σ 2 -ıω)(σ 3 -ıω) [ û2 h 2 ∂σ 3 ∂ξ 2 + û3 h 3 ∂σ 2 ∂ξ 3 ] = ρc 2 û1 [ - (-ıω) σ 1 -ıω (κ 2 + κ3 ) + (-ıω) σ 2 -ıω κ2 + (-ıω) σ 3 -ıω κ3 ] + ρc 2 (-ıω)σ (σ 2 -ıω)(σ 3 -ıω) [ κ2 3 κ 2 3 û2 h 2 ∂κ 3 ∂ξ 2 + κ2 2 κ 2 2 û3 h 3 ∂κ 2 ∂ξ 3 ] = ∑ i ( 1 -σi σi-ıω ) (q i + σr i ),
where we have introduced the additional fields

q1 = -ρc 2 (κ 2 + κ3 )û 1 , r1 = 0, q2 = ρc 2 κ2 û1 , r2 = ρc 2 1 σ 3 -ıω κ2 3 κ 2 3 û2 h 2 ∂κ 3 ∂ξ 2 , q3 = ρc 2 κ3 û1 , r3 = ρc 2 1 σ 2 -ıω κ2 2 κ 2 2 û3 h 3 ∂κ 2 ∂ξ 3 .
Defining the additional fields pi = -

1 σ i -ıω [ ρc 2 1 ∏ k h k ∂ ∂ξ i (( ∏ k̸ =i h k ) ûi ) + qi + σr i ] , for i = 1, 2, 3, system (8)
can be rewritten as

                                     -ıω p + ρc 2 1 ∏ k h k ( ∑ i ∂ ∂ξ i (( ∏ k̸ =i h k ) ûi ) ) = - ∑ i (σ i pi ) - ∑ i (q i + σr i ), -ıωû i + 1 ρ 1 h i ∂ p ∂ξ i = -σ i ûi , for i = 1, 2, 3, -ıω pi + ρc 2 1 ∏ k h k ∂ ∂ξ i (( ∏ k̸ =i h k ) ûi ) = -(σ i pi + qi + σr i ), for i = 1, 2, 3, -ıωr 2 = ρc 2 κ2 3 κ 2 3 û2 h 2 ∂κ 3 ∂ξ 2 -σ 3 r2 , -ıωr 3 = ρc 2 κ2 2 κ 2 2 û3 h 3 ∂κ 2 ∂ξ 3 -σ 2 r3 .
One equation can be removed by using p = ∑ i pi , as well as one term of the first equation since

∑ i qi = 0.
Using an inverse Fourier transform in time, we finally obtain the time-dependent PML system

                                   ∂p ∂t + ρc 2 ∇ • u = -σ 1 p 1 -σ 2 p 2 -σ 3 (p -p 1 -p 2 ) -σ(r 2 + r 3 ), ∂u ∂t + 1 ρ ∇p = - ∑ i σ i u i e i , ∂p i ∂t + ρc 2 ∇ i • u = -σ i p i -q i -σr i , for i = 1, 2, ∂r 2 ∂t = ρc 2 u 2 κ2 3 κ 2 3 (e 2 • ∇κ 3 ) -σ 3 r 2 , ∂r 3 ∂t = ρc 2 u 3 κ2 2 κ 2 2 (e 3 • ∇κ 2 ) -σ 2 r 3 . ( 13
)
with

u i = e i • u, q 1 = -ρc 2 (κ 2 + κ3 )u 1 , q 2 = ρc 2 κ2 u 1 and ∇ i • u def. = ∇ • (u i e i ). (14)
With this system, the governing equations of p and u are identical to the original ones, but with additional non-differential terms. The original equations are recovered if σ is equal to zero. This formulation has four supplementary differential equations: two with the differential operator of equation ( 14) and two ordinary differential equations. The latter involve the spatial variation of each principal curvature (κ 2 and κ 3 ) in the other principal direction (respectively, e 3 and e 2 ). For a sphere, the curvature is constant and the fields r 1 and r 2 are equal to zero. In two dimensions, only one supplementary partial differential equation is required.

We now consider system [START_REF] Bécache | On the analysis of Bérenger's perfectly matched layers for Maxwell's equations[END_REF], which can be rewritten as

           -ıω p + ρc 2 1 ∏ k h k ( ∑ i ∂ ∂ξ i (( ∏ k̸ =i h k ) ûi ) ) = -ıω ( 1 - ∏ k hk ∏ k h k ) p, -ıωû i + 1 ρ 1 h i ∂ p ∂ξ i = -ıω ( 1 - hi h i ∏ k̸ =i h k ∏ k̸ =i hk ) ûi .
For the sake of clarity, the superscript ⋆ has been removed. Using equation ( 12), the right-hand sides of these equations become

-ıω ( 1 - ∏ k hk ∏ k h k ) p = -ıω ( 1 - ( 1 - σ 1 ıω ) ( 1 - σ 2 ıω ) ( 1 - σ 3 ıω
)) p = -( ∑ i σ i ) p - ( ∑ i ∏ k̸ =i σ i ) p (-ıω) -( ∏ i σ i ) p (-ıω) 2 , -ıω ( 1 - hi h i ∏ k̸ =i h k ∏ k̸ =i hk ) ûi = -ıω ( 1 - 1 -σ i /(ıω) ∏ k̸ =i (1 -σ k /(ıω))) ) ûi = -ıω ∏ k̸ =i (ıω -σ k ) -ıω(ıω -σ i ) ∏ k̸ =i (ıω -σ k ) ûi = -ıω ıω ( σ i - ∑ k̸ =i σ k ) + ∏ k̸ =i σ k ∏ k̸ =i (ıω -σ k ) ûi = - ( σ i - ∑ k̸ =i σ k ) (v i + ûi ) - ( ∏ k̸ =i σ k ) ŵi ,
where vi and ŵi have been introduced such that

vi + ûi = (ıω) 2 ∏ k̸ =i (ıω -σ k ) ûi = (ıω) 2 (ıω) 2 -ıω ( ∑ k̸ =i σ k ) + ( ∏ k̸ =i σ k ) ûi , ŵi = ıω ∏ k̸ =i (ıω -σ k ) ûi = ıω (ıω) 2 -ıω ( ∑ k̸ =i σ k ) + ( ∏ k̸ =i σ k ) ûi , which give
-ıωv i = - ( ∑ k̸ =i σ k ) (v i + ûi ) + ( ∏ k̸ =i σ k ) ŵi , -ıω ŵi = -(v i + ûi ).
Defining supplementary fields for p/(-ıω) and p/(-ıω) 2 , and using an inverse Fourier transform in time, we finally obtain the system

                           ∂p ∂t + ρc 2 1 ∏ k h k ( ∑ i ∂ ∂ξ i (( ∏ k̸ =i h k ) u i ) ) = -( ∑ i σ i ) p - ( ∑ i ∏ k̸ =i σ i ) p (1) -( ∏ i σ i ) p (2) , ∂u i ∂t + 1 ρ 1 h i ∂p ∂ξ i = - ( σ i - ∑ k̸ =i σ k ) (u i + v i ) - ( ∏ k̸ =i σ k ) w i , ∂p (1) ∂t = p, ∂v i ∂t = - ( ∑ k̸ =i σ k ) (u i + v i ) + ( ∏ k̸ =i σ k ) w i ,
∂p (2) ∂t = p (1) ,

∂w i ∂t = -(v i + u i ).
The tensor form of this system reads

                       ∂p ∂t + ρc 2 ∇ • u = -(σ 1 + σ 2 + σ 3 ) p -(σ 1 σ 2 + σ 1 σ 3 + σ 2 σ 3 ) p (1) -(σ 1 σ 2 σ 3 ) p (2) , ∂u ∂t + 1 ρ ∇p = -A(u + v) -Bw, ∂p (1) ∂t = p, ∂v ∂t = -C(u + v) + Bw, ∂p (2) ∂t = p (1) , ∂w ∂t = -(u + v), (15)
with the second-order tensors

A = (σ 1 -σ 2 -σ 3 ) I 1 + (-σ 1 + σ 2 -σ 3 ) I 2 + (-σ 1 -σ 2 + σ 3 ) I 3 , B = (σ 2 σ 3 ) I 1 + (σ 1 σ 3 ) I 2 + (σ 1 σ 2 ) I 3 , C = (σ 2 + σ 3 ) I 1 + (σ 1 + σ 3 ) I 2 + (σ 1 + σ 2 ) I 3 ,
where I i = e i ⊗ e i . As in system [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF], p and u are governed by the original equations where source terms are added. By contrast, this formulation involves two scalar and two vectorial supplementary fields, which are governed by ordinary differential equations. Only one scalar field and one vector field are required in two dimensions.

In the remainder, systems [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] and [START_REF] Bérenger | Perfectly Matched Layer (PML) for Computational Electromagnetics[END_REF] are respectively called the PML-PDE system and the PML-ODE system.

Interpretation and extension for domains with corners

Though the PML systems have been derived for convex domains with regular boundary, they can straightforwardly be adapted for squared and cuboidal domains with edge and corner regions where layers overlap. Formula ( 12) is classically used to stretch scale factors in the Cartesian directions to derive PML for squared and cuboidal domains (see e.g. [START_REF] Teixeira | Unified analysis of perfectly matched layers using differential forms[END_REF][START_REF] Teixeira | Complex space approach to perfectly matched layers: a review and some new developments[END_REF]). Therefore, the PML systems ( 13) and ( 15) can be used as is for Cartesian PMLs, but each absorption function σ i (x i ) is chosen independently and corresponds to the damping of waves in the corresponding Cartesian direction e i .

An analogy with the Cartesian PMLs provides a nice interpretation for the absorbing functions σ 2 and σ 3 expressed in equations ( 10) and [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]. They can be considered as absorption functions associated with the principal directions, but they are due to the curvature of the domain boundary. They are equal to zero if the boundary is planar.

The validity of the PML systems ( 13) and [START_REF] Bérenger | Perfectly Matched Layer (PML) for Computational Electromagnetics[END_REF] can also be extended to other convex domains having corners with right angles. For instance, they can deal with cylindrical domains having convex cross-sections. Let us consider a cylinder with the axis e 3 . In the lateral PML that surrounds the cylinder, the PML systems are then simply obtained by considering that κ 2 is the curvature of the cross section and κ 3 → 0 in relations ( 10) and ( 11); e 1 remains the radial direction and e 2 is tangent the lateral surface and perpendicular to e 3 . At the two extremities of the cylinder, the PMLs are planar and only σ 3 is non-zero. At the corner, σ 1 and σ 2 are defined as in the lateral PML, and σ 3 is defined as in the PML of the border.

Classification and mathematical properties

The derived PML-PDE system (13) and PML-ODE system (15) can be categorized into two well-known families of PMLs. The first family contains the Bérenger-like PML systems that are built using the splittingfield technique of Bérenger [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] or the complex coordinate stretch technique [START_REF] Chew | A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates[END_REF][START_REF] Chew | Complex coordinate stretching as a generalized absorbing boundary condition[END_REF]. These systems involve additional fields governed by differential equations with spatial partial derivatives (see e.g. [START_REF] Bérenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Collino | The perfectly matched layer in curvilinear coordinates[END_REF][START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Hu | On absorbing boundary conditions for linearized euler equations by a perfectly matched layer[END_REF][START_REF] Hu | Development of PML absorbing boundary conditions for computational aeroacoustics: A progress review[END_REF][START_REF] Hu | Absorbing boundary conditions for nonlinear Euler and Navier-Stokes equations based on the perfectly matched layer technique[END_REF][START_REF] Navon | A perfectly matched layer approach to the linearized shallow water equations models[END_REF]). The second family corresponds to PML systems based on frequency-dependent complex material properties [START_REF] Gedney | An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices[END_REF][START_REF] Sacks | A perfectly matched anisotropic absorber for use as an absorbing boundary condition[END_REF][START_REF] Zhao | GT-PML: Generalized theory of perfectly matched layers and its application to the reflectionless truncation of finite-difference time-domain grids[END_REF] and written in the time-domain with supplementary ordinary differential equations [START_REF] Abarbanel | On the construction and analysis of absorbing layers in CEM[END_REF][START_REF] Chen | Long-time stability and convergence of the uniaxial perfectly matched layer method for time-domain acoustic scattering problems[END_REF][START_REF] Rahmouni | An algebraic method to develop well-posed PML models: Absorbing layers, perfectly matched layers, linearized Euler equations[END_REF]. Such PMLs are generally referred as uniaxial PMLs.

The mathematical properties of PMLs from both families have been studied in many works and are still an active field of research (see e.g. [START_REF] Abarbanel | A mathematical analysis of the PML method[END_REF][START_REF] Abarbanel | On the construction and analysis of absorbing layers in CEM[END_REF][START_REF] Bécache | On the analysis of Bérenger's perfectly matched layers for Maxwell's equations[END_REF][START_REF] Bécache | Remarks on the stability of cartesian PMLs in corners[END_REF][START_REF] Diaz | A time domain analysis of PML models in acoustics[END_REF][START_REF] Duru | The role of numerical boundary procedures in the stability of perfectly matched layers[END_REF][START_REF] Halpern | Hyperbolic boundary value problems with trihedral corners[END_REF][START_REF] Halpern | The analysis of matched layers[END_REF][START_REF] Kaltenbacher | A modified and stable version of a perfectly matched layer technique for the 3-d second order wave equation in time domain with an application to aeroacoustics[END_REF][START_REF] Rahmouni | An algebraic method to develop well-posed PML models: Absorbing layers, perfectly matched layers, linearized Euler equations[END_REF]). For wave problems described with symmetric hyperbolic systems (e.g. acoustic system (1) and Maxwell's equations), the PML systems derived with the original strategy of Bérenger are weakly hyperbolic, and lead to weakly well-posed problems [START_REF] Abarbanel | A mathematical analysis of the PML method[END_REF][START_REF] Bécache | On the analysis of Bérenger's perfectly matched layers for Maxwell's equations[END_REF]. By contrast, the uniaxial PML systems are symmetric hyperbolic and lead to strongly well-posed problems [START_REF] Abarbanel | On the construction and analysis of absorbing layers in CEM[END_REF][START_REF] Halpern | The analysis of matched layers[END_REF][START_REF] Rahmouni | An algebraic method to develop well-posed PML models: Absorbing layers, perfectly matched layers, linearized Euler equations[END_REF]. To our knowledge, only proofs of weak stability have been proposed for standard PML systems from both families [START_REF] Bécache | On the analysis of Bérenger's perfectly matched layers for Maxwell's equations[END_REF][START_REF] Bécache | On the long-time behavior of unsplit perfectly matched layers[END_REF][START_REF] Chen | Long-time stability and convergence of the uniaxial perfectly matched layer method for time-domain acoustic scattering problems[END_REF]. This means that no exponentially growing modes are supported by the equations, but linear modes could pollute the solution in long-duration simulations. Since the PML systems studied in [START_REF] Bécache | On the analysis of Bérenger's perfectly matched layers for Maxwell's equations[END_REF][START_REF] Bécache | On the long-time behavior of unsplit perfectly matched layers[END_REF][START_REF] Chen | Long-time stability and convergence of the uniaxial perfectly matched layer method for time-domain acoustic scattering problems[END_REF] correspond to planar versions of systems ( 13) and ( 15), the PML formulations proposed in this work will also suffer from the same long-time instability.

Let us note that modified PML formulations have been proposed to avoid any linear growing mode, and to improve the accuracy for problems involving evanescent and grazing waves. For instance, the so-called complex frequency shifted PML formulations [START_REF] Kuzuoglu | Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers[END_REF][START_REF] Roden | Convolutional PML (CPML): An efficient FDTD implementation of the CFS-PML for arbitrary media[END_REF] are derived using a modified coordinate stretch, where new parameters are introduced. In the time domain, these formulations can be seen as generalizations of standard PML formulations with supplementary terms and equations. The formulations proposed in this work could be enriched using the modified coordinate stretch, but at the cost of larger PML systems and more expensive computations.

Numerical simulations

Both time-dependent PML formulations are tested by means of a three-dimensional benchmark where a spherical wavefront is propagated in a truncated domain. The computational domain is the half of an ellipsoid of revolution, whose shell is extended with a PML. The simulations are performed using discontinuous finite element schemes that are described in section 3.1. After the description of the benchmark in section 3.2, accuracy results are proposed to validate our implementation and performance results are discussed (section 3.3). An illustration of application is finally proposed in section 3.4.

For the PML-PDE formulation [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF], both fields r i are set to zero and the corresponding equations are not computed. Since the domain has an axis of revolution, one of the principal curvatures is constant, and the corresponding field is equal zero. Setting the other field to zero is an approximation that corresponds to neglecting the effect of the spatial variation of the other principal curvature. With our specific benchmark, we have observed a posteriori that the performances of the formulation remain very good in comparison to the PML-ODE formulation, which is implemented without any approximation.

Discontinuous finite element schemes

We propose numerical schemes based on a nodal discontinuous Galerkin (DG) method for the PML systems. The schemes are derived following a classical strategy for hyperbolic systems [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: algorithms, analysis, and applications[END_REF], with some further developments to discretize the supplementary equations.

The unknown fields are approximated on a spatial mesh of Ω made of non-overlapping polyhedral cells, Ω = ∪ k=1,...,K D k , where K is the number of cells and D k is the k th cell. With the nodal DG method, each scalar field and each Cartesian component of vector fields is approximated by piecewise polynomial functions, where discontinuities correspond to interfaces between elements. The discrete unknowns correspond to the values of fields at the nodes of elements. Over each cell D k , the approximation of any field q(x, t) can then be written as

q k (x, t) = N k ∑ n=1 q k,n (t) ℓ k,n (x), ∀x ∈ D k ,
where N k is the number of nodes over element k, q k,n (t) is the value of the field at node n, and ℓ k,n (x) is the multivariate Lagrange polynomial function associated with node n.

For both the original system (1) and the PML systems ( 13) and ( 15), the governing equations of p and u are discretized in the same manner. Multiplying these equations by test functions p and û, integrating the resulting equations over any cell D k , and using integration by parts, we get the weak form

( ∂p ∂t , p) D k - ( ρc 2 u, ∇p ) D k + ⟨f p , p⟩ ∂D k = (s p , p) D k , ( ∂u ∂t , û) D k - ( p ρ , ∇ • û) D k + ⟨f u , û⟩ ∂D k = (s u , û) D k ,
with the fluxes

f p = ρc 2 n k • u, f u = p ρ n k ,
where n k is the outward unit normal on the boundary ∂D k . For both PML systems, s p and s u contain the right-hand side terms of the governing equations of p and u. The discrete fields being discontinuous at the interface between elements, the fluxes are replaced by numerical fluxes which involve values of fields from both sides of interfaces. We choose the upwind fluxes

f num p = ρc 2 n k • { u } -c p , ( 16
)
f num u = 1 ρ n k { p } -c n k ( n k • u ) , ( 17
)
where the mean { q } and the semi-jump q of any field q are defined as

{ q } = q + + q - 2 , q = q + -q - 2 .
The plus and minus subscripts denote values on the side of the neighboring and current cells, respectively. These fluxes are classically obtained by using a Riemann solver (see e.g. [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: algorithms, analysis, and applications[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]). At the external boundary, we use a first-order absorbing boundary condition (ABC) based on characteristics, which is simply implemented by taking the external values p + and u + equal to zero in the numerical fluxes ( 16)- [START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF].

For the PML-PDE system (13), the partial differential equations that govern the supplementary fields p 1 and p 2 are discretized by analogy with the governing equation of p. Multiplying these equations by test functions, integrating over an element D k of the layer Σ, and integrating by parts lead to the weak form

( ∂p 1 ∂t , p1 ) D k - ( ρc 2 e 1 (e 1 • u), ∇p 1 ) D k + ⟨f p1 , p1 ⟩ ∂D k = (s p1 , p1 ) D k , ( ∂p 2 ∂t , p2 ) D k - ( ρc 2 e 2 (e 2 • u), ∇p 2 ) D k + ⟨f p2 , p2 ⟩ ∂D k = (s p2 , p2 ) D k ,
where f p1 and f p2 are defined as and s p1 and s p2 contain the source terms written in the right-hand side of equations. Since system ( 13) is no longer strictly hyperbolic, we cannot use the Riemann solver to derive upwind fluxes. Noting the similitudes between f p and both f p1 and f p2 , we propose to use the Lax-Friedrichs fluxes

f p1 = ρc 2 (n k • e 1 )(e 1 • u), f p2 = ρc 2 (n k • e 2 )(e 2 • u),
f num p1 = ρc 2 (n k • e 1 )(e 1 • { u } ) -c p 1 , f num p2 = ρc 2 (n k • e 2 )(e 2 • { u } ) -c p 2 .
These fluxes are valid inside the layer Σ, but not at its external boundary Γ and at the interface Υ with the domain Ω. Indeed, because the supplementary fields are defined only inside the layer, we cannot use their semi-jumps p 1 and p 2 on Υ and Γ. Fortunately, we note that n k = -e 1 on Υ and n k = e 1 on Γ. We also have n k • e 2 = 0 on both surfaces. Therefore, we have f p1 = f p and f p2 = 0 in both cases. We then use the numerical flux f num p instead of f num p1 , and f num p2 is simply set to zero.

Since the supplementary equations of the PML-ODE system (15) do not involve spatial differential operators, the derivation of the scheme is straightforward. In addition to the above scheme for the governing equations of p and u, ordinary differential equations will be only local, at each node of each element, for the supplementary fields p (1) , p (2) , v and w. We use the same time-stepping scheme for both partial and ordinary differential equations.

Description of the benchmark

The computational domain Ω is the half of an ellipsoid, with axes of lengths 330 m (x-direction) and 120 m (y-and z-directions), respectively. A symmetry boundary condition is used in the section of the semiellipsoid, and the shell is extended with a PML Σ (Figure 3). To generate the wavefront, a Gaussian pulse is prescribed as initial condition on the pressure, p(x, 0) = e -r 2 /R 2 with r = ∥x -x 0 ∥, x 0 = (-122.5 m, 0 m, 20 m) and R = 7.5 m, while the other fields are initially set to zero. With this setting, the initial pressure is negligible inside the PML. We use the physical parameters c = 1.5 km/s and ρ = 1 kg/m 3 .

The mesh is made of tetrahedra in the truncated domain Ω and prisms in the layer Σ, with element sizes between 3 m and 5 m. The mesh of the layer is generated by extrusion. The main reason is that it can be automatically done for truncated domain of any shape with a regular enough boundary Υ. It does not require the analytic representation of the external boundary Γ of the layer. With an ellipsoidal domain, the external boundary Γ is not an ellipsoid (see e.g. [START_REF] Teixeira | Complex space approach to perfectly matched layers: a review and some new developments[END_REF]), and, to the best of our knowledge, no analytical representation is known. Ideally, a mesh with curvilinear elements should be used in order to exactly follow the interface Υ and the boundary Γ that are curved. Elements with planar faces could be used with a lower computational cost, but the approximate polynomial representations of the surfaces introduce a geometric error. To reduce the global error, it is not clear if using curvilinear elements is computationally more efficient than using elements with planar faces and a thicker layer. In this paper, we consider only elements with planar faces.

We consider different thicknesses for the layer (1, 3 and 5 mesh cells) and both first-degree and seconddegree basis functions. All the meshes are generated with Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. The time-stepping is made with the fourth-order Runge-Kutta scheme. The time step ∆t is equal to 0.05 ms and 0.025 ms for first-degree and second-degree basis functions, respectively. Two alternative absorption functions σ(ξ 1 ) are considered: the hyperbolic function σ h and the shifted hyperbolic function σ sh , defined as

σ h (ξ 1 ) = c δ -ξ 1 , σ sh (ξ 1 ) = c δ ξ 1 δ -ξ 1 .
It has been shown in previous works [START_REF] Bermúdez | An exact bounded perfectly matched layer for time-harmonic scattering problems[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF] Modave | An optimum PML for scattering problems in the time domain[END_REF][START_REF] Modave | Optimizing perfectly matched layers in discrete contexts[END_REF] that both these functions do not require the tuning of parameters. In particular, we have shown [START_REF] Modave | Optimizing perfectly matched layers in discrete contexts[END_REF] that σ sh is a good choice for a discontinuous Galerkin implementation with a Cartesian PML formulation similar to the PML-PDE system [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. Here, both functions are tested with both PML systems ( 13) and [START_REF] Bérenger | Perfectly Matched Layer (PML) for Computational Electromagnetics[END_REF].

The PML formulations require to know the coordinate ξ 1 , the principal curvatures κ 2 and κ 3 , and the local orthonormal frame (e 1 , e 2 , e 3 ). Both ξ 1 and e 1 are computed numerically by finding, in the initialization of the simulation, the closest node belonging to the mesh of the surface Υ. To define the principal curvatures and the other vectors, we use the fact that the domain is an ellipsoid of revolution. Then, e 2 belongs to the plane tangent to Υ and is oriented towards one pole, and e 3 is perpendicular to both e 2 and e 3 . κ 2 corresponds to the curvature of an ellipse [START_REF] Weisstein | Ellipse[END_REF],

κ 2 = a x a y [p 2
x (a y /a x ) 2 + (p 2 y + p 2 z )(a x /a y ) 2 ] 3/2 , and κ 2 is obtained from the Gaussian curvature G of an ellipsoid [START_REF] Weisstein | Ellipsoid[END_REF],

G = κ 2 κ 3 = (a x a y a z ) 2 [(p x a y a z /a x ) 2 + (p y a x a z /a y ) 2 + (p z a x a y /a z ) 2 ] 2 .
In these formulas, a x , a y and a z are the semi-axes of the ellipsoid and p x , p y and p z are the Cartesian coordinates of p(ξ 2 , ξ 3 ). Let us note that, in our case, a y = a z .

Validation and comparison

Figure 4 shows snapshots of the pressure p at different instants with a three-cells layer and the PML-PDE formulation. A spherical wave is generated by the initial pulse. As time goes by, the spherical wavefront reaches the PML with different angles of incidence, and propagates along the interface at a nearly grazing incidence. We observe that the spherical wave is not deformed near the interface, as we wished for, and the pressure is damped inside the layer.

The simulation is performed with different settings over the period [0 ms, 500 ms]. This includes a first phase where the wavefront is traveling in the truncated domain (before t = 200 ms) and a second phase where the wavefront has left the domain (after t = 200 ms). For all the settings, the difference between the numerical solution and the exact solution is quantified using the relative error defined as where the exact solution is given by

Error(t) = ∫ Ω ( 1 2ρc 2 (p exact (t, x) -p num (t, x)) 2 + ρ 2 ∥u exact (t, x) -u num (t, x)∥ 2 ) dx ∫ Ω ( 1 2ρc 2 (p exact (0, x)) 2 + ρ 2 ∥u exact (0, x)∥ 2 ) dx
p exact (t, x) = 1 2 ( r -ct r e -(r-ct) 2 /R 2 + r + ct r e -(r+ct) 2 /R 2 ) , u exact (t, x) = 1 2cρ [( R 2 2r 2 + r -ct r ) e -(r-ct) 2 /R 2 - ( R 2 2r 2 + r + ct r ) e -(r+ct) 2 /R 2 ] x -x 0 r .
The error corresponds to the ratio of the energy associated with the error fields to the energy associated with the initial fields in the truncated domain Ω.

Figure 5 shows the time-evolution of the error when the absorbing boundary condition (ABC) and threecells PMLs are used as truncation technique. Both first-degree and second-degree basis functions are tested. During the period [0 ms, 200 ms], all the PMLs exhibit a similar error for each polynomial degree. This error dramatically decreases when increasing the degree from one to two. By contrast, the error corresponding to the ABC is larger and stays at the same level (between 10 -2 and 10 -1 ) for both polynomial degrees. These observations are coherent with the properties of these truncation techniques. The considered ABC is only a first-order approximation that cannot simulate oblique outgoing waves without any reflection. The observed error then is dominated by the modeling error associated with spurious reflected waves already present at the continuous level. For the PMLs, the significant decay of error indicates that numerical errors dominate.

In the second phase of the simulation (after t = 200 ms), the ABC error curves are decreasing and similar with both polynomial degrees. The error is dominated by the modeling error corresponding to the multiple reflections of the initial wavefront that are leaving the domain. With the PMLs, the error curves are in a quite large range (between 10 -7 and 10 -5 ) and have slightly different behaviors, which confirms that the kind of PML formulation and the absorption function impact the error. For both degrees, the error is far higher with the ABC than with the PMLs at the beginning of the second phase, but the ABC error goes under the error obtained with some PMLs at the end of the simulation. While the ABC error is dramatically decreasing during the second phase, the PML errors are only slightly decreasing or even slightly increasing. This phenomenon could be explained by the linear growing modes that can be generated inside the PML and that can pollute the long-duration solution (see the end of section 3.3).

In order to compare the PML implementations, we now consider the mean value of the relative error over the period [250 ms, 500 ms], where the PMLs exhibit different errors. Figure 6 shows the mean error as function of the total runtime when using 2 Dual-Core Intel Xeon Processors L5420. All the one-cell PMLs provide similar errors for similar computational costs. For three-cells and five-cells layers, the PML-PDE implementation is more accurate with the shifted hyperbolic absorption function σ sh than with the hyperbolic one σ h . This result is coherent with the conclusions of previous works with similar formulations [START_REF] Modave | An optimum PML for scattering problems in the time domain[END_REF][START_REF] Modave | Optimizing perfectly matched layers in discrete contexts[END_REF]. For the PML-ODE implementation, both functions provide quite similar results. σ h is slightly better than σ sh with a three-cells layer, while the converse is true with a five-cells layer. Finally, the PML-PDE implementation (with the shifted hyperbolic function) is more efficient than the PML-ODE. For a five-cells layer, it even provides a more accurate result with a smaller runtime.

Illustration of application

Our strategy is now illustrated by the scattering of a submarine placed in the center of the ellipsoidal domain. The submarine is approximately 120 m in length, and the dimensions of the ellipsoid are not changed. The PML-PDE is used together with the shifted hyperbolic absorption function σ sh and the first-order ABC to terminate the layer. The mesh is made of 183,707 tetrahedra (truncated domain) and 62,650 prisms (five-cells layer). First-degree basis functions are used. There are 7,450,112 degrees of freedom.

Figure 7 shows snapshots of the pressure p at different instants. The spherical wave generated by the initial pulse hits the front of the submarine, and creates perturbations in the pressure field. The wavefront propagates along the submarine and is nearly grazing at the boundary of the domain. The perturbations, as well as the primary spherical wave, have correctly left the front zone at t = 200 ms. These results have been obtained using 8 Dual-Core Intel Xeon Processors L5420. The total runtime was 1 h 40 min, corresponding to the period [0 ms, 200 ms] and to 8000 time steps.

Conclusion

This paper is dedicated to generalizing PML formulations for acoustic wave propagation with generallyshaped convex domains, which offers flexibility when choosing the shape of the computational domain. Our strategy is based on the complex stretch of a specific curvilinear coordinate in the time-harmonic equations, and on the use of supplementary differential equations to write the time-domain PML formulations.

Two time-dependent PML formulations have been derived for the pressure-velocity system. One formulation involves supplementary PDEs, while only ODEs are required for the other. Both have been implemented in a discontinuous Galerkin finite-element solver and tested with a three-dimensional benchmark. The best accuracy is observed with the PML-PDE formulation, despite an approximation made during the implementation. That formulation is faster than the other, but the complete version requires the knowledge of supplementary geometrical information. In a more applicative context, we believe that all the required The approach is quite general and offers a wide range of extensions. Since the PML equations are written at the continuous level, they can be implemented with other numerical methods. As explained at the end of section 2.3, the formulations can be adapted to deal with domains having corners, though they have been derived assuming that the domain border was regular enough. The formulations can also be improved using a complex coordinate stretch with a frequency shift [START_REF] Kuzuoglu | Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers[END_REF][START_REF] Roden | Convolutional PML (CPML): An efficient FDTD implementation of the CFS-PML for arbitrary media[END_REF].

Figure 1 :

 1 Figure 1: Curvilinear coordinates and local frame associated with the boundary Υ in two dimensions. The curves of iso-value coordinates are represented in Fig. (a). Gray curves are parallel. Red lines are straight and perpendicular to Υ. Fig. (b) shows the local frame and the radial coordinate ξ1.

Figure 2 :

 2 Figure 2: Complex coordinate stretch of the radial coordinate ξ1.

Figure 3 :

 3 Figure 3: Geometry and mesh for the three-dimensional benchmark with a three-cells PML. This mesh is made of 58,396 tetrahedra in the truncated domain and 3 × 6, 998 prisms in the layer.

Figure 4 :

 4 Figure 4: Snapshots of the solution of the three-dimensional benchmark at different instants. Colored surfaces are iso-surfaces of p(x, t). The inner and outer white surfaces correspond to the interior and exterior boundaries of the PML, respectively.

  with σh PML-PDE with σsh PML-ODE with σh PML-ODE with σsh with σh PML-PDE with σsh PML-ODE with σh PML-ODE with σsh

Figure 5 :

 5 Figure 5: Time-evolution of the error inside the domain Ω with the absorbing boundary condition (ABC) and threecells PMLs. First-degree (a) and second-degree (b) basis functions are considered with both PML formulations and both absorbing functions.

  with σh PML-PDE with σsh PML-ODE with σh PML-ODE with σsh with σh PML-PDE σsh PML-ODE with σh PML-ODE with σsh

Figure 6 :

 6 Figure 6: Mean error over the period [250 ms, 500 ms] as a function of the runtime. The three points of each curve correspond to one-cell (left), three-cells (middle) and five-cells (right) layers.

Figure 7 :

 7 Figure 7: Snapshots of the solution of the three-dimensional benchmark at different instants. Colored surfaces are iso-surfaces of p(x, t). The inner and outer white surfaces correspond to the interior and exterior boundaries of the PML, respectively

Acknowledgements

This work was supported in part through the ARC grant for Concerted Research Actions (ARC WAVES 15/19-03), financed by the Wallonia-Brussels Federation, and by the Belgian Science Policy Office under grant IAP P7/02. Axel Modave was partially supported by an excellence grant from Wallonie-Bruxelles International (WBI). He was a Postdoctoral Researcher on leave with the F.R.S-FNRS. Jonathan Lambrechts was a Postdoctoral Researcher with the F.R.S-FNRS.