Bruno Lombard
email: lombard@lma.cnrs-mrs.fr

Jean-François Mercier
email: jean-francois.mercier@ensta.fr

Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz resonators

Keywords: nonlinear acoustics, solitons, Burgers equation, fractional derivatives, diffusive representation, time splitting, shock-capturing schemes

Acoustic wave propagation in a one-dimensional waveguide connected with Helmholtz resonators is studied numerically. Finite amplitude waves and viscous boundary layers are considered. The model consists of two coupled evolution equations: a nonlinear PDE describing nonlinear acoustic waves, and a linear ODE describing the oscillations in the Helmholtz resonators. The thermal and viscous losses in the tube and in the necks of the resonators are modeled by fractional derivatives. A diffusive representation is followed: the convolution kernels are replaced by a finite number of memory variables that satisfy local ordinary differential equations. A splitting method is then applied to the evolution equations: their propagative part is solved using a standard TVD scheme for hyperbolic equations, whereas their diffusive part is solved exactly. Various strategies are examined to compute the coefficients of the diffusive representation; finally, an optimization method is preferred to the usual quadrature rules. The numerical model is validated by comparisons with exact solutions. The properties of the full nonlinear solutions are investigated numerically. In particular, the existence of acoustic solitary waves is confirmed.

Introduction

Propagation of linear acoustic waves in lattices has already been the subject of a large body of theoretical and experimental works. Floquet-Bloch band gaps are present in ordered lattices [START_REF] Brillouin | Propagation des Ondes dans les Milieux Périodiques[END_REF], whereas localization occurs in disordered cases [START_REF] Richoux | Propagation of mechanical waves in a one-dimensional nonlinear disordered lattice[END_REF]. Nonlinearities, when they are present, are usually considered at discrete points [START_REF] Richoux | Acoustic wave dispersion in a one-dimensional lattice of nonlinear resonant scatterers[END_REF].

The propagation of nonlinear acoustic waves in lattices was addressed by Sugimoto and his coauthors in a series of theoretical and experimental studies, whose original purpose was the reduction of shock waves generated by high-speed train into a tunnel [START_REF] Sugimoto | Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves[END_REF][START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF][START_REF] Sugimoto | Acoustic solitary waves in a tunnel with an array of Helmholtz resonators[END_REF][START_REF] Sugimoto | Experimental demonstration of generation and propagation of acoustic solitary waves in a air-filled tube[END_REF][START_REF] Sugimoto | Verification of acoustic solitary waves[END_REF]. The configuration under study was made up of a tube connected with Helmholtz resonators (figure 1). These resonators induce dispersion that competes with the nonlinear effects and may prevent from the emergence of shocks.

More fundamental questions are also raised, concerning well-known nonlinear waves called solitons [START_REF] Toda | Nonlinear waves and solitons[END_REF][START_REF] Tao | Why are solitons stable ?[END_REF]. Those are stable solitary waves that maintain their shape while traveling at constant speed. Solitons are caused by the cancellation of nonlinear and dispersive effects in the medium. Many models have soliton solutions, for instance the Korteweg-de Vries equation, the nonlinear Schrödinger equation, and the sine-Gordon equation. In acoustics, dissipation is largely greater than dispersion, so that it was commonly thought that it was impossible to generate those solitary waves. Thanks to the Sugimoto's configuration, it was shown that acoustic solitons can exist and propagate in place of shock waves [START_REF] Sugimoto | Verification of acoustic solitary waves[END_REF].

The model proposed by Sugimoto involves two coupled equations: a nonlinear Partial Differential Equation (PDE) PDE describing the propagation of finite amplitude acoustic waves in the tube, and a linear Ordinary Differential Equation (ODE) describing the oscillations in the Helmholtz resonators. The dissipative effects in the tube and in the necks of the resonators are modeled by fractional derivatives [START_REF] Matignon | An introduction to fractional calculus[END_REF], that introduce convolution products. A good numerical modeling relies on the following three specifications:

• accurate computation of nonlinear non-smooth waves;

• efficient computation of fractional derivatives, without storing the previous values of the solutions;

• stable algorithm under a Courant-Friedrichs-Lewy (CFL) condition, whatever the physical parameters and the amplitude of the waves are. The first specification, based for example on shock-capturing schemes, has been well-known for decades [START_REF] Leveque | Numerical methods for conservation laws[END_REF]. The second specification is much less standard. A diffusive representation of fractional operators is used instead of a discretization of convolution operators [START_REF] Yuan | A numerical scheme for dynamic systems containing fractional derivatives[END_REF][START_REF] Lu | Wave field simulation for heterogeneous porous media with singular memory drag force[END_REF][START_REF] Galucio | An adaptation of the Gear scheme for fractional derivatives[END_REF][START_REF] Diethelm | An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives[END_REF][START_REF] Deü | Simulation of fractionally damped mechanical systems by means of a Newmark-diffusive scheme[END_REF][START_REF] Birk | An improved non-classical method for the solution of fractional differential equations[END_REF]. In doing so, the fractional derivatives are replaced by a set of memory variables that satisfy local-in-time linear differential equations. Determining the quadrature coefficients of the diffusive representation is a crucial issue for the accuracy and the efficiency of the method. Various strategies, usually based on orthogonal polynomials, have been proposed in the literature. We propose here a more efficient strategy, where the coefficients are optimized with respect to the dispersion relation of Sugimoto's model. Lastly, the stability specification requires an adequate coupling between the PDE and the ODE. A naive coupling between these equations usually ensures an increase of discrete energy.

On the contrary, we obtain here a stable scheme under the optimal CFL condition.

The paper is organized as follows. The Sugimoto's model is presented in section 2. Dispersion analysis in the linear case is developed, and degeneracy towards Korteweg-de Vries equations is specified. The diffusive representation of fractional derivatives, leading to a first-order system of PDE, is described in section 3. The numerical methods are detailed in section 4: a splitting procedure to ensure an optimal CFL condition and to take advantage of efficient methods; a Total Variation Diminishing (TVD) scheme for the advection-Burgers PDE; and an exact integration of the diffusive part. The determination of the weights and nodes of the diffusive representation is discussed in section 5. Numerical experiments are proposed in section 6. Four tests are presented, concerning successively nonlinear acoustic waves in the tube, oscillations in the resonators, and the coupling (linear and nonlinear) between the two subsystems. The numerical simulations reveal acoustic solitary waves. A conclusion is drawn and future directions of research are outlined in section 7. The configuration under study is made up of an air-filled tube connected with Helmholtz resonators (figure 2). The cylindrical resonators are uniformly distributed along the tube. The geometrical parameters are the radius of the guide R; the axial spacing between resonators D; the radius of the neck r; the length of the neck L; the radius of the cavity r h ; and the height of the cavity H, which may vary depending on the resonator. Hence, the cross-sectional area of the guide, A, is π R 2 and that of the neck, B, π r 2 , the volume of each resonator, V , is π r 2 h H, and the reduced radius is:

Physical modeling

R * = R 1 - R 2 D B A = R 1 - r 2 2 D R . ( 1
)
The physical parameters are the ratio of specific heats at constant pressure and volume γ; the pressure at equilibrium p 0 ; the density at equilibrium ρ 0 ; the Prandtl number Pr; the kinematic viscosity ν; and the ratio of shear and bulk viscosities µ v /µ; hence the linear sound speed, a 0 , the sound diffusivity, ν d , the dissipation in the boundary layer, C, the natural angular frequency of the resonator, ω 0 and the natural angular frequency of the tube coupled with the resonator, ω 1 , are:

a 0 = γ p 0 ρ 0 , ν d = ν 4 3 + µ v µ + γ -1 Pr , C = 1 + γ -1 √ Pr , ω 0 = a 0 B L V = a 0 r r h 1 √ L H , ω 1 = ω 0 1 + V 2 A D . (2)
Under a one-dimensional assumption, the variables are the axial velocity of the gas u and the excess pressure in the cavity p (denoted p ′ 2 = p 2p 0 in the papers by Sugimoto, p 2 being the pressure in the cavity of the resonators). The wavelength of the initial disturbance is λ; hence the acoustic Mach number, M, the parameter of nonlinearity, ε, the characteristic angular frequency, ω, the excess pressure in the tube (in the linear theory), p ′ , and the intensity of sound (in dB), I, are:

M = u a 0 , ε = γ + 1 2 M, ω = 2 π a 0 λ , p ′ p 0 = γ u a 0 , I = 20 log 2 γ γ + 1 p 0 p ref ε , (3)
where p ref = 2 10 -5 Pa.

Governing equations

The main assumptions underlying Sugimoto's model are [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF]:

• low-frequency (ω < ω * = 1.84 a 0 R ), so that only the plane mode propagates and the 1D approximation is valid [START_REF] Chaigne | Acoustique des Instruments de Musique[END_REF];

• weak acoustic nonlinearity in the tube (ε ≪ 1) [START_REF] Hamilton | Nonlinear Acoustics[END_REF];

• continuous distribution of resonators (λ ≫ D); • no interactions between neighboring resonators ( V A D = r h R 2 H D ≪ 1);
• linear response of the resonators, no turbulence.

Under these hypotheses, the right-going simple wave is modeled by the following coupled PDE-ODE system

         ∂u ∂t + ∂ ∂x au + b u 2 2 = c ∂ -1/2 ∂t -1/2 ∂u ∂x + d ∂ 2 u ∂x 2 -e ∂p ∂t , (4a)
∂ 2 p ∂t 2 + f ∂ 3/2 p ∂t 3/2 + g p = h u, (4b)
with the parameters

a = a 0 , b = γ + 1 2 , c = C a 0 √ ν R * , d = ν d 2 , e = V 2 ρ 0 a 0 A D , f = 2 √ ν r , g = ω 2 0 , h = ω 2 0 γ p 0 a 0 . (5)
PDE (4a) models nonlinear acoustic waves in the tube (coefficients a and b).

Viscous and thermal losses in the boundary layer of the tube are introduced by the coefficient c [START_REF] Chester | Resonant oscillations in closed tubes[END_REF]. The diffusivity of sound in the tube is also modeled by the coefficient d. ODE (4b) models the air oscillation in the neck of the resonators (coefficients f and g) [START_REF] Monkewitz | The response of Helmholtz resonators to external excitation. Part 1. Single resonators[END_REF][START_REF] Monkewitz | The response of Helmholtz resonators to external excitation. Part 2. Arrays of slit resonators[END_REF]. The coupling between the two equations is done by the coefficients e and h. If the resonators are suppressed (H → 0 and thus V → 0), then the coefficient e → 0: no coupling occurs, and the classical Chester's equation is recovered [START_REF] Menguy | Weakly non-linear gas oscillations in airfilled tubes ; solutions and experiments[END_REF].

Fractional operators of order -1/2 and 3/2 are involved in the system (4). These operators model the viscous and thermal losses in the tube and in the resonators. These losses are respectively proportional to 1/(i ω) 1/2 and (i ω) 3/2 in the frequency domain. In (4a), the Riemann-Liouville fractional integral of order 1/2 of a function w(t) is defined by

∂ -1/2 ∂t -1/2 w(t) = H(t) √ π t * w = 1 √ π t 0 (t -τ ) -1/2 w(τ ) dτ, ( 6
)
where * is the convolution product in time, and H(t) is the Heaviside step function [START_REF] Matignon | An introduction to fractional calculus[END_REF]. The fractional derivative of order 3/2 in (4b) is obtained by differentiating [START_REF] Chester | Resonant oscillations in closed tubes[END_REF] twice with respect to t. Figure 3: phase velocity in the linear regime, with resonators (V = 0 ⇔ e = 0). Inviscid (ν = 0) and viscous (ν = 0) case; see table 1 for the parameters. The vertical dotted lines denote ω 0 and ω 1 . The horizontal dotted lines denote υ < a 0 and a 0 (12).

Dispersion analysis

In this section we present the dispersion analysis of (4) in the linear case b = 0. The results obtained will be useful to adjust the terms in the diffusive representation of the fractional derivatives (section 5).

Let us define the Fourier transforms in time and space

w(ω) = +∞ -∞ w(t) e -i ω t dt, w(k) = +∞ -∞ w(x) e i kx dx, ( 7
)
where ω is the angular frequency and k is the wavenumber. Applying these transforms to (4) provides a system of two linear equations whose determinant must be null, which yields the dispersion relation between ω and k.

Defining the symbol of the half-order integral ( 6)

χ(ω) = 1 (i ω) 1/2 , (8)
and setting the coefficients

       D 2 (ω) = i d[g -ω 2 (1 + f χ)], D 1 (ω) = (a -c χ)[g -ω 2 (1 + f χ)], D 0 (ω) = ω[ω 2 (1 + f χ) -(g + e h)], (9)
the dispersion relation takes the form

D(k, ω) = D 2 (ω) k 2 + D 1 (ω) k + D 0 (ω) = 0. ( 10
)
Let us now describe this dispersion relation. First, we will discuss the general case in which the guide is coupled with resonators (e = 0). Neglecting the diffusivity of sound (d = 0), one obtains k = -D 0 (ω) / D 1 (ω). Otherwise, [START_REF] Flannery | Numerical Recipes in C: the Art of Scientific Computing[END_REF] has two roots k 1 and k 2 , and the root k with minimal modulus is selected, hence the phase velocity υ = ω / ℜe(k) and the attenuation -ℑm(k). In the inviscid case (c = f = 0 adding to d = 0), the explicit expressions are obtained

υ(ω) = a 0 ω 2 -ω 2 0 ω 2 -ω 2 1 and ℑm(k) = 0, (11)
where we use (2) and [START_REF] Chaigne | Acoustique des Instruments de Musique[END_REF]. The basic properties of (11) are deduced:

                 υ = υ(0) = a 0 1 + V 2 A D < a 0 , υ(ω 0 ) = 0, lim ω→ω ± 1 υ(ω) = ∓∞, lim ω→+∞ υ(ω) = a 0 , υ ′ (ω) = -2 a 0 ω 2 0 V 2 A D ω (ω 2 -ω 2 1 ) 2 < 0. ( 12
)
Hypothetically, V A D ≪ 1 and hence ω 1 is close to ω 0 (2). Far from ω 0 and ω 1 , the curve of υ is quite flat. In the viscous case (c = 0 and f = 0), the phase velocity does not vanish anymore at ω 0 . Two vertical asymptotes of υ occur near ω 0 and ω 1 . When ω → +∞, the horizontal asymptote of υ is still a 0 . On the other hand, the maximum attenuation occurs near ω 0 . These properties are illustrated in figures 3 and 4. Viscous case (ν = 0); see table 1 for the parameters.

In the limit-case without resonators (e = 0), the linear dispersion relation [START_REF] Flannery | Numerical Recipes in C: the Art of Scientific Computing[END_REF] simplifies to

[g -ω 2 (1 + f χ)][i d k 2 + k(a -c χ) -ω] = 0. ( 13
)
In this equation, we focus on the part corresponding to (4b), which leads to

i d k 2 + a - c (i ω) 1/2 k -ω = 0. ( 14
)
Neglecting the diffusivity of sound (d = 0) leads to the following phase velocity and attenuation

υ = a 2 ω -a c √ 2 ω + c 2 a ω -c ω/2 , α = c √ 2 ω 3/2 a 2 ω -a c √ 2 ω + c 2 . ( 15
)
In the inviscid case (c = d = 0), the phase velocity is equal to a and the attenuation is null. Otherwise, the following properties are deduced:

       υ(ω) ∼ 0 -c 2 ω , lim ω→+∞ υ(ω) = a 0 , α(0) = 0, α(ω) ∼ +∞ c a 2 ω 2 . ( 16
)
These properties are illustrated in figure 5.

Regimes of propagation

An analysis of wave regimes is presented in [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF]. Under the hypothesis of weak nonlinearity, ∂u/∂x in (4a) is replaced by -(1/a 0 ) ∂u/∂t in the terms with coefficients b, c and d. The resulting system is written in the (T, X) coordinates, where T is a non-dimensional retarded time, X is a non-dimensional slow space variable and ε is defined in (3):

T = ω t - x a 0 , X = ε ω x a 0 . (17)
The reduced variables F = O(1) and G = O(1) are introduced:

F = 1 ε γ + 1 2 u a 0 = 1 ε γ + 1 2 γ p ′ p 0 , G = 1 ε γ + 1 2 γ p p 0 = p p ′ F, (18)
leading to          ∂F ∂X -F ∂F ∂T = -δ R ∂ 1/2 F ∂T 1/2 + β ∂ 2 F ∂T 2 -K ∂G ∂T , (19a)
∂ 2 G ∂T 2 + δ r ∂ 3/2 G ∂T 3/2 + Ω G = Ω F, (19b)
with the new sets of parameters

δ R = C ε R * ν ω , β = ν d ω 2 ε a 2 0 , K = V 2A D ε , δ r = 2 r ν ω , Ω = ω 0 ω 2 . ( 20
)
The effect of dissipative terms δ r , β and δ R is analyzed in [START_REF] Sugimoto | Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves[END_REF][START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF]; in particular, β is negligible. The dynamics of the system is mainly governed by the two parameters K and Ω. K is the ratio of geometrical dispersion terms to the nonlinear terms. It is assumed that the resonators are coupled to the guide (V = 0), hence ω 0 = ∞ and 0 < Ω < +∞. Three limit-cases are then obtained in the inviscid case:

• first, if K ≪ 1 or Ω ≪ 1, one gets the evolution in a tube without resonators, leading to a shock

∂F ∂X -F ∂F ∂T = 0; (21)
• second, if K ≫ 1, ( 19) degenerates into a linear dispersive equation

∂F ∂X + K ∂F ∂T + 1 Ω ∂ 3 F ∂T 2 ∂X = 0; (22)
• third and last, if Ω ≫ 1, it is possible to obtain the Korteweg-de Vries equation

∂F ∂X + K ∂F ∂T -F ∂F ∂T = K Ω ∂ 3 F ∂T 3 . ( 23
)
This case is of particular interest, since it yields to solitons.

Mathematical modeling

Diffusive representation of fractional derivatives

Here we focus on the fractional terms in (4). The fractional integral ( 6) is non local in time, and relies on the full history of w(t). A much more efficient formulation relies on a diffusive representation of this operator, which can be written equivalently [9]

∂ -1/2 ∂t -1/2 w(t) = +∞ 0 φ(θ, t) dθ, (24)
where the diffusive variable φ defined by

φ(θ, t) = 2 π t 0 e -(t-τ ) θ 2 w(τ ) dτ (25)
satisfies the local-in-time differential equation

   ∂φ ∂t = -θ 2 φ + 2 π w, φ(θ, 0) = 0. ( 26
)
To compute the derivative of order 3/2, we differentiate the fractional integral twice [START_REF] Chester | Resonant oscillations in closed tubes[END_REF]. A first differentiation leads to the derivative of order 1/2:

∂ 1/2 ∂t 1/2 w(t) = d dt ∂ -1/2 ∂t -1/2 w(t) , = d dt 1 √ π t 0 (t -τ ) -1/2 w(τ ) dτ , = 1 √ π t 0 (t -τ ) -1/2 w ′ (τ ) dτ, (27)
where w(0) = 0 has been used. In the rest of the paper w(•) = p(x, •) and a zero initial condition p(x, 0) = 0 will always be chosen in the numerical tests. Thus, proceeding as previously for the fractional integral, a diffusive representation is introduced

∂ 1/2 ∂t 1/2 w(t) = +∞ 0 ξ(θ, t) dθ, (28)
where the diffusive variable ξ satisfies

   ∂ξ ∂t = -θ 2 ξ + 2 π w ′ , ξ(θ, 0) = 0. ( 29
)
The derivative of order 3/2 is immediately deduced [START_REF] Deü | Simulation of fractionally damped mechanical systems by means of a Newmark-diffusive scheme[END_REF]:

∂ 3/2 ∂t 3/2 w(t) = d dt ∂ 1/2 ∂t 1/2 w(t) , = +∞ 0 ∂ξ ∂t (θ, t)dθ, = +∞ 0 -θ 2 ξ + 2 π w ′ dθ. (30)

First-order system

To approximate the integral ( 24), a quadrature formula on N points is used, with weights µ ℓ and abscissae θ ℓ :

∂ -1/2 ∂t -1/2 w(t) ≃ N ℓ=1 µ ℓ φ(θ ℓ , t) = N ℓ=1 µ ℓ φ ℓ (t), (31)
where the φ ℓ satisfiy the ODE [START_REF] Rekik | Optimization of the collocation inversion method for the linear viscoelastic homogenization[END_REF]. Similarly, the integral ( 30) is written

∂ 3/2 ∂t 3/2 w(t) ≃ N ℓ=1 µ ℓ -θ 2 ℓ ξ(θ ℓ , t) + 2 π w ′ = N ℓ=1 µ ℓ -θ 2 ℓ ξ ℓ + 2 π w ′ , (32)
where the ξ ℓ satisfy the ODE [START_REF] Sousa | The controversial stability analysis[END_REF]. The determination of weights and nodes µ ℓ and θ ℓ is discussed in section 5. Then equations ( 4), (26, (29)), ( 31) and ( 32) are written as a first-order system

                                   ∂u ∂t + ∂ ∂x au + b u 2 2 = c N ℓ=1 µ ℓ φ ℓ + d ∂ 2 u ∂x 2 -e q, ∂p ∂t = q, ∂q ∂t = h u -g p -f N ℓ=1 µ ℓ -θ 2 ℓ ξ ℓ + 2 π q , ∂φ ℓ ∂t - 2 π ∂u ∂x = -θ 2 ℓ φ ℓ , ℓ = 1 • • • N, ∂ξ ℓ ∂t = -θ 2 ℓ ξ ℓ + 2 π q, ℓ = 1 • • • N, (33)
in an unbounded domain. The initial conditions are (ℓ = 1

• • • N) u(x, 0) = u 0 (x), p(x, 0) = p 0 (x) ≡ 0, ∂p ∂t (x, 0) = q 0 (x) ≡ 0, φ ℓ (x, 0) = 0, ξ ℓ (x, 0) = 0. ( 34
)
Taking the vector of (3 + 2 N) unknowns

U = (u, p, q, φ 1 , • • • , φ N , ξ 1 , • • • , ξ N ) T , (35)
the system (33) can be written in the form

∂ ∂t U + ∂ ∂x F (U) = S U + G ∂ 2 ∂x 2 U, ( 36
)
where F is the nonlinear flux function

F = a u + b u 2 2 , 0, 0, - 2 π u, • • • , - 2 π u, 0, • • • , 0 T , (37)
G is the (3 + 2 N) × (3 + 2 N) diagonal matrix diag(d, 0, • • • , 0). S is the (3 + 2 N) × (3 + 2 N)
diffusive matrix, as it contains the diffusive representation:

S =                        0 0 -e c µ 1 • • • c µ N 0 • • • 0 0 0 1 0 • • • 0 0 • • • 0 h -g - 2 π f N ℓ=1 µ ℓ 0 • • • 0 f µ 1 θ 2 1 • • • f µ N θ 2 N 0 0 0 -θ 2 1 . . . . . . . . . . . . 0 0 0 -θ 2 N 0 0 2 π -θ 2 1 . . . . . . . . . . . . 0 0 2 π -θ 2 N                        .
(38) Three properties are deduced from (36):

• the eigenvalues of the Jacobian matrix J = F ′ in (37) are real: a + b u, and 0 with multiplicity 2 N + 2. These eigenvalues do not depend on the coefficients of the diffusive representation;

• if the viscosity is neglected in the resonators (f = 0), then the 3 + 2 N eigenvalues of S are (ℓ = 1, • • • , N):

Sp(S) = ±i ω 1 , 0, -θ 2 ℓ , (39)
with ω 1 = √ e h + g; see (2) and [START_REF] Chaigne | Acoustique des Instruments de Musique[END_REF]. In particular, the real parts in (39) are null or negative, which implies a decrease in energy.

• a linear dispersion analysis can be performed as in the original model [START_REF] Brillouin | Propagation des Ondes dans les Milieux Périodiques[END_REF]. The formula [START_REF] Flannery | Numerical Recipes in C: the Art of Scientific Computing[END_REF] with coefficients (9) still holds, replacing χ in (8) by

χ(ω) = 2 π N ℓ=1 µ ℓ θ 2 ℓ + i ω . ( 40
)
4. Numerical scheme

Splitting

In order to integrate the system (36), a grid was introduced, with a uniform spatial mesh size ∆x and a variable time step ∆t n , which for the sake of simplicity will be noted ∆t. The approximation of the exact solution U(x j = j ∆ x, t n = t n-1 + ∆t) is denoted by U n j . Unsplit integration of ( 36) is not optimal, because the time step stability condition involves the spectral radius of S which increases with N. Moreover, it requires building an adequate scheme for the coupled system.

A more efficient strategy, based on a splitting method, was adopted here. Instead of integrating the original equation ( 36), a propagative equation

∂ ∂t U + ∂ ∂x F (U) = G ∂ 2 ∂x 2 U, (41)
and a diffusive equation

∂ ∂t U = S U, (42)
were considered successively. The discrete operators to solve (41) and (42) were denoted by H a and H b , respectively. Strang splitting [START_REF] Leveque | Numerical methods for conservation laws[END_REF][START_REF] Holden | Operator splitting for the KDV equation[END_REF] was then used between t n and t n+1 , solving successively (41) and ( 42) with adequate time increments:

• U (1) j = H b ( ∆ t 2 ) U n j , • U (2) j = H a (∆ t) U (1) j , • U n+1 j = H b ( ∆ t 2 ) U (2) j .
(

) 43
Provided that H a and H b are second-order accurate and stable operators, the time-marching (43) gave a second-order accurate approximation of the original equation [START_REF] Toda | Nonlinear waves and solitons[END_REF].

Propagative part of the system

Equation ( 41) is solved by any standard scheme for nonlinear hyperbolic PDE:

u n+1 j = u n j - ∆t ∆x F 1 j+1/2 -F 1 j-1/2 + d ∆t ∆x 2 u n j+1 -2 u n j + u n j-1 , φ n+1 j,ℓ = φ n j,ℓ + 1 π ∆t ∆x u n j+1 -u n j-1 , ℓ = 1, • • • , N, (44)
where F 1 j±1/2 is the numerical flux function of the advection-Burgers equation in [START_REF] Yuan | A numerical scheme for dynamic systems containing fractional derivatives[END_REF]. In practice, a second-order TVD scheme with MC-limiter was used in our numerical experiments [START_REF] Leveque | Numerical methods for conservation laws[END_REF]. Stability analysis of (44) provides the necessary and sufficient condition [START_REF] Sousa | The controversial stability analysis[END_REF][START_REF] Dehghan | Weighted finite difference techniques for the onedimensional advection-diffusion equation[END_REF]

α -α 2 2 ≤ δ ≤ 1 -α 2 , (45)
with the adimensionalized parameters α and δ and the discrete velocity a

(n) max α = a (n) max ∆t ∆x , δ = d ∆t ∆x 2 , a (n) max = a + b max j (u n j ). ( 46
)
Condition (45) was proven rigorously in the case of the advection equation and the upwind scheme, but numerical experiments indicated that it still holds for the nonlinear advection (modifying a into a (n) max ) and for the TVD scheme. Solving (45)-( 46) gives the condition

α ≤ min   1 + 1 Pe , 1 1 + 1 Pe   , (47)
where Pe = α/2δ = a

(n) max ∆x/2 d is the discrete Péclet number. In our configuration, Pe ≈ 10 5 which leads to the restriction on the time step

a (n) max ∆t ∆x ≤ 1 + 1 Pe -1 ≈ 1 - 1 Pe ≈ 1. ( 48
)
Therefore despite the explicit discretization of d ∂ 2 u/∂x 2 , the optimal CFL condition is maintained.

Diffusive part of the system

Since the physical parameters do not vary with time, the diffusive part (42) can be solved exactly. This gives

H b ∆ t 2 U j = e S ∆ t 2 U j . (49)
In the inviscid case ν = 0, only the unknowns p, q and u are involved, and since N = 0, the exponential can be computed analytically. Using ω 1 (2) and defining τ = ∆t/2, one obtains:

e Sτ =         1 ω 2 1 (g + e h cos ω 1 τ ) e g ω 2 1 (1 -cos ω 1 τ ) - e ω 1 sin ω 1 τ h ω 2 1 (1 -cos ω 1 τ ) 1 ω 2 1 (e h + g cos ω 1 τ ) 1 ω 1 sin ω 1 τ h ω 1 sin ω 1 τ - g ω 1 sin ω 1 τ cos ω 1 τ         . (50)
In the general case N > 0, the exponential is computed numerically using a (6, 6) Padé approximation in the "scaling and squaring method" [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later[END_REF]. If the physical parameters are constant, the computation is done only once at each time step, leading to a negligible computational cost. Even in the case N = 0, using the numerical evaluation of e Sτ is twice as fast as computing (50), because of the numerical evaluations of trigonometric functions. This part of the splitting is unconditionally stable, so that the global stability requirement is (48) and is not penalized by the diffusive part. In other words, the time step depends only on the advection and on the Burgers coefficient in (4). In particular, ∆t does not depend on the fractional parameters c and f or on the coupling parameters e and h.

Coefficients of the diffusive representation

The 2 N coefficients of the diffusive representation µ ℓ and θ ℓ in ( 38) have yet to be determined. These coefficients are derived from [START_REF] Monkewitz | The response of Helmholtz resonators to external excitation. Part 1. Single resonators[END_REF], [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF] and [START_REF] Sugimoto | Acoustic solitary waves in a tunnel with an array of Helmholtz resonators[END_REF], and they are used to approximate improper integrals of the form

+∞ 0 φ(θ) dθ ≃ N ℓ=1 µ ℓ φ(θ ℓ ), (51)
where time t has been omitted for the sake of simplicity. This issue is crucial both for the accuracy of the modeling and for the computational efficiency of the method. Many strategies exist for this purpose. We will begin by recalling three known methods based on orthogonal polynomials, and then we will propose another method based on optimization. To be consistent with the literature, the following notations are introduced: α is the order of the fractional derivative; ⌈α⌉ is the ceiling function that rounds up to the next integer; and lastly

α = 2 α -2⌈α⌉ + 1. ( 52
)
Since α = 1/2 or 3/2 in (4), then α ≡ 0.

Method 1: Gauss-Laguerre

This algorithm is proposed in [START_REF] Yuan | A numerical scheme for dynamic systems containing fractional derivatives[END_REF]. The improper integral (51) is evaluated with the Gauss-Laguerre quadrature:

+∞ 0 θ γ e -θ ψ(θ) dθ ≃ N ℓ=1 w ℓ ψ(z ℓ ), ( 53
)
where γ is a parameter, w ℓ are the weights and z ℓ are the nodes [START_REF] Flannery | Numerical Recipes in C: the Art of Scientific Computing[END_REF]. It implies

+∞ 0 φ(θ) dθ ≃ N ℓ=1 w ℓ z -γ ℓ e z ℓ φ(z ℓ ), (54)
and consequently the desired coefficients in (51) are

µ ℓ = w ℓ z -γ ℓ e z ℓ , θ ℓ = z ℓ . (55)
Very slow convergence, with the usual value γ = 0, was observed by many authors. In [START_REF] Diethelm | An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives[END_REF], two problems were identified :

• φ(θ) ∼ 0 θ α , and since α ∈] -1, 1[, an integrable singularity may occur at 0 if α = 0. Taking γ = α in (54)-(55) eliminates this problem;

• the diffusive variable φ (25) decreases polynomially: φ(θ) ∼ +∞ θ α-2 ≡ 1/θ 2 , which is badly represented by Gauss-Laguerre exponential weight. This problem cannot be solved.

Method 2: Gauss-Jacobi

A more efficient approach has been proposed and analysed in [START_REF] Diethelm | An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives[END_REF]. The improper integral (51) is evaluated with the Gauss-Jacobi quadrature:

+∞ 0 φ(θ) dθ = +1 -1 (1 -z) γ (1 + z) β φ(z) dz, ≃ N ℓ=1 w ℓ φ(z ℓ ), (56)
where

φ(z) = 2 (1 -z) γ (1 + z) β+2 φ 1 -z 1 + z . ( 57
)
The quadrature coefficients in (51) are deduced:

µ ℓ = w ℓ 2 (1 -z ℓ ) γ (1 + z ℓ ) β+2 , θ ℓ = 1 -z ℓ 1 + z ℓ . (58)
In [START_REF] Diethelm | An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives[END_REF], it is proposed to take γ = α ≡ 0 and β = -α ≡ 0 in (56)-( 58). The nodes of Gauss-Jacobi quadrature (58) cover a much wider interval in the θ-space than those of the Gauss-Laguerre quadrature (55), which explains qualitatively why the slowly decreasing diffusive variables φ (25) are better approximated.

Method 3: modified Gauss-Jacobi

In [START_REF] Birk | An improved non-classical method for the solution of fractional differential equations[END_REF], an improvement is proposed for method 2 , which consists in widening the range of nodes. Based on a modified Gauss-Jacobi quadrature, the function φ in ( 56) is now

z(θ) = 4 (1 -z) γ-1 (1 + z) β+3 φ 1 -z 1 + z 2 . (59)
The quadrature coefficients in (51) are deduced:

µ ℓ = w ℓ 4 (1 -z ℓ ) γ-1 (1 + z ℓ ) β+3 , θ ℓ = 1 -z ℓ 1 + z ℓ 2 . (60)
In [START_REF] Birk | An improved non-classical method for the solution of fractional differential equations[END_REF], it is proposed to take γ = 2 α + 1 ≡ 1, β = -(2 α -1) ≡ 1 in (59)-(60).

Method 4: optimization

Finally, we propose a fourth and last method based on the dispersion relation [START_REF] Flannery | Numerical Recipes in C: the Art of Scientific Computing[END_REF]. The original problem (4) and the first-order system (36) differ only in their symbol χ(ω): (8) in the first case, (40) in the second one. Adjusting them provides a way to estimate µ ℓ and θ ℓ . This technique is physically meaningful, and has proven its efficiency in a previous work about poroelastic waves [START_REF] Blanc | Biot-JKD model: simulation of 1D transient poroelastic waves with fractional derivatives[END_REF]. Let Q(ω) be the optimized quantity and Q ref the desired one:

Q(ω) = χ(ω) χ(ω) = 2 π N ℓ=1 µ ℓ θ 2 ℓ + i ω (i ω) 1/2 = N ℓ=1 µ ℓ q ℓ (ω), Q ref (ω) = 1. (61)
We implement a linear optimization procedure [START_REF] Groby | A time domain method for modeling viscoacoustic wave propagation[END_REF][START_REF] Blanc | Biot-JKD model: simulation of 1D transient poroelastic waves with fractional derivatives[END_REF] in order to minimize the distance between Q and Q ref in the interval [ω min , ω max ] containing the characteristic angular frequency ω of the initial pulse. The abscissae θ ℓ are chosen so that they are distributed linearly on a logarithmic scale

θ 2 ℓ = ω min ω max ω min ℓ-1 N-1 , ℓ = 1, ..., N. (62)
The weights µ ℓ are obtained by solving the system

N ℓ=1 µ ℓ q ℓ (ω k ) = 1, k = 1, ..., K, (63)
where the ωk are also distributed linearly on a logarithmic scale of K points ωk = ω min ω max ω min

k-1 K-1 , k = 1, ..., K. (64)
Since the q ℓ (ω) are complex functions, optimization is performed simultaneously on the real and imaginary parts

           N ℓ=1 µ ℓ Re(q ℓ (ω k )) = 1, N ℓ=1 µ ℓ Im(q ℓ (ω k )) = 0, k = 1, ..., K. (65)
A square system is obtained when 2K = N, whereas 2 K > N yields an overdetermined system, which can be solved by writing normal equations [START_REF] Flannery | Numerical Recipes in C: the Art of Scientific Computing[END_REF]. Higher accuracy is obtained with K = N, so we will make this choice in numerical experiments. The interval of optimisation [ω min , ω max ] depends on the configuration under study:

• for the coupled system with resonators, the attenuation is bounded (figure 4), and the existence of smooth solitary waves globally maintains the frequency content of the initial disturbance. Consequently, we chose a narrow interval centered around ω, by taking for instance ω min = ω / 2 and ω max = ω × 3/2;

• for the tube without resonators, shocks are expected. Consequently, higher harmonics are generated, and we proposed to use ω min = ω / 2 but ω max = ω × N , where N is the number of harmonics of interest in a Fourier decomposition of the wave [START_REF] Menguy | Weakly non-linear gas oscillations in airfilled tubes ; solutions and experiments[END_REF]. In numerical experiments, we took N = 20. [START_REF] Garnier | Asymptotic transmission of solitons through random media[END_REF]. Comparisons between the exact fractional model ( 4) and the approximate diffusive model [START_REF] Toda | Nonlinear waves and solitons[END_REF]. GL = Gauss-Laguerre (method 1); GJ = Gauss-Jacobi (method 2); GJm = Gauss-Jacobi modified (method 3); opti=optimisation (method 4).

The accuracy of the methods presented along section 5 is illustrated in Figure 7. The exact values of the phase velocity and of the attenuation are shown in solid red lines, and correspond to the cases shown in figures 3, 4 and 5. In the case with resonators (top), N = 6 diffusive variables are used, leading to 15 unknowns in [START_REF] Tao | Why are solitons stable ?[END_REF]. In the case without resonators (bottom), N = 12 diffusive variables are used, leading to 27 unknowns in [START_REF] Tao | Why are solitons stable ?[END_REF]. Increasing accuracy is observed passing successively from Gauss-Laguerre to Gauss-Jacobi, to Gauss-Jacobi modified, and then to optimization. Consequently, the optimization procedure was chosen in the numerical experiments presented in section 6. The physical and geometrical parameters are given in table 1. The physical data correspond to air at 15 • C, and the geometrical data are from [START_REF] Sugimoto | Verification of acoustic solitary waves[END_REF]. a 0 = 345.25 m/s, C = 1.478 and ν d = 3.92 10 -5 m 2 /s are obtained from [START_REF] Blanc | Biot-JKD model: simulation of 1D transient poroelastic waves with fractional derivatives[END_REF]. The parameters of ( 4)-( 5) are deduced. They are given in the upper part of table 2. The tube L = 80 m in length is discretized on N x = 8000 grid nodes. The maximal CFL number is 1-1/Pe, and the Péclet number is Pe = 1.75 10 5 (48). The CFL number is taken equal to 0.95. Lastly, a set of 10 receivers is put on the computational domain at abscissas x r = 15 + 5 (i-1), i = 1 • • • 10, where the time history of u was recorded at each time step.

Numerical experiments

Except in section 6.3, computations were initialized by a Gaussian pulse or by a rectangular force pulse on the velocity

u 0 (x) =          u m exp - x -x 0 σ 2 , u m H x -x 0 - λ 2 -H x -x 0 + λ 2 , ( 66
)
where H is the Heaviside function and x 0 = 7 m. λ is the width of the rectangular force pulse, and is also the width of the Gaussian pulse: taking σ = λ/2 √ ln 100 gives u 0 (x) = u m /100 at x = x 0 ± λ/2. All the other initial conditions in [START_REF] Sugimoto | Verification of acoustic solitary waves[END_REF] The key parameters governing the evolution of the system were K and Ω (20); see section 2.4. On the one hand, we took K = 0.5 which ensures a nonlinear regime of propagation, hence the amplitude u m ≈ 56.12 m/s (66), the parameter of nonlinearity ε = 0.19, the overpressure p ′ /p 0 = 0.22 and the sound intensity I = 181.1 dB (3). On the other hand, we took Ω = 1 or Ω = 16, yielding theoretically to dispersive waves and solitons, respectively. In the case Ω = 1, the linear dispersion analysis predicts a maximal attenuation. In the case Ω = 16, the natural frequency is ω = 411.4 rad/s and the central wavelength is λ = 5.27 m. The cut-off angular frequency is ω * = 15881 rad/s, so that the 1D approximation is justified (section 2.2).

Test 1: nonlinear acoustics

In a first test, the height of the cavity was H = 0, hence e = 0 and no coupling occured with the Helmholtz resonators. The dispersion relation was therefore [START_REF] Haddar | Efficient solution of a wave equation with fractional-order dissipative terms[END_REF], and the key adimentionalized parameters in section 2.4 were K = 0 and Ω = +∞. The number of diffusive variables was N = 12; since no memory variables ξ are required to model dissipative effects in the resonators, only 15 variables were involved in [START_REF] Tao | Why are solitons stable ?[END_REF]. Optimisation of the coefficients µ ℓ was performed between ω min = 185 rad/s and ω max = 8228 rad/s (section 5.4).

In the inviscid case, where only the coefficients a and b are non-null in (4), smooth initial data develop shocks in finite time, yielding a decrease in energy. The Gaussian initial pulse breaks at time

t * = exp(1) 2 σ u m b . (67)
It gives the time break t * = 0.021 s (67).

Figure 8 shows the initial values of the solution at the initial instant (iii) and at t = 0.16 s > t * (iii-iv), which corresponds roughly to 7000 time steps. Both the inviscid case and the viscous case are displayed, where the viscous boundary layer in the tube and the diffusivity of sound are accounted for. In the inviscid case, typical nonlinear phenomena are observed: shock on the right part of the Gaussian pulse (iii), rarefaction waves and rightgoing shock for the rectangular force pulse (iv). In this latter case, good agreement is obtained with the exact solution. In the viscous case, these phenomena are qualitatively maintained. A small decrease in amplitude is observed, together with a tail on the left part of the waves. The key issue is that viscous effects are not sufficient to prevent from the occurrence of a shock (iii) or to smooth an existing discontinuity (iv), which confirms the theoretical analysis performed in [START_REF] Sugimoto | Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves[END_REF].

Lastly, time evolution of the energy

E (n) = j (u n j ) 2 is displayed in (v-vi).
For the Gaussian pulse in the inviscid case (v), energy is conserved as long as the wave is smooth; at the scale of the figure, numerical diffusion is not seen. From t = t * where the shock appears, energy decreases. For the rectangular force pulse in the inviscid case (v), energy decreases linearly with time. The results are qualitatively the same in the viscous cases, with a greater decrease in energy.

Test 2: fractional oscillations in Helmholtz resonators

In this second validation test, we focused on the fractional oscillator of order 3/2. The coupling with nonlinear acoustics in the tube was neglected, and (4b) was solved with no coupling h = 0. The initial value of pressure was p 0 (x) = 1 and p 1 (x) = 0, leading to oscillations with damping. The analytical solution is obtained in terms of fractional power series on 200 modes: see equation ( 10) of [START_REF] Deü | Simulation of fractionally damped mechanical systems by means of a Newmark-diffusive scheme[END_REF]. The numerical solution is obtained by solving the following problem with N + 2 unknows

                 ∂p ∂t = q, ∂q ∂t = -g p -f N ℓ=1 µ ℓ -θ 2 ℓ ξ ℓ + 2 π q , ∂ξ ℓ ∂t = -θ 2 ℓ ξ ℓ + 2 π q, ℓ = 1 • • • N. (68)
Therefore the numerical solution is the exact solution of the system (49) with N memory variables; consequently, the only error is the error of model ε m , due to the quadrature of the diffusive representation.

Figure 9 shows the influence of the quadrature rule on the accuracy of the model. When optimization is used, only N = 2 memory variables are required to obtain an excellent agreement with the exact solution (left row). On the contrary, large errors are observed when Laguerre quadrature is used, even with N = 20.

Test 3: linear dispersive waves

In a third test, we considered the coupled system with resonators, in the linear regime: b = 0 in (4). The simulations were initialized by a Gaussian pulse, and they were performed over 8000 time steps. in [START_REF] Tao | Why are solitons stable ?[END_REF]. Optimisation of the diffusive coefficients µ ℓ is performed between ω min = 822 rad/s and ω max = 2468 rad/s if Ω = 1, and values 4 times smaller if Ω = 16 (section 5).

The case Ω = 1 is shown in figure 10. High dispersion is observed in the inviscid case, which confirms the dispersion analysis performed in section 2.3 (see figure 3, near the vertical dotted lines). The oscillations are highly damped in the viscous case, due to the large value of attenuation (see figure 4).

The case Ω = 16 is displayed in figure 11. Compared with figure 10, the dispersion is greatly reduced. In the inviscid case, an oscillating mode remains at the place of initialization; moreover, the energy is conserved (not shown here). The static mode is damped in the viscous case. Seismograms are built from the time signals stored at the receivers. Then, the celerity V of the highest amplitude is numerically measured. One obtains V = 312.02 m/s (if ν = 0) and V = 310.39 m/s (if ν = 0). These values are close to the zerofrequency limit υ = 314.53. The slight difference is due to the large-band of the intial pulse.

Test 4: acoustic solitary waves

Lastly, we considered the coupled system in the nonlinear regime. Computations were initialized by a Gaussian pulse, and simulations were performed during 8000 time steps, which corresponds roughly to 0.20 s of propagation. The case K = 0.5 and Ω = 1 is shown in figure 12, to be compared with the case presented in figure 10 obtained in the linear regime. A ballistic signal was observed, followed by a highly dispersive coda: no solitary wave emerged. If ν = 0, a large amount of attenuation was also introduced, which damped this coda. In the inviscid case, the energy began to decrease and then was almost conserved. This means that the initial smooth pulse had led to a shock, and then an equilibrium with dispersion had led to the emergence of a smooth structure. Similar conclusion is reached in the viscous case, except that energy always decreased.

The case K = 0.5 and Ω = 16 is displayed in figure 13, to be compared with the case Ω = 16 in figure 11 obtained in the linear regime. We recall that the theoretical analysis predicts the existence of solitary waves (section 2.4).

Compared with what can be seen in figure 12, the coda has disappeared. In the inviscid case, an oscillating mode remains at the place of initialization; moreover, the energy is conserved (not shown here), which indicates that no shock has been created. Two smooth structures are observed. Longer simulations show that these two components separate and propagate at different speeds. In the sequel, we will examine whether these solitary waves have the classical properties of solitons.

In the case K = 0.5 and Ω = 16, seismograms are built from the time signals stored at the receivers. The celerity V of the nonlinear wave with the highest amplitude is measured numerically. Similar measures are done for various amplitudes u m of the incident pulse, from 10 m/s to 100 m/s, or equivalently from K = 2.80 to K = 0.28. It is observed that V increases linearly with u m : a linear regression estimation yields V = 313.58+0.4435 u m (if ν = 0) and V = 312.17 + 0.3773 u m (if ν = 0). Waves propagate slightly faster in the inviscid case, because attenuation decreases the amplitude and consequently the celerity. The limit for u m = 0 is close to the value obtained in the linear case (test 3).

In figure 13, it can be seen that the original Gaussian pulse of u separated into two smooth structures. The taller one was thinner and traveled faster than the shorter one. At the last instant of simulation, we inverted these waves to initialize a new computation (top of figure 14). In the inviscid case, we observe that the two waves interact like classical solitons, exchanging their shape [START_REF] Leveque | Phase plane behavior of solitary waves in nonlinear layered media[END_REF][START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF]. After separation, the shape of the wave is the same as that of the original wave, though it was shifted in location from where it would be without interaction. These behaviors were qualitatively maintained in the viscous case, even if it is less clear due to the attenuation of these waves. Collision between two solitary waves. In green and red: location of the waves if they were alone.

Conclusion

We have considered nonlinear acoustic waves in a tube connected with Helmholtz resonators. Various challenging physical features are involved: nonlinearity due to the amplitude of the waves, dispersion induced by the resonators, and fractional derivatives of order -1/2 and 3/2 due to viscous losses.

Our original contribution was to propose an efficient and accurate numerical modeling of this configuration. Some tools are conventional (TVD scheme for the nonlinear hyperbolic part), some others are novel (diffusive representation of fractional derivatives). To our knowledge, it is the first time that a diffusive representation has been considered together with the advection-Burgers equation. Eventually, a splitting strategy has ensured an optimal CFL condition for an explicit scheme. The proposed approach is computationaly efficient: the CFL stability condition is only governed by the nonlinearity (as in the usual advection-Burgers equation), and a minimum number of supplementary arrays is required to discretize the fractional derivatives.

This work was motivated by the experimental configuration shown in figure 1, previously used in the linear propagation regime [START_REF] Richoux | Propagation of mechanical waves in a one-dimensional nonlinear disordered lattice[END_REF][START_REF] Richoux | Acoustic wave dispersion in a one-dimensional lattice of nonlinear resonant scatterers[END_REF], and currently investigated in the nonlinear propagation regime. Our objective was to provide an efficient and accurate numerical modeling, validating (or not) the model [START_REF] Brillouin | Propagation des Ondes dans les Milieux Périodiques[END_REF] and the underlying hypotheses. The numerical experiments showed that the viscous effects do not modify qualitatively the wave phenomena. Consequently, the theoretical predictions made in the inviscid case about the existence of acoustic solitons were also obtained in the viscous case [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF].

A first extension of this work concerns the coefficients of the diffusive representation. In some linear problems, is is possible to determine the time evolution of the energy and to prove that this energy decreases as soon as all the coefficients µ l > 0 (l = 1 • • • N) of the diffusive representations of the fractional derivatives are positive [START_REF] Haddar | Efficient solution of a wave equation with fractional-order dissipative terms[END_REF][START_REF] Deü | Simulation of fractionally damped mechanical systems by means of a Newmark-diffusive scheme[END_REF]. Therefore methods leading to positive µ l > 0 are usually prefered. As already pointed out, it is the case for the Laguerre method but not for the optimization method we used. One perspective is to develop an alternative method ensuring the positivity of the coefficients µ l . A first possibility is to use an analytical method consisting in appoximating the function χ(ω) by rational fractions. A second approach consists in using an optimization process with a positivity constraint, for instance a Shor algorithm [START_REF] Rekik | Optimization of the collocation inversion method for the linear viscoelastic homogenization[END_REF] initialized with the results from the Laguerre method. In the case of propoelasticity [START_REF] Blanc | Numerical modeling of anisotropic poroelastic waves in the full range of frequencies[END_REF], this method led to 10 up to 100 times more accurate results.

From a physical point of view, the dissipation effects also require further investigations. To get solitons, it is necessary to have Ω ≫ 1 [START_REF] Lu | Wave field simulation for heterogeneous porous media with singular memory drag force[END_REF], which implies ω ≪ ω 0 . In this regime, the dispersion analysis indicates that the attenuation is quite low (left part of figure 3). On the contrary, a high attenuation of waves is observed experimentally. A possible explanation of this mismatch is that some mechanisms of attenuation are not incorporated in the model. A good candidate is given by turbulence and nonlinear losses in the resonators. To account for these losses, a nonlinear theory for the response of the resonators has been proposed in the appendix of [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF]. Equation (4b) with notations (5) should be replaced by the nonlinear fractional ODE

where η is determined experimentally (η ≈ 0.82 r). The term ∂ 2 (p) 2 /∂t 2 models the nonlinearity due to the adiabatic process in the cavity, whereas the semi-empirical term depending on the sign of ∂p/∂t accounts for the jet loss resulting from the difference in inflow and outflow patterns [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF]. A more sophisticated numerical method must be developed to integrate (69).

A last extension of our work concerns the case where the height H of each resonator may vary with position, leading to variable coefficients in (4). Numerically, this requires smooth functions e(x), g(x) and h(x), for instance with cubic splines, to be built. The exponential of S in (38) and (49) needs to be computed at each grid node and at each time step, which increases the computational cost, but no other modifications are required. It will make it possible to investigate numerically the propagation of acoustic solitons in random media [START_REF] Garnier | Asymptotic transmission of solitons through random media[END_REF]. This topic is a subject of intense research in various fields of physics, with possible applications in the transport of information.

Figure 1 :

 1 Figure 1: Tube connected with Helmholtz resonators (courtesy given by O. Richoux).

Figure 2 :

 2 Figure 2: Sketch of the guide connected with Helmholtz resonators.

Figure 4 :

 4 Figure4: attenuation in the linear regime, with resonators (V = 0 ⇔ e = 0). Viscous case (ν = 0); see table 1 for the parameters. The vertical dotted lines denote ω 0 and ω 1 .

Figure 5 :

 5 Figure5: phase velocity and attenuation of the tube without resonators (V = 0 ⇔ e = 0). Viscous case (ν = 0); see table 1 for the parameters.

Figure 6 :

 6 Figure 6: error of model due to the optimization procedure, for the coupled system (tube with resonators, V = 0). Left: |Q(ω) -1| in (61), for various numbers of diffusive variables N ; the vertical dotted lines denote the upper and lower ranges of optimisation ω min and ω max . Right: ε m = ||Q(ω) -1|| L2 , in terms of N ; the slope of the linear regression is -5.

Figure 6 Figure 7 :

 67 Figure 6 illustrates how the number of diffusive variables influences the accuracy of the optimiaation procedure, in the case of the coupled system. In the left part, we show the error |Q(ω) -1| for various values of N. The vertical dotted lines represent the range of optimisation. By construction, the error vanishes at the abscissae ωk . As expected, the accuracy of the diffusive approximation increases with N. The right part of the figure displays the error of model ε m = ||Q(ω) -1|| L 2 on the range of interest, in terms of N and in log-log scale. The measured values are close to a straight line with slope -5, hence one can postulate a power-law ε m ≈ ε 0 (1/N) 5 .

Figure 8 :

 8 Figure 8: test 1. Nonlinear acoustic waves in the absence of oscillators, for a Gaussian pulse (left row) and a rectangular force pulse (right row). (i-ii): initial pulse. (iii-iv): numerical and exact solution at t = 0.16 s. (v-vi): time evolution of the energy; in (v), the vertical dotted line denote the time of shock t * .

Figure 9 :

 9 Figure 9: test2. Fractional oscillator of order 3/2. Comparisons between the numerical solution and the exact solution (right row: zoom).

Figure 10 :

 10 Figure 10: test3. Linear coupled system, with a Gaussian pulse and Ω = 1. Snapshots of u after 8000 time steps. Left: invicid case; right: viscous case.

Figure 11 :Figure 12 :Figure 13 :

 111213 Figure 11: test3. Linear coupled system, with a Gaussian pulse and Ω = 16. Left: invicid case; right: viscous case. Top: snapshot of u after 8000 time steps (the vertical dotted lines denote the receivers); bottom: seismograms.

Figure 14 :

 14 Figure 14: test4. Nonlinear coupled system, with a Gaussian pulse, K = 0.5 and Ω = 16.Collision between two solitary waves. In green and red: location of the waves if they were alone.

Table 1 :

 1 physical parameters of air at 15 • C, and geometrical data from[START_REF] Sugimoto | Verification of acoustic solitary waves[END_REF].

	6.1. Configuration
	γ	p 0 (Pa) ρ 0 (kg/m 3 ) P r	ν (m 2 /s) µ v /µ
	1.403	10 5	1.177	0.708	1.57 10 -5 0.60
	R (m) D (m)	r (m)	L (m) r h (m)	H (m)
	0.04	0.05	0.00355	0.0356 0.0125	0.1

Table 2 :

 2 are null. coefficients of the PDE (4). Lower part: natural angular frequencies (2), and low frequency limit of the phase velocity (11) in the coupled system.

	a (m/s)	345.25
	b	1.20
	c (m/s 3/2 )	50.42
	d (m 2 /s)	1.96 10 -5
	e (m/s 3 /Pa) 2.40 10 -4
	f (s -1/2 )	2.23
	g (s 2 )	2.70 10 6
	h (Pa/m/s) 1.10 10 9
	ω 0 (rad/s)	1645.67
	ω 1 (rad/s)	1724.16
	υ (m/s)	314.53

  )

		2.5E5
				inviscid	8.0E5
		2.0E5		viscous
					6.0E5
	E (J)	1.5E5		E (J)	4.0E5	inviscid
		1.0E5				viscous
		5.0E4			2.0E5
		0			0
		0.02 0.04 0.06 0.08	0.1	0.12 0.14 0.16		0.02 0.04 0.06 0.08	0.1	0.12 0.14 0.16
		t (s)				t (s)

Acknowledgments. This study has been initiated with Agnès Maurel (ES-PCI, France), Olivier Richoux and Vincent Pagneux (LAUM, France), and has received financial support from the Agence Nationale de la Recherche through the grant ANR ProCoMedia, project ANR-10-INTB-0914. Pierre Haldenwang (AMU, France) is thanked for his insights into hyperbolic equations.


Patrice Salzenstein
email: patrice.salzenstein@femto-st.fr

Michel Mortier

Hélène Sérier-Brault

Rémi Henriet

Aurélien Coillet

Kouomou Yanne

Alphonse Chembo

Yannick Rasoloniaina

Patrice Dumeige

Féron

Yanne K Chembo

Alphonse Rasoloniaina

Yannick Dumeige

Patrice Féron

Coupling of high-quality-factor optical resonators

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.

Introduction.

Optical resonators with high quality factor (Q-factor) present a great interest for several applications in optoelectronic [1], metrology or fundamental physics [2,3]. For instance they can be used in optical filtering, all optical switching, low threshold non-linear optics, etc. Main results have been obtained since 1996 [4]. Recently results have been obtained with mono-crystalline such as MgF 2 , Quartz or CaF 2 whispering gallery mode (WGM) resonators [5][6][7][8][9]. Rare earth element doped micro spheres are tiny lasers which intrinsically offer narrow line width emission [10]. A key feature of the characterization of the resonators is to perform a good coupling. Performing this characterization is quite challenging because it is a necessary but difficult technique. That's why we improve theoretically and experimentally the problem of the coupling between a high Q-factor resonator and its external coupler. They are characterized by different methods such as slow and fast frequency sweeping [11]. Resonators can be coupled to tapered fibres or prisms. We focus here on characterizations techniques of passive optical resonators or active (Erbium doped ZBLALiP [12,13]) WGM resonators in terms of Q-factor determination.

Properties of a doped micro resonator.

Figure 1 shows schematically the WGM micro resonator and its access line. It is possible to isolate each resonance, although the resonator is inherently multi mode. We can work with each resonance separately. The studied resonator is then considered single mode.

A mode with an amplitude a(t) is characterized by the time life of the field τ 0 limited by losses, by the coupling rate of the access line which is proportional to 1/τ e and by the angular frequency of the resonance ω 0 . The amplitude of the input signal is noted s in (t)

 = - +     (1)
In equation ( 1), τ =( 1/τ 0 + 1/τ e ) -1 is the time life of the field in the cavity related to factor Q=ω 0 τ /2.

We assume that input signal is harmonic s in (t)= S in e jωt and output signal is s out (t)= S out e jωt . As a consequence, the transfer function is given by :

( ) ( ) 2 2 2 2 0 2 2 2 0 1/ 1/ 4 ( ) 1/ 1/ 4 e out in e s T s τ τ π δ δ τ τ π δ - + = = + + (2)
In this expression δ is the frequency mismatch given by ω = ω 0 +2πδ. Different configurations are possible according to the relative value of τ e and τ 0 . If τ 0 → ∞ then resonator has no losses, so T(δ)= 1 as illustrated on figure 1(b). But if we suppose that τ e = τ 0 , then the resonator is in critical coupling and transmission at the resonance T(0) becomes zero. And if τ 0 < -τ e , resonator presents some gain (even if it stays less than the LASER threshold) and it behaves like an amplifier, then T(0) > 1 [14].

Figure 1. Scheme of the micro resonator and its access line.

Measurements.

Amplification measurements were performed in Erbium doped ZBLALiP micro spheres. Diameters of the micro spheres are between 50 µm and 150 µm. Theses spheres are obtained by fusion in a microwave inductively coupled plasma [10]. Usual measurement method consists in sweeping the frequency of a Laser emitting around 1550 nm, in order to obtain a description of one resonance for the micro sphere. Simultaneously, this micro sphere is pumped by a laser at 1480 nm emitting in opposition of propagation as described in Figure 2. Thanks to the pumping, the sphere delivers a Laser signal for one or several resonances far enough from the studying resonance: it helps to stabilize the gain. Optical Spectrum analyzr (OSA) helps for the control the micro sphere Laser emission. To illustrate how this technology can be applied, laser frequency is swept slowly enough to stay in stationary approximation, to describe thin resonance as illustrated in Figure 3. Our approach is different. The analysis begins by sweeping very quickly the frequency of the Laser. This process is analogous to what occurs when the speed is slow, but provides definitely an easier and unquestionably a precise determination of the Q-factor. This approach is implemented by utilizing the following process. This allows to work in the ringing regime of the micro sphere. It has already been demonstrated that it is possible to determine characteristic parameters of the micro sphere [11,15] if we compare the profile of the ringing oscillations to an analytical model obtained thanks to Eq.( 1) and to an incident field like ( )

( ) i S j V t t in in s t S e ω π + =
. The sweeping speed of the Laser V S , is often difficult to obtain precisely. But it is directly obtained by the method describes above. The results are illustrated in Figure 4, where two time profiles are obtained on the same resonance, with two different speeds (the second one is two time more quickly) for sweeping the Laser. While the deduced speeds are in the good ratio, adjustment provides almost the same results for τ 0 and τ e . Observed differences for τ 0 and τ e are mainly due to the fluctuations of the taper position. It has an influence on the signal coupling, and also on the pump. When performing the adjustment from the obtained ringing oscillation profiles, it has to be underlined that we deduce a quality factor as high as Q = 5.8×10 8 . The analysis proceeds by matching the theoretical and experimental curves. Q-factor estimation may be quickly obtained, and the efficiency of such a fact sweeping is really good. This is the desired result, which gives great insight into the resonator coupling and Q-factor estimation problem.

Coupling with passive resonator and integration into an oscillator .

An optoelectronic oscillator (OEO) with an optical resonator coupled to an optical fiber, delivers a free spectral range microwave signal [1]. To implement this OEO in terms of miniaturization, it is possible to use a high WGM resonator with magnesium fluoride or calcium fluoride for microwave photonics applications. We have fabricated and characterized 5 mm diameter toroid MgF 2 and 12 mm CaF 2 resonators with a good surface roughness (1 nm). Quality factor better than 3.4×10 8 for magnesium fluoride [7] and 1.4×10 8 for calcium fluoride were measured by cavity ring-down measurements. Even if an optical line can be used as a frequency selective component in the OEO loop, the use of the mini-resonator leads to a significant reduction of the oscillator size. Further work consists in measuring the phase noise on dedicated instruments developed at the laboratory [16,17]. We can also reduce or eliminate sensitivity to vibration by using a refractive index adapter gel. Q factor of resonators can still be improved by a better polishing of the crystal coupled to other treatments.

Conclusion.

The work performed here is significant in the sense that it appears that the problem of the coupling between a high Q-factor resonator and its external coupler has been theoretically and experimentally improved. The key is to recognize that oscillations of ringing induced by the sweeping of the excitation frequency of an active micro sphere are observed. Measured apparent gain is approximatively 27dB (up to 33.9 dB in stationary regime) in a 330 kHz band corresponding in the best case to a high quality factor Q = 5.8×10 8 . All of these approaches have been developed to address high quality factor resonator characterization. Addressing the combination of resonator polishing with other treatment, such as thermal annealing, to improve the surface roughness, is still an active area of research to achieve higher Q factors.

Figure 2 .

 2 Figure 2. Setup for measuring the selective gain. Probe : extended cavity laser diode with a 150 kHz spectral width, emitting around 1550 nm. OSA : Optical Spectrum Analyzer, CP: polarization controller, MUX : multiplexer 1550/1480, PWM : power meter, EDFA : optical amplifier, F : tunable filter with a 50 pm spectral width. 'pump' is a Laser diode emitting around 1480nm. The diameter of the micro sphere is 105 µm.

Figure 3 .

 3 Figure 3. Theoretical spectra obtained for 3 τ 0 .

Figure 4 .

 4 Figure 4. Obtained ringing oscillation profiles for two frequencies F driving the sweeping of the tunable Laser, PS = 125 nW, PP = 7 mW. (a) -left: F = 400 Hz, these values are deduced from the adjustment : τ0 = -40.9 ns, τe = 38.8 ns, VS = 2.2 MHz/µs, Q = 4.4 ×10 8 and T(0)= 31.2 dB. (b) -

  , we can then calculate the amplitude of

	the output signal	( ) s t out	= -	( ) s t in	+	2	τ	e	( ) a t	by integration :
				da dt		j ω	0		1 τ	( ) a t	2 e τ	( ) s t in

Acknowledgements.

Authors acknowledge the French National Research Agency (ANR) grant "ANR 2010 BLAN 0312".


Laura Oudjedi

A Nicholas

G Parra-Vasquez

Antoine G Godin

Laurent Cognet

Brahim Lounis

Metrological Investigation of the (6,5) Carbon Nanotube Absorption Cross Section

Single wall carbon nanotubes (SWCNTs) provide unique opportunities for applications in electronics, optoelectronics, photonics and photovoltaics [START_REF] De | Carbon Nanotubes: Present and Future Commercial Applications[END_REF] . They are indeed nearly ideal models of infinite π orbital conjugations where the sp² lattice provides exceptional charge carrier properties [START_REF] Avouris | Carbon nanotube electronics[END_REF] . Moreover, their true one-dimensional character confers them well-defined optical transitions from excitonic states 3, 4 5 .

While their Raman and photoluminescence properties have been extensively studied [START_REF] Jorio | Carbon Nanotubes: Advanced topics in the synthesis, structure, properties & applications Springer[END_REF] , the basic absorption properties of SWCNT are still poorly established in a quantitative manner, even if they are of prime importance for applications. Indeed, common absorption measurements on ensemble of SWCNT suffer from the problem of samples heterogeneities inherent to standard synthesis methods.

The large multiplicity of nanotube chiralities and lengths, as well as the presence of synthesis impurities impede to assess the exact amount of a SWCNT chirality contributing to an absorption signature. For (6,5) nanotubes, absorption cross-sections ranging from 0.5 10 -18 to 1.3 × 10 -17 cm² per C atom can be deduced from previous reports [START_REF] Islam | Direct measurement of the polarized optical absorption cross section of single-wall carbon nanotubes[END_REF][START_REF] Zheng | Solution Redox Chemistry of Carbon Nanotubes[END_REF][START_REF] Carlson | Fluorescence Efficiency of Individual Carbon Nanotubes[END_REF][START_REF] Schöppler | Molar Extinction Coefficient of Single-Wall Carbon Nanotubes[END_REF] considering an excitation light polarized along the nanotubes and resonant with the second order transition (S 22 ). The majority of the caveats mentioned above which are responsible for the disparity in determination can be lifted-off by performing measurements on individual SWCNTs.

Single SWCNT absorption measurements are challenging because the weak signals have to be extracted from laser intensity fluctuations (including the always present shot noise) and contributions from environment scattering [START_REF] Lefebvre | Polarized light microscopy and spectroscopy of individual singlewalled carbon nanotubes[END_REF] . Luminescence microscopy, which is an intrinsic dark field method, has been used to perform the first evaluation of for individual (6,5) nanotubes. The method was based on the study of the luminescence dependence with excitation intensity, which deviates from linearity at high excitation due to exciton-exciton annihilation processes [START_REF] Berciaud | Luminescence decay and the absorption cross section of individual single-walled carbon nanotubes[END_REF] . From the saturation intensity and the luminescence lifetime of the SWCNT was estimated to be ~ 10 -17 cm² per C atom [START_REF] Berciaud | Luminescence decay and the absorption cross section of individual single-walled carbon nanotubes[END_REF] . However, individual defects that are photoinduced upon increasing excitation intensity quench the luminescence and therefore can strongly alter the saturation behavior [START_REF] Georgi | Photoinduced Luminescence Blinking and Bleaching in Individual Single-Walled Carbon Nanotubes[END_REF][START_REF] Santos | All-Optical Trion Generation in Single-Walled Carbon Nanotubes[END_REF] . Rayleigh intensity scattering technique combined with AFM have been used to extract an estimate of ~ 2.5 10 -17 cm²/C for the resonant absorption cross section of nanotubes of different diameters [START_REF] Joh | Single-walled carbon nanotubes as excitonic optical wires[END_REF] . Recently, a value in the range of ~ 10 -17 cm²/C was reported for a large diameter (18,5) nanotube lying on opaque substrate excited around its S 33 transition using a reflective modulation imaging technique [START_REF] Christofilos | Optical Imaging and Absolute Absorption Cross Section Measurement of Individual Nano-objects on Opaque Substrates: Single-Wall Carbon Nanotubes on Silicon[END_REF] .

Given the large disparity in the reported absorption cross-section values, a metrological approach is essential to provide a definitive value. Here, we present a quantitative measurement of the absorption cross-section of chirality identified individual (6,5) nanotubes at their S 22 transition peak. We choose the (6,5) chirality as standard nanotubes since they can easily be identified by their luminescence spectra and be excited efficiently near their second order resonance S 22 using low noise CW solid-state lasers. Futhermore, they are widely studied because of their abundance in standard synthesis methods (CoMoCAT or HiPco for instance) and because their luminescence falls in the detection window of

σ C σ C σ C σ C
Silicon detectors. Our direct measurement of carbon nanotube absorption cross-section requires neither assumption on the sample content, nor on data interpretation. It was made possible by combining transmission microscopy for absorption determination with luminescence microscopy for chirality assignment and photothermal microscopy for assuring nanotube isolation.

We optimize the transmission spatial modulation technique [START_REF] Arbouet | Direct measurement of the single-metal-cluster optical absorption[END_REF] to detect the extinction cross section of small diameter SWNTs with high signal to noise ratios. For small diameter SWCNTs, the extinction cross section is dominated by the absorption as demonstrated by Rayleigh scattering experiments which reported scattering cross-section of the order of a few 10 -20 cm²/C atom for (6,4) nanotubes [START_REF] Joh | Single-walled carbon nanotubes as excitonic optical wires[END_REF] .

One limiting factor of this technique was that only relatively low frequencies (< kHz) could be used for modulating the sample position due to mass loading [START_REF] Carey | Imaging and Absolute Extinction Cross-Section Measurements of Nanorods and Nanowires through Polarization Modulation Microscopyâ€[END_REF] . Thus, these measurements suffer from 1/f noise. Our approach is based on a high frequency (f = 100 kHz) spatial modulation of a tightly focused Gaussian beam across the nanoobject (Figure 1). In practice, the spatial modulation of a 561nm laser beam is performed by varying the first-order diffraction angle of an acousto-optic modulator. The for index matching (~1.50). In the case of small nano-objects the extinction cross section is dominated by the absorption [START_REF] Bohren | Absorption and scattering of light by small particles[END_REF] (negligible intensity scattering). For a small spherical gold NP (with dimensions much smaller than the wavelength) at the position (x, y), the transmitted beam power writes with the incident beam power, the NP absorption cross-section, I(x,y) the intensity profile at the sample, and δ the amplitude of the spatial modulation along the y-axis. Assuming a Gaussian beam profile and a small beam modulation amplitude ( ) the demodulated power is proportional to the beam profile first derivative along the y-axis and writes . As shown on figure 2c-d, individual gold nanoparticles are imaged with a high signal-to-noise-ratio (>20) with short integration times of 25 ms and excitation power of 650 µW.

P trans (t) = P inc -σ Au I(x, y +δ sin(2π ft))

P inc σ Au I(x ', y') = I 0 e -2( x ' 2 +y' 2 ) w 0 2 δ << w 0 P mod (x, y) = - 8y π w 0 4 σ Au δ × P inc e -2( x 2 +y 2 ) w 0
It follows that can be directly measured from the peak-to-peak amplitude Α pp of the profile (

) given the precise determination of the beam size w 0 and displacement δ.

For this purpose, the beam profile was determined by confocal fluorescence images of 45 individual 20 nm luminescent spheres deposited on a glass coverslip (Figure 1S, Supp. Info). We find w 0 = 269 ± 6 nm (mean ± standard deviation). The beam displacement at the sample is calibrated by measuring the displacement of individual quantum dots positions found in two confocal images acquired with two different RF waves driving the acousto-optic modulator (Figure 2S, Supp. Info). From 13 acquisitions, we found δ = 36.4 ± 1.9 nm (mean ± standard deviation).

After careful calibration of the demodulated signal, one obtain a mean absorption cross-section for 10 nm gold NPs of 41±16 nm 2 (mean ± half width at half maximum, N=191) at 561 nm excitation wavelength (Figure 2e). This is in excellent agreement with the value of 42 nm 2 predicted by the Mie theory [START_REF] Bohren | Absorption and scattering of light by small particles[END_REF][START_REF] Mie | Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen[END_REF] for a particle surrounded by a medium of refraction index of 1.50. Noteworthy, the dispersion in the measured NP absorption cross-section (33%) is mainly imposed by the NP size dispersion (30% in volume), since the measurement uncertainties play here a minor role (see below).

This measurement demonstrates that our experimental procedure allows quantitative determination of small nano-object absorption cross-sections.

We then consider (6,5) SWCNTs to achieve precise characterization of their absorption cross-section.

Raw HiPco SWCNTs were solubilized in 1% aqueous sodium deoxycholate (DOC) using brief tip sonication followed by bench-top centrifugation to remove un-solubilized material from the suspension. Nanotubes are deposited on cleaned glass coverslip and covered with immersion oil for further analysis. In order to identify long individual (6,5) SWCNTs, single-molecule wide-field photoluminescence (PL) microscopy is first performed with a detection window around their S 11 emission transition (985 nm) [START_REF] Duque | Environmental and Synthesis-Dependent Luminescence Properties of Individual Single-Walled Carbon Nanotubes[END_REF] . Figure 3a shows bright individual (6,5) nanotubes displaying uniform luminescence.

The selected nanotubes are then imaged using the modulated absorption method. Since the signal amplitude and profile depend on the nanotube orientation with respect to the beam modulation axis, the measurements are performed only on long nanotubes perpendicular to y-axis (Figure ,3b).

Furthermore, the laser polarization is set along the nanotube for maximum interaction [START_REF] Lefebvre | Polarized Photoluminescence Excitation Spectroscopy of Single-Walled Carbon Nanotubes[END_REF] . To ensure that the nanotubes are isolated from other nanotube species or from catalyst impurities, always abundant in nanotubes samples, we performed photothermal heterodyne microscopy on the same sample region [START_REF] Berciaud | Photothermal Heterodyne Imaging of Individual Nonfluorescent Nanoclusters and Nanocrystals[END_REF][START_REF] Berciaud | Absorption spectroscopy of individual single-walled carbon nanotubes[END_REF] . In this case, we use circularly polarized excitation to obtain images independent from tube orientations (Figure 2c and Figure 3S, Supp. Info). The perfect correlation between PL, photothermal and modulated absorption images allows assigning the individual (6,5) nanotubes without ambiguities.

σ Au

σ Au = πe 1/2 8 w 0 3 δ A pp P inc
In contrast to gold nanoparticles, which can be considered infinitely small compared to the beam profile, the 1D geometry of the nanotube has to be taken into account in the derivation of demodulated power:

where L is the SWCNT length and the number of carbon atoms per tube length ( = 88271 µm -1 for (6,5) nanotubes). Figure 3d displays the signal profile along the y-axis measured around the middle part of the tube at excitation power of 70 µW.

In the case of nanotubes significantly longer than the diffraction limit (i.e. L>1µm), can be expressed as a function of Α pp using: . The excitation laser wavelength is being not strictly resonant with the S 22 transition; one has to account for the dependence of the absorption on the excitation wavelength to provide a value for at the peak of this transition. For this purpose, we acquired photoluminescence excitation (PLE) of individual (6,5) nanotubes prepared following the same sample procedure as for absorption measurements because the optical spectra depend on the nanotube local environment. A tunable dye laser (emission range 540-590 nm) is used for excitation and the SWCNT luminescence signal was recorded at low excitation intensities to avoid any saturation effect. PLE spectra displayed systematically a Lorentzian profile [START_REF] Cambré | Luminescence Properties of Individual Empty and Water-Filled Single-Walled Carbon Nanotubes[END_REF] with a peak S 22 transition found at 569±2 nm (Figure 4S, Supp. Info). From 8 spectra acquired, a correction factor of 1.47 ± 0.12 is deduced allowing the determination of at the S 22 peak transition using the measurements performed with low noise solid-state laser at 561 nm. Figure 3e displays the histogram measured on 73 different (6,5) nanotubes. A unimodal distribution is obtained providing a clear signature that all data points stem from individual chirality identified nanotubes. The histogram is well adjusted with a Gaussian distribution centered at 3.2 10 -17 cm 2 with a half-width-at-half maximum of 0.5 10 -17 cm 2 . This value compares to the upper range of previous indirect estimations of (Figure 3f) [START_REF] Joh | Single-walled carbon nanotubes as excitonic optical wires[END_REF] . Importantly, this obtained value originates from direct measurements (model independent), therefore, we can determine both its precision and its exactitude in a straightforward manner. The relative error on the predetermination of w 0 (2%, which translate to 4% error on ), δ (5%) and the excitation wavelength correction factor (9%) might induce a systematic bias on the determination of thus affecting its exactitude by up to 18%. Furthermore, we estimate the precision of this measured value to be 15% stemming from as a 10% noise-to-signal ratio for Α pp , a 1% precision in the determination of the beam power and 2% fluctuations on w 0 due to sample-to-sample focusing. We

P mod (x, y) = - 8 π δ w 0 3 P inc n C σ C ye -2 y 2 w 0 2 erf 2 w 0 (x + L 2 ) " # $ % & '-erf 2 w 0 (x - L 2 ) " # $ % & ' " # $ $ % & ' ' n C n C P mod (x = 0, y) σ C σ c = 2π e 1/2 8 w 0 2 n c δ A pp P inc σ C σ C σ C σ C σ C σ C
thus conclude that the measured spread of the data (±16 %) displayed on figure 4b reflects the precision on the determination of while its exactitude lies within 18%.

Correlative imaging (luminescence versus absorption) of long carbon nanotube provide a direct evidence of exciton diffusion and a straightforward determination of the exciton diffusion length before its recombination [START_REF] Wang | Observation of rapid Auger recombination in optically excited semiconducting carbon nanotubes[END_REF][START_REF] Ma | Femtosecond spectroscopy of optical excitations in single-walled carbon nanotubes: Evidence for excitonexciton annihilation[END_REF] . Previous determinations rely on stepwise quenching of nanotube luminescence by chemical reactions [START_REF] Cognet | Stepwise quenching of exciton fluorescence in carbon nanotubes by single-molecule reactions[END_REF][START_REF] Cognet | Subdiffraction Far-Field Imaging of Luminescent Single-Walled Carbon Nanotubes[END_REF] , by ensemble studies of the length dependence of PL efficiency [START_REF] Hertel | Diffusion Limited Photoluminescence Quantum Yields in 1-D Semiconductors: Single-Wall Carbon Nanotubes[END_REF] or on studying the luminescence spatial profile at nanotube ends [START_REF] Crochet | Disorder Limited Exciton Transport in Colloidal Single-Wall Carbon Nanotubes[END_REF] . In the latter case, determination of l D relies on fitting carefully the intensity profile at the nanotube end where local quenching occurs.

Here we show that by simultaneous acquisition of the luminescence and absorption profiles, l D can be directly determined by confronting the apparent nanotube ends between the two imaging modalities.

To this aim, individual nanotubes were imaged by confocal luminescence and spatial modulation absorption microscopies with the same confocal excitation beam (Figure 4a-b). Comparison of the images reveals that nanotube PL profiles appear systematically shorter than absorption ones (Figure 4c). This is a direct consequence of exciton diffusion. The apparent nanotube ends appear separated by ~ 150 nm between the two images. Using a simple 1D exciton diffusion model (Supp. Info) with quenching at nanotube ends, we simulated the expected PL and absorption profiles obtained in our confocal microscope for l D varying from 0 to 500 nm. The best agreement with the experiments is found for l D ≈ 200 nm (figure 4d). This value is in agreement with previous determinations for nanotube suspended in DOC [START_REF] Crochet | Disorder Limited Exciton Transport in Colloidal Single-Wall Carbon Nanotubes[END_REF][START_REF] Siitonen | Surfactant-Dependent Exciton Mobility in Single-Walled Carbon Nanotubes Studied by Single-Molecule Reactions[END_REF] .

In conclusion, we measured the absorption cross-section of (6,5) carbon nanotubes at their second order optical transition and obtained a value of 3.

  incident beam is focused onto the sample by means of a high NA microscope objective (x60, NA 1.45) and its transmission collected by a second identical objective is sent to a balanced photodetector connected to a lock-in amplifier for demodulation. Images of the samples are obtained by raster scanning the sample mounted on a piezo-scanning stage. The sensitivity and metrological potential of the method is first validated on a standard nanosphere sample. This sample consists of gold nanoparticles (NP) with mean diameter 10.0±0.8 nm (measured on 191 NPs by TEM, figure 2 a-b) deposited on cleaned glass coverslip, embedded in immersion oil

Figure 1 :

 1 FIGURES CAPTIONS

Figure 2 :

 2 Figure 2: a) TEM image of 10nm individual gold NP used in this study. b) Histogram of the diameters measured on 191 NP and its Gaussian fit. c) Modulated absorption image of a single 10 nm gold NP excited at 561 nm (P 0 = 650 µW, 25ms/pixel). d) Corresponding profile averaged over 4 lines of c). e) Histogram of the absorption cross-sections deduced for 100 individual 10 nm gold NPs with its Gaussian fit (see text).

Figure 3 :

 3 Figure 3: Luminescence (a, 300 ms integration time), photothermal (b, 5 ms/pixel) and modulated absorption (c, 25 ms/pixel) images of the same individual (6,5) carbon nanotube. Scale bars 1µm. d) Absorption profile perpendicular to the nanotube axis averaged over 10 lines of c). e) Histogram of the absorption cross sections deduced for 71 individual (6,5) carbon nanotubes and its Gaussian fit (see text). The mean value of the distribution is 3.2 10 -17 cm 2 and its standard deviation 0.5 10 -17 cm 2 . f)Comparison with the values reported in the literature for(6,5) nanotubes (all values were normalized to an excitation at the S 22 transition with a polarization along the nanotube axis):  6 ;  7 ;  8 ;  11 ;  9 ; 14 ; this work.

Figure 4 :

 4 Figure 4: Confocal microscopy images ((a) modulated absorption, (b) luminescence, 50ms/pixel) of an individual (6,5) carbon nanotube using the same excitation beam. c) Corresponding absorption (black) and luminescence (red) profiles along the nanotube axis. d) Simulated profiles using a 1D diffusion equation (see text and Supp. info): absorption (black) and photoluminescence (red and blue) for varying diffusion lengths (100 to 500 nm). The best agreement with experiments is obtained for l D =200 nm (red).

  Figure'2'

  2 10 -17 cm 2 per carbon atom. This constitutes the first

	metrological investigation of the absorption cross-section of chirality-identified nanotubes. Beside its
	importance for fundamental investigations of nanotube physical properties (e.g. accurate determination
	of exciton rate generation), this value is also essential for carbon nanotube material science and a key
	parameter for exploiting the large application potential that nanotubes offer in optoelectronics and
	photovoltaics.
	ACKNOWLEDGMENT. This work was funded by the Agence Nationale de la Recherche, Région
	Aquitaine, DGA and the European Research Council.

§ present address: Los Alamos National Laboratory, Los Alamos, NM 87545