L Bretheau

Ç Ö Girit

M Houzet

H Pothier
email: hugues.pothier@cea.fr

D Esteve

C Urbina

Theory of microwave spectroscopy of Andreev bound states with a Josephson junction

Keywords: numbers: 74.45.+c, 73.23.-b, 74.50.+r

We present a microscopic theory for the current through a tunnel Josephson junction coupled to a non-linear environment, which consists of an Andreev two-level system coupled to a harmonic oscillator. It models a recent experiment [Bretheau, Girit, Pothier, Esteve, and Urbina, Nature 499, 312 (2013)] on photon spectroscopy of Andreev bound states in a superconducting atomic-size contact. We find the eigenenergies and eigenstates of the environment and derive the current through the junction due to inelastic Cooper pair tunneling. The current-voltage characteristic reveals the transitions between the Andreev bound states, the excitation of the harmonic mode that hybridizes with the Andreev bound states, as well as multi-photon processes. The calculated spectra are in fair agreement with the experimental data.

I. INTRODUCTION

The Josephson effect, predicted and observed 50 years ago in superconducting tunnel junctions [1][2][3], describes the non-dissipative supercurrent that results from the coherent coupling between superconductors separated by a thin insulating barrier. Since then the supercurrent has been observed in many other weak links such as point contacts, semiconducting nanowires, carbon nanotubes, graphene sheets and thin ferromagnetic layers [4]. Microscopically the Josephson coupling is established through fermionic states whose energies depend on the superconducting phase difference δ across the weak link. For a weak link shorter than the superconducting coherence length ξ, these so-called Andreev bound states (ABS) have energies smaller than the superconducting gap ∆ and are therefore localized at the weak link over a distance of the order of ξ [5][6][7][8]. In a single-channel weak link, there is only one pair of ABS, |-and |+ , with energies ∓E A (see Fig. 1(a)), where

E A = ∆ 1 -τ sin 2 (δ/2). (1)
Here, δ is the superconducting phase difference and τ is the transmission probability for electrons in the normal state. The phase dependence gives rise to opposite supercurrents ∓ (1/ϕ 0 ) (∂E A /∂δ) in the ground and excited states (with ϕ 0 = /2e the reduced flux quantum). This pair of states can be seen as a spin-1 /2 and introduces an internal degree of freedom to Josephson weak links. At zero temperature only the lower energy state |is occupied. This ground state has been probed through measurements of the current-phase relation in superconducting atomic contacts [9]. The direct observation of the higher energy state |+ has been achieved only recently [10,11]. In the experiment of Ref. [10], a voltagebiased Josephson junction was used as an on-chip spec-trometer (see Fig. 2). The dissipative subgap current through the Josephson junction is due to inelastic tunneling of Cooper pairs, the released energy being absorbed in the environment (see Fig. 1). Current peaks are then observed at energies corresponding to the eigenenergies of the environment. In practice, the environment was an asymmetric SQUID formed by a superconducting atomic contact in parallel with a second Josephson junction (see Fig. 2). This environment can be modeled as a spin-1/2 degree of arXiv:1406.6301v2 [cond-mat.mes-hall] 10 Oct 2014 freedom (the Andreev doublet) coupled with a harmonic mode (the plasma mode of the Josephson junction) (see Fig. 1). The goal of the present work is to reach a quantitative understanding of the measured spectra.

Andreev doublet E

A -E A Δ -Δ 0 (a) δ E J hν = 2E A 2π 2 1 Δ -τ - + (
The rest of the paper is organized as follows. Section II gives a Hamiltonian description of the circuit schematized in Fig. 2. The spectrometer Hamiltonian is then treated as a perturbation, and following previous work [12][13][14][15][16][17], the atomic contact is treated under the assumption of small phase fluctuations. In section III, the resulting spin-boson Hamiltonian of the environment is solved numerically, and the calculated current spectra are compared with the experimental ones. Finally, in section IV the Hamiltonian is solved analytically in the Jaynes-Cummings approximation. Some perspectives are discussed in the conclusion.

II. MODEL

A. Setup

The setup used in the experiment [10] is shown schematically in Fig. 2. On the left-hand side, a voltagebiased Josephson junction (critical current I 0 = 48 nA) is used as a spectrometer. It is biased with a voltage source V in series with a resistor R = 2 kΩ. It radiates microwaves at the Josephson frequency ν = 2eV /h [18], which can be absorbed by its electromagnetic environment, an atomic-SQUID formed by an atomic point contact in parallel with a second Josephson junction. The critical current I L = 1.06 µA of this second junction is much larger than those of the spectrometer junction and of a one-atom contact. An external superconducting coil is used to apply a dc flux φ = ϕ 0 ϕ through the SQUID loop. The capacitance C = 280 fF corresponds to the sum of the junctions capacitors.The atomic-SQUID and the spectrometer are coupled through a capacitor Σ ∼30 pF. Whereas it behaves as an open circuit from the dc point of view and ensures that the dc voltage V falls on the spectrometer, it can be considered as a short for the Josephson radiation in the explored frequency range (ν ν p (C/Σ) 1/2 ∼ 2 GHz, with ν p the plasma frequency of the SQUID junction).

B. Hamiltonian of the circuit

Neglecting the coupling capacitor and the bias resistor, the circuit represented in Fig. 2 can be described with the Hamiltonian H = H spec + H SQ . The first term H spec = -E J cos (γ+2eV t) corresponds to the spectrometer with Josephson energy E J = ϕ 0 I 0 . The voltage drop across the spectrometer junction induces a superconducting phase difference α = γ + 2eV t, with γ the phase at The Hamiltonian

H SQ = E C N 2 -E L cos (γ) + H A (δ) (2)
accounts for the SQUID. The operator N is the number of Cooper pairs that have crossed the tunnel junction; N and γ are conjugated operators:

[γ, N ] = i. E C = 2e 2 /C
is the charging energy for pairs and E L the Josephson energy of the SQUID junction. The last term in Eq. ( 2) describes the atomic contact. The phase δ = ϕ+γ across the atomic contact differs from γ by ϕ, the reduced flux. For a single channel contact, and neglecting excitations which involve quasiparticles in the continuum, the Andreev Hamiltonian is [13,14]

H A = -∆ (ImU σ y + ReU σ z ) , with U = cos( δ 2 ) + i √ 1-τ sin( δ 2 ) e -i √ 1-τ δ 2 (3)
where the Pauli matrices σ x,y,z act in a 2 × 2 subspace corresponding to the two ABS. The physics of the ABS is therefore analogous to that of a spin-1 /2 in a magnetic field whose magnitude ∆ |U | = E A and direction depend on the superconducting phase difference δ across the contact. The eigenstates of H A are the ABS |± , with eigenenergies ±E A (see Eq. ( 1) and Fig. 1(a)).

C. Andreev spin and plasma boson

In the experiment [10], due to the large asymmetry E L E A , the phase dynamics is essentially ruled by the SQUID junction. The size of the phase fluctuations is determined by the dimensionless parameter

z = (E C /2E L ) 1/2
1. Therefore, following Ref. [12], we may treat the SQUID junction as a linear inductor and retain only the lowest order coupling (∝ √ z) between the SQUID junction and the atomic contact in Eq. ( 2). Namely,

H SQ ≈ E C N 2 + E L γ 2 /2 + H A (ϕ) + ϕ 0 γC A (ϕ) , (4)
where C A = ϕ -1 0 ∂H A /∂δ is the Andreev current operator. Hence, the parallel combination of the SQUID junction and the capacitor C forms a harmonic oscillator of resonant frequency ν p = √ 2E L E C /h (see Fig. 1(c)). Switching to second quantization, its Hamiltonian reads hν p (a † a + 1 2 ), where the phase γ = √ z a + a † is linked to the annihilation and creation operators a and a † of the plasma mode of the SQUID.

In the basis of the Andreev states {|-, |+ }, the SQUID Hamiltonian finally reads

H SQ = hν p a † a + 1 2 -E A σ z + a + a † (Ω x σ x + Ω z σ z ) , (5) where
Ω z = ∆ √ z τ sin (ϕ) 1 -τ sin 2 ϕ 2 (6)
and

Ω x = Ω z √ 1 -τ tan ϕ 2 .
The phase-dependence of these coupling energies is represented in Fig. 3. The spin-boson Hamiltonian (5) yields a discrete eigenspectrum which results from the hybridization of the bosonic plasma mode in the SQUID junction and the Andreev spin-1/2 degree of freedom in the atomic contact (see Fig. 4(a)). The term ∝ Ω x allows transitions between Andreev states. It is minimum at ϕ = 0 and maximum at ϕ = π. Note that the flux modulation of the SQUID plasma frequency due to the contribution of the effective inductance of the atomic contact would appear as a higher order effect in z.

0 π 2 π 0 0.5 -2 1-1-τ 2 τ ϕ -0.5 Ω x Δ z Ω z Δ z

D. Incoherent Cooper pair current

The global circuit Hamiltonian H cannot be diagonalized analytically. Along the lines of P (E) theory for dynamical Coulomb blockade [20,[START_REF] Ingold | Single Charge Tunneling[END_REF], the spectrometer may be treated as a perturbation. The dc current I flowing through the spectrometer is calculated using the current operator I 0 sin(γ + 2eV t) and perturbation theory up to the second order in E J . At zero temperature and for V > 0,

I (V ) I 0 = π 2 E J k k| e iγ |0 2 δ ((E k -E 0 )-2eV ) , (7)
where |k are the eigenstates of H SQ with energy E k . Equation ( 7) accounts for incoherent Cooper pair tunneling through the junction's barrier at rate I/2e, provided that the energy 2eV matches an excitation energy of the junction environment (see Fig. 1(b)). As in P (E) theory for zero temperature, we assume that, before each tunneling event, the environment has relaxed to its ground state |0 .

In P (E) theory [20,[START_REF] Ingold | Single Charge Tunneling[END_REF], the environment is purely electromagnetic, linear and is described by an impedance Z(ω). Equation ( 7) then simplifies and gives the current as a function of ReZ (ω). In particular, this would yield current peaks at 2eV = nhν p (n integer), with a width related to the quality factor of the bosonic plasma mode of frequency ν p . Such an approach is not possible here due to the strong non-linearity of the Andreev degree of freedom. As the model does not include dissipation, the spin-boson Hamiltonian H SQ yields a discrete spectrum. Then, Eq. ( 7) predicts infinitely sharp dc current peaks at the excitation energies of H SQ . A finite broadening of each of the peaks may be introduced phenomenologically with the substitution δ

(E -E 0 ) → (Γ/π) / (E -E 0 ) 2 + Γ 2
in Eq. (7). While the model is unable to predict the amplitude of the linewidth Γ, it should be large enough for the incoherent Cooper pair current through the spectrometer junction to be small, I I 0 , so that perturbation theory is valid. Near voltage V = k /(2e), where k = E k -E 0 , this condition reads (E J /Γ)P k 1, with transition probability P k = k| e iγ |0 2 .

III. SPECTRA CALCULATION A. Numerical resolution

The spectrum of Hamiltonian ( 5) can be obtained numerically. To do so, we write it as a matrix in the basis {|σ, n }, where σ = ± accounts for the Andreev spin-1/2 and n is the plasmon occupation number, and truncate to the lower energy states. By numerical resolution of a 14 × 14 Hamiltonian matrix (n ≤ 6), we have computed all the resonances of energy smaller than 2∆ and the corresponding transition probabilities P k . Figure 4 shows the excitation spectrum (b) and the transition probabilities (c) for the transitions indicated in (a), for a channel of transmission τ = 0.98. In the experiment, the superconducting gap energy is ∆/h 43 GHz, the charging energy is E C /h 270 MHz and the Josephson energy of the SQUID junction is E L /h 900 GHz, leading to z 0.012. Note that E L is renormalized as E L = ϕ 0 (I L + 2ϕ 0 /L), where L/2 = 0.44 nH is a parallel inductor (not shown in Fig. 2) accounting for a 1.35 mmlong aluminum connecting wire present in the actual geometry of the sample.

B. Incorporating line broadening

The experimental data displayed in the right panels of Fig. 5 show broad lines. One notices that the transition solely concerning the plasma mode is already broad, which indicates that this mode is subject to dissipation and/or dephasing; note also that voltage fluctuations across the spectrometer junction can limit the energy resolution [START_REF] Devoret | Likharev Dynamics of Josephson Junctions and Circuits[END_REF]. To account for the broadening of the energy levels of H SQ one should at least include the coupling of the plasma mode to a dissipative bath. Treating the latter as an infinite collection of harmonic oscillators, as it is usually done, would lead to a complex spin-boson problem without analytical solution. Here we simply introduce in Eq. ( 7) a phenomenological Lorentzian broadening, cf. Sec. II D. By fitting the plasma resonance peak, we found Γ/h = 2 GHz, which we used for all the calculated lines shown in the left panels of Fig. 5.

C. Comparison with experiment

In experiments atomic contacts have several channels. We extend the model to contacts with multiple channels by adding an Andreev term (3) for each channel to the Hamiltonian. This Hamiltonian describes the physics of N spins coupled to the same bosonic mode but not directly coupled to each other. All the approximations made before are still valid [23]. Figure 5 AC1 (transmissions 0.942, 0.26) and AC2 (transmissions 0.985, 0.37). The theory is obtained with the numerical method sketched before, still using seven photon levels (28 × 28 matrix).

The model describes both the Andreev transitions |-, 0 → |+, 0 , of energy 2E A (V-shaped lines), and the plasma transition |-, 0 → |-, 1 , of energy hν p (red horizontal line at 0.51 ∆). It also describes the higher harmonic transitions: |-, 0 → |+, 1 , of energy 2E A + hν p (replica of the Andreev transition, shifted up by 0.51 ∆) and |-, 0 → |-, 2 , of energy 2hν p (white horizontal line at 1.02 ∆). These processes correspond to the tunneling of one Cooper pair and emission of two photons. They are less probable and result in fainter transitions as seen both in experiment and calculation. Theory also accounts for the anti-crossings arising from the coupling between the Andreev-spin and the plasma-boson. Only crossings of transition lines involving the same number of photons show significant hybridization, in good agreement with 3 in Ref. [10]. The grey regions at ν < 4 GHz and around 25 GHz are not accessible because there the biasing of the spectrometer is unstable. Left: The calculated spectra are obtained using Eq. ( 7) with a phenomenological damping parameter Γ/h = 2 GHz. The calculated currents have been multiplied by 0.6 to match the measured ones. The transition probabilities and excitation energies are computed numerically, using 7 photon levels in the plasma mode.

the data.

Extra phase-independent current at 2eV ∆/2 in the experimental data is attributed to the coupling to uncontrolled environmental modes outside of the SQUID.

Finally, the global weakening of the signal at high V is well captured by the model. This weakening is due to the lower impedance of the SQUID capacitance compared to its inductance at frequencies larger than the plasma frequency. However, the agreement for the amplitude and width of the different peaks is not quantitative. Having set the width to 2 GHz, a correct peak amplitude is only obtained when multiplying theory by a factor 0.6. A rigorous treatment of dissipation is needed. This could be achieved by adding an electromagnetic impedance with a dissipative component in parallel with the SQUID.

IV. ANDREEV-PLASMA MODE HYBRIDIZATION

We use the rotating wave approximation, valid in the vicinity of the degeneracy hν p = 2E A , to obtain analytical expressions for the excitation energies and transition rates. The Hamiltonian can be approximated by the Jaynes-Cummings model [START_REF] Jaynes | Proc. IEEE[END_REF][START_REF] Shore | [END_REF]:

H JC SQ = hν p a † a + 1 2 -E A σ z + Ω x aσ + + a † σ - (8)
where

σ ± = 1 2 (σ x ∓ iσ y ).
Then, by blockdiagonalization in the subspace {|+, n , |-, n + 1 }, one derives the eigenstates and eigenenergies and computes the excitations energies and transition probabilities.

To fourth order in √ z, the current (7) through the spectrometer displays four peaks

I(V ) I 0 = π 2 E J σ=±;n=1,2 P σ n δ (2eV -σ n ) (9)
at the bias voltages matching the excitation energies (see top panel of Fig. 6):

± n = E A + n - 1 2 hν p ± W n (10)
with

W n = (E A -1 2 hν p ) 2 + nΩ 2
x . The amplitudes of the current peaks are proportional to the transition probabilities (see bottom panel of Fig. 6):

P - n = 1 n! z n e -z cos 2 θn 2 P + n = 1 n! z n e -z sin 2 θn 2 (11)
with θ n = arctan (nΩ x /W 0 ). These peaks correspond to excitations towards composite states, resulting from the hybridization of the plasma and Andreev modes (see Fig. 4(a)). The first two resonances in Eq. ( 9) correspond to the excitation with a single photon of the hybridized Andreev (2E A ) and plasma mode (hν p ); the last two resonances correspond to the excitation of higher harmonic modes at 2hν p and hν p + 2E A . These two-photon processes are possible because the spectrometer, a Josephson tunnel junction, is a nonlinear emitter. They correspond to the tunneling of one Cooper pair and emission of two photons. Note that far from degeneracy, using perturbation theory with z 1, one finds for the plasma resonance and its harmonic the transition probabilities ze -z and 1 2 z 2 e -z . These amplitude are consistent with the Poisson distribution found at arbitrary z in the P (E) derivation with a bosonic mode [START_REF] Ingold | Single Charge Tunneling[END_REF]. It is worth mentioning that the Jaynes-Cummings model gives an excellent description, even far from degeneracy when |W 0 | hν p . There, the error made in energy scales as Ω x / (2E A + hν p ), and is in general negligible. However, this model fails to predict the anticrossing between the second harmonic of the plasma resonance and the Andreev resonance at 2hν p = 2E A (red and green lines in Fig. 6). In practice this correction is one order of magnitude smaller than the anti-crossing at hν p = 2E A and is not seen in the experiment [10]. Note however that the derivation in the case of multiple chan-nels is much more involved when more than one Andreev transition crosses the plasma resonance.

V. CONCLUSIONS AND PERSPECTIVES

In conclusion, we have treated inelastic Copper pair tunneling through the emitter junction using perturbation theory in its coupling with the environment, and calculated the environmental transition energies and probabilities by diagonalization of the SQUID Hamiltonian. This Hamiltonian contains a spin-like term describing the atomic contact and which is coupled to a harmonic term for the plasma mode of the SQUID junction. The calculated spectra are in fair agreement with the experimental ones. However this theory could benefit from several extensions.

First of all, the atomic contact Hamiltonian was restricted to a two-level system Hamiltonian, by neglecting excitations which involve quasiparticles in the continuum. In particular, quasiparticle trapping, which was measured in a similar experiment [11], leaves the Andreev doublet in a long-lived "odd" state for which both the Andreev states are either occupied or empty [26], with time scales in the millisecond range [27]. It limits the Andreev transition rate and therefore reduces the corresponding dc current. On the other hand, transitions from the Andreev bound states to the continuum were measured on the same contacts using another detection method [11]. They do not contribute to the dc current and were not seen in the experiment described in Ref. [10] due to the long lifetime of the odd states. To deal with all the Andreev excitations on the same footing, one must consider the full Hamiltonian of the atomic contact, as done in recent works [17,28].

Second, the description used here for the environment eigenstates and for the spectrometer is far from being complete. As discussed above, this is a difficult problem for two reasons: (i) the coupling of a two-level system to a dissipative resonator is not known analytically, and (ii) the I(V ) characteristic of a Josephson junction coupled to a resonator is not generically known when the perturbative expansion in E J is not valid, despite recent progress along this direction [29]. Although the second order expansion is valid for the data discussed here since the Cooper pair current is small (I I 0 ), higher order contributions corresponding to the transfer of more than one Cooper pair emitting several photons were observed in other experiments [11,[START_REF] Bretheau | Localized Excitations in Superconducting Atomic Contacts: Probing the Andreev doublet[END_REF].

Finally, it was also assumed in the model that the environment is in its ground state each time a Cooper pair tunnels. This requires the relaxation time of the excited environment states to be shorter than the inverse tunnel rates. When this is not the case, the environment modifies the tunneling process. For the harmonic oscillator bosonic mode, this could lead to stimulated emission and lasing. For the Andreev two-level system, this could saturate absorption and reduce the Cooper pair current through the spectrometer junction. A regime in which the Andreev two-level system would coherently exchange an excitation with the Josephson junction can also be envisioned.

Figure 1 :

 1 Figure 1: (Color online) Principle of the photon spectroscopy of Andreev bound states. (a) Phase (δ) dependence of the Andreev levels with energies ∓EA in a short conduction channel of transmission τ. ∆ is the superconducting gap. (b) A Cooper pair tunneling across the spectrometer junction releases energy 2eV as a photon of frequency ν, which is absorbed either in the atomic contact by exciting the Andreev transition at energy 2EA (a) or in a harmonic oscillator mode (c) present in the embedding circuit of the atomic contact (see Fig. 2).

Figure 2 :

 2 Figure 2: (Color online). Simplified diagram of the experimental set-up. A voltage-biased Josephson junction (orange cross) is used as a spectrometer. It emits microwaves in its environment, an atomic-SQUID formed by an ancillary Josephson junction (green cross in parallel with a capacitor C) and an atomic point contact (blue triangles) of channel transmission probabilities {τi}. The absorption of a photon by the environment is accompanied by the transfer of a Cooper pair through the spectrometer. The phases γ and δ across the SQUID junction and the atomic contact are linked by the external reduced flux ϕ threading the loop: δ = γ + ϕ.The phase across the spectrometer is given by α = γ + 2eV t.

Figure 3 :

 3 Figure 3: (Color online). Phase dependence of the coupling energies Ωx (red) and Ωz (blue) in units of √ z∆, for τ = 0.8.

Figure 4 :

 4 Figure 4: (Color online). (a) Energy spectrum diagrams of HSQ for a single channel: each state is labelled |-, n or |+, n for the Andreev pair in the ground (-) or excited (+) state and n photons in the plasma mode. At degeneracy 2EA = hνp, the plasma and Andreev modes hybridize, which leads to an avoided crossing visible in (b). (b, c) Excitation energies (in units of ∆) (b) and transition probabilities P (c) of the first resonances, as a function of the reduced flux ϕ, for a channel of transmission τ = 0.98. These lines are obtained by numerical resolution of a 14 × 14 Hamiltonian matrix. Each color encodes a transition from the ground state towards a different excited state, as represented in (a).

Figure 5 :

 5 Figure 5: (Color online). Comparison between calculated (left) and experimental (right) spectra I (ϕ, V ) for contact AC1 (transmissions 0.942, 0.26) (top) and AC2 (transmissions 0.985, 0.37) (bottom). The two insets show theoretical (black lines) and experimental (green lines) I (V ) cuts for ϕ = π, for each contact. Right: The experimental spectra are extracted from Figure3in Ref.[10]. The grey regions at ν < 4 GHz and around 25 GHz are not accessible because there the biasing of the spectrometer is unstable. Left: The calculated spectra are obtained using Eq. (7) with a phenomenological damping parameter Γ/h = 2 GHz. The calculated currents have been multiplied by 0.6 to match the measured ones. The transition probabilities and excitation energies are computed numerically, using 7 photon levels in the plasma mode.

Figure 6 :

 6 Figure 6: (Color online). Excitation energies (in units of ∆) (top) and transition probabilities P (bottom) of the first resonances, as a function of the reduced flux ϕ, in the Jaynes-Cummings approximation for a channel of transmission τ = 0.98. Each color encodes a transition towards a different state, as shown in Fig. 4(a).

VI. ACKNOWLEDGMENTS

We thank Sebastian Bergeret, Juan Carlos Cuevas, Andrew C. Doherty, Alfredo Levy Yeyati, Fabien Portier, Yuli Nazarov and Vitaly Shumeiko for enlightening discussions. We gratefully acknowledge help from other members of the Quantronics group. Work partially financed by ANR through projects ANR-09-BLAN-0199-01 DOC-FLUC, ANR-12-BS04-0016 MASH and ANR-11-JS04-003-01 GLASNOST. The research leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme (FP7/2007-2013) under REA grant agreement n°PIIF-GA-2011-298415.


Pothier

G Urbina

D G Olivares

A Levy Yeyati

L Bretheau

C ¸ Ö Girit

H Pothier

C Urbina

Dynamics of

I. INTRODUCTION

There are several external mechanisms that undermine the quantum coherence of superconducting circuits being explored for quantum information processing 2 . Their influence has been reduced over the years by new designs that minimize the coupling with external degrees of freedom. However, a fundamental intrinsic decoherence process arises from the coupling of the qubit variables to superconducting quasiparticles tunneling through the Josephson junctions of the circuits. Although in principle the superconducting gap ∆ provides an inherent protection against low energy excitations at low temperatures, in practice there are residual nonequilibrium quasiparticles that can rule the behavior of the circuits [3][4][5][6][7][8][9] . As shown in a recent experiment 10 , this is particularly true for weak links containing channels of high transmission, where localized excitations occupying Andreev levels of energy below ∆ become possible. This has important consequences for the corresponding proposed qubits designs [11][12][13][14] . Furthermore, single quasiparticle trapping in localized levels could be detrimental in experiments proposed to detect "Majorana bound states" in condensed matter systems since their topological protection relies on parity conservation 15 . Understanding the dynamics of relaxation of quasiparticles in superconducting weak links is therefore an important present-day issue.

We report here on a theory highlighting the role of the electromagnetic environment in this dynamics. In experiments, Josephson junctions or weak links are very often embedded in electrical circuits having electromagnetic modes at frequencies lower than the superconducting gap. The environment can be a resonator intentionally coupled to the junction like in Ref. 10, or the plasma mode of another junction placed in parallel like in Ref. 1. We show that if the mode impedance is large enough, it rules the quasiparticle dynamics when the sum of the An-dreev level energy and of the energy of the mode exceed the superconducting gap.

The rest of the paper is organized as follows: In Sec. II we describe the model considered for a superconducting one channel contact coupled to a generic electromagnetic environment; Sec. III is devoted to the analysis of the transition rates between different quasiparticles states induced by quantum phase fluctuations; in Sec. IV we focus on the experimental situation of Ref. 1 and compare the theoretical results for the transition rates and the stationary probability for quasiparticles trapped in the subgap states with the corresponding experimental results. In Sec. V we present our main conclusions. The more technical details on our calculations are described in appendices A, B, C and D.

II. THEORETICAL FRAMEWORK

We consider the situation illustrated in Fig. 1(a) with a single superconducting channel (SC) coupled to an arbitrary impedance Z env . The excitation spectrum of the arXiv:1309.7283v3 [cond-mat.mes-hall] 12 Feb 2014 SC contains a discrete, spin-degenerate Andreev level, with an energy E A (δ) = ∆ 1 -τ sin 2 δ/2, where δ is the superconducting phase difference across the contact and τ the transmission probability for electrons 16 (see Fig. 1(b)). The Andreev level is completely empty when the channel is in its ground state, which has a phase dependent energy -E A and carries a supercurrent I = -(∂E A /∂δ) /ϕ 0 , where ϕ 0 = /2e is the reduced flux quantum. The lowest-energy excitations correspond to the occupation of the Andreev level by a single quasiparticle (of either spin), the global energy and the supercurrent of these "odd" configurations being then zero. There is also an excitation of energy 2E A with respect to the ground state, where the Andreev level is occupied by two quasiparticles of opposite spins. This "even" configuration can be seen as a localized excited "Andreev pair" 17 , and carries a supercurrent opposite to that in the ground state.

The system Hamiltonian can be written as Ĥ = ĤSC ( δ) + Ĥenv (γ), where the first term describes the SC and the second one the electromagnetic environment, modelled by the impedance Z env . The phases δ and γ across the SC and the impedance are related by δ -γ = Φ/ϕ 0 = ϕ, where Φ is the magnetic flux through the loop.

The population of the SC electronic states becomes then sensitive to the effects of quantum phase fluctuations. Assuming that Re(Z env ) R Q , we treat the fluctuations to lowest order in perturbation and write the Hamiltonian Ĥ = Ĥenv (γ) + ĤSC (ϕ) + ϕ 0 γ Î(ϕ), where Î = ϕ -1 0 ∂ ĤSC /∂δ is the current operator in the contact region.

To describe the unperturbed single-channel SC we use a one-dimensional SNS junction model with a Dirac delta potential barrier (to account for non-perfect transmission) inside a normal region of negligible length. Details of the diagonalization of this model in terms of Bogoliubov fermion operators γ α,σ , where σ indicates spin, are given in Appendix A. Two types of states are obtained, α ≡ k with energy E k ≥ ∆ corresponding to the extended continuum states and α ≡ A corresponding to the localized Andreev states with energy E A (see Fig. 1 (b)). The SC ground state |Ψ 0 corresponds to the absence of excitations, i.e. γ α,σ |Ψ 0 = 0.

III. TRANSITION RATES

The coupling of the SC to the environment allows for transitions between different quasiparticle states. We shall first consider processes which permit the removal of a quasiparticle from the Andreev level. These processes allow in particular the relaxation of the lowestenergy excited states with one trapped quasiparticle back to the ground state |Ψ 0 [START_REF]Similar processes but for photons from an external driving field were considered in Ref[END_REF] . They consist either in the absorption of an environmental photon and transfer of the trapped quasiparticle into the continuum states, or in the recombination of a quasiparticle from the continuum with the trapped one into a Cooper pair while releasing the energy as a photon. These two processes are illustrated in panels (a) and (b) of Fig. 2

= 2π k k, σ ϕ 0 Î A, σ 2 P (E k -E A (δ)) × (1 -f FD (E k , T qp )) , (1)
where f FD (E, T qp ) is the Fermi population factor for quasiparticles in the continuum (assumed to be in equilibrium at a temperature T qp ) and P (E) is the probability of absorbing a photon of energy E from the environment. This probability is

P (E) = D (E) f BE (E, T env ),
where f BE (E, T ) is the Bose population factor, and

D (E) = Re {Z env (E) /E} /R Q , with R Q = h/4e 2
, is the density of states for the modes in the environment [START_REF] Ingold | Nazarov in Single Charge Tunneling[END_REF] .

The environment is assumed to be in equilibrium at a temperature T env which can be in general different from T qp . The numerical evaluation of this rate (and every other) for different transmissions, shows a rather universal dependence in the Andreev level energy position E A .

Simple analytical expressions can be derived in the perfect transmission limit τ → 1 and in the tunnel limit τ → 0, for which the wavefunctions have a considerably simpler form. In this limit one obtains (see Appendix B)

Γ (a) out = 8∆ h ∞ ∆ dED (E -E A ) g (E, E A ) ×f BE (E -E A , T env ) (1 -f FD (E, T qp )) , (2)
with g (E,

E A ) = (E 2 -∆ 2 ) (∆ 2 -E 2 A )/(∆ (E -E A )
). We also give in Appendix B the expression of g (E, E A ) in the tunnel limit τ → 0.

When the environment of the SC contains a single mode with infinite quality factor, and at low temperature, this expression simplifies to Γ

(a) out = 2∆ Z0 R Q g (E A + hν, E A ) exp (-hν/k B T env )
, where ν is the mode frequency. The function g is of order 1 when ∆ -hν < E A < ∆, so that the rate is simply determined by the impedance Z 0 of the oscillator. For aluminum and for Z 0 = 50 Ω, 2∆ Z0 R Q ≈ 1 GHz. In a similar way for the second relevant process we find

Γ (b) out = 8∆ h ∞ ∆ dED (E + E A ) g (E, -E A ) × (1 + f BE (E + E A , T env )) f FD (E, T qp ) . (3)
We show in Appendix D that for perfect transmission the matrix elements for electron-phonon coupling have the same functional form in terms of E and E A as those for the coupling with the electromagnetic modes. Therefore its inclusion leads to the same expressions for Γ

IV. COMPARISON TO EXPERIMENTS

We focus on the recent experiments on superconducting atomic contacts 1 that have analyzed in detail the quasiparticle trapping in Andreev levels and its dynamics. In these experiments, an atomic contact was embedded in a superconducting loop containing a Josephson junction, thus forming an asymmetric SQUID. It was found that there is a significant probability for the SC to get trapped in an odd state in which the highest transmitted channel carries no supercurrent. The experiments also showed that the relaxation rates for these states fall into a nearly universal behavior as a function of the energy E A regardless of the particular values of the transmission and phase difference. Trapping occurred essentially when the Andreev level energy was smaller than half the superconducting gap E A 0.5∆, with the lifetime of trapped quasiparticles exceeding 100 µs. For larger energies no significant trapping could be detected. The origin of this sharp energy threshold was a puzzle not explained in the paper reporting the experiment.

In Ref. 1, the SQUID Josephson junction had a Josephson energy much larger than the charging energy, and it can therefore be described as an harmonic oscillator. Spectroscopy measurements 17,[START_REF] Bretheau | [END_REF] on similar circuits as the one used in Ref. 1 showed that the plasma frequency of this mode can be significantly renormalized by parallel inductances and approach ∼ 0.5∆ (see Appendix C and Ref. 19). As explained in the following, our theory shows that the main relaxation mechanism for the trapped quasiparticle states is their excitation into the extended continuum states above the superconducting gap by absorption of photons from the plasma mode. This mechanism becomes inefficient when the energy difference between the Andreev level and the continuum exceeds the plasma energy, ∆ -E A > hν p , hence providing a simple explanation for the observed behavior.

The results for the transition rates Γ

in is large in this energy range because quasiparticles near the gap edge can relax in the Andreev level while emitting a plasma photon (see Fig. 3). For lower energies, the energy of the plasma photons is not sufficient and the rate drops abruptly. Other processes, like phonon absorption or emission start to play a role. Hence, both Γ in and Γ out are determined by phonon processes for E A < ∆ -hν p and by photon processes for E A > ∆ -hν p . It should be noticed that three different temperatures enter the calculation. We assume that the phonons in the Al films of Ref. 1 are at equilibrium with the substrate and therefore T ph is taken equal to the base temperature measured by the thermometers in the experiment (30 mK). The two other temperatures, T env and T qp , can be significantly larger due to incomplete filtering of radiation. To fit the results we have used T env ∼ 120 mK, similar to what is deduced from measurements of the switching probability of the SQUID [START_REF] Bretheau | [END_REF]21 and T qp ∼ 180 mK which simulates the presence of a few tens of out-of-equilibrium quasiparticles per µm 3 , as typically found in experiments with Al resonators and qubits [3][4][5]7,[START_REF] Visser | [END_REF] .

The different transitions which determine the population and relaxation of the Andreev level are illustrated in the inset of Fig. 4. They also involve the even excited state |even * = γ † A↑ γ † A↓ |Ψ 0 . The analysis is further simplified by symmetry relations: the rates connecting the even excited state and the odd states are equal to the ones connecting the odd states and the even ground state. This is indicated by the color code used for the arrows in the inset of Fig. 4. Notice that the full determination of the level populations requires also the evaluation of the rates Γ e * →e and Γ e→e * . In Ref. 1 it was assumed that the relaxation rate Γ e * →e to the ground state was very fast compared to all other ones and that Γ e→e * was negligible. For photonic or phononic environments these have been calculated in Refs. 11,19 and 23 respectively and reproduced with the present formalism, as discussed in Appendix B and D. The calculation does corroborate that Γ e * →e is much larger than all the other rates for the transmissions explored in the experiment.

E + E A ! E E A ( b ) E -E A ! E A E ( a )
A last step in our calculation is to obtain the stationary distribution of quasiparticles by solving the master equation involving all transitions indicated in the inset of Fig. 4. The result for the occupation probability of the odd states, P ∞ , is shown in Fig. 4 and compared with the experimental results from Ref. 1 for contacts with different transmissions. As can be observed, the theory qualitatively describes the decrease in P ∞ at fixed E A which is observed experimentally for increasing transmission in the slow relaxation regime. There is, however, some discrepancy in the quantitative values of P ∞ which is overestimated in our model calculations.

We show in Fig. 5 the comparison of the experimental and theoretical results for the total rates Γ in and Γ out as a function of E A for different values of the contact transmission. One should remark the quite good agreement which is obtained for Γ out in the fast relaxation regime (E A > 0.5∆) and for Γ in in the slow relaxation regime. The drop in Γ out by more than two orders of magnitude at E A ∼ 0.5∆ is also captured by our model. In the slow relaxation regime the model correctly describes the decrease of Γ in and the increase in Γ out which is observed at fixed E A with increasing transmission. The largest discrepancies between model and experiment are found for Γ in when E A 0.5∆ and for Γ out when E A 0.5∆. It should be noticed, however, that the experimental determination of Γ in is less precise for E A > 0.5∆ where P ∞

1. The decay of Γ in in theory is essentially related to the assumption of a thermal distribution of quasiparticles, which leads to Γ in ∝ exp(-(E A + hν p )/k B T qp ). A weaker dependence on energy would be obtained with a self-consistent description of this distribution 4 . On the other hand, the deviations found for Γ out in the slow relaxation regime are more significant and could indicate that some additional relaxation mechanism, like tunneling to vortices or quasiparticle traps in the vicinity of the contact, could be contributing for small E A . Since in the limit where Γ e * →e Γ in , Γ out Γ e→e * , the stationary occupation probability P ∞ is simply given by 2Γ in /(3Γ in + Γ out ) 21 , the deviations for Γ out explain why the theory overestimates P ∞ in this regime. Further combined experimental and theoretical work would be required to clarify this point.

V. CONCLUSIONS

In conclusion, we have presented a theory which describes the dynamics of trapping and untrapping quasiparticles in phase-biased superconducting weak links. It is shown that in realistic conditions this dynamics can be controlled by the coupling of the weak link to its electromagnetic environment. The results are in semiquantitative agreement with the experiments of Ref. 1, where the sharp jump observed in the trapping and untrapping rates is associated to the onset of the coupling to the environment plasma mode. The mechanisms described here can be relevant for controlling decoherence in superconducting qubits involving channels with nonnegligible transmission. In case of the Andreev qubits discussed in Refs. 11,12, where poisoning by trapped quasiparticles in the ABS should be avoided, the presence of a mode of energy larger than ∆ -E A would be beneficial. In contrast, for the proposals of Refs. 13,14 which are based on the manipulation of the odd states, a larger lifetime of the trapped quasiparticles is desirable. In this case one would need an electromagnetic environment containing no mode of frequency larger than ∆-E A . Finally, as discussed in Ref. 6, even in the case of qubits based on tunnel junctions changes in the occupation of the Andreev states make the Josephson coupling and hence the qubit frequency fluctuate thus giving rise to dephasing. Therefore, even in this case of Andreev levels very close to the gap edge, slowing down the dynamics of these occupations could have an influence on the qubit decoherence.

The point contact is modelled as a 1D SNS junction with a Dirac delta barrier in the normal region.

The normal region's length LN can be taken to the limit LN → 0 in the ballistic regime, effectively turning the scattering problem into the problem of a delta barrier in a superconducting system, with a well-defined phase bias between the left and right leads.

Continuum wavefuntions: For an homogeneous superconducting system, the Hamiltonian's eigenfunctions

have a momentum |k| = k F ± κ E , where κ E = ξ -1 0 (E/∆)
2 -1 (ξ 0 being the superconducting coherence length, given by hv F /∆ in the ballistic regime). The eigenfunctions with positive energy E k and spin-up take the following shape:

ψ (↑) k, φ (x) = 1 √ L u E v E e i φ e +ikx for |k| > k F (quasielectrons) ψ (↑) k, φ (x) = 1 √ L v E u E e i φ e -ikx for |k| < k F (quasiholes) (A1)
where φ is the phase of the superconducting order parameter, L is a length over which the freely propagating eigenfunctions are defined and u E , v E are, respectively, the electron and hole components of the plane wave: These plane waves are schematically shown in Fig. 6.

(u/v) E = 1 √ 2 1 ± 1 -(∆/E) 2 (A2) E κ k k F -k F Δ +κ -κ +κ -κ u! v! ( ) u! v! ( ) v! u! ( ) v! u! ( )
Summations in momenta such as the one in Eq. ( 1) may be rewritten, for the sake of convenience, as integrals over quasiparticle energies weighted by the superconducting density of states ρ SC

(E) = ρ F |E| √ E 2 -∆ 2
, where ρ F is the normal density of states at the Fermi level.

That being the case we shall brand the wavefunctions (and the states they refer to) not by using their momentum k as an index, but their energy E k → E, their associated quasielectron/hole character and the direction of their momentum.

With the plane waves from Eq. (A1) it is possible to construct solutions to the BdeG equations in an inhomogeneous system following a scattering approach 24 .

These wavefunctions are Nambu spinors of the form

ψ (η±) E↑ (x) = U (η±) E (x) V (η±) E (x) (A3)
where U E is the electron amplitude, V E the hole amplitude, and η denotes the electron/hole character of the quasiparticle state. Spin-down eigenfunctions can be easily obtained from spin-up ones by use of the electronhole symmetry in the system, through the substitution

U → V * , V → -U * .
For the sake of simplicity, we shall ommit the spin subindex in the wavefunctions that we discuss next.

Condensing all the phase difference δ in the right lead, the wavefunctions take the following shape

ψ (η±) E (x) = ψ (η±) Src (x)+ A (η±) ψ (e-) E,0 (x) + B (η±) ψ (h-) E,0 (x) Θ (-x)+ C (η±) ψ (e+) E,δ (x) + D (η±) ψ (h+) E,δ (x) Θ (x) , (A4) where ψ (η±)
Src is a source term of an quasielectron or quasihole impinging the contact from any of the leads

ψ (e+) Src = ψ (e+) E,0 (x) Θ (-x) , ψ (h+) Src = ψ (h+) E,0 (x) Θ (-x) , ψ (e-) Src = ψ (e-) E,δ (x) Θ (x) , ψ (h-) Src = ψ (h-) E,δ (x) Θ (x) .
(A5) The rest of the contributions to the wavefunctions are outgoing partial waves (as illustrated in Fig. 7). Imposing continuity for the wavefunction and its derivative (taking into account the effect of the delta barrier) the values of the partial wave coefficients are obtained for each incidence

A (e+) = -µ τ τ (µ τ + i) sinh 2 θ E • Q -1 B (e+) = iτ sin δ/2 sinh θ E + i δ 2 • Q -1 C (e+) = -iτ (µ τ + i) e -iδ/2 sinh θ E sinh θ E + i δ 2 • Q -1 D (e+) = -µ τ τ e -iδ/2 sinh θ E sin δ 2 • Q -1 (A6) A (e-) = C (e+) (-δ) B (e-) = D (e+) (-δ) C (e-) = A (e+) (-δ)
D (e-) = B (e+) (-δ)

A (h+) = B * (e+) B (h+) = A * (e+) C (h+) = D * (e+) D (h+) = C * (e+) A (h-) = C (h+) (-δ) B (h-) = D (h+) (-δ) C (h-) = A (h+) (-δ) D (h-) = B (h+) (-δ) , ( A7
)
where sinh θ E = ( E ∆ ) 2 -1, µ τ = R τ , τ is the normal transmission probability from the potential barrier, R = 1-τ its normal reflection probability, and Q = sinh 2 θ E + τ sin 2 δ 2 .

Β"

A" D"

C"

1"

H δ(x) Andreev Bound State wavefunction: The wave ampli-tudes for states with E < ∆ exhibit subgap poles at

(Δ = 0) Δ Δe iδ ( a ) Β" A" D" C" H δ(x) (Δ = 0) Δ Δe iδ ( b )
E = E A (δ) = ∆ 1 -τ sin 2 δ
2 , which signals the existence of a bound state at such an energy.

The wavefunctions for the Andreev Bound States (ABS) may be obtained in a similar way than those for states lying at E > ∆ taking into account that the quasiparticle momentum gains an imaginary component below the gap. The partial waves for the ABS are

ψ (e±) E A , φ (x) = 1 √ ξ0 u A v A e i φ e ±i(k F +iκ A )x ψ (h±) E A , φ (x) = 1 √ ξ0 v A u A e i φ e ∓i(k F -iκ A )x (A8)
These differ from the propagating partial waves because the u E , v E coefficients and the quasiparticle momentum

E,η±,σ Δ E A (φ) γ E,η±,σ + Φ 0 = Even; E, η±, σ γ Aσ + Φ 0 = Odd, σ γ A↑ + γ A↓ + Φ 0 = Even * Φ 0 = Even FIG. 8:
Single-particle excitations of the unperturbed SC Hamiltonian, with the parity notation for different ABS occupancies. The ground state of the system is noted as the even state. As single-particle excitations of the ABS have an odd parity, the ABS doublet excitation has an even parity, so it is branded as the even excited |Even * state.

κ E become complex for E < ∆:

u E → u A = 1 √ 2 e iθ A /2 v E → v A = 1 √ 2 e -iθ A /2 (A9) κ E → iκ A = iξ -1 0 sin θ A (A10) with sin θ A = √ τ sin δ 2 = ∆ 2 -E 2 A ∆ . ( A11
)
Only partial waves confined within a length κ -1 A , which diverges for E A → ∆, may appear in the wavefunctions

ψ A (x) = A A ψ (e-) E A ,0 (x) + B A ψ (h-) E A ,0 (x) Θ (-x) + C A ψ (e+) E A ,δ (x) + D A ψ (h+) E A ,δ (x) Θ (x) . (A12)
A linear homogeneous system of equations is obtained for the partial wave weights by applying the same conditions as in the case of the continuum states. The system exhibits a nontrivial solution for E = E A (δ). Eliminating the redundant equation and imposing the normalisation condition for the wavefunction, it is finally obtained that, in the bound states

                       A A = -i ÑA sin θ A -δ 2 e -i βτ e iδ/2 B A = -i √ R ÑA sin δ 2 • e i βτ e iδ/2 C A = √ R ÑA sin δ 2 e -i βτ D A = σ δ ÑA sin θ A -δ 2 e i βτ (A13)
where sin βτ = √ τ , σ δ = sign (δ) and

ÑA = -σ δ √ τ 2 cos θ A sin θ A -δ 2 .
(A14)

A similar derivation of this result can be found in 19.

The solutions of the BdeG equations allow us to express the electron field operators as

Ψ σ (x) = U A (x) γ Aσ -σV * A (x) γ † Aσ + E,η± U (η±) E (x) γ E,η±,σ -σV (η±) * E (x) γ † E,η∓,σ (A15)
where γ Aσ and γ E,η±,σ are the quasiparticle operators which diagonalise the SC Hamiltonian. The excitation spectrum of the system is represented in Fig. 8.

Appendix B: Current operator and transition rates

The current operator in the new basis defined by Eq. (A15) is

Î (x) = - e 2mi i,j,σ γ † iσ γ iσ   U * i dUj dx - dU * i dx U j σ U * i dV * -j dx - dU * i dx V * -j σ V -i dUj dx -dV-i dx U j V -i dV * -j dx -dV-i dx V * -j   γ jσ γ † j σ (B1)
The sum in i, j indices are a shorthand notation for all the different contributions appearing in Eq. (A15). The minus sign in front of some particular wavefunction subindices notes that such a component corresponds to the antiparall wavefunction (e.g., if i = E (e+), then

U -i = U (e-) E
), which only applies when the index corresponds to an excitation in the continuum.

Transitions involving the odd states

A relevant matrix element in the problem is the one associated to the |Odd, σ → |Even; E, e+, σ process, which is found to be

- ek F m √ Lξ 0 C A C * (e+) -D A D * (e+) (u E u A + v E v A ) . (B2)
The electron-hole symmetry in the field transformations ensures that this matrix element is the complex conjugate of the matrix element associated to the process |Odd, σ → |Even; E, h-, σ .

In the ballistic limit τ → 1, restricting δ to the interval [0, π] so as to establish a bijection between δ and E A (δ), one finds

|C A | → ∆ 2 -E 2 A ∆ 2 1/4 D A → 0 (B3) C (e+) → E 2 -∆ 2 E 2 -E 2 A D (e+) → 0. (B4)
The squared amplitude of these matrix elements in this limit is

2 e 2 k 2 F m 2 Lξ 0 ∆ 2 -E 2 A ∆ E 2 -∆ 2 E 2 -E 2 A 1 + E A E . ( B5
)
Matrix elements for the |Odd, σ → |Even; E, e-, σ and |Odd, σ → |Even; E, h+, σ processes vanish in the limit of perfect transmission. Conversely, the squared amplitude for the quasiparticle recombination processes |Odd, σ; E, e-, σ → |Even and |Odd, σ; E, h+, σ → |Even in the limit τ → 1 is found to be

2 e 2 k 2 F m 2 Lξ 0 ∆ 2 -E 2 A ∆ E 2 -∆ 2 E 2 -E 2 A 1 - E A E . ( B6
)
Whereas the amplitudes for the other two recombination processes, which are |Odd, σ; E, e+, σ → |Even and |Odd, σ; E, h-, σ → |Even , are zero in this same limit. The products of the terms that contain in these expressions the functional dependence in E and E A with the superconducting density of states yield the factors g (E, E A ) mentioned in the main article

g (E, E A ) = (E 2 -∆ 2 ) (∆ 2 -E 2 A ) ∆ (E -E A ) . ( B7
)
In the opposite tunnel limit τ → 0 all squared amplitudes tend to zero as τ 3/2 , with the leading term being the same for the four different processes

2 e 2 k 2 F m 2 Lξ 0 τ 2 ∆ 2 -E 2 A ∆ 1 + ∆ cos δ E . ( B8
)
The same applies to amplitudes of recombination processes, with a minus sign appearing inside the parenthesis instead of a plus.

From these we may define another g (E, E A ) factor for the tunnel regime.

g (E, E A ) = τ 2 ∆ 2 -E 2 A E 2 -∆ 2 E ∆ + sgn (E A ) cos δ .
(B9) These results are in agreement with the ones recently derived by Kos et al. using a different method 25 taking also into account the factor (E ± E A )

-1 that comes from the environmental density of states (see Eq. (C1)).

Transitions between the even states

An analytical expression for the amplitude that links the two even states can be derived for any value of τ

Even Î Even * = ek F mξ 0 C A D A u 2 A -v 2 A e iδ (B10) Even Î Even * 2 = e 2 ∆ 4 2 (1 -τ ) τ 2 sin 4 δ 2 E 2 A (δ)
. (B11)

These results coincide with the results from 12.

corresponding electromagnetic ("plasma") mode is modeled with a resistance R J in series with the capacitor. The SQUID loop formed by the contact and the junction is connected to a biasing circuit through an inductor L e and a capacitance C e . Dissipation in this circuit is modelled by a resistance R e in series with the capacitor (see Fig. 9(a)). The total impedance seen by the SC is

Z env (ν) R Q = 4 π E C hν p0 i ν νp0 a e (ν) b p0 (ν) a p0 (ν) a e (ν) -ν 2 ν 2 3 b p0 (ν) , (C2) with a (ν) = 1 + i 1 Q ν ν -ν 2 ν 2
where ≡ e, p0 and b p0 (ν) = 1 + i 1 Qp0 ν νp0 . The characteristic frequencies and quality factors inside these expressions are

ν p0 = 1 2π (L J C J ) -1/2 ν e = 1 2π (L e C e ) -1/2 ν 3 = 1 2π (L J C e ) -1/2 (C3) Q p0 = 1 R J L J C J Q e = 1 R e L e C e . ( C4
)
The equivalent circuit in Fig. 9 has two modes. A lowfrequency mode determined essentially by the on-chip LC filter connecting the SQUID to the outside world, and a high-frequency mode corresponding to the plasma oscillation of the junction "dressed" by the external circuit:

ν P = 1 2π L -1 J + L -1 e C J (C5)
Parameters for this equivalent circuit were obtained in the following way:

The Josephson junction inductance L J = 595 pH is determined by the critical current extracted from the switching probability measurements. The other five parameters of the equivalent circuit of Fig. 9 were adjusted so as to reproduce at best all the available experimental information.

a. The energy gap measured from the IV characteristics is ∆ = 194 µeV.

b. The dressed plasma frequency must be close to ∆/2h in order to explain the position of the sharp threshold observed in the rates Γ in and Γ out . This is compatible with what was reported in 17 .

c. A value of C e = 60 pF is expected from a measurement at very low frequency (100 kHz) on a larger test capacitor fabricated on the same run.

d. The DC sub-gap current of the JJ alone is I J = 4.4 nA at V J = ∆/4e. Imposing the power equality I J V J = R J I 0 between the DC injected power and the microwave power absorbed at ∆/2h by the junction's environment, we get R J = 0.3 Ω.

e. The low frequency mode of the environment was measured at 558 MHz in a separate microwave reflectometry experiment [START_REF] Bretheau | [END_REF] .

f. At this resonance frequency, the reflection amplitude S 11 shows a dip of -15 dB, from which we determine R e = 0.25 Ω.

The two capacitances, C J and C e , and the environmental inductance L e were adjusted so as to reproduce the two characteristic frequencies of the circuit. The chosen value C J = 168 fF is a 25% lower than what is expected from the nominal area of the junction (2.8 µm 2 ) and the typical specific capacitance for the junctions fabricated usually in our laboratory (75 fF/µm 2 ). C e = 68 pF is a 13% higher than what is expected from the test low frequency measurement. Finally, the nominal value L e = 600 pH for the environmental inductance is an 80% of what is expected from a crude geometrical estimation.

With these values we predict a dressed plasma mode frequency of 0.48∆/h and a quality factor Q = 116.

Appendix D: Relaxation due to phonons

The electron-phonon interaction in real space is [START_REF] Fetter | Quantum Theory of Many-Particle Systems[END_REF] Ĥe

-ph = γ dr σ Ψ † σ (r) Ψ σ (r) φ (r) (D1)
where φ (r) is the phonon field operator: φ (r) = q hν q 2V b q e iqr + b † q e -iqr . (D2)

The electron-phonon coupling constant γ is

Z 2 π 2 mk F n0 B 1 2
, with n 0 being the atomic density, B the adiabatic bulk modulus and Z the electron valence from the superconductor.

The SC density inside the interaction Hamiltonian takes a form similar to the current operator's in the quasiparticle basis (Eq. (B1)). But differently from the coupling with the EM environment, the phonon coupling depends on the geometrical spread of the SC wavefunctions, a feature characteristic of the coupling of phonons with localized states 23,[START_REF] Ivanov | [END_REF] .

After eliminating terms linear in e ±2ik F x , which vanish in the spatial integration due to their rapid oscillatory behavior, the matrix element associated to the process |Odd, σ → |Even; E, e+, σ is found to be

C * (e+) C A -D * (e+) D A (u E u A -v E v A )
e -(κ A +iκ E )x √ Lξ 0 .

In the limit τ → 1, the squared amplitude of the part of this matrix element that does not depend on x, as well as the similar quantity obtained from the matrix element for the process |Odd, σ → |Even; E, h-, σ , tend to the expression in Eq. (B6), except for the factor ( ek F /m) 2 . It can also be found in the same limit that the analogous x-independent quantity for the |Odd, σ; E, e-, σ → |Even and |Odd, σ; E, h+, σ → |Even processes is, on the other hand, the same as in Eq. (B5).

The spatial integrals are of the form ∞ 0 dx • e -(κ A ±iκ E )x sin q x x • F (q ⊥ , x) ,

where

F (q ⊥ , x) = A ⊥ (x) d 2 r ⊥ e iq ⊥ r ⊥ |Ψ ⊥ (r ⊥ )| 2 ,
Ψ ⊥ (r ⊥ ) being the axial spread of the SC wavefunctions on the leads, whose geometric details are enclosed in their cross section A ⊥ (x).

The momentum transfer Q to or from a phonon taking an active role in these relaxation processes is large compared to the inverse penetration length of the ABS: Q κ A , κ E in the region E A (δ) < ∆ 2 . Following the approximations detailed in 23 in the theoretical description of the phonon-mediated |Even * → |Even relaxation, the F (q ⊥ , x) factor introduces a cutoff L in the integral in the x direction

F (q ⊥ , x) → θ L -|x| (D3)
and the spatial integration may be easily evaluated in the limit L-1

Q κ A , κ E : 4Q sin Q L Q 2 + |κ A ± iκ E | 2 ≈ 4 L.
Combining the different contributions, the squared amplitude for the process |Odd, σ → |Even; E, e/h±, σ mediated by the emission of a phonon is

|M 1 | 2 = 8 hν Q V L2 γ2 Lξ 0 1 - E A ∆ 2 E 2 -∆ 2 E 2 -E 2 A E -E A E .
The transition rate for such a process is

2π V d 3 Q (2π) 3 4 |M 1 | 2 ρ SC (E) (1 -f FD (E, T qp )) f BE (hν Q , T ph ) δ (hν Q -E + E A ) , (D4)
where we have used the same notation as in the main text. We may rewrite the integral over momenta in Eq. (D4) as an integral over energies, with a density of states quadratic in E -E A that appears as a result of this transformation. The resulting total rate is then Γ . If this factor is of the order of 10 -2 , the phonon-induced relaxation rates are reduced to around 1 kHz, which coincide with other estimations in analogous systems 14 .

FIG. 1 :

 1 FIG. 1: (a) Schematic representation of a single superconducting channel coupled to an impedance Zenv. δ and γ indicate the phase drops through the channel and the impedance, respectively, and Φ is the magnetic flux through the loop. (b) Quasiparticle excitation spectrum in the single channel weak link, with the continuum above the gap ∆ and a subgap discrete spin-degenerate Andreev level of energy EA.

  and the corresponding rates are denoted by Γ (a,b) out . The Fermi golden rule for the first process yields Γ

  as in Eqs.(2,3) but with a quadratic density of states and T env replaced by the phonon temperature T ph (see Appendix D for more details).

FIG. 2 :

 2 FIG. 2: Left panel: Schematic representation of processes removing a quasiparticle from the Andreev level. Process (a) involves absorption of an environmental photon or phonon and corresponds to the transition rate Γ (a) out . Process (b) corresponds to the recombination into a Cooper pair of the quasiparticle trapped in the Andreev level with a quasiparticle from the continuum, with the emission of a photon or a phonon, and is characterized by a transition rate Γ (b) out . The right panel shows the rates Γ (a) out (red) and Γ (b) out (blue) resulting from the absorption and emission of environmental photons (full lines) or phonons (dashed lines), for the parameters of Ref. 1 (see Appendix C). We have set kBTenv = 0.06∆, kBTqp = 0.09∆, kBT ph = 0.015∆.

  (a,b) out obtained using parameters which are appropriate for the experimental situation of Ref. 1 (see Appendix C) are shown by the solid lines in Fig. 2. For E A ∆ -hν p ∼ 0.52∆, Γ (a) out is large because photons in the plasma mode can excite the trapped quasiparticle out into the continuum. Similarly, Γ

FIG. 3 :FIG. 4 :

 34 FIG. 3: Left: Schematic representation of processes adding a quasiparticle in the Andreev level. Right: the emission processes (a) rates, denoted by Γ (a) in in the text, are given by the red curves for photons (full line) and phonons (dashed lines). The processes (b) involving the breaking of a Cooper pair (blue line in the plot) are much less efficient. The one involving the phonons is below this scale. Same parameters as in Fig. 2.

FIG. 5 :

 5 FIG. 5: Total rates Γout (upper curves and symbols) and Γin (lower cuves and symbols) calculated by the present model and obtained experimentally in Ref. 1 for different contact transmissions. The same convention as in Fig. 4 is used.

FIG. 6 :

 6 FIG. 6: Quasiparticle energies and two-component eigenfunctions near the Fermi level for spin-up quasiparticles. There is a four-fold degeneracy of states for each energy E, with two quasielectron states (blue dots) with a momentum |k| = kF +κE and another two quasiholes with an absolute momentum of |k| = kF -κE.

FIG. 7 :

 7 FIG.7: Schematic representation of the scattering problem for the case of a quasi-electron impinging the contact and the resulting outgoing quasiparticle partial waves (a), either normally or Andreev-reflected or transmitted. The scattering problem for the ABS does not require a source term (b), but these outgoing partial waves exhibit an evanescent behaviour.

FIG. 9 :

 9 FIG. 9: (a) Electrical elements in the relevant neighborhood of the SQUID loop. (b) Equivalent electrical model for the Josephson junction in the SQUID. (c) Environmental impedance as seen by the SC.

3 g 3 gγ2 π 2 ∆ cs 3 is 8 π L ξ0 2

 332382 (E, -E A ) f BE (E -E A , T ph ) (1 -f FD (E, T qp )) . (D5)Repeating the same process for the process for the form |Odd, σ; E, ησ → |Even yieldsΓ (E, E A ) (1 -f BE (E + E A , T ph )) f FD (E, T qp ) . (D6)Using the material constants for aluminum, the factor 2π of the order of 10 GHz. The phonon rate is reduced because of the relevant thermal factors (which are of the order of e -β∆ ∼ 10 -5 ) and the geometrical factor

Acknowledgments

The authors acknowledge fruitful discussions with J. Martinis, R. Egger, A. Zazunov, D. Urban, J. C. Cuevas, F. S. Bergeret and A. Martín-Rodero. Financial support by EU FP7 SE2ND project, Spanish Mineco project FIS2011-26516, ANR contracts DOCFLUC and MASH, C'Nano and by the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme (FP7/2007-2013) under REA grant agreement no. PIIF-GA-2011-298415 is acknowledged.

Appendix C: Modelling the EM environment for the experiment in Ref. [1] The density of environmental modes D (hν) is, following the formalism presented in Ref. 20,

where Z env (ν) is the electric impedance as seen from the SC, represented in Fig. 9(c).

In the experiment shown in Ref. 1, a SC is placed in parallel with a Josephson junction having a critical current much larger than that of the atomic contact (see Fig. 9(b)). This junction is perceived by the atomic contact as the parallel combination of an inductor L J = ϕ 0 /I 0 and a capacitor C J , I 0 being the critical current of the junction. The finite quality factor of the


J.-P Gras

J.-Y Delenne

M S El

Study of capillary interaction between two grains: a new experimental device with suction control

Keywords: Capillarity, Liquid bridge, Tensile test, Suction, Experimental device, Rupture criterion J.-P. Gras (B)

We investigated the behavior of a water liquid bridge formed between two grains. We mainly focused on tensile tests with suction control (capillary pressure). Theo-retical and experimental studies are compared. A new exper-imental device involving suction control of the liquid bridge was developed specifically for this kind of test. Most of the liquid bridge variables and characteristics were measured by image analysis (gorge radius, volume, contact angles, fill-ing angles). Capillary force was measured by differential weighting. Experimental conditions allows us to avoid vis-cous effects. Our experimental results were close to Young-Laplace equation solutions. The "gorge method", commonly used for calculating the capillary force, was also validated by our experiments. Liquid bridge rupture was studied and a new rupture criterion is proposed. This criterion depends on the grain radius, contact angle, surface tension and suction and was in agreement with the experimental results.

Introduction

Understanding the role of capillary interactions in the overall behaviour of a wet granular material is of considerable technological importance for applications in many fields, such as civil engineering, pharmaceutical research, agronomy, etc. Some properties, such as tensile strength and water reten-tion, strongly depend on the capillary interactions between grains. For example, sand castles may be sculpted when sand is mixed with water [START_REF] Schiffer | Granular physics: a bridge to sandpile stability[END_REF]. Liquid bridges between grains generate attractive forces when the water content is low enough. These forces are due to the surface tension γ and the difference between the air pressure p a and the liquid pressure p w across the liquid/air interface, also called suction. Know-ing the capillary force that generates such liquid bridges is important to determine the local origins of the mechanical behaviour of unsaturated granular media [START_REF] Richefeu | Force transmission in dry and wet granular media[END_REF].

A stable liquid bridge can be obtained by controlling the suction or liquid volume. We decided to control the suc-tion in this study. One advantage is that the suction of every liquid bridge in a sample is equal to the macroscopic suction imposed upon a laboratory sample (at equilibrium). We put forth the hypothesis that the equilibrium of suction and liquid bridge volumes is achieved by the formation of a thin film of water at the grain surface, thus enabling a quite rapid return to equilibrium. For example, after mixing water with a gran-ular medium, the relaxation time of liquid bridge volumes is around 5 min [START_REF] Scheel | Liquid distribution and cohesion in wet granular assemblies beyond the capillary bridge regime[END_REF].

Several theoretical studies on the behaviour of liquid bridges between two spherical grains have been carried out using numerical solution of the Young-Laplace equation or a method based on the approximation of liquid bridge pro-files via circles (toroidal approximation). In most cases, these studies involve tensile tests on liquid bridges with an imposed volume. Some expressions which link the capillary force, volume and geometrical and physical characteristics of the liquid bridge have been proposed for monodisperse [START_REF] Mikami | Numerical simulation of cohesive powder behavior in a fluidized bed[END_REF][START_REF] Pitois | Liquid bridge between two moving spheres: an experimental study of viscosity effects[END_REF] and polydisperse [START_REF] Soulié | Capillary cohesion and mechanical strength of polydisperse granular materials[END_REF][START_REF] Richefeu | A model of capillary cohesion for numerical simulations of 3d polydisperse granular media[END_REF] cases. Recently, some authors have stud-ied the relation between suction and liquid bridge volume [START_REF] Molenkamp | Interactions between two rough spheres, water bridge and water vapour[END_REF][START_REF] Likos | Hysteresis of capillary stress in unsaturated granular soil[END_REF][START_REF] Lechman | Capillary force and water retention between two uneven-sized particles[END_REF][START_REF] El Shamy | Micromechanical aspects of the shear strength of wet granular soils[END_REF] for given interparticle distances. Monodisperse [START_REF] Molenkamp | Interactions between two rough spheres, water bridge and water vapour[END_REF][START_REF] Likos | Hysteresis of capillary stress in unsaturated granular soil[END_REF][START_REF] El Shamy | Micromechanical aspects of the shear strength of wet granular soils[END_REF] or polydisperse [START_REF] Lechman | Capillary force and water retention between two uneven-sized particles[END_REF] grain-pairs were considered in these studies. The toroidal approximation has been validated in monodisperse and polydisperse cases for given interparticle distances [START_REF] Molenkamp | Interactions between two rough spheres, water bridge and water vapour[END_REF][START_REF] Lechman | Capillary force and water retention between two uneven-sized particles[END_REF]. Recently, tensile tests were briefly investigated with suction control [START_REF] Scholtès | Modélisation Micro-Mécanique des Milieux Granulaires Partiellement Saturés[END_REF], with a numerical resolution of the Laplace-Young equation.

Experiments have been carried out to study tensile tests on liquid bridges for a given volume. One of the first experimental studies was conducted by Mason and Clark [START_REF] Mason | Liquid bridges between spheres[END_REF], who measured the capillary force as a function of the interparticle distance between two monodisperse spheres or between a sphere and a plate. Since this experiment, several experimental studies have been carried out on polydisperse and monodisperse spheres [START_REF] Pitois | Liquid bridge between two moving spheres: an experimental study of viscosity effects[END_REF][START_REF] Hotta | The capillary binding force of a liquid bridge[END_REF][START_REF] Soulié | Cohésion par capillarité et comportement mécanique de milieux granulaires[END_REF][START_REF] Willet | Capillary bridges between two spherical bodies[END_REF][START_REF] Willet | Effects of wetting hysteresis on pendular liquid bridges between rigid spheres[END_REF], with close agreement with Young-Laplace equation solutions. We should underline the experiments of Willet et al. [START_REF] Willet | Capillary bridges between two spherical bodies[END_REF] which seem very close to modelling results because of the use of spheres with an exceptional surface state. To our knowledge, however, no experimental studies have focused on liquid bridges with suction control.

In this paper, we present an original experimental device which allowed us to conduct tensile tests on liquid bridges with suction control. The results of this experiment are compared to Young-Laplace equation solutions, and the "gorge method" is tested to calculate the capillary force. Then, we present some results concerning the rupture distance.

Experimental setup

An experimental device was designed to perform displacement controlled tensile tests on capillary bridges between grain-pairs for a given suction (Fig. 1). A "grain-pair" is an assembly of two spherical grains linked by a liquid bridge. This device allows the measurement of capillary force as a function of displacement, suction control, and measurement of the main liquid bridge characteristics by image analysis. The capillary force f cap is calculated by differential weighting :

Monitoring of tensile force as a function of displacement

f cap = (m d -m a )g (1)
where m d and m a are respectively the dry mass and the apparent mass (in the presence of a liquid bridge) of the lower bead

The experimental device is shown in Figs. 1 and Distilled water is introduced between two beads of 10 mm diameter via a channel in the upper bead connected to the water supply. The lower bead is attached to a high precision balance (10 -3 g) which monitors variations in mass during the tensile test. The output of the balance is connected to a Fig. 3 Details on the application of suction, which is deduced from the i measure and its support. The bridge volume in our experiment was estimated to be in the 0.1μl -4μl range. It corresponds to a weight in the 10 -6 N -0.4 * 10 -4 N ) range, which is relatively weak in comparison to the capillary force (10 -4 N -10 -3 N ) but could be substantial in rare cases. Its effect contributes to decreasing the measured capillary force. The force that we measure would thus be an underestimation of the real capillary force. The liquid bridge volume and its weight was determined (using image analysis) and taken into account to compare the theoretical to the experimental results.

Suction control

The air pressure is the atmospheric pressure. The water pressure in the bridge is given by the difference h i between the level of the liquid bridge neck and the level of the water in the water tank (Fig. 3). A hole in the upper bead ensures continuity between the liquid bridge water and the water of the water tank. The height of the water in the water tank is regulated by a laboratory scissor jack and measured by a graduated burette.

The Young-Laplace equation can be applied for the liquid bridge:

p a -p w = γ 1 ρ ext - 1 ρ int (2)
where ρ ext and ρ int are respectively external and internal absolute values of radii of curvatures of the liquid bridge surface. At the level of the liquid bridge gorge, the internal radius of curvature is called the gorge radius h.

In the water tank and buret, water and air pressures are equal at the air-water interface. Suction s is thus defined by: Suction was not completely constant during the tensile test. As the upper bead rose, h i changed during the test. Nevertheless, these changes could be easily evaluated and in most cases they were too weak in comparison to the initially imposed suction value because of the low rupture distance D rupt of the liquid bridge.

s = p a -p w = ρgh i (3)
In order to avoid air bubbles which may alter the water pressure control, a drain valve allowed the elimination of air bubbles (Fig. 2). A temperature of 23 • C was kept throughout the experiments.

Image analysis

For given geometrical characteristics (grain radius, interparticle distance) of the grain-pair (Fig. 4), and given physical characteristics (contact angles, liquid/air surface tensions), we examined the effect of suction on the liquid bridge behaviour. The image analysis allowed us to measure the main liquid bridge characteristics.

A digital camera with macrozoom allowed us to accurately record the shape of the liquid bridge during the tensile test at the rate of one picture every 5 sec. It also allowed us to control the alignment of the two beads.

An image analysis program was developed to extract the main liquid bridge characteristics from the photos. Each bead boundary was fitted by a circle which enabled us to determine the scale of the photos (Fig. 5). First, the left and right profiles y(x) of the liquid bridge were detected and estimated by a fourth order polynomial. It was then possible to localize the contact point between the liquid bridge and the beads. Then the main geometrical characteristics of the liquid bridge were evaluated (Fig. 6):

-Gorge radius. We localized the minimum A and B of the right and left profiles (Fig. 5). The gorge radius is half of the distance between these minimums. -Filling angles. These are determined according to the centers of the beads and the contact points. -Contact angles. Contact angles are the angles between bead boundaries and polynomial representations of the liquid bridge (Fig. 6). -Volume. The volume is obtained by calculating the surface covered by the liquid bridge while assuming that it is a surface of revolution (see Fig. 7):

V = π x c2 y 2 (x)dx -V 1 -V 2 (4)

Experimental results

The tensile tests were carried out with suction control. The granular material used for the experiments were monodisperse beads of 10 mm grain diameter. Figure 8 shows the variations in the liquid bridge during three steps of the tensile test. We noted that the liquid bridge formed was quite symmetrical during the test.

Force and volume as a function of the interparticle distance at constant suction

Several tests were conducted for different suction values (Fig. 9). The results were highly repeatable (very low standard deviation). The behaviour was very close to the solutions of the Young-Laplace equation:

-In the case of tensile tests with suction control, the function representing the relation between the capillary force and the interparticle distance was monotically decreasing and convex (Fig. 9) instead of concave as reported in tensile tests with volume control [START_REF] Mason | Liquid bridges between spheres[END_REF][START_REF] Willet | Capillary bridges between two spherical bodies[END_REF][START_REF] Soulié | Influence of liquid bridges on the mechanical behaviour of polydisperse granular materials[END_REF]]. -The rupture distance decreased with increased suction (Fig. 9). -For weak interparticle distances, the capillary force increased with increased suction (Fig. 9). -For large interparticle distances, the capillary force decreased with increased suction (Fig. 9).

The influence of suction on variations in the capillary force as a function of the interparticle distance in tensile tests with suction control was found to be the opposite of the influence of volume in tensile tests with volume control. No peaks were noted in tensile tests with suction control, as generally observed in tensile tests with volume control [START_REF] Mason | Liquid bridges between spheres[END_REF][START_REF] Soulié | Influence of liquid bridges on the mechanical behaviour of polydisperse granular materials[END_REF]. Such peaks have been previously attributed to changes in contact

x c1

Experimental protocol

The first step involves filling the hydraulic system with water. For this, the hole of the upper bead is stopped up. Then water is added via the water tank until water flows off the drain valve. Then the drain valve is closed and the hole of the upper bead is opened. Water flows through the hole until h i = 0. Then a liquid bridge is formed for h i = 0 by positioning the lower bead. By moving the water tank vertically and setting it at a given height, suction is applied to the liquid bridge. The tensile test can begin by applying a vertical velocity to the upper bead. [START_REF] Soulié | Influence of liquid bridges on the mechanical behaviour of polydisperse granular materials[END_REF]. Moreover, we noted that variations in volume during the tensile test were not monotonic (Fig. 10). The volume first increased with the interparticle distance and finally decreased until rupture of the liquid bridge.

Experimental validation of the gorge method

The "gorge method" is conventionally used to calculate the capillary force. In this method, the capillary force f cap is defined as a function of the gorge radius h, suction s and surface tension γ , and is expressed as:

f cap = 2πγ h + π sh 2 (5)
Surface tension is first measured at the outlet of the hydraulic system by wrenching (γ = 0.051 N /m), suction is applied, and the gorge radius can be measured by image analysis. The low water surface tension value is attributed to impurities in the hydraulic system.

Figure 11 presents variations in the gorge radius as a function of the interparticle distance. Experimental gorge radius, water surface tension and suction values were used to determine the capillary force using the "gorge method". Capillary forces determined by the "gorge method" were very close to the experimental capillary forces (Fig. 12). For the comparison, we took into account the weight of the liquid bridge, which was calculated from its volume V measured by image analysis. The "gorge method" appears to be a good approximation to evaluate capillary forces for a symmetrical liquid bridge.

Force displacement and rupture modelling

Experimental capillary forces were compared to the solutions of the Young-Laplace equation. The "grain-pair" characteristics are known: grain radius (R = 5 mm), surface tension (γ = 0.051 N /m) and variations in the contact angle θ with the interparticle distance D [START_REF] Gras | Approche micromécanique de la capillarité dans les milieux granulaires: rétention d'eau et comportement mécanique[END_REF]. Some effects are due to gravity: upper top and lower bottom contact to predict the rupture distance of the liquid bridge. The Lian criterion of rupture is expressed as:

D rupt = (1 + 0.5θ) V 1 3 (6)
Firstly, the validity of this criterion in tensile tests with volume control was tested experimentally with our device. For this particular test, the upper bead did not have a hole and the volume was introduced between the two beads using a microsyringe. In this case, the Lian criterion gave a good prediction of the rupture distance as a function of the volume (Fig. 14): experimental data (square) are close to the Lian criterion. Secondly, the validity of this expression in tensile tests with suction control was tested experimentally. In this case, we noted as expected that the experimental data (triangles) were not in agreement with the Lian criterion but fit well the estimation of rupture from Young-Laplace equation (Fig. 14). These results clearly confirm that the rupture distance depends on the loading path of the liquid bridge.

D (mm)

Fig. [START_REF] Scholtès | Modélisation Micro-Mécanique des Milieux Granulaires Partiellement Saturés[END_REF] Comparison of experimental results obtained with the "gorge method" (Eq. 5) for different suctions in the case of tensile tests with suction control angles do not have the same value. We noted that the difference between both decreased with increased suction. In order to simulate the experimental behaviour, experimental variations in the upper and lower contact angles were fitted and then taken into account in the model. Then, by using the grain-pair characteristics during the tensile test, it was possible to compare the Young-Laplace solutions, associated with the "Gorge method", to the experimental data in terms of capillary forces or rupture distances [START_REF] Gras | DEM and experimental analysis of the water retention curve in polydisperse granular media[END_REF]. We noted that the experimental and theoretical capillary forces were very close (Fig. 13). The approximation was better for high suctions. Otherwise, we noted that the gap between the theoretical approximation and the experimental data increased as a function of the interparticle distance. This increase could be explained by the gravity effect, which is greater for high interparticle distances (Appendix 1).

In tensile tests with controlled volume, an empirical expression proposed by Lian et al. [START_REF] Lian | A theoretical study of the liquid bridge forces between two rigid spherical bodies[END_REF], is generally used A rupture criterion would be necessary for tensile tests with suction control. In this part, the results are presented with dimensionless variables. For grains with different radii R 1 and R 2 , we used the arithmetic mean radius R 12 = 2R 1 R 2 R 1 +R 2 , the interparticle dimensionless distance D * = D R 12 and the dimensionless suction H 12 = s R 12 γ . We noted that the relation between the dimensionless rupture distance and the dimensionless suction hardly depended on the polydispersity ratio r = R 1 R 2 for 0, 1 ≤ r ≤ 1 in tensile tests with suction control (Fig. 15a) for null contact angles. We only plotted curves for r = 0.1 and r = 1, but the other curves for 0, 1 < r < 1 are included in the space between the two plotted curves. From this observation, the dimensionless rupture distance for a null contact angle D * rupt θ 0 can be approximated by:

D * rupt θ 0 = (H 12 + 3) -0.93 (7)
On Fig. 15b, we note that the relation between the ratio

D * rupt D *
rupt θ 0 and the contact angle θ hardly depends on the polydispersity ratio r for 0.1 ≤ r ≤ 1 and on the dimensionless suction H 12 and follows a cosine curve. From this observation, we formulated the following equation:

D * rupt D * rupt θ 0 = (cos θ ) 1.7 (8)
Finally, it is possible to establish a rupture criterion of the bridge which gives the dimensionless rupture of the liquid bridge as a function of the dimensionless suction H 12 and contact angle θ from Eqs. 7 and 8. Its expression is:

D * rupt = (H 12 + 3) -0.93 * (cos θ) 1.7 (9)
This rupture criterion fit well the estimation of rupture from Young-Laplace equation. We noted that, for a mean contact angle of 20, as measured by image analysis, the rupture criterion was very close to the experimental data for rupture (Fig. 16).

Conclusion

An innovative experimental device has been developed for tensile tests at constant suction and allows measurement of the main liquid bridge characteristics (capillary force, interparticle distance, liquid bridge volume, gorge radius, contact angles, filling angles) with excellent repeatability. This device also enables validation of an expression commonly used to calculate the capillary force, i.e. the "gorge method". A model based on the Young-Laplace equation coupled with the gorge method for calculation of the capillary force is proposed. This model specifically takes into account differences in contact angle between the liquid and two different grains. Our comparison between Young-Laplace solutions and experimental results concerning capillary forces indicated very close values. The validity of the Lian criterion for rupture did not seem to agree in the case of tensile tests with suction control. A criterion which depends on the dimensionless suction and contact angle is proposed for tensile tests with suction control. This result is very useful for studying the effect of capillarity on the overall cohesion of a granular material because the density of liquid bridges is directly linked to the rupture distance of liquid bridges.

Appendix: Effect of gravity

We consider a vertical liquid bridge under gravity. The effect of gravity decreases the rupture distance, creates distortion in the liquid bridge and changes the capillary forces. Under gravity, a difference of pressure p across the liquid/air interface is: p = p 0 -(ρ w -ρ a )gx [START_REF] Rotemberg | Determination of surface tension and contact angle from the shapes of axisymmetric fluid interfaces[END_REF][START_REF] Zeinali | Capillary holdup between vertical spheres[END_REF]. ρ w et ρ a are respectively densities of water and air, g is the acceleration due to gravity, p 0 is the pressure difference between water and air at the liquid bridge gorge and x is the vertical position from the gorge section. Then, for a vertical liquid bridge with x defining the vertical axis, consideration of the gravity on the liquid bridge shape modifies the Young-Laplace differential equation, which becomes:

-p 0 γ + (ρ w -ρ a )g γ x = ÿ (x) 1 + ẏ2 (x) 3 2 - 1 y(x) 1 + ẏ2 (x) (10)
From this equation, we note that the effect of gravity increases with the interparticle distance [START_REF] Mazzone | The effect of gravity on the shape and strength of a liquid bridge between two spheres[END_REF]. The dimensionless equation becomes: . In the experimental studies presented, the maximum ratio between the two left terms of the former equation was noted when x was high. Hence, the variation in water pressure inside the bridge increases as the interparticle distance increases. For example, in the above experiments, the maximum value of this ratio (at rupture) was respectively under 2, 7 and 15 % for suctions of 140 Pa, 70 Pa and 40 Pa. This means that the variation in water pressure inside the bridge due to gravity increases when suction decreases. This explains why the model was more accurate at high suction values. Considering that errors due to gravity are not very marked, we do not have to take gravity into account in our model. Nevertheless, its effects are of interest when analyzing the results.

Fig. 1 Fig. 2

 12 Fig. 1 Photo of the experimental device

Fig. 4

 4 Fig.[START_REF] Mikami | Numerical simulation of cohesive powder behavior in a fluidized bed[END_REF] Geometrical model of a liquid bridge between two grains of different sizes: h is the gorge radius, D is the interparticle distance, δ 1 and δ 2 are the filling angles, θ 1 and θ 2 are the contact angles, R 1 and R 2 are the grain radii with R 1 ≤ R 2

Fig. 5 Fig. 6

 56 Fig. 5 Detection of the beads and liquid bridge contours by image analysis

Fig. 7

 7 Fig. 7 Details of the liquid bridge

Fig. 8

 8 Fig. 8 Photos of the liquid bridge. From the top to bottom photos, the interparticle distance is 0.0 mm, 0.3 mm and 0.67 mm (near rupture), respectively; R = 5 mm and s ≈ 40 Pa

Fig. 9 Fig. 10

 910 Fig. 9 Variations in capillary force as a function of the interparticle distance in the case of tensile tests with suction control. Each point represents the average of 6 points

Fig. 11 Fig. 13 Fig. 14

 111314 Fig. 11 Experimental variations in the gorge radius measured by image analysis as a function of the interparticle distance and for different suctions in the case of tensile tests with suction control

Fig. 15

 15 Fig. 15 Tensile tests with suction control. a Variations in D * rupt as a function of H 12 for θ = 0. b Variations in the ratio D * rupt D * rupt θ0 for several r and H 12 values as a function of the contact angle θ

Fig. 16

 16 Fig. 16 Rupture criterion for tensile tests with suction control (θ = 20 • )

H = H 0 12 + 2 - 1 y

 1221 B 0 x * = ÿ * (x * ) 1 + ẏ * 2 (x * ) 3 * (x * ) 1 + ẏ * 2 (x * ) (11) with y * (x) = y(x) R 12 , x * = x R 12 , H 0 12 = -p 0 R 12 γ and B 0 = (ρ w -ρ a )g R 2 12 γ