Alice Sinatra

Jean-Christophe Dornstetter

Yvan Castin

Spin squeezing in Bose-Einstein condensates: Limits imposed by decoherence and non-zero temperature

Keywords: numbers: 03.75.Gg, 42.50.Dv, 03.75.Mn. Contents

We consider dynamically generated spin squeezing in interacting bimodal condensates. We show that particle losses and non-zero temperature effects in a multimode theory completely change the scaling of the best squeezing for large atom numbers. We present the new scalings and we give approximate analytical expressions for the squeezing in the thermodynamic limit. Besides reviewing our recent theoretical results, we give here a simple physical picture of how decoherence acts to limit the squeezing. We show in particular that under certain conditions the decoherence due to losses and non-zero temperature acts as a simple dephasing.

I. INTRODUCTION

Spin squeezing is about creating quantum correlations in a many-body system that can be useful for metrol-ogy. An example is that of atomic clocks, where the ultimate signal to noise ratio, once all the technical noise has been eliminated, can be improved by manipulating and controlling the system at the level of its quantum fluctuations.

A. Spin squeezing and atomic clocks

The aim of an atomic clock is to measure precisely the energy difference between two atomic states a and b that are for example two hyperfine states of an alkali atom. To explain how the clock works and to introduce spin squeezing, we shall describe the ensemble of N atoms used in the clock using the picture of a "collective spin" that evolves on the so-called Bloch sphere. The collective spin is simply the sum of the effective spins 1/2 that describe the internal degrees of freedom of each atom. In the second quantized formalism the three hermitian spin components S x , S y and S z are defined by:

S x + iS y = a † b , S z = a † a -b † b 2 , (1)
where a † and b † are creation operators of particles in the internal states a and b. The spin operators are dimensionless and obey the commutation relations [S x , S y ] = iS z and cyclic permutations. For the moment we do not care about the external degrees of freedom of the atoms. The component S z of the collective spin is half the population difference between states a and b, while S x and S y describe the coherence between these states. If the N atoms are prepared in a coherent superposition of states a and b with relative phase 2φ:

|φ N = 1 √ N ! e iφ a † + e -iφ b † √ 2 N |0 , ( 2
)
where |0 is the vacuum, the collective spin lies on the equatorial plane of the Bloch sphere, pointing at an angle -2φ with respect to the x axis. For non-interacting atoms, the further evolution is ruled by the Hamiltonian and the spin precesses around the z axis with the Larmor frequency ω ab . The atomic clock measures the phase accumulated by the collective spin during a long precession time τ . From this phase the frequency ω ab is deduced. Atomic clocks are nowadays so precise that they are sensitive to the quantum noise of the collective spin. If for example the spin is initially prepared along the x axis (in an eigenstate of S x ), it has necessarily fluctuations in the transverse components S y and S z such that:

H 0 = ω ab S z (3)
S x = N 2 , ∆S y ∆S z ≥ 1 2 | S x | . (4)
In particular fluctuations of S y introduce a statistical variance on the accumulated phase during the precession time and on the measured frequency ω ab . In the uncorrelated state [START_REF] Kitagawa | Squeezed spin states[END_REF] with φ = 0, root mean square fluctuations of S y and S z are equal: ∆S y = ∆S z = √ N /2. In a Ramsey measurement with interrogation time τ , these quantum fluctuations introduce the root mean square fluctuations of the measured frequency equal to [START_REF] Wineland | Squeezed atomic states and projection noise in spectroscopy[END_REF]:

∆ω unc ab = 1 √ N τ . (5)
This noise coming from quantum fluctuations, intrinsic to the initial state where each atom is in a superposition of a and b, is known in clocks as the "partition noise". The idea of spin squeezing [START_REF] Kitagawa | Squeezed spin states[END_REF] is that the Heisenberg relation (4) allows to reduce ∆S y provided that ∆S z is increased. This idea is illustrated in Fig. 1. To quantify the spin squeezing we use the parameter ξ 2 introduced in [1]:

ξ 2 = N ∆S 2 ⊥,min | S | 2 (6)
where N is the total atom number, ∆S 2 ⊥,min is the minimal variance of the spin orthogonally to its mean value S . The state is squeezed if and only if ξ 2 < 1. As explained in [START_REF] Wineland | Squeezed atomic states and projection noise in spectroscopy[END_REF], ξ directly gives the reduction of the statistical fluctuations of the measured frequency ω ab with respect to uncorrelated atoms, for the same atom number N and the same Ramsey time τ :

∆ω sq ab = ξ∆ω unc ab = ξ √ N τ . ( 7
)
The parameter ξ in Eq.( 6) is in fact the properly normalized ratio between the "noise" ∆S ⊥,min and the "signal" | S |. In experiments ∆S ⊥,min is directly measured by measuring S z after an appropriate state rotation and | S | is separately deduced from the Ramsey fringes contrast.

B. State of the art

On one hand the most precise atomic clocks using microwave transitions in cold alkali atoms have already reached the quantum partition noise limit with atom numbers up to N = 6 × 10 5 [3]. On the other hand, very recently a significant amount of spin squeezing, up to -8 dB (ξ 2 = 10 -0.8 ) [START_REF] Gross | Nonlinear atom interferometer surpasses classical precision limit[END_REF] was measured in dedicated, proof-of-principle experiments. In [START_REF] Leroux | Implementation of cavity squeezing of a collective atomic spin[END_REF] squeezing was created in a large sample of N = 5 × 10 4 atoms with a feedback mechanism in a resonant optical cavity, while in [START_REF] Gross | Nonlinear atom interferometer surpasses classical precision limit[END_REF] and in [START_REF] Riedel | Atom-chip-based generation of entanglement for quantum metrology[END_REF] the squeezing was created in smaller samples, of order N = 10 3 , using atomic interactions in bimodal condensates. The ultimate limits of the different paths to spin squeezing are still an open question. Here we concentrate on a dynamical scheme using interactions in bimodal condensates [START_REF] Gross | Nonlinear atom interferometer surpasses classical precision limit[END_REF][START_REF] Riedel | Atom-chip-based generation of entanglement for quantum metrology[END_REF][START_REF] Sorensen | Manyparticle entanglement with Bose-Einstein condensates[END_REF][START_REF] Uffe | Positive-P simulations of spin squeezing in a two-component Bose condensate[END_REF] and analyze in particular the influence of dephasing, decoherence and non-zero temperature on this squeezing scheme.

C. Two-mode scalings without decoherence

We consider for simplicity a bimodal condensate with identical interactions in the components a and b with coupling constants g aa = g bb = g and no crossed a-b interactions [23]. We assume that the initial state is the factorized state (2) with φ = 0 and a fixed total number of atoms N . In a two-mode picture, interactions introduce a Hamiltonian that is non-linear in the spin operator:

H nl = χS 2 z . (8)
The quadratic form ( 8) is obtained expanding the system Hamiltonian to second order around the average numbers of particles in components a and b, Na and Nb , both equal to N/2 for the initial state (2) [START_REF] Castin | Relative phase of two Bose-Einstein condensates[END_REF][START_REF] Sinatra | Phase dynamics of Bose-Einstein condensates: Losses versus revivals[END_REF]. χ is thus the derivative of the chemical potential with respect to the particle number in each component

χ = dµ a /dN a = dµ b /dN b evaluated in N a = Na , N b = Nb .
The general expression of the expanded Hamiltonian including drift terms for non-symmetric interactions, nonsymmetric splitting or fluctuations in the total particle number can be found in [START_REF] Sinatra | Phase dynamics of Bose-Einstein condensates: Losses versus revivals[END_REF][START_REF] Yun | Spin squeezing in a bimodal condensate: spatial dynamics and particle losses[END_REF]. We are interested in the best squeezing that can be obtained in the thermodynamic limit for a spatially homogeneous system: we therefore explicitly write χ in terms of the interaction constant g and the volume V of the system:

N → ∞ , ρ = N V = constant, (9)
µ a = gN a V , µ b = gN b V , χ = g V . ( 10
)
We can consider the non-linear Hamiltonian (8) as a Hamiltonian of the form (3) with a Larmor frequency ω ab that depends itself on S z . As explained in [START_REF] Kitagawa | Squeezed spin states[END_REF] and shown in Fig. 2 (left), H nl "twists" the transverse spin fluctuations and generates spin squeezing. However, in order to minimize ξ 2 , the evolution should not go to far. E.g. when the fluctuations become too much distorted and start to wrap around the Bloch sphere, the "signal" | S | in the denominator of ξ 2 decreases and the squeezing parameter increases again. A sketch of the time dependence of ξ 2 when the state (2) evolves under the influence of H nl [START_REF] Uffe | Positive-P simulations of spin squeezing in a two-component Bose condensate[END_REF] is given in Fig. 2 (right). We name "best squeezing time" t min the time that minimizes ξ 2 (t) and "best squeezing" ξ 2 min the corresponding squeezing. The analytical expression for the squeezing as a function of time is given in [START_REF] Kitagawa | Squeezed spin states[END_REF]. Introducing an appropriate rescaling of the time variable as in [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF], this gives the following scalings of the best squeezing and the best squeezing time for N ≫ 1 and ρ,g constants:

(ξ 2 min ) H nl ≃ 3 2/3 2 1 N 2/3 ; ρgt min H nl ≃ 3 1/6 N 1/3 . (11)
This differs from the original prediction in [START_REF] Kitagawa | Squeezed spin states[END_REF] by numerical factors. Subsequent studies [START_REF] Sorensen | Manyparticle entanglement with Bose-Einstein condensates[END_REF][START_REF] Søndberg | Bogoliubov theory of entanglement in a Bose-Einstein condensate[END_REF] gave indications that the achievement of large squeezing in condensates should be possible even in presence of decoherence, but were not able to confirm or disprove the scalings [START_REF] Yun | Spin squeezing in a bimodal condensate: spatial dynamics and particle losses[END_REF].

D. New scalings in presence of decoherence

We will show in sections III and IV that the scalings [START_REF] Yun | Spin squeezing in a bimodal condensate: spatial dynamics and particle losses[END_REF] are disproved when decoherence coming from particle losses or non-zero temperature is included in the description. To summarize, instead of tending to zero when N → ∞ in the thermodynamic limit, ξ 2 min tends to a positive constant

ξ 2 min lim.therm. → constant > 0 . ( 12
)
Concerning the best squeezing time t min we distinguish two cases. In the case of particle losses t min is finite in the thermodynamic limit and scales as

ρgt min ∝ 1 ξ 2 min . ( 13
)
In the case of finite temperature t min cannot be calculated within our analytical treatment that neglects interactions among Bogoliubov modes. Nevertheless, we introduce a "close-to-best" squeezing time t η defined in equation (37), at which ξ 2 approaches ξ 2 min with a finite precision η, that scales as

ρgt η ∝ 1 ηξ 2 min ( 14
)
provided that t η remains smaller than the typical collision time among Bogoliubov modes. We will show that these new scalings in presence of decoherence are quite general and that the physics of how decoherence acts is caught by a very simple dephasing model that we shall solve exactly and study in detail in the next section.

II. DEPHASING MODEL

In this section we consider a dephasing Hamiltonian model of the form

H = χ(S 2 z + DS z ) ( 15
)
where D is a Gaussian real random variable of zero mean. We assume here that D is time independent, but it varies randomly from one experimental realization to the other mimicking a stationary random dephasing environment. We also assume that D has a variance of the order of N for N large [24]:

D 2 N → ǫ noise , N → ∞. (16)
Finally ǫ noise , finite in the thermodynamic limit, is a small parameter of the theory and we limit ourselves in general to first order in this quantity. An exception is made in subsection II B where expressions to all orders in ǫ noise are given.

Starting with the initial state (2) with φ = 0, we will show that this minimal model reproduces the scalings ( 12) and ( 14). In the subsequent sections III and IV we will detail an analogy between the dephasing model and microscopic models accounting for the effect of particle losses or of non-zero temperature on squeezing. The parameter ǫ noise introduced here (16) will then be related to the lost fraction of particles or the populations of thermally excited modes, respectively.

A. Squeezing in the thermodynamic limit For the symmetric case we consider, the mean spin is always aligned along x. The minimum transverse spin variance is

∆S 2 ⊥,min = 1 2 S 2 y + S 2 z - ( S 2 y -S 2 z ) 2 + {S z , S y } 2 , ( 17
)
where the expectation values . . . represent the average over the quantum state and over the random variable D.

The notation {, } stands for the anticommutator. Introducing quantities A and B,

A = S 2 y - N 4 (18)
B = {S z , S y } (19) ∆S 2 ⊥,min = 1 2 N 2 + A -A 2 + B 2 . ( 20
)
To derive the scalings ( 12) and ( 14), and to have a physical insight, it is convenient to reason in terms of the phases of the operators a and b:

a = e iθa N a , [N a , θ a ] = i (21) b = e iθ b N b , [N b , θ b ] = i , (22)
where N a = a † a and N b = b † b. This is a legitimate representation as long as the condensate modes have a negligible probability of being empty [START_REF] Carruthers | Phase and angle variables in quantum mechanics[END_REF]. By neglecting the fluctuations of their modulus, the collective spin components S x , S y are simply given by

S x ≃ Re N 2 e -i(θa-θ b ) , (23)
S y ≃ Im N 2 e -i(θa-θ b ) . ( 24
)
At t = 0, the phase difference (θ aθ b ) has zero mean and root mean square fluctuations that scale as 1/ √ N . Both S y and S z scale as √ N . As a consequence, for N large we can expand the exponentials in (23)-(24). To lowest order we then have

S x ≃ N 2 , S y ≃ - N 2 (θ a -θ b ) , S z = N a -N b 2 . (25)
S y and S z are then simply proportional to the position operator Q and the momentum operator P of a fictitious free particle. As we will see, the expansions (25) remain valid for times (ρgt/ ) ≪ √ N . The squeezing occurs because in a given realization of the experiment S y becomes an enlarged copy of S z . Indeed after the pulse, for t > 0, from the Heisenberg equations of motions for the phase operators, with χ = g/( V ), one has

(θ a -θ b )(t) = (θ a -θ b )(0 + ) - gt V [2S z + D] . ( 26
)
As the squeezing dynamics goes on, S y that was initially of the same order as S z , grows linearly in time while S z stays constant. Correspondingly A ∝ t 2 and B ∝ t.

Since {S y (0), S z } = 0, one actually has at all times:

S y = S y (0) + S lead y t ; (27) A N = αt 2 ; (28) B N = βt , (29)
where we have introduced the time independent operator and coefficients

S lead y = ρg 2 [2S z + D] (30) α = (S lead y ) 2 N = ρg 2 2 (1 + ǫ noise ) (31) β = {S lead y , S z } N = ρg 2 . ( 32
)
Using the fact that from (25) S x ≃ N/2 for (ρgt/ ) ≪ √ N , and expanding the expression (20) for t ≫ /(ρg), we have in the thermodynamic limit

ξ 2 (t) = t→∞ 1 - β 2 α + β 4 4α 3 1 t 2 + O(t -4 ) . (33)
Using equations ( 30)-( 32), with ǫ noise ≪ 1, we finally obtain in the long time limit

ξ 2 (t) = ǫ noise + ρgt 2 [1 + O (ǫ noise )] + O 4 (ρgt) 4 . ( 34
)

Best squeezing and close-to-best time

According to (34), the best squeezing in the thermodynamic limit to leading order in ǫ noise is

ξ 2 min = ǫ noise = lim N →∞ D 2 N . (35)
Remarkably, the best squeezing (35) only involves the part of the phase difference D that is not proportional to S z . To understand physically this result we rewrite

ξ 2 min = 1 - β 2 α = (S lead y ) 2 S 2 z -{S lead y , S z }/2 2 (S lead y ) 2 S 2 z .
(36) Due to interactions, through the phase difference (26), S lead y is proportional to (2S z + D). In the absence of the term D this allows a perfect cancellation between the correlation S lead y S z 2 and the product (S lead y ) 2 S 2 z in (36) leading to ξ 2 min = 0 in the limit N → ∞. In presence of D, this is not possible and ξ 2 min has a non-zero limit. Considering the next to leading order in the time expansion of the squeezing parameter Eq.(34), the best squeezing is reached in an infinite time (in the thermodynamic limit). However as we will see ξ 2 (t) is quite flat around its minimum, and it suffices to determine a "close-to-best" squeezing time t η defined as

ξ 2 (t η ) = (1 + η)ξ 2 min , η > 0 . (37)
Then, according to (34), t η is given by

ρg t η = 1 ηξ 2 min . ( 38
)
The close-to-best squeezing time t η is thus very simply related to the best squeezing ξ 2 min . The important point is that t η is finite (non-infinite and non-zero) in the thermodynamic limit.

Geometrical interpretation

We give here a geometrical and pictorial interpretation to the squeezing process in the thermodynamic limit. Let us introduce the rescaled transverse spin components

Y = S y S 2 y , Z = S z S 2 z ( 39
)
that verify Y = Z = 0 and Y 2 = Z 2 = 1. At short times, the Wigner function representing the probability distribution of Y and Z is approximately Gaussian

W (y, z) ∝ exp - 1 2 (y, z) M -1 y z ( 40
)
where M is the covariance matrix

M = Y 2 1 2 {Y, Z} 1 2 {Y, Z} Z 2 . ( 41
)
We can represent graphically the fluctuations of Y and Z by drawing isocontours of W (y, z). Let us introduce the eigenvalues of M

λ 1,2 = 1 ± 1 2 {Y, Z} . (42)
and a rotated coordinate system y ′z ′ aligned with the eigenvectors of M :

y ′ = (z + y)/ √ 2 ; z ′ = (z -y)/ √ 2 . ( 43
)
The points in the y ′ z ′ plane such that

y ′2 λ 1 + z ′2 λ 2 = 1 (44)
form an ellipse whose semi-axis gives the mean square fluctuations of Y ′ and Z ′ that are time dependent linear combinations of S y and S z . The ellipse surface S ellipse = π(λ 1 λ 2 ) 1/2 divided by π is equal to the square root of the determinant of M and is thus asymptotically equivalent to ξ min according to (36):

ξ min = S ellipse π for t ≫ ρg . ( 45
)
In Fig. 3 (left) we show the time dependence of the squeezing parameter ξ 2 in presence and in absence of decoherence. On the same plot we show the determinant of M that asymptotically gives the value of ξ 2 min (36). In Fig. 3 (right) we show the isocontours of W (y, z) defined by (44) at different times for the case with decoherence. As the dynamics goes on, Y and Z become more and more correlated and the ellipse shrinks. In the absence of dephasing and in the thermodynamic limit the ellipse would collapse into a segment in the y = z direction. In the presence of decoherence the process is "blocked" and the ellipse keeps a finite width with a limit area S ellipse = πǫ noise .

B. Exact solution of the dephasing model

The dephasing model ( 15) is exactly solvable. One first writes the Heisenberg equations of motion for a and b, e.g.

i ȧ = χ 2 2S z + D + 1 2 a . (46)
Then one uses the fact that S z is a constant of motion to integrate the equations [START_REF] Sinatra | Phase dynamics of Bose-Einstein condensates: Losses versus revivals[END_REF][START_REF] Wright | Collapses and revivals in the interference between two Bose-Einstein condensates formed in small atomic samples[END_REF]. One obtains (see also [START_REF] Ferrini | Effect of phase noise on useful quantum correlations in Bose Josephson junctions[END_REF]):

ξ 2 (t) = N 2 N 2 + A - √ A 2 + B 2 C 2 (47)
with

A = N (N -1) 8 1 -e -2(χt) 2 D 2 (cos 2χt) N -2 (48) B = N (N -1) 2 sin χt e -1 2 (χt) 2 D 2 (cos χt) N -2 (49) C = S x = N 2 e -1 2 (χt) 2 D 2 (cos χt) N -1 . (50)
A first application of the exact solution (47)-( 50) is to determine the best squeezing ξ 2 min in the thermodynamic limit, to all orders in the dephasing parameter ǫ noise . To this aim we take the limit N → ∞ in (47)-(50) at fixed time t, density ρ and noise parameter ǫ noise = D 2 /N . We find that A/N , B/N and C/N have a finite limit and that

ξ 2 (t) lim.therm. → 1-   1 2 (1 + ǫ noise ) + 1 4 (1 + ǫ noise ) 2 + ρgt 2   -1 . (51)
From this solution one gets the best squeezing and the close-to-best squeezing time t η (for η ǫ noise < 1):

ξ 2 min lim.therm. → ǫ noise 1 + ǫ noise (52) ρg t η lim.therm. → 1 -η ǫ noise (1 + ǫ noise ) √ η ǫ noise . ( 53
)
Note that one can obtain (35) and ( 38) from ( 52) and (53) by linearizing for small ǫ noise .

In Fig. 4 we show ξ 2 (t) as a function of time for a large atom number. The curve is indeed quite flat around the best squeezing time t min . There are two solutions to equation (37): t η < t min and t ′ η > t min . When N → ∞, t η is finite and given by (53). On the other hand, as we show in Appendix A, t min diverges as N 1/4 and t ′ η diverges as N 1/2 . Knowing the asymptotic behavior of t min , by introducing appropriate rescalings of the time variable as in [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF], it is also possible to obtain the first finite size correction to ξ 2 min . This is given in equation (A4) of Appendix A.

C. Squeezing in the weak dephasing limit

We can use the exact solution of the dephasing model (47)-(50) to investigate the weak dephasing limit that we define as the limit where the noise D remains bounded for N → ∞ (see the footnote before equation ( 16)):

D 2 → constant N → ∞. ( 54
)
In this case the scaling of the best squeezing time with N is still given by ρgt min / ∝ N 1/3 as in the case without decoherence [START_REF] Yun | Spin squeezing in a bimodal condensate: spatial dynamics and particle losses[END_REF] and, using the same rescaling of time as in [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF], one has

ξ 2 min = 3 2/3 2 1 N 2/3 + 3 2 + D 2 N + o 1 N , (55)
ρgt min = 3 1/6 N 1/3 - √ 3 4 + o(1) . (56)

III. PARTICLE LOSSES

In this section we consider particle losses that are an intrinsic source of decoherence in condensed gases. Among those, one-body losses are due to collisions of condensate atoms with residual hot atoms due to imperfect vacuum. More fundamental in dense samples are three-body losses where, after a three-body collision, two atoms form a molecule and the third atom takes away the energy to fulfill energy and momentum conservation. After such a collision event the three atoms are lost. Three-body losses are present due to the metastable nature of ultra cold gases, whose real ground state at such low temperatures would be a solid and whose gaseous phase is maintained because the sample is very dilute. Finally two-body losses can also be present, caused by two-body collisions that change the internal state of the atoms. For a trapped gas, we have shown theoretically [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF] that the best achievable squeezing within a two-mode model at zero temperature in presence of one, two and three-body losses can in principle be very large (squeezing parameter ξ 2 of the order of 10 -4 ) provided that the harmonic trapping potential is optimized and a careful choice of the internal state of the atoms is made. To realize such conditions that minimize losses remains however an experimental challenge.

In this section we recall the main results of [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF] concerning the squeezing in presence of particle losses, and we use these results to show an analogy between the effect of the losses and the effect of the dephasing Hamiltonian (15) in the thermodynamic limit.

A. Monte Carlo wave functions

We consider two spatially separated, symmetric condensates. For a more general treatment, please refer to [START_REF] Yun | Spin squeezing in a bimodal condensate: spatial dynamics and particle losses[END_REF]. Initially the system is in the eigenstate of S x with maximal eigenvalue N/2. Besides the non-linear Hamiltonian for the two bosonic modes a and b given by

H nl = χS 2 z with χ = (∂ Na µ a ) Na / , (57)
we include one, two and three-body losses. Due to the losses, the system is "open" and we shall describe it with a density operator that obeys a Master Equation of the Lindblad form [START_REF] Sinatra | Phase dynamics of Bose-Einstein condensates: Losses versus revivals[END_REF]. In the interaction picture with respect to H nl :

dρ dt = 3 m=1 ǫ=a,b γ (m) c m ǫ ρc †m ǫ - 1 2 {c †m ǫ c m ǫ , ρ} (58)
where c ǫ = ã, b for ǫ = a, b and:

ã = e i H nl t a e -i H nl t ; b = e i H nl t b e -i H nl t . ( 59
)
Here the operators a and b are in the Schrödinger picture. The m-body loss rates γ (m) are defined in terms of the so-called rate constants K m as

γ (m) = K m m d 3 r|φ(r)| 2m (60)
where φ(r) is the condensate wave function in mode a or b for the initial atom number (weak loss approximation), so that for example:

d dt N a = -K 1 + K 2 N a d 3 r|φ(r)| 4 + K 3 N 2 a d 3 r|φ(r)| 6 N a . ( 61
)
It is convenient to rephrase the Master Equation (58) in terms of Monte Carlo wave functions [START_REF] Mølmer | A Monte-Carlo wave function method in quantum optics[END_REF]. In this picture pure states evolve deterministically under the influence of an effective Hamiltonian H eff acting during time intervals τ i = t it i-1 separated by random quantum jumps (described by the jump operators S ǫ ) occurring at times t i as illustrated in Fig. 5:

H eff = - ǫ=a,b i 2 γ (m) c †m ǫ c m ǫ , S ǫ = γ (m) c m ǫ . ( 62
)
More precisely the evolution of the non-normalized state vector |ψ(t) in between quantum jumps is given by

i d dt |ψ(t) = H eff (t)|ψ(t) (63)
and the effect of a quantum jump in the component ǫ i at time t i is Quantum averages of any atomic observables O are obtained by summing over all the possible trajectories of the non-normalized state vector [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF]:

|ψ(t + i ) = S ǫi (t i )|ψ(t - i ) . (
Ô = k 0<t1<t2<•••t k <t dt 1 dt 2 • • • dt k {ǫj } ψ(t)| Ô|ψ(t) . ( 65
)

B. Losses randomly kick the relative phase

Let us consider the action of a quantum jump over a phase state (2) with N for example the loss of one particle in state a or b at time t:

c a (t)|φ N = N 2 e -i χ 4 t e iφ |φ -χt/2 N -1 (66) c b (t)|φ N = N 2 e -i χ 4 t e -iφ |φ + χt/2 N -1 . ( 67
)
Under the action of a jump a phase state remains a phase state. On the other hand, the relative phase is shifted by a random amount that depends on the time of the jump and has a random sign depending on whether the jump was in a or b. This behavior is illustrated in Fig. 6 taken from [START_REF] Sinatra | Phase dynamics of Bose-Einstein condensates: Losses versus revivals[END_REF] where we plot the modulus squared of a relative phase distribution amplitude c(φ, t) at t = 0 and t = 2π/χ for three single Monte Carlo realizations. At this particular time (second revival time) the coherent evolution due to H nl has no effect and we can isolate the action of the losses. In the case of Fig. 6, as χt is of the order of unity, the shift in the relative phase due a single jump is large. In the case of squeezing with N ≫ 1, χt ≪ 1, the shifts due to single quantum jumps are on the contrary very small. These shifts nevertheless limit the maximum squeezing achievable as we shall see.

C. Spin squeezing limit and the lost fraction

Let us consider the case of one-body losses only with a loss rate constant γ equal in states a and b. In this case the effective Hamiltonian does not depend on time; referring to the time sequence in Fig. 5, we have explicitly

|ψ(t) = e -iH eff (t-t k )/ S ǫ k (t k )e -iH eff τ k / S ǫ k-1 (t k-1 ) . . . . . . S ǫ1 (t 1 )e -iH eff τ1/ |ψ(0) , ( 68
)
and there is an explicit analytical solution for the generated spin squeezing as a function of time [START_REF] Yun | Spin squeezing in a bimodal condensate: spatial dynamics and particle losses[END_REF][START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF]. Here we use this solution to find the best squeezing in presence of losses in the thermodynamic limit with N → ∞, ρ and γ constant. For simplicity we also assume that the fraction of lost particles at the relevant time t remains small in the thermodynamic limit:

γt ≡ ǫ loss . ( 69
)
We proceed similarly as we did to derive Eq.(51), taking the thermodynamic limit in the exact solution in presence of losses. Restricting for simplicity to the leading order in γt and to the case ρgt/ ≫ 1, we obtain

A N e -γt ≃ ρgt 2 2 1 - 5 3 γt (70)
B N e -γt ≃ ρgt 2 (1 -γt) (71)
ξ 2 (t) ≃ γt 3 + ρgt 2 [1 + O(γt)] , (72)
where A and B are defined in ( 18)- [START_REF] Sinatra | Limit of spin squeezing in finite-temperature Bose-Einstein condensates[END_REF]. Equation (72) shows that for long times, the squeezing parameter is asymptotically equivalent to one third of the lost fraction of atoms. Minimizing (72) with respect to time we obtain

ξ 2 min = 3 4 4 3 γ ρg 2/3 (73) ρg t min = 3 4 ρg γ 1/3 . (74)
In Fig. 7 we show the squeezing in presence of losses and we compare the exact solution [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF] with the approximate expression (72) valid at long times and in thermodynamic limit.

D. Analogy between dephasing and losses

Our goal here is to make an analogy between the dephasing model of section II and the present model with losses. Indeed the relative phase is perturbed by the losses. In presence of one-body losses, a non normalized Monte Carlo wave function after k jumps has the form:

|ψ(t) = N 1/2 e i φ|φ + χt 2 D N -k , (75)
N = k i=1 e -γti e -γ(N -k)t γ k × N (N -1) . . . (N -k + 1) 2 k , (76)
D = 1 t k l=1 t l (δ ǫ l ,b -δ ǫ l ,a ) . (77)
The factor N given by ( 76) is the norm squared of the wave function that is needed to calculate quantum averages, the first line in (76) is due to the effective Hamiltonian evolution, φ is an irrelevant phase and D given by (77) is a random perturbation of the relative phase 2φ that plays the role of the quantity D in Eq.( 26) of the dephasing model of section II, as appears from the fact that exp(-iχtDS z )|φ N = |φ -χtD/2 N . Contrarily to D, D given by ( 77) is time dependent. As detailed in Appendix B, using (65) we can calculate D 2 . To first order in γt we obtain:

D 2 N ≃ γt 3 . (78)
We have thus shown that ξ 2 (t) is asymptotically equivalent to D 2 /N as it the case in the dephasing model. We can then establish an analogy between the model with losses and the dephasing model as summarized in the first two columns of the Table I. In Fig. 8 we show a comparison of D 2 /N obtained from a Monte Carlo simulation, from the exact expression (B7), and from the approximate expression (78) valid to first order in ǫ loss = γt.

E. Optimum squeezing in a harmonic trap

A similar analysis can be performed for a trapped system, where we now consider the two components a and b in identical and spatially separated harmonic traps. In particular equation ( 19) of [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF] is very similar to (72) except that in [START_REF] Sinatra | Limit of spin squeezing in finite-temperature Bose-Einstein condensates[END_REF] of [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF] we have N χ instead of ρg in (72). Another difference is that the more general equation [START_REF] Sinatra | Limit of spin squeezing in finite-temperature Bose-Einstein condensates[END_REF] in [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF] that includes two and three-body losses besides one-body losses, is derived performing an approximation on the effective Hamiltonian: the constant loss rate approximation [25]

H eff = - ǫ=a,b i 2 γ (m) c †m ǫ c m ǫ ≃ - ǫ=a,b i 2 γ (m) N m ǫ . ( 79
)
By using [START_REF] Sinatra | Limit of spin squeezing in finite-temperature Bose-Einstein condensates[END_REF] of [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF] and the Thomas-Fermi profiles of the condensate wave functions, one can optimize the squeezing with respect to time, trap frequency and atom number. This optimum squeezing in a trap in presence of one, two and three-body losses has a simple expression as a function of the s-wave scattering length a aa = a bb = a and the rate constants K m , m = 1, 2, 3:

ξ 2 opt = 5 √ 3 28π M a 2/3 7 2 (K 1 K 3 ) + K 2 2/3 . (80)
In Fig. 9 we show the spin squeezing minimized over time ξ 2 min as a function of N in presence of one, two and threebody losses. The trap frequency is optimized for each value of N [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF]:

ω opt = 2 19/12 7 5/12 π 5/6 15 1/3 m a 1/2 N 1/3 K 1 K 3 5/12 . ( 81
)
We note that the minimum squeezing, instead of going to zero as in two-mode model without decoherence (red line) tends to a finite non-zero value for N → ∞ given by (80). min as a function of N in presence of one, two and three-body losses, obtained in [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF] in the constant loss rate approximation. The trap frequency is optimized for each value of N according to (81). Dashed blue line: Optimum squeezing for large N given by the analytical prediction (80). Red line: Best squeezing without losses for comparison. Parameters: a = 5.32nm, K (1) = 0.1s -1 , K (2) = 2 × 10 -21 m 3 s -1 , K (3) = 18 × 10 -42 m 6 s -1 . This figure is taken from [18].

IV. FINITE TEMPERATURE

The multimode nature of the atomic field and the population in the excited modes at non-zero temperature have important consequences on the squeezing and change the scaling laws with respect to the two-mode case [START_REF] Yun | Spin squeezing in a bimodal condensate: spatial dynamics and particle losses[END_REF]. In [START_REF] Sinatra | Limit of spin squeezing in finite-temperature Bose-Einstein condensates[END_REF] we use a powerful formulation of the Bogoliubov theory in terms of the time dependent condensate phase operator [START_REF] Sinatra | Nondiffusive phase spreading of a Bose-Einstein condensate at finite temperature[END_REF][START_REF] Sinatra | Coherence time of a Bose-Einstein condensate[END_REF] to perform a multimode treatment of the squeezing generation in a condensed gas in the homogeneous case. We show that the best squeezing ξ min has a finite non-zero value in the thermodynamic limit and we calculate this value analytically.

A. Multimode description

We consider a discretized model on a lattice with unit cell of volume dV , within a volume V with periodic boundary conditions [START_REF] Sinatra | Coherence time of a Bose-Einstein condensate[END_REF]. The Hamiltonian after the pulse for component a (and similarly for b) reads:

H a = k 2 k 2 2m a † k a k + g 2 dV r ψ † a (r)ψ † a (r)ψ a (r)ψ a (r) .
(82) The fields have commutators

[ψ µ (r), ψ † ν (r ′ )] = δ rr ′ δ µν dV (83)
with µ, ν = a or b, and a k (b k ) is the amplitude of ψ a,b over the plane wave of momentum k. We assume identical interactions in states a and b with a coupling constant g = 4π 2 a/m where a is the s-wave scattering length while states a and b do not interact (g ab = 0). Note that the coupling constant in H a,b should actually be a bare coupling constant g 0 different from the effective coupling constant g, but this difference can be made small in the present weakly interacting regime (ρa 3 ) 1/2 ≪ 1 by choosing a lattice spacing much larger than a but still much smaller than the healing length ∝ 1/ √ ρa [START_REF] Mora | Extension of bogoliubov theory to quasicondensates[END_REF]. In terms of the fields, the collective spin components are

S x + iS y = dV r ψ † a (r )ψ b (r ), S z = N a -N b 2 (84)
with N ν = dV r ψ † ν (r)ψ ν (r), ν = a, b. For low or high values of k B T /ρg we find the asymptotic behaviors

k B T ≪ ρg : ξ 2 min ρa 3 ≃ ξ 2 (T =0) min ρa 3 = 0.0234 . . . (85) k B T ≫ ρg : ξ 2 min ≃ N nc N . (86)

B. Best squeezing and close-to-best time

Performing a double expansion for large N and small non-condensed fraction ǫ Bog ,

ǫ Bog ≡ N nc N ≪ 1 (87)
we find that to first order the system effectively behaves as in the dephasing model presented in section II. Indeed the component S y of the spin develops a term that is proportional to the condensate relative phase θ aθ b , and the relative phase evolves as

θ a -θ b = (θ a -θ b )(0 + ) - gt V [(N a -N b ) + D th ] . (88)
The dephasing parameter D th in (88) is related to the population in the excited modes:

D th = k =0 (U k + V k ) 2 (n ak -n bk ) , (89)
where n ak = c † ak c ak and n bk = c † bk c bk are the occupation number operators of the Bogoliubov modes. Note that these modes are in a non-equilibrium state since the zero relative phase state between the condensates in a and b is prepared at t = 0 + by applying a sudden π/2 pulse to the gas (that was initially at thermal equilibrium in state a). U k and V k are the usual Bogoliubov functions

U k + V k = E k E k + ρg 1/4 ; E k = 2 k 2 2m , (90)
ρ/2 = N/2V being the spatial density in each single component a or b after the pulse. D th fluctuates from one realization to the other because n ak and n bk depend (in Heisenberg picture) on the creation and annihilation operators of the Bogoliubov modes before the pulse, and the initial state of the gas has thermal fluctuations. To first order in N and in ǫ Bog , the best squeezing parameter and the close-to-best squeezing time are given by [START_REF] Sinatra | Limit of spin squeezing in finite-temperature Bose-Einstein condensates[END_REF]

ξ 2 min = D 2 th N ; ρgt η = 1 √ ηξ min . (91)
An explicit calculation [START_REF] Sinatra | Limit of spin squeezing in finite-temperature Bose-Einstein condensates[END_REF] gives:

ξ 2 min = d 3 k (2π) 3 s 4 k 2ρ n (0) k + 1 2 (s (0) k ) 2 s 4 k + s 4 k (s (0) k ) 2 -1 , (92)
Particle Losses Dephasing model Multimode T = 0

|ψ(t) = N 1/2 e i φ|φ + χt 2 D (θa -θ b )(t) = (θa -θ b )(0 + ) -χt [2Sz + D] (θa -θ b )(t) = (θa -θ b )(0 + ) -χt [2Sz + D th ]
D from quantum jumps (77) D from a dephasing Hamiltonian (15) D th from excited modes population (89) where now n

ξ 2 (t) ≃ ρgt/ >1 D 2 N ξ 2 (t) ≃ ρgt/ >1 D 2 N ξ 2 (t) ≃ ρgt/ >1 D 2 th N D 2 N = γt 3 = ǫ loss 3 D 2 N = ǫnoise D 2 th N = ρa 2 F (kBT /ρg) ∼ k B T >ρg ǫBog
U (0) k + V (0) k = E k E k + 2ρg 1/4 ; E k = 2 k 2 2m . (94)
A consequence of (92) is that the squeezing divided by ρa 3 is an universal function of k B T /ρg.

C. Spin squeezing and non-condensed fraction

Within our treatment, valid for large N and T ≪ T c , we find that the best squeezing ξ 2 min (92) is always lower than the before-pulse non-condensed fraction. This is shown in Fig. 10.

Finally in Table IV A we can complete the correspondence table between dephasing noise, losses and non-zero temperature effects. This gives the possibility to tune the non-linearity [START_REF] Yun | Spin squeezing in a bimodal condensate: spatial dynamics and particle losses[END_REF] by either controlling the spatial overlap between the two species as done [START_REF] Riedel | Atom-chip-based generation of entanglement for quantum metrology[END_REF] or by using a Feshbach resonance changing the inter-species coupling constant g ab as done in [START_REF] Gross | Nonlinear atom interferometer surpasses classical precision limit[END_REF]. A fundamental source of decoherence in these experiments is two-body losses in the state F = 2. A theoretical analysis [START_REF] Yun | Spin squeezing in a bimodal condensate: spatial dynamics and particle losses[END_REF] shows that in typical experimental conditions these losses limit the squeezing to about ξ 2 min ≃ 6 × 10 -2 . For cold samples this limit is above the limit imposed by non-zero temperature. In the case of [START_REF] Riedel | Atom-chip-based generation of entanglement for quantum metrology[END_REF] we could explain in detail the squeezing results using a zero-temperature model including spatial dynamics and including particle losses and technical noise (dephasing noise) as sources of decoherence, the latter being dominant [START_REF] Riedel | Atom-chip-based generation of entanglement for quantum metrology[END_REF]18].

The perspective of studying experimentally the scaling of squeezing in a controlled decoherence environment, e.g. preparing the sample at different temperatures, is fascinating and challenging.

VI. CONCLUSIONS

We have considered a scheme to create spin squeezing using interactions in Bose-condensed gas with two internal states [START_REF] Kitagawa | Squeezed spin states[END_REF][START_REF] Sorensen | Manyparticle entanglement with Bose-Einstein condensates[END_REF][START_REF] Uffe | Positive-P simulations of spin squeezing in a two-component Bose condensate[END_REF]. The squeezing is created dynamically after a π/2 pulse applied on the system initially at equilibrium in one internal state. We have reviewed the ultimate limits of this squeezing scheme imposed by particle losses and non-zero temperature based on our recent works [START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF] and [START_REF] Sinatra | Limit of spin squeezing in finite-temperature Bose-Einstein condensates[END_REF] and we have extracted a simple physical picture of how decoherence acts in the system. An important result is that contrarily to the case without decoherence [START_REF] Kitagawa | Squeezed spin states[END_REF] the squeezing parameter minimized over time ξ 2 min has a finite non-zero value in the thermody-namic limit, that we determine analytically. Finally we have shown that the physics of spin squeezing in presence of losses or at non-zero temperature can be caught by a simple dephasing model (also considered in [START_REF] Ferrini | Effect of phase noise on useful quantum correlations in Bose Josephson junctions[END_REF]) that we have solved exactly and studied in details.

FIG. 1 :

 1 FIG. 1: (Color online) Uncorrelated state and squeezed state represented on the Bloch sphere.

FIG. 2 :

 2 FIG. 2: (Color online) Left: Deformation of quantum fluctuations due to the non-linear Hamiltonian (8). Right: Best squeezing time.

2 FIG. 3 :

 23 FIG. 3: (Color online) Left: Squeezing parameter ξ 2 as a function of time for ǫnoise = 0.05 (full line) and ǫnoise = 0 (dotted line). Determinant of the covariance matrix M as a function of time (dashed line). Simulation with 10 4 realizations. N = 4 × 10 4 , χ = 0.022577s -1 . Right: Isocontours of W (y, z) for t = 0, t = 2 × 10 -3 s and t = 2.5 × 10 -2 s.

2 FIG. 4 :

 24 FIG.4: Squeezing as a function of time for the dephasing model[START_REF] Wright | Collapses and revivals in the interference between two Bose-Einstein condensates formed in small atomic samples[END_REF] as given by the exact solution (47). N = 10 9 , D 2 /N = 0.02. Horizontal dashed lines: ξ 2 min given by (A4) and (1 + η)ξ 2 min with η = 0.1. Vertical dashed lines from left to right: Close-to-best squeezing time tη (53), best squeezing time tmin (A3) and t ′ η deduced from (A5).

FIG. 6 :

 6 FIG.6: Relative phase probability distribution at time t = 0 and at time t = 2π/χ in three single realizations Monte Carlo realization. From upper left to lower right the wave function has experienced 0, 3, 1 and 0 quantum jumps respectively. Three-body losses corresponding to N = 301 87 Rb atoms in F = 1, m f = -1 state in separated identical harmonic traps with ω = 500Hz. Figure taken from[START_REF] Sinatra | Phase dynamics of Bose-Einstein condensates: Losses versus revivals[END_REF].

2 FIG. 7 :

 27 FIG. 7: (Color online) Spin squeezing as a function of time in presence of one-body losses: γ = 0.005ρg/ , N = 10 6 . Solid black line: Exact analytical solution. Dashed red line: Approximate solution (72). Horizontal and vertical dashed lines: predictions (73) and (74) respectively.

FIG. 8 :

 8 FIG. 8: (Color online) Main plot: D 2 /N given by (77) as a function of time for γ = 0.005ρg/ , N = 10 6 . Black circles: Numerical simulation with 5 × 10 5 realizations. Dashed blue line: Approximate expression (78). Solid blue line: Exact expression (B7). Inset: Squeezing as a function of time. Solid red line: Exact analytical solution. Black circles: Numerical simulation with 5 × 10 5 realizations.

2 FIG. 9 :

 29 FIG.9:(Color online) Blue line: Best squeezing ξ 2 min as a function of N in presence of one, two and three-body losses, obtained in[START_REF] Li | Optimum spinsqueezing in Bose-Einstein condensates with particle losses[END_REF] in the constant loss rate approximation. The trap frequency is optimized for each value of N according to (81). Dashed blue line: Optimum squeezing for large N given by the analytical prediction (80). Red line: Best squeezing without losses for comparison. Parameters: a = 5.32nm,K (1) = 0.1s -1 , K (2) = 2 × 10 -21 m 3 s -1 , K(3) = 18 × 10 -42 m 6 s -1 . This figure is taken from[18].

2 FIG. 10 :

 210 FIG.10:(Color online) ξ 2 min (solid line) and before-pulse noncondensed fraction Nnc /N (dashed line), both divided by ρa 3 , as functions of kBT /ρg. Both quantities are obtained analytically within a Bogoliubov framework, see (92) for ξ 2 min .

are

  Bose mean occupation numbers of Bogoliubov modes in state a before the pulsen (0) k = 1 e ǫ k /kB T -1 , ǫ k = E k (E k + 2ρg) . (93)s k = U k + V kand s functions in internal state a before the pulse

V

  . A FEW WORDS ABOUT EXPERIMENTS Two recent experiments demonstrated spin squeezing in Bose-condensed bimodal condensates of rubidium atoms in internal states |a = |F = 1, m F = ∓1 and |b = |F = 2, m F = ±1 [4, 6]. Due to the fact that the three scattering lengths characterizing the interactions between atoms are very close g aa ≃ g ab ≃ g bb (95) the effective two-mode nonlinearity χ in H nl (8) is very small when the condensates a and b overlap spatially.

  . 5: Time sequence of deterministic evolution periods and quantum jumps for a single Monte Carlo wave function |ψ .

	τ 1 τ 2 ...	τ k t-t k
	0 t 1 t 2	...	t k-1 t k	t
	64)

FIG

TABLE I :

 I Correspondence table among the different physical models. The results are valid in the thermodynamic limit and to first order in the small parameters ǫ loss (69), ǫnoise (16) and ǫBog (87).

Appendix A: Times tmin and t ′ η in the dephasing model

The exact solution (47)-(50) allows to determine how the best squeezing time t min diverges in the thermodynamic limit. We first found numerically that it diverges as N 1/4 . We then introduce the rescaled time θ such that ρg t = θN

By minimizing (A2) over θ one obtains ρg t min lim.therm.

With a similar technique we can determine the divergence of t ′ η that is the second solution of equation (37) t ′ η > t min . When N → ∞, t ′ η diverge as N 1/2 . To calculate the prefactor we introduce again a rescaled time ρgt/ = θ ′ [N/(1 + ǫ noise )] 1/2 and take the large N limit in (47) to obtain

Solving the transcendental equation (1 + η)ξ 2 = ξ 2 min one finds the large N approximation to t ′ η . For η ≪ 1 and no constraint on the ratio η/ǫ noise , we obtain

Appendix B: Calculation of D 2 in the lossy model

From the definition of D (77) and the expression of a quantum average (65) in the Monte Carlo wavefunction method, using the expression of the norm squared of |ψ(t) of (76), we obtain

We have introduced the random variables η i = +1 for ǫ i = b and η i = -1 for ǫ i = a, and we have used the that the integrand is a symmetric function of the jump times t i to extend time integration from the ordered domain 0 < t 1 < . . . < t k to the hypercube [0, t] k (also dividing by k!). The notation N k represents the usual binomial coefficient N !/[k!(Nk)!]. We first sum over the variables η i :

then we perform the temporal integration to obtain

Taking the derivative with respect to u of the binomial identity k i=0 N k u k = (1 + u) N , we get the final expression

Expanding I 2 for small γt gives as expected


Ruth Sayer

Camden Pct

Mortimer Market

Clinic Miller

Trevor ; Baglin

Nhs Addenbrookes

Haematology Trust

Beverley ; Hunt

Gstt

Sylvia ; Benjamin

Raj ; Patel

D P Hart

R F Miller

S G Edwards

A Kelly

B J Hunt

S J Machin

Dr Marie Scully

Human immunodeficiency virus associated thrombotic thrombocytopenic purpura (TTP) -favourable outcome with plasma exchange and prompt initiation of Highly Active Antiretroviral Therapy (HAART)

Keywords: TTP, HIV, HAART TTP, HIV, ADAMTS 13, HAART, plasma exchange

Background: Thrombotic thrombocytopenic Purpura is an acute prothrombotic disorder. HIV is an identified precipitant. Design and methods Thirty episodes of HIV associated TTP in 24 patients from the South-East England Apheresis units, over past 10 years, are reviewed. Results: All patients were heterosexual Black Africans. First presentation of TTP revealed a new diagnosis of HIV in 8 patients. TTP relapse occurred on 6 occasions (in 4 patients) as a result of non-adherence to HAART. Prompt initiation/re-initiation of HAART in parallel with PEX +/-steroid led to prompt remission. Adjunct immunomodulatory agents (e.g.Rituximab) was required in 10% of cases. Once-daily HAART regimens are recommended, being compatible with PEX requirement, maximising drug exposure between PEX. High viral loads (>500,000copies/ml) require more PEX to remission. ADAMTS 13 activity was reduced (<5%) as detected by CBA and Anti-ADAMTS 13 IgG antibodies raised in 80%. Continued HAART-adherence ensured a durable TTP remission with associated viral control resulted in no evidence of relapse. Conclusion: PEX and HAART are associated with replenishment of ADAMTS 13 and viral suppression.

More PEX is required in cases with higher viral loads. Continued HAART maintains remission. In a small proportion of cases, further immunomodulatory therapy may be required.

INTRODUCTION

Thrombotic thrombocytopenic purpura (TTP) results from a deficiency of von Willebrand factor-cleaving protease (vWF-cp) activity. ADAMTS-13 (a disintegrin and metalloprotease with thrombospondin type 1 motif 13) is the enzyme identified as mediating this cleaving function [START_REF] Levy | Mutations in a member of the ADAMTS gene family cause thrombotic thrombocytopenic purpura[END_REF]. The loss of cleaving activity results in circulation of abnormal ultra large vWF multimers (ULvWF) provoking widespread microvascular thrombosis. Urgent plasma exchange (PEX) and adjunctive steroid therapy is the cornerstone of treatment for acquired TTP and has resulted in dramatic improvements in survival, from 10% in the pre-PEX era to 80-90% with current PEX-based regimens [START_REF] Rock | Comparison of plasma exchange with plasma infusion in the treatment of thrombotic thrombocytopenic purpura. Canadian Apheresis Study Group[END_REF]. Further improvements in clinical outcome of F o r P e e r R e v i e w 3 idiopathic-TTP have been reported using anti-CD20 antibody (Rituximab) to reduce relapse risk in antibody mediated idiopathic disease, accounting for the majority of cases [START_REF] Scully | Remission in acute refractory and relapsing thrombotic thrombocytopenic purpura following rituximab is associated with a reduction in IgG antibodies to ADAMTS-13[END_REF]. Simultaneous in-vitro assay development has improved quantitation of ADAMTS13 activity and detection of anti-ADAMTS13 antibodies [START_REF] Peyvandi | ADAMTS-13 assays in thrombotic thrombocytopenic purpura[END_REF]. Human immunodeficiency virus (HIV) associated TTP is a subgroup (5) defined as secondary TTP, as there is a precipitant identified. In this paper we describe patients with HIV-associated TTP registered in SE England TTP registry and subsequently the UK registry (but within the SE England), over a 10 year period.

PATIENTS, MATERIALS AND METHODS

Twenty four HIV positive patients, eighteen of which were female, mean age 36 years (range 1 -64) were identified from five South East England units offering apheresis between 1998 and 2009. We have included eight patients previously published [START_REF] Miller | Thrombotic thrombocytopaenic purpura in HIV-infected patients[END_REF]. A total of 30 episodes of acute TTP were identified from this population. Four patients had relapses (three experiencing a single relapse and the remaining patient had three relapses). Clinical case notes were accessed when necessary to complete data sets.

Diagnosis of TTP was made using BCSH guidelines [START_REF] Allford | Guidelines on the diagnosis and management of the thrombotic microangiopathic haemolytic anaemias[END_REF], on clinical grounds, including thrombocytopenia and microangiopathic anaemia characterised by blood film red cell fragmentation, reticulocytosis, normal coagulation screen, negative direct antiglobulin test (DAT) and elevated lactate dehydrogenase (LDH), supported by ADAMTS 13 assays. On diagnosis of TTP all patients received PEX and steroids (either pulsed methylprednisolone or high dose oral steroids) as per local protocol.

Only four patients did not receive steroids on admission, typically if they were already known to have HIV infection.

Since 2004, ADAMTS13 activity has been undertaken by a modification of collagen binding assay (CBA) [START_REF] Gerritsen | Assay of von Willebrand factor (vWF)-cleaving protease based on decreased collagen binding affinity of degraded vWF: a tool for the diagnosis of thrombotic thrombocytopenic purpura (TTP)[END_REF], normal 66-126%, and before 2004, by multimer analysis [START_REF] Furlan | Deficient activity of von Willebrand factor-cleaving protease in chronic relapsing thrombotic thrombocytopenic purpura[END_REF] (normal range 80 -120%. Presenting haematological parameters revealed universal thrombocytopenic (mean 20 x 10 9 /L, range 5-107x10 9 /L ) a raised lactate dehydrogenase (LDH) (mean 2,285, range 464-8,000 IU/L) and variable levels of anaemia (mean presenting haemoglobin 6.8 g/dL (range 3 -13.4g/dl). Nineteen cases presented with neurological symptoms ranging from headache and confusion to seizures and localising neurology Table 1. Ten patients had renal impairment at presentation, correcting to within the normal range with treatment. There was no deterioration in renal function on starting HAART. Four cases had cardiac symptoms/raised Troponin at presentation. Daily plasma exchange and concurrent initiation of HAART resulted in an average of 13.5 plasma exchange (PEX) (range 3 -33), to remission, when calculated for the whole group.

In seven cases, despite previous anti-retroviral treatment, there was an increase in viral load associated with the initial presentation of TTP. One was receiving suboptimal single agent antiretroviral therapy and three episodes were thought to be adherent with HAART but had HIV resistance. One other case had demonstrable HIV-resistance mutations and rising viral load despite multiple HAART regimens.

Rising viral load attributed to resistance was observed in a second case and no data was available for the third. Patients with HIV resistance mutations had not previously been reported as presenting with TTP at HIV diagnosis. HIV viral load was >10,000 copies/ml at TTP diagnosis in nineteen cases: three had viral loads >500,000 copies/ml; five between 100,000-500,000 copies/ml and nine between 10,000-100,000 copies/ml. Two cases presented with viral loads <10,000 copies/ml. Cases (n=10) with viral loads between 10,000-500,000 copies/ml required fewer PEX to remission (median 7.5 (range: 6-25). In those with viral loads >500,000 (n=3), the median PEX to remission was 27 [START_REF] Novitzky | Thrombotic thrombocytopenic purpura in patients with retroviral infection is highly responsive to plasma infusion therapy[END_REF][START_REF] Park | ADAMTS13 activity levels in patients with human immunodeficiency virus-associated thrombotic microangiopathy and profound CD4 deficiency[END_REF](14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30). The number of cases were too small to Only one relapsed case and one de novo TTP required additional therapy; Rituximab. In these cases, despite PEX, steroids and HAART with a reduction in viral load, there was prolonged thrombocytopenia and ADAMTS 13 activity remained <5%. One case received 2 infusions and the other 4 infusions (375mg/m 2 ), given weekly. This resulted in prompt sustained remission, with no adverse effects (specifically no infections). There was one death in this cohort, an atypical presentation for HIV-associated TTP. The patient was already receiving HAART; viral load was undetectable and remission induction and maintenance was difficult despite PEX, steroids, Rituximab and continued HAART. The cause of death was not determined. Although not formally tested, viral resistance and/or HAART noncompliance was unlikely as the viral load was undetectable. A further case had undetectable viral load, low ADAMTS 13 activity (<5%) and a raised Anti-ADAMTS 13 IgG (105%). They were already receiving HAART which was not altered on admission and responded to PEX and high dose steroids.

Patients who continued HAART associated with lowered, mainly undetectable viral loads, have not relapsed. However, relapse was associated with stopping HAART and re-presentation of TTP within 3 months. In twenty-four cases, pre-PEX ADAMTS 13 activity data was available -in 16 patients was <5% by CBA and was <50% in the remaining cases by the multimer assay. In eighteen cases, anti-ADAMTS 13 IgG levels were available, median 24% (4-105%). Remission was associated with normalisation of ADAMTS 13 activity; relapse resulted in a further decreased ADAMTS 13 activity to <5%.

DISCUSSION

Survival of HIV-associated TTP in the pre-HAART era was rarely more than 2 years despite improvements in acute TTP management with PEX and steroid treatment [START_REF] Gadallah | Disparate prognosis of thrombotic microangiopathy in HIV-infected patients with and without AIDS[END_REF]. Our data confirms prompt, concurrent treatment with both PEX and HAART results in a favourable outcome in most cases, with a low mortality of 4% of our cohort. This is in contrast to acquired, idiopathic TTP, where mortality is 10-20% (2;10). Indeed, when the patients are adherent with HAART therapy, the apparent relapse rate in this cohort is lower than historical relapse rates seen in idiopathic TTP treated with PEX/steroid alone. TTP relapses were seen in patients with virological rebound attributed primarily to non adherence with HIV therapy. One case was repeatedly non-compliant and as a consequence has relapsed on multiple occasions. One of these non-compliant relapses required adjunctive Rituximab to reinduce a remission. Only one patient appeared refractory to PEX/steroids/HAART at first presentation and consequently received Rituximab to achieve remission. Indeed, Rituximab has been used successfully in HIV-TTP [START_REF] Ross | HIV, thrombotic thrombocytopaenic purpura and rituximab in a violent noncompliant patient[END_REF]. Rituximab therapy has impacted significantly on refractory and relapsing TTP in non-HIV TTP (3). However, in the context of HIV associated TTP, we believe immunomodulatory therapy such as Rituximab, should only be considered in the unlikely scenario of primary resistance to PEX/steroid/HAART and after HAART non-adherence or HIV resistance have been excluded. There appears to be no association with the Anti-ADAMTS 13 IgG titre and the clinical severity of HIV-associated TTP, which is in contrast to idiopathic TTP, where antibody associated TTP has a more aggressive clinical presentation [START_REF] Scully | Regional UK TTP registry: correlation with laboratory ADAMTS 13 analysis and clinical features[END_REF]. Therefore, the immune component in the pathophysiology of this subgroup appears generally less relevant than the direct effect of the retrovirus. Why these specific patients develop TTP rather than other HIV associated complications, either as a presenting event or with resistant disease remains unanswered.

ADAMTS 13 assays are useful in helping to confirm a diagnosis of TTP. CBA and more recently FRETS assay are more sensitive and specific, compared to multimeric analysis. All patients achieving remission had a normal ADAMTS 13 activity that was maintained while patients continued HAART.

Fifteen cases in this cohort were known to be HIV positive and already on HAART, and all but one had a change made to their HAART regimen upon presenting with TTP. Indeed optimisation of HIV therapy with new drugs was required to achieve virological suppression and utilising combination therapy advance remission of TTP.

Of interest, those with a viral load between 10,000 and <500,000 copies/ml at presentation, regardless of a new presentation or previously diagnosed HIV, required fewer PEX to remission than those with higher viral loads (>500,000 copies/ml).

Cases treated promptly with PEX/steroid/HAART therapy achieved a robust remission if they continued HAART and achieved suppression of the viral load. A single patient in this cohort achieved prompt TTP resolution with steroid and plasma infusion alone, prior to transfer for plasma exchange. Plasma infusion has been reported as being sufficient to achieve remission either in the absence of available PEX facilities [START_REF] Novitzky | Thrombotic thrombocytopenic purpura in patients with retroviral infection is highly responsive to plasma infusion therapy[END_REF] or in the context of patient non-compliance with PEX [START_REF] Ross | HIV, thrombotic thrombocytopaenic purpura and rituximab in a violent noncompliant patient[END_REF]. However, administration of such large volumes of plasma causes considerable clinical problems. Plasma infusion should be used with steroid in the event of delay transferring to a PEX facility and plasma exchange should remain the gold standard therapy for HIV associated TTP when logistically possible.

The unwell, HIV-positive patient presenting with thrombocytopenia demands consideration of a broad differential diagnosis list (13;14). However, TTP remains one of the most important clinical diagnoses requiring immediate treatment rather than deferred treatment pending confirmatory testing or exclusion of other diagnoses. The cohort we have presented had no evidence of other HIV associated TMA, such as HUS or opportunistic infection. We did not see, as recently reported, increased mortality rates in HIV-TTP (15). ADAMTS 13 activity in our cohort could not be used as a predictor of outcome and only one of our cases had other AIDS defining complications. Patients with advanced HIV disease may have thrombocytopenia and red cell fragmentation alone, rather than displaying other thrombotic phenomena, in keeping with HIV thrombotic microangiopathy, rather than HIV associated TTP. ( 16).

The cohort consists of HIV cases presenting over a period of a decade, in which time there have been significant changes in anti-retroviral preparations and combinations.

Although no conclusions can be drawn as to the optimal HIV treatment regimen required to expedite remission, safe combinations enabling once daily administration are ideal. We believe that HIV treatment should be administered immediately with the aim to suppress circulating HIV virus as quickly as possible. As an HIV resistance result may not be available at the time of starting treatment, a once daily boosted protease inhibitor (which has a high genetic barrier to resistance) should be included in the regimen. In addition, HIV therapy including raltegravir has been shown to result in a more rapid decline in HIV viral load that standard triple combination therapy. HAART should be administered immediately after PEX, As an autoimmune phenomenon, women are more commonly affected by idiopathic TTP than men [START_REF] Scully | Regional UK TTP registry: correlation with laboratory ADAMTS 13 analysis and clinical features[END_REF]. This gender skew (70% female) is also evident in this HIVassociated TTP cohort, relative to 58% of new HIV diagnoses being women.

However, the striking skew to patients of African origin in the HIV-TTP setting (100%), compared to only 34.7% of new HIV diagnoses in the UK being of African origin, is more difficult to explain, but may suggest a possible underlying genetic predisposition.(17).

In conclusion, thrombocytopenia in HIV is relatively common and potential causes and underlying mechanisms may vary significantly, However, HIV-associated TTP remains an important diagnosis, particularly in heterosexual, mainly female Black African patients. We emphasise the importance of prompt identification of the HIVpositive status and urgent initiation of HAART in parallel to PEX. Those with the highest viral loads at presentation (>500,000 copies/ml) required more PEX to remission than those with viral loads <500,000 copies/ml. If TTP is diagnosed in an HIV positive patient already on HAART, adherence and/or resistance studies should be considered. Rarely, further immunomodulatory therapy such as Rituximab may be needed. A patient's adherence with HAART ensuring continued viral suppression prevents further TTP relapse.

F

4 RESULTS

 4 Inhibitory activity was determined by 1:1 mixing studies with pooled normal plasma for 1 hour incubation at 37 0 C. Anti ADAMTS13 IgG detected by a modified ELISA (3) Thirty episodes of HIV-related acute TTP were treated over the past 10 years. Twenty four cases, (seventeen women) represented a first episode of TTP, the remaining 6 episodes being TTP relapses. All patients were heterosexual Black Africans. In eight episodes, presentation with TTP resulted in a new diagnosis of HIV infection. The remaining twenty one episodes occurred in patients who were known to have HIV infection.

  statistical significance. Average CD4 count was 236 x 10 6 /l (n=24) (range 70 -540) on admission.

  and minimising chance of drug removal with subsequent PEX.

  Benjamin M, Terrell DR, Vesely SK, Voskuhl GW, Dezube BJ, Kremer Hovinga JA, et al. Frequency and significance of HIV infection among patients diagnosed with thrombotic thrombocytopenic purpura. Clin Infect Dis 2009 Apr 15;48(8):1129-37.

	17 (14) (15) Malak S, Wolf M, Millot GA, Mariotte E, Veyradier A, Meynard JL, et al. Human 42 F Black African N 14 Black 18 30 M African Y 10	2874 1500	6.8 5.5	380,000 120,000	170 240	c,b n,a,	9 9.5	da,r,ki da,r,tr
	immunodeficiency virus-associated thrombotic microangiopathies: clinical Platelets Black			Viral load			20 +
	o r P e F o r P e characteristics and outcome according to ADAMTS13 activity. Scand J Immunol age at New on 19 64 M African N 14 2008 Sep;68(3):337-44. (16) Becker S, Fusco G, Fusco J, Balu R, Gangjee S, Brennan C, et al. HIV-associated thrombotic microangiopathy in the era of highly active antiretroviral therapy: an observational study. Clin Infect Dis 2004 Nov 1;39 Suppl 5:S267-S275. (17) HPA: HIV in the United Kingdom: 2008 Report. http://www.hpa.org.uk/web/HPAwebFile/HPAweb_C/1227515298354 F o r P e e r R e v i e w LDH on Hb on 1668 6.9 Formatted: Font: (Default) Times <40 on New Roman, 12 pt Field Code Changed Formatted: Font: (Default) Times New Roman Formatted: Font: (Default) Times New Roman Patient presentation (years) gender Ethnicity diagnosis HIV admission (x10 admission admission admission Black 9 /L) (IU)(NR<480IU) (g/d L) (Copies/ml) CD4 Symptoms 393 n 1 30 F Black African Y 8 3078 9.4 854,000 70 n,f 2 27 M Black African Y 10 3317 8.9 n/a n/a f,r 3a 23 M Black African N 7 2390 6.7 67,800 170 n 4 45 F Black African Y 15 1967 8.1 19,500 400 n,b 5a 44 F Black African N 16 914 4.6 47,100 170 n,r 6 35 F Black African n/a 46 n/a 6.3 n/a n/a n 7 1 M Black African N 20 >8000 6.5 >500,000 260 n 8 40 F Black African Y 11 5009 6.0 51,900 130 n,a,r,c 9a 37 F Black African N 14 2196 8.5 n/a 190 n,f 10 38 F Black African N 9 483 7.8 48,000 490 n 20 25 M African Y 79 1199 10.5 79,091 540 n,r 21 35 F Black African N 18 1954 7.4 n/a 235 ? 22 29 F Black African N 27 1335 7.3 n/a n/a n,f 23 37 M Black African N 11 576 5.4 n/a n/a n,f,r Black F o r P e e r R e v i e w 24 42 F African N 6 1438 5.4 <50 74 r,c	Ritux No.of PEX 12 20 16 16 7 13 30 8 5 6 18 10 9 11 20	#t,d,k Commenced HAART n e,c l,lm,a dd,tc,k l,lm,d k,az z,lm,k l,tc,dd l *t,d4,tc,k,at n k,r,t N/A *nv, d4 N/A
	11 12	35 32	F F	e r Black e African Y Black r African N	12 22	975 3500	5.4 6.6	n/a 467,000	320 50	n,f n,r	11 5	l,lm,d z,lm,nv
	13 14 15a 16	37 41 39 33	F F F F	Black African Black African Black African Black african	R e v i e N 53 R Y 45 e N 9 v i e N 8	1234 1673 2124 3139	6.8 5.5 5.3 5.1	n/a 12,617 >500,000 230,000	209 250 200 420	n,f,c N/A n tired	9 6 27 25 + Ritux	*s,lm,nl s,lm,nl tr,l da,r,tr
						w w
				10

F o r P e e r R e v i e w

Table 1a : HIV-TTP Presenting features and treatment of first presentation Key:

 1a

	f:fever
	r:renal	#existing combination	d:didanosine
	a:abdo	* change of existing combination	az: azt
	n:neuro	n: commenced but details unavailable z: zidivudine
	c:cardiac	e: efavirenz	t:tenofovir
	b:bleeding	c: combivir	d4: d4t
		l: lopinavir	at: atazanavir
		lm: lamivudine	nv: nevirapine
		a: abacavir	s: stavudine
		dd:ddi	nl: nelfinavir
		tc: 3TC	tr:truvada
		k: kaletra	da:darunavir
			r: ritonavir
			ki:kivexa

British Journal of Haematology

F o r P e e r R e v i e w


Keywords: ParFlow, TLM, propagation des ondes, propagation indoor, simulation de canal, multi-résolution, domaine fréquentiel, OFDM

Optimisation de la méthode MR-FDPF pour simulation à multiple fréquences Résumé : Dans ce rapport nous proposons une nouvelle implémentation de la méthode ParFlow multi-résolution (MR-FDPF) permettant de simuler le canal de propagation radio dans une large bande de fréquences. Les performances de la méthode ParFlow implémentée ont été validées dans différents scenarios.

Introduction

The new generation of wireless communication systems, such as the IEEE 802.16 WiMAX and the 3GPP Long Term Evolution (LTE), has recently come to the stage of actual deployment. Broadband wireless access (BWA) technologies constitute indispensable part of these systems. The orthogonal frequency division multiplexing (OFDM) [START_REF] Jha | OFDM. Towards Fixed and Mobile Broadband Wireless Access[END_REF] is one of the key BWA techniques implemented in modern high data rate wireless communication systems.

In order to evaluate and optimize the performance of OFDM systems, reliable models that adequately describe radio propagation channels over the whole system frequency bandwidth are required. The deterministic channel models based on the numerical simulation of the Maxwell equations [START_REF] Daniel | Microwave Circuit Modeling Using Electromagnetic Field Simulation[END_REF], as well as the models based on the ray-optical methods [START_REF] Hoppe | Advanced ray-optical wave propagation modelling for urban and indoor scenarious including wideband properties[END_REF], are known for their relatively high accuracy [START_REF] Molisch | Wireless Communications[END_REF]. One of the main drawbacks of these modeling approaches is the large computational load.

During the last decade, a great deal of research has been focused on finding an optimal trade-off between the computational cost and the achievable accuracy of the deterministic channel modeling approaches. In [START_REF] Gorce | The Adaptive Multi-Resolution Frequency-Domain ParFlow (MR-FDPF) Method for Indoor Radio Wave Propagation Simulation. Part I : Theory and Algorithms[END_REF], the multi-resolution frequency domain ParFlow (MR-FDPF) method has been proposed. The method has proved to be efficient for modeling indoor and indoor-like environments. However, similar to the majority of the frequency domain methods [START_REF] Berini | A new frequency domain symmetrical condensed tlm node[END_REF][START_REF] Johns | New frequency-domain tlm method for the numerical solution of steady-state electromagnetic problems[END_REF][START_REF] Jin | A new Frequency-Domain TLM symmetrical condensed node derived directly from Maxwell's equations[END_REF], the MR-FDPF algorithm has to be repeated for every frequency in the range of interest. Consequently, the computational cost associated with applying the method quickly becomes excessively large as the number of frequencies increases. This is particularly the case for OFDM wireless communication systems, where data is transmitted in parallel over several sub-carriers.

further development of the MR-FDPF method. The proposed modification of the MR-FDPF method allows computationally efficient radio coverage prediction at multiple frequencies. In the traditional implementation of the MR-FDPF method (see, e.g., [START_REF] Gorce | The Adaptive Multi-Resolution Frequency-Domain ParFlow (MR-FDPF) Method for Indoor Radio Wave Propagation Simulation. Part I : Theory and Algorithms[END_REF]), simulations of the electromagnetic wave propagation at every single frequency consists of three steps: the preprocessing step, the upward step, and the downward step. The preprocessing step is by far the most demanding step of the MR-FDPF method in terms of the computational resources. The computational efficiency of the modified MR-FDPF method is achieved by conducting the preprocessing step only once during simulations of the wave propagation in a range of frequencies.

The rest of the report is organized as follows. Section 2, gives a general overview of the MR-FDPF method. In Section 3, we describe the proposed modification of the MR-FDPF method for channel modeling in OFDM systems. The results of the performance analysis for the proposed implementation of the MR-FDPF method are presented in Section 4. The concluding remarks are given in Section ??.

MR-FDPF Method

In this section, a general overview of the MR-FDPF method is provided in order to facilitate the presentation of the material in Section 3. For additional details, an interested reader is addressed to [START_REF] Gorce | The Adaptive Multi-Resolution Frequency-Domain ParFlow (MR-FDPF) Method for Indoor Radio Wave Propagation Simulation. Part I : Theory and Algorithms[END_REF][START_REF] Gorce | Deterministic approach for fast simulations of indoor radio wave propagation[END_REF][START_REF] De La Roche | On predicting in-building WiFi coverage with a fast discrete approach[END_REF].

ParFlow Formulation

The ParFlow (or, more exactly, the time domain ParFlow) method is based on the cellular automata formalism [START_REF] Chopard | Lattice-Boltzmann method for wave propagation in urban microcells[END_REF][START_REF] Chopard | Cellular Automata Modeling of Physical Systems[END_REF][START_REF] Luthi | Lattice wave automata: From radio waves to fractures propagation[END_REF][START_REF] Chopard | Lattice-Boltzmann computations and applications to physics[END_REF]. In this method, the scalar electrical field strength is obtained by summing the fictitious flows traveling along a regular grid of connection lines and experiencing scattering at the nodes of the grid (see Fig. 1).

At every node r of the grid, the discrete time evolution of the flows is described by the equation

-→ F (r, t) = W(r) ← - F (r, t -△t) + -→ S (r, t) (1)
where the incident (inward) flows vector ← -F (r, t), the reflected (outward) flows vector -→ F (r, t), and the source flows vector -→ S (r, t) are defined, respectively, as follows

← - F (r, t) = [ ← - f E (r, t) ← - f W (r, t) ← - f S (r, t) ← - f N (r, t) f0 (r, t) ] T -→ F (r, t) = [ -→ f E (r, t) -→ f W (r, t) -→ f S (r, t) -→ f N (r, t) f0 (r, t) ] T -→ S (r, t) = [ -→ s E (r, t) -→ s W (r, t) -→ s S (r, t) -→ s N (r, t) 0] T . ( 2
)
The operator {•} T signifies matrix transposition. The inner flow f0 (r, t) in (2) models dielectric media with the relative permittivities ε r ̸ = 1. Note that the reflected flows at the grid node r correspond to the incident flows at its immediate neighbor nodes, that is ← -

f d (r + △r, t) = -→ f d (r, t -△t), d ∈ E, W, S, N (3)
where △r is the spatial grid step, which is assumed to be the same along all the coordinate axes. The step △r satisfies the following conditions

△r = c 0 √ 2△t, △r ≪ λ ( 4
)
where c 0 denotes the speed of an electromagnetic wave in a propagation medium and λ is the smallest wavelength in the spectrum of the source signal. The scattering matrix W(r) in (1) reads as follows

W(r) =       1 α r 1 1 Y r α r 1 1 1 Y r 1 1 1 α r Y r 1 1 α r 1 Y r 1 1 1 1 β r       (5)
RR n°7726

where α r = 1 -2n 2 r ; β r = 2n 2 r -4. The parameters n r and Y r = 4n 2 r -4 denote the refraction index and the local admittance, respectively. In order to satisfy the boundary conditions, the scattering matrix W(r) must be modified as described, e.g., in [START_REF] Chopard | Lattice-Boltzmann method for wave propagation in urban microcells[END_REF].

The complex scalar electrical field strength Ψ(r, t) is given by

Ψ(r, t) = 1 n 2 r   Y r f0 + ∑ d=E,W,S,N -→ f d (r, t)   . ( 6
)

Multi-Resolution Frequency Domain (MR-FD) Approach

Taking the Fourier transform of both sides of the local scattering equation (1) leads to

-→ F (r, ν) = W(r)e -j2πν△t ← - F (r, ν) + -→ S (r, ν) (7)
where ν denotes the frequency. By aggregating the incident flow vectors ← -F (r, ν) for all nodes into a single incident flow vector ← -F (ν) and similarly gathering the reflected flow vectors -→ F (r, ν) and the source flow vectors -→ S (r, ν) into the global flow vectors -→ F (ν) and -→ S (ν), respectively, we come to the global scattering equation

-→ F (ν) = We -j2πν△t ← - F (ν) + -→ S (ν) (8)
where W is the global scattering matrix obtained by properly combining the local scattering matrices W(r) defined in [START_REF] Gorce | The Adaptive Multi-Resolution Frequency-Domain ParFlow (MR-FDPF) Method for Indoor Radio Wave Propagation Simulation. Part I : Theory and Algorithms[END_REF].

Taking into account the relationship between the incident and the reflected flows of the adjacent nodes expressed by (3), equation ( 8) can be rewritten in terms of the reflected flow vector 1

-→ F (ν) as -→ F (ν) = We -j2πν△t -→ F (ν) + -→ S (ν) (9)
where the matrix W signifies a permutated global scattering matrix W.

The solution to equation ( 9) can now be found as

-→ F (ν) = ( I -We -j2πν△t ) -1 -→ S (ν) (10)
where I denotes the identity matrix.

The multi-resolution approach [START_REF] Gorce | The Adaptive Multi-Resolution Frequency-Domain ParFlow (MR-FDPF) Method for Indoor Radio Wave Propagation Simulation. Part I : Theory and Algorithms[END_REF] provides a computationally efficient way of solving the equation ( 9). The approach is based on combining all the grid nodes into a single large head node that encompasses the considered propagation environment. After imposing prespecified boundary conditions, the head node is decomposed back into the original grid nodes. The process of applying the multi-resolution approach is divided into three steps: preprocessing step, upward step, and downward step. Further details related to the implementation of the MF-FD approach can be found in [START_REF] Gorce | The Adaptive Multi-Resolution Frequency-Domain ParFlow (MR-FDPF) Method for Indoor Radio Wave Propagation Simulation. Part I : Theory and Algorithms[END_REF][START_REF] Gorce | Deterministic approach for fast simulations of indoor radio wave propagation[END_REF]. Here, it is sufficient to mention that the highest computational load is associated with the preprocessing step. For a propagation environment of dimensions M × M nodes the computational loads for the preprocessing step, the upward step, and the downward step are estimated as [START_REF] Gorce | Deterministic approach for fast simulations of indoor radio wave propagation[END_REF]: 2 The computational cost is given per signal source.

• preprocessing step: O{52 • M 3 }; • upward step 2 : O{3 • M 2 }; • downward step: O{10 • log 2 (M ) • M 2 }.

RR n°7726

3 Optimization of MR-FDPF method for multifrequency simulations A major limitation of the MR-FD approach is that it provides a solution -→ F (ν)(or, equivalently, ← -F (ν)) of the global scattering equation ( 9) only at one frequency ν. That means that the MR-FDPF method has to be applied repeatedly at every frequency point in the range of interest B, e.g., for every sub-channel carrier frequency in OFDM communication systems.

In this section, a computationally efficient implementation of the MR-FDPF method is presented. The proposed implementation allows for a significant reduction in the computational load as compared to a straightforward repetition of the MR-FDPF method at multiple frequencies at the cost of acceptable degradation of the accuracy in predicting the electromagnetic wave propagation.

Without any loss of generality, we consider that only one signal source is present. Thus, we can also assume

-→ S (ν) = -→ S ̸ = 0), for all ν ∈ B.
If it is not the case, both sides of equation ( 10) can always be multiplied by a correction scalar factor.

A solution (10) of equation ( 9) at the frequency ν 0 ∈ B is given by

-→ F (ν 0 ) = ( I -W(ν 0 ) ) -1 -→ S (11)
where W(ν 0 ) = We -j2πν0△t . Similar expression for the frequency ν 0 + △ν reads as

-→ F (ν 0 + △ν) = ( I -W(ν 0 + △ν) ) -1 -→ S . ( 12
)
Note that W(ν 0 + △ν) = e -j2π△ν△t W(ν 0 ). Equation ( 12) can be rewritten as

-→ F (ν 0 + △ν) = ( I -W(ν 0 ) + W(ν 0 ) -e -j2π△ν△t W(ν 0 ) ) -1 -→ S (13)
After inserting [START_REF] Chopard | Lattice-Boltzmann method for wave propagation in urban microcells[END_REF] into ( 13) we obtain

-→ F (ν 0 + △ν) = ( I - ( I -W(ν 0 ) ) -1 W △ν ) -1 -→ F (ν 0 ) (14)
where the matrix W △ν = (e -j2π△ν△t -1) W(ν 0 ). Using the Neumann series expansion (see, e.g., [START_REF] Meyer | Matrix Analysis and Applied Linear Algebra[END_REF]) we can write

-→ F (ν 0 + △ν) = ( ∞ ∑ n=0 [ ( I -W(ν 0 ) ) -1 W △ν ] n ) -→ F (ν 0 ) ( 15
)
where it is implied that A 0 = I for any square matrix A.

According to [START_REF] Meyer | Matrix Analysis and Applied Linear Algebra[END_REF], an N -th order approximation of ( 12) is given by

-→ F (ν 0 + △ν) ≈ -→ F (ν 0 ) + -→ F (1) (△ν) + . . . + -→ F (N ) (△ν) (16)
where -→ F (1) (△ν) =

( I -W(ν 0 ) ) -1 W △ν -→ F (ν 0 ) (17) RR n°7726 and -→ F (n+1) (△ν) = ( I -W(ν 0 ) ) -1 W △ν -→ F (n) (△ν) (18)
It is worth noticing that (18) (as well as ( 17)) is similar to [START_REF] Chopard | Lattice-Boltzmann method for wave propagation in urban microcells[END_REF]. The only difference is that the source vector -→ S in ( 11) is substituted by the vector W △ν -→ F (n) (△ν) in (18). Thus, the MR-FDPF method can readily be applied in order to find the vector -→ F (n+1) (△ν) in ( 18).

An important question is the convergence speed of the approximation (16). A sufficient condition for the convergence of the Neumann series is that the spectral radius of the matrix

( I -W(ν 0 )
) -1 (e -j2π△ν△t -1) W(ν 0 ) is less than 1 [START_REF] Meyer | Matrix Analysis and Applied Linear Algebra[END_REF]. Based on the fact that in practice, the eigenvalues of the scattering matrix W(ν) must be real and smaller than 1, it can be induced that the smaller the spectral radius of the matrix W(ν) (or, in different words, the larger the attenuation due to the propagation environment), the faster is the convergence of the approximation (16). The advantage of using the approximation (16) lies in the fact that the scattering matrix W(ν) is calculated only for the frequency ν 0 . That means that the preprocessing step (see [START_REF] Gorce | The Adaptive Multi-Resolution Frequency-Domain ParFlow (MR-FDPF) Method for Indoor Radio Wave Propagation Simulation. Part I : Theory and Algorithms[END_REF][START_REF] Gorce | Deterministic approach for fast simulations of indoor radio wave propagation[END_REF]), which is the most computationally expensive one, is conducted only once. On the downside, however, equation (18) implies that while calculating the terms -→ F (n) (△ν), n = 1, 2, . . ., each node is a signal source. Thus, the computational cost of the upward step for each additional term -→ F (n) (△ν) becomes O{log 2 (M )M 2 } (up to a constant scaling factor). Of course, the additional computational loads associated with the upward step and the downward step are linearly proportional to the number of terms N in the approximation (16). Nevertheless, for large propagation environments (large M ), the computational cost of the preprocessing step is dominating, hence the approximate solution (16) promises significant computational savings.

In Section 4, we assess the performance of the proposed MR-FDPF method implementation based on the simulation results obtained for two propagation scenarios.

Experimental Results

In this section, we present the results of applying the proposed implementation of the MR-FDPF approach to predicting radio coverage in two realistic propagation scenarios. Unfortunately, we are not able to provide any experimental results for the MR-FDTLM method described in Section ?? as, presently, no TLM based prediction tool is at our disposal.

The first propagation scenario is inside the CITI building, INSA university campus in Lyon. The building dimensions are approximately 40 × 60 m. The electrical field strength Ψ(r, ν) inside the CITI building and its immediate surroundings has been computed with a spatial resolution of 2 cm. In Fig. 2, the power P (r, ν) of the electrical field predicted by the MR-FDPF method at the frequency ν = 2.345 GHz is depicted. The position of the omnidirectional transmitting antenna is indicated in Fig. 2.

We are to compare the electric field power P (r, ν) predicted with the MR-FDPF method optimized for multifrequency simulations as described in Section 3 to the electric field power P (r, ν) predicted with the traditional implementation of the MR-FDPF method, i.e., the implementation where the preprocessing step takes place at every frequency ν.

In the following, the elapsed time it takes to complete the preprocessing step is referred to as the preprocessing time. In a similar manner, the time period spent on the radio coverage prediction, i.e., the upward step together with the downward step, is called the prediction time. For the traditional MR-FDPF implementation, the total preprocessing time T prepr and the total prediction time T pred are K •620 s and K •40 s, respectively, where K is the number of frequencies. For the proposed implementation of the MR-FDPF method optimized for multifrequency simulations, the total preprocessing time is T prepr = 620 s, while the total prediction time is given by

T pred = K • 40 • k.
The factor k depends on the chosen number of terms N in the approximation (16). For example, k = 5.7, 10.2, and 15.66 for N = 1, 2, and 3, respectively. 16) is shown in Fig. 3. In Figs. 4 and5, examples of the cumulative distributions functions (CDFs) of the absolute approximation error are shown for several frequencies ν and several numbers of terms N , respectively. As expected, a larger deviation △ν from the central frequency ν 0 causes the approximation error to expand. Also, as the number of terms N in (16) increases, the approximation error decreases, or, alternatively, we observe expansion of the frequency band, in which the error is bounded by a certain level. For example, the mean approximation error is bounded by 3 dB in the frequency range of 3 MHz, 4 MHz for N = 1, N = 2 terms, respectively. Based on the electrical field strength Ψ(r, ν) predicted with the MR-FDPF implementations at the frequencies ν = 2.345, 2.346. . . . , 2.355 GHz, the delay spread has been estimated at 100 randomly selected positions. The chosen positions satisfy the condition that the approximation error is lower than 3 dB at the frequency ν = 2.345 GHz. In order to improve the estimate of the delay spread, the assessed power delay profiles (PDPs) have been averaged over 7 × 7 -point regions centered at the selected positions. The estimated delay spreads for the considered realizations of the MR-FDPF method are demonstrated in Fig. 6. In the figure, the red dashed line corresponds to the perfect match between the estimated delay spreads. Note that for larger delay spreads, the values obtained with the MR-FDPF method optimized for multifrequency simulations are slightly overestimated as compared to the delay spreads determined with the traditional implementation of the MR-FDPF method. The average relative absolute error is equal to 0.048.

RR n°7726

Figure 7 demonstrates the spatial distribution of the approximation error at the frequency ν = 2.347 GHz provided by the N = 2 term approximation. Note that the extreme values of the approximation error appear only at a limited number of points at the border of the propagation environment. It can also be seen from Fig. 7 that there exists a correlation between the predicted power P (r, ν) and the approximation error. Indeed, the smaller the power P (r, ν), the larger is the approximation error | P (r, ν)[dBm] -P (r, ν) [dBm]|. This observation is supported by the plot in Fig. 8, which illustrates the distribution of the observation points with respect to (w.r.t.) the predicted power P (r, ν) and the absolute approximation error.

The second propagation scenario is an indoor office environment in Stanford University. The dimensions of the environment are approximately 16×34 m. The positions of 8 transmitters equipped with omnidirectional antennas are depicted in Fig. 9. The electrical field strength Ψ(r, ν) has been calculated with the 2 cm spatial resolution.

The power P (r, ν) of the electrical field predicted by the MR-FDPF method at the frequency ν = 2.345 GHz is shown in Fig. 10. The delay spreads estimated with the MR-FDPF realizations for the considered propagation environment are depicted in Fig. 13. Again, the delay spreads have been RR n°7726 Fig. 14. Similar to the CITI propagation scenario considered above (although, less evident), the approximation error is lower at the points where the predicted power P (r, ν) is relatively larger (cf. the results presented in Fig. 7). The relationship between the RR n°7726

Conclusion

In this report, a new implementation of the MR-FDPF method optimized for multifrequency simulations has been presented. The described implementation can reduce the computational cost of the radio coverage prediction at multiple frequencies as compared to the traditional MR-FDPF method.

Experimental results provided in this paper demonstrate that the proposed implementation of the MR-FDPF method predicts the radio coverage in a range of frequencies with acceptable accuracy and without excessive computational overload. In general, the accuracy of the prediction or, alternatively, the frequency bandwidth, in which a specified accuracy is sustained, is dependent on the modeled propagation environment. It has been noticed that the prediction accuracy degrades as the distance from the transmitter increases.

Figure 1 :

 1 Figure 1: Incident and reflected flows.

Figure 2 :

 2 Figure 2: Radio coverage prediction inside the CITI building at ν = 2.345 GHz.

Figure 3 :

 3 Figure 3: CITI. Average approximation error (ν 0 = 2.35 GHz).

Figure 4 :

 4 Figure 4: CITI. CDF of the absolute approximation error N = 2 (ν 0 = 2.35 GHz).

Figure 5 :

 5 Figure 5: CITI. CDF of the absolute approximation error ν = 2.347 GHz (ν 0 = 2.35 GHz).

  Delay spread for proposed MR-FDPF (µs) Delay spread for traditional MR-FDPF (µs)

Figure 6 :

 6 Figure 6: CITI. Estimated delay spread N = 2 (ν 0 = 2.35 GHz).

Figure 7 :

 7 Figure 7: CITI. Spatial distribution of the approximation error at ν = 2.347 GHz and N = 2 (ν 0 = 2.35 GHz).

Figure 8 :

 8 Figure 8: CITI. Distribution of the observation points w.r.t. the predicted power and the absolute approximation error (ν = 2.347 GHz, N = 2, ν 0 = 2.35 GHz).

Figure 9 :Figure 10 :

 910 Figure 9: Stanford. Propagation environment.

Figure 11 :

 11 Figure 11: Stanford. Average approximation error (N = 2, ν 0 = 2.35 GHz).

Figure 12 :

 12 Figure 12: Stanford. CDF of the absolute approximation error (N = 2, ν = 2.345 GHz, ν 0 = 2.35 GHz).

  Delay spread for traditional MR-FDPF (µs) Delay spread for proposed MR-FDPF (µs)

Figure 13 :

 13 Figure 13: Stanford. Estimated delay spread N = 2, (ν 0 = 2.35 GHz).

Figure 14 :

 14 Figure 14: Stanford. Spatial distribution of the approximation error for TX 6 at ν = 2.347 GHz and N = 2 (ν 0 = 2.35 GHz).

Figure 15 :

 15 Figure 15: Stanford. Distribution of the observation points w.r.t. the predicted power and the absolute approximation error (ν = 2.347 GHz, N = 2, ν 0 = 2.35 GHz).

RR n°7726

Similarly, equation (8) can be reformulated in terms of the incident flow vector ← -F (ν).RR n°7726