Alexis Poncet

Olivier Bénichou

Vincent Démery

Gleb Oshanin

Bath-Mediated Interactions between Driven Tracers in Dense Single-Files

Single-file transport, where particles cannot bypass each other, has been observed in various experimental setups. In such systems, the behaviour of a tracer particle (TP) is subdiffusive, which originates from strong correlations between particles. These correlations are especially marked when the TP is driven and leads to inhomogeneous density profiles. Determining the impact of this inhomogeneity when several TPs are driven in the system is a key question, related to the general issue of bath-mediated interactions, which are known to induce collective motion and lead to the formation of clusters or lanes in a variety of systems. Quantifying this collective behaviour, the emerging interactions and their dependence on the amplitude of forces driving the TPs, remains a challenging but largely unresolved issue. Here, considering dense single-file systems, we analytically determine the entire dynamics of the correlations and reveal out of equilibrium cooperativity and competition effects between driven TPs.

The motion of particles in narrow channels in which particles cannot bypass each other is known as singlefile diffusion. Such systems have been studied in various experimental setups with zeolites [1,2], micro- [3,4] and nano-channels [5], and simulations of carbon nanotubes [6]. Key features, on the theoretical side, involve the existence of a subdiffusive scaling for the mean position of a given particle [7,8], and strong correlations between particles [9].

The basic phenomenology of the single-file transport is well-captured by the symmetric exclusion process (SEP). In this paradigmatic model of crowded equilibrium systems, particles perform symmetric random walks on a one dimensional lattice with the constraint of at most a single occupancy of each lattice site. Different facets of the SEP have been scrutinised (see Refs. [8,[10][11][12][13][14][15]), including several important extensions to out-of-equilibrium situations. In particular, the mean displacement [16], as well as all higher-order cumulants [17] of an unbiased tagged particle (TP) placed initially at the shock point of a step-like density profile have been determined. Moreover, for a SEP with a single biased TP (due to either an energy consumption or an external force), the mean displacement of the latter [18,19] and the higher-order cumulants in the dense limit [20] have been calculated, and shown to grow sublinearly as √ t. Here, the particles accumulate in front of the TP and are depleted behind it, which results in an inhomogeneous, non-stationary spatial distribution of particles.

A general open question concerns situations when several biased TPs are introduced in an otherwise quiescent medium of bath particles. The TPs are then expected to entrain the bath particles in a directional motion, which brings the system out-of-equilibrium and gives rise to effective bath-mediated interactions (BMIs) between the TPs. Such BMIs potentially lead to self-organisation, as observed in systems as diverse as colloidal solutions [21][22][23][24][25][26][27][28][29][30], nearly-critical fluid mixtures [31], dusty complex plasmas [32] pedestrian counter flows [33].

Quantifying the emerging interactions between biased TPs and the ensuing collective behaviour is thus a key issue which however remains largely unresolved. Here, modeling a host medium as a dense SEP, we analytically determine the temporal evolution of all correlation functions, reveal intrinsically out-of-equilibrium cooperativity and competition effects between multiple TPs and quantify BMIs.

The quiescent host medium is modelled as a SEP which involves a high density ρ (ρ → 1) of hard-core bath particles performing symmetric random walks (with unit jump rate) on a one-dimensional lattice. We then tag N particles at initial positions X 0 j (see Figure 1). These TPs are biased: the j-th TP jumps to the left (resp. right) with probability (1s j )/2 (resp. (1 + s j )/2). The bias s j ∈ (-1, 1) may either be due to an "activity" of the particle or to an external force f j , in which case one has the detailed balance condition: e βfj = 1+sj 1-sj , β being the reciprocal temperature.

We aim at determining the correlations between the TPs, embodied in the so called cumulant-generating function ψ(k, t) ≡ ln e ik•Y(t) , where k = (k 1 , . . . , k N ), Y = (Y 1 , . . . , Y N ) and Y j (t) = X j (t)-X 0 j is the displace-ment of the j-th TP. The cumulants are denoted by • c and defined by the expansion

ψ(k, t) = ∞ p1,...,p N =0 (ik 1 ) p1 . . . (ik N ) p N (p 1 + • • • + p N )! Y p1 1 . . . Y p N N c . (1)
In densely populated single-files (ρ → 1), the dynamics of the system can be reformulated in terms of independent vacancies. In essence, this amounts to neglecting events where two vacancies interact simultaneously with any TP [20,34]. Our approach is based first on considering an auxiliary problem involving a single vacancy, initially at position U . By counting all the interactions of this vacancy with the TPs, we determine the probability p U (Y, t) that the TPs have displacements Y at time t. Then, for a density of vacancies ρ 0 = 1ρ → 0, the cumulant-generating function reads [35]

lim ρ0→0 ψ(k, t) ρ 0 = U / ∈{X 0 i } [p U (k, t) -1] (2)
where pU (k, t) is the Fourier transform of p U (Y, t) and can be expressed in terms of first-passage quantities of simple random walks with or without absorbing sites [START_REF] Hughes | Random walks and random environments[END_REF]. Analysis of the explicit expression of the cumulantgenerating function for densely populated single-files [35] allows us to draw a number of important conclusions, which we present below. Bath-mediated binding. As expected, at short times, the TPs move independently subject to their own biases. Our first finding is that, at large time and at high particle density, they are moving as a single TP. More precisely, in the large-time limit the N -TPs cumulants are given by Y q1 1 . . . Y q N N c = Z q1+•••+q N c for positive integer q j , where Z = N j=1 Y j /N is the displacement of the center of mass. At high density, even and odd cumulants of Z satisfy

Z(t) 2n c ρ 0 = Z(t) 2n+1 c ρ 0 S = 2t π , (3)
where S = tanh(βF/2) is the effective bias, F = N j=1 f j the effective force and the ratio • /ρ 0 is understood as the limit ρ 0 → 0. Equation (3) implies in particular that for any number of TPs and arbitrary forces Y j = Z , meaning that at large time all the TPs move like their center of mass with an effective force F (in agreement with the hydrodynamic analysis of Ref. [START_REF] Bénichou | [END_REF]).

In the following, we determine the full dynamics of the correlations between TPs, and focus for simplicity reasons on the case of two TPs.

Bath-mediated entrainment. We examine the case of a single biased TP (s 2 = 0) followed by an unbiased TP (s 1 = 0), initially separated by a distance L (see Figure 2, top), which allows us to quantify the perturbation induced by a biased tracer in a quiescent medium. While 2 c/s2 for the same set of parameters. The black line corresponds to Eq. (6). Inset of (a): variances for L = 10, s2 = 0.8, -0.8 (solid, dashes). Predictions in black [35]. Inset of (b): law of the variation of distance D = Y2 -Y1 at times 10, 10 2 , 10 3 , 10 4 (blue to red) for L = 10, s = 0.8, ρ0 = 0.05. The squares are the numerical results, the colored lines are the theoretical predictions, the black line is the asymptotic prediction [35].

τ = t/L 2 • c /(ρ 0 L) Y 1 Y 2 c Y 1 Y 2 2 c/s2 (b)
the behaviour of Y 2 is known [20], we unveil an interesting scaling behaviour of Y 1 beyond the large-time regime. Indeed, in the limit t → ∞ with t/L 2 constant, one finds

Y 2 (t) ρ 0 = s 2 2t π , Y 1 (t) ρ 0 = s 2 2t π g L √ 2t , (4)
g(u) = e -u 2 - √ π u erfc(u). ( 5
)
This provides the dynamics of the entrainment of the TP1 by TP2, which admits a typical time scale L 2 and leads to the final bound state discussed above (see Fig. 2). The evolution towards the final regime can further be quantified by the dynamics of the two-TPs even (resp. odd) cumulants In the case (ii), the velocity of the first TP changes sign at rescaled time τ * = t * /L 2 (gray square in ii.a). We plotted the prediction τ * as a function of s1/s2 for s2 = 0.8. The dashed line are the asymptotic behaviors from Eq. ( 12), the grey square corresponds to the one in (ii.a).

κ e = Y p 1 Y q 2 c (resp. κ o = Y p 1 Y q 2 c ) L L 0.1 1 10 100 1000 0.4 0.6 0.8 1 1.2 A A 1 A 2 (i.a) 0 5 10 0 2.5 Y /(ρ 0 L) 0.1 1 10 100 1000 0.01 0.1 1 10 τ = t/L 2 • c /(ρ 0 L) Y 2 1 c Y 2 2 c Y 1 Y 2 c Y 1 Y 2 2 c /S (i.b) 0.1 1
with p + q even (resp. odd), p, q ≥ 1. They obey

κ e ρ 0 = κ o ρ 0 s 2 = 2t π g L √ 2t . (6)
Several comments are in order. (i) Equations ( 4) and ( 6) are similar to the expressions found for the random average process [38,39], which points towards their universality. (ii) The same scaling function g is involved in the expressions of Y 1 (t) and Y 1 Y 2 (t) c (Eqs ( 4) and ( 6)). This leads to the generalized fluctuation-dissipation relation

lim f2→0 2 β Y 1 (f 1 = 0, f 2 ) f 2 = Y 1 Y 2 c (f 1 = f 2 = 0). ( 7
)
Note that this relation holds in the opposite limit of a dilute (ρ → 0) SEP [40]. (iii) Our approach provides the time dependence of all cumulants of individual particles and the law of the distance between TPs (insets of Fig. 2 and [35]). The time, initial distance between TPs, and driving force dependences from numerical simulations are unambiguously captured by our theoretical expressions (Fig. 2). Bath-mediated cooperatively and competition. We now turn to the general case in which both TPs are biased (see Figure 3, top). The dynamics of effective interactions between TPs at the level of averages is conveniently analysed by introducing the rescaled instantaneous velocities

A j (t) = √ 2πt ρ 0 d Y j dt , (8)
which satisfy A j (t) = s j at small time and

A j (t) = S = s 1 + s 2 1 + s 1 s 2 at large time. The full time dependency is found [35] to be given by A 1 = H s1,s2(1+s1),-s1s2 (L/ √ 2t) and A 2 = H s2,s1(1-s2),s1s2 (L/ √ 2t) with H β0,β1,β2 (u) = ∞ n=0 (-s 1 s 2 ) n 2 m=0 β m e -[(2n+m)u] 2 (9)
while higher order cumulants follow

κ e ρ 0 = κ o ρ 0 S = 2t π G s1s2 L √ 2t , (10)
G σ (u) = (1 + σ) ∞ n=0 (-σ) n g([2n + 1]u). ( 11
)
These results fully quantify the dynamics of the BMIs between two biased TPs and reveal striking behaviors. (i) In the case of same sign biases, the TPs cooperate [START_REF] Bénichou | [END_REF] (Fig 3 left). At large times a pair of biased TPs moves faster than a single TP-in agreement with the asymptotic result Eq. (3). Note that such an accelerated dynamics has been numerically observed in two-dimensional systems [22,24]. At intermediate times, we unveil an overshoot of the rescaled velocity of the trailing TP. (ii) When the biases act in opposite directions (say 0 < -s 1 < s 2 ), each TP starts to move in the direction of its own force, and eventually both TPs move in the direction of the largest force (Figure 3 left). This competing stage can be quantified from Eq. ( 9) by determining the U-turn time t * at which the velocity of TP1 changes its sign. This time t * vanishes when s 1 is small, and diverges when s 1 is close to -s 2 according to the scaling laws

t * L 2 ∼ s1→0 1 2 log(-s 2 /s 1 ) ; t * L 2 ∼ s1→-s2 γ 1 + s 1 /s 2 (12)
with γ = 2(1 + s 2 ) 2 /(1s 2 ). Figure 3 shows an excellent quantitative agreement between the analytical predictions and the numerical simulations.

Our approach can be extended to determine the dynamics of correlations in the case of an arbitrary number of driven TPs. Cooperativity and competition involve a complex cascade of time scales associated with the initial distances between TPs, fully captured by our approach, as exemplified in the case of 3 TPs in Figure 4 and [35].

Bath-mediated interactions. The BMIs between two bound biased TPs can be further analysed by associating the probability distribution of the variation of distance D = Y 2 -Y 1 to an effective potential U (D) via P (D) = exp(-βU (D)) [41]. For identical forces (f 1 = f 2 = f ) and D sufficiently close to its average value, the two TPs are effectively bound by an harmonic potential

U (D) ∼ κ 2 (D -D ) 2 , D ρ 0 L = 1 cosh (βf ) -1 , ( 13
)
where the constant κ is explicitly given by

κ = cosh (βf ) βρ 0 L (1 + cosh (βf )) . (14)
In the regime, D D , the potential displays a weaker dependence on the distance U (D) ∼ D(ln D + ν)/β (see [35] for the value of ν). We remark that this qualitative change of regimes has been observed in 2d in numerical simulations of two biased TPs in a quiescent colloidal bath [24]. 4: Rescaled velocity of three TPs for biases s1 = 0.9, s2 = 0.3, s3 = -0.6, total distance L = X 0 3 -X 0 1 = 60 with X 0 2 -X 0 1 = 45, ρ0 = 10 -2 . The colored circles correspond to the numerical simulations while the dashed lines are the theoretical predictions [35]. Inset: average displacements of the TPs in linear rescaled time. The behavior of the second TP displays several regimes: it first moves to the right, then under the influence of the third TP it goes to the left and finally the first TP pushed it back to the right. Note that this complex dynamics is captured by our theoretical approach.

L 0.001 0.01 0.1 1 10 100 -1 -0.5 0 0.5 1 τ = t/L 2 A TP1 TP2 TP3 0 0.5 -0.2 0 0.2 FIG.
Altogether, we determined the full dynamics of correlation functions in a paradigmatic model of non equilibrium statistical physics and entirely characterized the corresponding bath-mediated interactions. Our model is based on the well known Symetric Exclusion Principle (SEP). Let us consider particles on a discrete one-dimensional line (the sites correspond to the integers). The mean density is denoted ρ ∈ [0, 1]. The particles follow a symmetric random walk with hard-core exclusions in continuous time. The rate of jumps of these particles to the left / right is 1/2. In addition to these particles, we consider two biased tagged particles (TPs), initially at positions X 0 1 = 0 and X 0 2 = L. They can jump to the right with probability p j = (1 + s j )/2 and to the left with probability p -j = (1s j )/2, where s 1 and s 2 are the biases on the two TPs. See Fig. 1 on the article for a sketch of the system.

For simplicity reasons, we choose to present only the computations for two TPs. The treatement of three or more TPs is very similar but it leads to heavier computations. Results for three TPs will be presented in the next section.

b. Approximation and consequences

We focus on the limit of a dense system (ρ → 1) and we follow the evolution of the vacancies as was done in Ref. [20]. In this limit of small number of vacancies, one can approximate the motion of the TPs as being generated by the vacancies interacting independantly (i.e. not simultaneously) with them: in the large density limit the events corresponding to two vacancies interacting simultaneously with the TPs happen only with negligible probability.

Let us consider a system of size N with M vacancies and denote by

Y(t) = (X 1 (t) -X 0 1 , X 2 (t) -X 0
2 ) the vector of the displacements of the two tracers. The probability P (t) (Y|{Z j }) of having displacements Y at time t knowning that the M vacancies started at sites Z 1 . . . Z M is exactly given by :

P (t) (Y|{Z j }) = Y1,...,Y M δ Y,Y1+•••+Y M P (t) ({Y j }|{Z j }) (S1)
where P (t) ({Y j }|{Z j }) is the probability of displacements Y j due to the vacancy j, for all j, knowing the initial positions of all the vacancies. Making the approximation that we described above, we can link it to the probability p

(t)
Z (Y) that the tracers have moved by Y at time t due to a single vacancy that was initially at site Z:

P (t) ({Y j }|{Z j }) ∼ ρ→1 M j=1 p (t) Zj (Y j ) (S2)
We then have :

P (t) (Y|{Z j }) ∼ ρ→1 Y1,...,Y M δ Y,Y1+•••+Y M M j=1 p (t) Zj (Y j ) (S3)
We take the Fourier transform and we average over the initial positions of the vacancies:

p(t) (k) ≡ 1 N -2 Z =0,L Y p (t) Z (Y)e ik•Y (S4)
and mutatis mutandis for P (t) (k). We obtain

P (t) (k) ∼ ρ→1 p(t) (k) M =   1 N -2 Z =0,L p(t) Z (k)   M . ( S5
)
Furthermore we write p(t) Z (k) = 1 + q(t) Z (k) (q Z corresponds to the deviation from a Dirac centered in 0) to get

P (t) (k) ∼ ρ→1   1 + 1 N -2 Z =0,L q(t) Z (k)   M . ( S6
)
We now take the limit N, M → ∞ with ρ 0 ≡ 1ρ = M/N (density of vacancies) remaining constant. The second characteristic function reads

ψ (t) (k) ≡ ln P (t) (k) ≈ ρ 0 Z =0,L q(t) Z (k). (S7)
ψ gives all the cumulants : the cumulant associated to

Y q 1 Y r 2 is Y q 1 Y r 2 = i -(q+r) ∂ q+r ∂k q 1 k r 2 ψ k=0 .
Our goal is now to find an expresion for ψ (t) (k).

Resolution

The approximation (S3) that leads to (S6) tells us that the case of a vanishing density of vacancies can be deduced from the case of a single vacancy. We first focus on the latter before using (S7) to get the cumulants.

a. Single vacancy

We consider the situation with a single vacancy. In the following we shall denote by ν = -1, +1, -2, +2 the "special sites" respectively to the left of TP1, to the right of TP1, to the left of TP2, to the right of TP2. We use the vectors e ±1 = (±1, 0) and e ±2 = (0, ±1).

The key is to introduce (conditional) first-passage time quantities

F (t)
η,A denotes the probability that the vacancy that started from site A (A can be a special site) at time 0 arrives for the first time to the position of one of the tracers at time t conditionned on the fact that it was on the special site η at time t -1.

A subtelty is that this probability shall depend on a quantity called L * = L * η,Z which is the distance between the tracers when a vacancy can arrive at site η when it was at site Z at the very begining. Let us denote z = I(0 < Z < L), Z being the inital position of the vacancy. Example: we consider η = +1, the vacancy has to be between the two tracers to get to this special site before any other. Now if z = 0 (vacancy started outside of the tracers) the distance between the tracers is L * = L + 1 (if z = 1, it would be L * = L). In any case, we always have L * = L + α(z, η), α ∈ {-1, 0, 1}. In the following F (j),z η,A ≡ F (j) η,A (L + α(η, z)). Note that if Z is the initial position of the vacancy, there is no ambiguity:

F (j) η,Z = F (j) η,Z (L * = L).
One can partition over the first passage of the vacancy to the site of one of the tracers to get an expression for qZ which is the main quantity involved in (S7):

p (t) Z (Y) = δ Y,0   1 - t j=0 ν=±1,±2 F (j) ν,Z   + t j=0 ν=±1,±2 p (t-j) -ν (Y -e ν )F (j) ν,Z (S8) p(t) Z (k) = 1 - t j=0 ν=±1,±2 F (j) ν,Z + t j=0 ν=±1,±2 p(t-j),z -ν (k)e ik•eν F (j) ν,Z , (S9) q(t) Z (k) = - t j=0 ν=±1,±2 1 -1 + q(t-j),z -ν (k) e ik•eν F (j) ν,Z . (S10)
An exponant z to a quantity means that this quantity is computed taking into account z = I(0 < Z < L).

We now need an expression for qz η where η is a special site. To do so we decompose the propagator of the displacements over the successive passages of the vacancy to the position of one of the tracers:

p (t),z η (Y) = δ Y,0   1 - t j=0 µ F (j),z µ,η   + ∞ p=1 ∞ m1,...,mp=1 ∞ mp+1=0 δ t, i mi ν1,...,νp δ Y, i eν i   1 - mp+1 j=0 µ F (j),z µ,-νp   F (mp),z νp,-νp-1 . . . F (m2),z ν2,-ν1 F (m1),z ν1,η (S11)
the sums on µ and ν i run over the special sites (±1, ±2).

The discrete Laplace transform of a function of time g(t) is ĝ(ξ) ≡ ∞ t=0 g(t)ξ t . We can now take both the Laplace and Fourier transforms of (S11) to get:

pz η (Y, ξ) = 1 1 -ξ    δ Y,0 1 - µ F z µ,ν + ∞ p=1 ν1,...,νp δ Y, i eν i µ 1 -F z µ,-νp F z νp,-νp-1 . . . F z ν2,-ν1 F z ν1,η    (S12) qz η (k, ξ) ≡ pz η (k, ξ) - 1 1 -ξ = 1 1 -ξ µ,ν {[1 -T z (k, ξ)] -1 } νµ × 1 -e -i keν e i keµ F z µη (ξ) (S13)
The matrix T is defined by T z (k, ξ) νµ = F z ν,-µ (ξ)e i keν .

b. Multiple vacancies

Introducing (S10) into (S7), one gets an expression for the Laplace transform of the second characteristic function:

ψ(k, ξ) = -ρ 0 ν=-1,2 1 1 -ξ - 1 1 -ξ + q (z=0) -ν (k, ξ) e i keν h ν (ξ) + ν=1,2 1 1 -ξ - 1 1 -ξ + q (z=1) -ν (k, ξ) e i keν h ν (ξ) (S14) h ν (ξ) = Z =0,L Fν,Z (ξ) (S15)
where qz ν is given by (S13). At the end of the day, the only quantities that need to be computed are : F-1,-1 , F2,2 (returning to the left of TP1 / to the right of TP2), F1,1 , F-2,-2 (returning to the left of TP1 without touching TP2, and vice-versa), F 1,-2 , F -2,1 (arriving to the right of TP2 starting from the left of TP1 without returning, and vice-versa), and the sums h

-1 , h 2 , h 1 , h -2 .
a. Computation of the "outside" quantities We call "outside" quantities the ones related to the sites -1 and 2: F-1,-1 (ξ), F2,2 (ξ), h 1 (ξ) and h 2 (ξ). In this case, we do not care about the fact that there are two tracers: we simply consider a random walk on a half line in which site 1 is biased. Note that in this case, we do not care either about the distance between the tracers (ie the value of z).

We first compute F 2,2 : this amounts to the probability to return for the first time to the origin at time t and being on site 1 at time t -1 knowing that one was at site 1 at time 0. We can partition on the first step of the walk.

For η a "special site", we call pη = 1/(2p η + 1) the probability for the vacancy to jump to the left (ν = 1) or right (ν = -1) of a tracer when it is next to it.

F 2,2 (t) = (1 -p2 ) δ t,1 + p2 t j=1 f 1 (j -1)F 2,2 (t -j) (S16)
where f l (t) the first-passage time density at the origin at time t of a symmetric Polya walk starting from site l. One knowns fl (ξ) = α |l| with α = ξ -1 (1 -1ξ 2 ) [START_REF] Hughes | Random walks and random environments[END_REF].

We take the Laplace transform and obtain:

F2,2 (ξ) = (1 -p2 )ξ 1 -p2 αξ (S17) Similarly F-1,-1 (ξ) = (1 -p-1 )ξ 1 -p-1 αξ (S18)
h -1 and h 2 are easily computed from the former quantities:

F +2,Z (t) = Θ(Z -(L + 1)) t j=0 f Z-1 (j)F 2,2 (t -j) (S19) F+2,Z (ξ) = Θ(Z -(L + 1)) fZ-L-1 (ξ) F2,2 (ξ) = Θ(Z -(L + 1)) F2,2 (ξ)α Z-L-1 (S20) h 2 (ξ) = Z =0,L F+2,Z (ξ) = F2,2 ∞ Z =0 α Z = F2,2 (ξ) 1 -α (S21) h -1 (ξ) = F-1,-1 (ξ) 1 -α (S22)
b. Computation of the "inside" quantities We first solve the case with no bias s 1 = s 2 = 0. There are only two quantity to compute: F (=) ≡ F+1,+1 = F-2,-2 and F ( =) ≡ F+1,-2 = F-2,+1 . We also compute (for 0 < Z < L)

F • Z ≡ F+1,Z = F-2,L-Z .
(L is the distance between the tracers) We recall the formula for the first passage time (of a Polya random walk) at site s 1 starting from s 0 and considering s 2 as an absorbing site [START_REF] Hughes | Random walks and random environments[END_REF],

f † (s 1 |s 0 , ξ) = fs1-s0 (ξ) -fs1-s2 (ξ) fs2-s0 (ξ) 1 -fs1-s2 (ξ) 2 (S23)
with fl (ξ) = α |l| in our 1-dimensional case. This gives us:

F (=) (ξ, L) = α -α 2L-1 1 -α 2L (S24) F ( =) (ξ, L) = α L-1 -α L+1 1 -α 2L (S25) F • Z (ξ, L) = α Z -α 2L-Z 1 -α 2L (S26)
The case s 1 = 0, s 2 = 0 can be deduced from the unbiased case for a distance L -2 by partioning on the first passage of the vacancy to either +1 or -2:

F 1,1 (t, L) =(1 -p1 )δ t,1 + p1   t j=1 F (=) (j -1, L -2)F 1,1 (t -j, L) + t j=1 F ( =) (j -1, L -2)F 1,-2 (t -j, L)   (S27) F1,1 (ξ, L) =(1 -p1 )ξ + p1 ξ F (=) (L -2) F1,-2 (L) + F ( =) (L -2) F1,-2 (L) (S28)
The first term corresponds to a return to 1, the second to a passage to -2. This leads us to a system of 2 × 2 equations with 2 × 2 unknowns.

           ξ p1 F (=) (L -2) -1 F1,1 + ξ p1 F ( =) (L -2) F1,-2 +(1 -p1 )ξ = 0 ξ p-2 F (=) (L -2) -1) F1,-2 + ξ p-2 F ( =) (L -2) F1,1 = 0 ξ p-2 F ( =) (L -2) -1 F-2,-2 + ξ p-2 F ( =) (L -2) F-2,1 +(1 -p-2 )ξ = 0 ξ p1 F (=) (L -2) -1) F-2,1 + ξ p1 F ( =) (L -2) F-2,-2 = 0 (S29)
The solution is:

                           F1,1 = (1 -ξ p-2 F (=) )(1 -p1 )ξ (1 -ξ p1 F (=) )(1 -ξ p-2 F (=) ) -ξ 2 p1 p-2 ( F ( =) ) 2 F1,-2 = ξ 2 p-2 (1 -p1 ) F ( =) (1 -ξ p1 F (=) )(1 -ξ p-2 F (=) ) -ξ 2 p1 p-2 ( F ( =) ) 2 F-2,-2 = (1 -ξ p1 F (=) )(1 -q -B )ξ (1 -ξ p1 F (=) )(1 -ξ p-2 F (=) ) -ξ 2 p1 p-2 ( F ( =) ) 2 F-2,1 = ξ 2 p1 (1 -p-2 ) F ( =) (1 -ξ p1 F (=) )(1 -ξ p-2 F (=) ) -ξ 2 p1 p-2 ( F ( =) ) 2 (S30)
with F (=) and F ( =) evaluted for L -2. For 0 < Z < L, one can again partion on the first passage of the vacancy to +1 or -2.

F +1,Z (t, L) = t j=1 F • Z-1 (j, L -2)F 1,1 (t -j, L) + t j=1 F • L-Z-1 (j, L -2)F -2,-2 (t -j, L) (S31) F+1,Z (ξ, L) = F • Z-1 (L -2) F1,1 (L) + F • L-Z-1 (L -2) F-2,-2 (L) (S32)
Using (S26), we obtain:

h 1 (ξ) = L-1 Z=1 F+1,Z = (1 -α L-2 )(1 -α L-1 ) 1 -α 2(L-2) F+1,+1 + F+1,-2 (S33) h -2 (ξ) = L-1 Z=1 F-2,Z = (1 -α L-2 )(1 -α L-1 ) 1 -α 2(L-2) F-2,-2 + F-2,+1 (S34) 3. Results

a. Symbolic calculation

We use a symbolic calculation software (Mathematica) to put all the bricks together according to the following graph. An asterix ( * ) indicates that one should be careful about which value of L (distance between the TPs) should be used.

ψ(k, ξ) (S14) qz ν (S13) hν (S21, S22, S33, S34) F-1,-1, F2,2 (S17, S18) Fα,β (α, β ∈ {1, -2}) (S30) * F (=) , F = (S24, S25) * α = 1- √ 1-ξ 2 ξ pν = 1 2pν +1 b. Scaling limit
We define a continuous Laplace variable p = 1ξ. When p goes to zero, α = (1 -1ξ 2 )/ξ = 1 -√ 2p + O(p). We will be interested in the limit of large time t and large distance L, with τ = t/L 2 constant. This amounts to writing p = p/L 2 with L → ∞ and p constant. The inverse Laplace transform of a function f (p, L) = L 2 ĝ(p) is f (t, L) = g(τ ) with g the inverse Laplace transform of ĝ.

Finally, our scaling limit in Laplace space amounts to the following substitutions that can be performed with a numerical software (in this order).

1 -ξ → p/L 2 1 -α → 2p/L α L ≈ 1 -2p/L L → e - √ 2 p α → 1 (S35)
Appendix B: Detailed results

One TP

We start by recalling the results of Ref. [20] for a single TP of bias s in the limit of high density (fraction of vacancies ρ 0 → 0), at large time (small p).

lim ρ0→0 ψ(k, p) ρ 0 = 1 √ 2p 3/2 {cos k -1 + s sin k} (S1) lim ρ0→0 ψ(k, t) ρ 0 = 2t π {cos k -1 + s sin k} (S2)
The even (resp odd) cumulants of the displacement Y are equal.

Y 2n c ρ 0 = 2t π Y 2n+1 c ρ 0 = s 2t π (S3)
The fraction • is understood as a limit ρ 0 → 0.

Two TPs

We now focus on the case of two TPs with biases s 1 and s 2 , initially separated by a distance L. The Laplace transform of the cumulant generating function reads

lim ρ0→0 ψ(k, p = p/L 2 ) ρ 0 = L 3 √ 2p 3/2 Ke,2 (p)(cos(k 1 + k 2 ) -1) + Ko,2 (p) sin(k 1 + k 2 ) + 2 i=1 Ke,1 i (p)(cos k i -1) + Ko,1 i (p) sin k i . (S4)
with Kα,j given in the following. In real time, this gives

1 L lim ρ0→0 ψ(k, t = L 2 τ ) ρ 0 = 2τ π K e,2 (τ )(cos(k 1 + k 2 ) -1) + K o,2 (τ ) sin(k 1 + k 2 ) + 2 i=1 K e,1 i (τ )(cos k i -1) + K o,1 i (τ ) sin k i . (S5)
Note that K(τ ) is related but not equal to the inverse Laplace transform of K(p).

From the previous expression one can deduce all the cumulants (i = 1, 2), including those of the variation of distance

D = Y 2 -Y 1 . Y j 1 Y 2n-j 2 c ρ 0 L = K e,2 (τ ) 2τ π Y j 1 Y 2n+1-j 2 c ρ 0 L = K o,2 (τ ) 2τ π (S6) Y 2n i c ρ 0 L = [K e,2 (τ ) + K e,1 i (τ )] 2τ π Y 2n+1 i c ρ 0 L = [K o,2 (τ ) + K o,1 i (τ )] 2τ π (S7) D 2n c ρ 0 L = [K e,1 1 (τ ) + K e,1 2 (τ )] 2τ π D 2n+1 c ρ 0 L = [K o,1 2 (τ ) -K o,1 1 (τ )] 2τ π (S8)
Again, these results correspond to a limit ρ 0 → 0. We now give the expressions of the quantities K(p) and K(t) both in the cases s 1 = 0, s 2 = 0 and s 1 , s 2 = 0.

Case

s 1 = 0, s 2 = 0 Case s 1 , s 2 = 0 Ke,2 (p) v (1 + s 1 s 2 )v/d 2 Ko,2 (p) s 2 v (s 1 + s 2 )v/d 2 Ke,1 1 (p) (1 -v)(1 + s 2 v) (1 -v)(1 + s 2 v)/d 2 Ko,1 1 (p) 0 s 1 (1 -v)(1 + s 2 v)/d 2 Ke,1 2 (p) 1 -v (1 -v)(1 -s 1 v)/d 2 Ko,1 2 (p) s 2 (1 -v) s 2 (1 -v)(1 -s 1 v)/d 2 with v = e - √ 2 p and d 2 = 1 + s 1 s 2 v 2 .
s 2 . The expansion gives

A 1 = ∞ n=0 (-s 1 s 2 ) n s 1 (e -1/(2τ ) ) (2n) 2 + s 2 (1 + s 1 )(e -1/(2τ ) ) (2n+1) 2 -s 1 s 2 (e -1/(2τ ) ) (2n+2) 2 (S18) = ∞ n=0 s 2n 2 (1 -n )× -s 2 (1 -)(1 -(2n) 2 η ) + s 2 (1 -s 2 + s 2 )(1 -(2n + 1) 2 η ) + s 2 2 (1 -)(1 -(2n + 2) 2 η ) (S19) = 0 + ∞ n=0 s 2n 2 s 2 (1 + (2n) 2 η) + s 2 2 -s 2 (1 -s 2 )(2n + 1) 2 η -s 2 2 -s 2 2 (2n + 2) 2 η (S20) = ∞ n=0 s 2n+1 2 + η ∞ n=0 s 2n 2 s 2 (2n) 2 -s 2 (1 -s 2 )(2n + 1) 2 -s 2 2 (2n + 2) 2 (S21) = s 2 1 -s 2 2 -η s 2 (1 + s 2 ) (1 -s 2 ) 2 (S22) A 1 = 0 gives η = s 2 1 -s 2 2 (1 -s 2 ) 2 s 2 (1 + s 2 ) = 1 -s 2 (1 + s 2 ) 2 . (S23)
At the end of the day,

τ * ≈ 1 2η = (1 + s 2 ) 2 2(1 -s 2 ) 1 1 + s1 s2 . (S24)
Eqs (S17) and (S24) correspond to the asymptots of the inset of Fig. 3, (b.ii) of the main text.

Single trajectories

We note that for large time and large distances, our results for cooperativity and competition (Eqs. ( 4) and ( 9) of the main text) actually hold at the very level of a single trajectory. Numerical evidence is shown on Fig. S1).

Three TPs

We now turn to the case of three TPs with biases s 1 , s 2 , s 3 . We denote L 1 = X 0 2 -X 0 1 and L 2 = X 0 3 -X 0 2 and L = L 1 + L 2 . Our result for the Laplace transform of the cumulant-generating function is 4) and ( 9) in the main text. We see that our approach is relevant even for individual realizations.

lim ρ0→0 ψ(k, p = p/L 2 ) ρ 0 = L 3 √ 2p 3/2 Ke,3 (p)(cos(k 1 + k 2 + k 3 ) -1) + Ko,3 (p) sin(k 1 + k 2 + k 3 ) + 2 i=1 Ke,2 i,i+1 (p)(cos(k i + k i+1 ) -1) + Ko,1 i (p) sin(k i + k i+1 ) + 3 i=1 Ke,1 i (p)(cos k i -1) + Ko,1 i (p) sin k i . (S25) with K α,n (p) = 2 a,b=0 Q α,n (a, b)v a 1 v b 2 1 + s 1 s 2 v 2 1 + s 2 s 3 v 2 2 + s 1 s 3 v 2 1 v 2 2 (S26)
with v 1 = e -(L1/L) √ 2 p, v 2 = e -(L2/L) √ 2 p. For completeness, we give the 12 matrices Q of coefficients.

Q e,3 = (1 + s1s2 + s2s3 + s1s3)   0 0 0 0 1 0 0 0 0   Q o,3 = (s1 + s2 + s3 + s1s2s3)   0 0 0 0 1 0 0 0 0   (S27) Q e,2 1,2 = (1 + s1s2)   0 0 0 1 s3 -1 -s3 0 0 0   Q o,2 1,2 = s1 + s2 1 + s1s2 Q 2,e 1,2 (S28)
Q e,2 2,3 = (1 + s2s3)   0 1 0 0 -1 -s1 0 0 s1 0   Q o,2 2,3 = s2 + s3 1 + s2s3 Q 2,e 2,3 (S29)
Q e,1 1 =   1 0 s2s3 s2 -1 0 s3(1 -s2) -s2 0 -s3   Q o,1 1 = s1Q 1,e 1 (S30) Q e,1 2 =   1 s3 -1 -s3 -1 -s1 1 + s1 -s3 -s1s3 s3(1 + s1) s1 s1(s3 -1) -s1s3   Q o,1 2 = s2Q 1,e 2 (S31) Q e,1 3 =   1 -1 -s2 s2 0 0 0 s1s2 -s1(1 + s2) s1   Q o,1 3 = s3Q 1,e 3 (S32)
The Laplace transform (S25) can be inverted numerically to obtain the time evolution of the cumulants. One notices that the sum of all the coefficients of the matrices Q 1 and Q 2 is zero. This means that for p → 0, (S25) simplifies into lim We consider three TPs undergoing forces f 1 , f 2 and f 3 . We denote L 1 = X 2 (t = 0) -X 1 (t = 0) the initial distance between TPs 1 and 2. Let us call D 12 = Y 2 -Y 1 the variation of distance between TPs 1 and 2. From (S25), we show that its cumulant generating function at large time reads ψ D12 (q) ≡ ln e iqD12 (S58) ψ stat D12 (q) ρ 0 , = L 1 2 e βf1 + e -β(f2+f3) (e iq -1) + (e -iq -1) .

(S59)

This still corresponds to a Skellam law. The presence of a force f 3 modifies the stationnary law. We see that the interactions between TPs are not additive.

Similarly, if we denote L 2 = X 3 (t = 0) -X 2 (t = 0) and D 23 = Y 3 -Y 2 , we have ψ stat D23 (q) ρ 0 = L 2 2 e β(f1+f2) + e -βf3 (e iq -1) + (e -iq -1) .

(S60)

FIG. 2 :

 2 FIG. 2: Entrainment of the bath particles by a biased TP for ρ0 = 10 -2 . Only the right TP is biased. (a) Evolution of the mean displacements, different symbols corresponding to L = 10, 50, s2 = -0.2, 0.8. The black lines are the predictions from Eq. (4). (b) Cumulants Y1Y2 c and Y1Y 2 2 c/s2 for the same set of parameters. The black line corresponds to Eq. (6). Inset of (a): variances for L = 10, s2 = 0.8, -0.8 (solid, dashes). Predictions in black [35]. Inset of (b): law of the variation of distance D = Y2 -Y1 at times 10, 10 2 , 10 3 , 10 4 (blue to red) for L = 10, s = 0.8, ρ0 = 0.05. The squares are the numerical results, the colored lines are the theoretical predictions, the black line is the asymptotic prediction [35].

2 FIG. 3 :

 23 FIG. 3: Cooperativity and competition (ρ0 = 10 -2 ). (i) The two TPs have identical biases s1 = s2 = 0.8. (ii) The biases are in opposed directions s1 = -0.6, s2 = 0.8. The rescaled velocities A are plotted on (i.a) and (ii.a) with the displacements in the inset, two initial distances are plotted: L = 50 and 200. At short time, the rescaled velocities are s1 and s2, at large time they are both equal to S. The variances and some cumulants are plotted on (i.b) and (ii.b) (initial distances L = 10 and 20).In the case (ii), the velocity of the first TP changes sign at rescaled time τ * = t * /L 2 (gray square in ii.a). We plotted the prediction τ * as a function of s1/s2 for s2 = 0.8. The dashed line are the asymptotic behaviors from Eq. (12), the grey square corresponds to the one in (ii.a).

  FIG. S1: Single realizations of the TP trajectories (ρ0 = 0.1). (a) Entrainment: s1 = 0, s2 = 0.8. (b) Competition: s1 = -0.4, s2 = 0.8. Circles: L = 1000, triangles: L = 2000, black lines: theoretical predictions from Eqs. (4) and (9) in the main text. We see that our approach is relevant even for individual realizations.

3 1 Ix 2

 312 cos(k 1 + k 2 + k 3 ) -1 + s 1 + s 2 + s 3 + s 1 s 2 s 3 1 + s 1 s 2 + s 1 s 3 + s 2 s 3 sin(k 1 + k 2 + k 3 ) (S33)with the effective potential U (D) given by 1 + µ 2 ) = cosh(βf )βρ 0 L[1 + cosh(βf )] . (S46) c. Large deviation function and large distance scalingNow, we want to derive the behavior of U (D) at large D. We first determine the large deviation function of the problem. We haveP stat D (D) = e -ρ0L(μ1+1) μD/2 D (2ρ 0 L μ1 ), (S47) ln P stat D (D) = -ρ 0 L(μ 1 + 1) + D 2 ln μ1 + ln I D (2ρ 0 L μ1 ) (S48) with μ1 = µ 1 /(ρ 0 L) = cosh(βf ) -1 .We use the expansion of the Bessel function at large order and large argument (here large D ∼ ρ 0 L).I D (Dz) ≈ 1 √ 2πD e ηD (1 + z 2 ) 1/4 η = 1 + z 2 + ln z 1 + √ 1 + z 2 (S49) ln I D (Dz) ≈ ηD ln I 2ρ0Lx (2ρ 0 Lxz) ≈ 2ρ 0 Lxη (+ μ1 + x ln μ1x ln(x + x 2 + μ1 ) (S53) ≡ φ(x)(S54) φ(x) is a large deviation function. The large D behavior corresponds to the limit x → ∞. φ(x) ≈ -x ln x + (1ln 2 + ln μ1 )x, (S55) ln P stat D (D) ≈ D 1ln D ρ 0 Lμ 1 ≈ D 1ln D cosh(βf )

  d. Three TPs and non-additivity of the interactions

Case s 1 = 0, s 2 = 0 Case s 1 , s 2 = 0 K e,2 (τ ) g 1 (τ ) G 0,s1s2,0 (τ ) K o,2 (τ ) s 2 g 1 (τ ) G 0,s1+s2,0 (τ ) K e,1 1 (τ ) 1 + (s 2 -1)g 1 (τ ) + s 2 g 2 (τ )

with

The time evolution of the cumulants is given by Eqs. (S6), (S7). The rescaled velocities studied in the main text are obtained by taking the time derivative of the average displacements.

U-turn time for two TPs with opposite biases

We consider two TPs with biases s 2 > 0 and -s 2 < s 1 < 0. We focus on the rescaled velocity of particle 1 defined as.

From (S7),

From this we deduce that

There exists a rescaled time τ * such that A 1 (τ * ) = 0. One can solve numerically for τ * . Here we derive the asymptotic behavior of τ * when s 1 /s 2 → 0 and when s 1 /s 2 → -1.

When s 1 /s 2 is small, τ * is small so (τ * ) -1/2 is large. We can keep only the first two terms in the sum (S15),

.

(S17)

In the opposite limit, we write s 1 = -s 2 (1 -) with 1. We guess that (2τ * ) -1 = η with η depending only on This corresponds to a single TP with effective bias S = s1+s2+s3+s1s2s3 1+s1s2+s1s3+s2s3 . This behavior generalises to an arbitrary number of TPs, the effective bias being

with s i = tanh(f i /2).

6. Law of the distance between two TPs a. Two TPs

We denote D = Y 2 -Y 1 the variation of distance between two TPs. The cumulant-generating function of D is given by

with ψ given by Eq. (S5). Thus the Laplace transform of ψ D reads:

Ke,1

. This lead us to the following stationnary law ψ stat D at large time:

= L 2 e βf1 + e -βf2 (e iq -1) + (e -iq -1) .

(S41)

We recall that the forces are defined by e βfi = (1 + s i )/(1s i ). Eq. (S41) corresponds to a Skellam distribution P D with parameters µ 1 = 2ρ 0 L/(e βf1 + e -βf2 ) and µ 2 = ρ 0 L.

The average of D is µ 1µ 2 and the variance is

b. Effective potential and Gaussian limit

We focus on the case

In this case, we can associate an effective potential U (D) to the probability law as

For D close to its average value µ 1µ 2 , the Skellam law is well approximated by a Gaussian

We do simulations of the SEP as follow. We consider a periodic line of N sites. We first place n = 1, 2 or 3 tagged particles (TPs) at initial deterministic positions. Then, we place Mn bath particles uniformly at random on the remaining sites. For the figures of the main text we used N = 2500, M = 2475 (ρ 0 = M/N = 0.01).

At each time step we draw a particle at random and try to move it either to the left or to the right according to its probability p = (1 + s)/2 for the TPs, p = 1/2 for the bath particles. The continuous time is then incremented with a exponentially distributed random number with characteristic time 1/M . This accounts for the exponentially distributed jump rates of the particles. At equally spaced continuous times we compute the observables: the moments with 1 to n TPs.

We perform the average over multiple simulations to get good statistics. Typically we perform 50000 simulations (except for the Figure for a single realization).

Numerical effective biases

We now discuss how to obtain the numerical curves for the rescaled velocities defined as Our procedure is as follow: we first gather the data into bins of logarithmic size over which we average. This gives us B(w). We then compute the derivative dB/dw using a Savisky-Golay filter. We finally obtain A(w) that we can plot on our graphs.


Radulesco Thomas
email: thomas.radulesco@ap-hm.fr

Meister Lionel

Bouchet Gilles

Giordano Jérome

Dessi Patrick

Pierre Perrier

Michel Justin

Thomas Radulesco

Computational fluid dynamics and septal deviations-Virtual surgery can predict post-surgery results: A preliminary study including two patients

Keywords: CFD, nasal airway obstruction, virtual surgery, unsteady flow

Nasal airway obstruction is a public health problem resulting in discomfort, manifested as a sensation of insufficient airflow through the nose [START_REF] Osborn | Chapter 2: Nasal obstruction[END_REF]. Septal deviation (SD) is one of the most common causes of nasal obstruction and SD's standard treatment is surgery [START_REF] Teixeira | Nasal Septal Deviations: A Systematic Review of Classification Systems[END_REF].

CTscanner and rhinomanometry are currently most used examinations in the exploration of SD [START_REF] Vogt | Four-phase rhinomanometry: a multicentric retrospective analysis of 36,563 clinical measurements[END_REF]. Computational Fluid Dynamics (CFD) is a tool bringing new perspectives, allowing to model anatomy and airflows of sinonasal tract [START_REF] Leong | A review of the implications of computational fluid dynamic studies on nasal airflow and physiology[END_REF]. CFD aims to study airflow and air conditioning and provides better understanding of nasal physiology.

Many CFD studies have been conducted in recent years [START_REF] Leite | The clinical implications of computerised fluid dynamic modelling in rhinology[END_REF]. Correlations has been established between patient perception and CFD variables i.e. heat flux [START_REF] Radulesco | Correlations between Computational Fluid Dynamics and clinical evaluation of nasal airway obstruction due to septal deviation: an observational study[END_REF] [START_REF] Kimbell | Changes in nasal airflow and heat transfer correlate with symptom improvement after surgery for nasal obstruction[END_REF]. This reinforce the interest and validate the use of CFD. Moreover, 3D reconstruction tools make it possible to achieve architectural modifications and simulate surgical procedures: it's Virtual Surgery (VS) [START_REF] Rhee | Role of Virtual Surgery in Preoperative Planning: Assessing the Individual Components of Functional Nasal Airway Surgery[END_REF]. Thus, association of VS and CFD may predict surgical outcomes and benefits of surgery [START_REF] Do | Predicting Postsurgery Nasal Physiology with Computational Modeling: Current Challenges and Limitations[END_REF]. However, CFD and VS currently raise two questions. Firstly, to date, most CFD studies were conducted in steady flow, a model of infinite inspiration (or expiration) which is different from physiological breathing [START_REF] Quadrio | Review of computational fluid dynamics in the assessment of nasal air flow and analysis of its limitations[END_REF]. Unsteady flows are closer to reality but requires an increased computation time. Secondly, few CFD studies compared VS to postsurgery (PS) models [START_REF] Do | Predicting Postsurgery Nasal Physiology with Computational Modeling: Current Challenges and Limitations[END_REF].

Establishing correlations between VS and PS would increase VS and CFD validity.

The main objective of our study was to compare VS and PS models using unsteady computations in a CFD study conducted on patients with SD. We also compared VS and PS to presurgery models.

MATERIALS AND METHODS

Ethical standards

All patients gave written consent before participating in the study, which was conducted according to the Declaration of Helsinki. Authorization to conduct this study was obtained from local Ethical Committee (APHM N°2017-14-12-005).

Patient selection

We included two patients referred in our Department for nasal obstruction, requiring [START_REF] Radulesco | Geometric morphometric contribution to septal deviation analysis[END_REF]. We used NOSE scores to assess patient impairment 1 month before and 6 months after surgery [START_REF] Floyd | Systematic Review and Metaanalysis of Studies Evaluating Functional Rhinoplasty Outcomes with the NOSE Score[END_REF].

Preoperative CTscans were obtained using a low dose protocol on a multi-detector CTscanner (AquilionPrime, Canon Medical Systems, Otawara, Japan) using a cranio-caudal helical acquisition. In both patients, PS CTscans were performed at least 3 months after surgery and for another medical reason. PS CTscans did not show postoperative complication.

Creation of 3-Dimensional Models 3D reconstructions were obtained using ITK-Snap (3.6.0). The procedure was as follows: 1.

Importation of CTscan images (DICOM formats), 2. Segmentation process using the half maximum height protocol (ImageJ software version 1.44) to determine the boundaries of anatomical structures, and 3. Nasal surface extraction. The half-maximum height protocol offers an objective process to delimit a reproducible interface between two tissues [START_REF] Spoor | Linear measurements of cortical bone and dental enamel by computed tomography: applications and problems[END_REF]. This method included pixels from -1024 to -400 Houndsfield Units. Paranasal sinuses were removed manually during the segmentation process.

Virtual Surgery protocols

In both cases, VS were blinded to operative protocol to avoid influencing anatomical modifications. VS were performed by another surgeon than the one who performed surgery using a combination of ITK-Snap® and Star-CCM+®. By adding air voxels or performing millimetric translations, it was possible to create a new shape aiming to simulate a realistic surgical procedure (Figure 1). Time needed to perform the procedure were evaluated.

CFD protocol CFD was performed using Star-CCM+® software (Siemens®). The computational domain was the Boolean union of extracted surface and a 3 cm-diameter sphere centered on the tip of the nose. The aim of this sphere was to impose a pressure boundary condition far from the nostrils, thus allowing a non-imposed airflow close to the nasal entrance. Volume mesh was performed using a polyhedral mesher with parameters defined after a convergence mesh study: the prism layer (total thickness = 0.16 mm) included 10 prismatic cells with 1.1 prism layer stretching ratio.

We defined the following computational hypothesis: airflow to standard conditions for temperature (19°C for inspiration, 37°C for expiration) and atmospheric pressure (101 325 Pa).

Air was considered to be an incompressible Newtonian fluid with density ρ=1.225 kg/m3 and viscosity μ=1.8x10-5 Pa.s. We considered the flow to be laminar. Sinonasal surface was a nonslip wall. Wall temperature was 34°C [START_REF] Lindemann | Impact of menthol inhalation on nasal mucosal temperature and nasal patency[END_REF]. An imposed atmospheric pressure (reference value) was defined on the sphere surface.

In order to perform a time-dependent computation, a time-varying boundary condition applied to nasopharynx was previously recorded following a normal breathing cycle (Static pressure, Figure 2). t1 was maximum inspiration (∆P=-171 Pa, time=4.15sec), t2 was the transition between inspiration and expiration (∆P=0Pa, time=5.25sec) and t3 was maximum expiration phase (∆P=+126 Pa, time=6.2sec).

Comparison between VS and PS

Comparison between models were evaluated by searching differences and/or establishing correlations between CFD variables. We also compared nasal resistances, airflows and rhinomanometries.

CFD variables measurements

Heat flux, total pressure and wall shear stress were measurement averages within a circular wall region of 1cm 2 . All measurements in axial view were performed 1.5 cm from the nasal floor on the lateral wall on 3 locations: location n°1 was defined at the maximum SD, 2cm after the entry of the nose (for both patients); location n°2, 4cm after the entry of the nose; location n°3, at the choana (7cm after the entry of the nose). Location n°0 was ambient air values.

Temperatures were mean temperature values recorded on cross sectional areas regarding locations n°1, 2 and 3. All these CFD variables were measured at t1, t2 and t3 for each location (n°1, n°2, n°3) on the most obstructed side (right side). We also analyzed streamlines.

Nasal Resistances, airflows and rhinomanometries

Nasal resistances and airflows were calculated for a static pressure differential of -75 Pa on the most obstructed side. Nasal resistances were expressed as R=ΔP/V, where ΔP is the differential pressure between the nostril and the rhinopharynx and V is the airflow rate, expressed in Pa.s/mL.

Anterior Active Rhinomanometries (AA-RMM) (without vasoconstriction) were performed before and after surgery. The same rhinomanometer was used for both patients (Otopront®, Rhinon-sys©, Happersberger, Otopront GmbH2008). CFD-computed RMM curves (CFD-RMM) were performed before surgery, after VS and PS.

Statistical analyses

A t-test was used to assess statistical differences between preoperative, VS and PS CFD variables according inspiration (t1) or expiration (t3) phases. t2 measurements were not used for this statistical evaluation. p-values < 0.05 were considered to be statistically significant.

The Pearson correlation test (r) was used to assess correlations between VS and PS CFD variables using t1, t2 and t3 measurements. Correlations were considered perfect if |r|=1, very strong if |r|>0.8, strong if 0.5>|r|>0.8, moderate if 0.2>|r|>0.5, low if |r|<0.2 and null if r=0.

The 95% confidence intervals of Person's rank correlation coefficients (IC95%) were computed by bootstrapping using the RVAideMemoire package for R3.5 (R Foundation for Statistical Computing, Vienna, Austria). Correlations were statistically significant if p<0.05.

Statistical calculations were performed using PAST software (Oyvind Hammer, University of Oslo, Norway).

RESULTS

Patient management

A septoplasty with Killian incision through endoscopic approach was undergone in both cases.

There was no complication after surgery. For patient 1, NOSE score was 70/100 before and 0/100 after surgery. For patient 2, NOSE score was 90/100 before and 0/100 after surgery. Our VS protocol aimed to simulate the septoplasty (Figure 1). VS duration were 6 and 26 minutes for Patient 1 and 2, respectively.

Comparison between CFD variables (Table 1)

Regarding comparison between VS and PS models, we found no statistical difference regarding heat flux, wall shear stress, total pressure and temperature for inspiration and expiration phases (all p>0.05, Figure 3). Figure 4 shows evolution of temperatures for both patients.

Concerning comparison between presurgery and PS or VS models, preoperative measurements were significantly different from VS and PS measurements during inspiration for heat flux, total pressure and temperature, but not for wall shear stress. We found no difference for expiration phase between models for all CFD variable (Figure 5).

The study of correlations between CFD variables of VS and PS models found very strong correlations for heat flux (r=0.93, IC95=0.78;0.98), total pressure (r=0.96, IC95=0.87;0.99), temperatures (r=0.98, IC95=0.95;0.99) and wall shear stress (r=0.92, IC95=0.75;0.97) (all p<0.01).

Nasal resistances and rhinomanometries (Table 2, Figure 6) VS and PS models provided comparable CFD-NR. Great improvement were found compared to preoperative evaluation in both cases. CFD-NR were quantitively different from rhinomanometry in most cases, while RMM curves shapes were comparable (Figure 6). For Patient 1, the main difference was found regarding maximum pressures reached in One of the major challenge for the surgeon is to choose which structure to modify, aiming the best outcome while seeking to be as conservative as possible [START_REF] Pawar | Advances in Technology for Functional Rhinoplasty[END_REF]. Indeed, extensive surgeries can lead to surgical complications such as empty nose syndrome or healing problems [START_REF] Di | Numerical Simulation of Airflow Fields in Two Typical Nasal Structures of Empty Nose Syndrome: A Computational Fluid Dynamics Study[END_REF].

VS allows to visualize the impact of surgery. If computation time is still a major obstacle to realtime production of CFD models, in the near future, it will undoubtedly be possible to show to the patient which corrections would be feasible and their consequences on airflow. This could provide better preoperative information to the patient (Figure 7).

Virtual Surgery protocols: clinical applicability

Performing VS is fast and feasible whatever the SD. However, it requires learning the use of 3D reconstruction and CFD softwares [START_REF] Burgos | New CFD tools to evaluate nasal airflow[END_REF]. In our protocol, the correction of SD that did not come into contact with the lateral wall SD was easy using Star-CCM+®, applying progressive translations to walls (Patient 2). In case of contact between nasal septum and the lateral wall, it was more convenient to use the 3D-segmentation software by adding air voxels between structures (Patient 1).

Reliability of Virtual Surgery

Our study aimed to compare CFD variables between VS and PS. No differences were found for all CFD measurements between VS and PS. Moreover, we found great correlations between VS and PS for all CFD variables (Table 1). Finding close correlation coefficients appears logical since determination of heat flux, wall shear stress, temperature and total pressure are related to solving the same equations. This reinforce the validity of our measurement procedure.

Regarding inspiration phase, heat flux, total pressure and temperatures were greatly improved after surgeries. Surprisingly, differences for wall shear stress were not statistically significant although values of wall shear stress in location n°1 were frankly improved. This was related to wall shear stress measurements in Locations 2 and 3 that were almost equal for all models (Figure 8). As Li et al., focusing on maximal value of wall shear stress (Peak wall shear stress) could be more relevant when studying SD [START_REF] Li | Investigation of the abnormal nasal aerodynamics and trigeminal functions among empty nose syndrome patients: Abnormal nasal aerodynamics in ENS patients[END_REF].

Comparing VS and PS models, the main difference for Patient 1 appeared in the upper and anterior part of the nasal fossa, where an indentation could be visualized postoperatively (Figure 3, black arrow). This was related to the surgical procedure. A higher correction of SD would have required an extra-corporeal septoplasty with open approach, but we believed our procedure would be compliant for the postoperative patient outcome. Looked in a different way, we may consider that VS overcorrected SD. Indeed, there is no anatomical or technical limit when performing a VS protocol. This explains why results could be slightly better after VS than PS, especially regarding total pressure and wall shear stress for Patient 2. VS protocols require surgical experience to remain realistic and avoid strong overcorrections.

Regarding expiration phase, no statistically significant differences were found for all CFD variables, nor between VS and PS, neither comparing presurgery to VS and PS. However, pressures variations in location n°1 decreased importantly after surgeries, testifying the correction of the anterior SD. We found few differences before and after surgeries regarding heat flux because nasal wall (34°C) and expired air (37°C) temperatures were close leading to a diminution of thermal exchanges compared to inspiration phase (19°C). Otherwise, it appears logical that most of temperatures measurements in expiration phase were similar as SD were anterior. Studying CFD variables during expiration appeared less sensitive than inspiration values, matching with patient's complaint, almost only reported for inspiration.

Regarding rhinomanometries, we found good cohesions between VS and PS CFD-RMM. Most of differences were found regarding maximal pressures reached during AA-RMM, higher than in CFD-RMM. When performing AA-RMM, the patient had one-side nasal breathing and needed, because of their SD, an increased inspiratory effort (∆P). In CFD models, patient had two-sided nasal breathing, leading to a diminution of airflow in each of the two nostrils, this artificially increasing nasal resistance. Deeper breathing cycle and unilateral nasal breathing should be used for a better comparison between rhinomanometries. For Patient 2, we found high nasal resistances with rhinomanometry after surgery (1.04 Pa.s/mL), although she had no nasal complain (NOSE = 0/100). It has already been shown that quantitative comparison between rhinomanometry and patient reported outcome are low [START_REF] Lara-Sánchez | Assessment of nasal obstruction with rhinomanometry and subjective scales and outcomes of surgical and medical treatment[END_REF]. Otherwise, nasal cycle, that was not assumed in our model, can lead to air pathway modification and could explain part of discrepancies between CFD-RMM and rhinomanometries. Some authors proposed virtual correction of nasal cycle: Gaberino et al. showed that correlations with patient perception improved after virtual correction of the nasal cycle [START_REF] Patel | Simulating the Nasal Cycle with Computational Fluid Dynamics[END_REF] , [START_REF] Gaberino | Estimates of nasal airflow at the nasal cycle midpoint improve the correlation between objective and subjective measures of nasal patency[END_REF].

The fact that most data between VS and PS were comparable provides scientific evidence of VS reliability. Other studies must to be carried out on this topic to consolidate these data.

Contribution of unsteady computations

To date, most of studies used steady flows, which correspond to an infinite inspiration (or expiration). This breathing model is quite far from a normal breathing which decomposes four phases: accelerating inspiratory, decelerating inspiratory, accelerating expiratory, and decelerating expiratory [START_REF] Vogt | The new agreement of the international RIGA consensus conference on nasal airway function tests[END_REF].

The main advantage of unsteady computation is to respect the normal breathing on the one hand, but also to provide visualizations close to reality showing inspiratory and expiratory airflows on the other (Video 1). The main disadvantage of unsteady computations is, to date, the computation time required to carry them out (about 5 days for 15 seconds computed).

While waiting for computers improvement, steady computations, providing many information, remains possible. Differences between steady and unsteady computations must be explored in the field of nasal obstruction. Hörscher et al., in a preliminary study, found that major differences between steady and unsteady flows where found at transition from inspiration to expiration [START_REF] Hörschler | On the assumption of steadiness of nasal cavity flow[END_REF].

Comparison with other studies

Few studies evaluated VS in the field of nasal airway obstruction (8) , ( 9) , [START_REF] Garcia | Septal Deviation and Nasal Resistance: An Investigation using Virtual Surgery and Computational Fluid Dynamics[END_REF]. If many have explored nasal physiology with CFD, this is the first, to our knowledge, to compare VS to PS with CFD variables such as heat flux or wall shear stress and using unsteady flows (7) , ( 25) , ( 26) , [START_REF] Ozlugedik | Numerical Study of the Aerodynamic Effects of Septoplasty and Partial Lateral Turbinectomy[END_REF].

Among existing studies, some focused on virtual turbinectomy (middle or inferior) and showed strong modifications in airflow and air conditioning (28) , [START_REF] Dayal | Impact of Middle versus Inferior Total Turbinectomy on Nasal Aerodynamics[END_REF].

In the field of SD, Frank-ito et al. compared VS to presurgery and PS models studying nasal resistances, airflow and pressures in 10 patients [START_REF] Do | Predicting Postsurgery Nasal Physiology with Computational Modeling: Current Challenges and Limitations[END_REF]. They found no statistically significant differences between VS and PS regarding unilateral nasal resistances (p=0.059) and unilateral flow (p=0.074). Comparing presurgery to VS or PS during inspiration, they always found statistically significant differences (p<0.05). They also showed that their VS protocol had a tendency to overestimate surgeries outcomes on most obstructed side concerning airflows. In our study, VS overestimated anterior total pressure improvement (Locations n°1 and n°2).

Rhee et al. tested different VS in a patient presenting multifactorial nasal obstruction and

showed that most of the reduction in nasal resistance was accomplished with performance of septoplasty and inferior turbinate reduction (-50.2% versus -4.9% for nasal valve repair alone) [START_REF] Rhee | Role of Virtual Surgery in Preoperative Planning: Assessing the Individual Components of Functional Nasal Airway Surgery[END_REF]. Obviously, everything depends on the model used. However, this is interesting and shows how, in the future, VS can help the surgeon in her/his therapeutic decision or surgical strategy.

In a previously published study, Rhee et al. also showed that VS provided reliable results [START_REF] Rhee | Toward Personalized Nasal Surgery Using Computational Fluid Dynamics[END_REF].

Using VS, Garcia et al. showed that anterior SD are more involved in nasal obstruction than posterior's [START_REF] Garcia | Septal Deviation and Nasal Resistance: An Investigation using Virtual Surgery and Computational Fluid Dynamics[END_REF]. However, these studies included only one patient, evaluated only nasal resistances, airflow or pressure measures in steady flows. Other CFD variables i.e. heat flux, are important to take into account while performing simulations and some authors highlighted how repetition of the respiratory cycle during simulation is essential to understand nasal physiology and avoid information loss (6)(7) [START_REF] De Gabory | Numerical simulation of two consecutive nasal respiratory cycles: toward a better understanding of nasal physiology: Simulation of nasal physiology[END_REF].

Perspectives

Using unsteady flows brings CFD closer to reality. However, many progresses remain to be made. Some authors incorporated humidity in CFD models [START_REF] Balakin | Aerodynamic evaluation of the empty nose syndrome by means of computational fluid dynamics[END_REF]. Hygrometry must be taken into account toward humidification remains one of the first functions of the nose and could play a major role in nasal patency. CFD-calculated humidity may be correlated to patient perception, especially in case of empty nose syndrome.

The most important challenge is the creation of deformable model. During breathing, architectural modifications appears during inspiration and expiration. More or less important according to patients, cartilages and soft tissues deformability has already been tested in vivo and is characterized by the Young's modulus [START_REF] Bonaparte | Assessment of Pliability and Elasticity of the External Nasal Skin in Patients With Unilateral Nasal Valve Collapse: A Static Biomechanical Evaluation[END_REF] [START_REF] Griffin | Biomechanical Characterization of Human Soft Tissues Using Indentation and Tensile Testing[END_REF]. We are currently working on a deformable CFD model.

Otherwise, this technology needs automating to be achievable in clinical routine. Deep learning techniques could help for the model reduction of the Navier-Stokes equations for unsteady flow problems [START_REF] Miyanawala | An Efficient Deep Learning Technique for the Navier-Stokes Equations: Application to Unsteady Wake Flow Dynamics[END_REF]. When incorporating these advancements, combination of VS and CFD could be, in the near future, an unavoidable tool for the ENT surgeon.

Limits

Our study included only two patients. This was due to the lack of postoperative CTscans. Indeed, all CTscans were performed for another reason after surgery, as we do not perform postoperative imaging for nasal airway obstruction.

Moreover, computation time are still long. However, steady computations are much faster and can be realized, upon request, in case of challenging situation. This is our practice today. CONCLUSION VS and PS models provided comparable results, matching with the clinical evaluation of our patients. In the near future, CFD and VS could be important tools to master for the ENT surgeon.

Nowadays, we believe that CFD allows to better understand nasal physiology and to anticipate the surgical strategy to adopt before surgery for selected cases. After surgery, in inspiration, most of the airflow was on middle pathway. However, in expiration, air pathway was more balanced (upper, middle and inferior pathway).

  septoplasty. Patient 1 (23-years-old men) was concerned with a major right-sided posttraumatic SD. Patient 2 (27-years-old female) had a congenital SD. Both were anterior rightsided and C-shaped SD according to Teixeira et al. or Radulesco et al. (2)

  preoperative AA-RMM compared to CFD-RMM (-600 Pa versus -171 Pa). For Patient 2, rhinomanometry showed still high nasal resistance postoperatively (1.04 Pa.s/mL). DISCUSSION Synopsis of new findings CFD offers new perspectives in nasal airway analyses by visualization and measurements on airflow. With VS, the surgeon can try different corrections facing her/his technical possibilities.

Figure 1 .

 1 Figure 1. Comparison of surgeries for Patient 1. Three-dimensional anterior and right inferior

Figure 2 .

 2 Figure 2. Breathing cycle used for unsteady computation: time-varying static pressure

Figure 3 .

 3 Figure 3. Computations during inspiration phase (t1) for Patient 1 (right nasal fossa). We

Figure 4 .

 4 Figure 4. Evolution of temperatures for Patient 1 and Patient 2 before and after surgeries

Figure 5 .

 5 Figure 5. Right nasal fossa of Patient 2. Representation of total pressure before surgery and

Figure 6 .

 6 Figure 6. Comparison between CFD-RMM (top) and AA-RMM (bottom), before surgery and

Figure 7 .

 7 Figure 7. Visual comparison of airflow during expiration (t3) before and after surgeries for

Figure 8 .

 8 Figure 8. Visualization of wall shear stress for patient 2. WSS values for Location n°2 and 3 are

Table 1 .

 1 Statistical analyses comparing CFD Variables according to inspiration or expiration

					Pearson correlation coefficient
	Statistical test		t test			(VS / PS)
	Data		p values		r	IC95%	p values
			Inspiration (t1) Expiration (t3)
	Heat flux	VS / PS	0.7	0.9	0.9361	0.78;0.98	<0.01*
		Pre-op / VS	0.01*	0.2
		Pre-op / PS	0.03*	0.17
	Total pressure	VS / PS	0.17	0.9	0.96	0.87;0.99	<0.01*
		Pre-op / VS	0.009*	0.19
		Pre-op vs PS	0.014*	0.14
	Temperature	VS / PS	0.07	0.6	0.98	0.95;0.99	<0.01*
		Pre-op / VS	0.007*	0.6
		Pre-op / PS	0.008*	0.1
	Wall shear	VS / PS	0.6	0.054	0.92	0.75;0.97	<0.01*
	stress	Pre-op / VS	0.2	0.23
		Pre-op / PS	0.1	0.2
	phases and locations. Statistical tests included CFD values of Patient 1 and 2. Pre-op=
	preoperative. VS = Virtual Surgery; PS = Postsurgery;
	* statistically significant test

Table 2 .

 2 Values of inspiratory nasal resistances and airflows on the obstructed side obtained with anterior active rhinomanometry (top) or CFD (bottom), for a pressure differential of -75 Pa after VS and RS. Nasal resistances in Pa.s/mL, Airflows in mL/s

2019]; Disponible sur: https://arxiv.org/abs/1710.09099v3