This article presents a new numerical scheme designed to solve for any scalar equation coupled with a Lattice-Boltzmann solver (in so-called hybrid methods). Its most direct application is to solve an energy equation, in parallel with a Lattice-Boltzmann solver dealing with mass and momentum conservation.

The numerical scheme is specifically designed to compute the energy flux consistently with the mass and momentum flux (as is done, for instance, using Riemann solvers). This scheme effectively lifts a major limitation of current compressible hybrid Lattice-Boltzmann, in which the energy conservation is tackled under nonconservative form, leading to discretization errors on jump conditions across shocks.

I. INTRODUCTION

Lattice Boltzmann Methods (LBM) has received a fastly growing interest in the past two decades. They are now recognized as a powerful tool for Computational Fluid Dynamics, and as a potential candidate for a breakthrough toward revolutionary CFD [START_REF] Stolnick | Cfd vision 2030 study: A path to revolutionary computational aerosciences[END_REF] . Most published works dealing with LBM are related to low Mach number continuum flows, covering a broad range of flow physics, among which multiphase flows, thermal flows, combustion and non-Newtonian flows (see [START_REF] Guo | The Lattice Boltzmann Method and its applications in engineering[END_REF][START_REF] Krüger | The Lattice Boltzmann Method. Principles and Practice[END_REF] for examples). Therefore, designing an efficient Lattice Boltzmann Method well suited for viscous high speed subsonic and supersonic flows remains a challenging issue, since the Asymptotic Preservation (AP) as defined by Filbet and Jin 4,[START_REF] Jin | Asymptotic preserving (ap) schemes for multiscale kinetic and hyperbolic equations[END_REF] for this case, i.e. the recovery of the compressible Navier-Stokes equations is tricky (see, e.g. [START_REF] Bennoune | Uniformly stable numerical schemes for the boltzmann equation preserving the compressible navier-stokes asymptotics[END_REF][START_REF] Chen | A comparative study of an asymptotic preserving scheme and unified gas-kinetic scheme in continuum flow limit[END_REF][START_REF] Kolobov | Unified solver for rarefied and continuum flows with adaptive mesh and algorithm refinement[END_REF][START_REF] Yang | Development of discrete gas kinetic scheme for simulation of 3d viscous incompressible and compressible flows[END_REF] for examples dealing with AP kinetic schemes). It is reminded here that a Lattice Boltzmann Method has the AP property, if, at fixed mesh size and time step, a consistent discretization of the Navier-Stokes equations is recovered in the limit of vanishing Knudsen number.

Since LBM can be derived as a discretization of the Boltzmann equation [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the navier-stokes equation[END_REF] , it can inherit from the later the capability to handle compressible flow physics. Here, the problem lies in the efficiency of LBM rather than in an intrinsic inability to capture compressible flows (assuming that the conservativity of the method has been enforced), since to recover accurately compressible flow physics a direct and "rough" discretization of the Boltzmann equation in velocity, space and time is known to lead to the definition of numerical methods with non-standard lattices that may involve a large number of discrete velocities [START_REF] Krüger | The Lattice Boltzmann Method. Principles and Practice[END_REF][START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the navier-stokes equation[END_REF] .

A common way to decrease the computational cost and the memory storage is to use a segregated approach, in which the energy equation is treated separately from the mass and momentum conservation equations. Doing so, it is expected to that it is possible to drastically reduce the number of unknowns per cell (i.e. the number of discrete velocities) along with a decrease of the algorithmic complexity of the equation associated to each unknown (i.e. the computational complexity of the collision model). On the opposite, in monolithic approaches all equations are treated simultaneously in the same way, e.g. [START_REF] Watari | Finite difference lattice boltzmann method with arbitrary specific heat ratio applicable to supersonic flow simulations[END_REF][START_REF] Zhang | A lattice boltzmann model for the compressible euler equations with second-order acuracy[END_REF][START_REF] Gan | Flux limiter lattice boltzmann scheme approach to compressible flows with flexible specific-heat ratio and prandtl number[END_REF][START_REF] Kataoka | New lattice boltzmann model for the compressible navierstokes equations[END_REF][START_REF] Gan | Discrete boltzmann trans-scale modeling of high-speed compressible flows[END_REF][START_REF] Atif | Higher-order lattice boltzmann model for thermodynamics[END_REF][START_REF] Fares | Validation of a lattice-boltzmann approach for transonic and supersonic flow simulations[END_REF] . Segregated Lattice Boltzmann methods can be grouped into two main classes. The first one corresponds to Double Distribution Function (DDF) approach in which a second set of distribution functions is defined to solve the energy equation via a kinetic advection-relaxation approach, e.g. [START_REF] Feng | A compressible lattice boltzmann finite volume model for high subsonic and transonic flows on regular lattices[END_REF][START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF][START_REF] Frapolli | Entropic lattice boltzmann model for compressible flows[END_REF][START_REF] Saadat | Lattice boltzman model for compressible flows on standard latices: variable prandtl number and adiabatic exponent[END_REF] . The second one is related to hybrid methods, in which the macroscopic energy equation is solved using a classical Finite Volume/Difference/Element method, e.g. [START_REF] Feng | Regularized thermal lattice boltzmann method for natural convection with large temperature differences[END_REF][START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Feng | Hybrid recursive regularized lattice boltzmann simulation of humid air with application to meteorological flows[END_REF] . Both approaches have shown promising results in the recent past, and are being actively developed by various research groups.

Segregated methods, as all methods based on a discretization of the Boltzmann equation, should be Asymptotic Preserving (AP) methods and therefore should be conservative. The conservativity issue has been addressed by several authors for methods originating in Boltzmann kinetic equation (e.g., see [START_REF] Aristov | The conservative splitting method for solving boltzmann's equation[END_REF][START_REF] Titarev | Conservative numerical methods for model kinetic equations[END_REF][START_REF] Mieussens | Discrete-velocity models and numerical schemes for the boltzmann-bgk equation in plane and axisymmetric geometries[END_REF][START_REF] Pieraccini | Implicit-explicit schmes for bgk kinetic equations[END_REF] ), and it appears that several features of usual LBM can lead to the loss of conservativity:

• Since LB methods are discrete ordinate/velocity methods with fixed discrete velocities, their collision models must be adapted in order to get an exact/conservative evaluation of the moments of interest and the preservation of the associated collision invariants, while still recovering the targeted macroscopic equations. Adaptation of the collision kernels can be done in several ways, among which supplementing a truncated polynomial expansion of the Maxwellian equilibrium function found in gas kinetic theory by forcing terms to balance truncation errors (e.g. [START_REF] Feng | A compressible lattice boltzmann finite volume model for high subsonic and transonic flows on regular lattices[END_REF][START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF][START_REF] Saadat | Lattice boltzman model for compressible flows on standard latices: variable prandtl number and adiabatic exponent[END_REF][START_REF] Feng | Regularized thermal lattice boltzmann method for natural convection with large temperature differences[END_REF][START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF] ), or building a synthetic collision model that explicitly enforces some constrains (usually solving a nonlinear problem at every grid point and time step), e.g. [START_REF] Watari | Finite difference lattice boltzmann method with arbitrary specific heat ratio applicable to supersonic flow simulations[END_REF][START_REF] Zhang | A lattice boltzmann model for the compressible euler equations with second-order acuracy[END_REF][START_REF] Kataoka | New lattice boltzmann model for the compressible navierstokes equations[END_REF][START_REF] Gan | Discrete boltzmann trans-scale modeling of high-speed compressible flows[END_REF][START_REF] Atif | Higher-order lattice boltzmann model for thermodynamics[END_REF][START_REF] Frapolli | Entropic lattice boltzmann model for compressible flows[END_REF][START_REF] Latt | Efficient supersonic flows through higherorder guided equilibrium with lattice boltzmann[END_REF][START_REF] Qiu | A double-distribution-function lattice boltzmann model for high-speed compressible viscous flows[END_REF] . It is worth noting that in the former case, many authors refer to conservativity error as Galilean invariance error or symmetry-breaking errors, since the error scales as a power of the velocity.

• Time integration, including the splitting between the streaming and the collision steps, must be performed in a carefully chosen way, e.g. [START_REF] Aristov | The conservative splitting method for solving boltzmann's equation[END_REF] .

To the knowledge of the authors, a complete theoretical analysis of the AP and conservation properties of segregated LBM for compressible viscous flows is still missing. Segregated hybrid LBM for thermal compressible flows renders the problem theoretically more difficult, since it mixes a mesoscopic approach for mass and momentum conservation with a macroscopic approach for energy. Numerical experiments have shown that it is a powerful approach, but recent theoretical analyses and numerical experiments [START_REF] Guo | Evaluation on formulation of energy conservation equation in hybrid compressible lattice boltzmann method[END_REF] have shown that some conservativity issues may arise that can lead to limitations in terms of numerical stability and convergence toward the correct shock jump relation for strong shock waves. More pre-cisely, a lack of consistency in the mass flux ρu evaluation between the mass and momentum equation on one hand, and the energy equation on the other hand has been identified. This problem can be partially cured by using a non-conservative form of the energy equation and choosing a thermodynamic variable whose evolution equation does not explicitly exhibit a pressure work term p∇•u, e.g. the entropy as in [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF] . While increasing the stability of the method, this choice doesn't lead to a fully conservative method, and Rankine-Hugoniot shock jump relations are not exactly recovered for strong shock waves. This can be understood considering the advection of a conserved scalar quantity φ in parallel of the mass conservation equation:

∂ρ ∂t I + (∇ • (ρu)) I = 0 ( 1
)
∂ρφ ∂t II + (∇ • (ρuφ)) II = 0 (2)
where subscripts I and II refer to LBM and macroscopic discretization, respectively. The last conservative equation can be rewritten as a combination of a non-conservative equation and the mass conservation equation

ρ ∂φ ∂t + φ ∂ρ ∂t II + (φ∇ • (ρu) + ρu • ∇φ) II = 0 (3)
from which one obtains

ρ ∂φ ∂t + u • ∇φ II = -φ ∂ρ ∂t + ∇ • (ρu) II (4)
showing that the scalar is conserved and the conservative and non-conservative forms of the scalar equation are equivalent if the right hand side of Eq. ( 4) is null. Since the mass flux is estimated via Eq. ( 1), one can see that an error proportional to φ

[(∇ • (ρu)) I -(∇ • (ρu)) II ]
appears in hybrid methods based on a conservative form of the scalar equation.

The present paper aims at designing a fully conservative hybrid LBM approaches for thermal compressible flows that will satisfy the AP condition for the continuum compressible regime, including solutions with shock waves. To this end, a new discretization of the macroscopic thermodynamic variable equation is introduced, that will allow for the use of the conservative form of the total energy, resulting in a dramatic improvement of both the conservativity and stability of the method. It relies on the idea that the macroscopic and the kinetic equation must be coupled in a deeper way than in previous hybrid LBM approaches, in order to alleviate the inconsistency in the mass flux evaluation between the two parts of the global method. To this end, it is proposed to evaluate the fluxes in the macroscopic equation for the thermodynamic variable directly from the distribution functions of the kinetic part.

The idea of using a kinetic approach to evaluate numerical fluxes in a macroscopic evolution equation is not new. As a matter of fact, it is the cornerstone of Gas Kinetic Schemes (GKS) proposed during the 1970's for Euler and Navier-Stokes equations (see [START_REF] Sanders | The possible relation of the 3-kiloparsec arm to explosions in the galactic nucleus[END_REF][START_REF] Pullin | Direct simulation methods for compressible inviscid ideal-gas flow[END_REF][START_REF] Benoit Perthame | An Introduction to Kinetic Schemes for Gas Dynamics[END_REF][START_REF] Xu | Numerical navier-stokes solutions from gas kinetic theory[END_REF][START_REF] Xu | Gas-kinetic relaxation (bgk-type) schemes for the compressible euler equations[END_REF][START_REF] Kun | Direct modelling for Computational Fluid Dynamics[END_REF] ) and their sequels, among which UGKS [START_REF] Chen | A comparative study of an asymptotic preserving scheme and unified gas-kinetic scheme in continuum flow limit[END_REF][START_REF] Xu | A unified gas-kinetic scheme for continuum and rarefied flows[END_REF][START_REF] Xu | An improved unified gas-kinetic scheme and the study of shock structures[END_REF][START_REF] Chen | A unified gas-kinetic scheme with moving mesh and velocity space adaptation[END_REF] , DUGKS [START_REF] Guo | Discrete unified gas kinetic scheme for all knudsen flows: Low-speed isothermal case[END_REF][START_REF] Guo | Discrete unified gas kinetic scheme for all knudsen flows. ii. thermal compressible case[END_REF][START_REF] Zhu | Discrete unified gas kinetic scheme on unstructured meshes[END_REF] , DGKS [START_REF] Yang | Development of discrete gas kinetic scheme for simulation of 3d viscous incompressible and compressible flows[END_REF][START_REF] Zhang | Double distribution function-based discrete gas kinetic scheme for viscous incompressible and compressible flows[END_REF] , FDGKS [START_REF] Xuan | An efficient high-order finite difference gas-kinetic scheme for the euler and navier-stokes equations[END_REF] , LBMFS [START_REF] Yang | A moment conservation-based non-free parameter compressible lattice boltzmann model and its application for flux evaluation at cell interface[END_REF][START_REF] Wang | Development of lbgk and incompressible lbgk-based lattice boltzmann flux solvers for simulation of incompressible flows[END_REF][START_REF] Wang | From lattice boltzmann method to lattice boltzmann flux solver[END_REF] and SLBM [START_REF] Chen | On improvements of simplified and highly stable lattice boltzmann method: Formulations, boundary treatment, and stability analysis[END_REF][START_REF] Chen | Three-dimensional simplified and unconditionally stable lattice boltzman method for incompressible isothermal and thermal flows[END_REF][START_REF] Chen | A simplified thermal lattice boltzman method without evolution of distribution functions[END_REF][START_REF] Chen | A simplified lattice boltzman method without evolution of distribution function[END_REF] . All these approaches are monolithic ones, and therefore they cannot be used straightforwardly to improve segregated hybrid LBM.

In a separate way, some authors have proposed so-called Lattice Boltzmann schemes for scalar transport [START_REF] Onishi | A lattice boltzmann model for polymeric liquids[END_REF][START_REF] Osmanlic | Lattice boltzmann method for oldroyd-b fluids[END_REF][START_REF] Ve Küng | Comparison of passive scalar transport models coupled with the lattice boltzmann method[END_REF][START_REF] Krivovichev | Analysis of the parametric models of passive scalar transport used in the lattice boltzmann method[END_REF][57] , with application to the transport of stresses in non-Newtonian fluid models and passive scalar transport. In these schemes, the advection term for the scalar quantity is directly computed by combining the distribution functions and the scalar quantity, rather than computing the macroscopic velocity and then applying a usual Finite Difference/Finite volume Scheme to solve the macroscopic equation. An important point is that these schemes have all been derived to recover the non-conservative form of the scalar equation, assuming flow incompressibility.

In order to improve conservativity and recover AP property for the continuum compressible flow regime, it is proposed in the present paper to improve the previously proposed hybrid LBM methods by extending the so-called Lattice Boltzmann scalar schemes to energy equation. A main novelty of the present method is that mass and momentum equation are solved in a classical LBM way, while a macroscopic evolution equation is still used for energy, the fluxes of the later being evaluated from distribution functions. To the knowledge of the authors, such a method has not been proposed yet.

This paper is organized as follows: the scalar transport schemes with LBM populations will be presented in Section II. Conservation properties and stability are assessed using the classical Taylor-Green Vortex (TGV) test case in Section IV. Compressible test cases with shocks are then addressed in Section V to demonstrate the performance of the schemes when discontinuities are encountered. Conclusions are drawn in Section VI.

II. SCALAR TRANSPORT USING LBM DISTRIBUTION FUNCTIONS A. Scalar flux schemes using collision populations

The present idea of using LBM populations for scalar transport comes from the Taylor analysis of LBM schemes. The streaming step of almost any LBM scheme can be expressed as

f i (x, t + ∆t) = f col i (x -c i ∆t, t) (5)
where f i is the population and f col i the collide population to be streamed [START_REF] Krüger | The Lattice Boltzmann Method. Principles and Practice[END_REF] , and c i are the discrete velocities of the considered DnQm lattice. Their detailed expressions depend on the LBM scheme, but Eq. ( 5) always gives the transport equation for the moments of

f i Π f α 1 •••αn ≡ i c iα 1 • • • c iα n f i (6)
through Taylor expansion, where the vector c i denotes the lattice velocity of direction i.

Multiplying Eq. ( 5) by the solution-independent term c iα 1 c iα 2 • • • c iα n and then writting the Taylor series expansion of its R.H.S around x, one obtains

Π f α 1 •••αn (x, t + ∆t) = i c iα 1 • • • c iα n f col i (x -c i ∆t, t) =Π f col α 1 •••αn -∆t ∂ ∂x α n+1 Π f col α 1 •••α n+1 + ∆t 2 2 ∂ 2 ∂x α n+1 ∂x α n+2 Π f col α 1 •••α n+2 + O(∆t 3 ) , (7)
where the variables without any space or time indicators denote the corresponding values at (x, t) for the sake of simplicity. This equation indicates a numerical scheme for the divergence of moments

∇ L • (Π f col α 1 •••α n+1 ) ≡ 1 ∆t Π f col α 1 •••αn - i c iα 1 • • • c iαn f col- i = ∂ ∂x α n+1 Π f col α 1 •••α n+1 - ∆t 2 ∂ 2 ∂x α n+1 ∂x α n+2 Π f col α 1 •••α n+2 + O(∆t 2 ) , (8)
where the upper script -denote the variables at (xc i ∆t, t). For instance, the spatial derivative of the first order moments of the collide population, usually the momentum Π α ≡ i c iα f col i = ρu α , can be evaluated as

∇ L • (ρu) ≡ 1 ∆t i f col i -f col i -= ∂ ∂x α ρu α - ∆t 2 ∂ 2 ∂x α ∂x β Π f col αβ + O(∆t 2 ), (9)
seemingly a first-order approximation. Yet, combining this equation with Eq. ( 7) for next order moment Π αβ yields

Π f α 1 •••αn (t + ∆t, x) =Π f col α 1 •••αn - ∆t 2 ∂ ∂x α n+1 Π f col α 1 •••α n+1 + Π f α 1 •••α n+1 (x, t + ∆t) + O(∆t 3 ) , (10)
showing that the numerical operator ∇ L is second-order accurate for the LBM mass flux evaluation, without any artificial viscous terms

∇ L • (ρu) = 1 2 ∂ ∂x α ρu α + ∂ ∂x α ρu α (t + ∆t, x) + O(∆t 2 ) . (11)
It is worth noting that the above equation just serves as a demonstration of the order of accuracy in the evaluation of the mass flux in LBM. In practice, values at the next time step t + ∆t are never used, i.e. the LBMs considered in the current study are always explicit.

More illustrative deductions of the above equations can be found in appendix A.

Let us now consider the multiplication of a scalar field φ and the collide population f col , denoted by

f φ i ≡ φf col i , (12)
whose moments are

Π f φ α 1 •••αn ≡ i c iα 1 • • • c iα n φf col i = i φΠ f col α 1 •••αn (13)
Inserting the first moments of f φ i into Eq. ( 8), one obtains the following numerical scheme for the scalar flux

∇ U • (ρuφ) ≡ 1 ∆t i f φ i -f φ i - = ∂ ∂x α ρu α φ - ∆t 2 ∂ 2 ∂x α ∂x β Π f φ αβ + O(∆t 2 ) (14)
Unfortunately, this scheme is only first order accurate with a significant artificial viscous term proportional to the second moments of f φ . This dissipation term can not be canceled by time marching like the mass flux. Indeed, unlike the mass flux, the scalar flux ρu α φ is not transported by the population f φ anymore. So that Eq. ( 10) can not be applied to cancel the artificial viscous terms, as discussed in Appendix A.

However, Some modifications can be performed to make the flux more "central" in order to reduce the numerical dissipation. Instead of using fully upwind reconstruction of φ, one can use a centered reconstruction the scalar, leading to the following numerical divergence operator

∇ C • (ρuφ) ≡ 1 ∆t i f col i φ + + φ 2 -f col i -φ + φ - 2 ( 15
)
where the upper script + denotes the variables at (x + c i ∆t, t). The Taylor expansion of this equation given in Appendix B shows that this centered scheme does not bring any artificial viscous terms with respect to the scalar field.

It is important noting that the terms "upwind" or "centered" in this context are defined with respect to the lattice velocities. In practice, both ∇ U and ∇ C operators use the values of the scalar and collision populations from all neighbor nodes in the lattice, which is quite different from classical upwind schemes.

It is worth mentioning that in both numerical divergence operators ∇ U and ∇ C , the mass flux part in the scalar flux is numerically identical to the mass flux resolved by the LBM. These numerical schemes are computationally affordable under the hybrid LBM-FD framework, because the collision population is already available. Indeed, these scalar flux operators are just special forms of finite difference method (FDM) except that, the coefficients are now accounting for the mass transport. In other words, the operators proposed here are just special forms of FD spatial derivation schemes which intend to inherit some numerical properties of double distribution approaches, with a reduced computational cost.

B. Bridging with other LBM scalar transport schemes

The idea of using LBM populations to compute the scalar transport can be found in several papers, see [START_REF] Onishi | A lattice boltzmann model for polymeric liquids[END_REF][START_REF] Osmanlic | Lattice boltzmann method for oldroyd-b fluids[END_REF][START_REF] Ve Küng | Comparison of passive scalar transport models coupled with the lattice boltzmann method[END_REF] . Generally speaking, these works focus on the transport of a passive scalar in non-conservative form, by evaluating the u • ∇φ term using the LBM populations from neighboring cells. For instance, Onishi et al. [START_REF] Onishi | A lattice boltzmann model for polymeric liquids[END_REF] consider that the mass transport from cell x to x + c i ∆t by streaming is equal to f col i (x); and that the one from cell x + c i ∆t to x is given by f col ī (x + c i ∆t), where the index ī corresponds to the lattice whose direction is opposite to lattice i. Thus, the scalar transport is modeled as

(u • ∇φ) O ≡ 1 ∆t i φf col i ρ (x) - φf col ī ρ (x + c i ∆t) , (16)
According to [START_REF] Onishi | A lattice boltzmann model for polymeric liquids[END_REF][START_REF] Ve Küng | Comparison of passive scalar transport models coupled with the lattice boltzmann method[END_REF] , it is also possible to reduce the numerical dissipation by modifying the population used to transport the scalar :

(u • ∇φ) OC ≡ 1 ∆t i φ f col i ρ -ω i (x) -φ f col ī ρ -ωī (x + c i ∆t) , (17)
Other authors [START_REF] Osmanlic | Lattice boltzmann method for oldroyd-b fluids[END_REF] noticed that the net mass exchange through the interface x + 1 2 c i ∆t into cell x is equal to

∆m i ≡ f col ī (x + c i ∆t) -f col i (x) , (18)
so that the total scalar exchange can be expressed by only considering the "upstream" values, leading to

(u • ∇φ) OS ≡ i -∆m i ρ(x, t + ∆t)    φ(x, t) if ∆m i ≤ 0 φ(x + c i ∆t, t) otherwise. (19)
It is worth noting that these numerical schemes for non-conservative scalar flux cannot be applied directly in compressible flows. As a matter of fact, one can demonstrate that these schemes rely on the incompressibility assumption, i.e. ∇ • u ≈ 0, to work properly. In order to compare these schemes with the ∇ U and ∇ C operators in compressible configurations, their conservative extensions given in Appendix C are used as references in the present work.

C. Applications in the energy transport for HRR-p LBM

As mentioned in the introduction, the scalar transport methods designed in the current work aim at curing the incompatibility between mass flux evaluations in the energy and mass conservation equations in hybrid LBM-FD approaches. For instance, in HRR-p LBM proposed in [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] , the governing equations for the flow fields are

∂ρ ∂t + ∂ρu α ∂x α = 0 (20a) ∂ρu α ∂t + ∂(ρu α u β + pδ αβ + T αβ ) ∂x β = 0 (20b)
supplemented by an entropy equation in non-conservative form for the energy transport

∂s ∂t + u α ∂s ∂x α = 1 ρT T αβ ∂u α ∂x β - ∂q β ∂x β (21)
where the viscous tensor reads

T αβ = µ ∂u α ∂x β + ∂u β ∂x α -δ αβ 2 3 ∂u γ ∂x γ (22)
with µ as dynamic viscosity, and the heat flux

q β = -λ ∂T ∂x β ( 23
)
where λ is the heat conductivity. The entropy field is connected to the flow fields through the ideal gas equation of state (EOS)

p = ρrT , s = C v ln p ρ γ (24)
with C v being the constant volume specific heat capacity, r the specific gas constant and γ ≡ C p /C v the adiabatic exponent.

In the original HRR-p method [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] , the flow fields (ρ and ρu) are computed by lattice Boltzmann method, while the entropy equation is updated using an explicit Euler method for time-integration supplemented by MUSCL scheme for advection along with a second order centered scheme for viscous terms. In practice, the density field evolves according to

ρ(x, t + ∆t) = ρ - p c 2 s + i f col i - = ρ - i f col i -f col i - , (25)
putting in evidence an identical mass flux (

f col i -f col i -) compared to classical density
based LBMs. The momentum fields are updated as

[ρu α ](x, t + ∆t) = ρu α - i c iα f col i -f col i - . (26)
The collision populations f col i are evaluated through a hybrid recursive regularized (HRR) collision procedure, yielding

f col i = f eq i + (1 - ∆t τ )f neq i + ∆t 2 F E i ( 27
)
where the equilibrium part (f eq i ), the non-equilibrium part (f neq i ) and the forcing terms for viscous tensor (F E i ) are evaluated from their projections onto the D3Q19 rotational symmetry basis of Gauss-Hermite polynomials up to third order [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann method-based large eddy simulation[END_REF] ,

f eq i = ω i p c 2 s + H (1) iα c 2 s ρu α + H (2) iαβ 2c 4 s ρu α u β + H (3r) iγ 6c 6
s a (3r),eq γ (28a)

f neq i = ω i H (2) iαβ 2c 4 s a (2),neq αβ + H (3r) iγ 6c 6 s a (3r),neq γ (28b
)

F E i = ω i H (2) iαβ 2c 4 s a (2),E αβ (28c)
where ω i is the Gaussian weight associated to the discrete velocity c i , p the pressure, and

c s ≡ ∆x/( √ 3∆t
) the lattice speed of sound. Further details about this collision operation can be found in Appendix D. The role of forcing term F E i is to enforce conservativity (in the sense given in Refs. [START_REF] Aristov | The conservative splitting method for solving boltzmann's equation[END_REF][START_REF] Titarev | Conservative numerical methods for model kinetic equations[END_REF] ) by balancing errors coming from both truncation of the polynomial expansion of the polynomial expansion of the Maxwellian and the quadrature error associated to the lattice.

The objective of the current work is to couple the conservative form energy equation with the LBM via this hybrid approach. To this end, instead of using the non-conservative formulation of entropy equation (21), the conservative form of total energy or entropy equations are discretized with the proposed numerical schemes, i.e.

∂ρE ∂t + ∂ρHu α ∂x α = ∂T αβ u α ∂x β - ∂q β ∂x β (29a) ∂ρs ∂t + ∂ρsu α ∂x α = 1 T T αβ ∂u α ∂x β - ∂q β ∂x β (29b)
where H ≡ E + p/ρ = e + 1 2 u α u α + p/ρ is the total enthalpy, and e = C v T the inner energy for idea gas. It is worth noting that when the scalar transport schemes are applied to the total energy equation, the scalar transported should be total enthalpy (H) instead of total energy (E). In the current work, the viscous terms in the R.H.S of Eq. ( 29) are resolved by second order central FD schemes same as in our previous works [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF][START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF] .

Using conservative formulation in energy transport should lead to conservation of the energy scalar. This is numerically guaranteed at the discrete level by the proposed schemes , i.e. the sum of all the scalar fluxes over any periodic domain is zero:

x ∇ U • (ρuφ) ≡ x 1 ∆t i f col i φ -f col i -φ - = 0 . x ∇ C • (ρuφ) ≡ x 1 ∆t i f col i φ + + φ 2 -f col i -φ + φ - 2 = 1 2∆t i x f col i φ + -f col- i φ + i x f col i φ -f col- i φ - = 0 . ( 30
)
The major benefit of the proposed schemes is the compatibility between the mass flux evaluation in the scalar transport equation and the one resolved by LBM. This feature is highlighted in the centered scheme ∇ C , granting it a skew-symmetry property, since it also conserves the associated quadratic invariant ρφ 2 (see [START_REF] Coppola | Numerically stable formulations of convective terms for turbulent compressible flows[END_REF] for a discussion about quadratic invariant conservation by numerical schemes). This can be demonstrated by examining the semi-discrete form of the scalar equation:

x δ t (ρφ 2 ) ≡ x 2φδ t (ρφ) -φ 2 δ t ρ ≡ - 1 ∆t x 2φ∇ C • (ρuφ) -φ 2 i f col i -f col- i = - 1 ∆t x φ i f col i φ + -f col- i φ + f col i φ -f col- i φ --φ 2 i f col i -f col- i = - 1 ∆t x i φf col i φ + -φf col- i φ - = - 1 ∆t i x φ + (φf col i ) -φ(φ -f col- i ) = 0 (31)
Here the δ t denotes the semi-discretized time derivation operator which only consider the convective part of the scalar evolution [START_REF] Coppola | Numerically stable formulations of convective terms for turbulent compressible flows[END_REF] . It is clear that invariant conservation is due to the identical mass flux evaluations in both the scalar flux and those resolved by LBM in Eq. (25). To the best of the authors' knowledge, none of the scalar transport schemes in the literature has such feature.

III. GLOBAL ACCURACY OF THE HYBRID LATTICE BOLTZMANN METHOD

Although the analysis based on Taylor expansion displayed in Appendix A gives the order of accuracy of the proposed spatial operators, the global order of segregated solver should be studied using numerical test cases. To this end, an one-dimensional entropy spot advection is considered [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] . This test case is initialized with an uniform initial velocity field corresponding to u = 1 (approximately at Ma 0.2) on a periodic one-dimensional domain of

L = 1.
Other variables are initialized according to

ρ = ρ 0 1 + δ exp -(x -x 0 ) 2 R 2 (32)
T = T 0 1 -δ exp -(x -x 0 ) 2 R 2 (33)
In practice, the parameters are chosen as T 0 = 20, ρ 0 = 1, x 0 = 0.5, R = 0.05 and δ = 0.001. The time step is given such that the CFL number is around 0.5. The energy scalar transported in the current test is the total energy following Eqn. (29a), without any viscous term to mimic a Euler calculation. Flow fields after one flow-over time are compared to their reference solutions (which are their initial values for this specific test). The global order of the solver is studied by the relative error defined as

(φ) ≡ (φ -φ ref ) 2 dx φ 2 ref dx (34)
Relative errors for density, velocity and energy are plotted versus the number of points in the calculation domain as figure 1. It is quite clear that the global convergence rate of the hybrid LBM solver with the total energy equation discretized by the upwinded operator ∇ U achieves first order, and the global order of the center operator ∇ C simulation is around two. Note that the flat tails in Fig. 1(b) merely indicate that the velocity fields reach the range of the machine tolerance with the chosen grid.

IV. TEST CASE: COMPRESSIBLE TAYLOR-GREEN VORTEX

The inviscid Taylor-Green vortex (TGV) test-case is addressed to assess the conservation properties of the proposed numerical schemes. It is calculated in a 2π × 2π × 2π periodic domain with only 32 grid points in each direction to emphasize the stability feature of the schemes. The viscosity is set to 10 -25 kg/m/s to mimic an inviscid condition for LBM solvers [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF] . The initial fields are given by Eqn (35), corresponding to constant density initial conditions (CDIC) at given Ma numbers [START_REF] Shu | Numerical convergence study of nearly incompressible, inviscid taylor-green vortex flow[END_REF][START_REF] Peng | Effects of the mach number on the evolution of vortex-surface fields in compressible taylor-green flows[END_REF] .

u x = sin x cos y cos z (35a) u y = -cos x sin y cos z (35b) u z = 0 (35c) ρ = 1 (35d) p = 1 γMa 2 +
∂s ∂t + u α ∂s ∂x α = 0 (36d)
The second order accurate isotropic spatial derivation operator used in [START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF] or MUSCL schemes employed in [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] are applied for the energy scalar gradients in primitive form. The scalar flux in conservative form are discretized using either the ∇ C and ∇ U operators, or the low-dissipation scalar transport scheme proposed in [START_REF] Onishi | A lattice boltzmann model for polymeric liquids[END_REF] in its conservative form given by Eq. C3. The original scheme proposed in [START_REF] Onishi | A lattice boltzmann model for polymeric liquids[END_REF] given by Eq. ( C1) is numerically equivalent to our ∇ U operator, therefore it is not tested separately. Also, although the operator ∇ OS proposed in [START_REF] Osmanlic | Lattice boltzmann method for oldroyd-b fluids[END_REF] is tested in its extended conservative form given by Eq. (C4), the results are not reported in the current work because they are very close to those obtained by operator ∇ U . The time integration of the scalar equation is performed using a first-order explicit Euler scheme to achieve a two-way coupling between energy scalar solvers and the LBM flow solver [START_REF] Feng | A lattice-boltzmann model for low-mach reactive flows[END_REF] , with an acoustic CFL number close to 0.5.

The variables that are conserved numerically are listed in Table II. It is worth noting that the HRR-p LBM flow solver always conserves density and momentum, as proved in Appendix E. As demonstrated in Eq. ( 31), the case scons should conserve both the linear invariant (ρs) and the quadratic invariant ρs 2 . However, in the case Econs, ρE 2 is only globally conserved due to presence of the pressure work pu α in the flux other than only ρEu α , as pointed out in [START_REF] Coppola | Numerically stable formulations of convective terms for turbulent compressible flows[END_REF] . Two distinct Ma numbers, Ma=0.08 and Ma=0.8, are chosen to cover both nearly incompressible and compressible regimes.

A. Nearly incompressible case

The temporal evolution of the concerned quantities when Mach number is 0.08 are displayed in Fig. 2. The results are reported by the normalized displacement from their original

f ≡ f -f 0 |f 0 | ≡ f (t)dV -f (0)dV | f (0)dV| . ( 37
)
All numerical configurations succeed in conducting this test up to t = 10 except Econs-Onishi and s-cons-Onishi, where the low-dissipation operators ∇ OC proposed in [START_REF] Onishi | A lattice boltzmann model for polymeric liquids[END_REF] is employed. The simulations s-cons-Onishi becomes unstable almost immediately after the sinitial time and therefore the results are not plotted. The simulation E-cons-Onishi remains stable up to t 1. This is not surprising because it can be shown from a Taylor analysis of the operator ∇ OC that, its numerical dissipation is proportional to ρuu, which is too small to stabilize the simulation at this Mach number. It is seen from Fig. 2 that in this low Mach configuration, the density and momentum are conserved numerically in all numerical recipes at a tolerance of 10 -10 . The kinetic energy from different schemes are decaying at almost the same rate as shown in Fig. 2(c), which is reasonable for hybrid LBM solvers because at low Mach number, the variation of kinetic energy should be dominated by the dissipation properties of the LBM flow solver which is identical in all setups. The kinetic energy maintains approximately 70% of its original value at t = 7, which is even better than the performance of some high-order compressible schemes on much finer grids [START_REF] Shu | Numerical convergence study of nearly incompressible, inviscid taylor-green vortex flow[END_REF] .

Scalar fields listed in table II are all conserved as expected, which can be seen in a zoomed view of the corresponding curves in Figure 3. The entropy quadrature ρs 2 in the s-cons configuration plotted in Figure 3(c) conserves with a tolerance of 10 -9 , which proves that the scalar operator ∇ C in equation ( 15) is quadrature preserving. The total energy quadrature ρE 2 resolved from configuration E-cons and E-cons-UP are given in Figure 3(d). Although this quantity is only globally conserved 61 , its variation is very small (at a tolerance of 10 -6 ), because the efforts of the pressure work should be very small at low Mach number.

It is also worth noting that the primitive formulation of the entropy equation (Eq. 36d)

equipped with an isotropic spacial operator (s-prim) performs quite well in preserving energy quantities under this nearly incompressible case. It effectively preserves the entropy (ρs) and the associated quadratic invariant ρs 2 . Even the internal energy and its quadratic invariant (ρe, ρe 2 ) are nearly conserved as seen in Fig. 3. The primitive formulation of entropy with MUSCL scheme (s-prim-MU) applied in our previous works [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF][START_REF] Feng | Solid wall and open boundary conditions in hybrid recursive regularized lattice boltzmann method for compressible flows[END_REF] for compressible flows fails to conserve any energy scalar due to the complex numerical treatment involving much larger stencils, but it gives a stable solution.

B. Higher Mach number case

The compressible simulations of TGV are performed at Ma=0.8. At this Mach number, operator ∇ OC leads to stable simulations. The conservative variables are displayed in Fig. 4 using the ratio between their spatial integration at time t and the initial values, i.e.

f f 0 ≡ f (t)dV f (0)dV (38)
The conservative properties of the schemes are the same as in Table II quantities in TGV simulation from different configurations, Ma=0.8.

V. COMPRESSIBLE CASES WITH DISCONTINUITIES

Having confirmed the stability and conservative properties of the proposed schemes, we now assess their performances in the presence of shock waves. To deal with discontinuities, we combine the schemes given in Eqs. ( 14) and ( 15) dynamically using a shock sensor already to enhance stability (at the cost of additional numerical diffusion).

Substituting the corresponding term in the the LBM mass flux scheme Eq. ( 9) will give a second order accuracy of the mass flux evaluation as in equation (A3)

∇ L • (ρu) = ∂ ∂x α ρu α - ∆t 2 ∂ 2 ∂x α x β Π f col αβ + O(∆t 2 ) = ∂ ∂x α ρu α - ∆t 2 1 ∆t ∂ ∂x α [ρu α -(ρu α ) (t + ∆t, x)] + O(∆t) = ∂ ∂x α ρu α + (ρu α ) (t + ∆t, x) 2 + O(∆t 2 ) (A3)
It is worth noting that although equation ( 10) is generally valid for any LBM distributions, there is one essential requirement to make the flux evaluation second order as indicated in the above deductions: the next order moments, e.g. mass flux ρu α here, must be also transported by the same population. In contract, the population f φ introduced in Eq. (12) only servers to transport the scalar field, the scalar flux ρu α φ are not governed by the same population anymore. As a matter of fact, the velocity components in the scalar flux are updated by the population f in practice. Thus, the upwinded divergence operator ∇ U can not achieve second order by time marching.

Appendix B: Taylor expansion of the center scalar operator Thus, although this scheme is first order corresponding to the flow part ( ∂Π f col αβ ∂xα∂x β ), it does not introduce any artificial diffusion corresponds to the scalar field. As there is no second order derivative of φ in this expression.

∇ C • (ρuφ

(D1j)

In the collision operation (eq.( 27)), the projection of the equilibrium distribution on the Hermite basis reads a (0),eq ≡ H a [START_REF] Stolnick | Cfd vision 2030 study: A path to revolutionary computational aerosciences[END_REF],eq α ≡ H

(1) i,α f eq i = ρu α (D2b) a

(2),eq αβ ≡ H

(2) i,αβ f eq i = ρu α u β (D2c) a

(3),eq αβγ ≡ H

(3) i,αβγ f eq i = ρu α u β u γ (D2d) ≡ a (3),eq yyza (3),eq xxz (D2j)

The second order projections of the non-equilibrium part are evaluated as a

(2),neq αβ = H

(2)

iαβ f col i (x -c i ∆t, t -∆t) -f eq i . ( D3
)
The corresponding third order projections are calculated through a recursive procedure a

(3),neq αβγ = u α a neq βγ + u β a neq αγ + u γ a neq αβ .

(D4)

The forcing term F E i in the collision step is added to obtain a correct viscous tensor. In the current study, it is evaluated as described in [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] .

Appendix E: Numerical conservation of mass and momentum in HRR-p

In any periodic calculation domains, the spacial integration of density at next time step evaluated by HRR-p can be written as which is also identical to the corresponding integration in the current time step.

Figure 1 :

 1 Figure 1: Error on density (a), velocity (b) and internal energy (c) as a function of the spatial discretization for the numerical schemes ∇ C (E-cons) and ∇ U (E-cons-UP), for the entropy spot convection test case.

Figure 2 :Figure 3 :

 23 Figure 2: Temporal evolution of (a)ρ, (b) ρu, (c) ρκ, (d) ρs, (e) ρe, (f) ρE, (g) ρs 2 , (h) ρe 2 , (i) ρE 2 quantities in TGV simulation from different configurations at Ma=0.08.

Figure 4 :

 4 Figure 4: Temporal evolution of statistical (a) ρ, (b) ρu, (c) ρs, (d) ρE, (e) ρs 2 , (f) ρE 2

Figure 5 :

 5 Figure 5: Density, velocity, pressure and internal energy evolution (resp. a, b, c, d) in SOD shock tube at t = 0.2 resolved by different schemes. (e, f, g, h) correspond to the same plots, zoomed in the region of interest.

Figure 6 :

 6 Figure 6: Density iso-contours of 2D Riemann case No. 3 in 68 at t = 0.3 obtained with (a) E-cons-dyn, (b) s-cons-dyn and (c) s-prim-MU schemes.

  same as the current time step. The momentum components at next time step x ρu α (x, t + ∆t) =

  equations is considered to emphasize the impact of the flux schemes, i.e.

	∂ρE ∂t	+	∂ρHu α ∂x α	= 0	(36a)
	∂e ∂t	+ u α	∂e ∂x α	+	p ρ	u α,α = 0	(36b)
	∂ρs ∂t	+	∂ρsu α ∂x α	= 0	(36c)
		1 16	(cos 2x + cos 2y) (cos 2z + 2)	(35e)

Numerical configurations tested in the current work are listed in Table

I

. The density and momentum fields are all resolved by the HRR-p LBM flow solver

[START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF]

while the energy scalar transport is modeled using several different solutions: internal energy, total energy and entropy are used in either conservative or primitive form. Only the Euler part of the energy transport

Table I :

 I Nomenclature in Taylor-Green vortex simulations

	Configuration Energy scalar equation	Flux scheme
	E-cons	Equation (36a)	Equation (15)
	e-prim	Equation (36b)	2 ed order isotropic 60
	s-cons	Equation (36c)	Equation (15)
	s-prim	Equation (36d)	2 ed order isotropic 60
	s-prim-MU	Equation (36d)	MUSCL 58
	E-cons-UP	Equation (36a)	Equation (14)
	E-cons-Onishi	Equation (36a)	Equation (C3) 53
	s-cons-Onishi	Equation (36c)	Equation (C3) 53

Table II :

 II Variable conserved locally or globally ( ) by tested numerical schemes

	ρ
	ρu
	ρE
	ρs
	ρE 2 ( )
	ρs 2
	values 61 , i.e.

E-cons e-prim s-cons s-prim s-prim-MU E-cons-UP E-cons-Onishi s-cons-Onishi

  , although the total energy quadratic invariant ρE 2 which is globally conserved in E-cons configuration varies significantly because of the non-negligible pressure work effects at this higher Mach number.

	Nevertheless, the proposed schemes are still stable on this coarse grid, and the conservation
	properties are still observed.
		1.01																			1.2
												1.000
		1.00
																					1.1
												0.998
		0.99
	ρ/ρ0	0.98										0.996 0 ρu/(ρu)								ρs/(ρs) 0 1.0
		0.97										0.994								0.9		E-cons e-prim
																							s-cons
																							s-prim
		0.96		E-cons		s-prim			E-cons-Onishi		0.992		E-cons e-prim		s-prim-MU E-cons-UP			0.8		s-prim-MU E-cons-UP
				e-prim		s-prim-MU		s-cons-Onishi					s-cons		E-cons-Onishi					E-cons-Onishi
				s-cons		E-cons-UP							s-prim		s-cons-Onishi					s-cons-Onishi
		0.95	0	1	2		3	4	5	6	7	0.990	0	1	2	3	4	5	6	7	0	1	2	3	4	5	6	7
							t										t								t
							(a)									(b)								(c)
		1.005											E-cons	s-prim-MU
												2.5		e-prim s-cons	E-cons-UP E-cons-Onishi			1.01
		1.000											s-prim	s-cons-Onishi
																					1.00
												2.0
	0.990 0 ρE/(ρE) 0.995		E-cons							1.5 0 ρs 2 /(ρs 2 )									ρE 2 /(ρE 2 ) 0 0.99
					e-prim s-cons																0.98
		0.985		s-prim s-prim-MU						1.0											E-cons		s-prim-MU
					E-cons-UP																	e-prim		E-cons-UP
					E-cons-Onishi															0.97	s-cons		E-cons-Onishi
					s-cons-Onishi																	s-prim		s-cons-Onishi
		0.980	0	1	2	3	4	5	6	7	0.5	0	1	2	3	4	5	6	7	0	1	2	3	4	5	6	7
							t									t								t
							(d)									(e)								(f)

Acknowledgements

The French Space agency (CNES) is acknowledged for supporting Song Zhao at M2P2. Centre de Calcul Intensif d'Aix-Marseille and GENCI-TGCC/CINES (Grant 2018-A0032A07679) are acknowledged for granting access to their high performance computing resources. We also acknowledge support from Labex MEC (ANR-10-LABX-0092) and the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the "Investissements d'Avenir" and ANR Industrial Chair ALBUMS (grant ANR-18-CHIN-0003-01). The work was performed using the ProLB solver.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

used in our previous works [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] . When a discontinuity is detected by the sensor, the fully upwind operator ∇ U given by Eq. ( 14) is applied, while the centered scheme ∇ C defined in Eq. ( 15) is applied elsewhere. Artificial viscosity obtained by the Jameson shock sensor [START_REF] Jameson | Numerical solution of the euler equations by finite volume methods using runge kutta time stepping schemes[END_REF] is also added to the LBM flow solver for better stabilization [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] . The conservative formulations of total energy and entropy are tested with the dynamic upwinding, denoted by E-consdyn and s-cons-dyn in the following. A one-dimensional shock-tube and a two-dimensional Riemann configuration are simulated.

A. SOD shock-tube case

This classical one-dimensional Sod shock tube case is set initially as (ρ L , u L , p L ) = (1, 0, 1) and (ρ R , u R , p R ) = (0.125, 0, 0.1). The computational domain is discretized by 400 points with an acoustic CFL number around 0.5. The density, velocity and pressure fields at t = 0.2 are shown in Fig. 5. Apart from the numerical methods E-cons-dyn and s-cons-dyn proposed in the current work, the s-prim-MU scheme 58 mentioned in the TGV test case is also examined. A reference solution obtained from an Riemann solver on a much finer grid (2000 points) is also provided. It can be seen that all numerical configurations lead to stable solutions at t = 0.2, however, only the E-cons-dyn scheme gives the correct inner energy profile. This is due to the fact that, under discrete form, only the conservative total energy equation can lead to the appropriate Hugoniot relations through a strong discontinuity [START_REF] Eleuterio | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF] .

B. Two-dimensional Riemann problems

Two-dimensional Riemann problems proposed in [START_REF] Lax | Solution of two-dimensional riemann problems of gas dynamics by positive schemes[END_REF] are considered to investigate the multidimensional performance of the proposed schemes. Configurations No. 3, No. 4, No. 6 and No. 12 in 68 are selected. All these simulations are stable using E-cons-dyn and s-cons-dyn schemes with a grid revolution of 400×400, but as for the SOD test case, only total energy equation leads to correct solutions when strong shocks are encountered. Figure 6 displays the density iso-contours of configuration No.3 resolved by different schemes. The initial condition of this configuration is given in Table III. Compared to the reference solutions reported in [START_REF] Lax | Solution of two-dimensional riemann problems of gas dynamics by positive schemes[END_REF] , only the E-cons-dyn case recovers the correct 2D pattern.

Table III: Riemann 2D conf. 3 68

VI. CONCLUSIONS AND DISCUSSIONS

We have presented a new numerical scheme for use in segregated hybrid LB models, i.e. models consisting of a LB solver for mass and momentum conservation coupled with a FD solver for the macroscopic energy conservation equation (or another scalar). The numerical The result, applied to the HRR-p model proposed [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] , is a numerical method that not only numerically conserves ρ, ρu, but also ρE (and ρE 2 globally), a characteristic which was never obtained, to the authors knowledge, using a hybrid LBM formulation.

Compared with the MUSCL scheme used in our previous work [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] , another advantagebesides leading to adequate jump conditions -is that it only requires the nearest-neighbor information, thus reducing MPI exchanges in an HPC framework.

The scheme can also be adapted easily to take into account external forces work, e.g. gravity, by simply modifying the flux to be computed (ρH in the case of total energy -see Eq. (29a)).

Lastly, its simple implementation allows to dynamically include upwinding where required

Appendix A: Second order accuracy of the LBM mass flux evaluation

To demonstrate the order of accuracy of the mass flux evaluation in LBM, let first expand the momentum transport using equation ( 7)

This equation leads to an expression of the cross derivation of the second order moments of

Appendix C: Conservative extensions of previously existing scalar transport schemes

The original scalar transport proposed by Onishi et al. [START_REF] Onishi | A lattice boltzmann model for polymeric liquids[END_REF] in Eq.n ( 16) can be modified for conservative formulation as

It is exactly the same as the proposed upwind operator ∇ U , because when considering

The corrected low-dissipating counterpart of this scheme introduced in [START_REF] Onishi | A lattice boltzmann model for polymeric liquids[END_REF][START_REF] Ve Küng | Comparison of passive scalar transport models coupled with the lattice boltzmann method[END_REF] in Eq. ( 17) can be expressed as

with ω i the Gaussian weights of the lattice. It is worth noting that the ρ is substituted by j f col j to be more general, thus can be cooperated with the pressure based HRR-p solver.

Finally, in the context of conservative formulations, the scheme proposed by Osmanlic et al. [START_REF] Osmanlic | Lattice boltzmann method for oldroyd-b fluids[END_REF] given in Eq. equation ( 19) can be written with the definition of the mass exchange ∆m i in Eq. ( 18)

i,5 ≡ H

i,xzz -H

i,yyz -H

(3) i,xxz 1998.


This Kélian

Le Brusquet Laurent

Frigerio Adrien

Colas Sébastien

Bondon Pascal

Contribution to continuum estimation in gamma spectrum by observation of local minima

Keywords: gamma spectrometry, continuum estimation, baseline correction method, peak characterization, local minima

This paper presents a method to estimate the continuum of a gamma rays spectrum through the observation of local minima. The method is simple, automatable and has a large field of use. Indeed, it is not limited by the peaks width, and consequently is usable as well with GeHP as with scintillators spectra. As the method exploits signal properties, its operation is easily explainable. It involves a limited set of meaningful parameters for which an adjustment is proposed. The potential of this method is demonstrated through simulations as well as real gamma spectrometry measurement.

I. INTRODUCTION

Gamma spectrometry is a common nuclear measurement technic which can be used for the detection of radioactivity, identification of radionuclides, and quantification of radioactive material. Eventhough other methods exist, in practice, the gamma spectrometry often constitutes the only possible and effective technic. As a consequence, gamma spectrometry has become essential in the nuclear field.

One will find in [START_REF] Gilmore | Practical Gamma-Ray Spectrometry[END_REF] a complete description of gamma rays Physics as well as a number of details on the measurement device. The result of a measurement is a histogram, called spectrum, which spreads detected photons by channels each corresponding to an interval of energy. All spectra have the same structure, that is to say a superposition of a background with peaks specific to some radionuclides, covered by an observation noise. Peaks are mathematically described by a mixture model, usually Gaussian, which contains a great deal of useful information. On the opposite, the background, also called continuum which is rather regular and smooth (at least, in view of the peaks), contains few information.

The purpose of the spectrum analysis is to estimate the mixture parameters from the data. Consequently, continuum is of little interest and the big step of the spectrum analysis will be to isolate the mixture from the background. This operation is known as baseline correction, and the objective of this study is to propose a Baseline Correction Method (BCM) adapted to gamma spectrometry. An ideal method enables the automation of the analysis. Therefore, it shall deal with various peaks shapes and widths, with various radiation detectors technologies, i.e. GeHP (Germanium Hyper Pur, or High Purity GermaniumRadiation) detectors as well as scintillators. The method shall admit a reduced set of parameters independent from the observation. Spectrum analysis is a recurrent topic in gamma spectrometry, and was a highly topical issue in the 70s where many propositions were made [START_REF] Schulze | Investigation of selected baseline removal techniques as candidates for automated implementation[END_REF]. This period brought to light two distinct strategies. On one hand, local analysis relying on the prerequisite estimate of the peaks positions [START_REF] Ma Mariscotti | A method for automatic identification of peaks in the presence of background and its application to spectrum analysis[END_REF] [4] [START_REF] Robertson | An automatic peakextraction technique[END_REF] [6] [START_REF] Connelly | Automatic location and area determination of photopeaks[END_REF] [8] which enables to estimate pieces of continuum beneath the peaks. One the other hand, global analysis which aims to estimate the whole continuum [START_REF] Inouye | Application of fourier transforms to the analysis of spectral data[END_REF] [10] [11] [12] [13] [14] [START_REF] Morháč | An algorithm for determination of peak regions and baseline elimination in spectroscopic data[END_REF] [16] [START_REF] Barat | Nonparametric bayesian inference in nuclear spectrometry[END_REF]. Nowadays, those propositions were naturally ranked by the operating experience, and the local approach based upon the second derivative method [START_REF] Ma Mariscotti | A method for automatic identification of peaks in the presence of background and its application to spectrum analysis[END_REF] is commonly used [START_REF][END_REF] and officially recommended [START_REF]Guidance for gamma spectrometry measurement of radioactive waste[END_REF]. This method chiefly draws its success from its simplicity and explainability. However, it remains difficult to automate, because it may fail in the presence of Compton edges or multiplets, i.e. mixtures of close overlapping peaks. On the other hand, a large pan of proposed global methods involve a model for the continuum (splines, Gaussian processes etc), which introduce an improper regularity prior: continuum often contains discontinuities which are difficult to model. Thus, spectrum analysis is still an active research topic.

The central idea of the present study relies on the following empirical observation: local minima rarely appear on peaks. Thus, it would be possible to estimate the continuum from local minima. Several authors [START_REF] Inouye | Application of fourier transforms to the analysis of spectral data[END_REF] [10] [START_REF] Tj Kennett | Automated analysis for high energy gamma ray spectra[END_REF] have beat around this idea, but the work of Tervo et al. [START_REF] Rj Tervo | An automated background estimation procedure for gamma ray spectra[END_REF] is the most accomplished: it enables to simply estimate the continuum without any prior nor parameters. However, as illustrates the application on two representative spectrum shown in Fig. 1, this estimator only works for thin peaks which quickly limits its use for real applications.

The present paper proposes a new method which covers a much wider range of spectra configurations as shown Fig. 1. The resulting process is in itself simple to describe, but maybe not easy to understand at first glance. The focus is put upon the statistical phenomenon which enables the method to work.

Section II, on one hand, gives a definition of the spectrum. On the other hand, it deducts a number of inherent signal properties on which is built the continuum estimation procedure in section III. Section IV comments the Fig. 1 real spectra application, and section V concludes this work. Confrontation of continuum estimation methods on real spectra. Blue line represents observations. On top is a GeHP spectrum, on the bottom a scintillator spectrum. Both are plotted with a log scale. Red lines represent the reference estimates of Tervo's method. Black lines represent the new method results, applied with parameters t break = 1.5, wopt = 15 and w = 10 for GeHP, w = 80 for scintillator.

II. SPECTRUM SIGNAL PROPERTIES

This section aims at formalizing the problem and derives some general properties about a gamma spectrum and its components.

A. Basic assumptions

Let y y y denote the observed gamma spectrum of n channels such that y y y = (y 1 , ..., y n ). Let m m m = (m 1 , ..., m n ) denote the peaks mixture and c c c = (c 1 , ..., c n ) the continuum. Denoting P the Poisson's distribution, Physics states [START_REF] Gilmore | Practical Gamma-Ray Spectrometry[END_REF] y y y is a sample from a random vector Y Y Y such that:

Y k ∼ P(m k + c k ) = P(µ k )
We introduced here µ µ µ = m m m +c c c, the noiseless signal. Poisson distribution is not practical to handle in literal calculations. Denoting N the normal distribution, the following approximation is possible [START_REF] Saporta | Probabilités, analyse des données et statistique[END_REF]:

Hypothesis 1:

Y k ∼ N (µ k , σ 2 k ) σk 2 = y k
The issue can now be specified: knowing y y y, how to estimate c c c? Because m m m is also unknown, the problem is unsolvable at this state: a constraint is missing. In the paragraphs below, one is looking through the definitions of peaks and continuum for a discrimination criterion that may play the role of the missing constraint.

B. Signal characterization

Let introduce the differential operator ∆x

k = x k -x k-1 .
The continuum is characterized by its low variations. Thus, continuum variations are majorated:

∃β, ∀k, |∆c k | ≤ β
A peak has characteristic areas. A top, at the center, has high values and low variations. Two flanks, uprising and downrising on both sides of the top have high variations, especially in view of continuum variations. Two flats at the borders have low values and low variations. Let denote F the set of all flanks in the spectrum. Then, F contains all tops and flats. The borders of the areas are thereby defined by means of an arbitrary threshold α such that:

β ≤ α, ∀k ∈ F, α ≤ |∆m k |
The unfixed threshold α is a necessary scaling variable, and its choice is a matter of convention. Indeed, what could be considered as a peak in a certain context could be considered as a continuum contribution in another. Fig. 2 shows a mono peak signal with a constant continuum. Consequently 0 is the smallest possible value for β. Choosing α = 50, resulting F areas are represented with grey bands.

As a consequence of the previous definitions, one can derive a lower bound for the variations of the signal:

Property 1:

∀k ∈ F, α -β ≤ |∆µ k | 40

C. Counter variations

Let denote respectively F + and F -the set of increasing flanks and the set of decreasing flanks:

F + = {k ∈ F|α ≤ ∆m k } F -= {k ∈ F|∆m k ≤ -α} ∀k ∈ F let F k
be the probability to have a counter-variation in y y y at k. More explicitly, F k is the probability for y y y to decrease where m m m is increasing, or to increase where m m m is decreasing:

F k = P (∆Y k ≤ 0|k ∈ F + ) = P (0 ≤ ∆Y k |k ∈ F -)
Notice that F k is almost the repartition function of ∆Y k ∼ N (∆µ k , y k + y k-1 ) evaluated at 0. Let denote Φ(.) the cumulative distribution function (CDF) of the standard normal distribution. Thanks to property 1 and noticing Φ is an increasing function, one has an upper bound for F k :

Property 2:

∀k ∈ F, F k ≤ A k = Φ -(α -β) √ y k + y k-1
On Fig. 2 is plotted A A A = (A 1 , ..., A n ) for α = 50 and β = 0. Because this signal is a simulation, one exactly knows the value of ∆m k , which enables to eval the limit admissible values for A A A as follows:

A lim k = Φ -∆m k √ y k + y k-1

One notes through A lim

A lim A lim = (A lim 1 , ..., A lim n ) that countervariations probabilities are close to zero on high variations areas. This observation is confirmed Fig. 3 where the value of A k quickly decreases.

D. Focus on local minima

Let introduce ξ, the set of indexes of y y y local minima:

ξ = {k|y k < y k-1 , y k < y k+1 } (1)
1) Local minima bias: As shown in Fig. 4, the local minima set is biased because local minima's expectation is not equal to the signal expectation. Let ϕ(.) denote the probability density function (PDF) of the standard normal distribution.

Property 3: μξi is an unbiaised estimator of µ ξi such that: where

     μξi = y ξi - C 1 C 0 √ y ξi var( μξi ) = 1 C 0 (C 0 y 2 ξi + (2C 1 -C 0 + C 2 )y ξi -C 1 √ y ξi ) 0
C i = +∞ -∞ u i ϕ(u)(1 -Φ(u)) 2 du
Proof of property 3 is given in the appendix. Note that C 0 is the density of ξ for a stationary signal. Numerical integration results in the following values:

   C 0 = 1/3 C 1 ≈ -0.28209479 C 2 ≈ 0.42522148
2) Occurrence of local minima: One reported on the Fig. 5 the evaluation by simulation of the probability P min that a point of a linear signal, with a slope γ and a noise variance σ, is a local minima. Notes that P min (γ/σ = 0) = C 0 , and P min quickly decreases.

Additionally, one can derive from property 2 an upper bound on the probability that a local minima belongs to F: Property 4:

P (k ∈ F|k ∈ ξ) ≤ A k P (k ∈ ξ)
Proof of property 4 is given in the appendix. One notices P (k ∈ ξ) can not be too small, because in practice there are always a non negligible portion of local minima in a measurement. Moreover, Fig. 5 testifies that ∀k ∈ F, A k is dramatically low. Therefore, P (k ∈ F|k ∈ ξ) is majorated by a constant close to zero, which explains a remarkable phenomenon easily noticeable through data: local minima are absent from the flanks. It is thus possible to identify points in F by observing ξ:

Hypothesis 2: Because it is easy to observe local minima in a given spectrum, hypothesis 2 is a convenient criterion upon which one may build an estimator for the continuum.

ξ ⊂ F

III. CONTINUUM ESTIMATION

A. Intruders filtering

In the previous section, we identified points from F. However, this is not exactly what we were looking for, i.e. c c c points (where m m m is negligible). Some undesirable intruders are present in ξ, as shown in Fig. 6. Indeed, it contains top points which must be removed. Moreover, local minima may accidentally appear on the flank of a significant peak. Anyway, all intruders values are substantially higher than those of the points attached to the continuum. This gives us an opportunity to filter them.

Let define the null hypothesis H0 : there is no discontinuity between ξ i-1 and ξ i . Let t break be the 1 -α/2 order quantile of N (0, 1) and:

z i = |∆y ξi | y ξi + y ξi-1
z i is a z-score for H0, i.e. if t break ≤ z i , one can reject H0 with a confidence α. By selecting a threshold t break for this hypothesis testing, one detects discontinuities in ξ, and forms groups of continuous ξ sets.

Observing the sign of ∆y ξi at the borders between the groups reveals groups which level is higher than those of their direct neighbours. These are intruders groups to be filtered as on Fig. 6.

B. Large peaks issue

Previous outlier filtering is able to deal with GeHP thin peaks spectra. But when faced to peaks acquired by a scintillator, peaks are very large with respect to F variations, and the estimator fails as shown on top plot on Fig. 7. Intruders filtering with t break = 1.5. Blue points represent observations y y y, red crosses represent remaining ξ after intruders filtering, orange points represent intruders.

A simple solution is to subsample the signal before the intruders filtering as shown on the middle plot on Fig. 7. It means that from the considered spectrum, one keep one point out of p, starting at point s. Parameter p is the subsampling step, s the subsampling offset such as 0 ≤ s < p. In this manner, the variation rate between two points is multiplied by p, whereas the noise level has not changed, which enables to fix large peaks issue.

To limit the information loss due to subsampling, one filters subsamples of ξ successively with all possible values of s for a given p, then merge the results as illustrated on the bottom plot of Fig. 7.

Actually, subsampling is a trick which allows to fall back on thin peaks analysis issue when faced to large peaks. An optimum choice for p depends on w, the actual full width at half maximum (FWHM) of the peaks, and w opt , a fix ideal FWHM that one strives to retrieve. This offers an alternative meaningful parametrization for the estimator:

p = max 1, w w opt

C. Full continuum estimation procedure

For a given spectrum, the proposed continuum estimation procedure is the following: 1) Subsample the signal as in part III-B. 2) Observe local minima defined equation 1. 3) Filter intruders as in part III-A. 4) Correct the local minima bias by property 3. 5) Merge remaining points from each subsamples. 6) Apply a linear interpolation to extend the signal to all canals.

Parameters t break and w opt enable to define the balance between peaks and continuum. A small values for t break or w opt will tend to attribute parts of the continuum as peaks, and big values will tend to attribute peaks to the continuum. The mean width of the peaks w is constant for a given measurement device, but still have to be specified (roughly) in order to adapt the method to the spectrometer. An extensive empirical campaign suggests to set t break = 1.5 and w opt = 15. The resulting full estimation is plotted on Fig. 8. One may note that the peak is essentially extracted, even if a little bit underestimated. A part of the basis of the peak is indeed attributed to the continuum which bias the estimation of the area. This is inevitably caused by the bordering areas where it is hard to attribute a point to F or F.

Ideally, one would have filtered the noise before the interpolation, and used a smooth interpolator. This would require to introduce a regularity prior. However, continuum often contains discontinuities which are difficult to take into consideration. Consequently, to find such a prior is an awkward task. For a first approach, one can simply set aside the noise and limit oneself to the padding of the signal with a linear interpolation. in both cases as almost all significant peaks that one may visually identify are extracted. The method seems to conveniently handle GeHP's Compton edges and multiplets where Tervo's method failed. Furthermore, the large peaks width difference between both spectra shows the adaptability of the method. To the knowledge of the authors, there exist no relevant quantitative performance metrics for such estimators, certainly because continuum are by essence unknown quantities. One would rather extend the method to derive peaks positions and areas which are much more convenient to compare.

IV. EXPERIMENTAL DEMONSTRATION

V. CONCLUSION

As a conclusion, this paper derives the properties of local minima under an original aspect, i.e. counter variations, and proposes a totally new estimation procedure, for which key parts are intruders filtering and subsampling.

This method enables a quick, adaptive and automated estimation of the continuum of a gamma spectrum. A full optimization procedure of the parameters, as well as a study of the robustness of the method would be appreciated. Nevertheless, the proposed empirical adjustment should be good enough for a general purpose, and the opinion of the authors is that little performances improvement is to be expected in this direction. Naturally, this method is showing its limitations when continuum variations are close to peaks variations.

As a conclusion, even with ideal settings, one can not ensure the goodness of the estimation, and consequently one should not use this method for quantification purpose. However, the present work stands as an improvement of the method of Tervo et al. [START_REF] Rj Tervo | An automated background estimation procedure for gamma ray spectra[END_REF]. Moreover, it would be possible to derive from this work a peaks detection method.

APPENDIX

Proof of property 3: Let consider a stationary signal of three channels simulated with X 1 , X 2 , X 3 ∼ N (µ, σ 2 ). Let f Xi (.) denote their PDF and g(x) =

x -µ σ . Consequently: f X (x) = g (x)ϕ(g(x)) P (X i ≤ x) = Φ(g(x))

One is looking for f , the conditional PDF of X 2 when it turned out to be a local minima: f (x) =f X2|X2<X1,X2<X3 (x) ∝f X2 (x)P (x < X 1 )P (x < X 3 ) ∝f X2 (x)(1 -P (X 1 < x)) 2 ∝g (x)ϕ(g(x))(1 -Φ(g(x))) 2 Let assess the normalization constant: (σ 2 g(t) 2 + 2µg(t) + µ 2 )f (t)dt

=σ 2 C 2 C 0 + 2µ C 1 C 0 + µ 2 V(Z) = E(Z 2 ) -E(Z) 2 = 1 C 0 (σ 2 C 2 -σC 1 + µ(2C 1 -C 0 ) + µ 2 C 0 )
Proof of property 4: One has P (k ∈ F) = P (k ∈ F + ) + P (k ∈ F -) Furthermore:

P (k ∈ ξ|k ∈ F + ) = P (∆Y k < 0, ∆Y k+1 > 0|α ≤ ∆m k ) ≤ P (∆Y k < 0|α ≤ ∆m k ) = F k ≤ A k
In the same way, one finds P (k ∈ ξ|k ∈ F -) ≤ A k .

Therefore:

P (k ∈ ξ|k ∈ F) = P (k ∈ ξ|k ∈ F + ) P (k ∈ F + ) P (k ∈ F) + P (k ∈ ξ|k ∈ F -) P (k ∈ F -) P (k ∈ F) ≤ A k P (k ∈ F) (P (k ∈ F + ) + P (k ∈ F -)) = A k
Finally:

P (k ∈ F|k ∈ ξ) = P (k ∈ ξ|k ∈ F) P (k ∈ F) P (k ∈ ξ) ≤ A k P (k ∈ ξ)

  Fig.1. Confrontation of continuum estimation methods on real spectra. Blue line represents observations. On top is a GeHP spectrum, on the bottom a scintillator spectrum. Both are plotted with a log scale. Red lines represent the reference estimates of Tervo's method. Black lines represent the new method results, applied with parameters t break = 1.5, wopt = 15 and w = 10 for GeHP, w = 80 for scintillator.

Fig. 2 .

 2 Fig. 2. Identification of peak's areas on a Gaussian example. Dark full line represents m m m + c c c signal, blue points represent observations y y y, grey bands indicates F, orange crosses represent A lim A lim A lim , red crosses represent A A A.

Fig. 3 .

 3 Fig. 3. A k values on a log scale.

Fig. 4 .

 4 Fig. 4. Local minima's expectation bias for µ = 100. Black is associated with the full signal, red is associated with its local minima. Vertical lines indicate expectations of the distributions.

PminFig. 5 .

 5 Fig. 5. Blue points represent P min (γ/σ) on a grid of 20 values of γ linearly spaced on [0, 100] and 20 values of σ spaced evenly on a log scale on [1, 1000]. Black line represents C 0 .

  Fig.6.Intruders filtering with t break = 1.5. Blue points represent observations y y y, red crosses represent remaining ξ after intruders filtering, orange points represent intruders.

Fig. 7 .

 7 Fig. 7. Subsampling effect on outlier filtering with t break = 1.5. Top figure uses no subsamplings, middle figure uses subsampling (p = 5, s = 0), bottom figure uses merged subsamplings (p = 5). Blue line represents observations y y y, red crosses represent remaining ξ after intruders filtering.

Figure 1 Fig. 8 .

 18 Figure1presents real spectra continuum estimations in order to permit a comparison. Undoubtedly, the new estimator has enhanced performances and qualitatively meet our expectations

g

  (t)ϕ(g(t))(1 -Φ(g(t))) 2 dt = +∞ -∞ ϕ(u)(1 -Φ(u)) 2 du = C 0 Thus f (x) = 1 C 0 g (x)ϕ(g(x))(1 -Φ(g(x))) 2 . Let Z arandom variable with f as PDF. Thus, involving variable mutation:

Commissariat à l'énergie atomique et aux énergies alternatives (CEA), Valduc, France

Laboratoire des Signaux et Systèmes (L2S), Université Paris-Saclay, CNRS, CentraleSupélec, Gif-sur-Yvette, France