Andrii Kleshchonok

Geneviève Fleury

Jean-Louis Pichard

Scanning Gate Microscopy of Kondo Dots: Fabry-Pérot Interferences and Thermally Induced Rings

Keywords: numbers: 07, 79, -v, 72, 10, -d 73, 63, Kv 72, 15, Qm

We study the conductance of an electron interferometer formed in a two dimensional electron gas between a nanostructured quantum contact and the charged tip of a scanning gate microscope. Measuring the conductance as a function of the tip position, thermally induced rings may be observed in addition to Fabry-Pérot interference fringes spaced by half the Fermi wavelength. If the contact is made of a quantum dot opened in the middle of a Kondo valley, we show how the location of the rings allows to measure by electron interferometry the magnetic moment of the dot above the Kondo temperature.

Scanning gate microscopy (SGM) is a new tool which allows to probe by electron interferometry [1] the properties of nanostructures created in a two-dimensional electron gas (2DEG). The nanostructures are made with charged gates deposited on the surface of a semiconductor heterostructure, allowing to divide the 2DEG beneath the surface in two parts connected via a more or less simple contact region. This can be a quantum point contact [2,3] (QPC), a quantum dot [4][5][START_REF] Grobis | Handbook of Magnetism and Magnetic Materials[END_REF], a double dot setup [START_REF] Van Der Vaart | [END_REF] or more complex nanostructures. With five gates, the contact between the right and left leads (left and right parts of the 2DEG) can be made of a quantum dot with a tunable gate (see Fig. 1(a)). With the charged tip of an atomic force microscope above the surface of the heterostructure, a depletion region can be capacitively induced in the 2DEG below the surface at a distance r from the contact. SGM consists in studying the conductance g of the electron interferometer formed by the contact and the depletion region. Scanning the tip outside the contact, one can record SGM images giving g as a function of the tip position. These images exhibit Fabry-Pérot interference fringes spaced by λ F /2, as observed by Topinka et al [8] for a QPC opened on its first conductance plateau. Using high-mobility 2DEGs, the SGM images of a single QPC have been more systematically investigated later at low temperatures in Refs. [9][10][11][12] for different points of a conductance plateau as well as between plateaus. This has led to revisit the theory of electron interferometers [13][14][15][16][17][18] which include a quantum point contact. Recently, the spacing between the interference fringes at a certain distance from the contact was found [12] to differ by more than 50% from the expected value λ F /2 when the QPC is biased on the second, third and fourth conductance plateaus. This difference gives rise to a large ring of radius ≈ 1µm visible in the SGM images of Ref. [12] where different scenarii were proposed for explaining this unexpected ring. This leads us to study if interference mechanisms other than those responsible for the λ F /2-oscillations can occur in the limit where the electron motion is purely ballistic between the contact V I FIG. 1: (Color online) (a): Metallic gates (yellow) create depletion regions (green) in the 2DEG (blue) beneath the surface. This makes a dot with a gate in the contact between left and right 2DEGs, while a charged tip yields a scannable depletion region near the contact. The SGM images give the conductance as a function of the tip position. (b): Two semiinfinite square lattices are contacted via an Anderson impurity (site I of coordinates (0, 0), Hubbard repulsion U , potential VI, coupling tc). Adding a potential V at a site of coordinates (x, y) gives rise to an electron interferometer of size r. and the tip. We show in this letter that this can indeed occur for a contact characterized by a series of transmission peaks, if it is opened between two peaks. Moreover, for a dot with an odd number of electrons and biased in the middle of the Kondo valley, this opens the possibility to measure by electron interferometry the magnetic moment [4][5][START_REF] Grobis | Handbook of Magnetism and Magnetic Materials[END_REF] induced by electron-electron interactions above the Kondo temperature.

In the Coulomb blockade regime, the gate voltage can be tuned for having an odd number of electrons in a lithographically defined quantum dot. This gives rise to an unpaired spin. Our study is based on the model routinely used [START_REF] Grobis | Handbook of Magnetism and Magnetic Materials[END_REF] for describing the Kondo effect due to this unpaired spin: An Anderson impurity coupled to two semi-infinite square lattices, the many body effects coming from the presence of an Hubbard repulsion U in the contact (see Fig. 1(b)). The Hamiltonian without tip reads H 0 = H I + H c + H leads where

H I = V I σ n Iσ + U n I↑ n I↓ (1)
H c = t c σ c † Iσ c (0,1)σ + c † Iσ c (0,-1)σ + H.C. (2)
H leads = iσ -4tn iσ + t j c † iσ c jσ + H.C. (3)
c iσ (c † iσ ) is the destruction (creation) operator of an electron of spin σ at site i and n iσ = c † iσ c iσ . H I describes the Anderson impurity located at the site I of coordinates (0, 0) making the contact and H c the coupling to the leads (hopping t c between I and the two neighboring sites of the two leads). H leads describes two semi-infinite square lattices making the right and left leads (nearest neighbor hopping t). The energy scale is defined by taking t = -1 and site potentials equal to -4t yield energy bands [0, 8] for the conduction electrons of the leads. Hereafter, we study the continuum limit and consider only small energies E. To this model for the contact, we add a term H tip = σ V n Tσ , assuming that the depletion region induced by the charged tip modifies only the potential of a single site T of coordinates (x, y) located at a distance r from the contact. The interferometer Hamiltonian reads

H = H 0 + H tip .
Dot without interaction: When U = 0, this model can be analytically solved [17,19]. Let us summarize the main results (partly given in Ref. [17] and with more details in Ref. [19]) which are necessary for understanding how a pattern of interference rings with a period different from λ F /2 can be seen in the SGM images. The contact being reduced to a single site I coupled to another single site per lead, the lead self-energies [START_REF] Datta | Electronic transport in mesoscopic systems[END_REF] are

only two complex numbers Σ l,r = R l,r (E) + iI l,r (E) = t 2 c < ±1, 0|G R l,r (E)| ± 1, 0 >, G R l,r (
E) being the retarded Green's function of the left and right leads evaluated at the sites directly coupled to I. Using the method of mirror images [START_REF] Molina | [END_REF], G R l,r (E) can be expressed in terms of the Green's function G R 2d (E) of the infinite 2d lattice [START_REF] Economou | Green's functions in quantum physics[END_REF]. Without tip (V = 0), the transmission T σ 0 (E) of an electron of spin σ through the dot making the contact reads:

T σ 0 (E) = 4I r I l (E -4 -V I -R r -R l ) 2 + (I r + I l ) 2 . (4)
If the variation of Σ l,r (E) can be neglected when E varies inside the resonance (typically t c < 0.5 in the continuum limit where the Fermi momentum k F ≪ 1), this is a Lorentzian of width Γ = -(I r + I l ) = -2I and center 4

+ V I + 2R if R l = R r = R.
If one adds a tip potential V = 0 in the right lead, the effect of the tip can be included by adding an amount ∆Σ r (E) = ∆R r (E) + i∆I r (E) to Σ r (E) (see Refs. [17,19,[START_REF] Darancet | [END_REF]). The interferometer transmission T σ (E) is still given by Eq. ( 4), once R r + ∆R r and I r + ∆I r have been substituted for R r and I r . Moreover, when r is sufficiently large, ∆Σ r (E, r) becomes small and one can expand T σ (E) in powers of ∆Σ r (E). The effect of the tip being restricted to a single site T, ∆Σ r can be obtained from Dyson's equation. In the continuum limit and for distances r ≫ k -1 F , one finds [17,19] where ρ and φ are the modulus and the phase of V /(1 -V 0, 0|G R 2d (E)|0, 0 ). Expanding T σ (E) to the leading order ∝ x 2 /r 3 in ∆Σ r , the effect of the tip upon the conductance at a temperature T can be obtained:

∆Σ r t 2 c ρ ≈ - k F x 2 2πr 3 exp[i(2k F r + π/2 + φ)] + O( x 3/2 r 3 ), ( 5
)
∆g = g -g 0 = σ dE(T σ (E) -T σ 0 (E))(- ∂f ∂E ), ( 6
)
where f is the Fermi-Dirac distribution. Assuming [13] -∂f 2 and approximating Σ l,r , ρ, and φ by their values at the Fermi energy E F , one eventually gets [19]:

/∂E ≈ (1/4k B T ) exp -[ √ π(E -E F )/(4k B T )]
∆g(T ) ≈ 2A(T ) cos(2k F r + Φ(T )), (7)
A(T ) = ρx 2 lT √ πkF r 3 exp -[(1 + v 2 )( lT lΓ ) 2 + r lΓ ] (8)
Φ(T ) = φ + v r lΓ -2v( lT lΓ ) 2 . ( 9
) v ≡ (V res I -V I )/Γ gives the energy shift of V I from the resonance V res I ≡ E F -4 -2R in units of Γ. l T = ( √ πk F )/(4k B T
) and l Γ = k F /Γ are two length scales associated respectively to T (Fermi-Dirac statistics) and to Γ (resonant transmission). To obtain Eq. ( 9), we have used asymptotic expansions valid when r > r * ≡ 2l T [1 + l T (1 + |v|)/l Γ ]. The factor 2 in ∆g(T ) comes from the spin degeneracy.

Zeeman splitting without interaction: Let us consider now the case where the spin degeneracy is removed by a Zeeman term ±h due to a parallel magnetic field applied in the dot. This removal is illustrated in Fig. 2 (left), g ↑ 0 (T ) having a peak shifted by an amount v ↑ = -h/Γ (electron with parallel spin) while the shift is v ↓ = h/Γ for the antiparallel spin. Hereafter, we study the effect of the tip at the value V res I (h = 0) for V I indicated by an arrow in Fig 2 left (middle of the valley). When h = 0, ∆g(T ) = σ A(T ) cos(2k F r + Φ σ (T )). The amplitudes A depending on v 2 are independent of σ, but the phases Φ σ depend on the sign of v, and hence of σ. This gives rise to a beating effect, the oscillations of ∆g ↑ (T ) and ∆g ↓ (T ) canceling each other and the tip having no effect on g(T ) when the distance r takes values

r D (n) = 2k F Γ ( l T l Γ ) 2 + (nπ + π 2 ) k F h , ( 10
)
where n is a positive integer. Conversely, the oscillations of ∆g ↑ (T ) and ∆g

↓ (T ) add if r(n) = r D (n) + πk F /(2h).
The SGM image giving ∆g(T ) as a function of the tip position is characterized by a first ring at a distance r D (n = 0) followed by other rings spaced by πk F /h where ∆g(T ) = 0. To optimize the contrast in the images, we calculate for a given value of h the temperature T * and the width Γ * for which ∆g(T , r = r D (n = 0)+πk F /(2h)) is maximum. The extrema are given by two coupled nonlinear algebraic equations which can be solved numerically, yielding T * ≈ 0.73h and Γ * ≈ 0.25h. In Figs. 2 (a) and (b), ∆g(T * , Γ * ) is shown as one varies the tip at the right side of the contact, keeping the tip coordinate y = 0. The figures correspond to two values of h chosen using Eq. ( 10) and the conditions r D (h, n = 0) = 75 and 100. The numerical results give rings at the expected distances r D (n = 0) and r D (n = 1). Though the period of the oscillations is λ F /2, the oscillations around r ≈ r D (n) become so small that one can easily miss a few of them, and draw the conclusion that the period of the oscillations exceeds λ F /2.

Anderson impurity making the contact: Let us now consider the case where the resonance peak of g 0 (E) is not split in two peaks by an applied magnetic field, but by an Hubbard repulsion of strength U acting in the contact. This case is of particular interest, since it describes a quantum dot with an odd number of electrons. This splitting was observed by measuring [4][5][START_REF] Grobis | Handbook of Magnetism and Magnetic Materials[END_REF] the dot conductance as a function of an applied gate voltage V g . The interval between the two conductance peaks is called the Kondo valley. It vanishes below the Kondo temperature, the dot becoming transparent (unitary limit) in the interval of V g where the number of electrons remains odd. Hereafter, we study the interferometer conductance when the dot is biased in the middle of the Kondo valley, describing the resonant level of the dot by the Anderson model. The middle of the Kondo valley corresponds to the symmetric case where

V I + (4 + 2R) = E F -U/2.
When U > Γ, it exhibits three fixed points as the temperature T decreases [24,25]. At large temperatures (T > √ 2U Γ/π), the Anderson impurity coupled to the conduction electrons (left and right 2DEGs) is described by the excitations of the free orbital fixed point. For an intermediate range of temperature ( √ 2U Γ/π > T > T K ), the impurity has a local magnetic moment and the system excitations become different (local moment fixed point). Below the Kondo temperature T K = ( √ 2U Γ/π) exp -[πU/(8Γ)], the local moment is screened by the conduction electrons and the excitations are those of the strong coupling limit. If the occurrence of a magnetic moment can be detected using the Hartree-Fock (HF) approximation [26], more involved many-body methods as the numerical renormalization group (NRG) algorithm [24] or the Bethe ansatz [25] are necessary to describe the Kondo screening of the magnetic moment. Hereafter, we study the effect of the tip for temperatures above T K using the HF approximation. Though this can be questionable in the magnetic region, where a HF-description can give artifacts, it is usually believed that the HF-behaviors are very suggestive, and, when suitably reinterpreted, indicate what we can expect an exact treatment to yield [27,28]. Notably, the HFapproximation does not give the spin-flip processes [29]. This point is extremely important. If the local moment has a finite time of flip τ sf , the conductance oscillations cannot persist beyond a coherence length l φ = k F τ sf , and a mean-field theory breaks down on scales r > l φ . Nevertheless, the oscillations should be given by the HFapproximation on shorter scales. τ sf can be estimated for the Anderson model using NRG [24] or Bethe ansatz [25]. From the expression giving the magnetic susceptibility as a function of the magnetic moment µ, T and T /T K [30][31][32], one can estimate τ sf as a function of T . For the temperature range

√ 2U Γ/π > T ≫ T K which we shall consider, l φ ≫ (l T , l Γ ) (l φ ≈ 1000 when T /T K ≈ 30).
In a mean-field approximation, the impurity potential V I is corrected by real Hartree potentials Σ H σ , given by the self-consistent solutions of coupled equations:

Σ H ↓ = -U π dEf (E, E F , T )ℑG R ↑ (E) (11)
Σ H ↑ = -U π dEf (E, E F , T )ℑG R ↓ (E), ( 12
)
where f is the Fermi distribution at temperature T and Fermi Energy E F and

G σ (E) = [E -V I -4 -Σ l -Σ r - ∆Σ r -Σ H σ ] -1
is the HF Green's function of an electron of energy E and spin σ. Without tip, this gives a splitting of the resonance by an amount U |n I↑ -n I↓ | proportional to the interaction induced magnetic moment. Such a magnetic moment can be modified by the tip. However, the effect of the tip upon Σ H σ remains negligible for the values of r and T used in this study. It was shown previously for a similar model that the tip modifies the HF potentials by an amount ∆Σ HF which decays as 1/r 2 at zero [15]. In the Anderson model, we obtain [START_REF] Kleshchonok | [END_REF] that ∆Σ H σ ∝ (1/r α )(exp -[r/l T ]) at a temperature T . Though the exponent α of the decay can slightly differ when the resonance is narrow (α → 1 when t c → 0), the exponential damping makes the correction quickly negligible. Neglecting the effect of ∆Σ r upon Σ H σ , the dot biased in the middle of a Kondo valley can be described by the theory previously developed for the non interacting dot, with a magnetization which is not now a free parameter induced by an external field, but takes a self-consistent value which depends on T , U and Γ. The corresponding SGM images should also exhibit interfer- ence rings, characterized by radii r D (n) given by Eq. (10) where one has substituted the self-consistent value of U |n I↑ -n I↓ | for 2h. The rings are now spaced by a distance 2πk F /(U |n I↑ -n I↓ |). This makes possible to measure the magnetic moment by electron interferometry. To optimize the contrast of the SGM images, we calculate for a given value of U the temperature T * and the width Γ * for which ∆g(T , r = r D (n = 0) + πk F /(U |n I↑ -n I↓ |)) is maximum. The extrema are now given by two coupled self-consistent differential equations, which can be solved numerically, yielding T * ≈ 0.17U and Γ * ≈ 0.06U . In Fig. 3, the SGM image of a dot biased in the middle of a Kondo valley is shown for a set of optimized values T * and Γ * . One can see two rings in addition to the Fabry-Pérot oscillations. Near the rings, the apparent increase of the spacing between the interference fringes (see lower part of Fig. 3) is reminiscent of the effect reported in Ref. [12] obtained with a QPC instead of a Kondo dot. In Fig. 4, the radii r D (n) of the rings are given as a function of T (4(a)) and U (4(b)) for a set of fixed values of U and Γ or T and Γ respectively. The right color scale gives a visibility parameter [START_REF] Born | Principles of Optics.7th[END_REF] equal to 0 without contrast and to 1 with a perfect contrast. The radii r D (n) obtained from the numerical self-consistent solutions of Eqs. ( 11) and (12) (without neglecting the effect of ∆Σ r in the HF-potentials) turn out to be well approximated by our simplified theory (which neglects it). When T → T K , the radii r D (n) become too large for observing rings. When T → √ 2U Γ/π, the magnetic moment vanishes and r D (n) → ∞. Hence, the observation of the rings requires a fine tuning of the temperature T and of the dot-leads couplings t c . In summary, a set of thermally induced interference rings can be seen when a quantum dot biased around the middle of a Kondo valley is studied above T K with a scanning gate. We have studied a case where the interaction effects are taken into account and which looks realistic enough for being amenable to experimental checks. We believe that the rings can be observed when the contact is biased between two resonances. This belief is supported by the numerical check that a contact made of two coupled sites in series (double dot setup) exhibits similar rings when it is biased between its two resonances. One motivation of this work comes from the interference ring observed [12] using a QPC biased between two channel openings. We believe likely that this is due to similar interference effects, two sharp consecutive channel openings of a QPC playing [17,19] a similar role than the two consecutive resonance peaks considered in this work.

Numerical studies are in progress for confirming this hypothesis. For the Kondo dot, the extension of the study below T K and beyond the mean-field approximation is also in progress.

This research has been supported by the EU Marie Curie network "NanoCTM" (project no.234970). Discussions with B. Brun, K. Ensslin, M. Sanquer and H. Sellier about SGM experiments are gratefully acknowledged.

FIG. 2 :

 2 FIG. 2: (Color online) Left: Conductance g0 of a non interacting dot without tip (V = 0, tc = 0.2, and Γ = 0.003) as a function of VI for T = 0 (blue) T = Γ/2 (green) and T = Γ (red). Cases with (h = 0, solid line) and without spin degeneracy (h = 4Γ, dashed line). The arrow gives the potential VI used in the right figures. Right: ∆g/g0(T * , Γ * , y = 0) as a function of x for a tip potential V = -2. (a): h = 0.0091, T * = 0.0066, Γ * = 0.0023; (b): h = 0.0068, T * = 0.005, Γ * = 0.0017. In all figures, EF = 0.1542 and λF /2 = 8.

FIG. 3 :

 3 FIG. 3: (Color online) ∆g/g0 (upper color scale) as a function of the tip coordinates (x, y) when an Anderson impurity with U = 0.0316 and tc = 0.1592 makes the contact. EF = 0.1542, V = -2 and g0 = 0.2607. The temperature (T * ≈ 0.0054 ≫ TK ) and the resonance width (Γ * ≈ 0.0019) have been calculated for having the largest amplitude ∆g/g0 after r D (n = 0) = 75. The apparent period between the visible maxima of ∆g/g0 is given in the lower figure. It is equal to λF /2 = 8, but becomes larger than λF /2 around r D (n = 0) = 75 and r D (n = 1) = 195.

FIG. 4 :

 4 FIG. 4: (Color online) Radii r D (n) of the rings for n = 0, 1 and 2 as a function of the temperature T (left, U = 0.022) and of the interaction strength U (right, T = 0.0032) with Γ = 0.003. The curves are in the regimes TK ≪ T < √ 2U Γ/π and U > Γ. The dots are obtained numerically, their colors corresponding to a visibility scale [33] indicated at the right. The solid lines are the analytical values of r D (n) derived assuming r D (n) > r * (dashed lines). One can see the ranges of temperature and interaction where the first ring is sufficiently close to the contact for being visible in the SGM images.


Alain Burgisser

Louis Forestier-Coste

François James

Simona Mancini

An expansion-coalescence model for bubble growth in viscous fluids

Keywords: AMS Classification: 82C22, 35Q90, 65C20

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Volcanic eruptions are classified in two families: effusive and explosive. The first ones are characterized by the separate emission of magma and gas. The resulting lava flows or lava domes affect the immediate neighborhood of the volcano. The second ones are more dangerous and are characterized by columns of cinders which may measure several kilometers. When falling on the ground these cinders may accumulate and break human infrastructures as well as pollute air. Finally, both eruptions may lead to pyroclastic avalanches that sweep the volcano flanks at high speeds and cause devastation up to dozens of kilometers from the source. Thus, the understanding of the mechanisms driving volcanic eruptions is useful for both scientific and societal reasons.

It is nowadays clear that the evolution of gas bubbles in very viscous magmas such as rhyolite affect the kind of eruption the volcano may have. This is why we are interested in the modeling of bubbles growth in viscous magmas during an eruption.

While ascending to the surface through volcano conduits, bubbles growth in viscous magmas is governed by decompression and exsolution of volatiles (mainly water) from the magma surrounding the bubbles to the bubbles. This phenomena is usually modeled considering that bubbles are mono-dispersed, i.e. they all evolve in the same way and have the same initial radius (or volume) and mass (or pressure). The most commonly used model consists of a system of two ordinary differential equations describing the time evolution of the radius and mass of a single bubble and coupled with an advection-diffusion equation modeling the space-time evolution of the volatile concentration in the magma surrounding the bubble, see for example [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF], [START_REF] Lensky | Bubble growth during decompression of magma: experimental and theorical investigation[END_REF] and [START_REF] Chouet | Dynamics of diffusive bubble growth and pressure recovery in a bubbly rhyolitic melt embedded in an elastic solid[END_REF]. Two previous works, [START_REF] Forestier-Coste | Numerical resolution of a mono-disperse model of bubble growth in magmas[END_REF] and [START_REF] Lensky | Bubble growth during decompression of magma: experimental and theorical investigation[END_REF], analyzed the dimensionless associated system of equations. This is characterized by two relaxation parameters, namely Θ V and Θ D , which vary by several order of magnitude and account for viscous (Θ D ≪ Θ V ), diffusive (Θ V ≪ Θ D ) or equilibrium (Θ V ∼ Θ D ) growth regimes. In particular, in the diffusion regime, for Θ V ≪ 1, the magma-bubble system is at the mechanical equilibrium; in the viscous regime, for Θ D ≪ 1, the magmabubble system is at the chemical equilibrium; finally, in the equilibrium regime, for Θ V , Θ D ≪ 1, the magma-bubble system is at the thermodynamical equilibrium. In [START_REF] Forestier-Coste | Numerical resolution of a mono-disperse model of bubble growth in magmas[END_REF] we numerically show that when both relaxation parameters converges to zero computational times are very long, and we propose some simplified models of the system to be used for these extremes values of Θ V and Θ D .

Although this mono-disperse bubbles population model gives very good results, as far as porosity is concerned with, when compared to experimental data (see [START_REF] Forestier-Coste | Numerical resolution of a mono-disperse model of bubble growth in magmas[END_REF] and [START_REF] Lyakhovsky | Bubble growth in rhyolitic melts: experimental and numerical investigation[END_REF]), it assumes that bubbles are isolated, not taking into account their interaction, i.e. the coalescence of two (or more) bubbles. Leaving aside coalescence is a severe limitation because of its large impact on the size distribution and because coalescence can create an interconnected network of bubbles from which the gas can escape from magma. Such gas run away is a main control on whether an eruption is effusive or explosive. There have been attempts to overcome the absence of coalescence in models. For instance, in [START_REF] Lovejoy | Bubble distributions and dynamics: The expansion-coalescence equation[END_REF] is presented a kinetic model describing the evolution of the distribution function of a set of bubbles that grow by decompression and coalescence. The distribution function f = f (t, v) represents in [START_REF] Lovejoy | Bubble distributions and dynamics: The expansion-coalescence equation[END_REF] the probability to find at time t a bubble of volume v, and the model is homogeneous in space. The volume growth rate is defined by an exponential function, corresponding to a growth at the thermodynamical equilibrium. The bubbles coalescence is described by a Smoluchowski type operator. Up to our knowledge, this is the only work presenting a statistical description of poly-disperse bubble growth in magma. Nevertheless, this model cannot describe the evolution of a bubble population when considering non-equilibrium growth, as for example growth of bubbles at the chemical equilibrium and not at the mechanical one, Θ V ≫ Θ D → 0. Therefore our aim is to propose a more general kinetic model describing the evolution of bubbles populations also for non-equilibrium growths. The abstract writing of the kinetic model is rather classic by means of Liouville theorem. The main goal is to define coalescence operator and rate depending on the volume and mass variables and not on position and velocity, as usual, and to validate it on a set of experimental data. Our coalescence operator turns out to be two-dimensional and is in some sense a direct generalization of the Lovejoy model [START_REF] Lovejoy | Bubble distributions and dynamics: The expansion-coalescence equation[END_REF]. Due to their complexity, two-dimensional coalescence operators have been rarely studied and simulated in the past, see for example [START_REF] Qamar | Solving population balance equations for twocomponent aggregation by a finite-volume scheme[END_REF] and [START_REF] Kumar | Numerical methods for solving two-dimensional aggregation population balance equations[END_REF]. In [START_REF] Forestier-Coste | A finite volume preserving scheme on nonuniform meshes and for multidimensional coalescence[END_REF], a numerical scheme is developed and validated by several numerical tests for the two-dimensional continuous Smoluchovski equation. We point out that in our modeling we have to overpass one main difficulty: to simplify the system of equations of the mono-disperse model in order to avoid the coupling with the advection-diffusion equation, so to diminish computational costs. This paper is organized as follows. In section 2 we present the abstract kinetic model and recall the main properties of its solution. Volume and mass rates are described in section 3, while section 4 is devoted to coalescence kernels. Next, in section 5 we describe the numerical discretization of the model. Finally, numerical results, together with their comparison with experimental data are gathered in sec-tion 6. In section 7 we discuss the proposed model and results and present some research directions to be further investigated.

The kinetic model

We want to describe at a statistical level the evolution of gas bubbles by decompression, exsolution and coalescence in a very viscous fluid. In this section we introduce the kinetic equation we will use to model this phenomena and recall the main properties of its solution.

We note that our kinetic model is defined on non-standard independent variables, since there is no dependence on bubbles velocity and position. In fact, we assume that bubbles are entrained by the magma which raises in the conduit with a laminar and constant velocity u, see also [START_REF] Lovejoy | Bubble distributions and dynamics: The expansion-coalescence equation[END_REF]. This is a first approximation since of course the magma fluid satisfies the Navier-Stokes equations for an incompressible fluid. By this assumption it follows that bubble positions are only related here to their depth in the conduit (and not to their distance from the conduit borders), which in his turn is linked to the magma pressure. Hence bubble positions are implicitly defined by the value of the magma pressure. This is consistent with experimental data which are obtained considering small samples of magma, where bubbles are homogeneously distributed in space and where their displacements can be neglected.

Therefore the distribution function f = f (t, v, m) represents the probability to find at time t ≥ 0 a bubble of volume v ∈ [v min , ∞] and mass m ∈ [m min , ∞]. In the following, we will then consider the domains

Ω T = [0, T ] × Ω, and Ω = [v min , ∞] × [m min , ∞],
describing the time interval and the phase space in which the bubbles evolution takes place, and, for brevity, we will denote by ω the phase space variables (v, m), we shall call ω the size of a bubble, implicitly meaning the bubble volume and mass; dω will denote the infinitesimal element dv dm. This choice of variable is consistent with the mono-disperse model described in [START_REF] Lensky | Bubble growth during decompression of magma: experimental and theorical investigation[END_REF].

Since bubbles grow by expansion and coalescence, the distribution function f (t, v, m) must satisfy the following kind of kinetic equation:

∂ t f + ∂ v (V f ) + ∂ m (M f ) = Q(f ).
(

) 1
The left-hand side in (1) is a conservative equation and describes the evolution by decompression and exsolution of the set of bubbles, where the two growth rates for the volumes and masses are defined by the non-linear functions V = V (t, v, m) and M = M (t, v, m), respectively:

V (t, v, m) = 3 Θ V (m -V ), (2)
M (t, v, m) = 1 Θ D (M -m), (3)
with V = V (t, v) and M = M (t, v) two equilibrium manifolds. Their expression will be given when introducing the volcanology model, see section 3. We note here that this kind of rates implies incoming drift on the boundaries.

The right-hand side in (1) describes bubbles coalescence by means of the Smoluchowski continuous coalescence operator:

Q(f )(t, v, m) = 1 2 ω 0 H(t, ω ′ , ω -ω ′ ) f (t, ω ′ ) f (t, ω -ω ′ ) dω ′ -f (t, ω) ∞ 0 H(t, ω ′ , ω) f (t, ω ′ ) dω ′ , (4)
where the integrals have to be understood in the two-dimensional setting (v, m).

Moreover, H(t, ω, ω ′ ) is the coalescence kernel which must be positive and symmetric, H(t, ω, ω ′ ) > 0 , H(t, ω, ω ′ ) = H(t, ω ′ , ω).

Its expression will be given later on, see section 4. We recall here that the first term on the right-hand-side of (4) is the gain term and counts the bubbles created by the coalescence of a bubble of size ω ′ and one of size ωω ′ . The second term on the right-hand-side of (4) is the loss term and counts the bubbles of size ω disappearing by coalescence with all other bubbles.

To conclude the description of the coalescence phenomena, we must describe how the sizes of two coalescing gas bubbles define the size of the created single bubble. Let us denote by (v ′ , m ′ ) and (v ′′ , m ′′ ) the sizes of the two coalescing bubbles and by (v, m) the size of the created one. It is physically natural to assume that masses sum up, but this is not the case for volumes. We recall that we have a poly-disperse population of bubbles and that in our model volume and mass are independent variables. This is rather different with respect to the modeling usually done when considering polymers, in which mass and volume are equivalent, see for example [START_REF] Escobedo | On a kinetic equation for coalescing particles[END_REF] or [START_REF] Carrillo | Exponential decay towards equilibrium for the inhomogeneous Aizenman-Bak model[END_REF] and references therein. Nevertheless, we have no way to know how volumes combine for two gas coalescing bubbles and therefore, in a first approach, we shall assume here that they do actually sum up. In conclusion volumes and masses verify:

v = v ′ + v ′′ , m = m ′ + m ′′ ,
and, as a consequence of the perfect gas law holding in each bubble, the pressure p in the created bubble may be given as a convex combination of the pressures in the two coalescing bubbles p ′ and p ′′ :

p = δp ′ + (1 -δ)p ′′ , δ = v ′ v ′ + v ′′ .
We associate to (1) an initial condition f 0 = f 0 (v, m) ≥ 0 describing the distribution of bubbles at the initial time 0 in the phase-space Ω. The initial condition f 0 given by the experimental measures, is defined by a sum of Dirac masses. Concerning the boundary conditions, we assume that there do not exist bubbles with volume or mass equal to zero :

f (t, 0, m) = f (t, v, 0) = 0,
and that for large values of v or m the distribution function f decays to zero :

lim v→∞ f (t, v, m) = lim m→∞ f (t, v, m) = 0.
The numerical phase space domain will be bounded both in v an m, and we will then consider no flux conditions for the larger values of v and m. On the other hand for v = 0 or m = 0 we note that the respective drifts are coming into the domain (see [START_REF] Carrillo | Exponential decay towards equilibrium for the inhomogeneous Aizenman-Bak model[END_REF] and (3)), so that if the initial condition f 0 = 0 for v = 0 or m = 0 are set (which seems natural), then f = 0 at each time.

Let us recall now some mathematical properties of the expansion-coalescence equation [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF]. Let r and s be two integers, then the moments M r,s = M r,s (t) of order p = r + s of the distribution function f are given by:

M r,s = ∞ 0 v r m s f (t, ω) dω . ( 5
)
For instance, the zero order moment:

M 0,0 = ∞ 0 f (t, v, m) dv dm , (6)
represents the number of bubbles N (t) = M 0,0 ; the first order moment (for r = 1 and s = 0):

M 1,0 = ∞ 0 v n(t, v, m) dv dm , (7)
represents the volume occupied by the bubbles, V(t) = M 1,0 ; and the first order moment (for r = 0 and s = 1):

M 0,1 = ∞ 0 m n(t, v, m) dv dm , (8)
represents the mass of the bubbles, M(t) = M 0,1 .

Since the left-hand side of equation ( 1) is a conservation equation, when neglecting coalescence, the number of bubbles must be conserved, that is N (t) = N (0). This can be proved integrating over the interval [0, t]. On the other hand, when neglecting the expansion terms in (1), the first order moments, M 1,0 and M 0,1 , must be conserved at least until the gelation time, if any, is reached, see [START_REF] Filbet | Numerical simulation of the Smoluchowski equation[END_REF], [START_REF] Forestier-Coste | A finite volume preserving scheme on nonuniform meshes and for multidimensional coalescence[END_REF] and reference therein. This means that the bubbles volume and the bubbles mass must be conserved : V(t) = V(0) and M(t) = M(0). We finally remark that if the first order moments are constant in time, then any linear combination of them will be constant in time, too. We note here that at a discrete level it will not be possible to conserve both first order moments and we will then choose to conserve only the bubbles mass, see [START_REF] Forestier-Coste | A finite volume preserving scheme on nonuniform meshes and for multidimensional coalescence[END_REF].

Growth rates

In order to define the volume and mass growth rates, we have to come back to the microscopic, or mono-disperse, model. We recall here the coupled system of equations governing the growth of a single bubble as it is presented in its dimensionless form in [START_REF] Lensky | Bubble growth during decompression of magma: experimental and theorical investigation[END_REF]. The function F (t, 1) is the volatile flux incoming the bubble from the influence region, ρ m is the magma density, R, S, M and C(t, r) respectively are the radii of the bubble and of the influence region, the bubble mass and the volatile concentration in the influence region.

Assuming that a bubble is surrounded by an influence domain (that is the region of magma exchanging gas with the bubble) represented by a sphere of radius S, see figure 1, the radius R and mass M of a volatile bubble growing in a viscous magma by decompression and exsolution are governed by the following system of two ordinary differential equations:

Ṙ = R η ef f Θ V P -P a (t) - Σ R , (9)
Ṁ = 3ρ m Θ D F (t, R), (10)
where η ef f is the magma viscosity, ρ m the magma density, Σ the surface tension between the gas and the magma, P a (t) the ambient (magma) pressure which will be detailed later on, and Θ V and Θ D are two relaxation parameters accounting for viscosity and diffusion, see section 6 for more details.

Equation [START_REF] Forestier-Coste | Numerical resolution of a mono-disperse model of bubble growth in magmas[END_REF] states that the time evolution of a bubble radius is governed by the difference between the bubble pressure and the magma pressure augmented by the surface tension.

Equation [START_REF] Kumar | Numerical methods for solving two-dimensional aggregation population balance equations[END_REF] is the differential version of the conservation of mass, which states that the variation of the bubble mass is equal to the mass of volatile diffusing between the magma and the bubble. In [START_REF] Kumar | Numerical methods for solving two-dimensional aggregation population balance equations[END_REF] the flux F (t, R) is defined by means of the volatile concentration in the magma surrounding the bubble, C = C(r, t), satisfying the following advectiondiffusion equation, for t > 0 and r ∈ [R, S]:

∂ t C + R 2 Ṙ r 2 ∂ r C = 1 Θ D 1 r 2 ∂ r r 2 D∂ r C , (11)
with D the diffusion coefficient and the following boundary conditions:

C(t, R) = C H √ P , (r 2 D∂ r C) r=S = 0, (12)
In equation ( 11), the left-hand side is an advective term that takes into account the fact that the border of the bubble moves with a radial velocity given by ṘR 2 /r 2 , thanks to the incompressibility of the magma. The right-hand side describes the radial diffusion of the volatile.

The left boundary condition in [START_REF] Lovejoy | Bubble distributions and dynamics: The expansion-coalescence equation[END_REF] asserts that the bubble boundary and the magma are at a chemical equilibrium, which is described by the solubility condition, with C H = K H √ P i , P i the initial pressure and K H the Henry solubility constant. The right boundary condition in [START_REF] Lovejoy | Bubble distributions and dynamics: The expansion-coalescence equation[END_REF], states that at the outer border (r = S) of the influence region the concentration flux is null, i.e. there is no volatile loss or gain. Finally, the initial concentration is assumed to be homogeneous in space C(0, r) = C i , which is equal to C H if the initial state is at the chemical equilibrium.

The incoming volatile flux in the bubble F (t, R) is then defined by:

F (t, R) = r 2 D∂ r C r=R . (13)
Following [START_REF] Lensky | Bubble growth during decompression of magma: experimental and theorical investigation[END_REF], assuming that the magma raises at a constant velocity, implies that magma pressure is linearly decreasing. The magma pressure eventually reaches a final value, which is the pressure at the surface of the volcanic conduit normalized to 0. The fact that pressure cannot become negative implies that time t belongs to the bounded interval [0, 1]. Hence, we assume that the ambient pressure P a (t) satisfies in dimensionless form:

P a (t) = 1 -t. (14)
Therefore, we may consider the magma pressure as a dimensionless time, recovering information on the depth of the bubbles in the conduit.

Following [START_REF] Lensky | Bubble growth during decompression of magma: experimental and theorical investigation[END_REF], the volume of the influence region S 3 (t) is constant in time so that it satisfies:

S 3 = R 3 + S 3 0 , (15)
with S 3 0 the dimensionless volume occupied by the magma surrounding a bubble of radius zero, see figure 1.

We note that (9) can be written in the form

Ṙ = R Θ V (P -P ), (16)
where the equilibrium manifold P = P (t, R) is given by:

P (t, R) = P a (t) + Σ R , (17)
which is the mechanical equilibrium Θ V → 0.

In order to diminish computational costs we want to avoid the coupling with the advection-diffusion equation ( 11) in [START_REF] Kumar | Numerical methods for solving two-dimensional aggregation population balance equations[END_REF]. By analogy with the radius growth, we introduce the equilibrium manifold defined by the chemical equilibrium case, that is when Θ D → 0:

M = M (R) = -Y C H + Y 2 C 2 H + 4R 3 Y C 0 2 4R 3 , (18)
with Y = ρ m S 3 0 , C 0 the total volatile concentration when the bubble radius is zero, and A a scaling function. We replace [START_REF] Kumar | Numerical methods for solving two-dimensional aggregation population balance equations[END_REF] by:

Ṁ = A Θ D (M -m). (19)
Numerical comparisons show that the choice A = 1 is relevant. The accuracy, for different values of the relaxation parameters, of this simplification is shown in section 6.

We finally note that the macroscopic quantities measured in experiments are the bubble radius R (or their volume V ) and the porosity α = α(t), that is the fraction of volume occupied by the volatile with respect to the total volume. In the case of a single bubble surrounded by its influence region, the porosity may be computed by the following formula, see [START_REF] Lensky | Bubble growth during decompression of magma: experimental and theorical investigation[END_REF]:

α(t) = R S 3 . ( 20
)
When considering a population of bubbles, under the assumption that the volume of magma is constant, see [START_REF] Qamar | Numerical solution of population balance equations for nucleation growth and aggregation processes[END_REF], we can compute the porosity by:

α(t) = V V + W 0 , (21)
where V = V(t) is the volume of bubbles defined by V = M 1,0 , see [START_REF] Forestier-Coste | Croissance et coalescence de bulles dans les magmas : analyse mathématique et simulation numérique[END_REF], and W 0 is the magma volume. Knowing the initial porosity α 0 and the initial distribution f 0 = f (0, v, m), we can determine W 0 using (21):

W 0 = 1 -α 0 α 0 V 0 . ( 22
)
We note that, since the first order moment is conserved by the coalescence operator, the only variation for the porosity will be given by the left-hand side of equation ( 1), and so the numerical results for porosity obtained by means of the mono-disperse model and the one obtained by means of the kinetic equation must be identical, unless gelation occurs.

Coalescence kernel

We now consider the coalescence of two bubbles of sizes ω ′ = (v ′ , m ′ ) and ω ′′ = (v ′′ , m ′′ ) that yields to a single bubble of size ω = (v, m). We assume that the coalescence is instantaneous and that the relaxation time of the created bubble to get back to a spherical geometry is zero. Hence, we neglect the possibility of having a chain of coalescing bubbles, although this is an important and observed behavior, see [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF]. Bubble chains influence the way volatiles may escape from magma and is thus linked to the type of eruption the volcano has. Multiple coalescence and bubble chains may be studied in future works by considering for instance a new characteristic for the bubble population linked to the geometrical form of the created bubble.

We recall that bubbles have no relative velocity, hence the only way they may coalesce is by expanding until the film of magma in between two bubbles is too thin to hold and breaks up. We want to define a physical relevant coalescence kernel H = H(t, ω, ω ′ ) which is positive and symmetric. The coalescence kernel represents the rate of coalescence and it has the dimension of the inverse of a time. Under the assumption that the mass of the magma film contained in between the walls of two coalescing bubbles is conserved during the coalescence process, in [START_REF] Castro | Mechanisms and dynamics of bubble coalescence in silicic magmas[END_REF] are described and compared two coalescing times : the planar profile one, T p , and the stretch one, T s . The coalescence kernel is then defined by:

H p (t, ω, ω ′ ) = 1/T p , H s (t, ω, ω ′ ) = 1/T s .
In figure 2 we show the two configurations of coalescing bubbles, on the left the planar coalescence and on the right the stretch one, in gray the region in which magma mass is conserved. We detail here the computations leading to the definition of the planar coalescence kernel H p (t, ω, ω ′ ). For the discussion in the case of the stretching kernel (or time) we refer the reader to [START_REF] Castro | Mechanisms and dynamics of bubble coalescence in silicic magmas[END_REF]. The magma film mass contained in the flat plane between two merging bubbles is given by:

m m (t) = ρ m πz(t)a 2 (t),
where z(t) is the distance between the bubbles and a(t) corresponds to the radius of the disc of the flat section of the bubble (which is assumed to be the same for each bubble), see Figure 2. Assuming that m m (t) is constant in time it is possible to define a rate of coalescence H p (t, ω ′ , ω ′′ ) of two merging bubbles, which border approach with a speed c, as the inverse of the time needed to reach the breaking film distance, z f , see [START_REF] Castro | Mechanisms and dynamics of bubble coalescence in silicic magmas[END_REF].

H p (t, ω, ω ′ ) = c a 0 √ z f √ z - √ z f .
For instance, we assume that a(t) linearly grows in time, a(t) = a 0 + ct, with a 0 the initial radius given by

a 0 = ε((v ′ ) 1/3 + (v ′′ ) 1/3 )/2, ( 23
)
with 0 < ε ≪ 1, and c the rate at which the film thin-off, which is given here by the mean between the two bubbles radius growth rates,

c = ( Ṙ′ + Ṙ′′ )/2.
The factor ε and the final distance z f (before break out) are determined by the physical characteristics of the magma and by experimental measures. We finally obtain the following coalescence kernel:

H p (t, ω ′ , ω ′′ ) = ( Ṙ′ + Ṙ′′ ) + ε((v ′ ) 1/3 + (v ′′ ) 1/3 ) √ z f ( √ z - √ z f ) + , (24)
where the distance between the two bubbles z = z(t) may be expressed in terms of the volume of each bubble and of the porosity α of the magma:

z = (v ′ ) 1/3 + (v ′′ ) 1/3 1 -α 1/3 α 1/3 . ( 25
)
This kernel is positive and symmetric. Considering the homogeneous coalescence equation ∂ t f = Q p (f ), with the kernel defined by ( 24) and (25), we have the following: Proposition 4.1 For any function g(v, m), the functions

f s = f s (v, m) given by f s (m, v) = 0 if m > Σv 2/3 , g(m, v) otherwise, ( 26
)
are equilibrium states of the coalescence operator Q p (f s ).

Proof: We first note that if f = f (t, ω) is a non-trivial solution, that is if f = 0 and Q p (f ) = 0 then we also have

∞ 0 Q p (f ) dω ′′ = 0, which is equivalent to - 1 2 ∞ 0 ∞ 0 H p (t, ω ′ , ω ′′ ) f (t, ω ′ ) f (t, ω ′′ ) dω ′ dω ′′ = 0. ( 27
)
We seek for non-trivial solutions, thus (27) holds if and only if H p (t, ω ′ , ω ′′ ) = 0, that is if and only if ( Ṙ′ + Ṙ′′ ) + = 0. Let us define the function k(t) for t ∈ [0, 1] as

k(t) = Θ V ( Ṙ′ + Ṙ′′ ) = m ′ (v ′ ) 2/3 -(1 -t)(v ′ ) 1/3 -Σ + m ′′ (v ′′ ) 2/3 -(1 -t)(v ′′ ) 1/3 -Σ = t((v ′ ) 1/3 + (v ′′ ) 1/3 ) -((v ′ ) 1/3 + (v ′′ ) 1/3 ) + m ′ (v ′ ) 2/3 + m ′′ (v ′′ ) 2/3 -2Σ. Note that H p (t, ω ′ , ω ′′ ) = 0 if and only if k ≤ 0. Since k(t) is increasing in time, if k(1) ≤ 0 then k(t) ≤ 0 and H p (t, ω ′ , ω ′′ ) = 0 for all t ∈ [0, 1]. It is easily seen that k(1) is negative if m ′ (v ′ ) 2/3 + m ′′ (v ′′ ) 2/3 ≤ 2Σ. ( 28
)
Assume that (m ′ , v ′ ) and (m ′′ , v ′′ ) verify (28) so that Ṙ′ + Ṙ′′ < 0. We have to consider two possibilities.

i) The first one is m ′ > Σ(v ′ ) 2/3 and m ′′ < Σ(v ′′ ) 2/3 , then the growth rates are such that Ṙ′ > 0 and Ṙ′′ < 0. However, H p (1, ω ′ , ω ′ ) > 0, and f (t, ω ′ ) must be the zero function for this choice. In other words if one of the two growth rate is positive, the stationary solution must be the trivial one.

ii) The second possibility is m ′ < Σ(v ′ ) 2/3 and m ′′ < Σ(v ′′ ) 2/3 , Then Ṙ′ < 0 and Ṙ′′ < 0, so that H p (t, ω ′ , ω ′′ ) = 0. Moreover, H p (t, ω, ω) = 0, so that f (t, ω) may be non-zero for all time t. In other words, there is a non-trivial stationary solution. This can be resumed by defining the steady state f s (ω) by ( 26).

On the contrary, it is clear that ∂ t f s = 0 and replacing f s in (27), we also obtain that Q p (f s ) = 0.

From a physical point of view, Proposition 4.1 means that if both bubbles have negative radius growth rates, that is if m < Σv 2/3 for both bubbles, then they cannot coalesce.

Concerning the stretching coalescence kernel, following [START_REF] Castro | Mechanisms and dynamics of bubble coalescence in silicic magmas[END_REF], we can define H s (t, ω, ω ′ ) as, see [START_REF] Castro | Mechanisms and dynamics of bubble coalescence in silicic magmas[END_REF] for more details:

H s (t, ω, ω ′ ) = z 2a 0 |P ′ -P ′′ | ,
where a 0 and z are defined as for the planar kernel, (23) and (25), and |P ′ -P ′′ | is the difference of pressure between the two merging bubbles and can be computed knowing the volume and mass of each bubble, by means of the perfect gas law P = m/v:

H s (t, ω, ω ′ ) = v ′ v ′′ ǫ|v ′′ m ′ -v ′ m ′′ | 1 -α 1/3 α 1/3 (29)
From the expression (29) we see that, on the contrary of what happens for Q p , there is no equilibrium state for Q s , and that the coalescence kernel H s can be defined only for bubbles having different gas pressures.

Numerical resolution

In the previous sections we have introduced all the ingredients needed to build the kinetic model. We resume here the complete model and the first step is to write the preceding equations in terms of mass and volume. The expansion-coalescence of a population of bubbles is governed by the kinetic equation ( 1):

∂ t f + ∂ v (V f ) + ∂ m (M f ) = Q(f ),
where the density function f = f (t, v, m) represents the distribution of bubbles at time t in terms of their volumes v and masses m. We underline once more that this model is based on the non-standard variables volume v and mass m, and has no kinetic variable as position x or velocity v. The respective expansion rates V and M are defined by ( 2) and (3):

V (t, v, m) = 3 Θ V m -(1 -t)v -Σv 2/3 , (30)
M (t, v, m) = 1 Θ D    -C H ρ m w + C 2 H ρ 2 m w 2 + 4ρ m vwC 0 2 4v -m    , (31)
where w is the volume of the influence region defined in analogy to (22):

w = 1 -α α v , (32)
and where the coalescence operator takes the form of the Smoluchowski continuous coalescence operator (4), with the kernel given by expression (24) or (29).

By defining the influence region volume w by (32), we have reduced the number of variables for the distribution function. In fact, by considering relation [START_REF] Qamar | Numerical solution of population balance equations for nucleation growth and aggregation processes[END_REF] we should take into account f = f (t, v, m, w) with a zero growth rate in terms of the variable w. From a physical point of view, definition (32) says that bubbles occupying the same volume v have the same influence region w. This assumption may be relaxed, but this has computational costs, since the numerical resolution would involve three variables.

We turn now to the discretization of the model. We split the numerical resolution of (1) in two parts. At each time step we first solve explicitly in time the coalescence term, right-hand-side of (1), and then implicitly in time the expansion term, lefthand-side of (1). We note that this is equivalent to solve the whole equation (1) in a single time step. In fact, for n ∈ N, let us denote the discrete times by t n , the time steps by ∆t n = t n+1t n and the time discretization of the distribution function by f n = f (t n , v, m). Then the explicit discretization of the coalescence step reads:

f n+1 = f n + ∆t n ζ(t n )Q(f n ),
where f n+1 denotes the value of the numerical solution after the coalescence step and has to be used to compute the expansion term. Then, the implicit discretization of the expansion term is given by:

f n+1 + ∆t n ∂ v (V n+1 f n+1 ) + ∂ m (M n+1 f n+1 ) = f n+1 .
Combining the two discretizations we obtain:

f n+1 -f n ∆t n + ∂ v (V n+1 f n+1 ) + ∂ m (M n+1 f n+1 ) = ∆t n ζ(t n )Q(f n ),
which is a semi-implicit time discretization of (1).

We will not detail here the discretization of the coalescence term. This has been the subject of [START_REF] Forestier-Coste | A finite volume preserving scheme on nonuniform meshes and for multidimensional coalescence[END_REF], where a multi-dimensional numerical scheme is studied, defined on non-uniform meshes and such that the first order moment is conserved. We note that, in the multi-dimensional framework, it must be chosen which first moment to conserve, since numerically it is not possible to conserve simultaneously all the first moments. We choose to conserve the moment representing the total mass of water M(t) = M 0,1 (t). We finally recall that this model is among the firsts considering two-dimensional coalescence and that numerical schemes for this kind of problem have not been widely studied yet, see [START_REF] Qamar | Numerical solution of population balance equations for nucleation growth and aggregation processes[END_REF], [START_REF] Qamar | Solving population balance equations for twocomponent aggregation by a finite-volume scheme[END_REF], [START_REF] Qamar | On the solution of population balances for nucleation, growth, aggregation and breakage processes[END_REF].

The scheme solving the expansion term is implicit in time because we want to avoid a stability (CFL) condition for this time discretization. Indeed, solving explicitly the expansion term would imply a CFL condition on the time step ∆t n , which is bounded by the minimum of the relaxation parameters Θ V and Θ D . This causes very long computational times, in particular when considering the chemical and/or mechanical equilibrium (Θ V ≪ 1 and/or Θ D ≪ 1). By using an implicit scheme, the time step will be fixed by the positivity condition on the coalescence term which is only bounded by the Θ V parameter, see [START_REF] Forestier-Coste | A finite volume preserving scheme on nonuniform meshes and for multidimensional coalescence[END_REF].

The experimental measurements are usually defined on a uniform mesh in terms of the bubbles radii, conditioning the discretization in v and m. This distribution of initial radii will be chosen as the basis to write the grid mesh in volume and mass, assuming that measurements are done at the thermodynamical equilibrium. Hence, fixing N grid points, we have for i = 1 . . . N , the following definitions for the mesh points in v and m:

v i-1/2 = R 0 + i∆R R i 3 , m i-1/2 = v i-1/2 + Σv 2/3 i-1/2 ,
where R i is the mean radius of the initial data (used also for the scaling), R 0 is the minimum radius considered in experiments and ∆R is the step of the uniform grid defining the initial data. We also define the space mesh sizes:

∆v i = v i+1/2 -v i-1/2 , ∆m i = m i+1/2 -m i-1/2 ,
and the middle grid points:

v i = (v i+1/2 + v i-1/2 )/2 , m i = (m i+1/2 + m i-1/2 )/2 .
For n ∈ N and i, j = 1 . . . N , we can now define the discrete distribution function f n ij which is an approximation of f (t n , v i , m j ), and the discrete marginal with respect to the volume:

V(v i ) = N j=1 f n ij ∆m j .
We split the numerical resolution of the expansion part in two steps: first we consider the discretization implicit in time of the term ∂ v (V f ), then we solve the part concerning ∂ m (M f ). This is equivalent to a directional splitting, solving first equation:

∂ t f + ∂ v (V f ) = 0, (33)
and then equation:

∂ t f + ∂ m (M f ) = 0. ( 34
)
We report here the implicit finite difference scheme for (33), equation (34) being solved analogously. Let j be fixed and i = 1 . . . N , and define:

V n ij = 3 Θ V m j -v i (1 -t n ) -Σv 2/3 i
and its positive and negative parts:

V + n ij = max(0, V n ij ) , V -n ij = min(0, -V n ij ) .
Then the equation ( 33) is discretized implicitly in time by:

1 + ∆t ∆v i |v n+1 ij | f n+1 ij - ∆t ∆v i (V + f ) n+1 i-1,j + (V -f ) n+1 i+1,j = f n ij (35)
To solve this implicit bi-dimensional scheme we transform it by writing the unknown square matrix f n ij of size N × N , in a vector X n h of size N 2 : for each time step n and index j = 1 . . . N we define h = i + (j -1)N for i = 0 . . . N and X n h = f n i,j , that is:

X n = (f n 1,1 , f n 2,1 , f n 3,1 , . . . , f n N -2,N , f n N -1,N , f n N,N ) T .
Moreover, we define the vectors a n , c n of size N 2 -1 by, for j = 1 . . . N :

a n h = -∆t ∆v i (V -) n i+1,j , i = 1 . . . N -1 0 , i = N c n h = -∆t ∆v i (V + ) n i-1,j , i = 2 . . . N 0 , i = 1
and the vector b n of size N 2 by, for i, j = 1 . . . N :

b n h = 1 + ∆t ∆v i |V n ij | ,
and the tridiagonal matrix D n which the h th line reads:

D n h = (0 . . . 0 c n h b n h a n h 0 . . . 0).
The discrete equation (35) then reads:

D n+1 h,h-1 X n+1 h-1 + D n+1 h,h X n+1 h + D n+1 h,h+1 X n+1 h+1 = X n h , (36)
which is solved by a LU decomposition. We can solve in the same way the equation (34). The algorithm involves, for the expansion term, at each time step n the resolution of (33) and then (34), before returning the resolution of the coalescence part.

Numerical and experimental results

We briefly explain how experiments are done because our numerical simulations are compared with experimental data, see [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF]. Various samples of the same magma are placed with water in sealed capsules and heated under high pressure for several days in order to ensure water saturation. A sudden decompression step at constant temperature, nucleates bubbles, that is creates microscopic gas bubbles (the mean radius being of order of one micron). Everything is left at rest until the thermodynamical equilibrium is reached. Then several decompression steps are applied rhythmically until a certain ambient pressure is reached, which approximates a constant decompression rate. After the final step, samples are very rapidly cooled, in order to quench bubbles in the magma. Thanks to microscopic observations, the number of bubbles is evaluated and theirs radii are measured so as to obtain a table with bubbles separated in classes of radius. An experimental series is defined by several experiments quenched at progressively lower final pressures under identical starting conditions and decompression rate. The nucleation step, however, cannot be reproduced exactly from an experiment to another. As a result, samples have slightly different initial bubble size distributions. These differences, although visible on the tabulated data, are smaller than the changes induced by coalescence and growth. Experimental conditions fix the values of T , α 0 , η ef f , Σ, ρ m , D, ∆P , P i , C i . The scaling value R i is defined as the mean bubble radius over all the measured radii at initial time.

We first verify our expansion model by neglecting coalescence. We recall that the ambient pressure acts in this model as the inverse of time, so that when time varies from 0 to 1, ambient pressure varies from 1 to 0. We consider the following initial data, before scaling: f 0 ij ∆v i ∆m j = δ 20,20 with δ I,J a Dirac mass such that δ I,J = 1 if i = 20 and j = 20, zero otherwise. This function corresponds to a mono-disperse initial data, that is we consider a population of bubbles with all the same radius and mass. The physical initial conditions are the following: the porosity is α 0 = 7.79%, the diffusion coefficient is D = 10 -10 m 2 s -1 , the volatile concentration is C i = 4.21, the ambient pressure P a (0) = 10 8 P a, the surface tension is σ = 0.1 Jm -2 , the magma density is ρ m = 2400 kg m -3 , the temperature is T = 825 o C, the decompression rate is ∆P = 10 5 P a s -1 , the mean radius is R m = 10 -5 m, the uniform radius step is ∆R = 0.5 • 10 -6 m and we choose N = 250 discretization points.

The two relaxation parameters, Θ V and Θ D , in equations ( 31)-(30) are related to the experimental data by:

Θ V = 4η ef f ∆P P 2 i , Θ D = ∆P R 2 i D P i . ( 37
)
where the quantities indexed by i are computed from the initial values.

We note here that we assume that the temperature T , as well as the diffusion coefficient D, and the viscosity one, η ef f , are constant. This is of course not always the case as also shown in [START_REF] Lensky | Bubble growth during decompression of magma: experimental and theorical investigation[END_REF] and considered in [START_REF] Forestier-Coste | Croissance et coalescence de bulles dans les magmas : analyse mathématique et simulation numérique[END_REF], and generalization to variable temperature, viscosity and diffusion is possible, but goes beyond the aim of this article.

The validation is carried out by comparing the outputs of three different computations. The two firsts (edp and edo simpl) refer to the simplified mono-disperse model and the kinetic model. Both are based on the approximation of the flux term in the mass rate (19). A perfect match between these two models indicates the implementation of the kinetic model is correct. The third computation (edo gen) refers to the original mono-disperse model, with the coupling with the advectiondiffusion equation [START_REF] Lensky | Bubble growth during decompression of magma: experimental and theorical investigation[END_REF]. It is used as a reference so as to assess the impact of the simplified flux.

In figure 3 we show the behavior of the solutions for four couples of the relaxation parameters, (Θ V , Θ D ) = (10 -3 , 10), (Θ V , Θ D ) = (10 -5 , 10), (Θ V , Θ D ) = (10 -5 , 10 -3 ) and (Θ V , Θ D ) = (10 -1 , 10 -1 ). The first two couples (top) correspond to two physical situations far from the chemical equilibrium, Θ D ≫ 0, one of them being at the mechanical equilibrium (top-right), the smallest Θ V , and the other one not (top-left). In both cases, the difference between the simplified mono-disperse model and the kinetic equation is negligible. The third model deviates from the two approximate models at low pressures. These results suggest that the kinetic model implementation is correct, but that the simplified flux introduces discrepancies. Considering that the values of the relaxation parameter Θ D represent a worse-case scenario with respect to the modeling situations to which we apply our model, it appears that the flux approximation is acceptable. The other two couples (bottom) correspond to the thermodynamical equilibrium with the smallest relaxation parameters values (bottom-left), and to a case neither at the mechanical, or at the chemical equilibrium (bottom-right). We note the very good agreement for the equilibrium case, which includes a radius jump at the beginning of the computation due to the initial data not at the thermodynamical equilibrium. The kinetic solution of the non-equilibrium case seems to better approximate the coupled system than the simplified system solution. We now compare the numerical solution of the kinetic model with the experimental data given in [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF]. The initial condition corresponding to the data ABG1, reads in the kinetic framework: with δ I,J Dirac masses such that δ I,J = 1 if i = I and j = J, zero otherwise.

f 0 ij ∆v i ∆m j = 1 127 34 k=18 D k δ k,k (38)
We compare in figure 4 the size distributions computed by our code (curves) and the data P P E4, P P E1, P P E10, and P P E11 (areas) of [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF]. The physical values are the following: the decompression rate is ∆P = 5 • 10 5 P a s -1 , the viscosity is assumed constant at η ef f = 9 • 10 5 P a s, the temperature is T = 825 o C, the surface tension is σ = 0.1 N m, the magma density is ρ m = 2354 kg m -3 , the diffusion coefficient is D = 6.95076 • 10 -12 m 2 s -1 . Initial values are from experiment P P E4: the volatile concentration is C i = 0.021756, the porosity is α 0 = 51.2%, the ambient pressure is P a (0) = 40 M P a, and the scaling radius is R i = 111 µm. Experiments P P E1, P P E10, and P P E11 were quenched at 30, 28, and 24 M P a, which correspond to dimensionless times t = 0.25, 0.3, and 0.4, respectively. The initial size distribution from P P E4 was obtained by counting 127 bubbles, which were then sampled in regular intervals of 2.5 µm. Simulations kept that sampling interval with N = 100 discretization points. The other experimental distributions, however, were obtained by counting a very small number of bubbles because they became of a size comparable to the sample. This has two consequences. First, we found preferable to resample both the simulations results and the experimental data to intervals of 12.5 µm for comparison (see figure 4). Second, it precludes the comparison between experimental and simulated porosities. Coalescence led some large bubbles to connect with the sample outer surface, thereby leaving the sample and causing the porosity measure to be underestimated and unreliable. Sample P P E1, for instance, as a measured porosity of 30%, which is below the starting porosity of P P E4 (51.2%), and hence also below that predicted by the three simulations (∼ 61%). Dashed curves are runs with H = 0, thin solid curves are runs with H = H s , and thick solid curves are runs with H = H p . For each dimensionless time (0.25, 0.3, and 0.4), simulated distributions were scaled back to the number of bubbles counted in the respective experiment. The jagged appearance of the experimental distributions is due to the fact that intervals contain an integer number of bubbles. The presence of experimental bubbles smaller than the initial distribution at t = 0.25 can be explained by a smaller nucleation population for that sample rather than bubble breakup.

After a short while (around t ∼ 0.1), the time steps became vanishingly small (< 10 -10 ), which caused prohibitively long calculation times (on the order of several months). These small time steps are due to numerical diffusion from bubble-bearing bins (f n ≫ 0) into near-zero bins (in practice between machine precision and 10 -10 ). A limiting condition so that f n ≥ 10 -10 was thus added at the end of each iteration. This kept the time steps above 10 -5 and the run times to less than an hour on an ordinary PC. We tested the accuracy loss by graphically superimposing the solutions with and without the limiting condition and found them to be negligible. A more accurate estimate is hindered by the fact that time steps are not identical between the runs. It can nevertheless be stated that the percent errors on the marginal with respect to the volume are of the same order as the percent in time differences (i.e. 10 -6 -10 -7 ).

Three simulations are shown on figure 4. The first run has no coalescence (H = 0), in the second run coalescence occurs by planar film thinning (H = H p ), and in the third run coalescence occurs by stretching of inter-bubble films (H = H s ). It is apparent on figure 4 that planar coalescence is more effective than the stretching one because much larger bubbles are produced at all times. The match between planar coalescence and the experimental data is qualitatively the best, although it is clear that experimental distribution are based on too few bubbles to carry out a quantitative comparison. In the absence of more precise experimental data, it can be concluded that the implementation of kernels H based on different physical mechanisms yield distinct evolutions of the size distributions. Our kinetic model is thus a powerful tool that can help assessing how bubble coalescence occurs in magmas.

In Figure 5 we plot the 2D distribution function for the initial data (left) and the final time T = 0.4 (right) in the case of planar coalescence Q p and for the previous set of experimental data. We note that the distribution function is concentrated on the equilibrium manifold defined by the chemical and mechanical equilibria.

Conclusions

The growth of volatile bubbles in very viscous fluids has been widely described in the mono-disperse framework, which is based on the microscopic description of a single bubble.

Starting from the analysis of a mono-disperse bubble population done in [START_REF] Lensky | Bubble growth during decompression of magma: experimental and theorical investigation[END_REF] and [START_REF] Forestier-Coste | Numerical resolution of a mono-disperse model of bubble growth in magmas[END_REF], we build a kinetic model statistically describing the growth by decompression, exsolution and coalescence of a poly-disperse bubble population. The proposed model is homogeneous in space and the distribution function represents the probability to find a bubble of volume v and mass m at time t. The volume and mass growth rates are described in agreement with a simplification of the original microscopic model [START_REF] Lensky | Bubble growth during decompression of magma: experimental and theorical investigation[END_REF] that avoids the resolution of the diffusion of water in the melt. Coalescence is two-dimensional in the proposed model, its rate involving both volumes and masses of the coalescing bubbles, see [START_REF] Castro | Mechanisms and dynamics of bubble coalescence in silicic magmas[END_REF]. Numerical simulations based on a semi-implicit numerical scheme show a good agreement for the coalescence-free case with those of the mono-dispersed model. When coalescence is introduced, numerical results show that coalescence kernels based on different physical mechanisms yields distinct evolutions of the size distributions. In the case presented, experiments underwent extensive coalescence, which yielded poorly defined distribution. This nevertheless allowed us to show a qualitative match of one out of three possible kernels. Our kinetic model is thus a powerful tool that can help assessing how bubble coalescence occurs in magmas.

Several simplifications have been done in order to write this model. For instance, the coupling with an advection-diffusion equation describing the behavior of the volatile concentration in the fluid has been simplified considering a mass growth rate converging to the chemical equilibrium for small values of the relaxation parameter Θ D . Finally, bubbles coalescence has been assumed as instantaneous, that is the created bubble is instantaneously spherical. Of course, this is not real, and it is possible to see the formation of chains of bubbles during the decompression. This is a relevant factor for the modeling of bubbles growth in volcanoes, and need to be taken into account in future works. In fact, it is directly linked to the permeability and to the possibility for the volatile to exit the magma or to accumulate giving rise to an explosive eruption.

Finally, from a mathematical and modeling point of view, it may also be of interest to study the convergence of this kinetic model towards simplified ones. For instance in [START_REF] Forestier-Coste | Numerical resolution of a mono-disperse model of bubble growth in magmas[END_REF] we have deduced several limit situations when the relaxation param-eters Θ V and Θ D converge to zero or infinity, simultaneously or not. Performing an asymptotic analysis on the proposed kinetic model should give simplified equations corresponding to the above limits cases leading to faster computational times.

Figure 1 :

 1 Figure1: A single bubble (white) with its surrounding influence region (hatched). The function F (t, 1) is the volatile flux incoming the bubble from the influence region, ρ m is the magma density, R, S, M and C(t, r) respectively are the radii of the bubble and of the influence region, the bubble mass and the volatile concentration in the influence region.

Figure 2 :

 2 Figure 2: Two coalescing bubbles : planar (left) and stretch (right) profiles.

Figure 3 :

 3 Figure 3: Comparison between the solution of the mono-disperse coupled system (edo gen), the simplified one (edo simpl) and the solution of the kinetic equation (edp).

D = ( 1 , 4 , 7 ,

 147 18, 17, 14, 19, 12, 11, 10, 4, 3, 1, 2, 1, 1, 1).

Figure 4 :

 4 Figure 4: Comparison between numerical runs and experimental data from[START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF]. The initial distributions are shown in light gray and correspond to sample P P E4 and to the initial conditions of the runs. Experimental data are shown in dark gray areas. Dashed curves are runs with H = 0, thin solid curves are runs with H = H s , and thick solid curves are runs with H = H p . For each dimensionless time (0.25, 0.3, and 0.4), simulated distributions were scaled back to the number of bubbles counted in the respective experiment. The jagged appearance of the experimental distributions is due to the fact that intervals contain an integer number of bubbles. The presence of experimental bubbles smaller than the initial distribution at t = 0.25 can be explained by a smaller nucleation population for that sample rather than bubble breakup.

Figure 5 :

 5 Figure 5: Bi-dimensional solution: left -initial data; right -time T = 0.4.

* This work was partially founded by the ERC-starting grant DEMONS (n. 202844) under the European FP7.