Christophe Chalons

Rodney O Fox

Frédérique Laurent

Marc Massot

Aymeric Vié

F Laurent § ¶

M Massot § ¶

A Vi

Multivariate Gaussian extended quadrature method of moments for turbulent disperse multiphase flow

Keywords: kinetic equation, multiphase flow, quadrature-based moment methods, kineticbased moment methods, particle trajectory crossing, hyperbolic conservation laws AMS subject classifications. 76T10, 76N15, 35L65, 65D32, 65M08, 76M12, 82C40

1. Introduction. The physics of inertial particles in a carrier fluid (e.g., fluidized beds, sprays, alumina particles in rocket boosters) can be described by a number density function (NDF) satisfying a kinetic equation. Solving such a kinetic equation relies on either a sample of discrete numerical parcels through a Lagrangian Monte-Carlo approach or on a moment approach resulting in a Eulerian system of conservation laws on velocity moments. For the latter, the main difficulty for particles with high Knudsen numbers (i.e., weakly collisional flows) where the velocity distribution can be very far from equilibrium, is the closure of the free-transport term in the kinetic equation.

In the field of high-Knudsen rarefied gases, Grad [START_REF] Grad | On the kinetic theory of rarefied gases[END_REF] introduced a perturbation of the equilibrium distribution using Hermite polynomials whose coefficients are found by transporting additional moments, whereas Levermore introduced a naturally wellposed hierarchy of hyperbolic systems of conservation equations [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF]. More recently, Struchtrup and Torrilhon [START_REF] Struchtrup | Regularization of Grad's 13-moment equations: derivation and linear analysis[END_REF] derived a regularization of the Grad 13-moment system; however, such systems lose the hyperbolic character of the original kinetic equation (see for instance [START_REF] Torrilhon | Characteristic waves and dissipation in the 3-moment-case[END_REF]). To address this issue, methods based on a realizable presumed NDF have been suggested. Torrilhon [START_REF]Hyperbolic moment equations in kinetic gas theory based on multi-variate Pearson-IV distributions[END_REF] used a Pearson-IV distribution as a basis NDF, which extends the hyperbolic region of the 13-moment method, but still is not globally hyperbolic. A class of hyperbolic methods was suggested by Levermore and Morokoff using entropy maximization (EM) [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF]. For moments up to second order, i.e., multivariate Gaussian distributions, the algorithm is fast and easy to implement, and has been used in rarefied gases [START_REF] Mcdonald | Extended fluid-dynamic model for micron-scale flows based on Gaussian moment closure[END_REF][START_REF] Mcdonald | Application of Gaussian moment closure to microscale flows with moving embedded boundaries[END_REF] as well as in particulate flows [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF][START_REF] Sabat | On the development of high order realizable schemes for the Eulerian simulation of disperse phase flows: a convex-state preserving Discontinuous Galerkin method[END_REF][START_REF] Sabat | Fully Eulerian simulation of 3D turbulent particle laden flow based on the Anisotropic Gaussian closure[END_REF][START_REF] Sabat | Eulerian modeling and numerical methods for the description of turbulent polydisperse sprays[END_REF]. For higher-order moments and multi-dimensional systems, McDonald and coworkers [START_REF] Mcdonald | Affordable robust moment closures for CFD based on the maximum-entropy hierarchy[END_REF][START_REF]Towards realizable hyperbolic moment closures for viscous heat-conducting gas flows based on a maximum-entropy distribution[END_REF] propose CFD strategies based on EM, yielding a high level of accuracy but with the drawback of an expensive inversion algorithm; besides, as investigated in [START_REF] Hauck | Convex duality and entropybased moment closures: Characterizing degenerate densities[END_REF], the structure of the moment space leads to some difficulties. This class of methods is referred to as the class of Kinetic-Based Moment Methods.

Another way to proceed is to use Quadrature-Based Moment Methods (QBMM) where the higher-order moments required for closure are evaluated from the lowerorder transported moments using multi-dimensional quadratures. In our previous work, we developed multivariate quadratures in the form of a sum of Dirac delta functions in velocity phase space (see [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF][START_REF] Kah | Prise en compte des aspects polydispersés dans le contexte d'une approche couplée Eulérienne-Lagrangienne pour la modélisation d'un jet de carburant dans les moteurs à combustion interne d'écoulements diphasiques[END_REF] and the references therein). Such an approximation also allows for a well-behaved kinetic numerical scheme in the spirit of [START_REF] De Chaisemartin | Eulerian models and numerical simulation of turbulent dispersion for polydisperse evaporating sprays[END_REF][START_REF] Massot | Eulerian multi-fluid models: modeling and numerical methods[END_REF][START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF] inspired from [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF] where the fluxes in a cell-centered finite-volume formulation are directly evaluated from the knowledge of the quadrature abscissas and weights with guaranteed realizability conditions. Such a quadrature approach has been shown to capture particle trajectory crossing (PTC) in direct-numerical simulations (DNS) where the distribution in the exact kinetic equation remains at all times in the form of a sum of Dirac delta functions [START_REF] Kah | Prise en compte des aspects polydispersés dans le contexte d'une approche couplée Eulérienne-Lagrangienne pour la modélisation d'un jet de carburant dans les moteurs à combustion interne d'écoulements diphasiques[END_REF][START_REF] Kah | Eulerian quadrature-based moment models for polydisperse evaporating sprays[END_REF][START_REF] Fréret | On the role of preferential segregation in flame dynamics in polydisperse evaporating sprays[END_REF][START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF][START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF], but lead to weakly hyperbolic system of conservation equation and can develop artificial singularities when the number of crossing goes beyond the number of Dirac delta function allowed as well as difficulties at the boundary of the moment space [START_REF] Chalons | Beyond pressureless gas dynamics: Quadrature-based velocity moment models[END_REF].

In a turbulent fluid, the effect of the turbulence on the particles through transport and drag leads to dispersion in velocity phase space due to the large number of PTC occurring over a wide range of length and time scales. Nevertheless, classical large scale PTC still occurs for large enough Stokes numbers because the small-scale dispersion is not strong enough to "randomize" the particle velocities resulting from free transport and drag is not strong enough to prevent such large scale trajectory crossing. However, capturing both small scale leading to velocity dispersion as well as large scale PTC requires a large number of quadrature nodes using a delta-function representation. Moreover, such QBMM result in entropic weakly hyperbolic systems of conservation laws and to the formation of δ-shock singularities, the mathematical structure of which is studied in [START_REF] Chalons | Beyond pressureless gas dynamics: Quadrature-based velocity moment models[END_REF]. The purpose of the present contribution is to introduce an extended quadrature-based reconstruction of the NDF for the closure of the free-transport term in the kinetic equation and drag, in between QBMM and KBMM methods, which also allows us to naturally account for velocity dispersion and large scale trajectory crossing. The proposed continuous representation of the NDF allows us both to limit the number of unknowns in multi-dimensional configurations, and to regularize the resulting system of conservation equations, while still being able to capture large scale PTC and velocity dispersion.

One-dimensional (1-D) multivariate Gaussian moment methods were introduced in [START_REF] Chalons | A multi-Gaussian quadrature method of moments for gas-particle flows in a LES framework[END_REF][START_REF] Laurent | A hierarchy of Eulerian models for trajectory crossing in particle-laden turbulent flows over a wide range of Stokes numbers[END_REF] as Multi-Gaussian KBMM. A generalization of these methods is proposed in [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF] and called the extended quadrature method of moments (EQMOM). Here we analyze what we refer to as Gaussian-EQMOM, extend it to 2-and 3-D phase space, and apply it to the solution of the kinetic equation describing inertial particles in a continuous fluid phase. The remainder of the work is organized as follows. In §2 a brief introduction to 1-D Gaussian-EQMOM is provided. In §3, we provide an in-depth description of the application to 1-D kinetic equations and the mathematical properties of two-node Gaussian-EQMOM. In §4, we describe the extension of Gaussian-EQMOM to a 2-D phase space using the conditional quadrature method of moments (CQMOM). In §5 we describe the application of Gaussian-ECQMOM to 2-D kinetic equations. Example applications are provided in §6. Finally, conclusions are drawn in §7.

2. One-dimensional Gaussian-EQMOM. Consider an NDF f (v) defined for v ∈ S where S ⊆ R is a closed support. Let us assume that the first 2N + 1 integer moments of f , defined by

M k := S f (v)v k dv for k ∈ {0, 1, . . . , 2N }, (2.1)
are finite and known. The objective of Multi-Gaussian KBMM [START_REF] Chalons | A multi-Gaussian quadrature method of moments for gas-particle flows in a LES framework[END_REF], which we refer here as EQMOM [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF] is to provide a continuous approximation f a > 0 defined such that

M k = M a k := S f a (v)v k dv for k ∈ {0, 1, . . . , 2N }. (2.2)
More precisely, f a will be assumed to have the form:

f a (v) = N α=1 ρ α δ σ (v, v α ) (2.3)
where the given kernel density function (KDF) δ σ (v, v α ) "tends" to a Dirac delta function δ vα (v) when σ tends to zero and the N weights ρ α , the N abscissas v α and the "spread" parameter σ > 0 are determined from the first 2N + 1 integer moments of f a by (2.2). Let us remark that, for a given type of KDF, it is not always possible to find such parameters. However, if the moment vector is not on the boundary of the moment space, one can always find a set of parameters such that (2.2) holds for k < 2N and M a 2N is as close as possible of M 2N . In the following, we briefly describe Gaussian-EQMOM for S = R and the associated algorithm for computing f a from the moment set {M 0 , . . . , M 2N }.

Definition of Gaussian-EQMOM.

In Gaussian-EQMOM, δ σ (v, v α ) is chosen as a Gaussian PDF of variance σ and centered at v α . The approximate NDF is then

f G (v) := N α=1 ρ α σ √ 2π exp - (v -v α ) 2 2σ 2 (2.4)
and its moments are denoted M G k . In the case σ = 0, f is a weighted sum of n ≤ N delta functions (i.e., an n-point distribution on a finite support). This case can be treated with an adaptive algorithm, described in [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF], using the Chebyshev algorithm [START_REF] Wheeler | Modified moments and Gaussian quadrature[END_REF]. Thus, our focus in this work is on cases with σ > 0, which occur when f is a continuous NDF.

Using the following definitions:

µ k := 1 √ 2π R s k e -s 2 /2 ds, M * k := N α=1 ρ α v k α , (2.5)
(note that µ k = 0 if k is odd, otherwise µ k ≥ 1), the nonnegative integer moments of f G can be written as

M G k = k k1=0 k k 1 σ k-k1 µ k-k1 M * k1 , (2.6)
which is a lower-triangular linear system of the form M G = A(σ)M * where A has unit eigenvalues. Given σ and the moment vector M G , (2.6) can thus be inverted to find the moment vector M * . This observation leads to the following moment-inversion algorithm.

2.2. Moment-inversion algorithm for Gaussian-EQMOM. Given the moment set {M 0 , . . . , M 2N }, the moment-inversion algorithm for Gaussian-EQMOM uses the following iterative procedure:

Starting from σ = 0, 1. compute {M * 0 , . . . , M * 2N -1 } from (2.6) by setting M G k = M k . 2.
Use the Wheeler algorithm [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF] to compute {ρ 1 , . . . , ρ N } and {v 1 , . . . , v N }.

Compute

M * 2N from (2.5). 4. Use (2.6) to compute M G 2N , compare with M 2N . 5. If M G 2N < M 2N , increase σ and repeat from step 1 until M G 2N = M 2N .
In practice, this algorithm will converge if the moment set {M * 0 , . . . , M * 2N -1 } is realizable for the given value of σ used in step 1. Thus, if the moment set becomes unrealizable (which can be determined from the adaptive Wheeler algorithm [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF] in step 2), then the iterations are aborted and the largest value of σ giving realizable moments is used (i.e., the value from the previous iteration). By construction, this moment-inversion algorithm will yield M G k = M k for k ∈ {0, . . . , 2N -1} and M G 2N ≤ M 2N where the equality holds when the iterations are not aborted. In §3, we consider the case with N = 2 where an explicit formula is found for σ, and thus only steps 1 and 2 are required to compute the quadrature parameters.

Evaluating integrals with Gaussian-EQMOM.

In the kinetic equation describing multiphase flows, there are typically terms for collisions or breakage. In the moment transport equations, these terms lead to unclosed integrals of the form

B = R B(v)f (v) dv, C = R 2 C(v 1 , v 2 )f (v 1 )f (v 2 ) dv 1 dv 2 (2.7)
where B and C are the breakage and collision kernels, respectively. Using the Gaussian-EQMOM representation for f , these integrals can be rewritten as

B = N α=1 ρ α 1 √ π R B( √ 2σs + v α )e -s 2 ds, C = N α,β=1 ρ α ρ β 1 π R 2 C( √ 2σs 1 + v α , √ 2σs 2 + v β )e -(s 2 1 +s 2
2 ) ds 1 ds 2 .

(2.8)

The remaining integrals in (2.8) can then be approximated using an M -node Gauss-Hermite quadrature:

B = N α=1 M i=1 ρ α w i B( √ 2σs i + v α ), C = N α,β=1 M i,j=1 ρ α ρ β w i w j C( √ 2σs i + v α , √ 2σs j + v β ) (2.9)
where w i , s i are the corresponding Gauss-Hermite weights and abscissas [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF]. It is important to note that M can be chosen independently from N , and thus that the accuracy of the quadrature approximation of (2.8) by (2.9) does not depend on the number of moments. The dual-quadrature form of f G [START_REF] Lage | On the representation of QMOM as a weighted-residual method -the dualquadrature method of generalized moments[END_REF] is defined by

f G (v) = N α=1 M i=1 ρ αi δ(v - √ 2σs i -v α ).
(2.10)

where ρ αi = ρ α w i .

The reader can note that the Gaussian KDF used in (2.4) could be replaced by any other normalized, symmetric function of s = (v -v α )/σ. The moment-inversion algorithm in §2.2 would remain unchanged, but it would be necessary to find the M -node Gaussian quadrature corresponding to the chosen KDF in order to efficiently evaluate the integrals in §2.3. Alternatively, [START_REF] Cheng | A class of quadrature-based moment-closure methods with application to the Vlasov-Poisson-Fokker-Planck system in the high-field limit[END_REF] uses a two-node B-spline quadrature with a compact support S ⊂ R for which the integral can be evaluated analytically. The B-spline quadrature can be extended to N nodes using the methods described above.

3. Application of 1-D two-node Gaussian-EQMOM to kinetic equations. We first introduce the two-node Gaussian-EQMOM for the NDF f (t, x, v) in 1-D phase/real space for the kinetic equation:

∂ t f + v∂ x f + ∂ v (Af ) = 0, t > 0, x ∈ R, v ∈ R, (3.1)
with initial condition f (0, x, v) = f 0 (x, v). The acceleration A is a real-valued function of v. The exact solution for free transport (when A = 0) is given by f

(t, x, v) = f (0, x -vt, v) = f 0 (x -vt, v).
In this work, we seek an approximation of f (t, x, v) in the form of a two-node Gaussian-EQMOM with weights ρ 1 (t, x) > 0, ρ 2 (t, x) > 0, velocity abscissas v 1 (t, x), v 2 (t, x) and standard deviation σ(t, x) ≥ 0. These five parameters are determined from the 1-D moment transport equations.

3.1. 1-D moment transport equations. Defining the i th -order moment:

M i (t, x) = R f (t, x, v)v i dv, i = 0, . . . , K; K ∈ N;
the associated governing equations are easily obtained from (3.1) after multiplication by v i and integration over v:

∂ t M i + ∂ x M i+1 = A i , i ≥ 0,
where the (unclosed) 1 moment acceleration term is

A i = R iA(v)f (t, x, v)v i-1 dv. (3.2)
For simplicity, we will focus our attention on the five-moment model and its abstract form:

∂ t M 0 + ∂ x M 1 = 0, ∂ t M 1 + ∂ x M 2 = A 1 , ∂ t M 2 + ∂ x M 3 = A 2 , ∂ t M 3 + ∂ x M 4 = A 3 , ∂ t M 4 + ∂ x M 5 = A 4 . =⇒ ∂ t M + ∂ x F(M) = A, (3.3)
with M = (M 0 , . . . , M 4 ) t , F(M) = (M 1 , . . . , M 4 , M 5 ) t and A = (0, A 1 , . . . , A 4 ) t . This model is closed provided that M 5 and A are defined as functions of M. Here we propose to define these functions using two-node Gaussian-EQMOM.

3.2. Two-node Gaussian-EQMOM. The function f G has exact moments M G i of orders i = 0, ..., 5 given by (2.6). The moment closure for (3.3) then naturally consists in setting M 5 = M G 5 where the five unknowns ρ 1 , ρ 2 , v 1 , v 2 and σ are found by solving the nonlinear system M i = M G i , i = 0, . . . , 4; which is clearly equivalent to solving the system

M 0 = ρ 1 + ρ 2 , M 1 = ρ 1 v 1 + ρ 2 v 2 , M 2 -σ 2 M 0 = ρ 1 v 2 1 + ρ 2 v 2 2 , M 3 -3σ 2 M 1 = ρ 1 v 3 1 + ρ 2 v 3 2 , M 4 -6σ 2 M 2 + 3σ 4 M 0 = ρ 1 v 4 1 + ρ 2 v 4 2 .
(3.4)

It remains to prove that this system is well-posed in the following proposition.

Proposition 3.1 (Two-Node Gaussian-EQMOM). For M = (M 0 , M 1 , M 2 , M 3 , M 4 ) t such that M 0 > 0, define e = M 0 M 2 -M 2 1 M 2 0 , q = (M 3 M 2 0 -M 3 1 ) -3M 1 (M 0 M 2 -M 2 1 ) M 3 0 , and η = -3M 4 1 + M 4 M 3 0 -4M 2 0 M 1 M 3 + 6M 0 M 2 1 M 2 M 4 0 . System (3.4
) is well-defined on the phase space Ω given by

Ω = M, M 0 > 0, e > 0, η > e 2 +
q 2 e , and η ≤ 3e 2 if q = 0 .

1 The acceleration terms will be closed if A is affine:

A(t, x, v) = -a(t, x)v + b(t, x
), in which case the moment acceleration term can be written as A k = k(-aM k + bM k-1 ). In gas-particle flows, this limit corresponds to Stokes drag in a stationary fluid.

Setting U = (ρ 1 , ρ 2 , ρ 1 v 1 , ρ 2 v 2 , σ) t , the function U = U(M)
is one-to-one and onto when v 1 = v 2 , and for all v 1 and v 2 provided that we set

ρ 1 = ρ 2 in the case v 1 = v 2 .
Moreover, σ 2 is given by the unique real root in (0, e] of the third-order polynomial

P(σ 0 ) = 2σ 3 0 + (η -3e 2 )σ 0 + q 2 , σ 0 = σ 2 -e. Proof. By setting ρ 1 = ρ1 M0 , ρ 2 = ρ2 M0 , v 1 = v 1 -M1 M0 , v 2 = v 2 -M1 M0 , solving (3.4) is equivalent to solving 1 = ρ 1 + ρ 2 , 0 = ρ 1 v 1 + ρ 2 v 2 , e -σ 2 = ρ 1 v 2 1 + ρ 2 v 2 2 , q = ρ 1 v 3 1 + ρ 2 v 3 2 , η -6σ 2 e + 3σ 4 = ρ 1 v 4 1 + ρ 2 v 4 2 , with e = (M 0 M 2 -M 2 1 )/M 2 0 , q = ((M 3 M 2 0 -M 3 1 ) -3M 1 (M 0 M 2 -M 2 1 ))/M 3 0 , η = (-3M 4 1 + M 4 M 3 0 -4M 2 0 M 1 M 3 + 6M 0 M 2 1 M 2 )/M 4 0 .
Dropping the overlines for the sake of clarity, it is then a matter of uniquely solving the following nonlinear system in (ρ 1 , ρ 2 , v 1 , v 2 , σ 2 ):

ρ 1 + ρ 2 = 1, ρ 1 v 1 + ρ 2 v 2 = 0, ρ 1 v 2 1 + ρ 2 v 2 2 = e -σ 2 , ρ 1 v 3 1 + ρ 2 v 3 2 = q, ρ 1 v 4 1 + ρ 2 v 4 2 = η -6σ 2 e + 3σ 4 . (3.5)
For any given value of σ such that e > σ 2 (or e ≥ σ 2 if q = 0), it is proved in [START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF] that the first four equations allow us to find (ρ 1 , ρ 2 , v 1 , v 2 ). We will then focus on the last equation to find σ 2 .

In the case q = 0, the second and fourth equations yield

ρ 1 v 1 (v 2 1 -v 2
2 ) = 0, and

ρ 2 v 2 (v 2 1 -v 2 2 ) = 0, which gives v := v 1 = -v 2 .
We then get ρ 1 = ρ 2 = 1/2 and σ 2 = e -v 2 , 2v 4 = 3e 2 -η. Recall that our objective is now to uniquely determine v and σ > 0 such that σ 2 ≤ e. A necessary and sufficient condition is then clearly η ∈ (e 2 , 3e 2 ). Note that the case η = e 2 would lead to σ = 0, meaning that the Gaussian functions degenerate into two Dirac delta functions that correspond to the usual quadrature. The case η = 3e 2 gives v = 0 and both Gaussian functions coincide.

In the case q = 0, from (3.5), we observe by using the usual algebra of quadrature methods and by setting

σ 0 = v 1 v 2 and σ 1 = -(v 1 + v 2 ), that e -σ 2 + σ 0 = 0, q + σ 1 (e -σ 2 ) = 0, η -6σ 2 (e -σ 2 ) + 3σ 4 + σ 1 q + σ 0 (e -σ 2 ) = 0.
The last equation then gives that σ 0 = σ 2 -e is a root of the third-order polynomial P(σ 0 ) = 2σ 3 0 + (η -3e 2 )σ 0 + q 2 . Note that one must have σ 0 ∈ (-e, 0) to fulfill the condition e > σ 2 and to be able to reconstruct σ > 0 from σ 0 . First, since lim σ0→-∞ P = -∞, P(0) > 0 and P (σ 0 ) = 12σ 0 , there exists a unique root σ 0 < 0 of P. It then follows that σ 0 > -e if and only if P(-e) < 0, that is if and only if η > e 2 + q 2 /e. This concludes the proof.

3.3.

Calculation of σ 2 . The three roots of P(σ 0 ) can be found analytically (one is real, and two are complex conjugates). The real root yields

σ 2 = e 1 -c 3 + c 1 c 3 (3.6)
where

c 1 = 1 6 η e 2 -3 , c 2 = q 2 4e 3 and c 3 = c 3 1 + c 2 2 1/2 + c 2 1/3
.

Because e, q and η depend only on the moments and not on the weights and abscissas, σ 2 can be computed directly and thus the right-hand side of (3.5) is known.

In transition zones between nondegenerate cases with v 1 = v 2 and the degenerate case with v 1 = v 2 , the two-node Gaussian-EQMOM based on (3.6) can generate quadrature points with high velocity but negligible density. As the CFL condition is based on the maximum velocity, this behavior will drastically reduce the time step. To avoid the generation of high velocities, a regularization procedure is employed. This kind of strategy is often used for moment systems close to singular distributions, see for instance [START_REF] Junk | Domain of definition of Levermore's five-moment system[END_REF][START_REF] Schaerer | On singular closures for the 5-moment system in kinetic gas theory[END_REF]. The difference between the mean velocity of each Gaussian is given by ∆u = |u 2 -u 1 | = q/(e -σ 2 ), and thus is an increasing function of σ 2 . If σ 2 = 0, the velocity difference will be bounded [START_REF] Chalons | Beyond pressureless gas dynamics: Quadrature-based velocity moment models[END_REF]. Regularization consists of using a limiter on σ 2 to control the maximum velocity difference:

σ 2 =    σ 2 from (3.6) if ∆u < ∆u lim e - |q| l if ∆u ≥ ∆u lim (3.7)
where l = ∆u lim + (∆u max -∆u lim ) tanh ∆u -∆u lim ∆u max -∆u lim .

(3.8)

The user-defined parameters ∆u lim and ∆u max > ∆u lim control the velocity difference starting from which we begin to limit σ and the maximum velocity difference, respectively. In practice, these parameters should be close to the local fluid-phase velocity. The difference between ∆u lim and ∆u max permits a smooth transition between regularized and non-regularized zones. Note that in the regularized zones only the fourth-order moment is not satisfied (i.e., M G 4 < M 4 ). For M 0 > 0, the two-node Gaussian-EQMOM algorithm consists of the following three steps:

1. Given moments M in Ω, compute e, q and η.

2. Compute σ 2 from (3.7).

3. Solve (3.5) according to [START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF] to find ρ 1 , ρ 2 , v 1 , v 2 (i.e., variables without the overlines). In the case where

σ 2 = e, set ρ 2 = v 2 = 0 and ρ 1 = M 0 , v 1 = M 1 /M 0 . For M 0 = 0, set ρ 1 = ρ 2 = 0 and (without loss of generality) σ = v 1 = v 2 = 0.
3.4. Effect of the σ 2 limiter on the flux. When using a five-moment method, the moment space can be described through a set of normalized moments defined by Junk [START_REF] Junk | Domain of definition of Levermore's five-moment system[END_REF]. Then, introducing the three central moments e, q, η defined in Proposition 3.1 and a fourth one, occurring in the flux:

s = M 5 M 0 -5 M 4 M 1 M 2 0 + 10 M 3 M 2 1 M 3 0 -10 M 2 M 3 1 M 4 0 + 4 M 5 1 M 5 0 , the normalized moments, corresponding to the normalized NDF f * (w) = √ e M0 f
√ ew + M1 M0 , are (1, 0, 1, q * , η * , s * ) with q * = q e 3/2 , η * = η e 2 and s * = s e 5/2 . In order to observe the effect of the σ 2 limiter on the flux, the normalized flux s * is plotted in Figure 3.4, with and without the use of a limiter. As observed for the closure with entropy maximization [START_REF] Mcdonald | Affordable robust moment closures for CFD based on the maximum-entropy hierarchy[END_REF], without any limiter, the flux diverges quickly towards negative or positive infinity as the line q * = 0 with η * > 3 is approached. However, with the limiter, the flux is smooth and does not diverge. 3.5. Mathematical properties of moment system with two-node Gaussian-EQMOM. The following theorem addresses an important mathematical property of system (3.3), which is its hyperbolicity. This property was shown in the near thermodynamic equilibrium limit, characterized by v 1 ≈ v 2 , by Cheng and Rossmanith [START_REF] Cheng | A class of quadrature-based moment-closure methods with application to the Vlasov-Poisson-Fokker-Planck system in the high-field limit[END_REF]. Here, the general case is addressed in the following theorem.

Theorem 3.2 (Hyperbolicity). Assuming that the vector M = (M 0 , M 1 , M 2 , M 3 , M 4 ) t lives in the space Ω defined in Proposition 3.1, system (3.3) with the two-node Gaussian-EQMOM closure is hyperbolic.

Proof. It is shown in Proposition 3.1 that one can define U = (ρ 1 , ρ 2 , ρ 1 v 1 , ρ 2 v 2 , σ 2 ) t , the vector of the reconstruction variables, using the two-node Gaussian-EQMOM closure. Moreover, σ is not equal to zero since P(-e) = e e 2 + q 2 e -η < 0 with the notation of Proposition 3.1.

The Jacobian matrix of system (3.3) is

J =       0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 α β γ δ       where (α, β, γ, δ, ) = DM 5 DU . DM DU -1
. After some tedious algebra, one obtains

= -2s 1 + v, δ = -s 2 1 -4s 0 -2v 2 + 10σ 2 , γ = -3s 0 s 1 -s 1 v 2 + 2s 0 v + 6σ 2 (2s 1 -v), β = -2s 0 v 2 -3s 2 0 + 3σ 2 (s 2 1 + 2v 2 + 4s 0 ) -15σ 4 , α = s 2 0 v + σ 2 (3s 0 s 1 + s 1 v 2 -2s 0 v) -3σ 4 (2s 1 -v)
where s 0 , s 1 , v and v 2 are defined by

s 0 = v 1 v 2 , s 1 = -(v 1 + v 2 ), v = ρ 1 v 1 + ρ 2 v 2 ρ 1 + ρ 2 and v 2 = ρ 1 v 2 1 + ρ 2 v 2 2 ρ 1 + ρ 2 .
The corresponding characteristic polynomial is

P (X) = -X 5 + X 4 + δ X 3 + γ X 2 + β X + α.
We then want to prove that this polynomial has five distinct roots.

In the case v 1 = v 2 (= v), using the change of variables X = Y + v leads to

P (Y +v) = -Y 5 +10 σ 2 Y 3 -15σ 4 Y = -Y Y 2 -5 + √ 10 σ 2 Y 2 -5 - √ 10 σ 2 ,
and hence there are five distinct roots when σ = 0.

In the case v 1 = v 2 , using the change of variables

X = v + v 2 -v 1 ρ 1 + ρ 2 Z and defining Σ = ρ 1 + ρ 2 v 2 -v 1 σ, leads to P v + v 2 -v 1 ρ 1 + ρ 2 Z = - v 2 -v 1 ρ 1 + ρ 2 5 Q(Z)
with

Q(Z) = Z 5 +2(ρ 2 -ρ 1 )Z 4 + ρ 2 1 -4ρ 1 ρ 2 + ρ 2 2 -10 Σ 2 Z 3 -2(ρ 2 -ρ 1 ) ρ 1 ρ 2 + 6Σ 2 Z 2 + ρ 2 1 ρ 2 2 -3(ρ 2 1 -4ρ 1 ρ 2 + ρ 2 2 )Σ 2 + 15Σ 4 Z + 2(ρ 2 -ρ 1 ) ρ 1 ρ 2 + 3Σ 2 Σ 2 .
To prove that Q(Z) has five distinct real roots, we use Sturm's theorem, which yields the number of distinct real roots in a given interval.

Let us define the Sturm sequence of polynomials: P 0 = Q, P 1 = Q and, for any n ∈ {0, 1, 2, 3}, -P n+2 is the remainder of the Euclidean division of P n+1 by P n . With the use of the formal calculation software Maxima, one can compute this sequence and remark that the coefficient of the highest-order term of each P n is positive for n ∈ {0, 1, 2, 3, 4, 5}, since they are positive coefficient polynomial functions of ρ 1 , ρ 2 and Σ 2 .

Let m(ξ) denote the number of sign changes (ignoring zeroes) in the sequence {P 0 (ξ), P 1 (ξ), . . . , P 5 (ξ)}. For b large enough, one has m(b) = 0 since each P n tends to +∞ as ξ → +∞. Likewise, for a small enough, one has m(a) = 5. Sturm's theorem then shows that Q has five real distinct roots in [a, b]. This concludes the proof.

3.6. Kinetic-based flux. In our numerical implementation to solve (3.3), the spatial fluxes F(M) are computed using a kinetic-based definition:

F i (t, x) = ∞ 0 f (t, x, v)v i+1 dv + 0 -∞ f (t, x, v)v i+1 dv, i = 0, . . . , 4; (3.9)
where the decomposition into positive and negative directions is used to define the flux function as proposed in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF]. The numerical representation of the flux function is a critical point in moment transport methods [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF] because only realizable moment sets can be successfully inverted. Formally, we close (3.9) using the two-node Gaussian-EQMOM:

F i (t, x) = 2 α=1 ρ α v i+1 + α + v i+1 - α , i = 0, . . . , 4; (3.10)
where

v i + α := 1 √ π ∞ -vα √ 2σ (v α + √ 2σs) i e -s 2 ds, v i - α := 1 √ π ∞ vα √ 2σ (v α - √ 2σs) i e -s 2 ds, (3.11)
can be computed analytically. To design a first-order scheme, this decomposition is sufficient as it corresponds to an upwind scheme at the kinetic level. In this case, the transport scheme is shown to be realizable with the following CFL-like condition:

∆t ∆x max α (|v α | + 1.8σ √ 2) ≤ 1 (see Appendix A).
For a quasi-high-order scheme [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF], the spatial fluxes can be found from (3.10) by employing a high-order spatial reconstruction for ρ α and a first-order reconstruction for the abscissas v α and σ. In summary, the numerical fluxes are computed as follows:

1. Given moments M, compute ρ 1 , ρ 2 , v 1 , v 2 , and σ using the two-node Gaussian-EQMOM algorithm in §3.2.

2. Compute kinetic-based fluxes from (3.10) using (3.11).

3. Compute finite-volume numerical moment fluxes as described in [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF]. The case where σ = 0 (i.e., f is composed of two delta functions) is handled by the adaptive Wheeler algorithm [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF].

4. Extension of Gaussian-EQMOM to 2-D phase space. Consider a 2-D phase space with NDF f (v) for v = (u, v) t and define the bivariate moments

M i,j := R 2 f (v)u i v j dv, i, j = 0, . . . , K; K ∈ N. (4.1)
Assuming that these moments exist, in the following we propose a bivariate extension of Gaussian-EQMOM using ideas from CQMOM [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF].

4.1. Definition of 2-D Gaussian-ECQMOM. For clarity, we limit our discussion here to two-node quadrature in each direction of phase space (i.e. a total of four nodes). Nevertheless, the same methodology can be used to develop the formulas for more nodes. For the four-node quadrature, we define an approximate bivariate NDF by

f G 12 (v) := 2 α=1 ρ α g(u; u α , σ 1 )   2 β=1 ρ αβ g(v -V (u); v αβ , σ 2α )   (4.2)
where the Gaussian KDF is

g(u; µ, σ) := 1 σ √ 2π exp - (u -µ) 2 2σ 2 . (4.3)
The form in (4.2) is based on v conditioned on u. An analogous form f G 21 (v) with u conditioned on v is found by permuting u and v. For clarity, we describe the properties of Gaussian extended CQMOM (or Gaussian-ECQMOM) for f G 12 .

The function V (u) in (4.2) replaces the phase-space rotation introduced in [START_REF] Fox | A quadrature-based third-order moment method for dilute gas-particle flow[END_REF] and is defined to have the following properties:

2 α=1 ρ α R V (u)g(u; u α , σ 1 ) du = M 0,1 (4.4)
and

2 α=1 ρ α R uV (u)g(u; u α , σ 1 ) du = M 1,1 . (4.5)
Here, we set 2 V (u) = a 0 + a 1 u where

a 0 = M 2,0 M 0,1 -M 1,0 M 1,1 M 0,0 M 2,0 -M 2 1,0 = µ v -µ u a 1 (4.6)
and

a 1 = M 0,0 M 1,1 -M 1,0 M 0,1 M 0,0 M 2,0 -M 2 1,0 = ρσ v σ u (4.7) with µ u = M 1,0 /M 0,0 , σ 2 u = (M 2,0 /M 0,0 )-µ 2 u , µ v = M 0,1 /M 0,0 , σ 2 v = (M 0,2 /M 0,0 )-µ 2 v and ρ = [(M 1,1 /M 0,0 ) -µ u µ v ]/(σ u σ v ) is the correlation coefficient. Note that V (u)
is well defined if the variance in the first direction σ 2 u is nonzero. In fact V (u) is the conditional expected value of v given u for a bivariate Gaussian distribution. To simplify the notation, we introduce the following definition:

u i V j α := R u i V (u) j g(u; u α , σ 1 ) du. (4.8)
This integral involves Gaussian integer moments up to order i + j, which are known functions of u α and σ 1 .

The integer moments of (4.2) are

M G i,j = 2 α=1 ρ α R u i g(u; u α , σ 1 )   2 β=1 ρ αβ R v j g(v -V (u); v αβ , σ 2α ) dv   du (4.9)
or, by defining v = y + V (u), as

M G i,j = 2 α=1 ρ α R u i g(u; u α , σ 1 )   2 β=1 ρ αβ R [y + V (u)] j g(y; v αβ , σ 2α ) dy   du. (4.

10)

2 A more general choice that uses all of the moments is V (u) = 4 n=0 anu n with an found from the linear system This choice is valid if σ 1 > 0, but care must be taken if the coefficient matrix becomes poorly conditioned. Note that the quantity on the left-hand side in the square brackets is M G i+n,0 , i.e., the univariate integer moments from two-node Gaussian-EQMOM. Thus the coefficient matrix is a Hankel matrix.

For integer j, a binomial expansion leads to

M G i,j = 2 α=1 ρ α j j1=0 j j 1 u i V j-j1 α µ j1 α (4.11)
where we have introduced the conditional moment of (v -V ) j given u = u α defined by

µ j α := 2 β=1
ρ αβ R y j g(y; v αβ , σ 2α ) dy. (4.12)

It follows immediately from the properties of the Gaussian distribution and the definition of ρ αβ that µ 0 α = 1 and µ 1 α = 0. 3 The form of (4.11) leads to the following moment-inversion algorithm.

Moment-inversion algorithm for 2-D Gaussian-ECQMOM.

The first step uses the univariate moments M i,0 with the algorithm in §2.2 for 1-D Gaussian-EQMOM to find {ρ 1 , ρ 2 , u 1 , u 2 , σ 1 }. There are two possible cases: (1) a nondegenerate case with u 1 = u 2 , (2) a degenerate case with u 1 = u 2 . The degenerate case occurs when the univariate moments M i,0 are Gaussian.

Nondegenerate case.

In order to determine the parameters v αβ and σ 2α , we must solve (4.11) to find µ j α for j = 2, 3, 4. This requires six equations, which we define by taking i = 0, 1 and j = 2, 3, 4. It is straightforward to show that the conditional moments can be found from three 2×2 linear systems, which can be solved sequentially:

2 α=1 ρ α µ 2 α = M 0,2 - 2 α=1 ρ α V 2 α , 2 α=1 ρ α u α µ 2 α = M 1,2 - 2 α=1 ρ α uV 2 α ; (4.13) 2 α=1 ρ α µ 3 α = M 0,3 -3 2 α=1 ρ α V α µ 2 α - 2 α=1 ρ α V 3 α , 2 α=1 ρ α u α µ 3 α = M 1,3 -3 2 α=1 ρ α uV α µ 2 α - 2 α=1 ρ α uV 3 α ; (4.14) 2 α=1 ρ α µ 4 α = M 0,4 -4 2 α=1 ρ α V α µ 3 α -6 2 α=1 ρ α V 2 α µ 2 α - 2 α=1 ρ α V 4 α , 2 α=1 ρ α u α µ 4 α = M 1,4 -4 2 α=1 ρ α uV α µ 3 α -6 2 α=1 ρ α uV 2 α µ 2 α - 2 α=1 ρ α uV 4 α .
(4.15) The final step is to use the algorithm in §2.2 for 1-D Gaussian-EQMOM with the moment set {1, 0, µ 2 α , µ 3 α , µ 4 α } to find the quadrature parameters {ρ α1 , ρ α2 , v α1 , v α2 , σ 2α } for α ∈ {1, 2}.

3 µ 0 α = 1 is the unique solution to (4.11) with i = 0, 1 and j = 0. Likewise, µ 1 α = 0 is the solution for i = 0, 1 and j = 1. The latter makes use of the properties of V (u) in (4.4) and (4.5).

Degenerate case.

For the degenerate case, we define a 2-node Gaussian-ECQMOM by

f G 12 (v) = M 0,0 g(u; µ u , σ u ) 2 α=1 ρ α g(v -V (u); v α , σ 2 ) (4.16)
that can further degenerate4 to a bivariate Gaussian with a full covariance matrix when ρ 2 = 0. The moments for this case are

M G i,j = M 0,0 j j1=0 j j 1 u i V j-j1 µ j1 (4.17)
where

u i V j := R u i V (u) j g(u; µ u , σ u ) du (4.18)
with V = M 0,1 /M 0,0 = µ v and uV = M 1,1 /M 0,0 . The conditional moments are defined by

µ j := 2 α=1 ρ α R y j g(y; v α , σ 2 ) dy. (4.19)
By definition, µ 0 = 1 and µ 1 = 0. To complete the quadrature, we find µ 2 , µ 3 and µ 4 from (4.17) using i = 0 and j = 2, 3, 4:

µ 2 = M 0,2 M 0,0 -V 2 (4.20) µ 3 = M 0,3 M 0,0 -V 3 -3 V µ 2 (4.21) µ 4 = M 0,4 M 0,0 -V 4 -6 V 2 µ 2 -4 V µ 3 (4.22)
where the moments V j are known. The final step is to use the algorithm in §2.2 for 1-D Gaussian-EQMOM with the moment set {1, 0, µ 2 , µ 3 , µ 4 } to find the quadrature parameters {ρ 1 , ρ 2 , v 1 , v 2 , σ 2 }.

For the nondegenerate case, the moment-inversion algorithm described above is able to recover 13 of the 16 extended optimal moments5 [10] defined by

      M 0,0 M 0,1 M 0,2 M 0,3 M 0,4 M 1,0 M 1,1 M 1,2 M 1,3 M 1,4 M 2,0 M 2,1 M 3,0 M 3,1 M 4,0 M 4,1       .
The reader can note that the formulas developed in this section for four-node Gaussian-ECQMOM can be extended to N × N nodes in a relatively straightforward manner.

Evaluating integrals with 2-D Gaussian-ECQMOM.

In 2-D phase space, the evaluation of integrals using Gaussian-ECQMOM is very similar to §2.3.

Consider again the unclosed integral

B = R 2 B(v)f (v) dv. (4.23)
Using the 2-D Gaussian-ECQMOM representation f G 12 , this integral can be rewritten as [START_REF]Towards realizable hyperbolic moment closures for viscous heat-conducting gas flows based on a maximum-entropy distribution[END_REF]) where s = (s 1 , s 2 ) t . The remaining integrals in (4.24) can then be approximated using an M -node Gauss-Hermite quadrature:

B = 2 α,β=1 ρ α ρ αβ π R 2 B( √ 2σ 1 s 1 + u α , √ 2σ 2α s 2 + V ( √ 2σ 1 s 1 + u α ) + v αβ )e -s 2 ds (4.
B = 2 α,β=1 ρ α ρ αβ M i,j=1 w i w j B( √ 2σ 1 s i +u α , √ 2σ 2α s j +V ( √ 2σ 1 s i +u α )+v αβ ) (4.25)
where w i and s i are the Gauss-Hermite weights and abscissas [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF]. The extension of (4.24) to a 3-D (or higher) phase space is analogous to (4.25). The dual-quadrature representation is thus

f G 12 (v) = 2 α,β=1 M i,j=1 ρ αβij δ(u - √ 2σ 1 s i -u α )δ(v - √ 2σ 2α s j -V (u) -v αβ ) (4.26)
where ρ αβij = ρ α ρ αβ w i w j .

A second, and equally valid, evaluation of the integrals in (4.23) can be found using f G 21 . In general, if there exist P permutations of the CQMOM conditioning variables (e.g., P = 2 in 2-D and P = 6 in 3-D), then each of them can be used to evaluate the integrals. The arithmetic average value of the P permutations would represent the best overall estimate [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF].

Application of four-node Gaussian-ECQMOM to kinetic equations.

Consider a 2-D velocity phase space with NDF f (t, x, v) for x = (x, y) t and v = (u, v) t that satisfies the kinetic equation

∂ t f + v • ∂ x f + ∂ v • (Af ) = 0, t > 0, x ∈ R 2 , v ∈ R 2 , ( 5.1)
with initial condition f (0, x, v) = f 0 (x, v). The acceleration A = (A x , A y ) t is a real-valued function of v. With a 2-D velocity phase space, we approximate the solution to f using a Gaussian-ECQMOM for the bivariate moments. In this work, we will consider only the minimal Gaussian-ECQMOM that uses four nodes in the 2-D velocity phase space. Nonetheless, the extension to more than four nodes would be analogous to the algorithm presented here.

2-D moment transport equations. Defining the bivariate moments

M i,j (t, x) = R 2 f (t, x, v)u i v j dv, i, j = 0, . . . , K; K ∈ N;
the associated governing equations are easily obtained from (5.1):

∂ t M i,j + ∂ x M i+1,j + ∂ y M i,j+1 = A i,j , i, j ≥ 0;
where the (unclosed)6 moment acceleration term is defined by

A i,j = - R 2 iA x (v)f (t, x, v)u i-1 v j dv - R 2 jA y (v)f (t, x, v)u i v j-1 dv. (5.2)
We will consider in this work a 16-moment model:

∂ t M 0,0 + ∂ x M 1,0 + ∂ y M 0,1 = 0, ∂ t M 1,0 + ∂ x M 2,0 + ∂ y M 1,1 = A 1,0 , ∂ t M 0,1 + ∂ x M 1,1 + ∂ y M 0,2 = A 0,1 , ∂ t M 2,0 + ∂ x M 3,0 + ∂ y M 2,1 = A 2,0 , ∂ t M 1,1 + ∂ x M 2,1 + ∂ y M 1,2 = A 1,1 , ∂ t M 0,2 + ∂ x M 1,2 + ∂ y M 0,3 = A 0,2 , ∂ t M 3,0 + ∂ x M 4,0 + ∂ y M 3,1 = A 3,0 , ∂ t M 2,1 + ∂ x M 3,1 + ∂ y M 2,2 = A 2,1 , ∂ t M 1,2 + ∂ x M 2,2 + ∂ y M 1,3 = A 1,2 , ∂ t M 0,3 + ∂ x M 1,3 + ∂ y M 0,4 = A 0,3 , ∂ t M 4,0 + ∂ x M 5,0 + ∂ y M 4,1 = A 4,0 , ∂ t M 3,1 + ∂ x M 4,1 + ∂ y M 3,2 = A 3,1 , ∂ t M 1,3 + ∂ x M 2,3 + ∂ y M 1,4 = A 1,3 , ∂ t M 0,4 + ∂ x M 1,4 + ∂ y M 0,5 = A 0,4 , ∂ t M 4,1 + ∂ x M 5,1 + ∂ y M 4,2 = A 4,1 , ∂ t M 1,4 + ∂ x M 2,4 + ∂ y M 1,5 = A 1,4 , (5.3)
which requires a closure for the fourth-order moment M 2,2 , the four fifth-order moments M 5,0 , M 3,2 , M 2,3 , M 0,5 , the four sixth-order moments M 5,1 , M 4,2 , M 2,4 , M 1,5 , and the acceleration terms. We propose to define these closures by reconstructing f with four-node Gaussian-ECQMOM. If unclosed, the acceleration term A can be evaluated using the dual-quadrature forms of f G 12 and f G 21 , and taking the arithmetic average. In our numerical examples, Stokes drag is used so that A is closed in terms of the transported moments, and operator splitting is used for the fluxes and the acceleration.

Mathematical properties of 2-D moment system with four-node

Gaussian-ECQMOM . Because we use a dimensional splitting to solve (5.3), let us consider the transport part of the system in the x-direction for the moments used in the Gaussian-ECQMOM reconstruction conditioned on the u velocity:

∂ t M + ∂ x F(M) = 0 (5.4) with M = (M 0,0 , M 1,0 , M 2,0 , M 3,0 , M 4,0 , M 0,1 , M 1,1 , M 0,2 , M 1,2 , M 0,3 , M 1,3 , M 0,4 , M 1,4 ) t
and

F(M) = (M 1,0 , M 2,0 , M 3,0 , M 4,0 , M 5,0 , M 1,1 , M 2,1 , M 1,2 , M 2,2 , M 1,3 , M 2,3 , M 1,4 , M 2,4 ) t .
Let us remark that the moments M 2,1 , M 3,1 and M 4,1 are not considered here since their equations are redundant with those of (5.4). Moreover, for the sake of simplicity, we use the reconstruction with V (u) = 0 in (4.2) to prove the hyperbolicity of system (5.4). Theorem 5.1 (Hyperbolicity). Assuming that the moment-inversion algorithm for 2-D Gaussian-ECQMOM with V (u) = 0 for the vector M is nondegenerate, system (5.4) with this closure is hyperbolic.

Proof. The Jacobian matrix J 2D of the flux is block triangular:

J 2D =       J 0 0 0 0 X 1 A 0 0 0 X 2 0 A 0 0 X 3 0 0 A 0 X 4 0 0 0 A      
where J is the 5 × 5 Jacobian matrix given in Theorem 3.2, corresponding to the 1-D system. The X i are 2 × 5 matrices and A is a 2 × 2 matrix given by

A = 0 1 ν ξ with (ν, ξ) = DM 2,1 D(M 0,1 , M 1,1 ) = DM 2,2 D(M 0,2 , M 1,2 ) = DM 2,3 D(M 0,3 , M 1,3 ) = DM 2,4 D(M 0,4 , M 1,4 )
given by ν = -u 1 u 2 + σ 2 1 and ξ = u 1 + u 2 . From Theorem 3.2, the matrix J is diagonalizable with five distinct eigenvalues. Matrix A is also diagonalizable with the following two eigenvalues:

λ ± = u 1 + u 2 2 ± (u 2 -u 1 ) 2 4 + σ 2 1 .
Moreover, neither of these eigenvalues is an eigenvalue of J. It follows from the block-diagonal structure of the submatrix, found by eliminating the first five rows and columns of J 2D , that J 2D is diagonalizable. This concludes the proof.

5.3. Kinetic-based flux. The moment transport system (5.3) has the form

∂ t M + ∂ x • F(M) = A
with flux vector F = (F x , F y ) t for the 16-moment vector M. In our numerical implementation, the components of the fluxes for moment M i,j are computed using a kinetic-based definition:

F x;i,j = R ∞ 0 f (t, x, v)u i+1 v j du dv + R 0 -∞ f (t, x, v)u i+1 v j du dv, (5.5)
F y;i,j = R ∞ 0 f (t, x, v)u i v j+1 dv du + R 0 -∞ f (t, x, v)u i v j+1 dv du. (5.6)
Here, we describe how F x;i,j is computed by conditioning on the u component of velocity. The treatment of F y;i.j is done analogously by conditioning on the v component of velocity.

Assuming that 1-D Gaussian-EQMOM has been applied to the u direction, we can assume that ρ α , u α and σ 1 are known. Thus, for the nondegenerate case with u 1 = u 2 , (5.5) can be written as

F x;i,j = 2 α=1 ρ α j j1=0 j j 1 u i+1 V j-j1 + α + u i+1 V j-j1 - α µ j1 α (5.7)
where µ j α are the conditional moments of (v -V ) j given u α , and

u i V j + α := 1 √ π ∞ -uα √ 2σ 1 (u α + √ 2σ 1 s) i V (u α + √ 2σ 1 s) j e -s 2 ds, u i V j - α := 1 √ π ∞ uα √ 2σ 1 (u α - √ 2σ 1 s) i V (u α - √ 2σ 1 s) j e -s 2 ds, (5.8)
can be computed analytically. For the degenerate case (ρ 2 = 0), the same expressions are used but with the conditional moments µ j defined in §4.2.2.

Numerical examples.

As example applications, we consider a Riemann problem with 1-D velocity phase space and a frozen turbulence problem with a 2-D velocity phase space. In each case, we solve the moment transport equations in (3.3) and (5.3), respectively.

6.1. 1-D Riemann problem. The initial conditions are defined on the real line with a step in the mean u velocity at x = 0:

U m = M 1 M 0 = 1 if x < 0, -1 otherwise.
For all x, the initial density is unity and the velocity distribution function is Maxwellian with energy σ 2 = 1/3. The velocity distribution is assumed initially to be in equilibrium (i.e., e 0 = σ 2 ). However, the discontinuous nature of the mean particle velocity quickly leads to particle trajectory crossing and a strongly non-equilibrium velocity distribution function.

For 1-D phase space, a measure of the degree of non-equilibrium is the ratio σ 2 /e, which is unity for an equilibrium distribution and zero when the distribution is composed entirely of Dirac delta functions. In order to identify clearly deviations of the higher-order moments from their equilibrium values, we will use the following normalized moments: e * = e/e 0 (energy), q * = q/e 3/2 (skewness), η * = η/e 2 (kurtosis); whose equilibrium values are e * = 1, q * = 0, and η * = 3.

In order to solve the moment equations numerically, the 1-D computational domain -2 < x < 2 is discretized into 402 finite-volume cells. The spatial fluxes are treated using the first-order kinetic-based approach. The time step is chosen based on the largest magnitude of the abscissas v α used to define the spatial fluxes with a CFL number of 0.5. For comparison, results are shown for σ 2 found with (3.6) and with the velocity limiter in (3.7).

Simulation results for the 1-D Riemann problem are presented in Figures 6.1-6.3 at time t = 0.5. Note that due to the equilibrium initial conditions, only one velocity abscissa is used when σ 2 /e = 1 (i.e., v 1 ) and the other (v 2 ) is set to zero automatically using the 1-D quadrature algorithm described in §3.2. At t = 0.5, one can observe from Figure 6.3 that the equilibrium condition is still present on the left and right sides of the computational domain. In the center of the domain, σ 2 /e ≈ 0.2, indicating that the overall distribution is composed of two Gaussian distributions with very little overlap. Also, note that unlike in a pure PTC problem where the velocity abscissas remain at their initial values (i.e., 1 and -1), in Figure 6.2 the abscissas have their largest magnitudes just behind the "shock" in density at the edge of the equilibrium domain. This behavior is a direct result of the definition of the spatial fluxes in terms of the underlying two-node Gaussian-EQMOM distribution. Indeed, the outer tails of the Gaussian distribution have higher velocity than the value at the peak density and thus penetrate faster into the equilibrium domain, resulting in a higher local flux velocity. The strong deviations from equilibrium are also clearly observed in the normalized energy, skewness, and kurtosis in Figure 6.3. Note, however, that using the velocity limiter in (3.7) does not compromise the prediction of the moments in Figure 6.1 because the high velocities correspond to negligibly small weights as can be seen in Figure 6.2.

Except at the edges of the equilibrium domain, we see from Figure 6.1 that the transported moments and σ 2 are smoothly varying functions of x. More importantly, the singularities appearing in the solution do not belong to the class of δ-shocks but to the less singular class of shocks encountered with hyperbolic systems of conservations laws, thus revealing a potentially well-behaved system. Moreover, due to the kinetic-based definition of the spatial fluxes, the moments are always realizable, and the moment-inversion algorithm always computes a well-defined quadrature from the updated moments. Overall, the two-node Gaussian-EQMOM reconstruction of the velocity distribution yields a robust numerical algorithm using a minimum number of moments. In comparison to the high-order delta function reconstruction described in [START_REF]Optimal moment sets for multivariate direct quadrature methods of moments[END_REF], the two-node Gaussian-EQMOM provides a higher fidelity flux representation for a fixed number of transported moments. Moreover, because the moments of the Gaussian-EQMOM distribution can be computed to any desired order, the flux representation described in §3.6 can be systematically improved. This advantage becomes even more significant for 2-D and 3-D phase spaces where the number of transported moments needed for the delta-function reconstruction increases rapidly with the order of the moments [START_REF]Optimal moment sets for multivariate direct quadrature methods of moments[END_REF]. 6.2. 2-D particle-laden turbulent flow. The proposed test case is a particleladen gas phase represented by 2-D frozen homogeneous isotropic turbulence (HIT) generated with the ASPHODELE code of CORIA [START_REF] Reveillon | DNS of Spray Combustion, Dispersion Evaporation and Combustion[END_REF][START_REF] Reveillon | Effects of the preferential segregation of droplets on evaporation and turbulent mixing[END_REF], which solves the 2-D and 3-D low-Mach-number Navier-Stokes equations. The turbulence is generated following the Pope spectrum [START_REF] Pope | Turbulent Flows[END_REF] with parameters p 0 = 4, c L = 0.013, c η = 0.105 and β = 5.2. The particle phase is placed homogeneously in the domain at t = 0 with the same velocity as the gas phase. The droplet number density is uniform and the Stokes number range based on the Kolmogorov time scale (t η = 0.3172 s) is St = τ p /t η ∈ [START_REF] Bouali | Numerical analysis of the influence of two-phase flow mass and heat transfers on n-heptane autoignition[END_REF][START_REF] Laurent | A hierarchy of Eulerian models for trajectory crossing in particle-laden turbulent flows over a wide range of Stokes numbers[END_REF], which is large enough to observe particle trajectory crossings (PTC). Predicting this type of flows is important, as it is expected to exhibit the main effects of a turbulence gas field on a disperse liquid spray, i.e., the preferential concentration of particles in low vorticity zones that greatly influences auto-ignition, which is of primary importance for turbulent combustion applications [START_REF] Bouali | Numerical analysis of the influence of two-phase flow mass and heat transfers on n-heptane autoignition[END_REF]. This effect is highly size dependent, so the test case is particularly interesting for quantifying the accuracy of the proposed methods.

As a reference, Lagrangian computations using a first order in time scheme with time step ∆t = 0.001 s are performed using 10 million particles, for which a satisfactory statistical convergence has been verified. The Eulerian projection of the Lagrangian particle statistics is done using a box filter at the grid-cell size. Example comparisons of the number density field M 0,0 found with the Lagrangian and Eulerian (MG) simulations are shown for four different Stokes numbers in Figures 6.4-6.7, respectively, at t = 4 s. Note that at the time shown in the figures, the Eulerian fields have nearly reached a steady-state condition during which the moments change very little with time. This is possible because of the frozen (time-independent) nature of the turbulent gas phase. At the smallest Stokes number shown in Figure 6.4, very little PTC is present. However, very strong preferential concentration of the particles is evident. In the absence of PTC, the Eulerian moment equations do not require a moment closure since the velocity distribution is monokinetic [START_REF] De Chaisemartin | Eulerian models and numerical simulation of turbulent dispersion for polydisperse evaporating sprays[END_REF] (i.e., the granular temperature is null). Thus, the comparison between the Lagrangian and Eulerian number density fields can be used to judge the ability of the finite-volume scheme to capture the steep density gradients caused by the hypercompressibility of the particle phase. From Figure 6.4, it is clearly evident that the Eulerian method is able to capture accurately the segregation of the particle phase, including regions in the domain where M 0,0 is null (i.e., vacuum zones), with very small numerical diffusion. For St = 5 shown in Figure 6.5, appreciable PTC is present and results in a decrease in the preferential concentration. Indeed, for this Stokes number, the velocity distribution has small (but nonzero) dispersion about the local mean velocity. In the Lagrangian simulation, the full structure of the velocity dispersion is captured, while in the Eulerian model it is only partially captured due to the moment closure needed to close the system. From Figure 6.5, we can observe that the Eulerian model captures much of the segregation structure seen in the Lagrangian number density field. Moreover, due to the hyperbolic nature of the moment closure, the Eulerian number density field is free from δ-shocks caused by not resolving the PTC, which inevitably arise in weakly hyperbolic systems that do not adequately account for velocity dispersion [START_REF] Chalons | Beyond pressureless gas dynamics: Quadrature-based velocity moment models[END_REF]. For St = 10 shown in Figure 6.6, the Lagrangian number density field becomes appreciably more uniform due to abundant PTC. Nevertheless, as clearly observed in the Eulerian field, near-vacuum zones still exist along side regions with preferential concentration of particles. Interestingly, due to the isotropic nature of the gas velocity field, the particle velocity dispersion becomes more isotropic with increasing Stokes number. Under these conditions, the Eulerian simulation becomes less dependent on the choice of the moment closure (albeit as long as the chosen moment closure can represent velocity dispersion). In fact, a full second-order moment closure based on a single Gaussian distribution [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF] (i.e., a degenerate case of Gaussian-ECQMOM) yields results very similar to Figure 6.6 for homogeneous isotropic turbulence. However, for more complicated configurations with large-scale PTC (e.g., crossing particle jets), the full multi-Gaussian represention will be needed to capture the particle phase. The last comparison for St = 20 is shown in Figure 6.7. Here, the Lagrangian number density field is free of vacuum zones since PTC is rampant. Indeed, for this Stokes number the particles are accelerated slowly by the gas and, once they achieve a given velocity, they move in a manner that is weakly coupled to the gas velocity. From the Eulerian field, we can observe that the moment closure again does a good job of reproducing the Lagrangian field. Furthermore, the Eulerian model has the advantage that the 'stochastic noise' present in the Lagrangian field due to the finite sample of particles is absent in the moment model. It would thus be possible to extract highfidelity information concerning, for example, the spatial gradients of moments (e.g., energy spectra) that would be impossible to obtain from a Lagrangian simulation using a tractable number of particles.

From a mathematical perspective, the 2-D simulations presented above, carried out with the kinetic-based fluxes described in §5.3, provide evidence that the proposed multivariate Gaussian-EQMOM closure is robust for simulating particle-laden turbulent flows. In these simulations, no evidence of unrealizable moments was observed (although we have not proved that the kinetic-based fluxes in the form of (5.7) are realizable even for first-order schemes), nor was any evidence of weakly hyperbolic behavior when V = 0 (i.e., outside the conditions of Theorem 5.1) observed.

7. Conclusions. The multivariate Gaussian extended quadrature method of moments and the related moment-inversion algorithms appear to be a very promising approach for the direct-numerical simulation of particle-laden turbulent flows [START_REF]Large-eddy-simulation tools for multiphase flows[END_REF]. The approach combines stability and a lower level of singularity compared to existing quadrature-based moment methods, see [START_REF] Chalons | Beyond pressureless gas dynamics: Quadrature-based velocity moment models[END_REF], and is able to capture both particle trajectory crossing (PTC) caused by the free-transport term and the effects of turbulent agitation. It is noteworthy that the Gaussian-EQMOM naturally degenerates toward the correct velocity distribution with the associated spatial fluxes in both the PTC and dispersion limits. Moreover, by relying on the recent advances in CQMOM [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF], the Gaussian-ECQMOM naturally adapts to the required number of nodes in even highly degenerate cases (e.g., in the absence of particles). As such, the Eulerian moment methods described in this work should offer an attractive alternative to Lagrangian particle tracking methods for simulating particle-laden flows.

where f n j is the two-node Gaussian-EQMOM reconstruction corresponding to the set of moments (M n k,j ) k∈{0,...,4} and given by the parameters (ρ j,α , v j,α , σ j ) α=1,2 . Let us denote by (u λ,α ) α=1,2,3 the three abscissas corresponding to the three-node quadrature of the measure exp(-x 2 )1 ]λ,+∞[ (x)dx. We then introduce the following conjecture, numerically checked for a large number of λ ∈ R:

[C] ∀λ ∈ R, ∀α ∈ {1, 2, 3}, |u λ,α | ≤ max α (u 0,α ) + max{0, λ} . Let us also remark that max α (u 0,α ) ≤ 1.8. The following proposition then gives a sufficient condition for the realizability of the scheme (A.1).

I j,k = M n k,j - ∆t ∆x j R v k |v|f n j (v) dv = R v k 1 -|v| ∆t ∆x j f n j (v) dv
It is sufficient to prove that (I j,k ) k∈{0,...,4} is a moment vector. For that, let us use the reconstruction of f n j and the formulas (3.10,3.11): Noticing that the integrals can be exactly computed by using the three-node quadratures of the measures exp(-x 2 )1 ] -v j,α √ 2σ j

I j,k = 2 α=1 ρ j,α √ π   ∞ -v j,
,+∞[ (x)dx and exp(-

x 2 )1 ] v j,α √ 2σ j
,+∞[ (x)dx, respectively, and using the conjecture concludes the proof.

Fig. 3 . 1 .

 31 Fig. 3.1. Normalized flux s * without (left) and with (right) the limiter on σ 2 .

u

  i+n g(u; uα, σ 1 ) du an = M i,1 for i ∈ {0, 1, 2, 3, 4}.

5 Fig. 6 . 1 .

 561 Fig. 6.1. Solution to 1-D Riemann problem at t = 0.5. Five transported moments (M 0,...,4 ) and reconstructed moment (M 5 ). The dashed line is found with (3.7) and the solid line with (3.6).

Fig. 6 . 2 .

 62 Fig. 6.2. Solution to 1-D Riemann problem at t = 0.5. Left: weights. Center: abscissas. Top right: σ 2 . Bottom right: Gaussian contribution to energy σ 2 /e (bottom right). The dashed line is found with (3.7) and the solid line with (3.6).

Fig. 6 . 3 .

 63 Fig. 6.3. Solution to 1-D Riemann problem at t = 0.5. Top left: density ρ. Top center: mean velocity Um. Top right: normalized energy e * -1. Bottom left: skewness q * . Bottom center: kurtosis η * -3. The dashed line is found with (3.7) and the solid line with (3.6).

Fig. 6 . 4 .

 64 Fig. 6.4. Number density in frozen turbulence with St = 1 at time t = 4 s for Lagrangian (left) and Eulerian (right) simulations.

Fig. 6 . 5 .

 65 Fig. 6.5. Number density in frozen turbulence with St = 5 at time t = 4 s for Lagrangian (left) and Eulerian (right) simulations.

Fig. 6 . 6 .

 66 Fig. 6.6. Number density in frozen turbulence with St = 10 at time t = 4 s for Lagrangian (left) and Eulerian (right) simulations.

Fig. 6 . 7 .

 67 Fig. 6.7. Number density in frozen turbulence with St = 20 at time t = 4 s for Lagrangian (left) and Eulerian (right) simulations.

Proposition A. 1 ( 2 ) ≤ 1 .

 121 Realizability). Let us assume that the conjecture [C] is true. Then the scheme (A.1) with flux (A.2) is realizable if ∀j, ∆t ∆x j max α (|v j,α | + 1.8σ j √Proof. The scheme can be written, for k ∈ {0, . . . , 4}:

If the bivariate moment set is (anisotropic) Gaussian, then the quadrature will degenerate to the exact distribution with w 1 = 1, v 1 = 0 and σ 2 2 = (1 -ρ 2 )σ 2 v .

If we define V (u) = 4 n=0 anu n , then all 16 moments can be recovered.

The acceleration terms will be closed if A is a linear function of the form (au, av) t , in which case the moment acceleration term can be written as A i,j = -a(i + j)M i,j . In gas-particle flows, this limit corresponds to Stokes drag in a stationary fluid.

Acknowledgments. ROF was supported by a grant from the U.S. National Science Foundation (CCF-0830214). The research leading to the results reported in this work has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement No. 246556.

Appendix A. Realizability of the first-order finite volume scheme with the kinetic-based flux.

The finite volume scheme corresponding to the transport part of (3.3) can be written:

where M n k,j is an approximation at time t n of the mean value of the k th -order moment on the cell ]x j-1/2 , x j+1/2 [ of size ∆x j . Using the kinetic-based definition (3.9), the flux is defined by

f n j+1 (v) dv, k = 0, . . . , 4 (A.2)


Guillermo Gallardo

Rutger Fick

William Wells

Rachid Deriche

Demian Wassermann

Groupwise Structural Parcellation of the Cortex: A Sound Approach Based on Logistic Models

Current theories hold that brain function is highly related with long-range physical connections through axonal bundles, namely extrinsic connectivity. However, obtaining a groupwise cortical parcellation based on extrinsic connectivity remains challenging. Current parcellation methods are computationally expensive; need tuning of several parameters or rely on ad-hoc constraints. Furthermore, none of these methods present a model for the cortical extrinsic connectivity. To tackle these problems, we propose a parsimonious model for the extrinsic connectivity and an efficient parcellation technique based on clustering of tractograms. Our technique allows the creation of single subject and groupwise parcellations of the whole cortex. The parcellations obtained with our technique are in agreement with anatomical and functional parcellations in the literature. In particular, the motor and sensory cortex are subdivided in agreement with the human homunculus of Penfield. We illustrate this by comparing our resulting parcels with an anatomical atlas and the motor strip mapping included in the Human Connectome Project data.

Introduction

The human brain is arranged in areas based on criteria such as cytoarchitecture or axonal connectivity. Current hypotheses attribute specialized functions to several of these areas. Hence, parcellating the cortex into such areas and characterizing their interaction is key to understanding brain function. Diffusion MRI (dMRI) enables the in vivo exploration of long-rage physical connections through axonal bundles, namely extrinsic connectivity. Current theories hold that extrinsic connectivity is strongly related to brain function, e.g. this has been shown in macaques [START_REF] Passingham | A survey of methods for analyzing clustered binary response data[END_REF]. Hence, parcellating the cortex based on its extrinsic connectivity can help to understand the internal organization of the brain. However, obtaining a whole-cortex groupwise parcellation based on extrinsic connectivity remains challenging [START_REF] Jbabdi | Long-range connectomics[END_REF]. Current extrinsic connectivity parcellation methods are computationally expensive; need tuning of several parameters or rely on adhoc constraints. For example, Clarkson et al. [START_REF] Clarkson | A Framework For Using Diffusion Weighted Imaging To Improve Cortical Parcellati[END_REF] propose to iteratively refine an anatomical parcellation using information from dMRI. This technique's main drawback is the strong dependence on the initial anatomical parcellation. Lefranc Fig. 1: Lower left corner: graphical model of the linear relationship between a tractogram ( Tsp ); the intra-cluster (˜ c ) and across-subject (˜ s ) variability. We transform the tractograms into a vectorial space while explicitly accounting for the variability, allowing us to propose a clustering technique in accordance. et al. [START_REF] Lefranc | Groupwise connectivity-based parcellation of the whole human cortical surface using watershed-driven dimension reduction[END_REF] calculate the average connectivity profile of regions using a watersheddriven dimension reduction, but they work parcellating predefined gyri. Parisot et al. [START_REF] Parisot | Tractography-Driven Groupwise Multi-scale Parcellation of the Cortex[END_REF] estimate a consistent parcellation across subjects using a spectral clustering approach, without averaging subject's connectivity profiles. Nevertheless, the method needs tuning of several parameters, including the expected number of parcels specified a priori. Moreno-Dominguez et al. [START_REF] Moreno-Dominguez | A hierarchical method for whole-brain connectivity-based parcellation[END_REF] present a parcelling method based on hierarchical clustering parcellation, in which it's not necessary to set an a priori number of clusters. However, they use the cosine distance to compare tractograms and the centroid of tractograms to represent their union. This can lead to an erroneous parcellation since the centroid criterion doesn't minimize the cosine distance between points. These examples show that an efficient groupwise parcelling technique alongside a sound model for the extrinsic connectivity is still needed.

In this work we present a parsimonious model for the cortical connectivity and an efficient parcelling technique based on it, both summarized in Fig. 1. Our model assumes that the cortex is divided in patches of homogeneous extrinsic connectivity. That is, nearby neurons in the cortex share approximately the same long-ranged physical connections, we call this the local coherence criterion.

Our assumption is based on histological results in the macaque brain [START_REF] Schmahmann | Fiber Pathways of the Brain[END_REF]. As in clustered data models in statistics [16] we allow intra-patch variability and across-subject variability in the patches.

Nowadays, the most common tool to estimate the extrinsic connectivity of a point on the cortex in vivo is dMRI-based tractography [START_REF] Jbabdi | Long-range connectomics[END_REF]. To frame tractography within our cortical connectivity model, we use Logistic Random Effects Models [START_REF] Schmahmann | Fiber Pathways of the Brain[END_REF]. This allows us to explicitly denote the relationship between tractography and the variability present in our model. Taking advantage of this, we propose an efficient clustering technique to create single subject and groupwise parcellations of the whole-cortex. Inspired by the method of Moreno-Dominguez et al. [START_REF] Moreno-Dominguez | A hierarchical method for whole-brain connectivity-based parcellation[END_REF], our technique creates a dendrogram: a structure that comprises different levels of granularity for the same parcellation. We also create the dendrogram while imposing the local coherence criterion using only one parameter: the minimum size of each parcel. Then, by choosing cutting criteria, we can explore different parcellation granularities without recomputing the dendrogram.

We validate our technique by taking advantage of the information available in the Human Connectome Project (HCP). Using our technique, we create single-subject and groupwise parcellations for 66 subjects. Then, we compare our purely structural parcellations against an anatomical atlas [5] and responses to functional stimuli [START_REF] Barch | Function in the human connectome: task-fMRI and individual differences in behavior[END_REF]. We show that our parcels subdivide some well-known anatomical structures in accordance with functional responses on the cortex.

Methods

Data and Preprocessing

In this work we used 66 male subjects aged 31-35 from the group S500 of the Human Connectome Project (HCP), all preprocessed with the HCP minimum pipeline [START_REF] Glasser | The minimal preprocessing pipelines for the Human Connectome Project[END_REF]. The main advantages of using this public data base are: each subject possess a dense mesh representing their cortical surface, we use it to create seedpoints for tractography; all the mesh's vertices are coregistered across subjects, property that we use to create the groupwise parcellation; each subject possess the Desikan Atlas [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF] parcellation already computed over their cortical mesh; for each cortical mesh there are also different z-score maps representing the response to different stimuli obtained with functional MRI (fMRI) [START_REF] Barch | Function in the human connectome: task-fMRI and individual differences in behavior[END_REF]. Finally, the group S500 contains z-score maps representing the average functional response to stimuli for 100 unrelated subjects (U100). These studies are used to validate our technique's results.

Cortical Connectivity Model and Tractography

Our model assumes that the cortex is divided in clusters of homogeneous extrinsic connectivity. That is, nearby neurons in the cortex share approximately the same long-ranged physical connections, we call this the local coherence criterion.

Our assumption is based on histological results in the macaque brain [START_REF] Schmahmann | Fiber Pathways of the Brain[END_REF]. As in clustered data models in statistics [16] we allow intra-cluster and across-subject variability. We formalize this concept as:

K = k i=1 K i , ∀ 1≤i,j≤k , i = j → K i ∩ K j = ∅ ∧ conn(K i ) = conn(K j ) (1)
where the set of points on the cortex K is the disjoint union of each cluster K i and conn(•) is the extrinsic connectivity fingerprint of a cluster. We will make the notion of variability explicit in eq. 3. In this work, the connectivity fingerprint of a seed-point in the brain is a binary vector denoting to which other seedpoints it is connected through axonal bundles. This is, the physical connections of a point p ∈ K i in the brain are represented by its connectivity fingerprint conn(p) = conn(K i ). Nowadays, the most common tool for estimating the extrinsic connectivity fingerprint of a point in vivo is probabilistic tractography [START_REF] Jbabdi | Long-range connectomics[END_REF]. Given a seed-point in the brain, probabilistic tractography creates a tractogram: an image where each voxel is valued with its probability of being connected to the seed through axonal bundles. One way of calculating these probabilities is with a Monte Carlo procedure, simulating the random walk of water particles through the white matter [START_REF] Behrens | Characterization and propagation of uncertainty in diffusion-weighted MR imaging[END_REF]. Each one of these paths is known as a streamline. If we think these streamlines as Bernoulli trials, were we get a value for the connection from our seed with other points (1 if they connected by the streamline, 0 if not) [START_REF] Behrens | Characterization and propagation of uncertainty in diffusion-weighted MR imaging[END_REF], then we can model the tractogram of the subject s in the seed-point p as:

T sp = [P ( Cspi = 1)] 1≤i≤n = [θ spi ] 1≤i≤n , Cspi ∼ Bernoulli(θ spi ) , (2)
where Cspi is a Bernoulli random variable3 representing "the point p of the subject s is connected to the voxel i". Each Bernoulli's parameter (θ spi ) represents the probability of being connected, and is estimated as the proportion of success in the Bernoulli trials of each seed.

To formulate the tractogram in accordance to our hypothesis of cortical connectivity, we model a tractogram as a vector of random variables. In our model, each element in a tractogram comes from a random variable depending of the point's cluster alongside its intra-cluster and across-subject variability:

p ∈ K c → Tsp = [P ( Cspi = 1| conn(K c ), ˜ ci , ˜ si )] 1≤i≤n , (3)
in this case, the point p belongs to the cluster c; ˜ ci represents the intra-cluster variability and ˜ si represents the across-subject variability for the connectivity to voxel i in the cluster c.

Since each Cspi follows a Bernoulli distribution (Eq. 2) it's difficult to find an explicit formulation for P ( Cspi = 1| conn(K c ), ˜ ci , ˜ si ) accounting for the variabilities. For this, we use the generalized linear model (GLM) theory. In this theory, the data is assumed to follow a linear form after being transformed with an appropriate link function [START_REF] Mccullagh | Generalized Linear Models[END_REF]. Using the following notation abuse:

logit( Tsp ) [logit(P ( Cspi = 1| conn(K c ), ˜ ci , ˜ si )] 1≤i≤n , (4)
we derive from GLM a logistic random-effects model [16] for each point p:

logit( Tsp ) = β c + ˜ c + ˜ s ∈ R n , ˜ c ∼ N (0, σ 2 c Id), ˜ s ∼ N (0, σ 2 s Id), (5)
where c and s represent the intra-cluster and across-subject variability respectively. According to GLM theory β c ∈ R n is the extrinsic connectivity fingerprint of cluster K c transformed:

logit -1 (β c ) = E( Tsp ) = conn(K c ) . (6)
The choice of logit as link function is based on the work of Pohl et al. [START_REF] Pohl | Using the logarithm of odds to define a vector space on probabilistic atlases[END_REF]. There, they show that the logit function's codomain is a Euclidean space, which allows us to transform and manipulate the tractograms in a well-known space.

Single Subject and Groupwise Parcelling Methodologies

In the previous section we hypothesised that the cortex is divided in clusters with homogeneous extrinsic connectivity, alongside intra-cluster and across-subject variability. In using the previous hypothesis, it is important to remark that we don't have a priori knowledge of the cluster's location or their variability. But, thanks to the proposed logistic random effects model, we formulated the problem of finding these clusters as a well-known clustering problem. This is because, after transforming the tractograms with the logit function as in eq. 4 they will be in a Euclidean space [START_REF] Pohl | Using the logarithm of odds to define a vector space on probabilistic atlases[END_REF]. Even more, eq. 5 states that the transformed tractograms come from a mixture of Gaussian distributions. This is known as a Gaussian mixture model.

To solve the Gaussian mixture model and find the clusters, we use a modified Agglomerative Hierarchical Clustering (AHC) algorithm. This was inspired by the method of Moreno-Dominguez et al. [START_REF] Moreno-Dominguez | A hierarchical method for whole-brain connectivity-based parcellation[END_REF]. To enforce the local coherence criterion we also modify the algorithm to accept one parameter: the minimum size of the resulting clusters. Clusters smaller than this size are merged with neighbors, i.e. physically close clusters in the cortex. As we are working in a Euclidean space, we use the Euclidean distance and the centroid as similarity and linkage functions, improving performance. Our technique's time complexity is O(n 2 log(n)), with n the number of tractograms to cluster [START_REF] Murtagh | Multidimensional Clustering Algorithms[END_REF]. AHC creates a dendrogram: a structure that comprises different levels of granularity for the same parcellation. This allows us to explore different parcellation granularities by choosing cutting criteria, without the need of recomputing each time.

The main advantage of the model we proposed in the previous section is that it allows us to create a groupwise parcellation using linear operations. Assuming direct seed correspondence across subjects, as in the HCP data set, our model lets us remove the subject variability of each seed's tractogram by calculating the expected value across subjects:

E s (g( Tsp )) = E s (β c + ˜ c + ˜ s ), = β c + ˜ c + E s (˜ s ) = β c + ˜ c . ( 7
)
where the last equality is due to E s (˜ s ) = 0 (Eq. 5). This allows us to create population-representative tractograms for each seed free of across-subject variability, which then can be clustered to create a groupwise parcellation. 3 Experiments and Results

Reliability of the Clustering Algorithm for the Model

In the previous sections we presented a model for the cortical extrinsic connectivity and a clustering technique to parcellate the brain. Our technique allows us to create single subject and groupwise parcellations, encoded with different levels of granularity in a dendrogram. However, is not immediate that the chosen clustering algorithm (AHC) solves a Gaussian mixture model (eq. 5) since it was not designed for this particular case [START_REF] Murtagh | Multidimensional Clustering Algorithms[END_REF]. That is, it's not immediate that the algorithm finds the clusters if they are stated as in the model (eq. 5). Now we show that the modified version of the algorithm (sec. 2.3) to enforce the local coherence criterion (sec. 2.2), solves the model for reasonable levels of variability. Moreover, it retrieves the right clusters using one of the simplest criterion to cut the dendrogram: the horizontal cut, i.e. cutting the dendrogram just by choosing the cut's height.

To test the technique, we started by creating synthetic data from the model (eq. 5). We randomly took 10 subjects from our chosen set alongside their extant Desikan parcellations. Then, we created synthetic connectivity fingerprints representing the connections between their Desikan areas. Next, for each vertex in their cortical surface we: replicated those fingerprints; transformed them with the logit function and added cluster-specific variability and across-subject variability as in our model. Finally, we grouped the vertex based on their connectivity using our clustering technique.

If our parcelling technique is able to solve the model, then the Desikan Areas should be encoded in the resulting dendrogram. To show that the Desikan areas were encoded in the resulting dendrogram, we calculated the best obtainable overlap between each Desikan area and the clusters in the dendrogram using Fig. 3: Groupwise dendrogram created by our technique. We can retrieve different granularities of the same parcellation by choosing cutting height, shown as a red dotted line. Motor and Sensory cortex appear at coarse levels of granularity.

the Dice coefficient. To evince the accuracy of the horizontal cut criterion, we compared every obtainable parcellation through cutting our dendrogram in the average subject case against the Desikan atlas using the corrected Rand index [START_REF] Hubert | Comparing Partitions[END_REF]. A Rand index of 1 means that the two parcellations were equal.

Parcelling Subjects From the Human Connectome Project

We next applied our parcellation technique to the HCP data. First, we performed Constrained Spherical Deconvolution (CSD) based tractography [START_REF] Tournier | Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF] from a dense set of points in the cortex. Then, we used our technique to parcellate the cortex by clustering the tractograms. Specifically, since each subject has a surface representing their gray-matter/white-matter interface, we used their vertices as seeds to create tractograms. To avoid superficial cortico-cortical fibers [START_REF] Reveley | Superficial white matter fiber systems impede detection of longrange cortical connections in diffusion MR tractography[END_REF], we shrank each of the 66 surfaces 3mm into the white matter. For each subject, we fitted a CSD model [START_REF] Tournier | Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[END_REF] to their diffusion data using Dipy (version 0.11) [START_REF] Garyfallidis | Dipy, a library for the analysis of diffusion MRI data[END_REF] and created 15000 streamlines per seed-voxel using the implementation of probabilistic tractography in Dipy. Later, we created a tractogram as in (Eq. 2) by calculating the fraction of particles that visited each white-matter voxel. Then, we transformed each tractogram with the logit function [START_REF] Pohl | Using the logarithm of odds to define a vector space on probabilistic atlases[END_REF] as in eq. 4. We clustered the tractograms of each subject using the modified AHC algorithm while imposing a minimum cluster size of 3mm 2 in the finest granularity. Fig. 4: Single subject parcellations created with our technique, alongside both anatomical (left) and functional (right) comparisons from sec. 3.3. Motor and Sensory Cortex (regions 4 and 5) appear to be found and their subdivision is consistent with the motor strip mapping in HCP.

To create the groupwise parcellation, we took advantage of the vertex correspondence across subjects in the HCP data set. Since we are in a vectorial space we calculated the average tractogram of each seed. Then, we created the groupwise parcellation by clustering the average tractograms with our proposed technique (sec. 2.3). The resulting dendrogram for the groupwise case, alongside some of the obtainable parcellations, are in fig. 3.

Functional and Anatomical Comparison

Here we present a proof of concept that our resulting dendrograms encode parcellations with both anatomical and functional meaning by comparing our results against an anatomical atlas and a functional study. To make the comparisons we extracted a parcellation from each subject's dendrogram using the horizontal cut criterion since it showed good accuracy in sec. 3.1. We manually searched for parcellations with a minimum of 36 parcels and a maximum of 150. This was made to get parcellations with coarse granularity while having at least the amount of parcels in the anatomical atlas of Desikan [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF].

Anatomical Comparison. To assess if some anatomical structures were present in the dendrogram and if our resulting parcels were subdividing them, we compared the extracted parcellation with the Desikan atlas [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF]. It's important to Fig. 5: Relation between our pure extrinsic parcellation and the anatomical atlas of Desikan [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF]. Some regions from Desikan projected over a parcellation with less than 150 parcels. Motor and Sensitive cortex appears to be found. remark that each subject in the HCP already has this atlas computed. We projected the Desikan regions over our parcels and then calculated: how many of our parcels were inside each projection and how many of them had more than 90% of its area inside the projection. Comparisons with 8 regions of the Desikan atlas for the single-subject and the groupwise case are in the fig. 4-left and the fig. 5 respectively. The table 1 shows the fraction of parcels inside of each projections that were well contained. Our division of the inferior frontal cortex differs from the Desikan's [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF] (regions 1 and 2 in fig. 5). However, it's similar to that of Anwander et al. [START_REF] Anwander | Connectivity-Based Parcellation of Broca's Area[END_REF]. In particular, our orange and cyan parcels inside regions 1 and 2 correspond with their blue and green parcels in Fig. 4, panels I, II, IV and V [START_REF] Anwander | Connectivity-Based Parcellation of Broca's Area[END_REF]. The Insula; Motor and Sensory cortex (regions 3, 4 and 5 in fig. 5) were well subdivided for both single subject and groupwise parcellations.

Functional Comparison. We studied if our parcels were related to brain function. To do so we calculated their overlap with thresholded maps representing responses to functional stimuli [START_REF] Barch | Function in the human connectome: task-fMRI and individual differences in behavior[END_REF]. These maps are also available for each subject in the HCP. The fig. 4-right shows the functional comparison for two single Fig. 6: Relation between our pure extrinsic parcellation and a functional study [START_REF] Barch | Function in the human connectome: task-fMRI and individual differences in behavior[END_REF]. Some of our groupwise parcels projected over responses to task-related functional activations. Our division of the Motor and Sensitive cortex appears to be consistent with the motor strip mapping. subject parcellations. The fig. 6 shows a projection from our groupwise parcels over averaged tFMRI activations. Each color encodes the response to a different stimulus thresholded with a z-score > 5. The comparison for the average-subject case was done against the average functional responses in the Unrelated100 population from the HCP. Table 2 shows the highest dice coefficient achieved by one of our regions for each task. Here it is important to remark than, in general, each stimulus will generate a functional response in both motor and sensory cortex. Our hypothesis is that this happens because, for example, while a subject is moving the hand he's also feeling it. Therefore, an overlapping of 0.5 suggests that at least one of our parcels in the motor or sensory cortex is having a good overlapping with the functional response, as confirmed visually. Table 2 shows high dice coefficients for each subject, but shows even higher coefficients for the average subject. We hypothesize that averaging tractograms in the Euclidean space is removing the across-subject variability.

Discussion and Conclusion

In this work we presented a parsimonious model for the long-range structural connectivity. Our model assumes that the cortex is divided in patches of homo- geneous extrinsic connectivity, with intra-patch and across-subject variability. Then, using Logistic Random Effects Models we formulate tractography in accordance to our model. This allowed us to transformed the tractograms into a Euclidean space. Working within this sound framework enabled us to easily manipulate and compare tractograms. Taking advantage of this we presented an efficient technique to parcellate the whole cortex in single subject and groupwise cases. The groupwise case assumes seed correspondence across-subject, we expect our model to account errors in seeds coregistration as across-subject variability.

Our technique creates a dendrogram using only one comprehensive parameter: the minimum size of each cluster. Then, different parcellation granularities can be obtained just by choosing at which height to cut the dendrogram. Also, as a preliminary validation of our results we showed that our pure structural parcellation had good agreement with anatomical and functional parcellations.

Fig. 2 :

 2 Fig. 2: Clustering performance for different levels of variability. SN R comprises both intra-cluster and across-subject variability, at smaller SN R more variability is present. Left, best overlap between synthetic regions and clusters in the dendrogram. Right, best corrected rand coefficient for horizontal cutting. The coefficients close to 1, specially in the average case (dark blue line), show the compatibility between our parcellation technique and the proposed model.

Fig. 2 -

 2 left shows the best dice coefficients obtained for every Desikan region under different levels of variability. The Signal-to-Noise-Ratio (SNR) in the figure represents the amount of variability added respect to the original variability of the synthetic connectivity fingerprint. The result in fig. 2-left shows that parcels were retrieved and well-encoded inside the dendrogram for reasonable levels of variability, specially in the average case (dark blue line) were we get ride of the across-subject variability by averaging. Fig. 2-right shows the best obtainable rand index using horizontal cut on the dendrogram under different levels of variability. The high Rand indices obtained show that we can solve our model by simply using the horizontal cut criterion.

Table 1 .

 1 Area Proportion Contained (see fig. 4-left)

	Variability

Table 2 .

 2 Best Dice Coefficient for Task (see fig. 4-right) Face Single Subj. 0.54 ± 0.10 0.46 ± 0.09 0.50 ± 0.1 0.43 ± 0.09 0.56 ± 0.19

	Foot Shape/Average Subj. Hand Tongue Story 0.62 0.60 0.58 0.65, 0.54

Table 1 :

 1 Area proportion contained for parcels inside Desikan projections. Table2: Best dice coefficient obtained for each functional task. Our parcels are mostly contained inside one anatomical region and overlap well with the activations.

For the sake of clarity we denote all random variables with a tilde, e.g. C.

Acknowledgment: This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (ERC Advanced Grant agreement No 694665 : CoBCoM).