Christophe Leclercq

F Gulminelli

T Furuta

O Juillet

Boundary conditions for star matter and other periodic fermionic systems

Keywords: numbers: 21, 60, -n 26, 60, -c 71, 10, Ca

Boundary conditions for star matter and other periodic fermionic systems

I. INTRODUCTION

Interacting fermionic systems in the bulk limit are a standard object of theoretical study in condensed matter physics. Electrons are subject to an external periodic potential in the presence of a cristalline ionic structure and the bulk limit can be seen as an infinite number of spatial replicas of a finite system within a specific geometry [START_REF] Ashcroft | Solid State Physics[END_REF]. In that case, the observables of the bulk system can be obtained from the modelization of one single elementary cell, provided adequate boundary conditions are applied to the many-body wave function. This amounts to introduce a Bloch phase or twist to each wave function in the single-particle basis, and average the twisted observables over the different phases within the first Brillouin zone [START_REF] Lin | [END_REF]3]. In practical applications, this technique has been applied to Quantum Monte-Carlo simulations of electron systems also in the absence of any external periodic potential [START_REF] Lin | [END_REF][3][4]. In this case, the introduction of Bloch phases has to be understood as a technique to accelerate the convergence towards the thermodynamic limit of these very expensive numerical calculations, which would otherwise become prohibitive in computation time. In the absence of an external potential, the periodic cell is just a computation cell with no physical meaning, and independence of the results respect to its size has to be checked.

Coming to the strongly interacting fermionic systems studied in nuclear physics, the bulk limit has not attracted much interest in the community since nuclei are finite. However there are physical situations where the propertis of fermionic matter composed of protons and neutrons in the thermodynamic limit are essential such as for corecollapsing supernova, and for the crust of the neutron stars which are left over by the explosion. This nucleonic stellar matter covers a very wide domain of densities ranging from ρ ≈ 10 8 g • cm -3 to a few times the normal saturation nuclear density ρ ≈ 10 14 g • cm -3 , temperatures between less than 1 and more than 20 MeV, and proton fractions varying between 0. and 0.5 In the sub-saturation density regime, it is well established that matter is charge neutral and mainly composed of neutrons, protons, electrons, positrons and photons in thermal and typically also chemical equilibrium [5,[START_REF] Haensel | Neutron stars: equation of state and structure[END_REF]. Depending on the thermodynamic condition, neutrinos and anti-neutrinos can also participate to the equilibrium.

A very large amount of literature exists on the microscopic modelizations of the neutron star outer and inner crust [START_REF] Sandulescu | [END_REF][8][9][10][11][12]. In this regime at zero temperature and relatively low density, matter consists of a lattice of Wigner-Seitz cells, each cell containing a spherical neutron-rich nucleus immersed in a sea of dilute gas of neutrons and relativistic electrons uniformly distributed inside the cell [13,14]. The linear size of the cell is of several hundreds of fermis in the outer crust. The size decreases going towards the center of the neutron star and the central nucleus becomes heavier and increasingly neutron rich. Typical values [13] in the inner crust range from about 200 particles in a cell of linear size of about 50 fm for ρ ≈ 10 -4 f m -3 to about 1500 particles in a cell of linear size of about 15 fm for ρ ≈ 8 • 10 -2 f m -3 . Most of the existing calculations are based on variational approaches (Hartree-Fock or Hartree-Fock-Bogoliubov) and employ mixed Dirichlet-Neumann boundary conditions at the edge of the cell [13], chosen to produce a flat density close to the cell border and thus to simulate a uniform neutron gas. The specific way of fixing these mixed boundary conditions is not completely clear, and discrepancies due to the choice of fixing these conditions increase with density [15]. A more conceptual problem with Dirichlet or Neumann boundary conditions, is that they neglect antisymmetrization correlations which extend beyond the cell size. This has been pointed out in refs. [16][17][18][19], where Bloch boundary conditions, similar to the ones used for the QMC modelizations of bulk electron systems, have been employed. In these works it was shown that both the neutron specific heat and the motion of unbound neutrons are affected by these effects.

In the intermediate region between crust and core, complex phases are predicted that can break the spherical symmetry of the Wigner-Seitz cell, and can violate the associated translational invariance. This is even more true at finite temperature, where the periodicity of the Wigner-Seitz cell is broken by thermal agitation. These conditions are met in stellar matter in the pre-and post-bounce supernova dynamics, as well as in the cooling process of proto-neutron stars. Even when periodicity cannot rigourously be assumed in these thermodynamic conditions, it can be kept as a practical working hypothesis allowing to address the thermodynamic limit. As already mentioned, the drawback of that is that convergence with respect to the cell size has to be systematically checked. Time-dependent variational microscopic calculations have been proposed to address this region, where statistical averages are calculated from time averages assuming ergodicity [20][21][22][23]. These works have shown that structures, though they can be degenerate in energy [23], are nonetheless approximately periodic in space even at finite temperature.

In these calculations simple periodic boundary conditions are employed, completely neglecting the antisymmetrization correlations beyond the calculation grid. The importance of properly accounting for these correlations was recently stressed in ref. [START_REF] Vantournhout | Nuclear pasta with a touch of quantum : Towards fully antisymmetrised dynamics for bulk fermion systems[END_REF]. In this work, a specific ansatz for the boundary conditions based on the localized Wannier representation has been proposed in the context of antisymmetrized molecular dynamics. It was shown that such boundary conditions allow obtaining a distribution similar to the one of a free Fermi gas, if gaussian wave packets of fixed width are periodically disposed on a two-dimensional grid, while an artificial Pauli potential is needed in order to obtain the same result with classical molecular dynamics. This result implies that, in the inner crust region where stellar matter contains an important component of quasi-free neutrons, properly accounting for the periodic character of the system may be of importance.

To conclude these introductory remarks, it appears that simple Dirichlet, Neumann, or periodic boundary conditions are not adapted to the description of bulk fermionic matter. Different solutions are proposed in the literature for different applications, namely the Twist Averaged Boundary Conditions (TABC) in QMC [START_REF] Lin | [END_REF]3], the Bloch method in mean-field calculations [16][17][18][19], and the Wannier replica method for molecular dynamics approaches [START_REF] Vantournhout | Nuclear pasta with a touch of quantum : Towards fully antisymmetrised dynamics for bulk fermion systems[END_REF], but the equivalence of the different techniques and/or their domain of validity is not completely clear.

In this paper, we will formally develop the link between the different methods and show some model applications for simple non-interacting nuclear systems described through the variational Fermionic molecular Dynamic (FMD) method. We will show that Bloch (or TABC) and the replica method are equivalent when they are applied to the exact eigenstates of the many-body Hamiltonian. When this is not the case, as for variational mean-field theories applied to interacting systems, the equivalence is broken. In this case, the Bloch or TABC technique can produce solutions which differ from the exact results more than if simple periodic boundary conditions are applied. Conversely, the replica method appears more powerful and can easily be applied to any variational based mean-field treatment with an affordable extra computational cost. Our numerical applications will concern one-dimensional systems, for which it has been recently argued [25] that additional drawbacks appear in the TABC method. However it is important to remark that one-dimensional modelizations can often be used in the stellar matter case where (except at very high density close to saturation) spherical symmetry is often a good approximation.

II. THE BLOCH THEOREM AND TWIST AVERAGED BOUNDARY CONDITIONS

We want to address the physical problem of an infinite system (specifically: the neutron star crust) constituted of an infinite number of spatial replicas of finite systems of linear size L (the Wigner-Seitz cells). In the absence of any periodic external potential, for the replicated system to be equivalent to the infinite system, each Wigner-Seitz cell is supposed to contain an integer number of structures (spherical nuclei or more exotic 'pasta' structures [20][21][22][23]) of size l, L i = n i l i , i = x, y, z. The minimal choice is to take a box containing exactly one structure, and in this case the symmetry of the box will follow the symmety of the physical system (a cylindrical box for cylindrical structures, etc).To simplify the notations we will work in one dimension only and write L = nl. The extension to three dimensions is straightforward.

A. The Bloch theorem at the N-body level

We are interested in the translational invariance properties of the Hamiltonian induced by the imposed periodicity. The Hamiltonian of the global system is obviously invariant respect to a simultaneous translation of all particle coordinates of the same arbitrary length r. This invariance physically corresponds to the general statement that the center of mass momentum is a good quantum number. This symmetry has no influence in an infinite system, and we will not consider it further. The fact of working in a finite box of linear size L which is replicated induces additionally an extra non trivial invariance, which we now discuss. The total Hamiltonian of the infinitely replicated system reads

Ĥ = ∞ m=1 H L ˆ x (m) (1)
where x (m) ≡ (x m 1 , . . . , x m N ) denotes coordinates of particles belonging to the m-th replica and the cell Hamiltonian is

ĤL = N i=1 ti + ∞ m=-∞ N i>j v (x i -xj -mL) , (2)
Note that the cell Hamiltonian depends on a finite number N of particles even if these particles are not necessarily confined into the specific volume of the cell. This Hamiltonian is invariant under the translation of any particle coordinate x k of a length mL, with m integer. This invariance can be expressed as

ĤL , Tk (m) = 0 (3)
where Tk (m) = exp -i mL • pk is the L-translational operator of particle k. Translational invariance implies that the eigenfunctions Ψ (x 1 , . . . , x N ) of the cell Hamiltonian eq.( 2) can be written as

Tk (m)Ψ θ k,m = exp (-iθ k,m ) Ψ θ k,m (4)
where θ k,m is the eigenvalue associated to the translation x k → x k -mL. Because of the property of translational operators Tk (m 1 ) Tk (m 2 ) = Tk (m 1 + m 2 ), the eigenvalues satisfy θ k,m1 + θ k,m2 = θ k,m1+m2 . This means that we can associate the translational invariance of periodicity L for particle k with an eigenvalue θ k such that

θ k,m = mθ k (5)
Let us define an auxiliary wave function as

Φ (x 1 , . . . , x N ) = exp -i θ k L x k Ψ (x 1 , . . . , x N ) (6)
Using the fact that Ψ is an eigenfunction of the translation operator eq.( 4), we get

Φ (x 1 , . . . , x N ) = exp -i θ k L (x k -mL) exp (-imθ k ) Ψ (x 1 , . . . , x N ) = Φ (x 1 , . . . , x k -mL, . . . , x N ) (7)
We have shown that, if Ψ is an eigenfunction, then the function Φ defined by eq.( 6) is a periodic function of period L. The same reasoning can be done for any particle k = 1, . . . , N . This means that the eigenfunctions of the translationally invariant Hamiltonian eq.( 2) can be written as

Ψ (x 1 , . . . , x N ) = exp i 1 L N k=1 θ k x k Φ (x 1 , . . . , x N ) ( 8
)
where Φ is a periodic function, that is invariant under the translation mL of any particle coordinate. Since the Hamiltonian eq.( 2) is invariant under the translation of any particle coordinate separately, one could consider in principle a different phase θ k for each particle. However the indistinguishability of nucleons imposes that all the phases must be equal, as we now show [4]. Let us consider a phase θ 1 for the L-translation of particle 1, and a phase θ 2 for the L-translation of particle 2:

Ψ (x 1 + L, x 2 , . . . , x N ) = exp (iθ 1 ) Ψ (x 1 , x 2 , . . . , x N ) (9) Ψ (x 1 , x 2 + L, . . . , x N ) = exp (iθ 2 ) Ψ (x 1 , x 2 , . . . , x N ) (10)
Applying the permutation symmetry to eq.( 9) gives

Ψ (x 1 , x 2 , . . . , x N ) = -exp (-iθ 1 ) Ψ (x 2 , x 1 + L, . . . , x N ) (11)
Applying a translation -L to the second coordinate gives

Ψ (x 1 , x 2 , . . . , x N ) = -exp (-i (θ 1 -θ 2 )) Ψ (x 2 , x 1 , . . . , x N ) (12)
Applying the permutation symmetry once again we get

Ψ (x 1 , x 2 , . . . , x N ) = exp (-i (θ 1 -θ 2 )) Ψ (x 1 , x 2 , . . . , x N ) (13)
which shows that the two phases must be equal, θ 1 = θ 2 = θ L . To be precise we could have θ 1θ 2 = 2nπ with any integer n, but the phase θ L can be taken without any loss of generality in the interval (0, 2π] or (-π, π], that is in the first Brillouin zone. Indeed if we consider a very large number N r of replicas of the WS cell, N r → ∞, the effect of the Bloch phase is negligible and we can write global periodic boundary conditions for the wave function

Ψ (x 1 + N r L, . . . , x N + N r L) = Ψ (x 1 , . . . , x N ) (14)
This implies exp(iθ

L N r ) = 1 or θ L = (2n -N r )π/N r with n integer. Let us write n = n ′ N r + n" with 1 ≤ n" ≤ N r and n ′ integer, then θ L = (2n ′ -1) π + θ (15)
with -π < θ ≤ π. Finally the wave function is

Ψ (x 1 , . . . , x N ) = exp i θ L N i=1 x i Φ (x 1 , . . . , x N ) (16)
This Bloch formulation is a possible way to represent the eigenfunction of the exact N-body Hamiltonian in the cell, accounting for the infinite range antisymmetrization correlations with the replicated system. It amounts to introducing an extra quantum number θ,

Ψ (x 1 , . . . , x N ) = x 1 , . . . , x N | n, θ (17)
where n denotes the ensemble of the other good quantum numbers. The system observables will thus explicitly depend on the phase θ.

The form of the periodic Hamiltonian eq.( 1) implies that the wave function of the global system can be expressed as an antisymmetrized product of N -body Bloch wave functions (16) according to

Ψ tot (x 1 , . . . , x ∞ ) = Â ∞ m=1 Ψ θm (x m 1 , . . . , x m N ) (18)
The Twist Averaged Boundary Condition method consists in assuming that the different angles θ m in the product are all different. Then the physical value of any observable Ô can be computed from eq.( 18) with the extra restriction

θ m = θ n giving Ô tot n = ∞ m=1 n, θ m | Ô| n, θ m (19)
We can see that within this decoherence hypothesis among the different phases, possible non-diagonal terms disappear and the observables per unit cell are obtained as simple averages over phases

Ô L n = 1 2π π -π dθ n, θ| Ô| n, θ (20)
We will explicitly show in the following that the parts of the Fock space corresponding to different phases are indeed disjoint in the case of Slater determinants. For more complex eigenstates, eq.( 20) has to be considered as an approximation which however appears very well verified in practical applications [3].

B. Application to Slater determinants

In the following we look for an ansatz for Φ (x 1 , . . . , x N ) to be introduced as a variational approximation to the full N-body problem. In particular mean field approaches, including Antisymmetrized Molecular Dynamics (AMD) and Fermionic Molecular Dynamics (FMD), are based on a Slater approximation. This means that the variational ansatz can be written as

Φ (x 1 , . . . , x N ) = Â N k=1 u k (x k ) ( 21
)
where  is the antisymmetrization operator among the N particles and u k (x + mL) = u k (x) ∀k. To fulfill this periodicity condition we can write

u k (x) = N ∞ m=-∞ g k (x -mL) ( 22
)
where g is an arbitrary function and N is a normalization factor. In the following we will show applications with the non-orthogonal single-particle FMD/AMD basis set [26,27] given by non-normalized gaussian wave packets

g Z k (x) ≡ exp - 1 2a k (Z k -x) 2 (23)
where Z k , a k are complex variational parameters. For this specific choice the translation of the argument is equivalent to a translation of the gaussian centroid

Z k u k (x) = lim Nr→∞ 1 √ N r Nr/2 m=-Nr/2 g Z k +mL (x) (24)
where the dependence on a k is implicit and omitted to simplify the notations. It is important to remark that the choice of the normalization in eq.( 24) guarantees that the norm of the wave function is finite, which will allow the numerical evaluations below. In the special case where the gaussian width is sufficiently small respect to the cell size, such that the overlaps between gaussians can be neglected this norm is readily evaluated as

< u k |u k >= π|a k |.
Let us consider a generic one-body operator Ô = N k=1 ôk . The matrix element in r-space representation reads

o jk =< u j |ô|u k >= ∞ -∞ dx ∞ -∞ dx ′ u * j (x ′ )o(x ′ , x)u k (x) = lim Nr→∞ Nr 2 L -Nr 2 L dx ∞ -∞ dx ′ u * j (x ′ )o(x ′ , x)u k (x) (25)
Because of the periodicity of the u k this simplifies to

o jk = lim Nr→∞ N r L/2 -L/2 dx ∞ -∞ dx ′ u * j (x ′ )o(x ′ , x)u k (x) (26)
where we have used the fact that the operator is translationally invariant. For a local one-body operator  = N k=1 âk the matrix element simplifies to

a jk = lim Nr→∞ N r L/2 -L/2 dx u * j (x)a(x)u k (x) (27)
This shows that, because of the periodicity, only integrals over one single box are needed. The expectation value of the observable Ô is given by (for the moment we are ignoring the Bloch phase, and considering expectations only over the auxiliary function Φ, which we indicate <> Φ ):

< Ô > Φ =< Φ| Ô|Φ >= N k=1 o kk (28)
if the u k constitute an orthonormal basis, and

< Ô > Φ = < Φ| Ô|Φ > < Φ|Φ > = N j,k=1 o jk B -1 kj ( 29
)
if not. Here B jk is the overlap matrix

B jk =< u j |u k > . (30)
Implementing the ansatz [START_REF] Vantournhout | Nuclear pasta with a touch of quantum : Towards fully antisymmetrised dynamics for bulk fermion systems[END_REF], the matrix element of a generic local operator is calculated by

a jk = lim Nr,N ′ r ,N ′′ r →∞ N ′′ r N r N ′ r Nr/2 m=-Nr/2 N ′ r /2 m ′ =-N ′ r /2 L/2 -L/2 dxg * Zj +mL (x)a(x)g Z k +m ′ L (x) (31)
Because of the finite norm of u k , the double sum is a convergent quantity with increasing number of replicas

lim Nr,N ′ r →∞ Nr/2 m=-Nr/2 N ′ r /2 m ′ =-N ′ r /2 L/2 -L/2 dxg * Zj +mL (x)a(x)g Z k +m ′ L (x) = (32) lim Nr,N ′ r →∞ 1 N r N ′ r Nr/2 m=-Nr/2 N ′ r /2 m ′ =-N ′ r /2 ∞ -∞ dxg * Zj +mL (x)a(x)g Z k +m ′ L (x) (33) =< u k |â|u k > (34)
This means that the limits in eq.( 31) can be performed separately giving:

a jk = L/2 -L/2 dx ∞ m,m ′ =-∞ g * Zj +mL (x)a(x)g Z k +m ′ L (x) (35)
In the hypothesis that the spatial extension of the wave function is negligible beyond a linear size L max = M L, lim x→±Lmax x|g Z = 0, we can write

a jk = L/2 -L/2 dx M m,m ′ =-M g * Zj+mL (x)a(x)g Z k +m ′ L (x) (36)
which is a feasible integral. Eq.( 36) is readily generalized to the case of a non-local operator.

Because of the translational invariance the wave function has to be multiplied by the Bloch phase factor eq.( 16). This means that the expectation values of observables have to be corrected respect to eq.( 28),(29)

< Ô > Ψ = N j,k=1 o jkθ B -1 kjθ = < Ψ θ | Ô|Ψ θ > < Ψ θ |Ψ θ > (37)
where the notation <> Ψ corresponds to a specific choice for the (arbitrary) Bloch phase -π < θ ≤ π,

Ψ θ (x 1 , . . . , x N ) = Â N k=1 ψ kθ (x k ), (38)
ψ kθ (x) = exp (iθx/L) u k (x) (39)
We can finally write the expression of the different matrix elements for the infinite periodic system, for a given choice of the Bloch phase θ. The overlap matrix is not influenced by the Bloch phase, and the same is true for the matrix element of any local operator:

a jkθ = M m,m ′ =-M L/2 -L/2 dxg * Z k +mL (x)a(x)g Z k +m ′ L (x) = a jk (40)
In particular the local density

ρ θ (x) = A j,k=1 < x|ψ k >< ψ j |x > B -1 kj = A j,k=1 < x|u k >< u j |x > B -1 kj = ρ(x)
is independent of θ.

The situation is different for non-local operators, because in this case the Bloch phase does not cancel any more

o jkθ = ∞ m=-∞ M m ′ =-M L/2 -L/2 dx ∞ -∞ dx ′ g * Zj +mL (x ′ )o(x, x ′ )g Z k +m ′ L (x) exp (iθ(x -x ′ )/L) (41)
Let us take the exemple of the momentum operator

kψ j (x) = -i M m=-M exp (iθx/L) ∂ ∂x + iθ L g Zj +mL (x) (42)
The expectation value of the total momentum in the cell K = N j=1 kj is given by

< K > Ψ = N ji=1 < ψ j | k|ψ i > B -1 ij = N ji=1 < u j | k|u i > B -1 ij + N θ L =< K > Φ +N θ L (43)
We can see that the Bloch phase physically represents a boost to the center of mass motion. It will cancel when performing the average between the different phases -π < θ ≤ π (TABC), but this will give a finite contribution to the kinetic energy. Indeed if we consider the square momentum

k2 ψ j (x) = -exp (iθx/L) M m=-M ∂ 2 ∂x 2 + 2iθ L ∂ ∂x - θ 2 L 2 g Zj +mL (x) (44)
This gives rise to a total kinetic energy expectation value

< ÊK > Ψ = 2 2m N jk=1 < ψ j | k2 |ψ k > B -1 kj + + 2 2m 2θ L < K > Φ + 2 2m N θ 2 L 2 = < ÊK > Φ + 2 θ mL < K > Φ + N 2m 2 θ L 2 (45)
This extra kinetic energy term linked to the Bloch phase explicitly enters in the energy variation (or in the equations of motion, in the case of dynamical models [22,23,26,27]). If we consider the typical size of a Wigner-Seitz cell as calculated by Negele and Vautherin [13], in the inner crust (L ≈ 15 fm) the phase effect gives rise to an energy contribution ∆E ≈ 11 MeV, which is far from being negligible if we compare to the corresponding Fermi energy E F ≈ 24 MeV [16][17][18].

C. One-dimensional model cases

As a first model case, we consider one-dimensional free particles. Even if this system is clearly very far from the correlated dishomogeneous solutions relevant for stellar matter, it has the advantage of being exactly solvable. Moreover, it is an especially interesting test case for molecular dynamics models [START_REF] Vantournhout | Nuclear pasta with a touch of quantum : Towards fully antisymmetrised dynamics for bulk fermion systems[END_REF]26,27] and more general models employing localized wave functions [23]: these models are optimized to treat density fluctuations but not necessarily adapted to treat the coupling to the continuum which is needed for the extreme neutron-proton ratio which is found in the inner crust of neutron stars.

In one dimension, the energy per particle of a degenerate ideal Fermi gas at density ρ = N/L = kF π is given by e F G = 2 π 2 ρ 2 6m . In the case of a finite system of N particles within a length L with periodic boundary conditions, the energy levels are given by

E i = 2 2m 2πn i L 2 , ( 46
)
where n i is an integer within the interval [-N 2 , N 2 ]. Due to this level degeneracy, the total cell momentum K L is equal to ± πN L for an even number of particles, while it is zero for an odd number. Because of the periodicity of the system, it is easier to work in Fourier space. All one-body observables can be expressed as a function of the one-body density in momentum space ρ

Ψ (k, k ′ ) = N i,j=1 k|ψ i ψ j |k ′ B -1 ij
, where the single-particle states are given in the momentum representation. These states can be obtained by taking the Fourier transform of the wave function. Considering that a periodic function can always be expressed as a Fourier series,

u n (x) = ∞ l=-∞ c nl e (-2πilx/L) (47)
we get

ψ nθ (k) = √ 2π ∞ l=-∞ c nl δ (θ/L -k -2πl/L) = lim Nr→∞ ∞ l=-∞ Nr -1 2 m=-Nr -1 2 2π N r 1 L L 2 -L 2 dx g Zn+mL (x)e 2iπlx/L δ (θ/L -k -2πl/L) = lim Nr→∞ ∞ l=-∞ 2π N r 1 L Nr L 2 -Nr L 2 dx g Zn (x)e 2iπlx/L δ (θ/L -k -2πl/L) = 2π L a n N r ∞ l=-∞ e -an 2 ( 2πl L ) 2 e i 2πl L Zn δ (θ/L -k -2πl/L) (48)
At variance with the spatial density, the momentum density thus explicitly depends on the value of the phase θ = Lk 0 :

ρ θ (k) = lim Nr→∞ 4π 2 N r L 2 N n,p=1 B -1 n,p a n a * p e - an +a * p 2 ( 2πl L ) 2 e i 2πl L (Zn-Z * p ) δ (θ/L -k -2πl/L) (49)
The momentum density is obtained by averaging over the different phase values, giving

ρ Ψ (k) = lim Nr→∞ 2π N r L Nr 2 n=-Nr -1 2 N n,p=1 B -1 n,p a n a * p e - an +a * p 2 ( 2πl L ) 2 e i 2πl L (Zn-Z * p ) (50)
We can see that the value of ρ in k is solely determined by the wave function Ψ corresponding to the phase value θ = k -2π m , which corresponds to a unique value in the first Brillouin zone. This means that the number of points chosen for the discretization of the integral over the Bloch phase defines the resolution of the momentum density: the Bloch phase generates new plane waves which 'fill up' the distribution as needed to obtain a continuous Fermi sea starting from the discrete levels of a finite system.

The total energy and momentum easily follow

Ê Ψ = 2 2m ∞ -∞ dkρ Ψ (k)k 2 , K Ψ = ∞ -∞ dkρ Ψ (k)k. ( 51
)
The result for the FMD model is shown in Fig. 1. We have chosen a particle density ρ = 1.37 π . The size L of the elementary cell is directly proportional to the number of particles in the cell according to L = N ρ . The FMD ground state is obtained minimizing the total energy with respect to the variational parameters z = {Z n , a n , n = 1, . . . , N }. This can be obtained using the gradient method, or alternatively by the numerical solution of the FMD equations of motion [26]

d Ĥ Ψ dz i = j C * ij dz * i dt , (52)
C ij = i ∂ 2 ∂z * i ∂z j ln Ψ|Ψ (53)
It is possible to show [26] that if a friction term is added as a multiplicative factor on the r.h.s. of eq.( 52)

d Ĥ Ψ dz i = A * j C * ij dz * i dt , (54)
it gives a dynamics of the variational parameters leading to a decrease in time of the total energy. The resulting dissipative dynamics can be written as:

d Ĥ Ψ dt = 2iℑ(A) ij C ij dz * i dt . dz j dt (55)
and this equation is numerically followed until convergence. As expected, Figure 1 shows that the periodicity in k-space is reproduced by the calculation. This is due to the implicit dependence of the expection values P Φ , Ĥ Φ on the choice of the Bloch quantum number k 0 = θ/L obtained through the variational procedure.

0 0.1 0.2 0.3 0.4 0.5 -3 -2 -1 0 1 2 3 ρ[f m -1 ] x[f m] 0 0.2 0.4 0.6 0.8 1 1.2 -1 0 1 <k| ρ|k> k F k 0 N k F -1 0 1 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 -6 -4 -2 0 2 4 6 K Ψ k F K Φ k F k 0 N k F 11.5 13
The exact results of a free Fermi gas are recovered in the model. This means that the FMD ansatz for the wave function appears adapted to describe the free particle system. This was not a-priori obvious because the choice of single-particle wave packets represents a restriction of the complete Slater determinant variational space, and the plane wave is obtained only as a mathematical limit a → ∞ of the wave function, which is not exactly accessible in numerical simulations. This implies that the FMD model, which has been introduced essentially to describe cluster degrees of freedom, can also address continuum states.

The other interesting point is that the exact Fermi gas result is obtained with any arbitrary number of particles in the simulation. For the 3-particles system shown in Figure 1, the triangles in the lower right part of the figure show the momentum values allowed for this finite system. These values are eigenvalues of the periodic part of the wave function Φ, and their spacing decreases with increasing number of particles. As shown by eq.( 50), the effect of the θ quantum number is to produce all the other missing values for the k quantum number which would be accessible to the infinite system, thus reproducing the continuous Fermi distribution.

In conclusion, the application of TABC to Bloch single particle wave functions appears as a very powerful method to recover the correct kinetic energies and wave functions at the thermodynamic limit with variational methods applied to small systems even in one dimension. A word of caution is however necessary: to demonstrate the Bloch theorem we have explicitly used the fact that the states are eigenvectors of the many-body Hamiltonian, which is in general not true when using variational methods. For this reason we will develop an alternative scheme in the next section.

III. THE REPLICA METHOD

The twist averaged boundary condition are currently applied in QMC calculations of correlated electron systems at the thermodynamic limit. In the specific framework of fermionic molecular dynamics, an alternative method to deal with the infinite range antisymmetrisation correlations has been proposed, based on the localized Wannier representation [START_REF] Vantournhout | Nuclear pasta with a touch of quantum : Towards fully antisymmetrised dynamics for bulk fermion systems[END_REF]. Taking advantage of the nested structure of the overlap matrix in a periodic system, an analytical method of inversion of this matrix is proposed in ref. [START_REF] Vantournhout | Nuclear pasta with a touch of quantum : Towards fully antisymmetrised dynamics for bulk fermion systems[END_REF]. We present here a slightly different derivation of the same equations, for an application to any arbitrary single particle basis.

A. General formulation

Let us consider that our WS cell contains already many different replicas of the physical system, N tot = N r N where finally we will let N r → ∞. Similarly, we now interpret the simulation cell L tot = nL as a huge length, n ≫ 1, such that working with a finite L tot can be considered as equivalent to the thermodynamic limit, meaning that the scaling has to be fulfilled for all observables Ô

< Ô >= N r < Ô > L (56)
with

< Ô > L = 1 N r Ntot j,k=1 o jk B -1 kj (57)
independent of the periodicity properties of the Hamiltonian. If we consider the N tot -body wave function, the invariance property of the system under the translation mL tot of any particle coordinate is still verified. The difference respect to the previous section is that the system composed of N tot particles has an additional symmetry which is not directly linked with the use of a finite box L but is rather due to the physical periodicity of the pasta structures over a lenght scale L = L tot /n. We have already observed that, because of the absence of an external field, the cell Hamiltonian eq.( 2), which now represents the full Hamiltonian eq.( 1),

Ĥ = Ntot i=1 ti + ∞ m=-∞ Ntot i>j v (x i -xj -mL tot ) ( 58
)
is invariant under the simultaneous translation of every particle coordinate of any arbitrary length. The translational invariance which will be relevant to us is the invariance under the simultaneous translation of every particle coordinate of the physical length L, defined as the periodic length of the one body density

ρ(x + mL) = ρ(x) (59)
If the Hartree-Fock field as an external field, this translational invariance would be the only one verified by the meanfield hamiltonian, and it would be verified for each particle coordinate separately. The situation is different in the case of an interacting system as described by the self-consistent Hartree-Fock approach. In this case, translational invariance with respect to any arbirary length is respected because of the self-consistency of the mean-field approach, U ij = δE HF /δρ ji where E HF is the interaction part of the Hartree-Fock energy. However, this translational invariance applies only to the simultaneous translation of all particle coordinates. The only special feature of the length L is that eq.( 59) is fulfilled, which will allow us to impose a simplified expression for the one-body wave functions, as we will see in a moment. Following the derivation of the previous section, these two invariance properties lead to the definition of a Bloch phase for the N tot -body wave function according to

Ψ (x 1 , . . . , x Ntot ) = exp i θ l LN tot Ntot k=1 x k Φ (x 1 , . . . , x Ntot ) ( 60
)
where Φ is a periodic function, that is invariant both under the translation mL of all N tot particle coordinates simultaneously, and under the translation mL tot of any particle coordinate.

The difference respect to the derivation of the previous section is that now we have a factor 1/N tot multiplying the Bloch phase. Since we are at the thermodynamic limit N tot → ∞, the correction on kinetic observables induced by this Bloch phase can be neglected and we can consider that the total N tot -body wave function Ψ shares the periodicity properties of Φ.

Let us now look for an ansatz for this L-and L tot -periodic Ψ (x 1 , . . . , x Ntot ) to be introduced as a variational approximation to the full N tot -body problem. We take the same Slater ansatz used in the previous section:

Ψ (x 1 , . . . , x Ntot ) = Â Ntot k=1 u k (x) (61) with u k (x + mL tot ) = u k (x).
Because of the periodicity of the one body density, each sub-cell of linear dimension L can be exactly associated to a finite number N = N tot /N r particles, and the wave function can be written introducing a sub-cell index m as:

Ψ (x 1 , . . . , x Ntot ) = Â N k=1 Nr m=1 u k,m (x) (62)
Moreover, the periodicity of ρ can be written as

ρ(x + L) = Nr m,m ′ =1 N k,j=1 u * k,m (x + L)u j,m ′ (x + L)B -1 jk = Nr m,m ′ =1 N k,j=1 u * k,m (x)u j,m ′ (x)B -1 jk = ρ(x) (63)
This condition can be satisfied introducing only N independent one body wave functions u k,m (x) = u k (x -mL) where in principle u k are arbitrary (non periodic) functions. If these functions would form an orthonormal basis

< j, m|k, m ′ >= δ jk δ mm ′ (64)
the problem would be reduced to a N body problem as with the ansatz eq.( 21). Indeed in the calculation (57) of an arbitrary one-body observable only the u k , k = 1, . . . , N functions would be needed

< Ô > L = 1 N r Ntot k=1 o kk = N k=1 o kk . (65)
B. Application to non-orthogonal single-particle basis

For practical applications we will not work with an orthonormal basis, but we want to take the FMD/AMD gaussian wave packet variational ansatz eq.( 23), namely

u k (x -ml) = g Z k +ml (x) (66)
Since the g Z k +ml do not form an orthonormal basis, we are left with a problem of computational size N 2 tot , that is impossible to solve (recall N tot → ∞).

A possibility would be again to impose the simplified expression:

Ψ (x 1 , . . . , x Ntot ) = Nr m=1 Â N k=1 g Z k +ml (x) = Nr m=1 Ψ (x 1 , . . . , x N ) (67)
which allows to work with only N functions.

The conceptual problem with this simple ansatz eq.( 67) is that it neglects the antisymmetrization correlations among different cells. We know that such correlations are important as they are at the origin of the band structure. For this reason we need to keep the antisymmetrization correlation among all the N tot = N r • N particles.

By introducing a linear transformation

u k = Û g k (68)
defined by

u k,n (x) = 1 √ N r Nr/2 m=-Nr/2 exp im 2πn N r g Z k +mL (x), (69)
we can reduce the computation size to N r N 2 keeping the antisymmetrization correlations among all the N tot = N N r particles. It is interesting to remark that, in the limit of N r → ∞, this transformation coincides with the definition of Wannier functions [START_REF] Ashcroft | Solid State Physics[END_REF]. This transformation preserves the determinant

 Nr m=1 g Z k +mL (x) =  Nr n=1 u k,n (x) (70)
so that we can write

Ψ (x 1 , . . . , x Ntot ) = Â Nr n=1 N k=1 u k,n (x) (71)
The advantage of the transformation (69) is that matrices are block-diagonal, that is functions corresponding to different phases γ = (2n -N r )π/N r are orthogonal

< u j,n |u k,n ′ >= δ n,n ′ B jk (n). (72)
We have thus shown that the decoherence hypothesis eq.( 19) of the TABC method is exactly verified in the case of variational approaches based on Slater determinants.

Having recovered the orthogonality on the level of the replica quantum number we can write

< Ô > L Ψ = 1 N r Nr n=1 N j,k=1 < u j,n |ô|u k,n > B -1 kj (n) (73)
The matrix element of any operator is given by

o jk (n) =< u j,n |ô|u k,n >= L/2 -L/2 dx ∞ m,m ′ =-∞ g * Zj +mL (x ′ )o(x, x ′ )g Z k +m ′ L (x) exp -i(m -m ′ ) 2πn LN r (74)
and explicitly depends on the n quantum number associated to the antisymmetrization among the different subcells.

In particular the overlap matrix reads

B jk (n) =< u j,n |u k,n >= L/2 -L/2 dx ∞ m,m ′ =-∞ g * Zj +mL (x)g Z k +m ′ L (x) exp -i(m -m ′ ) 2πn LN r (75)
Using the same argument that we had for the expectation (36), the series can be reduced to a finite sum

B jk (n) = M m,m ′ =-M L/2 -L/2 dxg * Zj +mL (x)g Z k +m ′ L (x) exp -i(m -m ′ ) 2πn LN r (76)
In practical application, a large value of M ≈ 10 is needed to describe non-interacting particles, but M = 1 turns out to be sufficient in the presence of the nuclear and Coulomb interaction for the particle densities relevant for the neutron star crust. From the computational viewpoint, M = 1 is equivalent to calculate a single gaussian overlap between each pair if periodic boundary conditions are applied to the gaussian, meaning that the simulation of the infinite system is not more expensive numerically than the one for the finite system. The analogous expression for the matrix elements reads

< u j,n |â|u k,n >= M m,m ′ =-M L/2 -L/2 dxg * Zj +mL (x)a(x)g Z k +m ′ L (x) exp -i(m -m ′ ) 2πn LN r (77)
Going to the N r → ∞ limit we can write lim

Nr→∞ 2πn LN r = γ (78)
with -π < γ ≤ π. The observables read

< Ô > L Ψ = 1 2π π -π dγ N j,k=1 < u j,γ |ô|u k,γ > B -1 kj (γ) (79)
which is expected to give equivalent results to the implementation where TABC are applied. The FMD ground state energy for a system of two free particles in one dimension is shown in Figure 2 as a function of the number of replicas, with the Wannier choice for the wave packets. For each calculation, the number M was chosen high enough to have convergent results. We have already observed that the FMD variational space is sufficiently flexible to describe a free-particles system. The FMD solution for each value of N r thus represents the exact energy of a finite system with periodicity L = 2/ρ and size L tot = L • N r . The thermodynamic limit, corresponding to the anlytical Fermi gas result, is obtained for about ten replicas, corresponding to ten values for the phase γ of eq.(78).

IV. COMPARISON OF THE TWO METHODS

In the previous sections, we have introduced the Bloch and the Wannier representation as two alternative ways to address the thermodynamical limit in variational theories. From a principle point of view the two approaches are very similar, however they lead to very different expressions in the calculation of observables. In particular, the fact of considering only diagonal terms in the phase quantum number naturally emerges as a consequence of the orthogonality of the Wannier basis in the replica method, while it appears as a simple working hypothesis in TABC. Moreover, the equivalence is expected in the ideal case of an infinite number of twists and under the condition that the wave function is an exact eigenvalue of the Hamiltonian. This first condition is numerically expensive, while the second is not assured in variational methods. For all these reasons, the equivalence of the two formalisms is not a-priori clear and will be discussed in the present section.

A. Finite number of replicas and finite number of twists

We now show that, in the hypothesis that the variational theory produces the exact eigenvectors of the many-body Hamiltonian, the two methods to treat the antisymmetrisation correlations can be mapped on each other for any system size.

In the Bloch formalism, if we impose that the condition of periodic boundaries eq.( 14) is verified for a finite translation of size N r L, the original problem of an infinite system with periodicity L is transformed into the problem of a finite system of size N r L, with periodic boundary conditions.

The allowed values for θ are then discrete θ = 2n-Nr Nr π, where n is an integer. This amounts to discretize the integral (20) with ∆θ = 2π

Nr . We expect then the two methods to give exactly the same results. This expectation is confirmed by Figure 3, which compares for the case of the non-interacting one-dimensional system the energy obtained considering a periodic system of N tot = ρN r L particles, with the calculation of a reduced N = ρL system averaged over N r phases with the two (Bloch and replica) proposed methods. We can see that the agreement is perfect for all sizes, and all calculations converge to the free Fermi gas result with increasing number of replicas (respectively: twists). FIG. 3: Energy of a free periodic one-dimensional system at a density ρ = 1.37 π as a function of the number of replicas. The upper abscissa represents the total linear size in the case of simple periodic boundary conditions (full line), the lower one is a number of replicas for the replica method (crosses), and a number of twists for the Bloch method (circles). Dash-dotted line: free Fermi gas in one dimension.

This shows the equivalence of the different methods when the wave vector is an eigenvalue of the Hamiltonian.

B. Limitations of the Bloch method in variational applications

In realistic applications of interacting fermion system, the variational solution is always an approximation of the exact eigenvalue. This may create a bias in the different method that we turn to discuss. For eq.( 4) to correctly account for the correlations due to the particles outside the cell L, it is necessary that the wave function Ψ θ is an eigenfunction of the hamiltonian. Conversely the Wannier representation can be obtained from the solution of the infinite periodic system via a simple linear transformation eq.( 68). Of course if the variational ansatz is not adequate to describe the exact state both representations will be false. However the error due to the inadequacy of the solution in the cell L will be propagated in an uncontrolled way in the Bloch method for the description of the global system, and this will not be the case for the replica method that directly addresses the replicated system of size N r L.

To illustrate this difference, we can once again take the simple example of the non-interacting Fermi gas. There is a popular microscopic theory in nuclear physics, the so-called AMD model [27], where the single particle wave packets are gaussians of fixed width. With this ansatz the variational space does not contain the exact solution of this problem at the bulk limit, that is an antisymmetrized combination of plane waves. In particular, if we choose a small value for the width, as it is variationally obtained to describe finite nuclei [27], it is clear that the states will be very far from the exact eigenvectors of the kinetic hamiltonian. The comparison of the two methods in this model case is shown in Figure 4. FIG. 4: Comparison of the different boundary conditions for the AMD model in the case of a non interacting system as a function of the system size. Dash-dotted line: free Fermi gas. Dashed line: AMD ground state energy for a system size large enough for finite size effects be negligible (L = 150).Circles: finite system with Bloch boundary conditions. Cross: replica method. Stars: periodic boundary conditions.

We can see in Figure4 that the AMD ansatz is not adapted for this problem, from the fact that asymptotically the energy of the free Fermi gas (full line) is not recovered by the AMD calculation (dashed line). This asymptotic value is perfectly reproduced by the replica method (crosses) for any system size, while the Bloch method (open circles) produces artificial fluctuations and attains convergence only when the Bloch phase is negligibly small, and the calculation is perfectly equivalent to simple periodic boundary conditions (stars).

The non-interacting system is certainly not the ideal application ground of molecular dynamics model. As a second model case which is closer to the physical case of the neutron star crust we take a one-body periodic potential, which is obtained as a periodic self-consistent mean field in realistic applications of mean-field variational models [START_REF] Sandulescu | [END_REF][8][9][10]12] if the two-body nuclear and Coulomb interaction is added, including the interaction with a uniform electron background. For this model one-dimensional application we only wish to discuss the effect of the boundary conditions, therefore we restrict to the simpler case of an external periodic potential. In particular, the choice of a gaussian potential has the advantage of making all calculations analytical. Let us take a potential form:

V gaus (x) = V 0 N 2 n=-N 2 ∞ m=-∞ exp - 1 2a (x -b n -mL) 2 (80)
where V 0 , a, b n are parameters. In particular, the choice a = 4 ρ , b n = n ρ for a system with average density ρ, corresponds to having a single particle per potential well. The average energy is readily calculated as

Vgaus Ψ = V 0 lim Nr→∞ 1 N r N i,j=1 B -1 ij 2π a i a * j a a * j a + a i a + a i a * j ∞ m1,m2,n=-∞ exp - a * j (Z i -b n + m 1 L) 2 + a i (Z * j -b n + m 2 L) 2 + a(Z i -Z * j -(m 1 -m 2 )L) 2 a * j a + a i a + a i a * j (81)
Results from the Bloch method for a case of 3 particles with an average density ρ = 1.37 π are presented in Figure 5. At variance with the free particle case, the periodic part of the wave function Φ now explicitly depend on the θ quantum number, leading to a θ dependence of the one body density.

The comparison between the different methods of treating the boundary conditions is shown in figure 6. This model is less trivial than the free Fermi gas, although still very schematic. The FMD model is expected to be a good approximation of this system, but there is no guarantee the FMD solution should be exact. Similar to the previous application, the replica method gives by construction the asymptotic result for any finite size, provided a sufficient number of phases is considered. Conversely for the Bloch method, which is calculated with the same number of phases as the replica method for the application of Fig. 6, a convergence is reached only when the effect of the phase can be neglected.

This model example shows that in realistic applications, where any variational ansatz might be simply an approximation of the exact energetics of the system, the different boundary conditions are not equivalent, and the employ of Bloch corrections to the single particle basis may increase the deviation respect to the exact solution.

V. CONCLUSIONS

In this paper we have introduced and critically discussed two different ways of addressing the boundary conditions in finite fermionic system in order to account for the infinite range antisymmetrisation correlations present at the thermodynamical limit. In the case of state vectors which are eigenstates of the many-body Hamiltonian, we have shown that the use of Bloch wave functions with twist averaged boundary conditions physically corresponds exactly to a bigger system constituted of a number of replicas of the system under study equal to the number of phases. In turn, this last formalism can be interpreted as the use of localized Wannier states. While these different representations are equivalent when dealing with eigenvectors, the same is not true in the case of approximate solutions of the many-body problem, as it is the case in the variational approach. In this case, the approach of the thermodynamic limit is not properly treated by the Bloch method, while the replica method provides a very accurate evaluation of the Fock energy. Even in the case of a non-orthogonal single particle basis this method can be numerically implemented with a moderate computational cost, opening the possibility of a completely quantum-mechanical treatment of nuclear matter with molecular dynamics approaches. FIG. 6: Comparison of the different boundary conditions for the FMD model in the case of a system subject to an external gaussian potential as a function of the system size. Dash-dotted line: free Fermi gas. Dashed line: FMD ground state energy for a system size large enough for finite size effects be negligible (L = 30). Circles: finite system with Bloch boundary conditions. Cross: replica method. Stars: periodic boundary conditions.

16 -FIG. 1 :

 161 FIG.1: FMD calculation of a system of three free particles in a one dimensional box of size L with Bloch boundary conditions. Upper part: behavior of the total linear momentum (left) and total energy (right) as a function of the Bloch wave number k0 = θ/L Full lines: total expectation values K Ψ, Ĥ Ψ (see text); dashed lines: averages over the variational part of the wave function K Φ, Ĥ Φ without the contribution of the Bloch phase (see text); dash-dotted line: total energy average over the Bloch phases. Lower part: spatial (left) and momentum (right) density.

FIG. 2 :

 2 FIG.2: Energy associated to a system of two free FMD particles in one dimension at a density ρ = 1.37 π as a function of the number of replicas. Dash-dotted line: free Fermi gas in one dimension.

FIG. 5 :

 5 FIG.5: FMD calculation of a system of three particles in a one dimensional box of size L subject to a periodic external potential with Bloch boundary conditions. Upper part: behavior of the total linear momentum (left) and total energy (right) as a function of the Bloch wave number k0 = θ/L Full lines: total expectation values K Ψ, Ĥ Ψ; dash lines: averages over the variational part of the wave function K Φ, Ĥ Φ without the contribution of the Bloch phase. Dashed-dotted lines: expectation values including the phase average. Lower part: dash-dotted lines: spatial (left) and momentum (right) density; full line: gaussian potential.

Acknowledgments

This paper has been partly supported by ANR under the project NEXEN. Discussions with Klaas Vanternhout and Karim Hasnaoui are gratefully acknowledged.


Adeline Augier
email: adeline.augier@math.u-psud.fr

François Dubois

Benjamin Graille
email: benjamin.graille@math.u-psud.fr

Francois Dubois
email: francois.dubois@math.u-psud.fr

Isotropy conditions for lattice Boltzmann schemes. Application to D2Q9

may come

Introduction

Numerous phenomena, in particular acoustic properties, are modelled by partial differential equations (PDEs) like Navier-Stokes equations. One of the main property of these PDEs is to take into account the isotropy of the environment. Namely, if the physical properties of the environment do not depend on the orientation, the associated PDEs are invariant by orthogonal transformations of the spatial frame. In this paper, we are interested by lattice Boltzmann schemes [START_REF] Humières | Generalized Lattice-Boltzmann equations[END_REF][START_REF] Dubois | Une introduction au schéma de Boltzmann sur réseau[END_REF][START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF] used by the industrialists because of their efficiency. As the directions of the lattice are privileged, they are not a priori numerically isotropic. Because of the criss-cross pattern of the scheme DdQq (where d represents the space dimension and q the number of discrete velocities), it is natural to look for parameters that give the same behaviour on each axes (i.e by a rotation of angle π/2 in 2D for example) and by extent we want to obtain this result for each orthogonal transformation of the spatial frame. Our purpose is then to understand and control the lack of isotropy for these schemes.

The first thing to do consists in defining the isotropy for lattice Boltzmann schemes. Actually, there is no uniform way to define isotropy: for example in [START_REF] Chen | Lattice Boltzmann computational fluid dynamics in three dimension[END_REF] the authors focus on the maintain of isotropy of the stress tensor required by fluid equations. An other point of view is described in [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF]: an isotropic fluid is simulated by eliminating as much as possible the effects due to the symmetry of the lattice and it means solving an eigenvalue problem. The problem of the space transformation is evoked in [START_REF] Van Der Sman | Diffusion lattice Boltzmann scheme on a orthorhonbic lattice[END_REF][START_REF] Yu | DNS and LES of decaying isotropic turbulence with and without frame rotation using the lattice Boltzmann method[END_REF]: in [START_REF] Yu | DNS and LES of decaying isotropic turbulence with and without frame rotation using the lattice Boltzmann method[END_REF] rotating frames are considered and in [START_REF] Van Der Sman | Diffusion lattice Boltzmann scheme on a orthorhonbic lattice[END_REF] the authors introduce a definition of isotropic diffusion that constraints the collision operator to be invariant under all isometries of the Bravais lattice. In this paper, we define isotropy as the invariance by the orthogonal transformations of the frame and we express it on the equivalent equations. Namely, the equivalent equations are PDEs that come from a Taylor expansion of the lattice Boltzmann scheme [START_REF] Dubois | Une introduction au schéma de Boltzmann sur réseau[END_REF] and a PDE is said isotropic if it is invariant by rewriting it in the "new" frame. Finally, a DdQq frame is isotropic at order M if its set of equivalent equations of order n, 1 ≤ n ≤ M , is isotropic.

In the first Section we quickly recall the lattice Boltzmann method and give some results on the equivalent equations.

In the second Section, we exploit the definition of isotropy for PDEs in an algebraic point of view. More precisely, we first define a new frame by orthogonal transformation. After redefining all the objects in the new frame, we identify both of the equivalent equations: the old one and the new one. The isotropy is guaranteed as soon as the various tensors in the PDEs appear as a fixed point of some group operation. Then, we obtain a very useful set of relations on the coefficients of the equivalent equations.

In the last Section, a lattice Boltzmann scheme is said isotropic at order M if its set of equivalent equations of order n, 1 ≤ n ≤ M , is isotropic in the sense explained in Section 2. Since the coefficients of the equivalent equations depend on the parameters of the scheme (relaxation times and equilibrium states), we are able to characterize, for the related scheme, the lack of isotropy on these parameters for each order. More precisely, we are able to determine order after order the relations between the parameters of the scheme that have to be satisfied to have isotropic equivalent equations. We determine these relations for the D2Q9 scheme and specify all the sets of parameters that give isotropy until the fourth order.

The lattice Boltzmann schemes considered in this paper are linear (for applications reasons) while the proposed definition of isotropy is more general and can be used as soon as the equivalent equations of the scheme are known.

Lattice Boltzmann method and equivalent equations

In this section, we recall some notations related to the lattice Boltzmann method. Then, we introduce the equivalent equations as well as their link with the physic. This link is essential because it allows us to decide what parameters have physical values and what parameters have to be chosen in a allowed set of parameters. We will see below (in section 3) what for a choice we have for this set of parameters in order to guaranty isotropy of the scheme.

Notations and description of the lattice Boltzmann method

We use the scheme proposed by d'Humières in [START_REF] Humières | Generalized Lattice-Boltzmann equations[END_REF]: we consider a regular lattice L with typical mesh size ∆x. The time step ∆t is determined thanks to the velocity scale λ: ∆t = ∆x/λ. For the DdQq scheme, we note v = (v j ) 1≤j≤q the set of q velocities and we assume that for each node x of L, and each v j in V , the point x + v j ∆t is also a node of the lattice L. The aim of the DdQq scheme is to precise the particle distribution f = (f j (t, x)) t 1≤j≤q (where t represents the transpose) for x ∈ L and discrete values of time t by solving a discretization of the Boltzmann's equation in two steps for every time: collision and transport. We finally have to propagate the so-called lattice Boltzmann scheme:

f j (x + v j ∆t, t + ∆t) = f * j (x, t), ∀ x ∈ L , 1 ≤ j ≤ q
, where f * j represents the jth component of the particle's distribution after the collision. In order to implement this scheme, as it is done in [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF], we introduce the orthogonalized moments m defined by m := M f where M is a given and invertible matrix. Concerning the moments m, we distinguish those who are preserved during the collision, called W ∈ R N , of the complementary ones denoted by Y ∈ R q-N . Since we are interested in linear applications, we assume that the collision operator is linear, so there exists a matrix J such that m * := J m, with the following structure:

m * = W * Y * = J m = Id N 0 S E Id q-N -S W Y
,

where E ∈ M q-N,N (R) and S is the diagonal matrix of the relaxation times s k = ∆t/τ k , N + 1 ≤ k ≤ q, with 0 < s k < 2, N + 1 ≤ k ≤ q for stability conditions. We introduce a last classical notation:

σ k = 1/s k -1/2 ∈ R + , 1 ≤ k ≤ q.
Finally a time step of linear lattice Boltzmann scheme scheme reduces to:

f (x + v j ∆t, t + ∆t) = M -1 J M f (x, t), ∀ x ∈ L, 1 ≤ j ≤ q.
(1)

Thus, the parameters of the lattice Boltzmann schemes are the coefficients of the matrices S and E. This coefficients are usually chosen thanks to the experiment (by the physicists and the industrialists). The aim of this paper is to precise these coefficients according to the property of isotropy.

Remark 1.1. There exist three types of lattice Boltzmann schemes according to the values of the relaxation times: either all of relaxation times are equal (BGK for Bhatnagar -Gross -Krook scheme) or there is only two different relaxation times (TRT for Two Relaxation Times) or there is no conditions on these (MRT scheme for Multiple Relaxation Times). Since we investigate the behavior of the scheme in the sense of the isotropy, it is natural to consider the third case: it is the most general case even if it is not the most often used. Thus, there are q -N independent unknowns parameters on S.

Focus now on the matrix E that is chosen full in order to give us as degree of freedom as possible: it gives N (q -N ) additional parameters. Actually, it is not an usual choice: the matrix E is usually chosen sparse by physicists, but in our knowledge, there is no mathematical justification to this type of assumption (based on the physical "good sense"). It is remarkable, that the isotropy condition for the DdQq scheme can give some justification of this set of coefficients.

Finally, both matrices E and S give (N + 1) × (q -N ) parameters to be precise in order to conserve physicals and isotropic properties.

For example, for the D2Q9 scheme, we consider the classical moments: the density ρ and the momentum = (q x , q y ) t that are conserved during the collision and the energy ε, the square of the energy ε 2 , the flux (ϕ x , ϕ y ) t , and the tensor (ϕ xx , ϕ xy ) t . Then, we orthogonalize them like it is done in [START_REF] Dubois | Quartic parameters for acoustic applications of lattice Boltzmann scheme[END_REF]. For this choice of moments we define both equilibrium states and relaxation times matrices by:

E t :=   E ρ ε λ 2 E ρ ε2 λ 4 E ρ ϕx λ 3 E ρ ϕy λ 3 E ρ ϕxx λ 2 E ρ ϕxy λ 2 E qx ε λ E qx ε2 λ 3 E qx ϕx λ 2 E qx ϕy λ 2 E qx ϕxx λ E qx ϕxy λ E

Equivalent equations

We first introduce the notation • i j for the element of • that is on the intersection of the ith line and the jth column. By convention, a Latin letter is an index related on the moments (of size N ) while Greek letter represents an index related on the space dimension (of size d). We then introduce the system of N equivalent equations [START_REF] Dubois | Une introduction au schéma de Boltzmann sur réseau[END_REF][START_REF] Dubois | Third order equivalent equation of lattice Boltzmann scheme[END_REF] of order M with which the linear lattice Boltzmann scheme is consistent at order M :

∂ t W i + M n=1 n A i : ∇ n W = 0, 1 ≤ i ≤ N, (2)
where n A i , 1 ≤ n ≤ M, 1 ≤ i ≤ N are tensors of order n that take into account the coefficients of E and S. By convention, the maximal contraction operator : is defined by:

( n A : ∇ n W ) i := n A i : ∇ n W := 1≤j≤N 1≤α1,••• ,αn≤d n A i,α1,••• ,αn j ∂ α1 • • • ∂ αn W j , 1 ≤ i ≤ N, 1 ≤ n ≤ M.
These equivalent equations come from a formal calculus explained in [START_REF] Dubois | Quartic parameters for acoustic applications of lattice Boltzmann scheme[END_REF] and an algorithm that is easy to use is described in [2].

Remark 1.2. The tensors n A belongs to the space R N 2 n d , 1 ≤ n ≤ M , but because of the Schwarz property, numerous of coefficients of n A i are equal (or are proportional between them). In order to count how much of these coefficients are independent, we have to calculate the number of d-uplet The equivalent equations ( 2) are also very useful to make the link with the physics. For example, it is wellknown (see for example [START_REF] Dubois | Third order equivalent equation of lattice Boltzmann scheme[END_REF]) that there exists a D2Q9 scheme (i.e. a choice of moments Y ∈ R 6 ) such that ( 2) is consistent with Navier-Stokes equations at order 2. Then, there exists a relation between both viscosities (shear viscosity ζ and bulk viscosity µ) and two relaxation times (see for example [START_REF] Dubois | Quartic parameters for acoustic applications of lattice Boltzmann scheme[END_REF]). Furthermore, equilibrium state of the energy given by E ρ ε ρ + E qx ε q x + E qy ε q y + E qz ε q z is such that E ρ ε only depend on the sound velocity. Finally, our purpose is to precise as much as possible the (N + 1)(q -N ) -3 parameters of the lattice Boltzmann scheme in order to understand the lack of isotropy of these schemes: for example, in the D2Q9 scheme, parameters that are given thanks to the environment are E ρ ε (it depends on the sound velocity), σ ε (it depends on the bulk viscosity ζ) and σ ϕxx (it depends on the shear viscosity µ) that means 21 parameters.

(β 1 , • • • , β d ) ∈ N d such that

Algebra for rotating the equivalent equations

We first recall the definition of an isotropic PDE:

Definition 2.1. A PDE ∂ t W i + M n=1 n A i : ∇ n W = 0, 1 ≤ i ≤ N, (3)
is said isotropic if it is invariant by orthogonal transformation of the frame.

We introduce the matrix r ∈ O d (R) that represents an orthogonal transformation of the frame (O, e). We then use the notation • for every element in the "new" frame: in this way, if the vectors x ∈ R d and the moments W ∈ R N are known in the frame (O, e), we have to define the news moments W , the new vector x (such that x := r x), the new basis e, the news tensors The aim of this section consists in giving the relation between objects that are defined before and after the transformation of the frame. First, we give a relation between both partial derivatives and we deduce a relation on both derivatives of the moments. Second, we introduce a group operation via an homomorphism that will show itself very useful for the definition of the isotropy. Lemma 2.2. For every 1 ≤ α ≤ d and 1 ≤ i ≤ N , the derivatives ∂ α W i at the ith component of the moment W reads:

n A i , 1 ≤ n ≤ M , 1 ≤ i ≤ N
∂ α W i = ∂ α (R(r) W ) i = j,β R i j (r)(r -1 ) β α ∂ β W j . (4)
Proof. As x = r x, for each 1 ≤ α ≤ d, we have the following relation on the partial derivatives ∂ α along the αth direction:

∂ α = 1≤β≤d ∂ x β ∂x α • ∂ β = 1≤β≤d r -1 β α ∂ β .
Then, it is easy to characterize the partial derivatives of the moments and we obtain the desired result (4).

We also define a group operation through the following Proposition.

Definition -Proposition 2.3. Let n be in N * and the map Φ

n : O d (R) × S(n, d) → S(n, d) defined for 1 ≤ i, j ≤ N , 1 ≤ α 1 , • • • , α n ≤ d, 1 ≤ n ≤ M
, by the relation:

(Φ n (r)( n A)) i,α1,••• ,αn j := 1≤β1,••• ,βn≤d 1≤k,l≤N (R(r)) i l n A l,β1,••• ,βn k r α1 β1 • • • r αn βn (R(r) -1 ) k j . (5)
It is an homomorphism that characterizes a group operation.

Proof. Let Id be the identity matrix of O d (R), then for every n A ∈ S(n, d) we have Φ n (Id)( n A) = n A. It remains to prove the following property:

∀ r, r ′ ∈ O d (R), ∀ n A ∈ S(n, d), Φ n (rr ′ )( n A) = Φ n (r) (Φ n (r ′ )( n A))
. First, by using the definition of the homomorphism Φ n given in (5), we get for 1 ≤ i, j ≤ N and 1

≤ α 1 , • • • , α n ≤ d: Φ n (r) (Φ n (r ′ )( n A)) i,α1,••• ,αn j = 1≤β1,••• ,βn,γ1,••• ,γn≤d 1≤k,l,m,o≤N R(r) i l R(r ′ ) l m n A m,γ1,••• ,γn o r α1 β1 (r ′ ) β1 γ1 • • • r αn βn (r ′ ) βn γn (R(r ′ ) -1 ) o k (R(r) -1 ) k j . (6)
Since R(rr ′ ) = R(r)R(r ′ ), and so R(rr

′ ) -1 = R(r ′ ) -1 R(r) -1 , (6) gives: Φ n (r) (Φ n (r ′ )( n A)) = Φ n (rr ′ )( n A)
and the homomorphism φ n truly defines a group operation.

Property 2.4. Let r be an orthogonal transformation of O d (R), then the system (3) reads in the new frame

∂ t W i + M n=1 n A i : ∇ n W = 0, 1 ≤ i ≤ N, (7)
where n A = Φ n (r)( n A), 1 ≤ n ≤ M , or equivalently n A = Φ n (r -1 )( n A), 1 ≤ n ≤ M .

Proof. Let us write (3) after the orthogonal change of frame r. Because of the equality W = R(r) W , ∂ t W becomes ∂ t R(r) W . Using the property (4), we obtain a relation that gives n A i : ∇ n W in function of the partial derivatives of the moments in the new frame, more precisely for 1 ≤ i ≤ N we get:

n A i : ∇ n W = 1≤β1,••• ,βn,γ1,••• ,γn≤d 1≤k,l,m,o≤N n A i,α1,••• ,αn j r -1 β1 α1 • • • r -1 βn αn R(r) j k β 1 • • • β n W k .
Thus, by identifying both (3) and R(r) × (7), we characterize n A, 1 ≤ n ≤ M , thanks

n A i : ∇ n W = 1≤j≤N R(r) -1 i j n A j : ∇ n W , ∀ W ∈ R N , 1 ≤ n ≤ M, 1 ≤ i ≤ N,
and that gives the desired result.

Finally, we obtain a characterization for a PDE to be isotropic:



  and S := Diag(s ε , s ε2 , s ϕx , s ϕy , s ϕxx , s ϕxy ).

d

  k=1 β k = n. This is a well-known combinatorics result: this number is the number of d combinations with repetition of a set which cardinal number is n and it reads: Γ n d := (d + n -1)!/(n!(d -1)!) and n A ∈ S(n, d), where S(n, d) is the space of the symmetrical tensors of R N 2 n d .

  and the new partial derivatives ∂ α , 1 ≤ α ≤ d. Finally, let be R(r) the matrix in O N (R) such that W := R(r) W . More precisely, R : O d (R) → O N (R) defines a morphism that gives a group representation. In the present contribution, mass and momentum are conserved thus W = (ρ, q) t , N = d + 1 and for every r ∈ O d (R), we get R(r) = Diag(1, r) ∈ O N (R).

Corollary 2 . 5 .

 25 Let (3) be a PDE of order M , it is isotropic if and only if n A = n A, 1 ≤ n ≤ M , ∀ r ∈ O d (R) and it means that the tensor n A is a fixed point of Φ n , 1 ≤ n ≤ M , that is Φ n (r)( n A) = n A, ∀ r ∈ O d (R), 1 ≤ n ≤ M .

* This work has been financially supported by the French Ministry of Industry (DGCIS) and the Region Ile-de-France in the framework of the LaBS Project [1].

Isotropy conditions for lattice Boltzmann schemes

In this section, we first give a definition of isotropy for lattice Boltzmann schemes and then we give an application on the D2Q9 scheme. Definition 3.1. A lattice Boltzmann scheme is said isotropic at the order M if the system of equivalent equations (2) at the order M is isotropic.

Furthermore, we denote by L N (r) := 1≤n≤N (Φ n (r)( n A)n A)△t n the lack of isotropy at order N of the orthogonal transformation r of the equivalent equations (2).

Then, we use this definition on the D2Q9 scheme in order to obtain all the sets of parameters that improve this scheme (in the sense of the isotropy) until the fourth order. In this case, the orthogonal transformation r is the well-known rotation matrix, d is equal to 2, N = 3 and M = 4. Proposition 3.2. Let L 4 (r) be the lack of isotropy for the D2Q9 scheme at fourth order for the orthogonal transformation r, then we get:

Thanks to [START_REF] Dubois | Quartic parameters for acoustic applications of lattice Boltzmann scheme[END_REF] this last condition gives the following relation on both viscosities: ζ = 6µ. Remark 3.3. Some of the properties for an isotropic D2Q9 scheme that are obtained here are well-known and naturally used in this scheme. In fact, the parameters that are taken null in order to obtain isotropy at first and second orders are also taken null for isotropic reasons considering the kinetic solution at equilibrium in the continuous environment. However, the results on the third and fourth orders are more surprising even if some usual D2Q9 schemes can be seen as particular cases of those we propose.

The proof of Proposition 3.2 is obtained thanks to a rigorous implementation of isotropy conditions with a software of formal calculus and [START_REF] Augier | Isotropy conditions for lattice Boltzmann schemes[END_REF] give details of it.


Nicolas Wiest-Daesslé

Olivier Commowick

Aymeric Stamm

Patrick Pérez

Christian Barillot

Romuald Seizeur

Sylvain Prima

Comparison of 3 diffusion models to track the hand motor fibers within the corticospinal tract using functional, anatomical and diffusion MRI

In this paper, we propose to compare three diffusion models to track the portion of the corticospinal tract dedicated to the hand motor function (called hand motor fibers hereafter), using diffusion, functional and anatomical MRI. The clinical diffusion data have few gradient directions and low b-values. In this context, we show that a newly introduced model, called diffusion directions imaging (DDI) outperforms both the DTI and the ODF models. This new model allows to capture several diffusion directions within a voxel, with only a low number of parameters. Two important results are that i) the DDI model is the only one allowing consistent tracking from the mesencephalon to the most lateral part of the cortical motor hand area, and that ii) the DDI model is the only model able to show that the number of hand motor fibers in the left hemisphere is larger than in the contralateral hemisphere for righthanded subjects; the DDI model, as the other two models, fails to find such a difference for left-handed subjects. To the best of our knowledge, this is the first time such results are reported, at least on clinical data.

Introduction

Diffusion MRI (dMRI) [START_REF] Le Bihan | Looking into the functional architecture of the brain with diffusion MRI[END_REF] allows in vivo and non-invasive imaging of tissue structure. It is based on the facts that i) the diffusion of water molecules is constrained by the micro-structure of the tissues (such as, typically, the white matter fibers in the brain), and that ii) MRI can be made sensitive to this diffusion, using specific MR pulse sequences. Diffusion models can then be devised, and their parameters can be estimated for further study and analysis of tissue architecture. The simplest model is that of a Gaussian diffusion function, which amounts to characterise the diffusion with a tensor (i.e. a 3×3 symmetric definite positive matrix), giving its name to diffusion-tensor imaging (DTI) [START_REF] Basser | Diffusion-tensor MRI: theory, experimental design and data analysis -a technical review[END_REF].

Fiber tracking, or tractography, has been developed to "reconstruct" or "dissect" the fiber tracts in vivo, and then infer brain anatomy [START_REF] Mori | Fiber tracking: principles and strategies -a technical review[END_REF]. Many of the association (e.g. uncinate fasciculus, cingulum) and commissural (e.g. transverse fibers of the corpus callosum) fiber tracts have been successfully reconstructed using i) clinical dMRI sequences (with few gradient directions and low b-values), ii) the simple Gaussian diffusion model and iii) simple deterministic streamline fiber tracking methods [START_REF] Mori | MRI Atlas of Human White Matter[END_REF]. On the contrary, it has proved much more difficult to reconstruct projection fibers, and especially the motor fibers of the corticospinal tract (CST) using such standard protocols. On one side, the anatomy of these fibers between the spinal cord and the internal capsule has been well-studied using DTI [START_REF] Kim | Corticospinal tract location in internal capsule of human brain: diffusion tensor tractography and functional MRI study[END_REF][START_REF] Hong | Somatotopic location of corticospinal tract at pons in human brain: a diffusion tensor tractography study[END_REF]. On the other side, the study of these fibers between the internal capsule and the cortex, and in particular of those dedicated to a specific motor function, is much more challenging using DTI, mostly due to the numerous crossings/kissing/merging/diverging fibers in the corona radiata. A particularly difficult fiber bundle to track within the CST is the portion corresponding to the motor hand area, because it is located laterally on the motor cortex [START_REF] Yousry | Localization of the motor hand area to a knob on the precentral gyrus. A new landmark[END_REF] (as shown by the homunculus of Penfield & Rasmussen), compared to the leg or the trunk areas for instance. In the following, we call this fiber bundle the HMFs, as "hand motor fibers".

The HMFs are a crucial fiber bundle to investigate, in the context of normal anatomy, within the more general study of handedness, cerebral dominance, and brain asymmetry [START_REF] Toga | Mapping brain asymmetry[END_REF]. More generally, the development of diffusion models and tractography methods for the CST in general, and the HMFs in particular, which could be used in clinical routine, is key for a better understanding of pathologies of the CST such as, typically, amyotrophic lateral sclerosis [START_REF] Iwata | Evaluation of corticospinal tracts in ALS with diffusion tensor MRI and brainstem stimulation[END_REF], Wallerian degeneration of the CST after ischemic stroke [START_REF] Yu | A longitudinal diffusion tensor imaging study on Wallerian degeneration of corticospinal tract after motor pathway stroke[END_REF], motor dysfunctions in infants [START_REF] Ludeman | Diffusion tensor imaging of the pyramidal tracts in infants with motor dysfunction[END_REF] or in patients with relapsing-remitting multiple sclerosis [START_REF] Lin | Diffusion tensor tractography-based group mapping of the pyramidal tract in relapsing-remitting multiple sclerosis patients[END_REF].

New types of image acquisition schemes (e.g. HARDI sequences), diffusion models (e.g. multiple tensors, ODF, etc.) [START_REF] Lenglet | Mathematical methods for diffusion MRI processing[END_REF] and tractography methods [START_REF] Mori | Fiber tracking: principles and strategies -a technical review[END_REF][START_REF] Mukherjee | Diffusion tensor MR imaging and fiber tractography: theoretic underpinnings[END_REF][START_REF] Mukherjee | Diffusion tensor MR imaging and fiber tractography: technical considerations[END_REF] have been introduced to account for intricate fiber configurations, but these techniques i) have been reported to often miss entirely the lateral portions of the CST, and thus the HMFs [START_REF] Behrens | Probabilistic diffusion tractography with multiple fibre orientations: What can we gain?[END_REF] and ii) are not applicable at hand in a clinical setting, mostly due to long acquisition times.

In this paper, we propose to investigate the usefulness of a new diffusion model, termed Diffusion Directions Imaging (DDI), to track the HMFs on clinical data, using a deterministic streamline tractography algorithm. This model allows to capture several diffusion directions within a voxel, with a low number of parameters [START_REF] Stamm | Diffusion Directions Imaging (DDI)[END_REF][START_REF] Commowick | Multifiber deterministic streamline tractography of the corticospinal tract based on a new diffusion model[END_REF].

The two goals of this paper are: 1) to evaluate the ability of this new diffusion model and of two other standard models (DTI/ODF) to track the left and right HMFs, using a common tractography algorithm, in a clinical setting, i.e. with few diffusion gradients (typically, less than 15) and low b-values (typically, less that 1000s/mm 2 ), and 2) to study whether the number of tracked HMFs in one hemisphere is different from that in the contralateral hemisphere, for right-handed and left-handed subjects. Note that we do not test multi-tensor models, as these have been shown to be unable to provide a unique solution in the context of single-shell (one unique b-value) acquisitions [START_REF] Scherrer | Why multiple b-values are required for multi-tensor models. Evaluation with a constrained log-Euclidean model[END_REF], as is the case here and in most standard clinical protocols.

In Section 2.1, we present the three tested models, and we outline our implementation thereof. In Section 2.2, we describe the common multiple fibers deterministic streamline tractography algorithm we use for these three models. The same algorithm is used to make sure that the subsequently reported results can be interpreted as differences in models, rather than differences in dMRI sequences or tractography algorithms. The data are described in Section 2.3, and we perform statistical tests and numerical evaluation in Section 3, before discussing these results, concluding and giving some perspectives in Section 4.

Material & Methods

Diffusion models

⊲ DTI model: The Gaussian model assumes that the diffusion process can be captured by a tensor (6 parameters), which is proportional to the covariance matrix of the unknown Gaussian pdf. The tensor T is parametrised as T = exp(M ), where M is an unknown 3 × 3 symmetric matrix, and its estimation is done using a least-squares (LS) fitting on the raw DWI intensities [START_REF] Fillard | Clinical DT-MRI estimation, smoothing, and fiber tracking with log-Euclidean metrics[END_REF]. The LS criterion is optimised numerically using the NEWUOA algorithm [START_REF] Powell | The NEWUOA software for unconstrained optimization without derivatives[END_REF]. Within a given voxel, the single putative fiber direction is considered to be aligned with the direction of the eigenvectors associated to the largest eigenvalue of the tensor. The tractography algorithm uses a log-Euclidean interpolation scheme [START_REF] Arsigny | Log-Euclidean metrics for fast and simple calculus on diffusion tensors[END_REF].

⊲ ODF model: The orientation distribution function (ODF) describes the orientational structure of the diffusion function [START_REF] Tuch | Q-ball imaging[END_REF]. The raw DWI intensities are modelled with a modified basis of spherical harmonics, whose c coefficients are estimated using a LS fitting including a Laplace-Beltrami regularisation term. The number of unknown coefficients depends on the order l of the basis: c = (l + 1)(l + 2)/2. This LS problem has a closed-form solution, from which the optimal ODF (or to be precise, an approximation thereof) can be computed analytically using the Funk-Hecke theorem [START_REF] Descoteaux | Regularized, fast, and robust analytical Q-ball imaging[END_REF]. Then, ODF sharpening is performed using spherical deconvolution to compute the fiber ODF [START_REF] Descoteaux | Deterministic and probabilistic tractography based on complex fibre orientation distributions[END_REF]. Within a given voxel, the putative fiber directions are selected as the local maxima of the normalised, sharpened ODF for which the ODF function value is above a user-specified threshold set here to 0.1. These local maxima are computed using NEWUOA (with starting points homogeneously distributed on the unit sphere), and they are sorted according to their ODF function value. The tractography algorithm uses a trilinear interpolation scheme [START_REF] Lenglet | Mathematical methods for diffusion MRI processing[END_REF].

⊲ DDI model: The diffusion function is modelled as a mixture of distributions with a common parametric form [START_REF] Stamm | Diffusion Directions Imaging (DDI)[END_REF]. The number of mixture components is that of the number m of different fiber directions within the voxel. In essence, the pdf of each of these distributions is defined as the convolution of a von Mises-Fisher pdf (which models the direction of the fiber) and of a centered cylindrical Gaussian pdf (which models the amplitude of the diffusion along the fiber). The covariance matrix of the Gaussian distribution is actually a function of the two parameters of the von Mises-Fisher distribution: the mean direction µ (unit vector) and the concentration parameter κ. A scalar parameter λ, assumed to be identical for all fiber directions, completes the model. Therefore, to allow for m fibers, the DDI model requires 3m + 1 parameters. An anisotropy value ξ, akin to the fractional anisotropy (FA) (resp. the generalised FA (GFA) [START_REF] Tuch | Q-ball imaging[END_REF]) in the DTI (resp. ODF) model, is also defined. The 3m + 1 unknown parameters are estimated using a LS fitting on the raw DWI intensities, and this optimisation is performed using NEWUOA. Within a given voxel, the putative fiber directions are a natural output of this model, sorted according to their diffusion function values. Note that, as of now, the number of fibers m in this model is set to 2, and automated estimation thereof will be a topic of future investigation. The tractography algorithm uses a trilinear interpolation scheme.

The common multiple fibers tractography algorithm

Our goal is to track the fibers linking two ROIs. Our deterministic streamline algorithm can be viewed as an extension of the original FACT method [START_REF] Mori | Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF], adapted to ODF and DDI models, using a breadth-first-type search. It must be made clear that for the DTI model, the tractography is led without considering multiple directions; we omit this important detail below for the sake of clarity.

Starting from one of the two ROIs, we define n starting points within each voxel of the ROI. The DTI/ODF/DDI models at these starting points are estimated using the previously described interpolation schemes. For each of these starting points, we compute the two principal putative directions (defined using the previously described sorting out procedures), we follow the first direction with a step size of l millimeters and we record the second direction for future use, as it can be indicative of crossing/kissing/merging/diverging fibers. We then reestimate the DTI/ODF/DDI models at this new spatial position (using the previously described interpolation schemes), and we compute all the putative directions for each model. Among these, we follow the one closest (i.e. with minimal angular difference) to the previously estimated first direction. A second direction, having the highest ODF/diffusion function value among the remaining putative fiber directions, is recorded for future use. The tracking of the main fiber is achieved when i) the angle between two successively estimated first directions is higher than α, or when ii) FA/GFA/ξ is lower than β, or when iii) the fiber reaches the border of a precomputed brain mask [START_REF] Smith | Fast robust automated brain extraction[END_REF]. Once this main fiber has been tracked, we perform the same tracking from all the possible crossing/kissing/merging/diverging points that have been recorded along its path. Importantly, for these trackings, the stepping rule and stopping criteria are identical as those for the main fiber, but we do not record any possible mixed fiber configuration along these secondary paths, for which we only follow the first direction at each step. The same tracking is then led on the second ROI, and only the tracts linking the two ROIS are kept for further analysis. In practice, we choose the parameters n = 1, l = 1, α = 60 degrees and β = 0.15.

Data

The data consist of dMRI, anatomical (aMRI) and functional (fMRI) MRI on 14 right-handed (8 males, 6 females) and 9 left-handed (6 males, 3 females) healthy volunteers. The mean age was 30.3 (21 to 45). Handedness was determined using the Oldfield questionnaire [START_REF] Oldfield | The assessment and analysis of handedness: the Edinburgh inventory[END_REF]. The aMRI, dMRI and fMRI data were acquired using standard sequences on a Philips Achieva 3T system:

-aMRI: T1-w 3D TFE, 184 sagittal slices of size 256×256 (1mm×1mm) and 1mm thickness. -fMRI: gradient echo EPI using BOLD contrast, 24 contiguous axial slices of size 128×128 (1.8mm×1.8mm) and 4mm thickness. The hand motor task consisted in opening and closing the hand, and was implemented in a standard block design. Motion correction, slice-timing and detection of the activation areas for both right and left hands were performed within SPM5. -dMRI: single shot EPI, 60 contiguous axial slices of size 128×128 (2mm× 2mm) and 2mm thickness. Diffusion gradients were applied in 15 non-collinear directions with b = 800s/mm 2 . Each diffusion-weighted MRI was corrected for eddy current-induced geometric distortions [START_REF] Netsch | Quantitative evaluation of image-based distortion correction in diffusion tensor imaging[END_REF] and denoised using the Rician non-local means algorithm [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio MRI: applications to DT-MRI[END_REF]. Given this low number of directions, the (modified) spherical harmonics basis of order 4 (15 parameters) was used for ODF estimation, while 6 parameters (resp. 7) were to be estimated for the DTI (resp. DDI) model.

For each subject, the aMR and fMR images were rigidly registered to the B0 image of the dMRI sequence [START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF]. A first ROI was manually delineated by an expert neuroanatomist in an axial slice on the aMRI data through the superior part of the mesencephalum, both on left and right sides. Tractography was then performed between these two anatomical ROIs and the two (left and right) cortical functional ROIs to reconstruct the HMFs using the three above-mentioned diffusion models and the previously described tractography algorithm.

Results

Connections between the ROIs

Our objective here was to evaluate whether the three diffusion models, coupled with the tractography algorithm, were able to connect fully, partially, or not at all, the anatomical and functional ROIs. First of all, we split each functional ROI into a medial and a lateral area, the latter corresponding to the extremity of the hand representation on the motor homunculus, i.e. the thumb. Then for each model DTI/ODF/DDI, each hemisphere, and each of the 23 subjects, we computed a discrete score of 0, 1 or 2 depending on the quality/strength (subjectively based on the number of fibers) of the connection between the anatomical ROI and the medial part of the functional ROI; in a word, we estimated 3 × 2 × 23 = 138 scores. Similarly, we computed another set of 138 scores for the connection with the lateral part of the functional ROI. At last, the overall number of fibers composing the reconstructed HMFs, i.e. linking the two (anatomical and functional) ROIs, was also computed.

The Pearson χ 2 test is particularly adequate to handle such qualitative, discrete scores. We performed pairwise Pearson χ 2 tests with a significance level of 0.05, corrected for multiple comparisons (Bonferroni) to compare DTI vs ODF, ODF vs DDI, and DDI vs DTI for the medial and lateral areas. To compare the overall number of fibers, we first showed that the data were not Gaussiandistributed using the Jarque-Bera test, and then we performed pairwise sign tests (which allows to test for differences in medians) with a significance level of 0.05, corrected for multiple comparisons (Bonferroni).

The p-values are reported in Tab. 1, left, and mainly show that i) the ODF model was able to track more medial fibers than the DTI model, but as many lateral fibers, and that ii) the DDI model did not track more medial fibers than the ODF model, but did track more lateral fibers, at the 0.05 significance level. These two results are confirmed by the sign test on the overall number of fibers.

Asymmetry of the hand motor fibers

Our objective here was to evaluate whether the three diffusion models, coupled with the tractography algorithm, were able to show significant differences (in terms of number of fibers) between the left and right reconstructed HMFs, in right-handed [START_REF] Ludeman | Diffusion tensor imaging of the pyramidal tracts in infants with motor dysfunction[END_REF] and left-handed (9) subjects. We pooled males and females for increased statistical power. After showing that the data were not Gaussiandistributed using the Jarque-Bera test, we performed pairwise sign tests with a significance level of 0.05, corrected for multiple comparisons (Bonferroni).

The p-values are reported in Tab. 1, right, and mainly show that i) none of the models is able to show differences in left-handed subjects, but that ii) the DDI model is the only model allowing to show that the bundle of HMFs is larger (in terms of number of fibers) in the left hemisphere than in the right for right-handed subjects.

Medial

Lateral # Fibers DTI vs ODF 2.6 × 10 -9 3.4 × 10 -1 3. : "Is there a significant difference (level=0.05) between the 3 models in recovering tracts reaching the medial and lateral areas of the functional ROI, in each hemisphere, for the 23 subjects? And in the overall number of tracked fibers?". Right table: "Are the 3 models able to show significant differences (level=0.05), in terms of number of fibers, between the left and right tracked HMFs, for the 14 right-handed (RH) and the 9 left-handed (LH) subjects?".

Discussion, Conclusion & Perspectives

In this paper, we showed that i) the DDI model outperforms both the DTI and the ODF models to track the HMFs (see Fig. 1), and, maybe more importantly, that ii) the DDI model is the only model able to show that the number of HMFs in the left hemisphere is larger than in the contralateral hemisphere for right-handed subjects, which seems to fit the intuitive idea that the hand motor fibers in the hemisphere contralateral to the dominant hand should be somewhat more developed than those in the other hemisphere. Interestingly, the DDI model, as the other two models, failed to find such a difference for left-handed subjects, which may suggest that left-handedness is not simply a mirrored righthandedness. This is already known from e.g. the notoriously differing patterns of hemispheric dominance for language between left and right-handed subjects, but to our knowledge, this is the first time such results are reported on white matter fibers, at least on clinical data. These results must now be further investigated in light of the huge literature on brain asymmetry and cerebral dominance [START_REF] Toga | Mapping brain asymmetry[END_REF]. In particular, recruiting more males/females and right/left-handed subjects for increased statistical power and population-specific analysis would be necessary to confirm these first results, of potentially important anatomical significance.

As expected, the DTI model performs very poorly when tracking the HMFs. Importantly, we stress here that we do not state that the DDI model outperforms the ODF model in general, but only in this particular experimental setting. Low angular resolution of the order-4 ODF model can explain why it is outperformed here. It would be of high interest to try to replicate our experiments on HARDI data using higher-order ODF models and improved (e.g. probabilistic) tractography methods, to further support our first findings.

Fig. 1 .

 1 Fig. 1. Tractography of the HMFs in the left and right hemispheres. From left to right: DTI, ODF and DDI. Note that we use the neurological convention, i.e. the left (resp. right) hemisphere is displayed on the left (resp. right). The DDI model is the only one able to consistently reach the lateral area within the functional ROI.

Table 1 .

 1 3 × 10 -7 ODF vs DDI 3.7 × 10 -2 7.6 × 10 -12 3.1 × 10 -7 DDI vs DTI 1.3 × 10 -14 5.5 × 10 -14 2.8 × 10 -14 p-values of the statistical tests. Left table

	RH	LH
	DTI 0.0923 1
	ODF 0.5811 0.1797
	DDI 0.0189 0.5078