Mathieu Ducousso
email: mathieu.ducousso@safrangroup.com

Eduardo Cuenca

Maxance Marmonier

Laurent Videau

François Coulouvrat

Laurent Berthe

Bulk probing of shock wave spatial distribution in opaque solids by ultrasonic interaction

Longitudinal shock wave propagation in opaque materials can be imaged only using high brilliance large X-ray facilities. [START_REF] Gupta | Realtime, high-resolution x-ray diffraction measurements on shocked crystals at a synchrotron facility[END_REF][START_REF] Schropp | Imaging Shock Waves in Diamond with Both High Temporal and Spatial Resolution at an XFEL[END_REF] Elastography, first proposed by Ophir et al., [START_REF] Ophir | Elastography: a quantitative method for imaging the elasticity of biological tissues[END_REF] is an acoustical-based technique developed to image the shear elasticity of biological tissues by mapping bulk elastic properties of materials through ultrasound (US). The general principle is to generate a relatively low frequency shear modulation, while a high-frame-rate longitudinal probe wave records the response to this modulation. [START_REF] Parker | Imaging the elastic properties of tissue: the 20 year perspective[END_REF][START_REF] Pinton | Adaptive motion estimation of shear shock waves in soft solids and tissue with ultrasound[END_REF] Magnetic Resonance Imaging can also be used instead of US. [START_REF] Mariappan | Magnetic resonance elastography: a review[END_REF] Elastography has already been used to monitor shear shock waves, [START_REF] Catheline | Observation of Shock Transverse Waves in Elastic Media[END_REF] even inside the brain. [START_REF] Espíndola | Shear shock waves observed in the brain[END_REF] [9,[START_REF] Itoh | Breast disease: clinical application of US elastography for diagnosis[END_REF] However, such imaging technique cannot be used to image laser-driven shock propagation in hard materials, as it requires a large ratio between longitudinal and shear wave velocities, of the order of 10 3 . [9,[START_REF] Itoh | Breast disease: clinical application of US elastography for diagnosis[END_REF] At the lab-scales, except this class of soft materials and also transparent ones, only shock induced surface displacement or velocities can be monitored, using optical [START_REF] Barker | Laser interferometer for measuring high velocities of any reflecting surface[END_REF] or piezo-element setups. [START_REF] Ageev | Noncontact ultrasonic acquisition of femtosecond laser-driven ablative Mbar-level shock waves on Ti alloy surface[END_REF] However, imaging at the lab-scale could efficiently bring new understandings into shock wave physics, such as caustic formation or Mach reflection in solids, into warm dense matter, for high pressure and temperature phase diagram exploration, [START_REF] Denoeud | Dynamic X-ray diffraction observation of shocked solid[END_REF] or into engineering applications of shock waves. [START_REF] Ducousso | Quantitative evaluation of the mechanical strength of titanium/composite bonding using laser-generated shock waves[END_REF][START_REF] Wu | The New Technologies Developed from Laser Shock Processing[END_REF] Burgers [START_REF] Burgers | On the transmission of sound waves through a shock wave[END_REF] and Brillouin [START_REF] Brillouin | Reflexion et refraction d'ondes acoustiques par une onde de choc[END_REF] introduced a pioneering model of the nonlinear interaction between an acoustic wave and a high amplitude shock in fluids. The detailed analysis of all possible waves (acoustic, entropy or vorticity ones) diverging from a shock, and of the distorsion of the shock itself, was achieved by Mc Kenzie and Westphal. [START_REF] Mckenzie | Interaction of Linear Waves with Oblique Shock Waves[END_REF] This theory found multiple applications in aerodynamics, [START_REF] Robinet | Shock wave instability and the carbuncle phenomenon: same intrinsic origin?[END_REF] aeroacoustics [START_REF] Mahesh | The interaction of an isotropic field of acoustic waves with a shock wave[END_REF] and magnetohydrodynamics. [START_REF] Spreiter | External aerodynamics of the magnetosphere[END_REF][START_REF] Axford | Acceleration of Cosmic Rays by Shock Waves[END_REF] Most of these applications involve strong shocks, with Mach number much above unity. In a solid medium, very few theoretical papers exist. [START_REF] Morro | Interaction of acoustic waves with shock waves in elastic solids[END_REF] The objective of this paper is to demonstrate the never explored possibility of using this interaction as a way to investigate shock propagation in an opaque solid. This is exemplified by the emission of a 15 MHz planar US elastic wave used as a probe to monitor a ns laser-driven shock wave of 0.65 GPa amplitude in a thick aluminum sample. The detection mechanism is based on the interaction of the probe wave and the shock.

Varying time-delay Δt between their respective emission times allows to reconstruct the shock propagation. Such elastic-based, shock detection scheme could complete optical-based detection and bring a new leap to the use of laser-driven shock waves.

We first introduce some physical basis related to interactions between shock and elastic plane waves. Second, 1D hydrodynamic numerical simulations are used to quantify laser-induced surface loading and the corresponding acoustic impedance variation due to the induced shock wave. [START_REF] Bardy | Development of a numerical code for laser-induced shock waves applications[END_REF] Third, experimental investigations are presented in an aluminum plate. The detection mechanism, driven by elastic scattering detection, is not straightforward and is interpreted using a 3D elastic simulation code. [START_REF] Cuenca | Propagation of laser-generated shock waves in metals : 3D axisymmetric simulations compared to experiments[END_REF] Similar experiments in titanium and water are presented in supplementary materials. [26] Finally, conclusions emphasizes on perspectives for the near future. Shock wave formation and propagation are irreversible processes. Therefore, the problem cannot be treated simply using an acoustical point of view, which is intrinsically limited to adiabatic phenomena. Equations of continuity of mass, movement and entropy must be considered, along with Rankine-Hugoniot (RH) jump relations through shocks. Their resolution gives rise to at least three remarkable points. First, an acoustic plane wave cannot reflect on the shock front, as the shock propagates at a velocity cs larger than the longitudinal speed of sound, ca. A shock front perturbation, viewed by Brillouin as an "accompanying wave", is also attached to the shock front. [START_REF] Brillouin | Reflexion et refraction d'ondes acoustiques par une onde de choc[END_REF] Second, the interaction gives rise to an entropy mode, convected by the upstream flow behind the shock. It must be considered to properly satisfy RH energy balance on the shock, [START_REF] Burgers | On the transmission of sound waves through a shock wave[END_REF] but is of negligibly small amplitude in case of a weak shock, as entropy jump then gets of third order only. In the case of oblique incidence (not considered here) a vorticity mode has also to be taken into account. Third, the transmission coefficient of the acoustic wave on the shock depends on several parameters, notably the Mach number and the nonlinear parameters of the materials. Thus, the amplitude of the scattered wave may be greater than the amplitude of the incident wave. [START_REF] Mckenzie | Interaction of Linear Waves with Oblique Shock Waves[END_REF] The laser-driven shock wave profile is quantified using ESTHER software, a 1D hydrodynamic code simulating surface sublimation due to laser-matter interaction. The resulting expanding dense plasma from the surface and, by reaction, the shock wave generation and propagation are solved by the resolution of the mechanical balance equations. The Mie-Grüneisen equation of state is used to account for non-linear elasticity. [START_REF] Bardy | Development of a numerical code for laser-induced shock waves applications[END_REF] The simulations are performed according to the experimental conditions. The calculated surface loading is presented in fig. 1(a). Lateral dimensions (20x20 cm²) are such it can be considered infinite inlateral directions.

The maximum pressure is around 0.65 GPa with a duration of around 10 ns at Full Width at Half Maximum (FWHM). Such dynamical loading is at the fringe of aluminum Hugoniot Elastic Limit (HEL), also estimated at 0.65 GPa in laser-driven loading. [START_REF] Cuenca | Propagation of laser-generated shock waves in metals : 3D axisymmetric simulations compared to experiments[END_REF] An illustration of the calculated relative variations of the density (dashed lines) and of the acoustic velocity (symbols) is presented fig. 1(b), at two times after laser illumination. At 300 ns, the shock traveled around 1.6 mm, less than the sample thickness and is a compressive shock wave. At 450 ns however, the shock has reflected on the back face, converting into an expansion or release wave. The nonlinear Mie-Grüneisen equation of state makes the velocity variation strongly dependent on the compressive (+0.3%) or release (-0.02%) nature of the shock, contrarily to the density variation (±0.5%). After 12 mm of propagation, the Mach number of the shock wave is around 0.002, with a rise time evaluated at 4 ns and a FWHM at 16 ns

[26].
Figure 2(a) presents the experimental setup. Shock waves are generated using 532 nm wavelength laser pulses of 7 ns, with energy ranging from 1 to 6 J. Fluencies up to 8 GW/cm² can easily be reached. A Diffractive Optical Element focusses the laser at the sample surface with an almost perfect circular top-hat profile of 3 mm diameter with 4.5 GW/cm 2 laser fluency (unless specified otherwise). [START_REF] Cuenca | Propagation of laser-generated shock waves in metals : 3D axisymmetric simulations compared to experiments[END_REF] In the near-field, the resulting laser-driven shock wave has a similar circular flat profile with sharp edges. Diffraction at these edges additionally generates shear waves of significant amplitude, [START_REF] Cuenca | Propagation of laser-generated shock waves in metals : 3D axisymmetric simulations compared to experiments[END_REF] that are not considered in our work as they arrive much later than the main shock. Compared to simulation, a thick sample (12 mm thickness) is used for the experiments, other parameters being identical. A phased-array (PA) is coupled by an acoustic coupling gel to the back-face of the sample, [START_REF] Drinkwater | Ultrasonic arrays for non-destructive evaluation: A review[END_REF] opposite to the laser focus. It contains 128 linearly aligned piezo-electric elements of 0.15x5 mm 2 elementary dimensions with flat surface. The central acoustic frequency is 15 MHz. The PA is used in both emission and reception modes. In emission, the US probe wave of lateral dimensions around 23x5mm² is contrapropagative with respect to the shock. The PA driver is locked onto the laser source with a time-delay driver. It allows to adjust a time delay Δt of several µs, with ps accuracy, between laser and US emission, changing the crossing point of the two waves inside the material. [START_REF] Ghrib | Generation of controlled delaminations in composites using symmetrical laser shock configuration[END_REF] Unless specified otherwise, Δt = 0. The PA driving voltage is raised to its maximum value (100 V) to emit a wave as intense as possible. The stress field is estimated (by manufacturer) at 3 MPa inside the sample, 216 times smaller than the shock amplitude. A minimum analogic gain of the PA driver is used (5 dB) in order to unsaturate the detected signal from the shock wave.

The performance analysis of the experimental setup is quantified over 30 identical sequences of generation/ detection of US and shock waves. The reproducibility of the PA generation/detection is evaluated for each element of the PA with a measured standard deviation of 10 -3 . The reproducibility of the shock generation/ PA detection is similarly measured, with a standard deviation of around 10 -2 . The difference can be explained mainly by shot-to-shot fluctuations of the laser intensity, which can be increased by the nonlinear optical processes of shock generation. [START_REF] Bardy | Development of a numerical code for laser-induced shock waves applications[END_REF] The accuracy of the overall experimental setup is therefore of the same order of magnitude as the relative impedance variation induced by the shock wave (around ±0.7%). The relative amplitude of the detected signal being around 0.1 (see later), the measured signal is one order of magnitude larger than the noise level of the detection. This relatively huge value will now be shown to allow a monitoring of shock propagation process. An alternative could have been to use the Doppler frequency shift resulting from the interaction between the counter-propagating probe wave and moving shock. It is quantified around 0.3 MHz for each round trip in the structure. However, this shift is too small with respect to the complex frequency spectrum of the probe and could not be clearly outlined experimentally. [26] A time-space representation of the involved waves, fig. 2(b), is used to infer the experimental strategy for the proposed shock elastography. The measured intensity field is , with pi the pressure field detected at a PA element (position x) for each i th back and forth travel (see fig. 2(c)). It can be decomposed into two contributions, a reference signal (noted r) and a scattered one (noted s) : u m = u r +u s . We define the relative variation ratio Qi = ui m /ui r = 1+ ui s /ui r , with superscript T (resp. A) when measuring the arrival of the transmitted probe pulse (resp. when measuring the arrival of the shock pulse and its perturbation). This is possible because probe and shock arrivals are well separated in time. In the experiments, u m = u r +u s is detected at the PA when both plane and shock waves are generated and detected. The reference ui r for Qi T (resp. Qi A ) evaluation is measured when emitting only the probe (resp. shock). The ratios Q1 A and Q1 T are presented on fig. 3(a) and 3(b), respectively, over 10 different, but identical, experimental acquisitions (blue dotted lines). Their mean values are also presented as the red thick lines. This averaging allows increasing the SNR. Detection Q1 T is strongly reproducible while more variations appear for Q1 A . This could be attributed to the particular nature of Q1 A (i.e. a shock perturbation, and not a reflected wave), and to waves superposed in the detection : probe for Q T versus shock for Q A . Both curves have a symmetric shape, centered on laser axis. The Q1 A curves feature two negative and symmetric peaks, while Q1 T presents a V shape with a clear central dip (10%). This amplitude is unexpectedly large. It is could be related to the supersonic motion of the shock discontinuity. A correct description should generalize the nonlinear approach of Burgers, Brillouin and McKenzie and Westphal. However, their simplified assumption of a static step shock has to be abandoned. Moreover, the non uniformity of the flow behind it should be considered. Indeed, a 100 ns shock duration, (Fig 1) comparable to the US probe period (66 ns at 15 MHz), may ensure a maximum interaction.

u m i (x) = ∫ p 2 i (x , t )d t
A preliminary model uses an elastic wave approach, [START_REF] Cuenca | Propagation of laser-generated shock waves in metals : 3D axisymmetric simulations compared to experiments[END_REF] replacing the shock by a 3 mm diameter inclusion with density and sound speed values provided by ESTHER (Fig. 1(b)). The resulting scattering profile (Fig. 3(b))

ideally reproduces the experimental V shape, including the strongly attenuated lateral oscillations likely due to the Fresnel diffraction of the probe on the shock circumference. [START_REF] Coulouvrat | Nonlinear Fresnel diffraction of weak shock waves[END_REF] However, a huge difference (3 decades) in amplitude remains. As already outlined, it is attributed to the over simplified linear acoustic model. Anyway, it can be attributed only to the interaction with the longitudinal shock wave as this characteristic shape is also observed in water.

[26] Signals Qi T , i = 1 to 3, for laser energies from 2 to 8 GW/cm² are reported on fig. 4(d) and compared to surface pressures simulated using ESTHER, fig. 4(c), increasing continuously with the laser intensity up to 0.75 GPa, above the material HEL (0.65 GPa). This means that plastic and elastic deformations simultaneously arise for highest energies. However, plastic deformation is confined on the very first µm under the surface while elastic deformations can propagate over several mm. [START_REF] Cuenca | Propagation of laser-generated shock waves in metals : 3D axisymmetric simulations compared to experiments[END_REF] Here, only the elastic deformation is probed at the center of the sample (Δt = 0), which explains the observed plateau. We also note that Qi T presents a non-monotonous variation : Q1 T < Q3 T <Q2 T (also observed on Qi A , [26]. This probably results from a balance between the attenuation of the different waves, and the cumulative scattering detection each time the probe wave is transmitted through the shock (twice per each round trip).

Finally, introducing a time delay between shock and probe emissions allows to tune the crossing point between them inside the bulk of the material and therefor to monitor the shock propagation. Figure 5 shows a first experimental implementation. The ratios Qi T , i = 1 to 3, are plotted as a function of space (horizontal axis) and depth (vertical axis). The shock widening is clearly imaged. Shock wavefront is here planar, but more complex one, related for instance to phase-induced transformations in materials or shock focalization (caustics) could be investigated. To conclude, we designed a new and efficient way to monitor longitudinal shock wave propagation using a US probe. The method relies on the detection of the probe transmission through the shock. We assume the observed high amplitude of the phenomenon results from a singular interaction between acoustic and shock waves, predicted in the middle of the XXth century but never observed before in solids. A numerical model is currently under development to confirm this. Doppler effect could also be used as a complementary approach to quantify bulk shock properties. A broad range of applications could be explored using PA versatility, that allows to generate plane waves of different angles with respect to the shock direction, or focusing probes for instance. [START_REF] Drinkwater | Ultrasonic arrays for non-destructive evaluation: A review[END_REF] This time-space acoustic monitoring of shock propagation is intrinsically complementary to optics-based detection, from the visible to the X-Ray range. It could bring new understandings to numerous domains, from the physics of weak shocks (here considered) to warm dense matter investigations, including laser-driven shock waves used to simulate accretion or giant and high-velocity impacts of astrophysical bodies and which are currently be imaged on X-FEL facilities only. Other applications include the study of caustics in wave physics or

Figure 1 :

 1 Figure 1 : (a) Surface pressure induced by a ns laser pulse in aluminum. (b) Simulations of the density (dashed lines) and velocity (symbols) relative variations in a 2 mm thick plate, 300 ns (blue) and 450 ns (red) after generation.

Figure 2 :

 2 Figure 2: (a) Principles of experiments. (b) Space time representation : shock (plane wave) propagation is represented by the purple (blue) arrows and diffracted signals by dashed arrows. (c) Time-space PA detection of the pressure field. The dotted line rectangles illustrate the temporal integration on the acoustic echoes for the investigation procedure.

Figure 3 :

 3 Figure 3 : Detection of ratio Q1 A (a) and Q1 T (b) from single shot detection (dot blue lines) and averaging of 10 detections (red continuous line). On (b), the signal presented in black (with its own amplitude scale) is obtained from linear simulations of a static inclusion.

Figure 4 :

 4 Figure 4 : (a) Average detection of Qi T (i=1 (red), 2 (green), 3 (blue)). (b) FWHM of Qi T for a laser focus of 3 mm (black stars) and 5 mm (blue points) diameter. (c) Laser pressure loadings as function of the laser intensity from ESTHER simulations. (d) Intensity of the detected signal Qi T , i = 1 to 3, as function of the laser intensity.

Figure 4 (

 4 Figure 4(a) presents Qi T for up to three round trips in the material. One observes a widening with the number of acoustic trips, linked to the broadening of the shock beam along its propagation and to the directivity of the probe scattering on the shock edge.Figure 4(b) compares the spatial width of Qi T , i = 1 to 3, for a 3 mm

Figure 4 (

 4 b) compares the spatial width of Qi T , i = 1 to 3, for a 3 mm (asterisks) and a 5 mm (dots) laser diameter (noted D). The observed increase is almost linear with respect to the number of round trips for D = 3 mm. For D = 5 mm, the size is nearly constant during the first two acoustic travels and shows beyond a growth similar to the previous case. Considering the PA frequency, Rayleigh distance, Z = D²/2λ where λ is the acoustic wavelength, is closed to 17 (48) mm with a laser diameter of 3[START_REF] Pinton | Adaptive motion estimation of shear shock waves in soft solids and tissue with ultrasound[END_REF] mm. Therefore, in the D = 3 mm case, the successive detections are all performed in the diverging far-field regime, while in the second case the near to far field transition occurs during the second round trip, explaining the observed widening. This close link between detection, laser spot diameter and longitudinal shock broadening confirms the process is only due to the probe/shock interaction.

Figure 5 :

 5 Figure 5 : Elastography of the shock wave propagation during three round trips in the material.


Lucien Vienne

Emmanuel Lévêque

Recursive Finite-Difference Lattice Boltzmann Schemes

Keywords: Lattice Boltzmann method, recursive formulation, finite difference, von Neumann analysis

The motivation of this study is twofold. First, a recursive mathematical formulation of the discrete-velocity Boltzmann equation (DVBE) under the Bhatnagar-Gross-Krook (BGK) approximation is introduced. This formulation allows us to formally express the solution of the DVBE as an infinite sum over successive particle derivatives of the distributions associated with local equilibrium. A Chapman-Enskog multiplescales expansion shows that this sum can be safely truncated beyond the second order if the Navier-Stokes level of description is requested. Therefore, the distribution functions depend only on the first and secondorder derivatives of the related equilibrium distributions. This alternative equation to the DVBE defines a basis to design kinetic schemes for the evolution of the distribution functions based solely on flow variables that are sufficient to define local equilibrium. Second, a family of mass-conserving numerical schemes is introduced from this kinetic equation by discretizing the particle derivatives by backward finite differences. Interestingly, a so-called "simplified Lattice Boltzmann method" introduced by Chen et al. in 2017 can be recast in this family. Numerical simulations highlight a level of numerical dissipation that is generally higher than the level obtained with a standard Lattice Boltzmann scheme, as expected by approximating derivatives by finite differences. Nevertheless, we show by using a von Neumann analysis that it is possible to parametrize our scheme, according to the relaxation coefficient of the DVBE, to reduce significantly its numerical dissipation and improve its spectral properties. We believe that this modeling can also be of interest to connect macroscopic and kinetic representations, e.g. when dealing with initial and boundary conditions or in hybrid simulations matching Navier-Stokes and Lattice Boltzmann schemes.

Introduction

The Lattice Boltzmann (LB) method follows from a discretization of the Boltzmann equation to solve weakly-compressible fluid dynamics [START_REF] Shan | Discretization of the velocity space in the solution of the Boltzmann equation[END_REF][START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. Using a kinetic approach to simulate continuum flows seems at first glance unreasonable but, nevertheless, has certain advantages. At first, the velocity space is amenable to a radical decimation so that only a small set of microscopic velocities, e.g. nine in two dimensions, is sufficient to reconstruct the Navier-Stokes dynamics at the macroscopic level. In addition, non-locality and non-linearity are disentangled in the kinetic equations, which facilitates numerical integration [START_REF] Succi | Lattice boltzmann[END_REF]. As a result, the so-called stream-and-collide LB algorithm is simple, accurate and formidably efficient in terms of computations [START_REF] Körner | Parallel Lattice Boltzmann Methods for CFD Applications[END_REF].

In the LB approach, the degrees of freedom, or nodal values, refer to the distribution functions (f 0 , • • • , f N -1 ) of particles with the microscopic velocities (c 0 , • • • , c N -1 ). The scheme governs the evolution in time of these distributions at each lattice node. At the macroscopic level, flow variables are recovered by summing the contributions from the different (microscopic) velocities so that

ρ = N -1 α=0 f α and ρu = N -1 α=0 f α c α , (1)
where ρ and u denote respectively the mass density and the velocity of the fluid. Even if the decimation in velocity space is drastic, the number of nodal values remains significantly higher that the number of reconstructed flow variables, e.g. in two dimensions nine distribution functions are required to reconstruct locally the mass density and the two components of the velocity. The difference is even more pronounced in three-dimensions where nineteen or twenty seven densities are required. This overload of information is problematic when a given macroscopic solution needs to be prescribed at the mesoscopic level, e.g. for initial or boundary conditions. In that case, the kinetic solution is indeterminate and ad-hoc constraints must be invoked [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. From a computational viewpoint, the LB algorithm is data intensive and memory bound, which can be detrimental to its portability on accelerators such as Graphical Processing Units (GPUs). In this article, a LB algorithm that relies on flow variables only, and therefore requires less degrees of freedom, is introduced. It stems from the truncation and the discretization of a recursive formulation of the discretevelocity Boltzmann equation (DVBE). Some variants of the original stream-and-collide algorithm, depending on flow variables only, have already been proposed. In particular, Inamuro suggested a Bhatnagar-Gross-Krook (BGK) collision with a relaxation time (towards statistical equilibrium) equal to the time step of the algorithm. Therefore, distributions are considered at equilibrium and depend only on flow variables (ρ, u) [START_REF] Inamuro | A lattice kinetic scheme for incompressible viscous flows with heat transfer[END_REF]. This strong assumption gives a nonphysical value to the fluid viscosity that can then be corrected by adding a contribution of the viscous stress tensor estimated by finite differences. Asinari et al. derived a link-wise artificial compressibility scheme that can be formulated in the LB framework, with a collision operator that is non-local but depends on distributions at equilibrium only [START_REF] Asinari | Link-wise artificial compressibility method[END_REF]. Let us mention that the scheme developed by Inamuro can be viewed as a special case of the link-wise artificial compressibility scheme. More recently, Chen et al. suggested to consider the equilibrium state at the previous time step to build a prediction of the solution, which is afterwards corrected [START_REF] Chen | Three-dimensional simplified and unconditionally stable lattice boltzmann method for incompressible isothermal and thermal flows[END_REF] by considering a first-order approximation (in Knudsen number) of the DVBE. This prediction-correction scheme only involves macroscopic quantities. We shall see that this last scheme fits into the general framework developed in the present study.

A key ingredient behind formulating a LB scheme that relies only on flow variables is to express the nonequilibrium component of the distributions, f neq i ≡ f i -f eq i , as a function of the space-and-time derivatives of the equilibrium component f eq i (ρ, u). In this regard, Holdych et al. [START_REF] Holdych | Truncation error analysis of lattice boltzmann methods[END_REF] showed that f neq i could be rewritten in a recursive manner as

f neq α = ∞ m=1 (-τ ) m (∂ t + c α • ∇) m f eq α . (2)
Their motivation was to derive the truncation errors of this expansion and to show that some errors cancel out for specific values of the relaxation time τ . This iterative approach is similar to the so-called "Maxwell iteration" for the LB method [START_REF] Zhao | Maxwell iteration for the lattice boltzmann method with diffusive scaling[END_REF], in which the macroscopic equations are derived without explicitly resorting to a multiple-scale Chapman-Enskog expansion [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. Finally, it has also been used to tailor the equilibrium function of the standard LB scheme to solve non-linear equations such as the Burgers, Korteweg-de-Vries or Kuramoto-Sivashinsky equations by Otomo et al. [START_REF] Otomo | Two complementary lattice-boltzmann-based analyses for nonlinear systems[END_REF]. In the present article, our contribution is to use Eq. ( 2) to derive original LB schemes which rely on flow variables only. This is made possible by first showing that the sum can actually be truncated at m = 2 to comply with the Navier-Stokes dynamics at the macroscopic level. Then, numerical schemes can be designed by finite-difference discretization of the first and second-order particle derivatives. The remainder of the paper is organized as follows. Section 2 introduces the recursive formulation of the discrete-velocity Boltzmann equation and its discretization by finite differences. An original finite-difference LB scheme is proposed. It is also shown that the aforementioned prediction-correction scheme by Chen et al. [START_REF] Chen | Three-dimensional simplified and unconditionally stable lattice boltzmann method for incompressible isothermal and thermal flows[END_REF] can be derived in this framework. In section 3, a detailed comparison with the classical LB scheme is performed on the double shear layer and Taylor-Green vortex two-dimensional flows. Section 4 gives further insight on the proposed scheme by exploring spectral properties with a von-Neumann analysis. Finally, concluding remarks and perspectives are drawn in Section 5.

2. On a recursive formulation of the discrete-velocity Boltzmann equation

Context and derivation of the recursive formulation

In the LB approach, the sums replace the integrals in the statistical moments, as expressed in Eq. ( 1). This discretization of the velocity space stems from expanding and truncating the solution of the continuum Boltzmann equation onto a finite basis of Hermite polynomials in velocity, and resorting to a Gaussian quadrature formula to express the statistical moments [START_REF] Shan | Discretization of the velocity space in the solution of the Boltzmann equation[END_REF][START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the navier-stokes equation[END_REF]. Therefore, the discrete set of velocities may be thought of as the nodes in the Gaussian quadrature formula. It can then be established that the distributions f α (x, t) evolve according to a discrete-velocity analogue of the Boltzmann equation

(∂ t + c α • ∇) f α = - 1 τ (f α -f eq α ) , (3)
where the BGK approximation to the collision [START_REF] Bhatnagar | A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems[END_REF] is used in the right-hand side. This approximation expresses as a relaxation of the distributions towards their values at equilibrium, f eq i (ρ, u), with a unique relaxation time τ = ν/c 2 s fixed by the kinematic fluid viscosity ν and the speed of sound c s . This is sufficient to ensure that the slowly varying solution of the hierarchy of statistical equations built from Eq. ( 3) satisfy the isothermal Navier-Stokes equations in the low-Mach-number limit.

The set of microscopic velocities (c 0 , • • • , c N -1 ) determine the lattice so that the particles are streamed from one node to a neighbouring node during exactly one time step. In the present study, numerical results will refer to the two-dimensional D2Q9 lattice associated with the set of velocities (c 0 , • • • , c 8 ) given by the Cartesian components (see Fig. 1)

c α = 0 1 0 -1 0 1 -1 -1 1 0 0 1 0 -1 1 1 -1 -1 T ( 4
)
in lattice units 1 . The equilibrium distribution in Eq. ( 3) is an essential ingredient of the kinetic representation. In theory, it refers to a (continuous) Maxwell-Boltzmann distribution for classical fluids. In the LB framework, it is replaced by its Hermite expansion, which is commonly truncated at the second order. For our D2Q9 lattice, it yields

f eq α (ρ, u) = w α ρ 1 + c α • u c s 2 + (c α • u) 2 2c s 4 - u 2 2c s 2 (5)
where the weights w α are related to the lattice connectivity with w 0 = 4/9, w 1...4 = 1/9 and w 5...8 = 1/36 (see Fig. 1).

Let us now return to the discrete-velocity Boltzmann equation (3) which can be rewritten as

f α = f eq α -τ (∂ t + c α • ∇) f α . (6)
The recursive formulation arises by replacing f α in the rhs of Eq. ( 6) by its expression and by repeating this procedure infinitely. This leads formally to the expansion

f α = f eq α + ∞ m=1 (-τ ) m (∂ t + c α • ∇) m f eq α ( 7
)
in substitution of Eq. ( 3). This formulation shows that the f α 's can be expressed in terms of flow variables (ρ, u), however, this comes at the price of an infinite sum of increasing-order derivatives. The BGK collision operator is hidden in the form of this equation. Interestingly, the distribution function expresses (in lattice units) as an expansion in powers of the small parameter τ = τ /δ t ,

f α = f (eq) α (ρ, u) -τ (∂ t + c α • ∇) f (eq) α (ρ, u) + τ 2 (∂ t + c α • ∇) 2 f (eq) α (ρ, u) + . . .
At this stage, one might wonder whether this infinite sum could be truncated. We will demonstrate through a multiple-timescales Chapman-Enskog expansion that it is sufficient to consider the first two terms in the sum to comply with the Navier-Stokes equations at the macroscopic level.

Multiple-timescales Chapman-Enskog expansion

Macroscopic fluid dynamics follows from seeking solutions (to the hierarchy of statistical equations) that vary on a much slower timescale than the collisional timescale τ . This is usually done by using a multiple-timescales Chapman-Enskog expansion in the small parameter , with

f α = f (0) α + f (1) α + + 2 f (2) α + • • • ∂ t = ∂ t1 + 2 ∂ t2 + • • • ∇ = ∇ 1 (8)
The parameter may be identified physically with the Knudsen number and, t 1 and t 2 with advective and diffusive timescales. Substituting these expansions into our recursive formulation, we obtain at successive orders in O( 0 ) :

f (0) α = f eq α (9) O( 1 ) : f (1) α = -τ (∂ t1 + c α • ∇ 1 ) f (0) α (10) O( 2 ) : f (2) α = -τ ∂ t2 f (0) α + τ 2 (∂ t1 + c α • ∇ 1 ) 2 f (0) α = -τ ∂ t2 f (0) α -τ (∂ t1 + c α • ∇ 1 ) f (1) α (11)
Consistently, this set of equations (up to O( 2)) is identical to that obtained by applying the Chapman-Enskog expansion to the original discrete-velocity Boltzmann equation [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. This set is sufficient to obtain the continuity equation and the Navier-Stokes equations by taking its statistical moments α (•) and α c α (•), as briefly recalled in Appendix A. An essential point is that Eqs. ( 9), [START_REF] Otomo | Two complementary lattice-boltzmann-based analyses for nonlinear systems[END_REF] and [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the navier-stokes equation[END_REF] come only from the first two terms of the recursive formulation, and that subsequent terms in the sum lead to higher orders in . This is evidenced by reformulating Eq. [START_REF] Chen | Three-dimensional simplified and unconditionally stable lattice boltzmann method for incompressible isothermal and thermal flows[END_REF] as

f α = f eq α + ∞ m=1 (-τ ) m (∂ t1 + c α • ∇ 1 + ∂ t2 + • • • ) m f eq α ( 12
)
where the m-th term in the sum contributes at least to O( m ) in the multiple-scale expansion2 . Therefore, it is enough to consider the first two terms to comply with the Navier-Stokes dynamics at the macroscopic level, i.e.

f α = f eq α -τ (∂ t + c α • ∇) f eq α + τ 2 (∂ t + c α • ∇) 2 f eq α . (13)
In the following, this second-order truncation will be the cornerstone of our numerical developments.

In practice, the LB scheme is obtained by integrating the discrete-velocity Boltzmann equation along the characteristics during a time interval δ t , and by approximating the integral of the collision term by a trapezium rule [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. A convenient change of variables for the distribution functions is usually operated to render the scheme fully explicit. Namely,

g α = f α + δ t 2τ (f α -f eq α ) ( 14
)
where the g α 's define a new set of distribution functions linked to the lattice discretization. This implies in particular that g eq α = f eq α and g neq α ≡ g α -g eq α = (τ g /τ )f neq α with the modified relaxation time

τ g = τ + δ t 2 . ( 15
)
In these new variables, the recursive formulation (restricted to second order) writes equivalently as

g α (x, t) = g eq α (x, t) -τ g (∂ t + c α • ∇) g eq α (x, t) + τ g (τ g - δ t 2 ) (∂ t + c α • ∇) 2 g eq α (x, t). (16)
In this article, the developments will be carried out with f α and τ but could be easily transposed to g α and τ g . The final step towards a numerical scheme is the discretization of space and time variables. In this respect, Eq. ( 13) offers a large panel of possibilities. Since the equation implies particle derivatives, a natural choice (taken later in this study) is to integrate it along particle trajectories and, therefore, resort to an underlying uniform lattice structure, in the spirit of the standard lattice Boltzmann method. However, time and space variables can also be disentangled by treating independently the temporal and spatial derivatives, leading to so-called off-lattice Boltzmann schemes [START_REF] Nannelli | The lattice boltzmann equation on irregular lattices[END_REF]. These latter allow for grid cells with arbitrary shapes and thus provide a higher flexibility in the mesh geometry. Recent studies (starting from the DVBE) reveal the attractive potential of such schemes, which have now overcome the accuracy and stability issues initially encountered [START_REF] Bardow | General characteristic-based algorithm for off-lattice boltzmann simulations[END_REF][START_REF] Patil | Chapman-enskog analysis for finite-volume formulation of lattice boltzmann equation[END_REF][START_REF] Zhu | Performance evaluation of the general characteristics based off-lattice boltzmann scheme and dugks for low speed continuum flows[END_REF]. The additional expense in computational costs is offset by the gain in the geometrical flexibility [START_REF] Di Ilio | Simulation of turbulent flows with the entropic multirelaxation time lattice boltzmann method on body-fitted meshes[END_REF][START_REF] Shu | Development of lattice boltzmann flux solver for simulation of incompressible flows[END_REF]. Combining several discretization methods is also possible to take advantage of the high efficiency (and low dissipation rate) of the stream-and-collide algorithm and the geometrical flexibility of off-lattice discretization for example, as recently demonstrated in [START_REF] Di Ilio | Hybrid lattice boltzmann method on overlapping grids[END_REF][START_REF] Horstmann | Hybrid simulation combining two space-time discretization of the discrete-velocity boltzmann equation[END_REF].

Space-and-time discretization by backward finite differences

At first, it would be tempting to recast Eq. (13) into a system of two coupled first-order differential equations reminiscent of the original discrete-velocity Boltzmann equation. Such system would write

(∂ t + c α • ∇) f eq α = h α ( 17
) (∂ t + c α • ∇) h α = 1 τ 2 (f α -f eq α ) + 1 τ h α . (18)
Continuing the analogy, one could integrate the left-hand side along the characteristics and approximate the right-hand side by a trapezium rule. However, this would ultimately lead to a system of two coupled streamand-collide algorithms with more degrees of freedom (two sets of distribution functions) than the original LB algorithm. Moreover, our preliminary simulations have shown that this double-distribution scheme is unstable.

It is rather proposed here to discretize the particle derivatives in Eq. ( 13) by backward finite-differences in time, namely

(∂ t + c α • ∇) f eq[0] α = 1 2δ t 3f eq[0] α -4f eq[-1] α + f eq[-2] α + O δ t 2 (19)
(∂ t + c α • ∇) 2 f eq[0] α = 1 δ 2 t f eq[0] α -2f eq[-1] α + f eq[-2] α + O (δ t ) (20)
where the short notation f

eq[i] α
stands for f eq α (x + ic α δ t , t + iδ t ). This results in the formula

f neq[0] α = τ (τ - 3 2 )f eq[0] α + 2(1 -τ )f eq[-1] α + (τ - 1 2 )f eq[-2] α with τ = τ /δ t . (21)
This equation allows us to update the macroscopic variables ρ(x, t) and u(x, t) from the values at neighbouring nodes at the previous time steps, by noting importantly that the non-equilibrium distributions f

neq[0] α must satisfy the solvability conditions 3 α f neq[0] α = 0 and α f neq[0] α c α = 0 (22)
and that α f

eq[0] α = ρ(x, t) and α f eq[0] α c α = (ρu)(x, t) by construction.
It is a priori desirable to have the same order of accuracy for the first and second differential operators. To do this, a four-point backward formula should be considered for the second-order operator

(∂ t + c α • ∇) 2 f eq[0] α = 1 δ 2 t 2f eq[0] α -5f eq[-1] α + 4f eq[-2] α -f eq[-3] α + O δ t 2 (23)
which this time yields

f neq[0] α = τ (2τ - 3 2 )f eq[0] α + (2 -5τ )f eq[-1] α + (4τ - 1 2 )f eq[-2] α -τ f eq[-3] α . (24)
In order to bridge the two discretizations ( 21) and [START_REF] Chen | On numerical diffusion of simplified lattice boltzmann method[END_REF], it is here proposed to leave the weight related to f eq[-3] α as a free parameter γ. After manipulating Taylor expansions, this leads to the parametric family of finite-difference schemes

f neq[0] α = τ (-γ + τ - 3 2 )f eq[0] α + (3γ + 2(1 -τ ))f eq[-1] α + (-3γ + τ - 1 2 )f eq[-2] α + γf eq[-3] α (25)
consistent with the original differential equation [START_REF] Nannelli | The lattice boltzmann equation on irregular lattices[END_REF]. One must mention that the sum of the coefficients in Eq. ( 25) being zero ensures that this scheme is mass-conserving. Eq. ( 25) also shares some similarity with the multi-step LB methods [START_REF] Wilde | Multistep lattice boltzmann methods: Theory and applications[END_REF], in which the distribution functions are expressed according to their value at previous times and neighbouring nodes. If γ = 0, the first and second differential operators are approximated respectively at second and first order and Eq. ( 21) is obtained. If γ = -τ , the first and second derivatives are approximated with a second-order accuracy and Eq. ( 24) is obtained. If γ = -τ + 1 3 , the first and second derivatives are approximated respectively with a third-order and second-order accuracy. Interestingly, we will see later that the parameter γ can be optimized to significantly reduce the discretization errors.

At this point, we can mention that the regularized LB method introduced by Latt & Chopard [START_REF] Latt | Lattice boltzmann method with regularized pre-collision distribution functions[END_REF] resembles our modeling in that both methods conform to the Navier-Stokes equations at the macroscopic level, and provide a reconstruction of the non-equilibrium component of the distributions by an alternative means to the standard LB scheme. However, the regularized LB is a local reconstruction (in (x, t)) that can be viewed as a projection on a restricted basis of Hermite polynomials, while our method is a nonlocal reconstruction based on the equilibrium distributions. If the source of numerical dissipation may be attributed to the projection operation (in the velocity space) in the regularized method, it will be rather related to the spatio-temporal discretization in our case. Concerning the link-wise artificial compressibility (LW-AC) method, it shares with the present modeling its finite-difference formulation. However, the asymptotic analysis of this model (developed in appendix B in [START_REF] Asinari | Link-wise artificial compressibility method[END_REF]) does not yield the high-order temporal and mixed derivatives that we obtain through the recursion. Therefore, the two methods cannot infer from each other. We suspect that the "empirical" substitution of the distributions by their equilibrium components (Eqs. [START_REF] Inamuro | A lattice kinetic scheme for incompressible viscous flows with heat transfer[END_REF][START_REF] Otomo | Two complementary lattice-boltzmann-based analyses for nonlinear systems[END_REF] in [START_REF] Asinari | Link-wise artificial compressibility method[END_REF]) is the main factor behind this difference.

From Eq. ( 25) and by resorting to the solvability conditions [START_REF] Latt | Lattice boltzmann method with regularized pre-collision distribution functions[END_REF], the flow variables are updated according to the sums

ρ (x, t) = α f * α , (26)
(ρu) (x, t) = α f * α c α (27)
with the effective distributions f * α given by

f * α = 3γ + 2(1 -τ ) γ -τ + 3 2 f eq[-1] α + -3γ + τ -1 2 γ -τ + 3 2 f eq[-2] α + γ γ -τ + 3 2 f eq[-3] α . (28)
Since, f eq α depends on flow variables only (see Eq. ( 5)), f * α can be evaluated from the previous values of ρ and u at neighbouring lattice nodes and previous time steps. It can be pointed out that the schemes not only advance in time the density and the velocity fields but also provide a direct estimation of the rate-ofstrain tensor through the second-order non-equilibrium moment: -2ρτ c 2 s S ≈ α f neq α c α c α . The first steps of the algorithm can be performed by using the standard (one-step) lattice Boltzmann scheme. The A-B memory access pattern, commonly referred to as "ping-pong buffering", is typically used in the LB scheme for better performance [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. It requires two sets of distributions at each node, and thus 18 and 38 values for the D2Q9 and D3Q19 lattices. These two numbers are reduced to 12 and 16 in the present modelling (with γ = 0). In terms of computations, our algorithm requires to upload the three previous values of ρ and u at neighbouring positions to reconstruct f In the following, our scheme will be referred to as recursive finite-difference lattice Boltzmann scheme with reference to the recursive formulation of the original DVBE.

Simplified prediction-correction Lattice Boltzmann algorithm

Interestingly, a so-called simplified LB method relying solely on the equilibrium distributions has recently been introduced by Chen et al. in [START_REF] Chen | Three-dimensional simplified and unconditionally stable lattice boltzmann method for incompressible isothermal and thermal flows[END_REF] and has proven to be relevant for many applications [START_REF] Chen | Simplified lattice boltzmann method for non-newtonian power-law fluid flows[END_REF][START_REF] Chen | On numerical diffusion of simplified lattice boltzmann method[END_REF][START_REF] Lu | The more actual macroscopic equations recovered from lattice boltzmann equation and their applications[END_REF]. In the following, we shall see that this method may be recast in the present framework by resorting to different levels of approximation in the recursive formulation of the discrete-velocity Boltzmann equation.

The simplified LB method proceeds in a two-step prediction-correction procedure. The density and the fluid momentum are first predicted as

ρ * (x, t) = α f eq[-1] α , ( 29
) (ρu) * (x, t) = α c α f eq[-1] α ( 30
)
which is equivalent to considering a recursive expansion restricted to m = 1 and a first-order backward finite-difference approximation of the derivative, i.e.

f neq[0] α -τ f eq[0] α -f eq[-1] α . (31)
Obviously, this dynamics does not conform to the Navier-Stokes equations and must be corrected. A Taylor expansion in time of

f eq[-1] α ≡ f eq α (x -c α δ t , t -δ t ) gives f eq[-1] α = f eq[0] α -δ t (∂ t + c α • ∇) f eq[0] α + δ t 2 2 (∂ t + c α • ∇) 2 f eq[0] α + o δ t 2 , (32)
which may be reformulated as

f eq[-1] α = f eq[0] α -δ t 1 - 1 2τ (∂ t + c α • ∇) f eq[0] α + δ 2 t 2τ 2 f neq[0] α + o δ t 2 (33)
by expressing the second-order derivative from the recursive expansion [START_REF] Nannelli | The lattice boltzmann equation on irregular lattices[END_REF]. The first and second moments then yield at (x, t)

ρ * = ρ -δ t 1 - 1 2τ (∂ t ρ + ∇ • (ρu)) + o δ t 2 , (34)
(ρu) * = ρu -δ t 1 - 1 2τ (∂ t (ρu) + ∇ • Π eq ) + o δ t 2 . ( 35
)
Therefore, no correction is needed for the density, i.e. ρ(x, t) = ρ * (x, t), to comply with the continuity equation, whereas the momentum must be corrected as

ρu = (ρu) * -δ t (1 - 1 2τ )∇ • Π neq + o δ t 2 with Π neq = α f neq α c α c α (36)
to be consistent with the Navier-Stokes dynamics. The divergence term needs to be discretized to complete the model. To do so, a Taylor expansion in space of f neq α (x, t) gives

f neq α x ± 1 2 c α δ t , t = f neq α (x, t) ± 1 2 δ t (c α • ∇)f neq α (x, t) + 1 4 δ 2 t (c α • ∇) 2 f neq α (x, t) + o δ t 2 (37)
which allows us to finally rewrite the correction for the momentum as

(ρu)(x, t) = (ρu) * (x, t) -(1 - 1 2τ ) α c α f neq α x + 1 2 c α δ t , t -f neq α x - 1 2 c α δ t , t + o δ t 2 . ( 38
)
The approach is so far general and results from manipulating the discrete-velocity Boltzmann equation and using Taylor expansions at second order. Interestingly, Chen et al. [START_REF] Chen | Three-dimensional simplified and unconditionally stable lattice boltzmann method for incompressible isothermal and thermal flows[END_REF] then proposed to simplify Eq. ( 38) by evaluating the correction at intermediate time t-δt/2 and by approximating the non-equilibrium distributions as

f neq α x + 1 2 c α δ t , t - δ t 2 -τ (f eq * α (x + c α δ t , t) -f eq α (x, t -δ t )) (39)
f neq α x - 1 2 c α δ t , t - δ t 2 -τ (f eq * α (x, t) -f eq α (x -c α δ t , t -δ t )) . (40)
This is equivalent to considering a recursive expansion restricted to m = 1 and a second-order central finitedifference approximation of the particle derivative. As expected in a prediction-correction procedure, the distributions at time t are evaluated from the predicted density and momentum, i.e. f eq * α = f eq α (ρ * , u * ). The resulting algorithm is indeed simple and may be summarized as

• prediction step: ρ * (x, t) = α f eq[-1] α (41) (ρu) * (x, t) = α c α f eq[-1] α (42) • correction step: ρ(x, t) = ρ * (x, t) (43) (ρu)(x, t) = (ρu) * (x, t) -(τ - 1 2 ) (ρu)(x, t -δ t ) - α c α f eq α (ρ * , u * ) (x + c α δ t , t) (44)
Therefore, this simplified LB scheme can be recast in the present modeling framework (under specific assumptions). It will be considered for comparisons in the following numerical tests.

Numerical simulations

In this section, two flow configurations are considered to assess the validity of our recursive finitedifference scheme (abbreviated as rfd, Eqs. ( 26), [START_REF] Coreixas | Impact of collision models on the physical properties and the stability of lattice boltzmann methods[END_REF], and ( 28)) against the prediction-correction scheme (precorr, Eqs. (41), (42), (43), and (44)) and the standard stream-and-collide LB algorithm (lbm, Eq. ( 45)). This latter expresses as

g α (x + c α δ t , t + δ t ) = g α (x, t) - 1 τg (g α (x, t) -g eq α (x, t)) (45)
with τg = τ + 1/2 and g eq α = f eq α . The rfd algorithm is initialized by performing the first two steps with the standard LB method. Two different values of the adjustable parameter γ are considered in the analysis, namely γ = 0 and an optimal value obtained by a trial-and-error procedure (detailed later). To improve the stability of the LB scheme, the equilibrium distributions have been developed at fourth order in Hermite expansion in the various schemes [START_REF] Chen | Fundamental conditions for n-th-order accurate lattice boltzmann models[END_REF][START_REF] Coreixas | Impact of collision models on the physical properties and the stability of lattice boltzmann methods[END_REF], i.e.

f eq α (ρ, u) = w α ρ 1 + c α • u c s 2 + (c α • u) 2 2c s 4 - u 2 2c s 2 + (c α • u) 3 6c s 6 - (c α • u)u 2 2c s 4 + (c α • u) 4 24c s 8 - (c α • u) 2 u 2 4c s 6 + u 4 8c s 4 (46)
In the remainder of the article, all quantities are expressed in lattice units obtained by normalizing space and time by the lattice spacing and the time step, respectively.

Double shear-layer flow

The double shear-layer flow is a classical test case for numerical methods. Two thin horizontal shear layers are located at y = L/4 and y = 3L/4 in a two-dimensional periodic domain of size [0, L] × [0, L]. A small vertical velocity perturbation is added in order to disturb the flow. A Kelvin-Helmholtz instability is expected to develop with two counter-rotating growing vortices. The thinning of the shear layers can lead to the development of spurious vortices and possibly the blow-up of the simulation.

The initial velocity field is

u x (x, 0) = U 0 tanh k y L - 1 4
for y ≤ L/2 (47)

u x (x, 0) = U 0 tanh k 3 4 - y L for y ≥ L/2 (48) u y (x, 0) = U 0 δ sin 2π x L + 1 4 ( 49
)
where U 0 is the characteristic velocity of the flow, k/L is inversely proportional to the thickness of the two shear layers and δ is the amplitude coefficient of the initial perturbation in velocity. Following previous works, the double shear-layer flow is simulated with k = 80 and δ = 5% [START_REF] Dellar | Bulk and shear viscosities in lattice boltzmann equations[END_REF][START_REF] Coreixas | Recursive regularization step for high-order lattice boltzmann methods[END_REF]. The reference velocity U 0 is fixed so that the Mach number is Ma = U 0 /c s = √ 3U 0 = 0.3. The length L corresponds to the resolution of the simulation (in lattice units) and is considered as a parameter. The Reynolds number is kept constant at Re = U 0 L/ν = 3 • 10 4 by adjusting the viscosity. This configuration corresponds to a challenging test case that can become unstable at low resolutions. Finally, let us mention that the initial density (or pressure) field is not known. To circumvent this difficulty, the iterative procedure proposed by Mei et al. [START_REF] Mei | Consistent initial conditions for lattice boltzmann simulations[END_REF] was used to generate the distribution functions in a consistent manner and avoid transient errors; the initial density was set to unity and the stopping error criterion was taken to be 10 -8 in the iterative procedure.

In Fig. 3, the vorticity for various resolutions is displayed at the same convective time t * = 1 with t * ≡ t U 0 /L for the different algorithms. For the rfd algorithm, the top shear layer refers to a simulation with γ = 0, whereas the bottom shear layer is simulated with an optimized γ = 0 obtained by a trial-anderror method. An analytical approach (based on a von-Neumann analysis) will be introduced later in the article to infer the optimal γ. For the moment, the optimal γ is found by incremental testing, i.e. by trying successive values of γ with an increment of 0.01. The optimal γ minimizes the error made on the averaged kinetic energy with respect to a reference solution in the interval 0 ≤ t * ≤ 2. In practice, we used the solution of the lbm simulation at the highest resolution (L = 512) as reference.

At low resolution (L = 128) the lbm simulation is unstable whereas the rfd and precorr simulations remain stable but the counter-rotating vortices are obviously over-damped. At intermediate resolution (L = 256) the lbm gains in stability but spurious vortices still appear where the shear layer is the thinnest. These vortices are not observed in the rfd and precorr simulations, however, the latter produce a larger numerical dissipation as evidenced by the diffusion of the shear layer. Interestingly, the situation can be greatly improved by selecting the optimal parameter γ = 0.18 in the rfd algorithm. Finally, all schemes give satisfactory results and converge towards the same solution at high resolution (L = 512) except for the rfd with γ = 0, for which the vortex is still too dissipated.

In summary, these first qualitative results show that the rfd and precorr algorithms allow us to simulate the double-shear-layer flow (at all resolutions) but this comes with an accentuated numerical dissipation as compared to the original LB scheme. Nevertheless, it is possible to reduce significantly the numerical dissipation in the rfd simulation by adjusting the value of the parameter γ. This is quantitatively confirmed in Fig. 4, where the mean kinetic energy is plotted at the intermediate resolution for the different algorithms. At this resolution, the rfd simulation with the optimized γ is stable and very close to the reference lbm solution.

Taylor-Green vortex flow

The two-dimensional Taylor-Green vortex flow is another meaningful test case as it provides an exact analytical solution for the Navier-Stokes equations (at low Reynolds number). The solution reads for a fully

x (x, t) = -U 0 k y k x cos(k x x) sin(k y y) exp(-t/T ν ), (50)
u y (x, t) = U 0 k x k y sin(k x x) cos(k y y) exp(-t/T ν ), (51)
ρ(x, t) = 1 - U 2 0 4c s 2 k y k x cos(2k x x) + k x k y cos(2k y y) exp(-2t/T ν ) ( 52
)
where U 0 is the initial characteristic velocity, k x = k y = 2π/L, and

T ν = 1/[(k 2 x + k 2 y )ν]
is the characteristic decay time. The initial distribution functions for the lbm algorithm are taken at equilibrium value, i.e. f α (x, 0) = f eq α (ρ(x, 0), u(x, 0)); here both ρ(x, 0) and u(x, 0) are known. At low Reynolds number, the initial vortices decay in time with an exponential rate T ν due to the viscous dissipation. As before, various simulations are performed with the rfd, precorr and lbm schemes and compared with the analytical solution. The relative error is quantified by using the L 2 -norm over the entire computational domain, i.e.

ε ϕ (t * ) = x (ϕ simulated (x, t * ) -ϕ analytical (x, t * )) 2 x ϕ 2 analytical (x, t * ) (53)
for the variables ϕ = ρ, u x and u y at the normalized diffusive time t * = t/T ν . In practice, only the error in ε ux will be examined, the error ε uy being of the same order of magnitude. At first, different domain sizes (or resolutions) are considered, namely L = 32, 64, 128 and 256, and a convergence analysis is carried out. Here, it is appropriate to adopt a diffusive scaling, meaning that δ x 2 /δ t is kept constant through all simulations. Consequently, the relaxation time (or the viscosity) is constant in lattice units but the Mach number varies inversely to the resolution. The similitude in the Mach number can be disregarded as long as acoustic features are ignored, i.e. the Mach number remains low. Such condition is met here. The characteristic velocity is given by

U 0 = Re ν L . (54)
The Taylor-Green vortex flow is simulated up to two decay times (t * = 2) at Re = 100 with ν = 0.01. The Mach number varies as Ma = √ 3/L 1 with the resolution L. The relative errors are evaluated at t * = 2 and shown in Fig. 5. For all schemes, a fourth and second-order convergence are found for the density and the velocity respectively. This is consistent with the second-order accuracy in space of the different algorithms. The difference of order between the density and the velocity agrees with the hypothesis that typically δρ ∼ ρ 0 u 2 /c 2 s in the weak-compressibility limit. Nevertheless, while sharing the same order of convergence, the initial error for the rfd algorithm with γ = 0 and the precorr algorithm is much greater. The optimal value γ = 0.15 in the rfd algorithm reduces this initial error at a level close to that of the lbm simulation. No variation of the optimal γ was found for the different resolutions with a constant relaxation time. The influence of the relaxation time, or viscosity, at constant Reynolds number Re = 100 and resolution L = 32 is now examined. Changing the relaxation time (τ = 3ν in lattice units) leads to a variable Mach number with Ma = (Re/ √ 3L)τ . It is observed in Fig. 6a that the error in density naturally increases with τ . Large errors on the velocity are observed in Fig. 6b for small relaxation times as compared to the lbm algorithm. The numerical dissipation of the rfd algorithm with γ = 0 and the precorr algorithm damages the numerical accuracy at this low resolution (L = 32). At the smallest viscosity, the numerical dissipation inherent to the finite-difference discretization exceeds the actual dissipation and the optimization of γ no longer helps to reduce ε ux in the rfd algorithm. This obviously points to an obvious limitation of this type of approach in the limit of vanishing viscosity, unless the resolution is considerably increased. On the other hand, at high viscosity the rfd algorithm with an optimized γ appears to provide the most accurate value for the velocity, the lbm solution being degraded by errors in Mach number. Such errors, which are a priori also present in the finite-difference approach, can be "compensated" by picking the optimal γ. Finally, let us notice that the optimal γ depends on τ . In the next configuration, both the viscosity (or relaxation time) and the resolution are kept constant with ν = 0.01 and L = 32. The impact of the characteristic velocity (or Mach number) alone is shown in Fig. 7. The relative error in density typically increases as Ma 2 whereas the velocity error remains mostly identical. For U 0 = 0.3125, the Mach number reaches Ma = 0.5 which is detrimental to the accuracy of the lbm algorithm. The value of the optimal γ remains almost constant. Finally, the evolution in time of the error is plotted for L = 32, Re = 100 and ν = 0.01 in Fig. 8. The relative error in density improves in time as the level of density fluctuations decreases. On the other hand, the velocity error deteriorates with time for the rfd (γ = 0) and precorr algorithms as spurious numerical dissipation accumulates. Corrections to limit this over-dissipation for the precorr scheme (as suggested in [START_REF] Chen | On numerical diffusion of simplified lattice boltzmann method[END_REF]) or γ = 0 for the rfd scheme appears to be mandatory. In summary, we first confirm that the precorr algorithm gives very satisfactory results in view of its simplicity, namely, its macroscopic character, its locality (restricted to first neighbours) and the absence of any adjustable parameter; it performs better than the rfd algorithm with γ = 0. However, adjusting γ in the rfd scheme provides a way to significantly reduce the numerical dissipation, and obtain results comparable with those given by lbm scheme, except for very low relaxation times (in lattice units) where the lbm remains superior in terms of accuracy. The selection of γ in the rfd scheme is important, however, a prediction of the optimal γ would be desirable to avoid any prior trial-and-error procedure. In the next section, we will see by using a von Neumann analysis that it is possible to infer (to a certain extent) the optimal value of the rfd scheme as a function of the relaxation coefficient τ .

Von Neumann analysis of the recursive finite-difference lattice Boltzmann scheme

Comparisons

The so-called von Neumann (stability) analysis is a procedure that is commonly used to investigate the stability and Fourier spectral properties of finite-difference schemes [START_REF] Hirsch | Chapter 7 -consistency, stability and error analysis of numerical schemes[END_REF]. It is here used to gain insight about the calibration of the parameter γ in the rfd scheme. At first, the distributions f α (x, t) are considered as small perturbations around a stationary and uniform reference state (ρ, U ): f α = fα + f α with fα = f eq α (ρ, U ). The fluctuating components are expressed as complex monochromatic plane waves

f α = A α exp (i(k • x -ωt)) (55)
where ω and k denote respectively the (complex) pulsation and the wavevector, A α is the amplitude of the wave. This decomposition is introduced in our scheme

f neq[0] α = τ (-γ + τ - 3 2 )f eq[0] α + (3γ + 2(1 -τ ))f eq[-1] α + (-3γ + τ - 1 2 )f eq[-2] α + γf eq[-3] α (56)
by assuming that f eq α fα + J αβ f β with the Jacobian matrix J αβ = (∂f eq α /∂f β ) being evaluated at the reference state. This eventually yields for each α

(δ αβ -J αβ )f β = τ (-γ + τ - 3 2 ) + (3γ + 2(1 -τ )) exp (i(ω -k • c α )δ t ) + (-3γ + τ - 1 2 ) exp (2i(ω -k • c α )δ t ) + γ exp (3i(ω -k • c α )δ t ) J αβ f β (57)
with an implicit summation over the index β. This dispersion equation may be recast into a generalized eigenvalue problem by first rewriting it in a matrix form as (Aλ 3 +Bλ 2 +Cλ+D)f = 0 with λ = exp (iωδ t ) and by introducing the new vectors g = λf and h = λg . Precisely, one obtains

  -D 0 0 0 I 0 0 0 I     f g h   = λ   C B A I 0 0 0 I 0     f g h   (58)
where I stands for a 9 × 9 identity matrix (in the case of a D2Q9 lattice) and

A αβ = τ γ exp (-3iφ α )J αβ , B αβ = τ (-3γ + τ -1 2 ) exp (-2iφ α )J αβ , C αβ = τ (3γ + 2(1 -τ )) exp (-iφ α )J αβ and D αβ = -δ αβ + (1 + τ (-γ + τ -3 2 ))J αβ with φ α = k • c α δ t .
The non-locality of the finite-difference stencil increases the size of the eigenvalue problem. The Jacobian of the equilibrium distributions has been derived analytically with the SymPy Python module for symbolic computation. The eigenvalue problem has been solved by using the SciPy Python module. The pulsations are eventually given by

ωδ t = -i log(λ) (59)
for a given reference state (ρ, U ) at wavevector k. By noting that exp(-iωt) = exp( (ω)t) • exp(-i (ω)t), it follows that the dispersion of the wave is associated with (ω) whereas the dissipation is linked to (ω).

To improve the readability of the plots, over-damped modes ( (ω α )/ν < -10 4 ) have been disregarded. The dispersion and dissipation properties have been compared with those of the standard LB scheme, for which the von Neumann analysis leads to the eigenvalue problem

f = λA I - 1 τg (I -J ) f with A αβ = exp (-ik • c α δ t ) δ αβ (60)
instead of Eq. ( 58). The pulsations obtained analytically from the isothermal Navier-Stokes equations [START_REF] Marié | Comparison between lattice boltzmann method and navier-stokes high order schemes for computational aeroacoustics[END_REF] will also be indicated in the figures. In Fig. [START_REF] Zhao | Maxwell iteration for the lattice boltzmann method with diffusive scaling[END_REF] and Fig. [START_REF] Otomo | Two complementary lattice-boltzmann-based analyses for nonlinear systems[END_REF], the dispersion and dissipation rates are displayed as a function of the wavenumber for the standard lbm scheme, the Navier-Stokes equations and the proposed rfd scheme with different values of γ. The lattice units are used, i.e. δ x = 1 and δ t = 1. The Mach number and the relaxation time are the two independent parameters of the analysis. The Mach number is fixed at Ma = 0.1 considered as a representative value. On the other hand, two relaxation times τ = 3 • 10 -3 (Fig. 9) and τ = 3 • 10 -1 (Fig. 10) are examined. In both cases, the reference velocity is horizontal and the analysis is restricted to waves propagating along the x-axis for simplicity. The wavenumber varies between 0 and π in lattice units.

First, we can notice that the LB approach, which treats the fluid at an underlying level, involves more degrees of freedom than the macroscopic Navier-Stokes approach and, consequently, more modes in the von Neumann analysis. It is interesting to understand to what extent the "macroscopic information" is carried by the mesoscopic modes [START_REF] Wissocq | An extended spectral analysis of the lattice boltzmann method: modal interactions and stability issues[END_REF]. In our case, the attention will focus on the dissipation properties, namely, how (ω) compares with the theoretical prediction -νk 2 at wavenumber k. For τ = 3 10 -3 in Fig. 9, increasing the accuracy of the second-order derivative from first (γ = 0) to second-order accuracy (γ = -τ ) has a negligible effect on the spectral properties. For γ = -τ + 1/3, the dispersion behaviour of the shear mode is greatly improved for k x > π/2, however, (ω) takes positive values involving an exponential growth detrimental to the stability. Therefore, improving the accuracy does not a priori help. For all γ, we observe that the level of dissipation of the rfd scheme is much higher than that of the lbm for k x > π/8, which agrees with the observations made in the test simulations of the previous section. In Fig. 10, the relaxation time is increased at τ = 0.3. According to τ = √ 3ν/c s δx, this can be achieved by decreasing δ x by a factor of 100 with all physical parameters kept constant. In that case, it is observed that the spectral behaviour of the rfd scheme is greatly improved and very comparable to the spectral behaviour of the lbm scheme for the dissipation, which is remarkable for a finite-difference scheme. Nevertheless, this comes with the price of increasing significantly the computational cost by reducing the grid spacing δ x . Interestingly, a non-trivial degraded behaviour is observed for γ = -τ even if this value corresponds to a higher order in accuracy than γ = 0. The optimization of γ is therefore not straightforward.

The strategy of optimizing the spectral properties of finite-difference schemes has already been widely used in computational acoustics with the so-called dispersion-relation-preserving schemes [START_REF] Tam | Dispersion-relation-preserving finite difference schemes for computational acoustics[END_REF]. In the lattice Boltzmann domain, the stability of the multiple-relaxation-time collision operator is optimized by carefully choosing the relaxation rates of moments through a von Neumann analysis [START_REF] Lallemand | Theory of the lattice boltzmann method: Dispersion, dissipation, isotropy, galilean invariance, and stability[END_REF][START_REF] Xu | Optimal low-dispersion low-dissipation LBM schemes for computational aeroacoustics[END_REF][START_REF] Chávez-Modena | Improving the stability of multiple-relaxation lattice boltzmann methods with central moments[END_REF]. The impact of various collision models on the linear stability has been recently reviewed in [START_REF] Coreixas | Impact of collision models on the physical properties and the stability of lattice boltzmann methods[END_REF]. In [START_REF] Krivovichev | Parametric schemes for the simulation of the advection process in finite-difference-based singlerelaxation-time lattice boltzmann methods[END_REF], parametric finitedifference schemes for the Boltzmann equation are proposed, and by selecting a well-suited parameter, the spurious numerical dissipation can be reduced. In the same vein, analytical developments are carried out in the next section to infer an optimal value of γ with respect to the dissipation dynamics.

Predicting the parameter γ

The trial-and-error method used previously to find the optimal parameter γ is efficient but can not be considered as an end in itself. A more straightforward way is desirable. In this subsection, we try to take . Only the dependence along the x-axis is here considered. For the analytical Navier-Stokes solution, the different waves correspond to two acoustic modes, ω ac± propagating at speed U ± cs and a shear mode, ωs,propagating at speed U ; the acoustic and shear modes are dissipated at rate -νk 2 . In the rfd, the value γ = 0 corresponds to an approximation of the first and second derivatives with second-order and first-order accuracy respectively. For γ = -τ the derivatives are approximated with second-order accuracy, whereas they are approximated with third-order and second-order accuracy for γ = -τ + 1 3 .

advantage of the von Neumann analysis to infer the optimal γ. To allow analytical developments, we limit ourselves to the case Ma = 0. This is not restrictive since the dissipation behaviour only weakly depends on the Mach number provided that this latter remains small; the optimal γ depends essentially on τ . The eigenvalue problem Eq. (57) can then be solved analytically by means of a computer algebra system. The shear mode is of particular interest in our case, and is related to the root of the equation to yield γ( kx , τ ) = 4A( kx ) (τ -1) exp(

2 4 cos 3 ( kx ) -3 cos( kx ) + 2 γ τ λ 3 + 4τ cos 2 ( kx ) -12γ cos 2 ( kx ) -2 cos 2 ( kx ) -6γ + 2τ -1 τ λ 2
k2 x τ 3 ) + B( kx ) (-2τ + 1) exp( 2 k2 x τ 3 ) -6τ + 9 6A( kx ) exp( k2 x τ 3 ) -6B( kx ) exp( 2 k2 x τ 3 ) + 2C( kx ) exp( k2 x τ ) -6 , (63)
with A( kx ) = cos( kx ) + 2, B( kx ) = 2 cos 2 ( kx ) + 1 and C( kx ) = 4 cos 3 ( kx ) -3 cos( kx ) + 2. The function γ( kx , τ ) may be interpreted as the optimal parameter γ associated with each wavenumber kx of the planewave decomposition of the solution, at a given τ . It a priori ensures that the shear mode of the rfd scheme matches exactly the shear mode of the (linearized) Navier-Stokes equations at all wavenumbers. The colormap of γ( kx , τ ) and several isocontours are displayed in Fig. 11a. In practice, only an effective parameter γ(τ ) is selected in a simulation. Therefore, the idea would be to identify a colour covering the widest range of wavenumbers at fixed τ . We see in Fig. 11a that the light-blue colour covers the range of wavenumbers 0 ≤ kx ≤ π/4 for τ ≤ 0.05. Nonetheless, a steep variation of γ with kx is also observed at larger wavenumbers in the same range of relaxation coefficients, which indicates that approximating γ( kx , τ ) by a constant effective value can not be relevant at all wavenumbers. As a rule of thumb, one may claim from Fig. 11a that the scheme can only be optimized to resolve wavenumbers in the range 0 ≤ kx ≤ π/4 (at least eight points per wavelength) for τ ≤ 0.05. The situation improves for larger relaxation coefficients, in particular when τ 1/2 for which γ = 0 appears to be optimal at all wavenumbers. In that particular case, the rfd scheme consistently reduces to the lbm scheme and simply writes

f neq[0] α = - 1 2 f eq[0] α -f eq[-1] α (64) with ρ(x, t) = α f eq[-1] α and (ρu)(x, t) = α f eq[-1
] α c α . This coincidence is not surprising since the contribution of the non-equilibrium component of the distributions in the lbm scheme cancels out for the specific value τ = 1/2 (τ g = 1). However, such favourable regime is only attainable for flows at very low Reynolds numbers or by greatly increasing the resolution.

As expected, the dissipation rate of the shear mode displays the correct behaviour in Fig. 11b. Nevertheless, we also observe that an acoustic mode triggers an instability at kx ≥ π/3 thus showing that optimizing the dissipation at all wavenumbers while preserving the stability turns out to be unrealistic. As already mentioned, a more accessible goal seems to find an effective γ to optimize the dissipation at low wavenumbers only. In this regard, we can mention that

lim kx→0 γ( kx , τ ) = 1 3 τ 2 - 1 2 τ + 1 6 (65)
provides a lower bound in the limit of vanishing kx = k x δ x . In general, the value of the optimal γ should depend on the spectral decomposition of the simulated flow. If k * represents the characteristic energy-carrying wavenumber of the solution, γ( k * , τ ) expected to be relevant for the optimal γ provided that 0 ≤ k * ≤ π/4. The values obtained by the trial-and-error method are compared with the present analytical modeling in Fig. 12. In the case of the Taylor-Green vortex simulation, one can consider that the kinetic energy is concentrated at wavenumber k * ≈ 2π/L from the analytical form of the initial condition. The resolutions L = 32, 16 and 8 thus give respectively k * ≈ π/16, π/8 and π/4 (in lattice units). We observe in Fig. 12 that the experimental values obtained by trial-and-errors method are in excellent agreement with the theoretical prediction γ( k * , τ ) by varying τ for each level of resolution. Nonetheless, we observe an increasing discrepancy with τ that we attribute to errors in Mach number; the prediction is made at Ma = 0. The optimal γ decreases slightly with the velocity. For the double-shear-layer configuration, the solution involves fluid structures at different scales (see Fig. 3) which implies a much broader spectral decomposition. Different resolutions are considered with L = 128, 256 and 512. The Mach and Reynolds numbers are kept constant so that the viscosity (or τ ) increases proportionally to L. We observe that the optimal γ obtained by the trial-and-error method lies in the interval 0.18 -0.20 corresponding to a spectral decomposition typically in the range of wavenumbers π/8 -π/4. This is consistent with the number of lattice nodes used to resolve the fluid structures. The expected decrease of the optimal γ with the resolution L is also captured. However, we are not able to predict finely the optimal γ but only to provide an order of magnitude related to the spectral domain (in wavenumbers) of the solution. The lines correspond to the theoretical (von Neumann analysis) optimal γ( kx, τ ) as a function of τ for various wavenumbers kx at Ma = 0. The symbols correspond to the optimal γ obtained by the trial-and-error method (with an increment of 0.001) for the Taylor-Green vortex (TGV) and double-shear-layer (DSL) simulations.

Conclusion and perspectives

This work introduces on a rigorous basis an alternative equation to the discrete-velocity Boltzmann equation to address fluid dynamics at a kinetic level, while conforming to the Navier-Stokes dynamics at the macroscopic level (in the low-Mach-number limit). This kinetic equation expresses the full distribution functions, i.e. including the equilibrium and non-equilibrium components, in terms of the flow variables (ρ, u) only. However, this comes at the cost of considering a second-order differential equation instead of the first-order discrete-velocity Boltzmann equation:

f α = f eq α (ρ, u) -τ ∂ ∂t + c α • ∇ f eq α (ρ, u) + τ 2 ∂ ∂t + c α • ∇ 2 f eq α (ρ, u).
From a numerical viewpoint, a first practical application that emerges from this modeling is the possibility of designing kinetic schemes that depend only on ρ and u, and which, hopefully, would inherit the advantageous spectral properties of the original LB method. The discretization in space and time opens a Pandora's box of possible strategies. However, since the differential equation involves particle derivatives, it was naturally chosen to integrate it along particle trajectories and therefore resort to an underlying uniform lattice structure. In addition, backward finite differences were considered to ensure an explicit character, allowing the use of solvability conditions to obtain a scheme that depends only on flow variables. A parametric family of mass-conserving finite-difference schemes was derived. Numerical simulations in two dimensions and a von-Neumann analysis indicate that the developed schemes are consistent and stable but suffer from a level of numerical dissipation which is higher than the standard Lattice Boltzmann scheme. The von-Neumann analysis shows that this deficiency is essentially linked to the approximation by backward finite differences. An optimization of the stencil is possible to approximate the physical dissipation of the Navier-Stokes equations in the range 0 < kx < π/4, i.e. with a resolution greater or equal to eight points per wavelength. However, that requires prior information about the spectral decomposition of the solution itself, which is not entirely satisfactory from a numerical point of view. The use of central differences might certainly improves the situation, however in that case, we would lose the possibility of using the solvability conditions to evolve the solution in terms of the flow variables only. The tracking of the entire distributions (including equilibrium and non-equilibrium components) would therefore be a necessary compromise.

Finally, while admitting that backward finite differences are too dissipative to propagate accurately the solution over long time intervals, the developed finite-difference scheme may nevertheless be useful to locally reconstruct some missing distributions due to a defect in the lattice connectivity, e.g. at the transition between two domains with different resolutions or near a solid boundary. It may also be of interest in hybrid simulations to match Navier-Stokes and Lattice Boltzmann schemes, by providing a map between the macroscopic variables (ρ, u) and the f α 's. This point is currently under investigation.

By gathering the equations at O 1 and O 2 and multiplying each order by and 2 respectively, one obtains

∂ t ρ + ∂ i (ρu i ) = 0 (A.8) ∂ t (ρu i ) + ∂ j (ρu i u j + pδ ij ) + ∂ j ( Π (1)
ij ) = 0 (A.9)

where Π

(1) ij = -τ ρc 2 s (∂ j u i + ∂ i u j ) + τ ∂ k (ρu i u j u k ) (A.10) arises from Eqs. (A.3), (A.4) and (A.5) with the equilibrium function f eq α given by Eq. ( 5). The first term in Eq. (A.10) suitably represents a Newtonian viscous stress with dynamic viscosity τ ρc 2 s . However, the second term indicates that the Navier-Stokes equations are not exactly recovered. This discrepancy stems from a lack of isotropy of the discrete set of microscopic velocities, here associated with the D2Q9 lattice. This lattice is not capable to represent the third-order moment Π eq ijk correctly [START_REF] Dellar | Lattice boltzmann algorithms without cubic defects in galilean invariance on standard lattices[END_REF]. Developing the equilibrium functions to an higher order in Hermite expansion does not solve the problem; only off-diagonal terms vanish. The only remedy would be to resort to a lattice with a much higher level of connectivity, e.g. the D2Q21 lattice with 21 velocities. In practice, the error term is assumed as negligible (compared to the viscous stress) provided that the Mach number remains sufficiently small.

Figure 1 :

 1 Figure 1: Set of velocities of the D2Q9 lattice.

  and evaluate f * α for all directions α. The updated values of ρ and ρu are eventually computed as the statistical moments of f * α . A schematic view of the algorithm is given in Fig.2.

Figure 2 :

 2 Figure 2: Schematic view of the algorithm.

Figure 3 :

 3 Figure 3: Double shear-layer simulations. Normalized vorticity at convective time t * ≡ t U 0 /L = 1 for resolutions L = 128, L = 256 and L = 512. The lbm simulation is unstable at L = 128.

  u

Figure 4 :

 4 Figure 4: Double shear-layer simulations. Normalized averaged kinetic energy according to the convective time t * for the intermediate resolution L = 256. The reference simulation refers to the lbm simulation at high resolution L = 512.

Figure 5 :

 5 Figure 5: Taylor-Green vortex flow simulations. L 2 -norm (relative) error at time t = 2Tν for the different resolutions L = 16, L = 32, L = 64 and L = 128.

Figure 6 :

 6 Figure6: Taylor-Green vortex flow simulations. L 2 -norm (relative) error at time t = 2Tν for various relaxation times, τ = 3ν in lattice units. The optimal γ in the rfd algorithm is γ = 0.17 for τ = 0.0003 and τ = 0.003, γ = 0.15 for τ = 0.03 and γ = -0.07 for τ = 0.3.

Figure 7 :

 7 Figure 7: Taylor-Green vortex flow simulations. L 2 -norm (relative) error at time t = 2Tν . The optimal γ in the rfd simulation is γ = 0.15 for U 0 = 3.125 • 10 -4 , 3.125 • 10 -3 , 3.125 • 10 -2 and γ = 0.14 for U 0 = 0.3125.

Figure 8 :

 8 Figure 8: Taylor-Green vortex flow simulations. L 2 -norm (relative) error according to the normalized time t * ≡ t/Tν at Re = 100, L = 32 and ν = 0.01.

Figure 9 :

 9 Figure9: Von Neumann analysis: Dispersion and dissipation rates given by the real and imaginary parts of the pulsation of plane-wave perturbations around a uniform and stationary reference state at Ma = 0.1 and τ = 3 • 10 -3 . Only the dependence along the x-axis is here considered. For the analytical Navier-Stokes solution, the different waves correspond to two acoustic modes, ω ac± propagating at speed U ± cs and a shear mode, ωs,propagating at speed U ; the acoustic and shear modes are dissipated at rate -νk 2 . In the rfd, the value γ = 0 corresponds to an approximation of the first and second derivatives with second-order and first-order accuracy respectively. For γ = -τ the derivatives are approximated with second-order accuracy, whereas they are approximated with third-order and second-order accuracy for γ = -τ + 1 3 .

+ 2 Figure 10 :

 210 Figure 10: Von Neumann analysis: Dispersion and dissipation rates given by the real and imaginary parts of the pulsation of plane-wave perturbations around a uniform and stationary reference state at Ma = 0.1 and τ = 3 • 10 -1 .

  (a) Colormap of γ( kx, τ ) given by Eq. (63). Some isocontours are represented. Von Neumann analysis at Ma = 0, τ = 3 • 10 -3 .

Figure 11 :

 11 Figure 11: Optimal γ parameter at various wavenumbers and relaxation times, and its influence on the spectral properties of the scheme.

Figure 12 :

 12 Figure 12: The lines correspond to the theoretical (von Neumann analysis) optimal γ( kx, τ ) as a function of τ for various wavenumbers kx at Ma = 0. The symbols correspond to the optimal γ obtained by the trial-and-error method (with an increment of 0.001) for the Taylor-Green vortex (TGV) and double-shear-layer (DSL) simulations.

Table 1

 1

	gathers the parameters

Table 1 :

 1 Independent parameters of the performed simulations. The characteristic velocity U 0 is given by U 0 = Re ν/L and the Mach number is Ma = √ 3 U 0 .

Lattice units are obtained by normalizing space and time by the lattice spacing and the time step of the algorithm, respectively.

Let us mention that m =

and m =

should be included in the sum to go beyond Navier-Stokes equations and recover the Burnett and super-Burnett equations, respectively.

The BGK collision operator conserves microscopic mass and momentum.

Funding

The authors acknowledge the funding of the french DGAC under the project OMEGA3 (DGAC/DTA/SDC n o 2018-16). This work part of a scientific collaboration including CS-group, Renault, Airbus, Ecole Centrale de Lyon, CNRS and Aix-Marseille University.

Appendix A. Navier-Stokes equations from Chapman-Enskog expansion

At O( 0 ), f

Taking the zeroth, first, and second moments of Eq. ( 10) yields

t Π eq ij + ∂

(1)

with Π eq ij = α c α,i c α,j f eq α = ρu i u j + ρc s 2 δ ij , Π eq ijk = α c α,i c α,j c α,k f eq α and Π

(1)

α . The first two equations correspond to the Euler equations with the pressure being defined by p = ρc s 2 . Similarly, computing the zeroth, first, and second moments of Eq. [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the navier-stokes equation[END_REF] gives ∂