Christopher Herzog

Maxim Olshanii

Yvan Castin

A liquid-gas transition for bosons with attractive interaction in one dimension

Keywords: unidimensional Bose gas, bright soliton, quantum liquid, ultracold atoms

A liquid-gas transition for bosons with attractive interaction in one dimension

Christopher Herzog, Maxim Olshanii, Yvan Castin

Introduction and motivations

In the last twenty years, ultracold atom experiments have been revealing their extreme flexibility in controlling the parameters and the geometry. By freezing with a laser trap the atomic motion in its vibrational ground state along one or two spatial dimensions, one can prepare systems of effective reduced dimensionality, which has led for example to the observation of the Berezinskii-Kosterlitz-Thouless transition [1] in a two-dimensional Bose gas. Furthermore, simply by applying a magnetic field close to a Feshbach resonance, one can adjust at will the s-wave scattering length, that is the amplitude of the atomic interactions [2]. The adiabatic use of both tools (transverse confinement followed by a switch of the scattering length from positive to negative values) to an almost pure three-dimensional Bose-Einstein condensate has allowed the experimentalists to observe the first N-mers in a one-dimensional atomic gas with attractive interactions [3,4], that is of N-body bound states whose existence was predicted by the purely quantum theory for all N > 1 [5], and that correspond for large N to the bright solitons of the classical field (Gross-Pitaevskii or non-linear Schrödinger) equation.

In the present work, we identify another mechanism of formation of those N-mers, which is purely at thermal equilibrium 1 and in a totally different physical regime. We indeed show that a one-dimensional spatially homogeneous attractive Bose gas in a non-degenerate regime can be thermodynamically less favoured (that is of a higher free energy) than an almost pure N-meric phase. In other words, the non-degenerate gaseous phase may exhibit a first-order phase transition to an almost fully bound phase, that we shall term a "liquid" in a thermodynamic limit that we shall specify.

It turns out that the concept of a liquid-gas transition in a system of ultracold bosonic atoms is quite timely in three dimensions, in the so-called unitary limit of an interaction with negligible range and infinite scattering length. After the numerical discovery of bound states that are probably of a liquid nature [6], from N = 3 (the Efimov trimers [7]) up to the largest accessible values of N (a few tens), such a transition was very recently observed in quantum Monte Carlo simulations with about one hundred particles [8]. The experimental realisation, however, is hampered by strong three-body losses induced by the Efimov effect [9,10], and may require a dedicated experiment.

In this context, our study of the one-dimensional case has both a theoretical and a practical interest. On one hand, the integrability of the theory allows one, under conditions to be specified in section 2, to construct in a controlled way a model that is simple to study numerically and to interpret analytically, as done in section 3, even in appropriate thermodynamic limits, see section 4. On the other hand, the smallness of the three-body losses, here away from the unitary limit, opens the way to a short term experimental realisation with ultracold atoms; we shall come back to this point in the conclusion.

Construction of the model

We consider N bosonic spinless quantum particles of mass m, living in one dimension on the x axis in the absence of an external potential, and with binary interactions via the attractive Dirac delta function potential V(x ix j ) = gδ(x ix j ), that is g < 0. A fundamental property of this hamiltonian model is its integrability, since one can determine the eigenstate wavefunctions with the Bethe ansatz [START_REF] Gaudin | La fonction d'onde de Bethe[END_REF], the price to pay with respect to the usual repulsive case [START_REF] Lieb | [END_REF] being that one must obviously include complex quasi-wavevectors [13,14] due to the existence of negative energy eigenstates, as it was already pointed out in [START_REF] Lieb | [END_REF] and explicitly done for N = 2 in an appendix of that reference. In free space-The Bethe ansatz solution is simple on the open line, that is in a pure scattering problem of the N bosons; in this case, the values of the quasi-wavevectors are known explicitly [13]. Under the condition that the Nbody wavefunction does not diverge at infinity, one finds that each eigenstate is composed of an arbitrary collection of indistinguishable n-mers, the special case n = 1 corresponding to an unbound bosonic particle, that is an atom. Each n-mer is characterized by a well-defined total momentum K, and elastically scatters off the other n-mers, since any chemical reaction (association or dissociation), or even any simple backscattering of two n-mers of different sizes, is forbidden by integrability. For all n ≥ 2, there exists one and only one possible bound state, of internal energy [5,13]

E 0 (n) = - mg 2 24 2 n(n 2 -1) (1)
In a given scattering state, there can be of course several n-mers with the same size, so that the total number of fragments may vary from one (all the particles are bound in the ground N-mer, as in a liquid) up to N (all the particles exist under an atomic, gaseous form). In terms of the size n i of the ith fragment and of its center of mass wave vector K i , the scattering state energy may be simply written as the following sum of internal and kinetic energies:

E = i       E 0 (n i ) + 2 K 2 i 2n i m       (2)
In the box -We shall here study the system at thermal equilibrium at temperature T , so we have to enclose it in a quantisation box of length L, with the usual periodic boundary conditions 2 . The quasi-wavevectors of the Bethe ansatz then solve a non-linear system of coupled equations [14] with no analytical solution, which makes it difficult to determine the discrete energy spectrum in the box, even numerically, at large N. The central idea of the present work is thus to restrict to a limiting case that may be handled in a simple way, such that the free space expression (2) remains approximately true in the box, with (as the only twist) the natural quantisation prescription:

K i ∈ 2π L Z (3)
A box that is larger than the n-mers -A first condition to get the expression (2) for the spectrum is that the structure and the internal energy of the bound state are weakly affected by the presence of the box. As the dimer state is the less strongly bound state, with the largest spatial diameter, this requires:

L ℓ = π |µ 0 | E F 1/2 ≫ 1 (4)
2. Even if it is possible to produce box potentials [15] it is useful to transpose our work to the harmonically trapped case, as was done in the pioneering reference [16].

where ℓ ≡ 2 2 /(m|g|) is the diameter of the dimer, µ 0 ≡ -mg 2 N 2 /( 82 ) is the chemical potential of the N-mer, and E F = 2 (πρ) 2 /(2m) is the Fermi energy of the fictitious one-dimensional gas of fermions with the same density ρ = N/L as the bosonic particles. As the n-mer wavefunction is an exponentially decreasing function of the sum of the distances between the particles [5,13], we expect that neglecting the effect of the box on the internal energy E 0 (n) introduces an exponentially small error in L. For N = 2, this expectation is confirmed by the appendix A of [START_REF] Lieb | [END_REF]: Two "bound states" are found, corresponding to our dimer state with total momenta K = 0 and K = 2π /L, and their energies deviate from (2) by ≃ ±4 exp(-L/ℓ) in relative value. In the non-degenerate regime -The condition (4) guides all the present work. Since |µ 0 | is the binding energy of a particle in the N-mer, one can reasonably expect that a liquid-gas transition (if there is any) takes place at a temperature of the order of a fraction of

|µ 0 |/k B ≫ E F /k B ,
that is in a regime where the atomic gaseous phase is highly non degenerate:

ρλ = 4E F πk B T 1/2 ≪ 1 (5)
where λ = [2π 2 /(mk B T )] 1/2 is the usual thermal de Broglie wavelength. This atomic phase is then also close to the continuous spectrum limit, meaning that any discrete sum over the K i in equation ( 3) can be replaced by an integral, since L ≫ Nλ implies L ≫ λ. In short, this phase can be considered as a classical gas. In what follows we consider this high temperature limit, where the calculation of the external partition function Z ext n-mères (that is of kinetic origin) of the n-mers associated to the spectrum (2) is greatly simplified. For a given internal configuration of the N particles, the subset of the n-mers forms a gas (supposed to be ideal for the moment) with N n indistinguishable (bosonic) elements of mass nm ≥ m, de Broglie wavelength λ/n 1/2 ≤ λ and density ≤ ρ, so that this gas is as close to the classical limit as the (n = 1) atomic gas, which leads to [START_REF] Diu | Physique Statistique[END_REF]

Z ext n-mères ≃ Ln 1/2 λ N n 1 N n ! (6)
Negligible interactions among the n-mers -A second validity condition for the use of the energy spectrum defined by Eqs. (2,3) is that the energy shifts in the box due to the elastic scattering, i.e. interaction among n-mers, are negligible as compared to the mean kinetic energy of each internal configuration 3 . In the purely atomic phase, the estimation of the interaction energy is made simpler by the following remark: at the considered temperatures of order |µ 0 |/k B , the kinetic energy per atom ≈ k B T/2 is much larger than the dimer binding energy |E 0 (2)| as soon as N ≫ 1, which leads to

λ ℓ = 1 N 4π|µ 0 | k B T 1/2 ≪ 1 (7)
so that the typical relative wavevector of two atoms ≈ 1/λ is much larger than m|g|/ 2 , and two-body scattering takes place in the Born regime, see equation (99) in the reference [13] 4 . The interaction energy per atom is then simply the one ρg of the mean field theory, including the factor two due to the bosonic bunching effect, and it can be neglected if

ρ|g| k B T = 4 π (E F |µ 0 |) 1/2 Nk B T = 1 π (ρλ) λ ℓ ≪ 1 (8)
This conditions seems to be generically satisfied, similarly to (7), in some sort of thermodynamic limit, since N appears in the denominator; this will be confirmed by section 4. Furthermore, this condition is implied by the product of conditions (5,7), so that it can be omitted in what follows.

To get an estimate of the interaction energy of an atom with a quasi-classical n-mer (n ≫ 1), let us treat the corresponding bright soliton with Bogoliubov theory, where the atom-n-mer scattering process at relative wave vector k is described by the quasi-particle mode functions u k (x) and v k (x). If u k (x) ∼ e ikx for k > 0 and x → -∞, the exact expression of the modal functions [START_REF] Kaup | [END_REF]19,20] leads to u k (x) ∼ e iθ(k) e ikx for x → +∞, with a phase shift θ(k) = -4 arctan(kℓ/n). The periodic boundary conditions then impose the quantization condition exp(ikL) exp[iθ(k)] = 1. Hence in our big box L ≫ ℓ/n, a spacing between successive wave numbers equals δk ≃ 2π/[L + θ ′ (k)] where θ ′ (k) is the derivative of θ(k). The effective quantization length at a given thermal kinetic energy becomes

L eff ≈ L + θ ′ (k = √ 2π/λ) (9)
We insert this effective length in the partition function (6) written for a single atom, and we approximate θ ′ (k) by its asymptotic expression ∼ -4n/(k 2 ℓ) at large k, which is justified by the fact that k B T ≈ |µ 0 |. We then find a length change that is very small in relative value, and that corresponds, for an atom-n-mer pair, to an interaction energy of order 2n|g|/L, that is of order n times the interaction energy of a pair of atoms 5 . Thus the same condition (8) obtained for interactions among atoms allows us to neglect interactions between atoms and n-mers.

The model -To conclude, for the validity conditions (4,5,7) given above, our high-temperature system of N bosonic particles with attractive interaction in a quantization box, may be assumed to exhibit, as internal configurations, all the possible partitions of the N particles into atoms and into two-body, . . . , N-body bound states, each subset of atoms or bound-states of a common size n forming an ideal classical gas of N n indistinguishable particles of mass nm, which corresponds to the partition function (6). In mathematical terms, the full partition function Z of our model is the sum over all internal configurations with a fixed total number N of particles, of configurational partition functions:

Z = conf Z conf , where conf ∈ {(N n ) 1≤n≤N ∈ N N such that N n=1 nN n = N}, ( 10
)
that can be written, with β = 1/(k B T ), as products of internal and external partition functions,

Z conf = exp        -β N n=1 N n E 0 (n)        N n=1        Ln 1/2 λ N n 1 N n !        (11)

Mesoscopic numerical investigations and interpretation in terms of a liquid-gas transition

We have performed a direct numerical study of the model (10,[START_REF] Gaudin | La fonction d'onde de Bethe[END_REF] simply by generating all possible internal configurations (N n ) 1≤n≤N . The numerical results are thus as exact as the model is. The number of configurations to explore is, however, equal to the partition number p(N) of the integer N, which increases exponentially with N [21]. This limits the simulations on a desktop computer to N 150. For a total number of N = 100 particles and a fixed geometry L/ℓ = 100, which lies well within the large box limit (4), we show in figure 1a and 1b the dependence on temperature of two observables, the total free energy F = -k B T ln Z and the mean fragmentation rate ν , given by a thermal average of the relative number of fragments in each configuration with weights Z conf ,

ν = 1 N N n=1 N n ∈ 1 N , 1 (12)
One observes in figure 1a a rather abrupt change of slope of the free energy, around T ≃ 105T F ≃ 0.1|µ 0 |/k B , which is well within the non-degenerate regime (5). Since the derivative -∂ T F gives the entropy of the system, this change of slope may be the precursor for our mesoscopic system of a first order phase transition with a non-zero latent heat. Around the same temperature, there is a rather abrupt change of the mean fragmentation rate, from a value close to unity at high temperature, to a value close to 1/N, see figure 1b: This change suggests the transition from a mainly atomic gaseous phase (i.e. an almost fully dissociated phase) to a liquid phase mainly formed of the ground N-mer. This expectation is confirmed by the shape of the distribution (n N n /N) 1≤n≤N of the number of particles per size of the n-mers, which exhibits at the transition two maxima, one at n = 1 and the other at n = N, at the top of narrow peaks, separated by infinitesimal occupation numbers for intermediate sizes, see figure 1c. If one increases (or reduces) the temperature, one observes that the n = N peak (the n = 1 peak respectively) rapidly collapses (not shown here). Let us try to provide a more quantitative support to our liquid-gas transition scenario. At a first approximation, let us keep in the partition function Z only two internal configurations, the one with N atoms, contributing as Z at , and the one of the pure N-mer, contributing as Z N-mère . The corresponding estimates for the free energy and the fragmentation rate,

F (0) = -k B T ln(Z at + Z N-mère ) and ν (0) = (N -1 Z N-mère + Z at )/(Z at + Z N-mère )
, are plotted as dotted lines in figure 1a and1b. This reproduces well the general behavior of the numerical results, up to about a 3.5% shift towards higher temperatures, as we have checked explicitly. The corresponding transition temperature T (0) c , obtained by equating the partition functions Z at and Z N-mère , which amounts here to taking a mean fragmentation rate equal to (1

+ 1/N)/2, solves e -β (0) c E 0 (N) L √ N λ (0) c = 1 N! L λ (0) c N ( 13
)
For the parameters of figure 1, this leads to T (0) c ≃ 108.90T F , whereas the mean fragmentation rate in the numerical simulation reaches the same value (1 + 1/N)/2 at T ≃ 105.1T F .

Note that a convenient approximate form of equation ( 13) is obtained by taking the logarithm and then dividing by N. One thus takes the large N limit, without precisely specifying yet how the thermodynamic limit must be taken for L (this will be specified in section 4). One only makes the weakly restrictive hypothesis that L/λ (0) c increases less rapidly than some power of N, so that ln(L/λ (0) c ) = O(ln N). Using Stirling asymptotic formula for N!, one gets

1 3 β (0) c |µ 0 | -ln       e ρλ (0) c       = O ln N N ( 14
)
Neglecting the terms that vanish at large N, and taking as unknown the dimensionless quantity

β (0) c |µ 0 |, one reaches the implicit equation 2 3 β (0) c |µ 0 | + ln 2 3 β (0) c |µ 0 | = ln πe 2 |µ 0 | 6E F ( 15
)
that can be solved in terms of the Lambert W function (inverse function of x → xe x for x ≥ -1)

as 6 2 3 β (0) c |µ 0 | = W πe 2 |µ 0 | 6E F ( 16
)
6. To the same level of approximation, one predicts a latent heat per particle given by According to the validity condition (4) of the model, the argument of W must be much larger than unity in equation (16). For the parameters of figure 1, this large N approximation to T (0) c reproduces the exact value within only 3% from below, even if the number of particles N = 100 is not very large.

Q (0) /N = 1 3 |µ 0 | + 1 2 k B T (0) c ,
Let us now go beyond the two pure phases model and try to understand why its prediction (dotted line in figure 1) is slightly shifted along the temperature axis, with respect to the numerical results. At this stage, we can point out another shortcoming of this model: at temperatures beyond the cross-over region, it assumes in practice that the system is purely atomic with a unit fragmentation rate, whereas the numerical results present a persistent deviation of about 10% from full fragmentation. This last observation is consistent with a non-negligible population of dimers at the transition, see the channel n = 2 in figure 1c, and gives the idea of a second order approximation, that includes, in addition to the pure atomic phase, arbitrary mixtures of atoms and dimers, so as to obtain an atomo-dimeric phase 7 .

To be symmetric, we also improve the competing liquid phase, even if this is quantitatively not really required for the figure 1a and 1b (the dotted line is much closer to the solid line for T ≤ 105T F than for T > 105T F ): to the N-mer we add the minimally dissociated configuration with one (N -1)-mer and one atom, whose partition function is ≈ e -β|µ 0 | L λ Z N-mère for N ≫ 1. Note here the occurrence of the expected thermal activation law exponential factor, which reveals the fact that the N-mer state is protected from dissociation by an energy gap |µ 0 | (in the center of mass frame) below the continuum of N-body excited states. On the contrary, the minimally dissociated phase is favored by the kinetic entropic factor ∝ L/λ provided by the relative motion of the (N -1)-mer and the atom. In the particular case of figure 1, the effect of the energy gap largely dominates over the kinetic entropic effect, e -β|µ 0 | L λ ≈ 0.06 at the transition, but this is not necessarily true for a large system where L/λ can diverge.

For the extended phases that we have just described (the atomo-dimeric phase and the up-to-an-atom pure N-meric phase), the free energy and the mean fragmentation rate are plotted as dashed lines on the figure 1a and1b. The dashed lines are significantly closer to the numerical results than the dotted line: including the dimers has lowered the free energy of the gaseous phase and explains most of the shift in temperature. The physics for this N = 100 particle system thus seems to be quantitatively well understood.

We cannot be sure however that the phase extension performed here is sufficient in the limit of a large system. In the gaseous phase, it remains to clarify analytically the role of dimers, or even trimers, quadrimers, etc; in the liquid phase, as we have seen, inclusion of more severe dissociations of the N-mer may become necessary. One can even imagine that a phase of totally different nature may emerge at large N. Can all this lead to large deviations of the critical temperature from the simple estimate (13,16), or even invalidate it? These questions that cannot easily be answered by our mesoscopic numerical calculations are treated in the next section.

Existence of the phase transition and various limits for a large number of particles

According to the numerical simulations of the previous section, limited to a hundred of particles, our onedimensional system of bosons with attractive interaction should exhibit, in some macroscopic limit, a first order phase transition between a gaseous phase and a liquid phase, with a critical temperature estimated within a two pure phases approximation, see equation (13). The goal of this section is first to give an analytical proof of the existence of such a transition, with an upper bound on the error on T c in (13); then to construct various types of thermodynamic limits, to see how they are influenced by finite size effects, and to determine to what extent the two competing phases are mainly atomic or mainly liquid, respectively.

Existence of the transition and bound on the error on T c

Appropriate bounds on the partition function Z conf of equation [START_REF] Gaudin | La fonction d'onde de Bethe[END_REF] will allow one to show in an elementary way that deviations from monoatomicity in the gaseous phase and impurity of the liquid phase do not destroy the scenario of a liquid-gas transition, provided that this transition takes place in the non-degenerate regime ρλ ≪ 1. This excludes the possibility that a third phase, that would have been missed by the mesoscopic simulations, may emerge and dominate in the large N limit. To this end, we shall define, in the configurational space, neighbourhoods of the

(z) = 2 F 0 (-N 2 , 1 2 -N 2 ; 4z).
purely atomic and N-meric phases, and we shall bound the corrections that they bring to the free energy per particle, and hence to T c , with respect to the two pure phases approximation of the previous section.

In the neighbourhood of the monoatomic phase -Let us first take as a reference the purely atomic phase, with a partition function Z at = (L/λ) N /N!, by forming the ratio Z conf /Z at and eliminating the atom number N 1 thanks to the conservation of the total number of particles N, N 1 = N -N n=2 nN n . This elimination is made simpler by the expected property that E 0 (1) = 0 in equation (1). The only non trivial remaining factor originates from the factorial due to atomic indistinguishability, that we shall simply bound,

N! N 1 ! =         N k=1+N 1 k         ≤ N N-N 1 = N N n=2 nN n ( 17
)
to obtain

Z conf Z at ≤ N n=2 [α n (N)] N n N n ! with α n (N) ≡ n 1/2 e -βE 0 (n) N Nλ L n-1 (18)
Intuitively, this inequality is almost saturated, and the integer numbers N n for n ≥ 2 almost obey a Poisson distribution with parameter α n (N) if a) the number of atoms in the considered phase only very weakly differs from N in relative value (in particular one should have |1 -N 1 /N| ≪ 1), and b) the truncation effect due to conservation of the total number of particles is negligible on the distribution of N n (in particular one should have α n (N) ≪ N/n).

A simple physical interpretation of the quantity α n (N) is obtained as follows: The partition function of a gas of Nn atoms in the presence of a single n-mer, with 2 ≤ n ≪ N, is approximately equal to α n (N)Z at , so that α n is the first correction (in relative value) to the atomic gas partition function due to the n-mers. In the expression for α n (N), the factor N n ≫ 1 corresponds to an increase in the counting entropy, induced by the substitution of n atoms by a distinguishable n-mer. For a fixed N, the larger the size n of the n-mer is, the smaller α n (N) and the less probable it is to include the n-mer, since ρλ ≪ 1; this already suggests that the trimers have a weaker effect on the free energy of the gaseous phase than the dimers, and so on.

Let us now define the neighbourhood of maximal clustering order n max for the atomic phase, by the set of all possible configurations composed of n-mers of arbitrary sizes n up to n max . We shall take n max less than N/2 so as to avoid overlap with the neighbourhood of the N-mer to be defined soon. Our neighbourhood of maximal order n max contains the purely atomic phase as a particular case, so that its partition function Z (n≤n max ) vois. at. is bounded from below by Z at . It can be bounded from above by applying [START_REF] Kaup | [END_REF] and extending the summation over the (N n ) 2≤n≤n max to N n max -1 :

Z at ≤ Z (n≤n max ) vois. at. ≡ conf/N n =0 ∀n>n max Z conf ≤ Z at e nmax n=2 α n (N) (19)
We shall use these inequalities to bound the change of critical temperature due to the presence of dimers, trimers, etc, in the gaseous phase.

In the neighbourhood of the N-mer -Let us now take as a reference the pure N-mer, with a partition function

Z N-mère = exp[-βE 0 (N)]N 1/2 L/λ. For k < N -k, that is k < N/2,
let us define its subneighbourhood of dissociation degree k as the set of configurations with one (N -k)-mer in presence of an arbitrary distribution of the remaining k particles among atoms, dimers, . . . or at most one k-mer. Its partition function will be noted as Z (k) ss-vois. N-mère . The particular case k = 0 corresponds to the pure N-mer.

Let us then apply the previous upper bound reasoning, up to slight modifications: as one necessarily has N N-k = 1, and elsewhere N n = 0 as soon as n > k, the number of atoms is now more conveniently eliminated by the relation

N 1 = k - k n=2 nN n so that k! N 1 ! =         k s=1+N 1 s         ≤ k k-N 1 = k k n=2 nN n ( 20
)
One also has (Nk) 1/2 ≤ N 1/2 . All this leads to an upper bound for Z (k) ss-vois. N-mère in the same spirit as [START_REF] Kaup | [END_REF], except that the first corrections α n (due to the n-mers) to the atomic phase surrounding the (Nk)-mer now apply to k particles only, rather than N, hence the substitution in equation ( 18) of α n (N) by α n (k), of Z at by (L/λ) k /k! and of the upper bound N by k:

Z conf Z N-mère ≤ e -β[E 0 (N-k)-E 0 (N)] (L/λ) k k! k n=2 [α n (k)] N n N n ! (21)
Note that the first (exponential) factor is here a barrier to dissociation, according to the expected activation law, corresponding to the minimal energy E 0 (Nk) -E 0 (N) > 0 required to extract k particles from the N-mer, an energy which is close to k|µ 0 | for k ≪ N. On the contrary, the powers of L/λ ≫ 1 favor dissociation due to a gain in kinetic entropy. Another point is that, as in the discussion that follows equation ( 18), the inequality ( 21) is expected to be saturated when the product of dissociation is almost purely atomic, in which case the n-mers (with n ≥ 2) are distributed according to Poisson laws of parameters α n (k). Let us now overbound the summation over all configurations (N n ) 2≤n≤k of the subneighbourhood of degree k by extending this summation to N k-1 , and let us define the N-mer neighbourhood of maximal degree of dissociation k max , by collecting all the subneighbourhoods of degree k ≤ k max , which results in the partition function

Z (k≤k max ) vois. N-mère ≡ k max k=0 Z (k)
ss-vois. N-mère [START_REF] Dunjko | [END_REF] We of course limit this definition to k max < N/2. We finally obtain the inequalities

Z N-mère ≤ Z (k≤k max ) vois. N-mère ≡ conf/∃n≥N-k max /N n 0 Z conf ≤ Z N-mère k max k=0 e -β[E 0 (N-k)-E 0 (N)] (L/λ) k k! e k n=2 α n (k) , (23)
that we shall use to bound the critical temperature change due to partial dissociation of the N-mer.

Implications for the phase transition -Let us test the two pure phases scenario of section 3, by extending the first phase to the gaseous phase containing n-mers of arbitrary size up to n max , and the second phase to the liquid phase containing the partially dissociated N-mer in the presence of an arbitrary phase with at most k max particles. The critical temperature T (e) c , where the exponent means "extended", is then obtained by requiring that the partition function Z (n≤n max ) vois. at. of the extended gaseous phase, see equation (19), be equal to the partition function Z (k≤k max ) vois. N-mère of the generalized liquid phase, see equation [START_REF] Landau | Physique statistique[END_REF]. The implicit equation ( 14) is replaced by

1 3 β (e) c |µ 0 | -ln       e ρλ (e) c       + O ln N N = 1 N        ln Z (n≤n max ) vois. at. Z at -ln Z (k≤k max ) vois. N-mère Z N-mère        (24)
In order to validate the two pure phases scenario, one has to show that the right-hand side of this equation is small as compared to one of the first two terms of the left-hand side. To this end, let us make the upper bounds more explicit in the sums over n of the α n (N) and the α n (k), using the fact that -βE 0 (n)/(n -1) = n(n + 1)β|µ 0 |/(3N 2 ), and that n(n + 1) is there at most n max (n max + 1) or k max (k max + 1):

n max n=2 α n (N) ≤ N n max n=2 n 1/2 X n-1 ≤ N √ 2X (1 -X) 2 with X = Nλ L e β|µ 0 |n max (n max +1)/(3N 2 ) ( 25
) k n=2 α n (k) ≤ k k max n=2 n 1/2 Y n-1 ≤ k √ 2Y (1 -Y) 2 with Y = k max λ L e β|µ 0 |k max (k max +1)/(3N 2 ) ( 26
)
For small enough values of n max and k max , we expect that X < 1 and Y < 1, a property that we used to overbound by extending the summations up to infinity 8 . Insertion of (25) in the inequalities (19) directly gives:

0 ≤ 1 N ln Z (n≤n max ) vois. at. Z at ≤ √ 2X (1 -X) 2 (27) 8. We used +∞ n=2 (n/2) 1/2 u n-1 ≤ +∞ n=2 (n/2)u n-1 = u(1 -u/2)/(1 -u) 2 ≤ u/(1 -u) 2 , ∀u ∈ [0, 1[.
To obtain a usable form of the upper bound [START_REF] Landau | Physique statistique[END_REF], we note that [E 0 (Nk) -E 0 (N)]/(k|µ 0 |) is a decreasing function of k and is smaller than its value η = 1 -(k max /N) + (k 2 max -1)/(3N 2 ) in k = k max . Then we use the upper bound in (26). Then we bound the resulting sum over k by the corresponding series over N:

0 ≤ 1 N ln Z (k≤k max ) vois. N-mère Z N-mère ≤ e -β|µ 0 |η ρλ e √ 2Y/(1-Y) 2 (28)
To confirm the existence of a liquid-gas phase transition at a temperature close to T (0) c defined in equation ( 13), in the case N ≫ 1, we maximally extend the two neighbourhoods, taking (for N odd) n max = k max = (N -1)/2. These maximal neighbourhoods are disjoint sets and they include as a whole all the possible internal configurations 9 . Let us define T c as the temperature of equality of the partition functions of these maximal neighbourhoods, and let us bound the supposedly small deviation of T c from T (0) c by linearizing the left-hand side of equation ( 24), and by evaluating the right-hand side of that equation for T = T (0) c , which allows us to eliminate β|µ 0 | in terms of ρλ in the quantities X and Y to obtain X = 2Y = e 1/4 (ρλ (0) c ) 3/4 . We find that the upper bound in equation ( 27) is larger than the one in (28). As √ 2e 1/4 < 2, we conclude that there exists a neighbourhood of ρλ (0) c = 0, that is of L/ℓ = +∞ according to (14,16), such that

|T c -T (0) c | T (0) c < 2(ρλ (0) c ) 3/4 ln[e/(ρλ (0) c )] → ρλ (0) c →0 0 ( 29
)
Therefore there is indeed a liquid-gas transition in such a neighbourhood, in thermodynamic limits that we now have to specify. We note en passant that the upper bound in (28), even if it is good enough in the present discussion on the existence of the transition close to T (0) c , has the bad property of diverging as 1/(ρλ) at very high temperature, so that it must be improved if one wishes to exclude the presence of a transition at temperatures arbitrarily higher than T (0) c . This is done in Appendix A.

A limit in which the model is exact

As is apparent in the caption of figure 1, the figure only depends on two dimensionless parameters, the number of particles N and the box length scaled by the diameter of the dimer, L/ℓ. To take the thermodynamic limit, the most natural way is to let N tend as usual to infinity for fixed interaction strength and mean density, which corresponds to

N → +∞ with ρℓ = constant (30)
This makes the first validity condition (4) of the model exactly satisfied. Let us first determine the critical temperature T (0) c of the simplest, two pure phases approximation of section 3. According to (4), the limit (30) implies that (|µ 0 |/E F ) 1/2 = N/(πρℓ) → +∞; this may look strange to the reader and we shall come back to this point in section 4.3. By iteration of the implicit equation (15), taking into account the fact that β (0) c |µ 0 | is much larger than its logarithm, we obtain the asymptotic behavior:

1 3 β (0) c |µ 0 | = N→+∞ ln N - 1 2 ln ln N + ln e ρℓ √ 12π
+ O ln ln N ln N (31)
When N diverges, the critical temperature T (0) c thus very slowly becomes arbitrarily smaller than the chemical potential of the N-mer. Furthermore, the last two validity conditions (5, 7) of the model become arbitrarily well satisfied at that temperature:

ρλ (0) c ρℓ = λ (0) c ℓ ∼ N→+∞ (12π ln N) 1/2 N → 0 (32)
In other words, the simple model (10,[START_REF] Gaudin | La fonction d'onde de Bethe[END_REF] is asymptotically exact. Last, due to the fact that ρλ (0) c → 0 essentially as a power law, the equation (29) gives the important result that the two pure phases approximation is also asymptotically exact for the calculation of T c :

(

k B T c -k B T (0) c )/|µ 0 | = O(1/N 3/4 ) → N→+∞ 0 (33) 9.
A configuration either has no n-mer with a size larger than N/2, in which case it belongs to the maximal neighbourhood of the monoatomic gas; or it has such a n-mer, and necessarily at most one, in which case the configuration belongs to the maximal neighbourhood of the N-mer.

This does not imply that the dimers have a negligible contribution to the gaseous phase partition function, or that their mean number in that phase tends to zero. From a study of the atomo-dimeric partition function, see section 3, or more simply by using the fact that the gaseous phase is here in the quasi-poissonian regime defined below equation [START_REF] Kaup | [END_REF], we find on the contrary that Z at-dim /Z at diverges as exp[α 2 (N)], and that the mean number of dimers diverges as α 2 (N), where α 2 (N) is defined in [START_REF] Kaup | [END_REF] and diverges as (ln N) 1/2 at the transition. This is however o(N) and too slow to asymptotically contribute to the free energy per particle or to the mean fragmentation rate. Even better, the inclusion of trimers brings, according to [START_REF] Kaup | [END_REF], a negligible contribution to the partition function as compared to Z at , since this contribution is bounded by exp[α 2 (N)]{exp[α 3 (N)] -1} → 0; also, the mean number of trimers tends to zero as α 3 (N) at the transition. Similarly, the liquid phase is asymptotically pure there, because the probability that it contains an unbound atom tends to zero as exp(-β c |µ 0 |)L/λ c , that is as ln N/N, see section 3. As a consequence, one can determine the scaling law of the thermodynamic limit simply using the two pure phases approximation. For example, one gets the following universal law for the mean fragmentation rate 10 :

ν → N→+∞ exp δ 1 + exp δ à δ ≡ T -T c T c N 1 2 + 1 3 β (0) c |µ 0 | fixé, (34)
This means that the relative temperature width of the cross-over region vanishes as 1/(N ln N).

A more usual thermodynamic limit

The large N limit of section 4.2 is inappropriate for an experimental realisation at fixed interaction strength, since |µ 0 | then diverges as N 2 and the condition for the one-dimensionality of the system (|µ 0 | small as compared to the transverse vibrational energy quantum ω ⊥ ) is asymptotically unreachable at the transition 11 . This gives the idea of fixing the value of Ng to a constant in the experiment, by tuning g with a Feshbach resonance. The limit (30) then corresponds to N 2 /L =constant, which is unusual. Let us rather consider here the limit

N → +∞ with L ℓ = constant (35)
that indeed leads to a constant density ρ at fixed Ng, and to a constant ratio |µ 0 |/E F as expected in a regular thermodynamic limit. Then β (0) c |µ 0 | has a finite limit, given by ( 16), and ρλ (0) c /e = exp(-β (0) c |µ 0 |/3). However, T (0) c no longer coincides with T c ; it is only a good approximation in a sufficiently non degenerate regime, as guaranteed by the mathematical result (29) 12 . The same conclusion holds for the universal relation (34), that now predicts a width of the cross-over region scaling as 1/N. The first two validity conditions (4,5) of the model are now only approximately obeyed; only the third condition ( 7) is exactly fulfilled when N → +∞ at the transition.

Even if the mathematical result (29) has to consider the maximal neighbourhoods of the two pure phases, with n max ∼ k max ∼ N/2, to be rigorous, it is physically expected that the first corrections to T (0) c actually originate from narrower neighbourhoods n max ≪ N/2 et k max ≪ N/2. As the close neighbourhood of the purely atomic gas is quasipoissonian at the transition, as defined below equation [START_REF] Kaup | [END_REF], we find that its leading contribution to the right-hand side of ( 24) is the one α 2 (N)/N → √ 2ρλ of the dimers 13 , the contribution of the trimers α 3 (N)/N being O(ρλ) 2 , etc. We realize however that there is a limitation to our model: although it is non degenerate, the purely atomic phase remains 10. One simply sets

T = T (0) c /(1 -ǫ), so that Z at /Z N-mère = exp[-ln(1 -ǫ)(N -1)/2] exp[-ǫβ (0) c E 0 (N)], with (N -1) ln(1 -ǫ) ≃ -Nǫ and β (0) c E 0 (N) ≃ -β (0) c |µ 0 |N/3.
Using (31) and a numerical opto-dimeric calculation (cf. the dashed line in figure 1), one finds for ρℓ = 1 that (34) is already almost reached for N = 1000.

11. This one-dimensionality condition also justifies the use of a Dirac δ function for the interaction potential, at least in the usual cold atom case where the real three-dimensional interaction V 3D has a negligible range as compared to the mean interparticle distance [3]. It was indeed shown in [START_REF] Dunjko | [END_REF] that, due to the transverse harmonic confinement, V 3D induces in one-dimension an effective interaction of non-zero effective range r e ≈ [ /(mω ⊥ )] 1/2 . This effective range is negligible in the Born regime under the sufficient condition k rel |r e | ≪ 1 where the relative wavevector is here k rel = O(N/ℓ).

12. The more precise forms (27,28) guarantee the existence of a solution T c to the equation Z (k≤N/2) vois. N-mère (T ) = Z (n≤N/2) vois. at. (T ) as soon as L/ℓ > 105.1. 13. To go beyond the poissonian approximation, we have applied the method of Laplace to the sum S (z) of the footnote 7, after having used Stirling's formula and replaced the sum over s by an integral over x = 2s/N. Then Z at-dim /Z at ∼ exp[Nu(x 0 )]/(1 + x 0 ) 1/2 when N → +∞, with u(x) = (x/2) ln[2 √ 2ρλ/(ex)] -(1x) ln(1x) and x 0 is the root of u ′ (x) in the interval [0, 1]. One always has u(x 0 ) < √ 2ρλ, in agreement with the upper bound (19) applied to n max = 2, and u(x 0 ) ∼ α 2 (N)/N if ρλ → 0. in reality a bosonic gas, and its partition function Z Bose at already deviates from the classical gas one Z at = (L/λ) N /N! to first order in ρλ [START_REF] Diu | Physique Statistique[END_REF], N -1 ln(Z Bose at /Z at ) = 2 -3/2 ρλ + O(ρλ) 2 , which brings a correction to the critical temperature to the same order in ρλ as the one of the dimers, but which is ignored by the model.

To estimate the contribution of the close neighbourhood of the N-mer, we take inspiration from the above discussion and assume, to leading order in ρλ, that the N-mer is weakly dissociated in a purely atomic gas with k particles, well in the poissonian regime. This amounts to saying that the second inequality in [START_REF] Landau | Physique statistique[END_REF] is almost saturated, that the α n (k) are negligible, and that the N-mer dissociates into a (Nk)-mer plus k particles with an essentially poissonian weight, since one takes k max ≪ N/2:

Z (k) ss-vois. N-mère /Z N-mère ≃ γ k k! where γ = N e -β|µ 0 | ρλ ( 36
)
At the transition, the mean dissociation rate of the liquid phase is then γ/N ≃ (ρλ (0) c ) 2 /e 3 , which brings to the righthand side of (24) a contribution negligible as compared to the one of the dimers of the gaseous phase. Beyond the model -In order to determine the true correction to T (0) c to first order in ρλ, let us apply to the gaseous phase 14 a systematic approach at high temperature, that is the quantum virial expansion, N -1 ln(Z Bose gaz /Z at ) = b 2 ρλ + O(ρλ) 2 , that includes at this order quantum statistics, atomic interactions and the presence of dimers [START_REF] Landau | Physique statistique[END_REF]. With the method of §2.3 of reference [START_REF] Castin | [END_REF] 15 we find for attractive interactions b 2 / √ 2 = (e u /2)[1 + erf (u 1/2 )] -1/4, where u ≡ β|E 0 (2)|; in the limit (35), u → 0 so that b 2 → 2 -3/2 , which reduces to the effect of quantum statistics, contrarily to the predictions of the model, and makes the temperature shift attributed to the dimers in figure 1 not meaningful 16 ! We interpret this as being due to the orthogonality of the diatomic unbound states with the dimer state, a fact not taken into account in the model; if the coupling constant g was changing from 0 + to 0 -, the dimer state would appear but a diatomic unbound state would disappear, at the energy scale |E 0 (2)| ≪ k B T , leading to a compensation in b 2 : the second cluster coefficient b 2 indeed varies continuously around g = 0. In our corrected model, we finally obtain the following implicit equation for the critical temperature of the liquid-gas transition in the limit (35):

Conclusion

We have here restricted our study of the unidimensional system of attractive (g < 0) bosons to the particular case of a big box, with a length L much larger that the diameter ℓ of the two body bound state. This means that, in a uniform distribution of the N bosons at the density ρ = N/L in the box, the mean interparticle spacing 1/ρ is much larger that the ground N-mer of diameter ξ ≃ 2 /(Nm|g|) and density profile [N/(4ξ)]/ cosh 2 [x/(2ξ)] at large N. In short, we are in the regime ρξ = ℓ/(2L) ≪ 1.

In this particular regime and at thermal equilibrium, we have shown that a liquid-gas transition occurs in the thermodynamic limit, at a temperature T c given to various orders of approximation by (16,37), and so high that the system lies in the non-degenerate regime ρλ c ≪ 1, where λ c is the corresponding thermal de Broglie wavelength. Over 14. In the liquid phase, one can apply Bogoliubov theory to the bright soliton if the Bogoliubov quasi-particles, of spectrum

ǫ k = |µ 0 |+ 2 k 2 /(2m), are a number δN = k [exp(βǫ k ) -1] -1 ≪ N, which imposes k B T ≪ |µ 0 | if ρλ ≪ 1, that is δN/N ≃ e -β|µ 0 | /(ρλ) ≪ 1: in the limit (35) for fixed Ng, the corresponding partition function obeys N -1 ln Z Bog = 1 3 β|µ 0 | + (ρλ) -1 s≥1 exp(-β|µ 0 |s)/ √
s, with a leading s = 1 term, for k B T ≪ |µ 0 |, that is indeed correctly given by (36); the contribution to N -1 ln Z Bog of the phase shift θ(k) = -4 atan(kℓ/N) due to the soliton [see above (9)] indeed tends to zero at large N, as -(4ℓ/L) exp(-β|µ 0

|)[1 + O(k B T/|µ 0 |)]/(Nρλ).
15. As in [25,26] we use a harmonic regulator of angular frequency ω, so that b 2 / √ 2 = 1/4 + lim ω→0 q≥0 (e -βǫqe -βǫ (0) q ) where (ǫ q ) q∈N is the spectrum of H = p 2 /m + mω 2 x 2 /4 + gδ(x) restricted to the even (bosonic) states, solving Γ( 1 4 -

ǫq 2 ω )/Γ( 3 4 - ǫq 2 ω ) = -2(2 ω/m) 1/2 /g [27] and ǫ (0) q = (2q + 1/2) ω is the corresponding g = 0 spectrum. If ω → 0 for a fixed g < 0, ǫ 0 → E 0 (2) and ǫ q≥1 = ω[2q -1/2 + 2∆(ǫ q ) + O(1/q)] ≥ 0, with ∆(ǫ) ≡ (1/π) arctan{[-ǫ/E 0 (2)] 1/2 }. If g > 0, ǫ q≥0 = ω[2q + 3/2 -2∆(ǫ q ) + O(1/q)]
≥ 0 with the same function ∆(ǫ). It remains to replace the sum over q by an integral over ǫ q , to obtain b 2 as given in the text for g < 0, and b 2 / √ 2 = e u 2 [1erf (u 1/2 )] -1 4 (in agreement with [28]) for g > 0.

16. On the contrary, the fraction of dimers predicted by the model in the gaseous phase agrees with the one

√ 2ρλ exp[-βE 0 (2)] deduced from b 2 .
an interval of temperature of width O(1/N)T c , the system changes from an essentially monoatomic uniform gaseous phase of density ρ [with a O(ρλ c ) n fraction of n-mers], to a liquid phase in the form of a single inhomogeneous drop of diameter ξ, surrounded by a small O(ρλ c ) 2 gaseous fraction which is mainly monoatomic and uniform.

On a theoretical point of view, our result seems to contradict a general argument due to Landau, that precludes the coexistence of two different phases in unidimensional systems with short range interactions, due to the fact that these phases have a tendency to experience arbitrarily high fragmentation and to form an emulsion [START_REF] Landau | Physique statistique[END_REF]. However the n-mers here do not constitute an ordinary extensive liquid with a density bounded from above: their energy E 0 (n) is not a linear function of the number of particles n at large n; it is rather a cubic function, so that the splitting of the liquid in mesoscopic or macroscopic fragments has a very high energy cost (for example (3/4)|E 0 (N)| if one cuts the N-mer into two halves), and it is not simply proportional to the number of fragments with an intensive proportionality factor as it was supposed in §163 of reference [START_REF] Landau | Physique statistique[END_REF].

At this stage, one may formulate the following objection: The lack of extensivity of the predicted liquid phase results from the Dirac delta function interaction potential, which should become unrealistic when N → +∞ since the liquid density then (presumably) diverges. In this limit, any realistic model for the interaction is expected to exhibit some finite range effect, such as a hard-core contribution. In the thermodynamic limit, however, it is the dimensionless quantity L/ℓ that matters, see equation ( 35), rather than the length ℓ. The objection can thus be circumvented by adjusting the interaction strength with the particle number N, thanks to a Feshbach resonance for ultracold atoms, in such a way that the central liquid density N/(4ξ) remains constant. This requires that N 2 /ℓ remains constant, so that the coupling constant g has to vanish as 1/N 2 . According to (35), the box length L then diverges as N 2 , and the ground liquid phase chemical potential µ 0 , the gaseous phase Fermi energy E F and the transition temperature T c all vanish as 1/N 2 . At the transition temperature, however, the quantum degeneracy parameter ρλ c remains fixed, and so does the entropy per particle in the gaseous phase, which shows that these low ∝ 1/N 2 temperatures can in principle be reached by adiabatic cooling.

The lack of extensivity of the liquid phase may have an important consequence on the metastability of the gaseous phase below the transition temperature. In general, it is known that a slowly cooled macroscopic gas in principle presents a very large delay to liquefaction, due to the extensive free energy barrier between the two phases, but in practice, this is avoided by local nucleation of droplets. Here this local mechanism may be inhibited by the strong non linearity of E 0 (n) with n. Let us indeed isolate in thought a small fraction of the gaseous phase with N ′ particles, with 1, ρℓ ≪ N ′ ≪ N, which is allowed by its homogeneity and the absence of long range interaction. As the liquefaction temperature T c of equation ( 14) is, within logarithmic corrections, proportional to the binding energy of a particle in the liquid that forms, it becomes obvious that the small fraction shall liquefy at a temperature T ′ c ∝ N ′2 much smaller than the one ∝ N 2 of the whole gas.

If it happens that this metastability constitutes a serious experimental handicap, a way out could be to measure the free energies of the two phases on both sides of the critical temperature, using a gradual heating of the already observed liquid phase at low temperature [3,4] and a slow cooling of the easy-to-prepare gaseous phase at high temperature, and then to realise that the free energies cross at T c . The free energy F(T ) can indeed be deduced from the mean energy E(T ) by integration of the relation E = ∂ β (βF), where β = 1/(k B T ). The mean energy, as it was shown in reference [29], is a known functional of the momentum distribution n(k) of the particles, that can be measured by time of flight after abrupt switch-off of the interactions thanks to a Feshbach resonance.

All this allows us to reasonably hope that an experimental test of our scenario is possible, maybe leading to the observation of a liquid-gas transition with ultracold atoms, which would constitute a "grande première". ) vois.at. and Z (k≤N/2) vois.N-mère , more precisely of their logarithm per particle N -1 ln Z vois. , as functions of ρλ for L/ℓ = 150 in the thermodynamic limit N → +∞. For the gaseous phase, the bracketing corresponds to the lower bound and the upper bound in (27) (red solid lines). For the liquid phase, the bracketing corresponds to the lower bound and the upper bound in (A.5) (black solid lines), or to the upper bound in (28) (black dashed line). The value of ρλ at the zeroth order estimate T (0) c of the transition temperature is shown by the vertical dotted line, at the intersection of the lower red and black solid lines. The improved upper bound (A.5), contrarily to the one of ( 28), excludes the existence of a phase transition at a temperature arbitrarily higher than T (0) c .

Figure 1 :

 1 Figure1: For N = 100 bosonic particles of mass m in one dimension, with an attractive interaction potential -|g|δ(x ix j ) in a box of length L = 100ℓ, where ℓ ≡ 2 2 /(m|g|) is the diameter of the two-body bound state: (a) free energy per particle F/N, in units of the Fermi energyE F = k B T F = 2 (πρ) 2 /(2m), and (b) mean fragmentation rate ν , as functions of the temperature in units of the Fermi temperature. Solid line: the numerical result obtained by generation of all possible configurations. Dotted line: when restricted to the two extreme configurations, the purely atomic one and the N-mer one. Dashed line: further including the gaseous mixtures (in arbitrary proportions) of atoms and dimers, as well as the minimally dissociated liquid phase of an atom plus a (N -1)-mer. In (c), the mean fraction of bosonic particles as a function of the size of the n-mers, n → n N n /N (for all n where it is > 10 -15 ), numerically obtained at the temperature ≈ 105.1T F where ν crosses (1 + 1/N)/2 [vertical line in (a) and (b)]. Note the logarithmic scale on the vertical axis in (b) and (c).

  since there is a transition from a N-mer of entropy k B O(ln N) to an atomic classical gas of entropy Nk B ln[e 3/2 /(ρλ)] + k B O(ln N).

7 .

 7 Its partition function is Z at-dim = Z at S (z), where z = exp[-βE 0 (2)]λ √ 2/L and the sum over the number of dimers s, S (z) = [N/2] s=0 z s N!/[s!(N -2s)!], can be expressed in terms of an hypergeometric function S

+ 2 - 3 / 2

 232 ρλ c + O(ρλ c ) 2 (37)

Figure A. 2 :

 2 Figure A.2: Bracketings of the partition functions Z (n≤N/2)vois.at. and Z (k≤N/2) vois.N-mère , more precisely of their logarithm per particle N -1 ln Z vois. , as functions of ρλ for L/ℓ = 150 in the thermodynamic limit N → +∞. For the gaseous phase, the bracketing corresponds to the lower bound and the upper bound in(27) (red solid lines). For the liquid phase, the bracketing corresponds to the lower bound and the upper bound in (A.5) (black solid lines), or to the upper bound in (28) (black dashed line). The value of ρλ at the zeroth order estimate T (0) c of the transition temperature is shown by the vertical dotted line, at the intersection of the lower red and black solid lines. The improved upper bound (A.5), contrarily to the one of (28), excludes the existence of a phase transition at a temperature arbitrarily higher than T (0) c .

Those elastic scattering events are for sure taken into account in the reference[13] in free space, but they do not impact the spectrum (2) since the various n-mers are then asymptotically free.

Contrarily to the three-dimensional case, the scattering of an incoming plane wave on a short range attractive potential in one dimension always tends to its total reflection at low energy, which is beyond the Born regime.

It may be useful to refer to footnote 14.

Acknowledgments

The group of Y.C. is also affiliated to IFRAF. One of us (C.H.) has received a Fulbright grant for his 1997-1998 stay in Paris, where he has accomplished most of this work, as well as National Science Foundation grants PHY-1316617 and 0844827; he also thanks the Sloan Foundation. M.O. was supported by grants from the Office of Naval Research (N0004-121-0400) and the National Science Foundation (PHY-1019197). We thank Christoph Weiss for his incitation to finalize this project.

Appendix A. An improved upper bound for the N-mer neighbourhood at very high temperature

In the main text, see equation (28), we have obtained the following upper bound of the partition function of the N-mer neighbourhood of maximal degree of dissociation k max < N/2:

This upper bound is useful close to the critical temperature but it is not in the limit of very high temperatures, in particular T ≫ T (0) c . When ρλ → 0 indeed, β|µ 0 | tends to zero, as well as Y, according to the equality

Then N -1 ln Z N-mère tends to zero as (ρλ) 2 , but θ diverges as 1/(ρλ), that is faster that N -1 ln Z at , which grows logarithmically only, as ln[e/(ρλ)]. The upper bound (A.1) is thus not sufficient to mathematically exclude the existence of a transition other than the one discussed in this paper, at a higher temperature.

The existence of such a competing transition can fortunately be excluded by a simple improvement of the upper bound in (A.1). One simply has to stop before the last step having led to (28), that is one keeps the finite sum over k:

where Γ is the incomplete Gamma function. The integral representation of this function leads to the following discussion, for N → +∞ at fixed θ and k max /N, which corresponds to the usual thermodynamic limit of section 4.3. If Nθ/k max < 1, an equivalent of the upper bound in (A.3) is indeed obtained by extending the sum over k to the entire N, and one recovers (28). If Nθ/k max > 1, an equivalent is obtained by the upper bound

where we have set k ′ = k maxk, we have used an upper bound of the type [START_REF] Diu | Physique Statistique[END_REF] and we have extended the summation over k ′ to N. The logarithm of the resulting geometric series is bounded and has a vanishing contribution in the thermodynamic limit. In other words, depending on the fact that Nθ/k max is less than or larger than one, the sum over k in (A.3) is dominated by its last or its first terms, at large N. We finally keep:

In the limit ρλ → 0, the upper bound in (A.5) diverges logarithmically only, with a prefactor k max /N < 1/2, that is at least two times less rapidly than N -1 ln Z at .

In the case of the maximally extended neighbourhoods of the atomic phase (n max = N/2) and of the N-mer (k max = N/2), we show, in figure A.2, the lower bound and the upper bound of N -1 ln Z vois. , where Z vois. is the partition function of each neighbourhood, as deduced from the bracketings ( 27), ( 28) and (A.5), as a function of ρλ and for L/ℓ = 150. Contrarily to (28), the improved bracketing allows to conclude without ambiguity that the liquid-gas transition takes place at a critical temperature close to T (0) c , even if one takes into account arbitrarily small values of ρλ, keeping in mind that the thermodynamically favoured phase is, at fixed temperature, the one with the largest partition function (that is, with the smallest free energy).


Flora Ponjou Tasse
email: flora.ponjou-tasse@cl.cam.ac.uk

Arnaud Emilien
email: arnaud.emillien@inria.fr

Marie-Paule Cani
email: marie-paule.cani@inria.fr

Stefanie Hahmann
email: stefanie.hahmann@inria.fr

Adrien Bernhardt
email: adrien.bernhardt@inria.fr

First Person Sketch-based Terrain Editing

Keywords: First person editing, terrain, sketch-based modelling, silhouettes I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-

We present a new method for first person sketch-based editing of terrain models. As in usual artistic pictures, the input sketch depicts complex silhouettes with cusps and T-junctions, which typically correspond to non-planar curves in 3D. After analysing depth constraints in the sketch based on perceptual cues, our method best matches the sketched silhouettes with silhouettes or ridges of the input terrain. A specific deformation algorithm is then applied to the terrain, enabling it to exactly match the sketch from the given perspective view, while insuring that none of the user-defined silhouettes is hidden by another part of the terrain. As our results show, this method enables users to easily personalize an existing terrain, while preserving its plausibility and style.

INTRODUCTION

Terrain is a key element in any outdoor environment. Applications of virtual terrain modelling are very common in movies, video games, advertisement and simulation frameworks such as flight simulators. Two of the most popular terrain modelling methods are procedural [START_REF] Fournier | Computer rendering of stochastic models[END_REF][START_REF] Miller | The definition and rendering of terrain maps[END_REF][START_REF] Lewis | Generalized stochastic subdivision[END_REF][START_REF] Musgrave | The synthesis and rendering of eroded fractal terrains[END_REF] and physics-based techniques [START_REF] Musgrave | The synthesis and rendering of eroded fractal terrains[END_REF][START_REF] Roudier | Landscapes synthesis achieved through erosion and deposition process simulation[END_REF][START_REF] Chiba | An erosion model based on velocity fields for the visual simulation of mountain scenery[END_REF][START_REF] Nagashima | Computer generation of eroded valley and mountain terrains[END_REF][START_REF] Neidhold | Interactive physically based fluid and erosion simulation[END_REF][START_REF] Kristof | Hydraulic erosion using smoothed particle hydrodynamics[END_REF]. The former are easy to implement and fast to compute, while the latter produce terrains with erosion effects and geologically sound features. However, the lack of controllability in these methods is a limitation for artists. Sketch-based or example-based terrains have been very popular recently in addressing these issues [START_REF] Cohen | Harold: A world made of drawings[END_REF][START_REF] Watanabe | A sketching interface for terrain modeling[END_REF][START_REF] Zhou | Terrain synthesis from digital elevation models[END_REF][START_REF] Gain | Terrain sketching[END_REF][START_REF] Hnaidi | Feature based terrain generation using diffusion equation[END_REF][START_REF] Tasse | Enhanced texture-based terrain synthesis on graphics hardware[END_REF][START_REF] Genevaux | Terrain generation using procedural models based on hydrology[END_REF]. However, many of these methods assume that the user sketch is drawn from a top view, which makes shape control from a viewpoint of interest very difficult. Others only handle a restricted category of mountains, with flat silhouettes. Lastly, terrains fully generated from sketches typically lack details. Dos Passos et al. [START_REF] Passos | Landsketch: A first person pointof-view example-based terrain modeling approach[END_REF] recently presented a promising approach where example-based terrain modelling and a first person point-of-view sketch are combined. However their method does not support local terrain editing and cannot handle typical terrain silhouettes with T-junctions. Moreover, terrain patches are often repeated which may spoil the plausibility of the results from other viewpoints. In this work, we address the problem of intuitive shape control of a terrain from a first person viewpoint, while generating a detailed output, plausible from anywhere. To achieve the intuitive shape control goal, we stick to the sketch-based approach, but allow the user to input complex silhouettes with cusps and T-junctions, as those typically used to represent terrains (see Figure 1). To get plausible, detailed results from anywhere, we focus on editing an existing terrain rather than starting from scratch. This approach captures the coherent small details from the existing terrain, while avoiding the patch blending and repetition problems that are typical of example-based methods. The use of an existing terrain also enables matches of sketched silhouettes with plausible, non planar curves on the terrain.

In practice, the user edits the input terrain by over-sketching it from a first person viewpoint. The user strokes, forming a graph of curves with T-junctions, represent the desired silhouettes for the ter-rain. The input terrain is then deformed such that its silhouettes exactly match the strokes in the current perspective view. This means that each stroke segment is to be some silhouette of the output terrain, and that no other part of the deformed terrain should hide them. Previous sketch-based modelling methods have successfully use feature curves to deform surfaces [START_REF] Singh | Wires: A geometric deformation technique[END_REF][START_REF] Zimmermann | Silsketch: Automated sketch-based editing of surface meshes[END_REF]. Our work explores the use of terrain features for sketch-based terrain editing.

First, we order the sketched strokes by inferring their relative depth from the height of their end-points and from the T-junctions detected in the sketch. Next, features of the input terrain such as silhouette edges and ridges are assigned to each stroke and extended if necessary, to cover the length of the stroke. This assignment is the solution of a minimization problem expressing the similarity between a terrain feature and a stroke in the drawing plane, and the amount of deformation caused by their matching. The selected features then become constraints for an iterative diffusion-based terrain deformation method. Our main contributions are:

• An algorithm for ordering strokes in a complex, perspective sketch with respect to their distance from the camera.

• A method for matching terrain features with user-specified silhouettes, drawn from a given first-person viewpoint.

• A deformation method for matching silhouette constraints while preventing them from being hidden by other parts of the terrain.

Related work is discussed in Section 2. Then we present an overview of our solution in Section 3. This is followed by a description of stroke ordering in Section 4, generation of feature constraints in Section 5 and terrain deformation in Section 6. We discuss results in Section 7 before concluding.

RELATED WORK

Most terrain modelling systems use one or a combination of the following: procedural terrain generation, physics-based simulation, sketch-based or example-based methods. See [START_REF] Natali | Modeling terrains and subsurface geology[END_REF] for a detailed survey.

Procedural terrain modelling methods are based on the fact that terrains are self-similar, i.e. statistically invariant under magnification [START_REF] Mandelbrot | The fractal geometry of nature[END_REF]. These methods are the popular choice for landscape modelling due to their easy implementation and efficient computation. They mainly consist of pseudo-randomly editing height values on a flat terrain by using either adaptive subdivision [START_REF] Fournier | Computer rendering of stochastic models[END_REF][START_REF] Miller | The definition and rendering of terrain maps[END_REF][START_REF] Lewis | Generalized stochastic subdivision[END_REF] or noise [START_REF] Miller | The definition and rendering of terrain maps[END_REF][START_REF] Musgrave | The synthesis and rendering of eroded fractal terrains[END_REF]. Adaptive subdivision progressively increases the level of detail of the terrain by iteratively interpolating between neighbouring points and displacing the new intermediate points by increasingly smaller random values. Noise synthesis techniques are often preferred because they offer better control. Superposing scaleddown copies of a band-limited, stochastic noise function generates noise-based terrains. For more information on fractal terrain generation methods, see Ebert et al. [START_REF] Ebert | Texturing and Modeling: A Procedural Approach[END_REF]. Fractal-based approaches can generate a wide range of large terrains with unlimited level of details. However, they are limited by the lack of user control or non-intuitive parameter manipulation, and the absence of erosion effects such as drainage patterns. To address the last issue, fractal terrains can be improved using physics-based erosion simulation [START_REF] Musgrave | The synthesis and rendering of eroded fractal terrains[END_REF][START_REF] Roudier | Landscapes synthesis achieved through erosion and deposition process simulation[END_REF][START_REF] Chiba | An erosion model based on velocity fields for the visual simulation of mountain scenery[END_REF][START_REF] Nagashima | Computer generation of eroded valley and mountain terrains[END_REF][START_REF] Neidhold | Interactive physically based fluid and erosion simulation[END_REF][START_REF] Kristof | Hydraulic erosion using smoothed particle hydrodynamics[END_REF]. Alternatively, river network generation can be incorporated in the procedural method [START_REF] Kelley | Terrain simulation using a model of stream erosion[END_REF][START_REF] Genevaux | Terrain generation using procedural models based on hydrology[END_REF]. In particular, Genevaux et al. [START_REF] Genevaux | Terrain generation using procedural models based on hydrology[END_REF] create procedural terrains from a hydrographically and geomorphologically consistent river drainage network, generated from a top-view sketch. However, this method only captures terrains resulting from hydraulic erosion, and there is no mechanism for controlling their silhouettes from a first person viewpoint.

Physically-based techniques generate artificial terrains by simulating erosion effects over some input 3D model. Musgrave et al. [START_REF] Musgrave | The synthesis and rendering of eroded fractal terrains[END_REF] present the first methods for thermal and hydraulic erosion based on geomorphology rules. Roudier et al. [START_REF] Roudier | Landscapes synthesis achieved through erosion and deposition process simulation[END_REF] introduce a hydraulic erosion simulation that uses different materials at various locations resulting in different interactions with water. Chiba et al. [START_REF] Chiba | An erosion model based on velocity fields for the visual simulation of mountain scenery[END_REF] generate a vector field of water flow that then controls how sediment moves during erosion. This process produces hierarchical ridge structures and thus enhances realism. Nagashima [START_REF] Nagashima | Computer generation of eroded valley and mountain terrains[END_REF] combines thermal and fluvial erosion by using a river network pregenerated with a 2D fractal function. Neidhold et al. [START_REF] Neidhold | Interactive physically based fluid and erosion simulation[END_REF] present a physically correct simulation based on fluids dynamics and interactive methods that enable the input of global parameters such as rainfall or local water sources. Kristof et al. [START_REF] Kristof | Hydraulic erosion using smoothed particle hydrodynamics[END_REF] propose fast hydraulic erosion based on Smooth Particle Hydrodynamics. The main drawback of all these methods is that they only allow indirect user-control through trial and error, requiring a good understanding of the underlying physics, time and efforts to get the expected results.

Sketching interfaces and more generally feature-based editing have been increasingly popular for terrain modelling. These methods can be combined with some input terrain data to generate terrains with plausible details.

Cohen et al. [START_REF] Cohen | Harold: A world made of drawings[END_REF] and Watanabe et al. [START_REF] Watanabe | A sketching interface for terrain modeling[END_REF] present the first terrain modelling interfaces that take as input a 2D silhouette stroke directly drawn on a 3D terrain model. They only handle a single silhouette stroke, interpreted as a flat feature curve. McCrae and Singh [START_REF] Mccrae | Sketch-based path design[END_REF] use stroke-based input to create paths which deform terrains. However user strokes are interpreted as path layouts and not as terrain silhouettes. Multi-grid diffusion methods enable generation of terrains that simultaneously match several feature curves, either drawn from a top view [START_REF] Hnaidi | Feature based terrain generation using diffusion equation[END_REF] or from an arbitrary viewpoint [START_REF] Bernhardt | Realtime Terrain Modeling using CPU-GPU Coupled Computation[END_REF]. The main limitation is that generated terrains typically lack realistic details.

In contrast, Zhou et al. [START_REF] Zhou | Terrain synthesis from digital elevation models[END_REF] use features (actually, sketch maps painted from above) to drive patch-based terrain synthesis from real terrain data. Closer to our concerns, Gain et al. [START_REF] Gain | Terrain sketching[END_REF] deform an existing terrain from a set of sketched silhouettes and boundary curves. The algorithm deforms the terrain based on the relative distance to the feature-curves in their region of influence, and on wavelet noise to add details to the silhouettes. In this work we rather use a diffusion-based deformation method to propagate feature constraints, avoiding the need for boundary curves. Lastly, Tasse et al. [START_REF] Tasse | Enhanced texture-based terrain synthesis on graphics hardware[END_REF] present a distributed texture-based terrain generation method that re-uses the same sketching interface. Unfortunately, all these methods interpret each sketched silhouette as a planar feature curve, which reduces the realism of the result. Dos Passos et al. [START_REF] Passos | Landsketch: A first person pointof-view example-based terrain modeling approach[END_REF] propose a different approach to address this issue. Given a set of sketched strokes drawn from a first person point-of-view, copies of an example terrain are combined such that the silhouettes of the resulting terrain match the strokes. This gives a realistic, varying depth to silhouettes. To achieve this, the algorithm assumes each stroke represents a terrain silhouette. A stroke is matched with a portion of a silhouette, selected from a set of silhouettes viewed from several standing viewpoints around the example terrain. Terrain slices representing portions of matched silhouette are cut from the example terrain and then combined through a weighted sum to produce a smooth terrain. A drawback of this method is that it does not handle complex sketches with T-junctions, which are common in landscape drawings. Moreover, the matching process may select the same silhouette portions for different strokes, thus producing unrealistic repeating patterns in the final result. Finally, the weighted sum function used for merging may fail to remove the boundary seams produced by combining different terrain slices. In this work, we address these issues by presenting a sketch-based method that handles T-junctions in complex sketches and deforms an input terrain to match the sketch rather than copypasting parts of it.

OVERVIEW

Let us describe our processing pipeline. As in many terrain modelling and rendering methods, our terrains are represented by a height field, namely a greyscale image storing elevation values. This representation cannot emulate features such as overhangs and caves, but it is the most prevalent format in terrain generation because of its simplicity and efficient use of storage space. For rendering purposes and silhouette detection, a 3D triangular mesh is constructed from the height field by connecting adjacent terrain points (x, y, altitude(x, y)). Users are able to navigate on a 3D rendering of the existing terrain, possibly flat, with a first-person camera always at a standing viewpoint. A sketch is created by drawing one or multiple strokes from the same camera position. The drawn strokes represent silhouettes that should be visible from that position. Our main goal is to deform the terrain such that these user constraints are respected. The following requirements should be satisfied:

• Every sketched stroke should be a terrain silhouette, in the current perspective view from the first-person camera viewpoint.

• Each of these terrain silhouettes should be visible, i.e. not hidden by any other part of the terrain.

• The deformed terrain should not have artifacts nor contain unrealistic deformations, from any other viewpoint.

Our solution consists of five main steps, illustrated in Figure 2: 1. We order strokes according to their depth, from front to back with respect to the camera position. This order is used when we generate constraints for terrain deformation, so that a curve constraint is not occluded by another, when viewed from the first-person viewpoint.

2. Terrain features such as silhouettes and ridges are detected. Deforming existing terrain features to match the desired silhouettes results in a more realistic terrain since no extra features are added and thus, the nature of the existing terrain is best preserved.

3. For each stroke, we select a terrain feature that will be deformed to fit the stroke, when seen from the camera position. These deformed features represent the positional constraints that we use in the diffusion-based terrain deformation. A key idea of our framework is the expression of this feature selection step as an energy minimization problem, in which we penalize features with large altitude differences compared to their corresponding stroke as well as features that would result in too large deformations.

4. We use a multi-grid Poisson solver for diffusion-based terrain deformation. It solves for altitude differences instead of absolute terrain positions, thus preserving the small-scale features of the input terrain.

5. After terrain deformation, other parts of the terrain may hide the user-specified silhouettes. To address this issue, we run the following iterative process: we detect terrain silhouettes that do not fit any user stroke and yet hide one of the sketched silhouettes. Extra deformation constraints are constructed to enforce lowering these protruding silhouettes until the usersketched silhouettes are no longer occluded. The terrain is deformed with a combination of previous constraints and the newly constructed constraints. We repeat this process until there is no longer protruding silhouette. In this section, we explain how depth ordering of silhouette strokes is extracted from the user sketch.

The different silhouette strokes in the input sketch first need to be ordered, in terms of relative depth from the camera viewpoint. This will enable us to ensure, when they are matched with features, that they will not be hidden by other parts of the terrain. Our approach to do so is based on two observations:

• If, in the viewing plane, a silhouette lies above another, it obviously corresponds to a mountain A farther away from the viewpoint than the other mountain B. Otherwise A would hide B. Using height coverage for ordering them in depth is however not sufficient, since some strokes may overlap in height, as for the green and blue strokes in Figure 3.

• Furthermore, the terrain being a height field, the projection of each stroke onto the horizon (x-axis of the viewing plane) is injective (no more than one height value per point). These two observations allow us to solve the relative stroke ordering problem thanks to a new sweeping algorithm (see Figure 3): We consider the projections of all the strokes onto the horizontal x axis (depicted in the bottom part of the Figure ) and sweep from left to right, examining the extremities (starting and endpoints in sweeping direction) and junction points of the silhouette strokes. While doing so, we label the strokes' extremities and the junction points in the following way: an extremity q s of stroke s is a Tjunction if its closest distance to another stroke r is smaller than a threshold. An endpoint q s is labelled (occluded-by, r) if the oriented angle, measured counterclockwise, between the tangent 1 of s 1 Strokes are always oriented clockwise. Hence, stroke tangents are independent of the direction in which the stroke was sketched. When labelling a starting point q s as T-junction, we flip its tangent.

at q s and the tangent of r at q s , (t s ,t r ) < 180 • . This indicates that s is occluded by, and thus behind, r. Otherwise, s is in front of r and we label q s as (in-front-of, r).

If a stroke s has no T-junctions, then it is behind a stroke r either if the projection of s completely contains the projection of r or if the smallest height value of s's endpoints is larger than the smallest height value of r's endpoints.

While scanning the sketch from left to right, we insert each stroke in a sorting structure, at a relative depth position determined by the cues above. This results in a relative ordering of the user strokes.

POSITIONING STROKES IN WORLD SPACE

The key idea of our approach is to create a 3D terrain that matches the user drawing, by deforming an existing one. More precisely, we deform the features of the existing terrain, like its ridgelines, to match the user silhouette strokes. Because a terrain has many features, we first have to compute to which one of them it is the most appropriate to apply a deformation. In this section, we detail how we compute the set of terrain features (Section 5.1), how we affect one of them to each of the user strokes (Section 5.2) and we present a feature completion algorithm that infers the hidden parts of the silhouettes, enabling a more realistic terrain deformation result (Section 5.3).

Feature detection: silhouettes and ridgelines

Silhouette detection on the existing terrain is based on a common and naive algorithm for computing the exact silhouettes of a 3D mesh. Silhouette edges are detected by finding all visible edges shared by a front face and a back face in the current perspective view. Neighbouring silhouette edges are then linked to form long silhouette curves.

Ridge detection is based on the profile-recognition and polygonbreaking algorithm (PPA) by Chang et al. [START_REF] Chang | A visual basic program for ridge axis picking on dem data using the profile-recognition and polygon-breaking algorithm[END_REF]. The PPA algorithm marks each terrain point that is likely to be on a ridge line, based on the point height profile. Segments, forming a cyclic graph, connect adjacent candidate points. Polygon-breaking repeatedly deletes the lowest segment in a cycle until the graph is acyclic. Finally, the branches on the produced tree structure are reduced and smoothed. The result is a graph where nodes are end points or branch points connected by curvilinear ridgelines. An improvement of the PPA algorithm connects all the terrain points into a graph using a heightbased or curvature-based weighting and computes the minimum spanning tree of that graph [START_REF] Bangay | Minimum spanning trees for valley and ridge characterization in digital elevation maps[END_REF]. Because we are mainly concerned with performance and detection of large-scale ridges, we simply connect candidate terrain points as in the original PPA algorithm and replace the polygon-breaking with a minimum spanning forest algorithm.

Stroke -Feature matching

In this section, we discuss a method for determining, for each stroke, the terrain features which can be used to construct deformation curve constraints. Viewed from the first person camera, these curve constraints should match the user-sketched strokes. To achieve this, we first construct a features priority list for each stroke and then select features for each priority list such that the sum of their associated cost is minimized.

Features priority list per stroke

For a stroke s, we project all terrain features on the sketching plane (i.e. we use the 2D projection of the feature from the first-person viewpoint) and select feature curves that satisfy the following condition: the x interval they cover matches the one of the stroke s. We deform the selected feature curves, and if necessary extend their endpoints, such that viewed from the camera position, they cover the length of s. This deformation is simply achieved by displacing the feature curve points according to their projection on the 2D stroke in the sketching plane, and their distance to the camera position. Let f be a terrain feature and f p its projection on the stroke plane. If a portion of f p is below or above s, f is truncated to that portion. Moreover, for each point q ∈ f , its altitude is modified as follows:

q.z = q.z + k * ||q p -q s p || * ||q -eye|| ||q p -eye||
where eye is the camera position, k = -1 if f p is below s and k = 1 otherwise, q p the projection of q on the stroke plane, and q s p the intersection of s and the vertical line passing at q p . We used this deformed version of the feature to associate the following cost E( f , s) to each feature f with respect to stroke s:

E = E dissimilarity + E de f ormation + E sampling + E extension (1) E dissimilarity ( f ) = w 1 curvelength( f p ) * f p h f p dt E de f ormation ( f ) = w 2 curvelength( f ) * f h f dt E sampling ( f ) = w 3 * longestedgelength( f ) max g∈prioritylist(s) longestedgelength(g) E extension ( f ) = w 4 * extendedcurvelength( f ) curvelength( f )
where, w i are weights, f p is the projection of f on the stroke plane, h f is the altitude difference between f and f 's projection on the terrain, and h f p is the altitude difference between f p and the stroke s. E dissimilarity represents the dissimilarity between f and s, E de f ormation expresses the amount of along f , E sampling penalizes features with long edges and E extension penalizes features that were extended to fully cover s when viewed from the camera position. All the results shown here were generated with w 1 = w 2 = w 3 = w 4 = 1.0. All features are sorted in a priority list according to their cost. Figure 4 illustrates this process for a single stroke (in this simple case, the feature of minimal cost is selected).

Energy minimization

The goal here is the selection of a feature curve f from the priority list of each stroke s i , to construct deformation constraints for terrain deformation. In addition to the feature order within the different priority lists, we need to take into account the depth ordering for silhouette strokes computed in Section 4.

Therefore, this selection process can be seen as a minimization problem. Let S = {s i : i = 1, ..., n} be the stroke list (ordered by depth) and f i denote a feature in the priority list L(s i ) = { f i k : k = 1, ..., m i } for a stroke s i . We are looking for { f i : i ∈ 1...n} such that f i < f j if i < j and ∑ E( f i ) is minimized. Here, f i < f j means that f i should not be occluded by f j , so that all deformation curve constraints are visible from the first person viewpoint. We process the ordered stroke list from front to back, and after each stroke, we remove from the priority list of the next strokes, features that will be occluded if selected. We chose to process strokes from front to back for two main reasons. Firstly, strokes that are closest to the eye are processed first and due to E de f ormation , the algorithm attempts to select constraints that will minimize the terrain deformation. Thus, features closer to the eye are more likely to be selected. Secondly, if all the features of interest for a given stroke s i were already selected, and therefore its priority list was empty, an arbitrary curve on the terrain would be used instead. If this ever occurs, we prefer it to be for background silhouettes. In practice, feature selections that cause any stroke to have an empty priority list are penalized with a very high cost. Thus, a configuration that guarantees a non-empty priority list for each stroke is always selected, if it exists. If no such configuration exists and s i has an empty priority list, we automatically compute a 3D embedding of the 2D stroke s i and use the resulting curve as a deformation constraint. To easily compute this 3D embedding, we take the two strokes lying just in front and just behind s i . Then we place s i halfway between the terrain features assigned to these two strokes. If there is no stroke restricted to lie behind s i , we place it behind the furthest stroke from the viewpoint. If there is no stroke restricted to lie in front of s i , we place it in front of the closest stroke to the viewpoint. With this approach, each stroke is represented by a deformation constraint even if it was not matched to a terrain feature during the energy minimization.

The energy minimization problem we have described so far is NP-hard. We use a branch-and-bound algorithm to efficiently discard all partial solutions that have a cost higher than the current best cost, without having to explore the whole solution tree. The branch-and-bound algorithm consists of two steps: a branching step and a bounding step. The branching step consists of exploring possible choices for s i+1 once we have made a feature selection for s i . In other words, we split the node (s i , f i ) into multiple nodes (s i+1 , f i+1 k ), where f i+1 k are features in the priority list of s i+1 . The bounding step allows the algorithm to stop exploring a partial solution if the total cost of features in the solution is higher than the cost of the best solution found so far.

It is possible for a feature to be the first choice in the priority lists for two or more strokes. To handle this, when exploring a possible solution, a feature curve assigned to a stroke is no longer considered for subsequent strokes. Our branch and bound algorithm will explore other solutions with the feature curve assigned to different strokes as long these solutions are guaranteed to have a smaller cost than the current best solution.

Stroke in world space

The previous minimization gives us, for each stroke s, an associated terrain feature f . However, the stroke s has its points in screen space, whereas the points of f are in the world space. Our goal is to place the stroke in the world space, in order to deduce terrain constraints, i.e. find the distance of their projection from the camera.

For each point of the stroke q s = (x s , y s ), we check if there exists a feature point q f whose projection on screen q p = P(q f ) = (x p , y p ) has the same x-coordinate as q s , i.e. x s = x p . If this point exists, we project the stroke point on the world space, using the distance of q f from the camera as a depth value.

The possible undetermined points depth, at the stroke borders, are set in world space to follow the stroke tangent, in the world space.

Completing selected 3D features

Using user-specified endpoints of an occluded stroke during the generation of deformation constraints would create silhouettes that appear to start exactly at these endpoints. This can look quite unnatural when viewed from a different position than the first person camera position used for sketching: indeed, the endpoint of the occluded stroke (a junction) is typically above the terrain and thus, a sharp deformation will be created at that point.

We address this problem by simply extending 3D features assigned to strokes at both endpoints along their tangents, until they reach the surface of the terrain. An example of this feature completion is presented in Figure 5. More sophisticated contour completion methods such as the one presented in SmoothSketch [START_REF] Karpenko | Smoothsketch: 3d free-form shapes from complex sketches[END_REF] could alternatively be used, but this simple method was sufficient in our case, and is proposed as an optional step in the editing process.

TERRAIN DEFORMATION

In the previous section, we analysed terrain features and used them to position the strokes in the world space. We present in this section how we use them as constraints to deform the existing terrain.

Diffusion-based equation solver

Our deformation algorithm relies on iterative diffusion of displacement constraints, which are computed from the 3D strokes positioned in the world space. The diffusion method, first introduced by the authors in [START_REF] Emilien | Interactive procedural modeling of coherent waterfall scenes[END_REF], consists in computing the difference of the curve height and the terrain height H , and to diffuse these differences (instead of absolute height values) using a multi-grid Poisson solver similar to the one in [START_REF] Hnaidi | Feature based terrain generation using diffusion equation[END_REF].

More precisely, for each point p = (x, y, z) of the stroke in the world space, we compute δ = z -H (x, y), and set it as a displacement constraint. The constraints are rasterised on a grid, whose resolution is equal to the terrain resolution. After having set the constraints of all strokes, we perform the diffusion, which gives the displacement map M .

The displacement is finally applied on the terrain height field H , whose feature line silhouettes are now matching the user strokes, when seen from the first-person viewpoint used for sketching. The deformation only consists of adding the two heights, H ′ (x, y) = H (x, y) + M (x, y), where H ′ is the resulting terrain. Because height differences are propagated, instead of absolute heights, the terrain preserves fine-scale details during the deformation.

Lowering protruding silhouettes

After deformation, the user-defined silhouettes may be hidden by other parts of the terrain. To address this issue, we detect the unwanted protruding silhouettes and constrain them to a lower position so that the user-defined silhouettes become visible.

Detecting most protruding silhouette edges

First, all visible silhouettes are detected, with the algorithm discussed in Section 5.1. These silhouettes are projected onto the sketching plane. Let s be a silhouette of the deformed landscape, inherited from the example terrain. The mountain of silhouette s hides a user-specified silhouette g if s is closer to the camera than g and the projection s p of s in the sketching plane has a higher altitude than g p , the projection of g. In this case, s is an unwanted protruding silhouette. Determining how much s should be lowered is done as follows: Let h be the maximum height difference between s and a silhouette g hidden by s. h is the minimum altitude by which s should be lowered to ensure the silhouettes it hides become visible. Our solution is simply to uniformly lower s by an offset h. This method is applied to all unwanted protruding silhouettes and we use the set of lowered silhouettes to form new deformation constraints.

Updating deformation constraints

The new deformation constraints from the lowered protruding silhouettes are added to the set of constraints associated to the sketched silhouettes, and the terrain is deformed once again using the method of section 6.1. This operation maintains the userspecified silhouettes while lowering areas around the unwanted protruding silhouettes, so that user specifications are satisfied.

The process of detecting protruding silhouettes and using this information to further constrain the terrain is repeated until protruding silhouettes are no longer detected. In practice, a single iteration is usually sufficient to make all user-specified silhouette strokes visible.

RESULTS

Validation examples : The examples below and the joined video illustrate the results of our method in a variety of cases. In particular, Figure 6 shows editing of a terrain with a complex sketch containing 5 T-junctions. Our method is also able to handle complex mountains where ridges are not as well-defined as they are on smooth landscapes. An example of this is presented in Figure 7. Our proposed approach differs from other sketch-based methods in that non-planar silhouettes can be generated from planar usersketched strokes. This is illustrated in Figure 8. Moreover, the method is robust enough to support terrains with few or no features, as shown in the example given in Figure 9. This complex sketchbased editing framework comes with interactive rates, as illustrated in the attached video, which makes it a very attractive alternative to other terrain generation/editing techniques discussed in Section 2.

User tests: We performed an informal user test with two experienced computer artists. The system was briefly introduced to the users, who had no prior knowledge about it. They were then asked to draw sketches to deform existing terrains. Both of them reported that our system was very easy to take in hand and use, and where able to quickly create new sceneries. Their feedback indicates that the approach is original, and seems a promising way to create a scene that matches their artistic intend. However, they had some difficulties for predicting the result of their sketch. This could have been improved by manually tuning the feature matching weights, which was not among their skills. These first users also asked for the ability to move within the scene and edit the terrain from multiple viewpoints. Note that this is not a technical issue: we simply did not implement it since this editing mode was already provided by previous methods. Lastly, the users insisted on the natural aspect of the resulting terrain, and noted that it matched their sketches in the expected way.

Limitations: Although our system succeeds in matching a complex user-sketch through a natural deformation of the terrain, based on its existing features, the lack of predictability of the strokefeature solver may be a problem. A full user study would be useful to understand the user intent when sketching over the existing terrain, and identify the ranges of weight values allowing us to better match this intent. We could also improve our matching method using extra error functions.

Another limitation comes from our deformation solver. The diffusion-based deformation method sometimes creates small declivities around the extremity of a constraint curve, when the slope of the curve is high and the extremity is located on the terrain: in this case, the terrain locally inflates, except at this end-point where the deformation is zero, which causes the problem. Using an inverse distance to deform a terrain [START_REF] Jenny | Interactive local terrain deformation inspired by hand-painted panoramas[END_REF] does not work either, because of our use of curves as constraints. Future work still needs to be done on terrain deformation, especially for curve-based deformations.

CONCLUSION

We presented the first sketch-based modelling method enabling the deformation of a terrain from a single viewpoint. The user sketches a few silhouette strokes forming a graph with T-junctions, similar to the silhouette representations used in artistic terrain sketching. A key feature of our method is that sketched silhouettes are matched with existing terrain features: this enables our technique to both match silhouette strokes with a non-planar curve, and produce a deformation that does not spoil plausibility, since the structure of ridges and valleys typically remains unchanged. This work could easily be extended to a multiple-view editing interface, where the user can move over the terrain and iteratively edit it from different points of view, while keeping his other sketches as constraints.

Figure 1 :

 1 Figure 1: A typical artist sketch (top left), is used to edit an existing terrain (right). Results are shown on the second row from the same two viewpoints. Note the complex silhouettes with T-junctions, matched to features of the input terrain. The bottom image shows a rendering of the resulting terrain, from a closer viewpoint.

Figure 2 :

 2 Figure 2: Overview of our terrain editing framework. In (b), stroke colour indicates stroke ordering: blue indicates that a stroke is closer to the camera position and red indicates that it is the furthest.(c) illustrates detected terrain features in white and (d) shows the subset of terrain features that have been assigned to user strokes. In (f) the terrain features are deformed so that they match the strokes from the user viewpoint. The final result in (h) is obtained after removing some residual artifacts.

Figure 3 :

 3 Figure 3: An input sketch (top) and the different steps of the sweeping algorithm used for scanning the sketch, labelling T-junctions and ordering strokes (bottom). As a result, stroke 3 is detected to be in front of stroke 2, which is itself in front of stroke 1.

Figure 4 :

 4 Figure 4: Computing possible features to match with a user stroke. Images (a) and (d) show the terrain from the first person viewpoint used for editing, while image (b) and (c) use a higher viewpoint to better show features on the input terrain. Feature colour indicates cost: blue for the lowest cost and red for the highest.

  (a) User input (b) Matched features (c) Extend the matched features (d) Resulting terrain

Figure 5 :

 5 Figure 5: Completing selected features: after matching 2D strokes to terrain features, we extend these features until they reach the surface of the terrain, to ensure a smooth transition from specified silhouettes to the terrain.

  (a) User input (b) Existing terrain (c) Deformed terrain (d) Result viewed from a different point

Figure 8 :

 8 Figure 8: Terrain editing produces non-planar silhouettes in the output, from 2D planar strokes.

Figure 9 :

 9 Figure9: Adding deformation constraints automatically: the stroke furthest away from the user did not have an assigned feature to it and so one was automatically generated, and positioned on a plane orthogonal to the view direction, such that stroke ordering is respected.

ACKNOWLEDGEMENTS

This work was conducted during an internship of Flora Ponjou Tasse at Inria Rhône-Alpes in Grenoble. It was partly supported by the ERC advanced grant EXPRESSIVE.