Paolo Visco

Frédéric Van Wijland

Emmanuel Trizac

Collisional statistics of the hard-sphere gas

Keywords: numbers: 05.40.-a, 51.10.+y, 02.50.-r

We investigate the probability distribution function of the free flight time and of the number of collisions in a hard sphere gas at equilibrium. At variance with naive expectation, the latter quantity does not follow Poissonian statistics, even in the dilute limit which is the focus of the present analysis. The corresponding deviations are addressed both numerically and analytically. In writing an equation for the generating function of the cumulants of the number of collisions, we came across a perfect mapping between our problem and a previously introduced model: the probabilistic ballistic annihilation process [Coppex et al., Phys. Rev. E 69 11303 (2004)]. We exploit this analogy to construct a Monte-Carlo like algorithm able to investigate the asymptotically large time behavior of the collisional statistics within a reasonable computational time. In addition, our predictions are confronted against the results of Molecular Dynamics simulations and Direct Simulation Monte Carlo technique. An excellent agreement is reported.

I. INTRODUCTION

Despite the great interest that the hard sphere gas has triggered since the early days of statistical physics, there are still, to this day, simple questions which have not been deeply investigated yet. In this paper we address one such question: what is the probability that a tagged particle suffers a given number of collisions in a time t? Apart from the simplicity of the question, relevant in its own right, the collisional statistics of the hard sphere gas turns out to be useful in several purposes, such as the estimation of transport coefficients or the characterization of the dissipation in the closely related granular gases.

In a low density hard sphere gas at equilibrium, the velocities of two colliding particles are uncorrelated just before the collision: this is the molecular chaos (strosszahlansatz) statement. This remark could naively lead to the conclusion that collisions are uncorrelated random events, implying therefore that the number of collisions is simply a Poisson random variable. However, even if molecular chaos is exactly verified, the collision number is not a Poisson random variable. The non-Poissonian nature of collisional statistics has already been noticed in the literature (see e.g. [1]) but to our knowledge, has resisted analytical investigations. It is our purpose here to fill this gap.

The reason for the non-Poissonian behavior alluded to above is that the probability for a collision to take place depends on the scattering cross section of the colliding pair (which depends on the relative velocity g of the pair). Namely, for the hard sphere gas, the probability of a collision behaves roughly as g • σΘ(g • σ), where σ is a unitary vector joining the centers of mass of the two particles at contact (directed along the apse line), and Θ is the Heaviside step function. For a gas made up of particles interacting through a potential other than the hard-sphere potential, the probability of having a "collision" is of course different, which affects the distribution of the number of collisions. Of particular inter-est is the gas of particle interacting through a pair potential V (r) ∼ 1/r 2(d-1) , where d is the space dimension. Such particles are known as Maxwell molecules [2] and lead to a velocity independent probability of having a collision. Within this model, it then appears that the collisions are truly uncorrelated random events, and hence that the number of collision is a Poisson random variable.

In the next section, the free flight time distribution of a hard sphere gas is addressed, together with the distribution of free path lengths. These results pertain to N = 0 collision properties. In section III, the analysis is extended to consider the full probability P (N , t) encoding the number of collisions statistics. The large time behavior and a complementary perturbative treatment are worked out analytically. Explicit expressions are obtained for the cumulants of the number of collisions N . In order to put the theoretical predictions to the test, three types of numerical simulations are performed. The first two, Molecular Dynamics and Direct Simulation Monte Carlo (DSMC) are routinely employed in the field and beyond. The third type, of the Monte Carlo family, is discussed in IV, and specifically devised to solve the particular problem at hand. It relies on a reinterpretation of the eigenvalue equation for the generating function for the cumulants of N , in terms of a population dynamics with cloning and annihilation. A Markov chain with the desired properties in then constructed, which allows for a direct measure of several key quantities involved in the analytical treatment. In doing so, we uncover a fruitful mapping with the probabilistic ballistic annihilation model proposed in Ref [3]. The three numerical methods provide results that are in excellent agreement with the analytical predictions. Preliminary accounts of part of this work has appeared in [4], where the numerical aspect of the work was limited to its Molecular Dynamics content, and where the theory is restricted to the zeroth order of the treatment put forward here.

II. FREE FLIGHTS TIME DISTRIBUTION

A. General results

We consider a hard sphere gas in d dimensions composed of N particles of equal mass m and equal diameter σ. The gas is thermalized at some temperature T 0 in a homogeneous state of (constant) density ρ, so that the one point distribution function of the system is a Gaussian:

φ(r, v) ≡ φ(v) = ρ(2πT 0 ) -d/2 exp - v 2 2T 0 . (1)
If one decides to follow the evolution of one particle among the N , then the velocity probability distribution function (pdf) of that particle will evolve according to the homogeneous linear Boltzmann equation:

∂ ∂t f (v 1 , t) = 1 ℓ dv 2 ′ d σ v 12 • σ [f * * 1 φ * * 2 -f 1 φ 2 ] , (2)
where ℓ = (σ d-1 ρχ) -1 is a length proportional to the mean free path denoted l below, and

v 12 = v 1 -v 2 .
The factor χ is the pair correlation function at contact and is the so called Enskog correction factor [5]. In the dilute limit where the molecular chaos assumption is justified, it tends to unity. Moreover, the primed integral means that the integration has to be performed in the domain for which v 12 • σ > 0. Finally, we are employing the short-hand notation

f * * i ≡ f (v * * i , t)
, where the two star superscript refers to the precollisional velocity:

v * * 1 = v 1 -(v 12 • σ) σ , v * * 2 = v 2 + (v 12 • σ) σ .
(3) In this framework the evolution of our tagged particle is exactly a Markov process. The probability of hopping from a state with velocity v to another state of velocity in a narrow interval dv ′ around v ′ and in a time interval dt is given by W (v ′ |v)dv ′ dt, where

W (v ′ |v) = 1 ℓ |v ′ -v| 2-d (2πT 0 ) -1/2 exp - 1 2T 0 × v ′2 -v 2 + v • v ′ -v |v ′ -v| (4)
is the transition rate density per unit time (see e.g. [6]). For any general Markov process, the probability of leaving a given configuration Γ in a time interval dt is r(Γ)dt, where

r(Γ) = dΓ ′ W (Γ ′ |Γ) . (5)
This is simply the loss term of the Linear Boltzmann equation (2), and it reads for the hard sphere gas:

r(v) = ω √ 2 v 2 dT 0 1 F 1 1 2 , 1 + d 2 , - v 2 2 T 0 + e -v 2 2 T 0 1 F 1 d -1 2 , d 2 , v 2 2 T 0 , (6)
where 1 F 1 is the confluent hyper-geometric function of the first kind [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs and mathematical tables[END_REF], and ω is the collision frequency of the gas:

ω = dvφ(v)r(v) = 2π d-1 2 √ T 0 Γ(d/2) . (7)
From the Markovian property, it follows that the probability of having a given velocity v for a given time t and making a collision exactly at time t is exponential:

P F F T (t|v) = r(v) exp (-r(v)t) . (8)
The probability density of the free flight time follows:

P F F T (t) = dvφ coll (v)P F F T (t|v) . (9)
Here φ coll (v) is the velocity pdf of the colliding particles (i.e. obtained sampling the velocity of the particle only on collision):

φ coll (v) = r(v) ω φ(v) , (10)
where the normalization factor ω is exactly the collision frequency of the gas. The free flight time distribution for hard spheres was already investigated in [START_REF] Wiegel | [END_REF], but the authors used the equilibrium velocity distribution φ instead of the on collision velocity distribution φ coll as a weight in [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF], and hence obtained an incorrect result (this issue is further commented in Appendix A). Two measurements of the free flight time distribution in an event driven molecular dynamic simulation of hard disks (see section IV for more details) and in a Direct Simulation Monte Carlo (DSMC) [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF], are shown in Fig. 1, together with the expression of Eq. ( 9). The agreement is excellent. The molecular dynamics simulations are performed at a low density (ρσ 2 = 10 -2 ) to ensure that molecular chaos holds, whereas the DSMC technique relies by construction on the molecular chaos approximation, and can therefore be considered as providing a ρ → 0 benchmark. Figure 1 also provides a comparison with the exponential law that would hold for a constant r(v), as is the case for Maxwell molecules. The latter expectation [ω exp(-ωt)] is seen to hold at early times, but significantly fails at large t, a time regime that we investigate in more details in the next section.

B. Large time and high dimension analysis

Apart from simplified models (see Appendix B), it does not seem possible to obtain a simple exact analytical expression for the free flight time distribution, mainly because of the particular behavior of the hyper-geometric function involved in the expression of the collision rate r(v). In order to get a simpler approximated result, one has to look at limiting cases. For large times, the integral involved in Eq. ( 9) can be estimated through the saddle-point approximation (a.k.a. method of steepest descent, see e.g. [START_REF] Morse | Methods of Theoretical Physics[END_REF]), yielding: MD simulations in all this work are shown at density ρ = 0.01σ -2 . The full line is the numerical integration of Eq. ( 9). The dot-dashed line shows the asymptotic large time behavior of Eq. [START_REF] Chapman | The mathematical theory of nonuniform gases[END_REF], and the dashed line is the exponential distribution, associated to the Poisson distribution. Time is measured in units of the mean free time (ω = 1).

P F F T (t) ∼ exp - ωt √ 2 ω 2 1 - 2 d + ωt √ 2d -d/2 . (11)
The above function is plotted in Fig. 1, and seen to be in very good agreement with our numerical data provided t is large enough. Note that expression [START_REF] Chapman | The mathematical theory of nonuniform gases[END_REF] does not provide a normalized probability (a feature visible on Fig. 1), but only an asymptotic expansion.

Besides, a way to get an explicit expression for the collision rate r(v) is to investigate the infinite dimension limit. In this case, the saddle-point approximation yields:

r(v) ≃ ω √ 2 1 + v 2 2T 0 d . (12)
A detailed derivation of this result is given in Appendix C.

C. Free path length distribution

From the knowledge of the free flight time distribution, it is straightforward to compute the distribution of the free path length. In fact, for a particle with a velocity v, the length covered in a time t is x = |v|t = vt. Therefore, the conditional distribution of free path length is obtained as:

P F P L (x|v) = 1 v P F F T x v v = r(v) v exp - r(v) v x . ( 13
)
The free path length distribution can be obtained averaging the above expression with the weight φ coll : FIG. 2: (Color online) Plot of the free path length distribution from MD simulations (at density ρ = 0.01σ -2 ) and from a numerical integration of Eq. ( 14). In the inset, the dashed line is the dominant exponential term, Eq. ( 16).

P F P L (x) = dvφ coll (v)P F P L (x|v) . (14)
Although the conditional probability of free path length is very similar to the conditional probability of free flight time, the step of averaging over the collisional velocity distribution may drastically change the shape of the distribution. A case worthy of attention is discussed in Appendix D, for Maxwell molecules: the free flight time distribution is a pure exponential, while the free path length distribution has a stretched exponential decay. In the case of hard particles, the dominant large path behavior remains exponential, as for the free flight time. Nevertheless, the leading exponential prefactor changes. The free flight time distribution is dominated by e -r(0)t , where r(0) is the minimum of the function r(v). In the case of free path length, the dominant exponential behavior is e Ax , where A is the minimum of r(v)/v, which is obtained in the limit v → ∞:

A = lim v→∞ r(v) v = 1 l √ 2 = π d-1 2 ℓ Γ( d+1 2 ) , (15)
where l = ω/ v is the mean free path. This yields:

P F P L (x) ∼ e x √ 2l , (16)
(here the ∼ should be understood as an equivalence between the logarithms in the large x limit). Unfortunately, this expression only gives the behavior of the dominating exponential term, which becomes visible at a length scales significantly larger than the mean free path. Hence, if one analyzes the result of MD simulations as in Fig. 2, one sees that there are sub-leading (algebraic) terms in the asymptotic behavior, which still play a role and that are responsible for the mismatch between Eq. ( 16) and MD data in Fig. 2. A similar feature holds for the free flight time distribution: upon neglecting the algebraic sub-dominant correction in Eq. ( 11), the agreement displayed at late times in the inset of Fig. 1 would be spoiled.

III. THE NUMBER OF COLLISIONS

A. Uncorrelated approach

After having shown how the probability of the time between two subsequent collisions behaves, our interest goes to the probability of the sum of many such times, and in particular, the probability of having N collisions in a time t. The time interval t in which the N collisions take place can be decomposed in a sum of N + 1 (correlated) times:

t = t in + t 1 + • • • + t N -1 + t f . (17)
We are considering that a given particle begins a trajectory with a given initial velocity, and then waits a time t in before it makes a collision and gets a new velocity. During a time t 1 , it then flies freely, before colliding again etc. The successive behavior of the particle's trajectory will be a sequence of collisions, spaced by other free flight times t 2 , t 3 , . . . , t N -1 , until the last N -th collision, which will be followed by a final time interval t f which is not ended with a collision. The free flight times t 1 to t N -1 are of course distributed following the probability P F F T given by Eq. ( 9). We will denote by P in and P f the pdf of t in and t f , although their precise statistics are irrelevant for our purposes, as becomes clear below.

If one would assume these N + 1 times to be uncorrelated, then the probability P (N , t) of having N collisions in a time t could be deduced from the knowledge of the free flight time probability as follows. We first introduce the Laplace transform:

P (z) = dt e -zt P (t) , (18)
where P can be either P F F T , P in , or P f . The transform of P (N , t) can hence be expressed as:

P (N , z) = P in (z) P F F T (z) N -1 P f (z) . (19)
Of course, since P F F T cannot be computed analytically, a closed form expression for P (N , z) is not available. One interesting limit to analyze is that of large times. To get an asymptotic form of the following integral

P (N , t) = dz e zt P in (z) P F F T (z) N -1 P f (z) , (20)
we make the remark that the number of collisions increases with time linearly on average, and we define a time intensive counterpart to N , a fluctuating collision rate, which we denote by n ≡ N /t. The integral of Eq. ( 20) then reads:

P (N , t) = dz exp t(z + n ln( P F F T (z))) + O(1) .
(21) Finally, when t → ∞, the above probability behaves as:

P (N , t) ∼ e π(n)t , (22)
where π(n) is a large deviation function, and is related to P F F T through the following relation:

π(n) = min z (z + n ln( P F F T (z))) . (23)

B. Large time behavior

The result of the previous subsection leads to an approximate evaluation of the large deviation function of the number of collisions that neglects temporal correlations. We will now investigate this large deviation function keeping into account these correlations. Let us define the joint probability f (v, N , t) of having a velocity v and having suffered N collisions up until time t. In a homogeneous state, the time evolution of the above defined probability is governed by a slightly different form of the linear Boltzmann equation:

∂ t f (v 1 , N , t) = dv 2 ′ d σ(v 12 • σ)× [f (v * * 1 , N -1, t)φ(v * * 2 ) -f (v 1 , N , t)φ(v 2 )] . (24)
The function φ(v) still represents the one point velocity pdf of the gas, which is in equilibrium. Moreover, since the particle whose collisions we are counting has the same velocity pdf, we enforce that:

∞ N =0 f (v, N ) = φ(v) . ( 25
)
It is useful to introduce the generating function of f as:

f (v, λ, t) = ∞ N =0 e -λN f (v, N , t) . (26)
It can be seen that f evolves according to

∂ t f (v 1 , λ, t) = dv 2 d σΘ(v 12 • σ)(v 12 • σ)× [e -λ f (v * * 1 , λ, t)φ(v * * 2 ) -f (v 1 , λ, t)φ(v 2 )] . (27)
The large time behavior of the solution of the above equation is dominated by the largest eigenvalue µ(λ) of the evolution operator:

µ(λ) f (v 1 , λ) = L λ f (v 1 , λ) = dv 2 d σΘ(v 12 • σ)(v 12 • σ)× [e -λ f (v * * 1 , λ)φ(v * * 2 ) -f (v 1 , λ)φ(v 2 )] , (28
) where f denotes the eigenfunction of L λ associated with µ:

f (v, λ, t) ∝ e µ(λ)t f (v, λ). (29)
Moreover, since P (λ) ∼ e µ(λ)t , one sees that µ(λ) is proportional to the cumulant generating function:

N p c = t(-1) p ∂ p µ ∂λ p λ=0 . (30)
Furthermore, P (N ) ∼ dλe λN P (λ), and hence, applying the saddle point method for t → ∞, one has that

P (N , t) ∼ exp(tπ(n)) , (31)
where n = N /t and the large deviation function π is related to µ through a Legendre transform:

π(n) = min λ (µ(λ) + λn) . ( 32
)
1. Large λ behavior Before trying to solve Eq. ( 28) extending the methods of kinetic theory, we shall first provide exact results, which can be extracted from the analysis of asymptotically large values of λ. The evolution operator appearing in Eq. ( 28) can be cast in the form:

L λ = L 0 + e -λ L 1 , (33)
where

L 0 f (v 1 , λ) = -dv 2 ′ d σ(v 12 • σ) f (v 1 , λ)φ(v 2 ) = -r(v 1 ) f (v 1 , λ) , (34a) L 1 f (v 1 , λ) = dv 2 ′ d σ(v 12 • σ) f (v * * 1 , λ)φ(v * * 2 ) .
(34b) For large values of λ, the coefficient e -λ plays the role of a small parameter, and therefore the eigenvalue equation ( 28) can be solved in perturbation theory. We will therefore try to get an expression of µ(λ) for large λ of the form:

µ(λ) = µ (0) + e -λ µ (1) + O(e -2λ ) . (35)
To zero-th order, the largest eigenvalue of L λ is given by the maximum of the function -r(v), whose expression is written in Eq. ( 6). The maximum of this function occurs at v = 0, and

µ (0) = -r(0) = - ω √ 2 . ( 36
)
The eigenfunction associated with this eigenvalue is indeed a delta function, centered in v = 0. In order to get the first order correction to the eigenvalue (36) one has to project the zero-th order eigenfunction on the operator proportional to the small parameter:

µ (1) = dv 1 L 1 δ(v 1 ) = ω √ 2 . ( 37
)
Finally one finds that for large λ, the largest eigenvalue µ behaves like:

µ(λ) = ω √ 2 (e -λ -1) + O(e -2λ ) . (38)
Hence, the probability of having N collisions behaves, for values of N small with respect to its average N = ωt, as a Poisson distribution with a frequency equal to ω/ √ 2:

P (N ) ∼ e -ωt √ 2 N ! ωt √ 2 N , for N ≪ ωt . (39)
For such a distribution, the large deviation function π easily follows:

π(n) = n -n log(n √ 2/ω) -ω/ √ 2 (40)
The result obtained in Eq. ( 39) is compatible with the large time behavior embodied in Eq. ( 11): the time derivative of the probability of having N = 0 collisions is of course (minus) the free flight time distribution function. We emphasize that both results (39) and ( 11) hold for large times. In addition, we note that the only space dimension dependence involved is through the collision frequency ω [see Eq. ( 7)].

An approximate solution

In order to get an approximate expression for the largest eigenvalue µ(λ), we will suppose the associated eigenvector f (v, λ) to be a Gaussian with a given temperature T (λ). We expect that this approximation will provide accurate results for small values of λ, given that for λ = 0, f (v, 0) is exactly a Gaussian with a temperature T 0 , but it is not a priori a systematic approximation for larger values of λ. Projecting the Boltzmann-like equation (28) onto the first two velocity moments of the eigenfunction we want to compute, we are left with two closed equations for µ(λ) and for T (λ):

µ(λ) = ν 0 , µ(λ)dT (λ) = ν 2 , (41)
where the ν n are the collisional moments (the expression of the first ones is given in Appendix E):

ν n = dv 1 dv 2 ′ d σv n 1 (v 12 • σ)× f (v * * 1 , λ)φ(v * * 2 ) -f (v 1 , λ)φ(v 2 ) . (42)
Solving simultaneously Eqs. (41), one obtains:

µ(λ) = - ω √ 2 (1 -e -λ ) 1 + T (λ) T 0 , (43a)
T (λ) = √ 2T 0 √ 1 + e λ . ( 43b
)
One may notice that this result satisfies some of the previous requirements obtained from the asymptotic large λ analysis.

In particular µ(∞) = -ω/ √ 2, as in the previous subsection. Moreover, the fictitious temperature T (λ) vanishes for infinite λ, meaning that the eigenfunction associated with µ does indeed tend towards a delta function. Nevertheless, while these features concerning the zero-th order perturbation results of the previous section are fulfilled, the behavior of Eq. ( 43) is different at the next order. In fact, one can see from Eq. (43) that, for λ → ∞:

µ(λ) ∼ - ω √ 2 -ω √ e -λ 2 + O(e -λ ) , (44)
to be compared with Eq. ( 38). This is a deficiency of the Gaussian approximation.

In order to get an expression of π(n) one should compute the Legendre Transform of µ(λ). This seems not feasible analytically, but can be achieved numerically. The function π(n) is plotted in Fig. 3, together with the result of uncorrelated estimation of Eq. ( 23), as well as two Poisson distributions, of average ω and ω/ √ 2. One can see that the two estimations carried out are in very close agreement, apart from a slight difference in the behavior of the right tail. Moreover, the two new results presented here are clearly different from the Poisson distribution with mean ω [denoted Poisson (1) since ω = 1 in Fig. (3)], which is narrower, and underestimates extreme events, characterized by either very few or many collisions with respect to typical realizations.

Sonine perturbation and cumulants

We now further exploit the property that for λ = 0 the solution of the Boltzmann equation ( 28) is exactly a Gaussian. This incites us, for λ close to 0, to search for a solution f (v, λ) as a small perturbation of a Gaussian distribution. One of the most useful expansions in kinetic theory is the Sonine polynomials expansion 1 . In practice this expansion consists in looking for solutions expressed as a Gaussian times a series of Sonine polynomials, denoted S n (x):

f (v, λ) = (2π) -d/2 exp - v 2 2T ∞ n=0 a n S n - v 2 2T .
(45) The first Sonine polynomials are:

S 0 (x) = 1 , (46a)
1 For a general treatment of Sonine polynomials in kinetic theory, see e.g. [START_REF] Chapman | The mathematical theory of nonuniform gases[END_REF] and references therein.

S 1 (x) = -x + d 2 , (46b)
S 2 (x) = 1 2 x 2 - d + 2 2 x + d(d + 2) 8 . (46c)
These polynomials have the property of being orthogonal with respect to a Gaussian measure in dimension d, and are therefore related to Laguerre polynomials. From this feature, it follows that the coefficient of the series (45) a 0 = 1 and that a 1 = 0. Hence, the first nontrivial correction to the Gaussian approximation comes from the term proportional to a 2 in the expansion (45). The procedure to get an estimate of the coefficient a 2 consists of solving a closed system of equations obtained projecting the equation (28) onto the first velocity moments :

µ(λ)m n = ν n , (47)
where:

m n (λ) = dv f (v, λ) , (48)
and ν n (λ) denotes the collisional moment of order n defined in Eq. ( 42). Truncating the expansions (45) up to the second Sonine polynomial, one gets for the moments m n :

m 0 = 1 , m 2 = dT (λ) , (49a)
m 4 = d(d + 2)T 2 (λ) (1 + a 2 (λ)) . (49b)
The expression of the first collisional moments in the Sonine approximation is given in the Appendix E. Hence, taking the moment equation (48) for n = 0, 2 and 4 gives a closed system of equations for µ(λ), T (λ) and a 2 (λ). This system can be solved perturbatively expanding its solutions in power series around λ = 0. The cumulants of N obtained by means of the latter expansion turn out to be remarkably accurate. The expansion of µ(λ) up to the 3 rd order in λ gives access to the first three cumulants, which read:

N c ωt = 1 , (50a)
N 2 c ωt = 9 64 8 + 1 4d + 3 , (50b)
N 3 c ωt = 28d(64d(320d + 729) + 35775) + 257391 8192(4d + 3) 3 , (50c)
As can be noted from these values, the variance of the number of collisions already deviates by more than 12.5% from its Poisson value (equal to unity).

IV. NUMERICAL RESULTS

In this section, we compare the theoretical results obtained in the previous sections against numerical simulations. In addition to the two aforementioned techniques of Molecular Dynamics and Direct Simulations Monte Carlo (see Fig. 1), we have shaped a third numerical tool, constructing a Monte Carlo algorithm in order to directly solve the eigenvalue equation (28). We start by setting the stage for the latter method, before briefly commenting on the Molecular Dynamics simulations used to measure the statistics of the number of collisions suffered by tagged particles.

A. Monte Carlo approach

In order to derive an algorithm for solving Eq. ( 28), it is useful to rewrite this equation in the form:

µ(λ) f (v, λ) = e -λ I[ f |φ] + (1 -e -λ )A[ f |φ] , (51)
where

I[f |g] = dv ′ d σ(v 12 • σ)[f * * 1 g * * 2 -f 1 g 2 ] , (52)
is the collision integral describing the elastic collision between two particles having respectively a velocity pdf f and g. We are using the short-hand notation

f * * i ≡ f (v * * i ). The functional A[f |g]
is the loss term of the above collision integral. It actually describes the statistics of hard spheres which annihilate after each collision:

A[f |g] = -dv ′ d σ(v 12 • σ)f 1 g 2 . ( 53
)
In the context of the nonlinear Boltzmann Equation, i.e. when φ = f in Eq. ( 51), and for positive values of λ, the rhs of Eq. ( 51) exactly describes the time evolution of the velocity pdf of the probabilistic ballistic annihilation process. This model, introduced in [3], consists in a system of N particles 2, which should asymptotically dominate for large λ (see Fig. 6). The inset shows a zoom near λ = 0, where the difference between Gaussian and Sonine orders can only be appreciated for the larger values of λ displayed. which move ballistically, and interact when at contact. The interaction may be an elastic collision, with probability e -λ , or an annihilation (with probability 1e -λ ). In the context of the linear Boltzmann equation, this process can be extended in the following way. Consider a set of N independent systems, one of each is just a single particle, characterized by its velocity v i (t) (i = 1, . . . , N ), assuming spatial homogeneity. If each of these systems (particles) evolves in a hard sphere gas (in the thermodynamic limit) thermalized at temperature T 0 , with a one point velocity pdf φ, and can both collide elastically (with probability e -λ ), or annihilate (with probability 1e -λ ), then, the reduced one point velocity pdf g(v, t) = Σ N i δ(vv i (t)) will verify the Boltzmann Equa- tion

∂g(v, t) ∂t = L λ g(v, t) . (54)
Since N (t) = dvg(v, t) is the number of particles at time t, and since we know that for long times g(v, t) ∼ e µt , it is clear that for long times one has that N (t) ∼ N 0 e µt . Note that the present interpretation implicitly assumes positive values of λ. Moreover, the particles can only annihilate or collide; then the total number of particles can only decrease, and hence µ must be negative. The above observations provide a numerical tool for measuring µ(λ) as the decay rate of the total number of particles. The main difficulty, which we have successfully addressed, with the above algorithm, is that the number of particles constantly decreases, and there is no steady state but the trivial state N = 0. To circumvent this practical difficulty (that would lead to somewhat noisy statistics in the simula-tions), we have introduced an external source of particles (systems), acting in such a way that the total number of particles is conserved. If every time that an annihilation takes place, a new particle is inserted as the clone of one of the N -1 remaining particles (chosen uniformly among this population), then the evolution equation of g(v, t) reads:

∂ t g(v, t) = L λ g(v, t) + sg(v, t) , (55)
where s is a constant rate. At late times, if the above equation has a steady state, the largest eigenvalue of the operator L λ +s vanishes, and hence one has that µ(λ) + s = 0. Finally one can measure µ(λ) simply as (minus) the steady state average of the number of particles injected by the external source. So far, we showed how to construct a Markov chain in order to simulate the eigenvalue equation ( 28) for positive values of λ. For negative values of λ the procedure is almost identical. Introducing a new time scale t = (2e -λ -1)t, Eq. ( 54) can be rewritten as:

∂ tg(v, t) = e -λ 2e -λ -1 I[g|φ] - e -λ -1 2e -λ -1 A[g|φ] . (56)
All the previous physical interpretations and remarks still hold, except that now instead of having annihilation (with probability (e -λ -1)/(2e -λ -1)), one has duplication (or cloning). Hence µ(λ) will be positive, and one can add an external source in order to remove particles when new particles are created. Summarizing, the algorithm proceeds as follows:

(o) The velocity of the N particles are stored in a N × d matrix. A scalar s is set to 0. A (small) time step dt is chosen.

(i) A particle is chosen randomly in the population with uniform probability. Its velocity is denoted v 1 .

(ii) An interaction is accepted with a probability ∝ v 12 • σΘ(v 12 • σ) × dt, where σ is a random direction in d dimensions, and v 2 a d-dimensional zero-mean Gaussian random variable of variance T 0 . In practice, dt has to be chosen in such a way that |v 12 • σ|dt is always smaller than (or equal to) one.

(iii) When the interaction is accepted, if λ > 0 (resp. λ < 0), the particle will have a post-collisional velocity v * 1 with probability e -λ resp. e -λ 2e -λ -1 ) or will be removed (resp. duplicated) otherwise.

(iv) If the particle has been removed (resp. duplicated) in step (iii), one of the N -1 remaining particles is chosen randomly and uniformly, and is duplicated (resp. removed). s is increased by 1.

(v) Time is increased by an amount dt (resp. d t).

(vi) Back to (i).

This algorithm bears some similarities with the approach proposed in [START_REF] Giardinà | [END_REF], also intended to directly measure large deviation functions. Nonetheless, the version proposed here is more inspired by some variants of the DSMC algorithm for systems which do not conserve the total number of particles [13,14,15]. In order to check the reliability of the algorithm, we have first performed simulations in the case of Maxwell molecules, where the number of collisions is exactly distributed following the Poisson distribution (cf. the Introduction). The measurements of µ(λ) for this particular model are shown in Fig. 4, together with the generating function of the Poisson distribution. The agreement is excellent. In the simulations the temperature scale is set by the temperature of the heat bath T 0 , which we set to unity. The time scale is set by the mean free time, which we also set to unity. In the simulation data, the mean collision frequency is therefore ω = 1.

In the case of the hard-sphere model, a measurement of the largest eigenvalue µ(λ) is shown in Fig. 5. One can see that the numerical results are in very good agreement both with the Gaussian and the Sonine approximations. Figure 6 shows the large λ behavior of µ(λ) divided by the prediction (38). In this case, a good agreement between the numerical data and the theoretical predictions is found. In the framework of the above described Monte Carlo algorithm, it is also possible to measure the stationary velocity pdf g(v, t = ∞, λ) which is equal to the eigenfunction f (v, λ) associated with µ(λ). Hence one can also compare the analytical predictions for the temperature T (λ) and for the Sonine coefficient a 2 (λ) with the Monte Carlo results, see Figs. 7 and8. In Fig. 7, one can see the temperature T (λ). Here the Gaussian approximation is already able to capture the general behavior of this "effective" temperature, and for small values of λ, the Sonine corrections compare very well with the simulation data. In Fig. 8, one can see the coefficient a 2 as a function of λ. As expected, for small values of λ the expansion carried out in the previous section correctly describes the result of the simulations.

B. Molecular Dynamics

The most direct way to measure the large deviation function of the number of collisions N is of course to count it in a (numerical) experiment, and then construct its probability distribution function. To this end, we have performed Molecular Dynamics (MD) simulations of a two dimensional harddisk gas of N = 10 3 particles of diameter σ = 1 at density ρ = N/V = 10 -3 σ -2 and 10 -2 σ -2 . The particles evolve in a square box with periodic boundary conditions, and the time is measured in mean free time units, in such a way that N = t. We measured the statistics of the number of collisions N suffered by each particle. The value of the first cumulants is reported in Table I, together with the Poisson prediction, as well as the results from the Gaussian and Sonine approximations. It seems that, when time increases, the cumulants converge towards a value which is close to the Sonine predictions. However, it must be noted that even if increasing time makes finite time corrections smaller, the statistics become poorer and poorer. The agreement with the results from the Sonine approximation is very good, while the Gaussian order already provides a reliable estimation. This is further shown in Figure 9, where two measurements of the large de- viation function of the number of collisions for two different times, t = 10 and t = 50 (with σ = 1 and T 0 = 1), are compared to the numerical inverse Legendre Transform of µ(λ) at Gaussian order (given by Eq. ( 43)). Here again, the Poisson prediction is distinctly off.

V. CONCLUSION

We have shown that the collisional statistics of the hard sphere gas exhibits clear deviations from the Poisson distribution. These deviations have been consistently quantified both analytically and numerically. In the analytical treatment, the cumulant generating function µ(λ) plays a pivotal role. The eigenvalue equation defining this object can be fruitfully interpreted in terms of population dynamics with annihilation and cloning events, which eventually leads to an efficient algorithm allowing to compute the various quantities involved in the theoretical analysis. The corresponding numerical method, of Monte Carlo type, should not be confused with the more conventional Direct Simulation Monte Carlo technique, which we also implemented, and that is intended to solve a different kinetic equation (the Boltzmann equation). Finally, a third numerical method (Molecular Dynamics) was employed. These three routes provide complementary and valuable results, that strongly support our predictions.

Interestingly, the present formalism can be extended to the study of out of equilibrium systems, such as granular gases, where the question of the collisional statistics has been the focus of recent interest [17,18,19,20,21]. In particular, the strong effect of dissipation on the distribution of free flight times reported in [19,20] calls for further investigations. In addition, if one considers a gas of inelastic smooth hard spheres, kept in a steady state by a velocity independent force (as e.g. a vibrating wall of the container, or a stochastic force acting independently on each particle), then the phase space volume of the system has been reduced, after a time t, by a factor (1α) N (t) , where α is the coefficient of normal restitution [START_REF] Pöschel | Granular Gases[END_REF]. Hence one sees that in this non-equilibrium system, the number of collisions can be exactly identified, up to a constant prefactor, with the integrated phase space contraction rate. This quantity has already been the subject of many works in non-equilibrium statistical mechanics, and is often intimately related with the irreversible entropy production (see e.g. [22] and references therein). . The full line is the numerical integration of Eq. ( 9). The dashed line is the numerical integration of P W F F T , and the dotted line its large time behavior described by Eq. (A8). The inset shows the same data in semi-logarithmic scale, where it is possible to note that at large times the prediction of P W F F T is always above the true distribution. Time is measured in units of the mean free time (ω = 1).

where

ω V HP = 2π d/2 √ T 0 ℓΓ(d/2) (B3)
is the collision frequency for VHP. Due to the simpler form of the collision rate, the free flight time distribution can be expressed analytically: where the time is here in units of the mean collision rate ω V HP . Note that here the large time behavior is ∼ e -ωt/2 , which is even slower than the Hard-Sphere gas behaviour (which is ∼ e -ωt/ √

P V HP F F T (t
2 ) with respect to the Poissonian case (∼ e -ωt ). A comparison between these three distributions is shown in Fig. 11. The observed large time behaviors suggest that Maxwell and VHP models provide a upper and lower bounds for the hard-sphere model, a phenomenon reminiscent to that observed in the long time behaviour of dynamics that do not conserve the density [15].

APPENDIX C: SADDLE-POINT APPROXIMATION FOR THE COLLISION RATE

In this appendix we show how to recover Eq. ( 12) with the saddle point method. We first note that r(v) can be expressed as the difference of two different integrals: where For the first integral it is more convenient to perform first the integration over the angle θ, which leads to:

r(v) = Ω d-1 (I 1 (v) -I 2 (v)) , (C1)
I 1 (v) =
I 1 (v) = v d -1 +∞ -∞ dv 2σ 1 - v 2σ v 2 d-1 2
φ(v 2σ ) .

(C4) Then, defining the rescaled variables ṽ = v/ √ 2T 0 d and ṽ2σ = v 2σ / √ 2T 0 d one finds, to leading order in d:

I 1 (ṽ) ≃ ṽ T 0 π +∞ -∞ dṽ 2σ × exp d -ṽ 2 2σ + 1 2 ln 1 - ṽ2σ ṽ 2 . (C5)
When d is very large, the above integral is dominated by the maximum of the function inside the exponential, which turns out to be located in ṽ2σ = 0. One can then perform a series expansion around ṽ2σ = 0 up to the second order. This results in a Gaussian integral which is easily integrated, and yields: All data are reported in units of mean free path.

I 1 (ṽ) ≃ T 0 d 2ṽ 2 √ 1 + 2ṽ 2 . ( C6
As for the second integral involved in the expression of r(v) the simplest is first to perform the integration over v 2σ ; hence one obtains, to leading order in d:

I 2 (ṽ) = -T 0 2π π 0 exp d -ṽ 2 cos 2 θ + ln sin θ .

(C7) The only maximum of the function in the exponential between 0 and π is in θ = π/2. Then, expanding as usual this function around the maximum up to second order, and extending the range of integration from -∞ to +∞, one finds:

I 2 (ṽ) ≃ T 0 d 1 √
1 + 2ṽ

FIG. 1 :

 1 FIG.1:(Color online) Free flights time distribution of a hard disc gas. The circles and squares correspond to the results of Molecular Dynamics (MD) and Direct Simulation Monte-Carlo (DSMC) methods. MD simulations in all this work are shown at density ρ = 0.01σ -2 . The full line is the numerical integration of Eq. (9). The dot-dashed line shows the asymptotic large time behavior of Eq. (11), and the dashed line is the exponential distribution, associated to the Poisson distribution. Time is measured in units of the mean free time (ω = 1).

FIG. 3 :

 3 FIG.3:(Color online) Plot of π(n) obtained in several fashions (ω = 1). The solid line has been obtained from the numerical Legendre Transform of µ(λ) from Eq. (43). The dashed line corresponds to the uncorrelated result obtained in Eq. (23), while the dot-dashed line is the large deviation function of the Poisson distribution of average 1. For completeness, we also show the large deviation function of the Poisson distribution with average 1/ √ 2 (dotted line).

FIG. 4 :

 4 FIG. 4: (Color online) The circles are the numerical measurements of µ(λ) for the Maxwell model. The solid line is the generating function of the Poisson distribution. The inset shows the same data, shifted vertically by one (in such a way that µ(λ) + 1 is always positive), in semi-logarithmic scale.

FIG. 5 :

 5 FIG.5:(Color online) Numerical measurement of µ(λ) for the hardsphere model. The circles show the results of Monte Carlo simulations. The solid line is the theoretical prediction within the Gaussian approximation. The dashed line is the theoretical prediction (obtained solving numerically (47-48) in the framework of the first Sonine correction). Finally the dotted line is the generating function of a Poisson distribution of parameter ω/ √ 2, which should asymptotically dominate for large λ (see Fig.6). The inset shows a zoom near λ = 0, where the difference between Gaussian and Sonine orders can only be appreciated for the larger values of λ displayed.

FIG. 6 :

 6 FIG. 6: (Color online) Same data as in Fig. 5, divided by the generating function of a Poisson distribution of parameter ω/ √ 2.

FIG. 7 :FIG. 8 :

 78 FIG. 7: (Color online) Temperature (defined as the variance) of the eigenfunction f (v, λ) measured in Monte Carlo simulations. The solid and dashed lines are the predictions of the Gaussian and Sonine approximations respectively.

FIG. 9 :

 9 FIG. 9: (Color online) Large deviation function π(n) of the number of collisions N suffered by each single particle. Symbols show the results of Molecular Dynamics simulations. The solid line is the result in the framework of the Gaussian approximation, and the dashed line is the large deviation function of the Poisson distribution, of average 1: π(n) = nn log n -1.

FIG. 10 :

 10 FIG.10:(Color online) Free flights time distribution of a hard disc gas. The circles corresponds to the results of Molecular Dynamics (MD) simulations at density ρ = 0.01σ -2 . The full line is the numerical integration of Eq. (9). The dashed line is the numerical integration of P W F F T , and the dotted line its large time behavior described by Eq. (A8). The inset shows the same data in semi-logarithmic scale, where it is possible to note that at large times the prediction of P W F F T is always above the true distribution. Time is measured in units of the mean free time (ω = 1).

FIG. 11 :

 11 FIG.11:(Color online) Free flights time distributions for three different interaction kernels, in two dimensions. The (red) dashed line shows the result for VHP (Eq. B4), the (blue) solid line and the (green) dotted line show respectively the results for Hard Spheres and Maxwell Molecules. Time is measured in units of the mean free time (ω = 1).

dθI 2 (dθ

 2 sin θ d-2 (v cos θ)× Θ(v cos θv 2σ )φ(v 2σ ) , (C2) sin θ d-2 v 2σ × Θ(v cos θv 2σ )φ(v 2σ ) . (C3)

FIG. 12 :

 12 FIG.12:(Color online) Free path length distribution for Maxwell molecules in two dimensions. The circles are the result of a DSMC simulation, while the solid line is the numerical integration of Eq. (D3). The dashed line is the large length approximation of Eq. (D4). All data are reported in units of mean free path.

TABLE I :

 I Cumulants for the number of collisions N from MD simulations, and comparison with the Gaussian and Sonine approximations (time is measured in units of the mean-free time).

	t = 10	N c/t 1.	N 2 1.1228 c/t	c/t 1.1282 N 3
	t = 50	1.	1.1354	1.1045
	Poisson	1	1	1
	Gaussian	1	1.125	1.1289
	Sonine	1	1.1377	1.1073

Acknowledgments

The authors acknowledge useful discussions with J. Piasecki, J. M. J. van Leeuwen, M. H. Ernst, D. Frenkel, H. van Beijeren and P. Krapivsky. This work was supported by the French Ministry of Education through a ANR-05-JCJC-44482 grant.

APPENDIX A: COMMENT ON THE FREE FLIGHT TIME DISTRIBUTION

We shall give here some additional arguments on the incorrectness of Eq. ( 9), when used with a Gaussian weight φ, instead of the weight φ coll defined in Eq. [START_REF] Morse | Methods of Theoretical Physics[END_REF]. We shall refer to this (erroneous) distribution as P W F F T :

A first argument bears on the inconsistency between the above relation and the definition of the collision frequency ω = r(v) , where the brackets denote an average over a Gaussian weight. Indeed, the average time τ between two subsequent collisions (mean free time) of a given particle is equal to the inverse of the collision frequency: τ = 1/ r(v) . Besides, from Eq. (A1) the mean free time τ is obtained as:

Hence we see that expression (A1) is in contradiction with the definition of the collision frequency. On the other hand, the counterpart of Eq. (A2) with the distribution provided by [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF] and [START_REF] Morse | Methods of Theoretical Physics[END_REF], yields

which is the required relation. Second, expression (A1) explicitly slightly differs from the correct free flight time distribution in two limiting cases.

• At small values of t, P W F F T behaves as:

while the true distribution behaves as

so that even the value in t = 0 is different, although this difference is numerically small (see also Fig. 10). For instance, in d = 2 one has:

and

• The large time behavior of the probability also slightly differs from Eq. ( 11):

.

(A8) In the above expression, the leading exponential term remains the same as in Eq. ( 11), since it is only determined by the minimum of the function r(v). Nonetheless, the subleading prefactor is slightly different, as it can be appreciated in Fig. 10.

APPENDIX B: FREE FLIGHTS TIME DISTRIBUTION FOR VERY HARD PARTICLES

For the sake of completeness, we report in this appendix the result for the free flight time distribution arising in the very hard particle (VHP) model. This framework allows for an explicit analytic computation of P F F T (t). The model, introduced by Ernst and Hendriks [23], consists in a choice of the collision rate proportional to the total energy of the system. In the case of a single tagged particle, the corresponding linear Boltzmann equation reads:

B1) The velocity dependent collision rate is thus quadratic in

APPENDIX D: FREE PATH LENGTH DISTRIBUTION FOR MAXWELL MOLECULES

In this appendix, we investigate the free path length distribution for Maxwell molecules. As already mentioned, in this case the collision rate is a constant, ω, independent of the velocity of the particle. It follows then that the free flight time probability is exponential:

Besides, the free path length distribution for a given velocity v reads:

and the free path length distribution is simply the average of the above probability over a Gaussian weight (for Maxwell molecules the on collision distribution φ coll is still a Gaussian):

This last expression can be expressed analytically in terms of Meijer G-functions [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs and mathematical tables[END_REF], but here we shall focus only on the large length behavior, for which simpler expressions are available. In particular, the saddle point approximation gives:

which leads to a stretched exponential behavior at large x (∼ exp -x 2/3 ) (cf Fig. 12).

APPENDIX E: COLLISIONAL MOMENTS

Here we provide the expressions of the first collisional moments, defined by Eq. ( 42), for f a Gaussian, and a Gaussian multiplied by a Sonine Polynomial. We denote by

where L λ is the operator appearing in (28). Hence in the Gaussian approximation one has ν n ≡ ν


Laurent Riffaud
email: laurent.riffaud@chu-rennes.fr

Gabriel Eustache

Chloé Rousseau

Anne Defontaine

Eric Wodey

Nicolas Nardi

Sylvie Aillet

Laurent Pr

Riffaud

Does the participation of a senior plastic surgeon improve perioperative care in craniosynostosis repair surgery?

Keywords: Craniosynostosis, blood loss, perioperative care, craniofacial surgery, surgeon collaboration

Introduction. The real impact of the participation of other surgical specialties together with neurosurgeons on perioperative care in craniosynostosis repair surgery has not been determined.

The purpose of this study was to determine whether the participation of a second senior surgeon (plastic surgeon) during surgical repair of pediatric monosutural craniosynostosis improved perioperative medical care.

Material and methods. The authors retrospectively reviewed 2 cohorts of patients who had consecutively undergone primary repair surgery for trigonocephaly and unicoronal craniosynostosis. Infants were operated on by a single senior pediatric neurosurgeon before December 2017, and with the collaboration of a senior plastic surgeon after January 2018.

Results. Overall, 60 infants were included in the study: 29 in group 1 (single surgeon, 2011-2017), and 31 in group 2 (pair of surgeons, 2018-2021). Median surgery time was significantly shorter in group 2 than group 1: 180 vs 167 minutes; p = 0.0045. There was no significant difference between the 2 groups in blood loss or intra/postoperative packed erythrocyte transfusion. Postoperative drain output was significantly lower in group 2. Median length of hospital and intensive care stay were significantly shorter in group 2, by 1 and 2 days, respectively; p  0.0001. Volume of infused solution, diuresis, immediate postoperative hemoglobin level, hematocrit, hemostasis (platelet count, fibrinogen, prothrombin time and activated partial thromboplastin time), and return to oral feeding did not differ from one group to the other.

Conclusion.

Results confirmed our impression of an improvement in perioperative medical care.

However, the role of surgical experience and the influence of the medical/nursing staff must not be minimized in these complex surgical procedures.

Introduction

Craniosynostosis repair surgery is a major surgical procedure. Practice has varied exceedingly from one center to another, even for management of isolated craniosynostosis [START_REF] Lee | Management options of non-syndromic sagittal craniosynostosis[END_REF][START_REF] Alperovich | Is craniosynostosis repair keeping up with the times? Results from the largest national survey on craniosynostosis[END_REF][START_REF] Baker | Management of sagittal synostosis in the Synostosis Research Group: baseline data and early outcomes[END_REF]. These variations concerned surgeon profiles, methods of patient evaluation, imaging, surgical technique and timing, and follow-up [START_REF] Kestle | Variation in the management of isolated craniosynostosis: a survey of the Synostosis Research Group[END_REF]. It is well-established that management of faciocraniosynostoses needs multidisciplinary expertise and collaboration of pediatric neurosurgeons with maxillofacial surgeons for surgical procedures, especially when midface distraction is needed [START_REF] Mathijssen | Updated guideline on treatment and management of craniosynostosis[END_REF]. Similarly, cooperation of pediatric neurosurgeons with plastic/maxillofacial surgeons is recommended when facial corrections or complex skull shape corrections are necessary in isolated single-suture craniosynostoses in order to optimize cosmetic results, even if there are no minimal requirements of team members for management of these craniosynostoses [START_REF] Mathijssen | Updated guideline on treatment and management of craniosynostosis[END_REF]. Since low surgical case volumes have been associated with increased perioperative transfusion volumes in this surgery [START_REF] Fernandez | Perioperative outcomes and surgical case volume in pediatric complex cranial vault reconstruction: A multicenter observational study from the pediatric craniofacial collaborative group[END_REF], the real impact on perioperative care of the participation of other surgical specialties together with neurosurgeons during the surgical procedure has not been yet demonstrated.

The purpose of this study was to determine whether the participation of a second senior surgeon (plastic surgeon) during surgical repair of pediatric monosutural craniosynostoses improved conditions of hospitalisation.

Material and methods

Patient population

We retrospectively reviewed 2 cohorts of patients who consecutively underwent primary repair surgery for trigonocephaly and unicoronal craniosynostosis in our Department. Infants who needed an fronto-orbital "bandeau" were operated on by a single senior pediatric neurosurgeon (LR) before December 2017, and those who were treated after January 2018 were operated on both by senior J o u r n a l P r e -p r o o f pediatric neurosurgeon (LR) and senior plastic surgeon (SA). We began the study after routinely introducing use of tranexamic acid (TXA) for craniosynostosis surgery in our institution in March 2011 [START_REF] Eustache | Importance of tranexamic acid in pediatric monosutural craniosynostosis surgery[END_REF]. The last child included was operated on in July 2021. We included in our study only craniosynostoses for which the surgical technique and timing did not change over time i.e.

trigonocephaly and unilateral craniosynostoses, and we also excluded brachycephaly and atypical cranionynostoses. This study was approved by the "Collège de Neurochirurgie" IRB#1: 2022/44 n° IRB00011687. The study was designed as a monocentric before-after study.

Surgical procedure

Since 2000, the same standardized procedures had been performed for trigonocephaly and unilateral plagiocephaly [START_REF] Eustache | Importance of tranexamic acid in pediatric monosutural craniosynostosis surgery[END_REF]. These procedures were open cranial vault reconstructions since endoscopic procedures are not performed in our institution. A senior plastic surgeon (SA) started to collaborate during surgery after December 2017 when fronto-orbital advancement was required.

The contribution of this surgeon was in observance of the conditions required by our French center of reference in surgical procedures involving the facial features. Duration of anesthesia was calculated from the induction of general anesthesia to extubation, and duration of surgical procedure from skin incision to the end of skin closure.

Perioperative management

During the 2 periods, anesthetic management was homogenized by the team of pediatric anesthesiologists for general anesthesia, patient monitoring, decisions regarding blood transfusions and TXA administration [START_REF] Eustache | Importance of tranexamic acid in pediatric monosutural craniosynostosis surgery[END_REF]. TXA was used in all the patients with a standardized protocol reported previously: a loading dose of 10 mg/kg over 15 minutes was followed by a 10 mg/kg per hour infusion until skin closure, and no TXA was used postoperatively. No patient with hematologic or that of leaving hospital were decided by both the pediatric neurosurgeon and the pediatric anesthetist.

Assessment of blood loss

Numbers of patients transfused intra-and postoperatively and volume of packed erythrocytes (PE) transfusion were examined in both groups. Hemoglobin, hematocrit, platelet count, fibrinogen, prothrombin time and activated partial thromboplastin time were also analyzed. The estimated blood volume lost (EBVlost) was calculated with a formula detailed in previous reports [START_REF] Eustache | Importance of tranexamic acid in pediatric monosutural craniosynostosis surgery[END_REF][START_REF] Kearney | Craniosynostosis: an assessment of blood loss and transfusion practices[END_REF]. The EBV for children younger and older than 12 months was respectively 80 and 75 mL//kg, and the hematocrit of transfused packed red blood cells (PRBCs) was 0.65 at our hospital.

Statistical analysis

Age at surgery, sex, weight, American Society of Anesthesiologists (ASA) physical status, type of craniosynostosis, duration of anesthesia and surgery, volume of solution infused during the procedure, use of dopamine, diuresis, drain output, length of stay in intensive care unit and in hospital and postoperative day on which the patient started oral food were also recorded for statistical analysis. Quantitative variables were described as number, mean  SD, and median interquartile range [Q1; Q3]. Groups were then compared using Student test (S) or Mann-Whitney Wilcoxon nonparametric test (W). Qualitative variables were presented as numbers and percentages, and groups were compared using the chi-square (K) or Fisher's exact (F) tests.

Association with duration of surgery was performed using Spearman coefficients.

J o u r n a l P r e -p r o o f

Analyses were achieved using SAS version 9.4 (SAS Institute, Inc.), and p < 0.05 was considered for statistical significance.

Results

Patient characteristics

Overall, sixty infants were included in the study: 29 infants in group 1).

Duration of surgery

The median duration of surgical procedures was significantly lower in group 2 (pair of surgeons) compared with group 1 (single surgeon): 167 vs 180 minutes; p = 0.0045 (Table 2).

Median timesaving was 5% in trigonocephaly procedures with the pair of surgeons (172 minutes in group 2 vs 182 in group 1, p = 0.0146) and 21% in unicoronal craniosynostosis (140 minutes vs 177, p = 0.0302) (Table 3).

Estimated Blood Volume lost and blood product transfusion

Median Hemocue values before transfusion were the same in both groups. There was no statistical difference between 2 groups regarding EBVlost or intra/postoperative packed erythrocytes J o u r n a l P r e -p r o o f transfusion (Table 4). Postoperative drain output was significantly lower in group 2 (pair of surgeons).

Intraoperative and postoperative data

Median lengths of stay in hospital and in intensive care were statistically lower in the group of pair surgeons: the median timesaving was estimated respectively as one and two days; p  0.0001 (Table 2).

Volume of solution infused, diuresis, immediately postoperative hemoglobin, hematocrit, hemostasis (platelet count, fibrinogen, prothrombin time and activated partial thromboplastin time), and the return to receiving food orally were not different from one group to the other.

Multivariate analysis

Spearman coefficient values demonstrated that duration of surgical procedures positively influenced drain output (coefficient > 0.5), but was slightly correlated to the median volume of PE transfusion or the length of stay in hospital and in intensive care (Table 5).

Discussion

During the last two decades, our impression was a progressive improvement in perioperative care and management of child craniosynostosis patients, which appeared to result not only from the surgeon's level of experience. Consequently, we initially conducted a study assessing the impact of the routine use of TXA in our practice after its initiation in March 2011: we observed that TXA not only significantly reduced blood loss and the need for transfusions, but also decreased the hospital length of stay and minimized medical care costs, thus having major impact on J o u r n a l P r e -p r o o f postoperative care. This result was in agreement with the subjective impression of the medical staff: pediatric anesthesiologists, intensive care unit doctors and the neurosurgeon.

A second question concerned the effect of cooperation with a senior plastic surgeon, which we initiated a few years ago in some cases of craniosysnostosis repair surgery: did this collaboration influence our surgical performance, and consequently perioperative medical care? Although retrospective, our study allowed a reliable comparison of the two practices: a single senior pediatric neurosurgeon or a pair of senior surgeons, in the particular context of two complex monosutural craniosynostoses. Results from statistical analysis suggested a positive impact of collaboration: duration of general anesthesia and of surgical procedures remained significantly shorter, both for unicoronal craniosynostoses and trigonocephaly. Despite the significant reduction in operating time observed in the group of patients operated on by the pair of surgeons, we did not observe any improvement in blood loss and transfusion. This finding was probably the consequence of several factors but we assume that the duration of surgical procedures was too limited to observe a positive influence on blood loss, given the important benefit brought about by the use of TXA. Another significant consequence was shorter hospital stay in the group of patients operated on by the pair of surgeons. This result needs to be seen in the context of multivariate analysis: a shorter surgical procedure mainly influenced output drainage, with little impact on other factors such as blood loss, PE transfusion, length of stay in hospital or in intensive care. By reducing blood loss and thus output drainage, routine use of TXA may already have changed our practice over the years: drains have progressively been removed earlier (on day 1) reducing postoperative pain and therefore hospital stay [START_REF] Eustache | Importance of tranexamic acid in pediatric monosutural craniosynostosis surgery[END_REF]. In other words, we suppose that the progressive adjustment of our habits may have induced a potential bias, and the positive role of TXA should not be underestimated in the interpretation of the results. Furthermore, the role of surgical experience must not be minimized in this 10-year overall study. We have already demonstrated the relationship between surgical practice J o u r n a l P r e -p r o o f and surgical experience in live surgery by recording the surgical processes: surgical expertise is synonymous with shorter procedures and greater efficiency in manual gesture [START_REF] Riffaud | Recording of surgical processes: a study comparing senior and junior neurosurgeons during lumbar disc herniation surgery[END_REF]. On the other hand, it has been demonstrated that the learning curve generally reached a plateau after about fifty cases in craniofacial procedures firstly [START_REF] Sakamoto | The learning curve of a craniofacial surgeon: A retrospective study of 100 consecutive cases[END_REF] and secondly the senior neurosurgeon in this report already had more than 10 years' experience of craniosynostosis surgery at the beginning of this study. For these reasons, we suppose that the influence of experience on results was potentially real but limited. In a previous study comparing monosutural craniosynostoses operated on by the same pediatric neurosurgeon, the median duration of operation was not statistically different between 2 successive groups of patients over a 7-year period [START_REF] Eustache | Importance of tranexamic acid in pediatric monosutural craniosynostosis surgery[END_REF].

A further unanswered issue is whether the contribution of a senior plastic surgeon influenced the cosmetic results in this population? This question remains scientifically very difficult to demonstrate, with a strong element of subjectivity. However, our opinion is largely positive: the emulation resulting from working in partnership, the intellectual stimulation of potential exchange looking for new solutions and concepts, in association with the pediatric anesthetic staff, have indubitably contributed to improving both cosmetic results and overall medical care. Furthermore, the presence of 2 senior surgeons was an important psychological factor generating confidence in the families and stimulating new ideas in the operating theatre.

Study limitations

Data were collected retrospectively, which suggested possibilities of inaccuracy. This being so, there was very little missing data, and all the recovery was done by a single researcher, which limited potential imprecision of data collection. Furthermore, this study compared a homogeneous population of patients operated on in a standardized manner in terms of surgical technique and physician team.

Conclusion

To the best of our knowledge, this study is the first to evaluate the added benefit of the presence of a second senior surgeon in craniosynostosis repair. The results confirmed our impression of an improvement in perioperative medical care. However, the role of surgical experience and the influence of the medical/nursing staff must not be minimized in these complex surgical procedures.

  154.0; 190.0] 140.0 [130.0; 165.0] 0.0302 (W) W: Mann-Whitney Wilcoxon test Median values are presented as the median [Q1; Q3]. Boldface type indicates statistical significance. Page 18 of 19 J o u r n a l P r e -p r o o f

Table 1

 1

	Patient characteristics
	Parameter	Group 1 (single surgeon)	Group 2 (pair of surgeons)	P-value
		2011-2017	2018-2021
	Patients, n	29	31
	Median age, months	8 [8; 10]	8 [8; 10]	0.7436 (W)
	Sex M/F, n	25/4	19/12	0.0292 ( 2 )
	Type of craniosynostosis, n			0.3152 ( 2 )
	Trigonocephaly	15	20
	Unicoronal	14	11
	Mean weight, kg	8.2  1.2	8.6  1.2	0.2636 (S)
	ASA class I/II, n	25/4	25/6	0.7324 (F)
	Preoperative biological values
	Mean hemoglobin, g/dL	11.7  0.7	11.7  0.7	0.8592 (S)
	Mean hematocrit, %	35.1  2.0	34.4  1.9	0.1567 (S)
	Mean ERCV, mL/kg	28.1  1.6	27.3  1.5	0.0878 (S)
	Median platelet count,	346.0 [305.0; 425.0]	405.0 [356.0; 458.0]	0.0357 (W)
	giga/L
	Mean thrombin time, %	99.8  12.3	102.5  11.9	0.4050 (S)
	Median aPTT	1.1 [1.0; 1.2]	1.1 [1.0; 1.1]	0.4221 (W)
	F: Fisher's exact test; S: Student's test;  2 : Chi-square test; pts: patients; W: Mann-Whitney Wilcoxon
	test; ERCV: estimated red cell volume; aPTT: activated partial thromboplastin time
	Median values are presented as the median [Q1; Q3]. Boldface type indicates statistical significance.

Table 3

 3 Duration of surgical procedure according to type of craniosynostosis

	Parameter	Patients, n	Group 1
			(single surgeon)
			2011-2017

Table 4

 4 Intra-and postoperative data regarding blood product transfusion, estimated blood loss and drain output

	Parameter	Group 1	Group 2	P-value
		(single surgeon)	(pair of surgeons)
		2011-2017	2018-2021
	Intraoperative transfusion, n (%)	27 (93%)	28 (90%)	1.0000 (F)
	Postoperative transfusion, n (%)	7 (20%)	10 (32%)	0.4855 (K)
	Median volume of PE transfusion,
	mL/kg
	Total	20.0 [17.0; 20.5]	20.0 [15.0; 30.5]	0.7883 (W)
	Intraoperative	20.0 [15.0; 20.0]	20.0 [11.0; 20.0]	0.6518 (W)
	Postoperative	0.0 [0.0; 0.0]	0.0 [0.0; 15.0]	NA
	Median EBVlost, mL/kg	38.3 [26.8; 51.8]	44.6 [32.2; 54.3]	0.2974 (W)
	Median drain output, mL/kg	13.9 [10.5; 21.6]	6.1 [4.5; 10.2]	0.0005 (W)
	Median hemocue pretransfusion,	8.4 [7.8; 9.3]	8.4 [7.1; 8.7]	0.1814 (W)
	g/dL

PE: packed erythrocytes; EBV: estimated blood volume F: Fisher's exact test; W: Mann-Whitney Wilcoxon test Median values are presented as the median [Q1; Q3]. Boldface type indicates statistical significance.

Acknowledgment

The authors thank Mrs. Deirdre McKeown for her help with the English language.

This work did not receive any specific grant from funding agencies in the public, commercial, or non-profit sectors.

Disclosure of interest

The authors report no conflict of interest concerning the materials or methods used in this study or the findings specified in this paper. J o u r n a l P r e -p r o o f

J o u r n a l P r e -p r o o f