Vincent Jacques

E Wu

Frédéric Grosshans

François Treussart

Philippe Grangier

Alain Aspect

Jean-François Roch

Delayed-choice test of complementarity with single photons

Keywords:

We report an experimental test of complementarity using clock-triggered single-photon pulses emitted by an individual N-V color center in a diamond nanocrystal. The single photons are sent into a Mach-Zehnder interferometer with an output beamsplitter of adjustable reflection coefficient R. In addition, the choice of introducing or removing this beamsplitter is random and relativistically space-like separated from the entering of the photon inside the interferometer, as required for the Wheeler's delayed-choice regime. Each set value of R allows us to observe interference with visibility V and to obtain incomplete which-path information characterized by the distinguishability D. The measured values of V and D are found to obey the complementarity relation V 2 + D 2 ≤ 1.

As emphasized by Bohr [1], complementarity lies at the heart of quantum mechanics. A celebrated example is the illustration of wave-particle duality by considering single particles in a two-path interferometer [START_REF] Feynman | Lectures on Physics[END_REF], where one chooses either to observe interference fringes, obviously associated to a wave-like behavior, or to know which path of the interferometer has been followed, according to a particle-like behavior [3]. Although interference has been observed at the individual particle level with electrons [START_REF] Tonomura | [END_REF], neutrons [5], atoms [6,7], molecules [8], only a few experiments with massive particles have explicitly checked the mutual exclusiveness of which-path information (WPI) and interference [9,10,11,12,13].

In the case of photons, it has been pinpointed that meaningful two-path interference experiments demand a single-photon source [14] for which full and unambiguous WPI can be obtained, complementary to the observation of interference [14,15,16]. In order to rule out a too naive view of complementarity, which would assume that the particle could "know" when entering the apparatus which experimental configuration has been set (record of interference or determination of WPI) and would then adjust its behavior accordingly [START_REF] Greenstein | The Quantum Challenge[END_REF], Wheeler proposed the "delayed choice" scheme where the choice between the two complementary measurements is made long after the particle entered the interferometer [START_REF] Wheeler | Quantum Theory and Measurement[END_REF]. Realizations of that gedanken experiment [START_REF] Hellmut | [END_REF]20,21,22] have confirmed that the chosen observable can be determined with perfect accuracy even if the choice, made by a quantum random number generator, is space-like separated from the entering of the particle into the interferometer [22].

In 1978, Wooters and Zurek [23] considered an intermediate situation in which interaction with the interferometer considered as a quantum device allows one to gain an imperfect -but significant -knowledge of WPI, without destroying the interference pattern, which re-mains observable with a good -although reduced -visibility. In 1988, Greenberger and Yasin noticed that in an unbalanced interferometer as used in some neutron interferometry experiments, one has partial WPI while keeping interference with limited visibility [24]. The complementary quantities WPI and interference visibility could then be partially determined simultaneously.

Consistent theoretical analysis of both schemes, independently published by Jaeger et al. [25] and by Englert [26] leads to the inequality [27]

V 2 + D 2 ≤ 1 (1)
which puts an upper bound to the maximum values of simultaneously determined interference visibility V and path distinguishability D, a parameter that quantifies the available WPI on the quantum system. The all-or-nothing cases (V = 1, D = 0) or (V = 0, D = 1) [START_REF] Tonomura | [END_REF]5,6,7,8,14,15,16] obviously fulfill inequality (1). Intermediate situations, corresponding to partial WPI and reduced visibility, have been investigated using atoms [START_REF] Durr | [END_REF], nuclear spins [29] and faint laser light [30]. However none of them has been realized in the delayed-choice scheme. We report here an experimental test of the complementarity inequality (1) in intermediate regimes using true single-photon pulses, and in the delayed-choice operation mode.

Following Englert [26], we point out that the distinguishablility D constrained by inequality (1) actually corresponds to two different notions. The a priori distinguishability, also called "predictability", refers to a WPI obtained by using an unbalanced interferometer with different particle flux along the two paths. Only the case where path distinguishability is introduced a posteriori, i.e. after the entering of the particles into the interferometer, offers the opportunity of a delayed choice test of complementarity. This a posteriori distinguishability can Theroretical sketch of the delayed-choice complementarity-test experiment. A single-photon pulse is sent into a Mach-Zehnder interferometer, composed of a 50/50 input beamsplitter (BS) and a variable output beamsplitter (VBS). The reflection coefficient is randomly set either to the null value or to an adjustable value R, after the photon has entered the interferometer. The singlephoton photodetectors P1 and P2 allow to record both the interference and the WPI.

be introduced either by creating entanglement between the particle and a which-path marker [13,31] or by using an interferometer with an unbalanced output beamsplitter [START_REF] Durr | [END_REF]. We have chosen the latter case by implementing the scheme depicted on Fig. 1, where a single-photon pulse is sent into a Mach-Zehnder interferometer with a variable output beamsplitter (VBS) of adjustable reflection coefficient R. When R is not 0.5, one can have some WPI by observing which detector (P 1 or P 2 ) is fired. The choice of introducing or removing this beamsplitter is random and relativistically space-like separated from the entering of the photon inside the interferometer, as required for the delayed-choice regime.

The experiment starts from a clock-triggered singlephoton source, based on the photoluminescence of a single N-V color center in a diamond nanocrystal [32]. The linearly polarized single-photon pulses are then directed to a polarization Mach-Zehnder interferometer described in Ref. [22]. The input polarization beamsplitter BS splits the light pulse into two spatially separated components of equal amplitudes, associated with orthogonal S and P polarizations. The two beams then propagate in free space for 48 m, which corresponds to a time of flight of 160 ns.

The variable output beamsplitter VBS is the association of a polarization beamsplitter (PBS) which spatially overlaps the two beams, an electro-optical modulator (EOM) which acts as an adjustable waveplate, and a Wollaston prism (WP) with its polarization eigenstates corresponding to the S and P polarized channels of the interferometer (Fig. 2). Given the relative orientation β of the EOM, the VBS reflection coefficient R depends on the voltage V EOM applied to the EOM, according to the relation:

R = sin 2 2β × sin 2 π 2 V EOM V π (2)
where V π is the half-wave voltage of the EOM. The pa-rameters β and V π have been independently measured for our experimental conditions, and found equal to β = 24 ± 1 • and V π = 217 ± 1 V at the wavelengh λ = 670 nm which is the emission peak of the negatively charged N-V color center [32]. This allows R to vary between 0 and 0.5 when V EOM is varied between 0 and 170 V. When R = 0, the VBS is equivalent to a perfectly transparent (or absent) beamsplitter. Then, each "click" of one of the two photodetectors (P 1 or P 2 ) placed on the output ports of the interferometer is associated to a specific path. It then gives access to the full WPI (D = 1) and no interference effect will be observed (V = 0). When R = 0, paths 1 and 2 are partially recombined by the VBS. The WPI is then partially washed out, up to be totally erased when R = 0.5. On the other hand, interference can be observed when the dephasing Φ between paths 1 and 2 is varied. The experiment will consist in checking the relation between D and V for a given value of R, controlled by the EOM voltage V EOM .

In order to perform the experimental test of complementarity in the delayed-choice regime, the choosen configuration of the interferometer, defined by R, has to be causally isolated from the entering of the photon into the interferometer. This condition is ensured by a relativistically space-like separated random choice. For each measurement, the value of the reflection coefficient of VBS is randomly chosen between 0 and a given value of R, using a quantum random number generator (QRNG) located at the output of the interferometer (Fig. 2) [22]. The random numbers are generated from the amplified shotnoise of a white light beam which is an intrinsic quantum random process. At the experiment clock-frequency, i.e. every τ rep = 238 ns, fast comparison of the amplified shotnoise to the zero level generates a binary random number 0 or 1 which changes the VBS reflectivity between 0 and R, by applying or not the corresponding voltage to the EOM (see Eq. ( 2)). In the laboratory framework, the random choise is realized simultaneously with the entering of the photon into the interferometer, ensuring the required space-like separation [22].

As meaningfull illustration of complementarity requires the use of single particles, the quantum behavior of the light field is first tested using the two output detectors feeding single and coincidence counters with no voltage applied to the EOM. In this situation of an absent output beamsplitter, we measure the correlation parameter α [14,16], which is equivalent to the second order correlation function at zero delay g (2) (0). For an ideal single-photon source, quantum optics predicts a perfect anticorrelation α = 0, in agreement with the particlelike image that the photon cannot be detected simultaneously in the two paths of the interferometer. With our source [32], we find α = 0.15 ± 0.01, a value much smaller than one, showing that we are indeed close to the pure single-photon regime [START_REF]α parameter is due to residual background photoluminescence of the diamond sample and to its two-phonon Raman scattering line, which both produce uncorrelated photons associated to Poissonian statistics[END_REF].

The delayed-choice test of complementarity with single-photon pulses is performed with the EOM randomly switched for each photon sent in the interferome-ter, corresponding to a random choice between two values 0 and R of the VBS reflectivity. The phase-shift Φ between the two arms of the interferometer is varied by tilting the polarization beamsplitter PBS of VBS with a piezoelectric actuator (see Fig. 2). For each photon, we record the chosen configuration of the interferometer, the detection events, and the actuator position. All raw data are saved in real time and are processed only after a run is completed. The events corresponding to each configuration of the interferometer are finally sorted. For a given value R, the wave-like information of the light field is obtained by measuring the visibility of the interference, predicted to be

V = 2 R(1 -R) .
(3)

The results, depicted in Fig. 3, show a reduction of V when the value of R randomly applied decreases.

To test inequality (1), a value of the distinguishability D is then required, to qualitatively qualify the amount of WPI which can be extracted for each value of R. We introduce the quantity D 1 and D 2 , respectively associated to the WPI on path 1 and path 2:

D 1 = |p(P 1 , path 1) -p(P 2 , path 1)| ( 4
)
D 2 = |p(P 1 , path 2) -p(P 2 , path 2)| (5)
where p(P i , path j) is the probability that the particle follows path j and is detected on detector P i . For a single particle arriving on the output beamsplitter, one obtains

1 = D 2 = 1 2 -R . (6)
The distinguishability parameter D is finally defined as [26]

D = D 1 + D 2 = 1 -2R . (7)
In order to test this relation, we estimate the values of D 1 and D 2 by blocking one path of the interferometer and measuring the number of detections N 1 and N 2 on detectors P 1 and P 2 , which are statisticallly related to D 1 and D 2 according to [22,[START_REF] Durr | [END_REF] :

D 1 = N 1 -N 2 N 1 + N 2 path 2 blocked (8)
D 2 = N 1 -N 2 N 1 + N 2 path 1 blocked . (9)
These measurements are also performed in the delayedchoice regime, using the procedure described above. We finally obtain independent measurement of D and V for different values of the reflection coefficient R, randomly applied to the interferometer. The final results, depicted on Fig. 4, leads to V 2 + D 2 = 0.97 ± 0.03, close to the ideal balance between V and D constrained by inequality (1), eventhough each quantity varies from 0 to 1.

The effects observed in this delayed-choice experiment are in perfect agreement with quantum mechanics predictions. No change is observed between a so called "normalchoice" experiment and the "delayed-choice" version. It demonstrates that the complementarity principle cannot be interpreted in a naive way, assuming that the photon at the input of the interferometer could adjust its nature according to the experimental setup installed. As Bohr pointed out [START_REF] Bohr | Albert Einstein: Philosopher Scientist[END_REF], "it obviously can make no difference as regards observable effects obtainable by a definite experimental arrangement, whether our plans of constructing or handling the instrument are fixed beforehand or whether we prefer to postpone the completion of our planning until a later moment when the particle is already on its way from one instrument to another". Such intriguing property of quantum mechanics forces one to renounce to some common-sense representations of the physical reality.

We warmly thank A. Clouqueur and A. Villing for the realization of the electronics, and J.-P. Madrange for all mechanical realization of the interferometer. This work is supported by Institut Universitaire de France.

FIG. 1 :

 1 FIG. 1:Theroretical sketch of the delayed-choice complementarity-test experiment. A single-photon pulse is sent into a Mach-Zehnder interferometer, composed of a 50/50 input beamsplitter (BS) and a variable output beamsplitter (VBS). The reflection coefficient is randomly set either to the null value or to an adjustable value R, after the photon has entered the interferometer. The singlephoton photodetectors P1 and P2 allow to record both the interference and the WPI.

FIG. 2 :

 2 FIG. 2: Variable output beamspltter (VBS) implementation.The optical axis of the polarization beamsplitter (PBS) and the polarization eigenstates of the Wollaston prism (WP) are aligned, and make an angle β with the optical axis of the EOM. The voltage VEOM applied to the EOM is randomly chosen accordingly to the output of a Quantum Random Number Generator (QRNG), located at the output of the interferometer and synchronized on the 4.2 MHz clock that triggers the single-photon emission.

FIG. 4 :

 4 FIG. 4: Delayed-choice test of complementarity with singlephoton pulses. (a)-Wave-like information V 2 and which-path information D 2 as a function of the EOM voltage corresponding to a given value R of the VBS reflectivity. The solid lines are the theoretical expectations, with β = 24 • and Vπ = 217, using Eqs (2), (3) and (7). (b)-V 2 + D 2 as a function of the EOM voltage.

D


Panayotis Spathis

Sylvain Colson

Cornelis Feng Yang

Jacominus Van Der Beek

Piotr Gierlowski

Takasada Shibauchi

Yuji Matsuda

Marat Gaifullin

Ming Li

Peter H Kes

Cornelis J Van Der Beek

Piotr Gier

c-axis coupling in underdoped Bi2Sr2CaCu2O8 with varying degrees of disorder

Keywords: numbers: 74.20.De, 74.25.Bt, 74.25.Dw, 74.25.Op, 74.25.Qt, 74.72.Hs

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.

I. INTRODUCTION

Significant controversy remains concerning an appropriate model description of high temperature superconducting cuprates (HTSC) in the underdoped regime, i.e. the regime in which the number of additional holes per Cu, p is smaller than the value 0.16 at which the critical temperature T c is maximum. 1 Whereas it is well established that the charge dynamics and transport properties in the normal-and superconducting states in the overdoped regime (p > 0.16) are, by and large, determined by well-defined quasi-particles, the role of quasiparticles in the underdoped regime is debated. The underdoped region of the (p, T ) phase diagram is characterized by several salient features. 2 At T * > T c , the wellknown "pseudo-gap" in the excitation spectrum opens up. This has been interpreted as either being related to the advent of another type of (spin-or charge-) order competing with superconductivity, driving T c down as p is diminished, or, alternatively, as signalling the formation of precursor Cooper pairs without long-range phase coherence. Then, T c is interpreted as the demise of long range superconducting phase order due to strong thermal 3,4,5 or quantum 6,7 phase fluctuations. Strong support for this scenario has come from the violation of the Glover-Tinkham-Ferrell conductivity sum-rule applied to the c-axis spectral weight; 5 also, the linear relation between T c and the superfluid density 8 has been interpreted as the result of T c being determined by phase fluctuations in a Kosterlitz-Thouless-type scenario. 9,10 A smoking gun for such a scenario would be an important reduction of the c-axis superfluid density ρ c s , which is determined by Cooper pair and quasiparticle tunneling between adjacent strongly superonducting CuO 2 layers through the weakly superconducting rocksalt-like blocking layers, with respect to the in-plane stiffness ρ ab s in underdoped cuprates below T c . However, apart from phase fluctuations, other mechanisms for the reduction of ρ c s , arising from disorder in the crystalline structure of underdoped cuprates, cannot be ignored. 11 First, scanning tunneling spectroscopy (STS) experiments 12,13,14 have revealed large variations of the magnitude of the gap maximum ∆ peak , as interpreted from conductance curves measured on the surface of Bi 2 Sr 2 CaCu 2 O 8+δ crystals. This has motivated recent interpretations of weakened c-axis superfluid response in this material 15,16 as well as in La 2-x Sr x CuO 4 17 in terms of finely dispersed 5 -20 nm-sized non-superconducting regions within the CuO 2 planes. Such regions may arise from the suppression of the superconducting order parameter by dopant atoms, 18 such as out-of-plane oxygen atoms in the Bi 2 Sr 2 CaCu 2 O 8+δ compound. 19 Moreover, the d-wave symmetry of the gap function is at the basis of several mechanisms by which pointlike disorder reduces the c-axis superfluid density. The appearance of quasiparticle (virtual-) bound states and their smearing by a finite defect density leads to an increase of the density of states (DOS) near the nodal (π, π) directions (the so-called lifetime effect). 20 The same pointlike defects increase the quasiparticle scattering rate Γ s . It was conjectured that in the case of coherent (in-plane momentum preserving) quasiparticle tunneling, the cancellation of these effects leads to disorder-independent low-temperature c-axis quasiparticle conductivity and superfluid density. 21 The approach of Ref. [21] neglects the crystal structure of the tetragonal HTSC, 11 which leads to the dependence of the interlayer hopping integral t ⊥ on the in-plane momentum (k x , k y ). For simple tetragonal structures, interlayer hopping occurs via Cu 4s orbitals in adjacent planes. Its momentum dependence t ⊥ = t 0 ⊥ [cos k x a -cos k y a] 2 is determined by the in-plane overlap of the bonding oxygen 2p level with the 4s level of the neighboring Cu atom. 22,23 As a result, c-axis tunneling occurs nearly exclusively for the anti-nodal directions at which quasiparticles are unlikely to be excited. In body-centered tetragonal structures such as Bi 2 Sr 2 CaCu 2 O 8+δ , hopping is also suppressed along the (k x , k y ) = (π, 0) and (0, π) lines, yielding 24 In either case, disorder is always relevant for the nodal directions. Then, from the lifetime effect, one expects a quadratic decrease with temperature of the reduced c-axis superfluid density 23,25 ρ c s (T )

t ⊥ = t 0 ⊥ [cos k x a -cos k y a] 2 cos 1 2 k x a cos 1 2 k y a.
ρ c s (0) ∝ 1 -α c 8π 3 Γ s ∆ 0 T ∆ 0 2 , (k B T ≪ Γ s ). (1)
Here α c is a dimensionless constant of order unity and the parameter ∆ 0 was assumed, in Refs. 23 and 25, to correspond to the maximum amplitude of the Bardeen-Cooper-Schrieffer d-wave gap. Eq. ( 1) essentially differs from that derived for the ab-plane superfluid density

ρ ab s (T ) ρ ab s (0) ∝ 1 -α ab ∆ 0 Γ s T ∆ 0 2 , (k B T ≪ Γ s ) (2)
in that the leading temperature-dependent term has a coefficient that is smaller by a factor (Γ s /∆ 0 ) 2 . 23 The presence of defects in the rocksalt-like (BiO) layers tends to break the d-wave symmetry of the hopping integral, and renders quasiparticle hopping possible for other values of the in-plane momentum, and notably along the order parameter nodes. 11,23,25,26 A condition for this "impurityassisted hopping" (IAH) to be effective is an anisotropic scattering matrix of the interplane defects. Notably, for strong forward scattering, the result 23,27 Here V 1 is the magnitude of the impurity scattering potential of the out-of-plane defects and N (E F ) is the density of states at the Fermi level in the normal state. The effect of impurities can be distinguished from that of boson-assisted interlayer hopping; for the latter, a very similar result is obtained, but with the leading temperature-dependent term proportional to T 3 . 25,28 Finally, direct hopping of quasiparticles was suggested to lead to a small, linearly temperature-dependent, reduction of ρ c s . 27 In this paper, we address the mechanism by which the c-axis superfluid density in underdoped Bi 2 Sr 2 CaCu 2 O 8+δ ( with p = 0.10 ) is reduced by using disorder, in the form of Frenkel pairs introduced by high energy electron irradiation, as an independent control parameter. Electron irradiation, the effects of which are taken to be similar to those of Zn-doping, 29 has previously been used to study the effect of pointlike disorder on the resistivity, critical temperature, 30 and Nernst effect of YBa 2 Cu 3 O 7 and YBa 2 Cu 3 O 6.6 . 31 In the latter material, electron irradiation eventually leads (at high fluences) to the breakdown of the well-known Abrikosov Gor'kov relation 32,33 ln

ρ c s (T ) ≈ 2πV 1 ∆ 0 N 2 (E F ) 1 -8 ln 2 T ∆ 0 2 , (3)
was obtained for Γ s ≪ k B T ≪ 1 2 2πV 1 ∆ 0 N (E F )/(t 0 ⊥ ) 2 1/3 T c .
T c T c0 = Ψ 1 2 -Ψ 1 2 + Γ 2πk B T c (4)
(with T c0 the critical temperature when the normal state scattering rate Γ is equal to zero, and Ψ the digamma function) as well as a significant increase of the fluctuation regime near T c . 31 Both effects were interpreted as the effect of strong superconducting phase fluctuations. 30,31 The in-plane and c-axis superfluid densities of Zn-doped YBa 2 Cu 3 O 7-δ were studied by Panagopoulos et al. 34 and by Fukuzumi, Mizuhashi, and Uchida. 35 The progressive inclusion of Zn leads to a rapid decrease of the in-plane superfluid density ρ ab s , corresponding to an increase of the in-plane penetration depth λ ab (0) ∝ (ρ ab s ) -1/2 , and a more modest decrease of ρ c s ∝ σ c (T c ) ∝ T c , that violates the c-axis conductivity sum-rule 5,35 Below, we report on c-axis coupling in the superconducting state measured through the Josephson Plasma Resonance (JPR), 40,41,42,43 which, in our underdoped Bi 2 Sr 2 CaCu 2 O 8+δ crystals takes place in the microwave frequency regime below 70 GHz. The JPR frequency f pl is sensitive to the value of ∆ 0 , as well as to fluctuations of the superconducting order parameter phase in the CuO 2 planes. 6 The evolution of f pl (T ) with temperature depends simultaneously on the quasiparticle dynamics and on the strength of fluctuations; the plasma resonance peak is broadened both by the quasiparticle tunneling rate and by crystalline disorder. 17,44 However, the dependence of f 2 pl ∝ ρ c s on the disorder strength is expected to be quite different, depending on which mechanism is predominant. In the following, we show that the disorder dependence of the c-axis plasma frequency is a sensitive probe, that allows one to identify in detail what physical mechanisms are at the basis of the reduction of the superfluid density in Bi 2 Sr 2 CaCu 2 O 8+δ . It turns out that, even in our heavily underdoped crystals, (incoherent) c-axis quasiparticle hopping is essential for a consistent description of the data. We find that the energy scale ∆ 0 , which turns out to be Delta 0 ≈ 2.5k B T c for all underdoped crystals, is to be interpreted as an energy scale governing nodal quasi-particle excitations.

II. EXPERIMENTAL DETAILS

The underdoped (T c = 65 ± 0.5 K, p ≈ 0.10) Bi 2 Sr 2 CaCu 2 O 8+δ single crystals, of typical dimensions 500 × 300 × 40µm 3 , were selected from the same boule, grown by the travelling solvent floating zone method at the FOM-ALMOS center, the Netherlands, in 25 mBar O 2 partial pressure. 45 The crystals were annealed for one week in flowing N 2 gas. We have also measured a set of optimally doped control samples ( T c = 86 K ). These were also grown by the travelling solvent floating zone technique, at 200 mbar oxygen partial pressure, and subsequently annealed in air at 800 • C. The crystals were irradiated with 2.3 MeV electrons using the Van de Graaff accelerator at the Laboratoire des Solides Irradiés. The beam was directed along the crystalline c-axis during the irradiation. To prevent recombination and clustering of point defects, the irradiation is carried out with the crystals immersed in a liquid hydrogen bath ( 22 K ). The electron flux is limited to 2 × 10 14 e -cm -2 per second. Crystals UD5-UD8 were irradiated to a total fluence of 0.53×10 18 , 3×10 18 , 7.7×10 18 , and 8.8×10 19 e -cm -2 respectively. After measurements, crystal UD5 was irradiated a second time to a total fluence of 6.0×10 19 e -cm -2 and was henceforth labeled UD5b. The high energy electron irradiation creates random atomic displacements in the form of Frenkel pairs, both in the CuO 2 bilayers and in the intermediate cation layers, throughout the samples.

The superconducting transition temperature T c was determined by ac susceptibility measurements using a driving field of amplitude 4.2 mOe and a frequency of 560 Hz, directed parallel to the c-axis. For all underdoped crystals, the superconducting transition is rather broad. The superconducting transition takes place in two steps: there is a slow increase of screening at high temperature, followed by a rapid step of the diamagnetic signal at lower temperature ( Figure 1 ). The high temperature screening vanishes when the excitation field amplitude is increased beyond 0.5 Oe, while the step at lower temperature is robust. This shows that doping is macroscopically inhomogeneous, and that the crystals are surrounded by a thin surface layer of higher doping. This layer could not be eliminated by cutting the crystals. The overall shape of the transition is unaffected by the electron irradiation. The transition widths are of the order of 4 K, which is usual for such low doping. After irradiation, the transition widths slightly increase. For all crystals, the transition to zero dc resistivity and bulk superconductivity occurs at the temperature at which the lower screening step takes place, see e.g. the Inset to Fig. 2. Therefore, the lower temperature feature was adopted as characterizing the bulk T c of the underdoped crystals. Crystals were further characterized by the measurement of the temperature variation of the in-plane penetration depth, λ ab (T )/λ ab (0) -1 ≡ ∆λ ab /λ ab . For this, a crystal is mounted on a sapphire rod, in the center of a superconducting (Pb) resonant cavity immersed in liquid 4 He, and operated in the TE 011 mode. The cavity resonant frequency was f ∼ 27.8 GHz, and the quality factor Q ∼ 4 × 10 5 . The crystal is mounted in such a way that the magnetic microwave field is perpendicular to its ab plane and solely in-plane screening currents are induced. From the shift ∆f of the cavity resonance frequency induced by the sample, we determine the surface reactance X s = 2πµ 0 G 2 ∆f (with µ 0 = 4π × 10 -7 Hm -1 ). The surface resistance R s = 2πµ 0 G 1 (∆Q) -1 was obtained from the change of the quality factor. The geometrical factors G 1 and G 2 were determined by comparing the surface impedance in the normal state, X s = R s = πµ 0 f δ s , to the value expected from the normal state resistivity, ρ = πµ 0 f δ 2 s . 46 It was retrospectively checked that all measurements were carried out in the skin effect regime, in which the normal state skin depth δ s is much smaller than the sample dimensions. The relative change of the penetration depth was determined from the behavior of the surface reactance at low temperature, at which X s ≈ 2πµ 0 f λ ab .

Figure 2 shows that the temperature at which the decrease of the surface resistance R s was observed corresponds to the (high temperature) onset of screening in the ac susceptibility measurement. This indicates that the surface skin depth of thickness ∼ 7µm probed by the microwave field contains patches with larger hole content p.

The JPR measurements were performed using the cavity perturbation technique, using the TM 01n modes (n = 0, ..., 4) of two Oxygen-Free High Conductivity (OFHC) copper resonant cavities ( Q ∼ 4000 -10000 ), mounted on a cryocooler cold head. The measurement frequencies ranged between 19.2 and 39.6 GHz. 47 Further measurements were made applying the bolometric technique, using waveguides in the TE 01 travelling wave mode. 39 In both measurement set-ups, the electrical microwave field is applied along the c-axis of the crystal. In contrast to the previously described surface impedance measurements, screening of the electric microwave field is very poor because of the high electronic anisotropy of the crystals. The underdoped samples are in the complete depolarization regime and thus the bulk electromagnetic response is probed. By monitoring the power absorption as a function or frequency for a fixed temperature, the Josephson Plasma Resonance is detected as a sharp absorption peak in the microwave response, see Fig. 3. We determine the JPR frequency at a given temperature, f pl (T ), as the measurement frequency at the temperature at which dissipation is maximum.

III. RESULTS

Figure 4 collects the values of the critical temperature as function of electron dose, for the set of underdoped samples as well as the optimally doped control samples. Both the T c of the underdoped and the optimally doped crystals decrease linearly with irradiation fluence. The derivative of T c with respect to fluence of the optimally doped crystals concurs with that measured by Behnia et al. 37 and Nakamae et al. 38 but is two times lower than that measured by Rullier-Albenque et al.. 36 The overlap between the variation with fluence of the screening onset temperature in the underdoped crystals and the T c of the optimally doped crystals shows that the thin surface layer on the underdoped Bi 2 Sr 2 CaCu 2 O 8+δ has p ∼ 0.16. The critical temperatures, normalized to the critical temperatures T c0 of the unirradiated crystals, can be superimposed on Eq. ( 4), yielding estimates of the normal state scattering rate Γ (see Fig. 5). This procedure supposes that T c0 corresponds to the critical temperature in the absence of disorder; we shall see below that this is not justified, so that the estimated Γ values are in fact lower limits for each crystal.

The relative change with temperature of the in-plane penetration depth λ ab is depicted in Fig. 6. For all underdoped crystals, including the unirradiated ones, λ ab (T ) -λ ab (0) varies quadratically with temperature at low T . Such a temperature dependence has been associated with quasiparticle scattering in the unitary limit by point defects situated within the CuO 2 planes of the d-wave superconductor, i.e. N (E F )V ≫ 1, Γ ∼ n d /πN (E F ), and Γ s ∼ 0.6(Γ∆ 0 ) 1/2 . 48 Here V is the scattering potential of the defects in the planes, with density n d . The quadratic temperature dependence of ∆λ ab (T ) is at odds with a possible important role of thermal phase fluctuations, for which a linear T -dependence was predicted. 49 As for the magnitude of the T 2 -contribution to λ ab , a very modest change is found for the lower irradiation fluences. Only for fluences exceeding 10 19 e -cm -2 does the in-plane penetration depth increase significantly with defect density. We now switch to the central results concerning the Josephson Plasma Resonance. Figure 7 shows the JPR frequency f pl (T ) of crystals UD5-UD7 as function of temperature ( measured in Earth's magnetic field ). The temperature at which f pl (T ) extrapolates to zero is welldefined and corresponds to the critical temperature of the bulk, underdoped portion of the crystals, i.e. the main transition in the ac susceptibility and zero resistance. This shows that the JPR probes the c-axis response in the heavily underdoped bulk and is insensitive to the surface quality of the samples. From Fig. 7 one sees that not only T c , but also f pl is strongly depressed by the electron irradiation. Fig. 8 collects values of the low-temperature extrapolated value f pl (0) versus the critical temperature, and reveals the proportionality between f 2 pl (0) and T c . This dependence is clearly different from the variation of f pl with oxygen doping. The same Figure recapitulates results for doping levels p = 0.13 39 and 0.11, 47 the evolution of which recalls the exponential f 2 pl (0)(T c )-dependence found by Shibauchi et al.. 16 The results are somewhat similar to those obtained by Fukuzumi et al. for Zn-doped and oxygen-deficient YBa 2 Cu 3 O 7-δ : 35 however, their Fig. 5 shows a linear dependence of f pl on T c for underdoped YBa 2 (Cu 1-x Zn x ) 3 O 6.63 , and a hyperbolic or exponential evolution of f pl (T c ) as function of δ.

As for the temperature variation of f 2 pl (T ) ∝ ρ c s , it turns out to be identical for all the underdoped crystals, including the pristine ones, and does not depend on the defect density. Josephson (c-axis) critical current,

j c c = 2π(h/2e)ǫǫ 0 s -1 f 2 pl = h 4πeµ 0 λ 2 c s (5)
normalized to its value for T → 0. Here s is the spacing between CuO 2 planes, the c-axis dielectric constant ǫ ≈ 11.5, 50 and ǫ 0 = 8.854 × 10 -12 Fm -1 .

IV. DISCUSSION

To distinguish between the different mechanisms responsible for the reduction of the c-axis plasma frequency as temperature is increased, we dispose of three tools. First, there is the variation of f 2 pl (0) with disorder strength, which manifests itself starting from the smallest electron fluences, and proportionally follows the evolution of T c with disorder. Second is the temperature variation f 2 pl (T )/f 2 pl (0), that follows a 1 -a(T /T c ) α dependence, with α ∼ 2 independent of the disorder strength. Finally, there is the comparison with the behavior of the in-plane penetration depth, λ ab (T )/λ ab (0) ∼ 1 + βT 2 . A successful model description should account for all three dependences correctly.

The theoretically expected low-temperature dependence λ c (T ) depends strongly on a number of circumstances. First is the question whether superconductive coupling is three-dimensional ( i.e. the c-axis momentum k z is a good quantum number ), 23 or whether it is mediated by Josephson tunneling between CuO 2 layers. Josephson coupling can be weakened by direct 21,27 or boson-assisted quasiparticle tunneling, 25,27,28 or by tunneling that involves intermediate defect-induced states between the layers (IAH). 23,27 For "direct" tunneling, e.g. occurring through direct overlap of the superconducting wave functions in the CuO 2 planes, one distinguishes the case of conserved in-plane momentum k ("coherent" tunneling) from the case where it is not preserved ("incoherent" tunneling 21 -this situation would yield a vanishing j c c in a d-wave superconductor). A prevailing effect of nodal quasiparticles leads to a more rapid decrease of j c c , see Eq. ( 1). Finally, in all cases the anisotropy of the transfer matrix t ⊥ is expected to have an important influence on the temperature dependence of λ c , 23,52 notably reconciling a weak T -dependence of λ c with a strong variation of λ ab (T ). Here, the experimentally observed temperature dependence of f 2 pl (T )/f 2 pl (0) actually allows one to discard a dominant role of a possible d-wave symmetry of the transfer matrix t ⊥ , 23 since ( for k B T ≫ Γ s ) this leads to a weak temperature dependence, O 8+δ compound will influence the maximum c-axis critical current. However, it will not change the expected 1 -ãT 5 temperature dependence, since this arises from the specific coincidence in k-space of the zero of t ⊥ (k x , k y ) with the nodal direction of the order parameter. Thus, models for superconductive coupling, 23 or direct Josephson coupling with a vanishing hopping integral along the nodal line 52 are in inadequacy with the data.

f 2 pl (T )/f 2 pl (0) ∼ 1 -ãT 5 ,
We next exclude a dominant role of direct quasiparticle tunneling. Even though the similar T 2 -dependences of the low-temperature ab-plane and c-axis penetration depths suggests such coupling, the disorder dependence is at odds with the experimental result. Radtke et al. 27 and Latyshev et al. 21 find that for direct coupling, the low temperature c-axis critical current

j c,direct c = πσ c q (0)∆ 0 es = 4πet 2 ⊥ N n (E F ) h (6)
is, for k B T ≪ Γ s ∼ 20 -30 K, to lowest order independent of the defect density due to the cancellation of the scattering-rate dependences of the quasi-particle conductivity σ c q and ∆ 0 . The model was further worked out by Kim and Carbotte, who find that to first order

j c,direct c ∝ 1 -α Γ s Γ 2 s + ∆ 2 0 ∼ 1 -α 1 1 + ∆ 0 /Γ (7)
both for the case of constant t ⊥ (where α ≈ 1) and angular-dependent t ⊥ (α = 16 9 ). 11 The nonlinear dependence ( 7) is at odds with the observed linear evolution of j c c with irradiation fluence. The temperature and disorder-dependence of the low temperature c-axis JPR frequency is more successfully described by a model for incoherent tunneling. According to Latyshev et al., an incoherent tunneling process yields

j c,i c ≈ j c,direct c ∆ 0 /E F ≈ 4πet 2 ⊥ N n (E F )∆ 0 /hE F .
The extra factor ∆ 0 then explains the linear relation between

f 2 pl (T )/f 2 pl (0) ∼ λ 2 c (0)/λ 2 c (T ) ∼ ρ c s (T )/ρ c s (0) ∼ j c c (T )/j c c (0)
. The drawn line is a fit to Eq. ( 3) with ∆0 = 2.5kB Tc. j c c (0) and T c , accepting that in a d-wave superconductor with impurity scattering in the unitary limit, ∆ 0 (Γ) is simply proportional to T c (Γ). 33 The linear dependence on ∆ 0 is found in the IAH model, see Eq. (3). 23,27,28 The latter expression consistently describes the fact that the temperature dependence of f pl does not change with defect density: the temperature dependent term writes (T /∆ 0 ) 2 ∝ (T /T c ) 2 . For the same reason, the "lifetime effect", Eq. (1), does not describe the data: in the unitary limit, the leading temperature-dependent term has an extra factor Γ s /∆ 0 ∼ (Γ/∆ 0 ) 1/2 and is therefore expected to strongly depend on defect density. The observed defect-density independence of ρ c s (T )/ρ c s would require the strength of the scattering potential of the individual irradiation defects in the CuO 2 planes to be in the Born limit, which contradicts the results on the temperature dependence of the ab-plane penetration depth. We note that the toy model for incoherent hopping of Ref. [11], which yields λ

-2 c ∝ [1 -5 14 (Γ s /∆ 0 ) 2 -. . .] ∼ [1 -3 14 (Γ/∆ 0 ) -. . .]
, also describes the initial linear decrease of j c c (T c ) (Eq. ( 31) of Ref. [11]). Given the success of the IAH-model in qualitatively explaining the temperature-as well as the disorder dependence of the JPR data, we perform a direct fit of f 2 pl (T )/f 2 pl (0) to Eq. (3), shown in Fig. 9. The only parameter is the value ∆ 0 = 2.5k B T c . This value not only means that Γ s > 30 K, comforting our interpretation of the decrease of j c c ∝ T c in terms of unitary scatterers induced in the CuO 2 planes, it is also remarkably close to characteristic energy scales found in recent Raman scattering 53 and STM experiments. 54 The first study finds that in underdoped HgBa 2 CuO 4 , the B 2g Raman mode, which directly couples to the same lowenergy nodal quasi-particle excitations that are responsible for the reduction of the c-axis superfluid density in the IAH model, is characterized by an energy scale ∼ 2k B T c (related to the nodal slope of the gap function). 53 We surmise that ∆ 0 is a closely related parameter. The second study finds that the total tunneling gap amplitude is determined by two energy scales, the smaller of which, ∆ φ ≈ 2.8k B T c , is related to superconductivity. 54 It is interesting that the existence of a second small energy scale describing nodal quasiparticle exitations may well be responsible for the observed suppression of j c c as function of doping (decreasing p or δ). 16 The proportionality of this decrease to the ratio of the gap maximum (at the antinode) and T c finds a logical explanation if the smaller energy scale ( ∼ 2 -2.5k B T c ) is what determines the magnitude of the c-axis critical current density.

Thus, the (incoherent) interlayer assisted hopping is a feasible candidate for the reduction of the c-axis superfluid density in underdoped Bi 2 Sr 2 CaCu 2 O 8+δ : it parametrically describes the data, and numerical values extracted for the relevant energy scale determining the quasiparticle excitations is the same as found in other experiments. The model does have several caveats: first, there is the above-mentioned puzzle that it requires the temperature dependence of the ab-plane superfluid density to be explained by the lifetime effect, whereas the same effect does not seem to play a role in the caxis superfluid density, other than providing the quasiparticles. The second is the identification of the defects in the rocksalt-like layers that are responsible for interlayer scattering. The model by Xiang et al. 23,25 requires these out-of-plane defects to be weakly scattering, with a strongly anisotropic potential that leaves the reflected wavevector close to the incident ("strong forward scattering"). Although candidates may be out-of-plane oxygen defects 19 or cation disorder, the constraints imposed on the scattering potential seem very strict. In the end, the agreement of defect-density independence of the experimental f pl (T )/f pl (0) with the IAH prediction 23 may be completely fortuitous. The coincidence of the temperature dependence of f pl of the pristine crystals with even the most heavily irradiated ones indicates that disorder plays an important role in all samples. Notably, we expect T c of pristine crystals to be substantially suppressed with respect to that of hypothetically "clean" underdoped Bi 2 Sr 2 CaCu 2 O 8+δ . Furthermore, impurity scattering is likely to completely suppress any role of quasi-particles in the c-axis electromagnetic response of this compound, leaving only pair tunneling. 6,55 The remaining mechanism for the reduction of the superfluid density in the presence of pair tunneling only is that of order parameter phase fluctuations in the CuO 2layers. The effect of quantum phase fluctuations on the ab-plane and c-axis superfluid densities was examined by Paramekanti et al., 6,7 who performed an analytical study of an XY model for the superconducting order phase φ on a two-dimensional (2D) lattice of spacing ξ 0 (representing the coherence length), within the self-consistent harmonic approximation. The authors conclude that, in contrast to thermal phase fluctuations, quantum phase fluctuations in quasi 2D high temperature superconductors are important at all temperatures. The low carrier density in these materials leads to inefficient screening of the Coulomb interaction between charge carriers, and a sizeable reduction of the magnitude of the ab-plane superfluid density (without change of its temperature dependence). Within this model, the JPR frequency is also renormalized downwards because of fluctuations of the phase in the layers: ⊥ ≈ e 2 /4πǫǫ 0 ξ 0 ε 0 (0)s, 7 with ε 0 s = h 2 s/16e 2 πµ 0 λ 2 ab (0) the in-plane phase stiff-We observe that, given the Uemura relation ε 0 s ∝ k B T c , 8 the resulting expression naturally describes the experimentally observed exponential depression of j c c as function of doping. 16 However, even if one explicitly develops the dependence of ε 0 s in terms of the variance of the in-plane phase φ 2 ∼ φ 2 ⊥ , Eq. ( 8) fails to describe the linear f 2 pl (T c )-dependence (Fig. 8). Therefore, the reduction of ρ c s with increasing disorder cannot be ascribed to quantum phase fluctuations only -pair-breaking in the CuO 2 layers must play a significant role.

f 2 pl (T ) = f 2 pl (0)e -1 2 φ 2 ⊥ ≈ f 2 pl (0) 1 - 1 2 φ 2 ⊥ -. . . . (8
A noteworthy prediction of the quantum fluctuation scenario is that the temperature evolution of the caxis superfluid density is entirely determined by that of the in-plane phase stiffness, i.e. the c-axis and abplane superfluid densities follow the same dependence at low temperature. Inserting the experimental result λ ab = λ ab (0)(1 + βT 2 ) into the prediction of Ref. [7], one would expect

∂ [f pl (T )/f pl (0)] 2 ∂(T /T c ) 2 = - C 1 4 βT 2 c 2πe 2 4πǫǫ 0 ξ 0 ε 0 (0)s , (9)
This is experimentally verified; taking the data of Fig. 6 and λ ab (0) ≈ 300 nm, 47 we find that Eq. ( 9) is obeyed with C 1 ≈ 0.3 ( C 1 should be order unity 7 ). The experimental independence of f pl (T )/f pl (0) on defect density demands that βT 2 c λ ab (0) is disorder independent. Adopting the generally accepted view that λ -2 ab (T ) is described by Eq. (2) with Γ s in the unitary limit, this would imply that βT 2 c λ ab (0) ≈ α ab λ ab (0)(T c /Γ) 1/2 (T c /∆ 0 ) 3/2 , and therefore that λ ab (0) -2 ∝ T c /Γ. This is as yet unverified, as the different sizes and aspect ratios of our crystals prohibit a direct comparison of the absolute values of λ ab .

Note that the case of screening by nodal quasiparticles was also studied in Ref. [7]. Then,

∂ [f pl (T )/f pl (0)] 2 ∂(T /T c ) 2 ≈ - 2βT 2 c σ q , (10)
where σ q is the ab-plane quasiparticle sheet conductivity, normalized to the quantum conductivity e 2 /h. This formula also describes the temperature dependence of the data, provided that σ q ≈ 3; moreover, the ratio β/σ q ∼ α ab ∆ 0 (T c /∆ 0 )m/N s (0)e 2 should be disorderindependent (N s (0) is the quasiparticle density of states and m the effective mass). Given the theoretical expectation N s (0) ∼ Γ 1/2 [48] and ∆ 0 ∝ 1 -Γ [33], this model again fails to describe the reduction of the zero temperature c-axis superfluid density in terms of quantum fluctuations only.

V.

SUMMARY AND CONCLUSIONS

We have cross-correlated the dependence of the c-axis Josephson Plasma Resonance frequency in heavily underdoped Bi 2 Sr 2 CaCu 2 O 8+δ on temperature and controlled disorder (introduced by high energy electron irradiation), and compared both with the behavior of the in-plane penetration depth. It is found that the c-axis critical current is depressed with increasing disorder strength, proportionally to the critical temperature T c . Both the in-plane-and out-of-plane superfluid densities follow a T 2 temperature dependence at low T . The temperature dependence of the c-axis response is independent of disorder, indicating that we are probing the superfluid density. The superfluid response of the pristine underdoped crystals is indistinguishable from that of heavily irradiated ones, suggesting that pristine underdoped Bi 2 Sr 2 CaCu 2 O 8+δ commonly contains sufficient disorder in the CuO 2 planes for the critical temperature to be significantly suppressed with respect to what the T c of the hypothetically "clean" material would be. The dominating in-plane disorder in as-grown crystals is likely to be of the same kind as introduced by the electron irradiation. Apart from unitary scatterers in the CuO 2 planes, this also encompasses the "order parameter holes" induced by dopant oxygen and cation disorder in the rocksalt-like layers. 18,19 The experimental data were confronted with a variety of theoretical models describing the reduction of c-axis superfluid density in terms of either quasiparticle dynamics or quantum fluctuations of the superconducting order parameter phase in the CuO 2 layers. We find that the quantum phase fluctuation description 6,7 yields excellent agreement as to the experimentally observed similar temperature dependences of λ ab and λ c , and quantitatively describes the temperature derivative ∂[f pl (T )/f pl (0)] 2 /∂(T /T c ) 2 . However, it fails to describe the dependence of the zero-temperature Josephson Plasma frequency on disorder strength.

We therefore surmise that the reduction of f pl (0) with increasing disorder must be due to pair-breaking within the CuO 2 layers. Data for λ ab (T ) and λ c (T, Γ) are in agreement with scattering in the unitary limit by the irradiation-induced point defects in the CuO 2 -planes. Only one model consistently describes all aspects of the reduction of the c-axis superfluid density with temperature and disorder strength. This is the Impurity Assisted Hopping model of Radtke et al., 27 elaborated upon by Xiang and Wheatley, 23,25 and by Kim and Carbotte. 11 The model supposes a reduction of ρ c s (T ) through hopping of nodal quasi-particles assisted by weak, highly anisotropic scattering by defects in the insulating SrO and BiO layers. Candidates for such impurities are out-of-plane oxygen defects 19 or cation disorder. From a fit of f pl (T ) to the IAH model, we extract the energy scale ∆ 0 ∼ 2.5k B T c characterizing nodal quasiparticle excitations. This is remarkably close to the number obtained by Le Tacon et al. from anisotropic Raman scattering, 53 giving further confidence in the IAH interpretation.

  [σ c (T c ) is the c-axis conductivity at T c ]. As for the low-T temperature dependence, a gradual change of both ρ ab s and ρ c s from T -linear to T -squared has been reported.

FIG. 1 : 2 FIG. 2 :

 122 FIG.1:(color online) Real part of the ac magnetic susceptibility of the investigated underdoped Bi2Sr2CaCu2O 8+δ crystals. The crystals were irradiated at 22 K with 2.3 MeV electrons to the indicated fluences. The ac field amplitude hac = 4.2 mOe, the ac frequency was 560 Hz. The curves show a screening onset determined by a surface layer containing optimally doped material; this screening is suppressed when hac 0.5 Oe. The steep drop at the lower end of the transitions corresponds to bulk screening by the underdoped material.

FIG. 3 :

 3 FIG.3: Microwave absorption, measured using the TM010 mode ( 19.2 GHz ) of one of the OFHC copper cavities, of three of the underdoped Bi2Sr2CaCu2O 8+δ crystals before and after irradiation with 5.3 × 10 17 , 3 × 10 18 , and 7.7 × 10 18 2.3 MeV electrons cm -2 , respectively.

FIG. 4 :FIG. 5 :FIG. 6 :

 456 FIG.4: Variation of the critical temperature with electron fluence for single-crystalline underdoped Bi2Sr2CaCu2O 8+δ (•). A comparison with the Tc of a series of optimally doped control crystals ( ) shows that the onset temperature of magnetic screening (•) corresponds to flux exclusion by a surface region containing optimally doped material.

FIG. 7 :

 7 FIG. 7: Temperature dependence of the JPR frequency of electron-irradiated underdoped Bi2Sr2CaCu2O 8+δ .

Fig. ( 9 )FIG. 8 :

 98 FIG.8: f pl (0) versus Tc (open symbols), as well as f 2 pl (0) ∝ ρ c s (0) ∝ j c c (0) versus Tc (closed symbols), for both the underdoped irradiated crystals (squares) and a set of crystals with different doping levels (p ≈ 0.13, 39 and p ≈ 0.1147 ).

FIG. 9 :

 9 FIG.9: Square of the JPR frequency, normalized to its low temperature extrapolation f pl (0), versus reduced temperature T /Tc. This plot is representative of the temperature dependence of the c-axis superfluid stiffness, as well as of the maximum c-axis Josephson current:f 2 pl (T )/f 2 pl (0) ∼ λ 2 c (0)/λ 2 c (T ) ∼ ρ c s (T )/ρ c s (0) ∼ j c c (T )/j c c(0). The drawn line is a fit to Eq. (3) with ∆0 = 2.5kB Tc.

  ) Paramekanti et al. estimate the phase difference between two points separated by a vector perpendicular to the superconducting layers as φ 2

  Sr 2 CaCu 2 O 8+δ .39

	34 Stud-
	ies on Bi 2 Sr 2 CaCu 2 O 8+δ are limited to electron irra-
	diation of the single crystalline optimally doped ma-
	terial, that show a linear decrease of T c with electron
	fluence. 36,37,38 The c-axis superfluid density in a under-
	doped pristine Bi 2 Sr 2 CaCu 2 O 8+δ single crystal has been
	previously studied by Gaifullin et al., who invoked the
	IAH model to explain the much stronger temperature de-
	pendence of ρ c s in underdoped with respect to optimally doped Bi 2

  observed in slightly overdoped Bi 2 Sr 2 CaCu 2 O 8+δ 39,50 and optimally doped HgBa 2 Ce n-1 Cu n O 2n+2+δ , 51 but not in the present data on underdoped Bi 2 Sr 2 CaCu 2 O 8+δ . The modification of t ⊥ arising from the body-centered Bravais lattice of the Bi 2 Sr 2 CaCu 2

Acknowledgements

This work was supported in part by the French-Japanese bilateral program SAKURA, Grant No. 122313UL. C.J. van der Beek wishes to thank the Department of Physics of Kyoto University, where the surface resistance measurements were performed, for its hospitality. We thank A.E. Koshelev for useful discussions and a thorough reading of the manuscript. One of us (P.G.) was partially supported by MNiSW Grant No. N202 058 32/1202.


Yann Brenier

Optimal Transport, Convection, Magnetic Relaxation and Generalized Boussinesq equations

We establish a connection between Optimal Transport Theory (see [Vi] for instance) and classical Convection Theory for geophysical flows [Pe]. Our starting point is the model designed few years ago by Angenent, Haker and Tannenbaum [AHT] to solve some Optimal Transport problems. This model can be seen as a generalization of the Darcy-Boussinesq equations, which is a degenerate version of the Navier-Stokes-Boussinesq (NSB) equations. In a unified framework, we relate different variants of the NSB equations (in particular what we call the generalized Hydrostatic-Boussinesq equations) to various models involving Optimal Transport (and the related Monge-Ampère equation [Br, Ca]). This includes the 2D semi-geostrophic equations [Ho, CNP, BB, CGP, Lo] and some fully non-linear versions of the so-called high-field limit of the Vlasov-Poisson system [NPS] and of the Keller-Segel for Chemotaxis [KS, JL, CMPS]. Mathematically speaking, we establish some existence theorems for local smooth, global smooth or global weak solutions of the different models. We also justify that the inertia terms can be rigorously neglected under appropriate scaling assumptions in the Generalized Navier-Stokes-Boussinesq

equations. Finally, we show how a "stringy" generalization of the AHT model can be related to the magnetic relaxation model studied by Arnold and Moffatt to obtain stationary solutions of the Euler equations with prescribed topology (see [START_REF] Arnold | Topological methods in hydrodynamics[END_REF][START_REF] Moffatt | Magnetostatic equilibria and analogous Euler flows of arbitrarily complex topology. I. Fundamentals[END_REF][START_REF] Moffatt | Relaxation under topological constraints, Topological aspects of the dynamics of fluids and plasmas[END_REF][START_REF] Schonbek | Decay of Solutions to non-oscillating Magneto Hydrodynamics equations[END_REF][START_REF] Vladimirov | On general transformations and variational principles for the magnetohydrodynamics of ideal fluids. IV[END_REF][START_REF] Nishiyama | Magnetohydrodynamic approaches to measure-valued solutions of the two-dimensional stationary Euler equations[END_REF]).

1 The Angenent-Haker-Tannenbaum model for Optimal Transport problems

In this section, we consider the model introduced by Angenent, Haker and Tannenbaum [AHT]. This model was designed in order to seek the solutions of some optimal transport problems as equilibrium states of a suitable dynamical system that could be efficiently solved on a computer. The concrete applications have been computer vision, image registration and image warping.

Optimal transport and rearrangements

Let us briefly recall some typical results in Optimal Transport Theory, such as the polar factorization of maps. More precisely, let D be the closure of a bounded connected open set in R d , with a boundary of zero d-dimensional Lebesgue measure. Up to a rescaling, we assume the Lebesgue measure of D to be 1. Given an L 2 map y : D → R d , we call image measure of the Lebesgue measure on D by y the unique nonnegative (Borel) measure µ defined by:

R d f (x)µ(dx) = D f (y(a))da, (1)
for all compactly supported continuous functions f on R d . We have

R d µ(dx) = 1, R d |x| 2 µ(dx) = D |y(a)| 2 da,
which means that µ belongs to the set P rob 2 (R d ) of all (Borel) probability measures µ on R d such that |x| 2 µ(dx) < ∞. In this space, we say that a sequence µ n converges tightly to µ in P rob 2 (R d ), if:

R d f (x)µ n (dx) → R d f (x)µ(dx)
for all continuous function f on R d such that sup

x∈R d |f (x)| 1 + |x| 2 < +∞.
Given two L 2 maps y and z from D to R d , we say that they are rearrangement of each other if they define the same image measure. When y is a rearrangement of the identity map, we say, in short, that y is Lebesgue measure preserving.

Next, we define the class of maps with convex potential:

Definition 1.1 We say that an L 2 map from D to R d belongs to the class C of maps with a convex potential if there is a lower continuous convex function p : R d →] -∞, +∞] such that, for Lebesgue almost every point x ∈ D, p is differentiable at x and its gradient ∇p(x) coincides with y(x).

Then, we get from [Br]:

Theorem 1.2 (Rearrangements with convex potentials) For each L 2 map y : D → R d there is a unique rearrangement map with a convex potential y * ∈ C. The map y * depends on y only through the associated measure µ defined by (1).

In addition, the nonlinear operator y ∈ L 2 → y * ∈ L 2 is continuous as well as the induced operator µ ∈ P rob 2 (R d ) → y * ∈ L 2 , with respect to the tight convergence.

We get more precise results if y is a non degenerate map, in the sense that the pre-image of every Lebesgue negligible set is also negligible: Theorem 1.3 (Polar factorization of maps [Br]) Let y be a non degenerate L 2 map from D to R d . Then, there is a unique "polar factorization" y = Y • X where Y belongs to C and X is a Lebesgue measure preserving map of D. In this decomposition, Y is the unique rearrangement y * of y in C and X is the unique measure preserving map of D that minimizes D |X(a) -y(a)| 2 da. In addition, X can be written: X(a) = (∇Φ)(y(a)), a.e. a ∈ D,

where Φ is a convex Lipschitz function defined on R d .

For (much) more results on optimal transport, we refer to Villani's textbook [Vi]. The expression "optimal transport" comes from the fact that y * , among all possible rearrangements y of y 0 , is the unique minimizer of the "transportation cost".

D |y(x) -x| 2 dx,
where | • | denotes the Euclidean norm. The name "map with convex potential" is due to Caffarelli [Ca]. The concept of polar factorization has been extended to Riemannian manifolds by McCann [Mc]. Examples of concrete applications of optimal transport techniques to natural and computer sciences can be found in [FMMS, HZTA].

The AHT model

The AHT model is an attempt to get the unique rearrangement y * of y 0 , with convex potential, as the equilibrium state at t = +∞ of the following set of evolution equations:

∂ t y + (v • ∇)y = 0, ( 3
)
Kv + ∇p = y, ∇.v = 0, (4)
where

y = y(t, x) ∈ R d , v = v(t, x) ∈ R d , p = p(t,
x) ∈ R depend on t ≥ 0 and x ∈ D, and K is a "dissipative" operator to be chosen, for instance K = I or K = -∆. In these "AHT" equations, we denote the inner product in R d by • and we use notations:

∇ i = ∂ ∂x i , v • ∇ = j=1,d v j ∂ ∂x j , ∆ = j=1,d ∂ 2 ∂x 2 j .
The boundary conditions for the AHT system (3,4) are: i) the initial value of y at t = 0, y(0, x) = y 0 (x), ii) v is parallel to the boundary ∂D if K = I and v = 0 along the boundary if K = -∆. Notice that neither p nor v need initial conditions. As a matter of fact, as K = I, the second AHT equation (4) just corresponds to the "Helmholz decomposition" of y as a sum of a gradient field and a divergence-free field parallel to the boundary ∂D. So, we can write: v = P y, where P is a linear singular integral operator bounded in all L p space for 1 < p < +∞, provided that the domain D is smooth enough. In the case K = -∆, in a similar way we can write v = P ∆ y, where P ∆ is a linear singular integral operator bounded from L p to the Sobolev space W 2,p , for all 1 < p < +∞, D being assumed to be smooth. Thus we can write the AHT system (3,4) in a more abstract form:

∂ t y + (P K y • ∇)y = 0, (5)
with

P K = P if K = I and P K = P ∆ if K = -∆.

Expected long time behaviour of the AHT model

Let us now explain why the AHT model is expected to solve the Optimal Transport (or rearrangement) problem, at least for a large class of data. First, we observe that equation ( 3) expresses, at least formally, that, at each time t, y(t, •) is a rearrangement of y 0 . Indeed, for any smooth compactly supported function f , we get:

d dt D f (y(t, x)) dx = - D (∇f )(y) • (v • ∇)y dx = - D v • ∇[f (y)] dx,
(using the chain rule) which is zero, since v is divergence free and parallel to ∂D and is therefore L 2 orthogonal to any gradient field. (Notice that this calculation can be made rigorous provided that v has enough regularity.

According to Ambrosio's recent improvement of the DiPerna-Lions theory on ODEs [Am, DL], it is enough that v belongs to L 1 loc (R + , BV (D, R d )).) Next, we get the following balance law for the AHT

d dt D 1 2 |y(t, x) -x| 2 dx = - D (v • Kv)(t, x) dx, (6)
which implies, at least formally, that y and v respectively belong to the functional spaces

L ∞ (R + , L 2 (D, R d )) and L 2 (R + , H K (D)).
Here H K (D) denotes the Hilbert space of all divergence free fields w(x) ∈ R d for which D w•Kw dx is finite with suitable boundary conditions (w parallel to ∂D as K = I or w = 0 on ∂D as K = -∆). The formal proof of ( 6) is as follows:

d dt D 1 2 |y -x| 2 dx = - D (y -x) • ((v • ∇)y) dx (using the first AHT equation (3)) = - D v • ∇( 1 2 |y -x| 2 ) dx - D v • (y -x) dx = - D v • (y -x) dx
(since v is divergence free and parallel to ∂D and is therefore L 2 orthogonal to any gradient field)

= - D v • (Kv + ∇(p - 1 2 |x| 2 )) dx (using the second AHT equation (4)) = - D v • Kv dx.
At this point, we can describe the expected long time behaviour of the AHT system through the following heuristics. Since Kv •v is space-time integrable, we first argue that, as t → +∞, v presumably tends to zero. Then we expect y to have a definite strong limit y ∞ in L 2 , which is, then, necessarily a rearrangement of y 0 . Passing to the limit in the second AHT equation (4), we conclude that y ∞ must be a gradient. Therefore, at the end of the process, y 0 has been rearranged as a map y ∞ with a potential. Observe that this potential needs not being convex. This is obvious in the special case when y 0 is itself a map with a potential which is not convex. Indeed, then y(t, x) = y 0 (x), v(t, x) = 0, is a trivial stationary solution to the AHT equations (3,4) and we get y ∞ = y 0 as a map with a non convex potential. So, we need further assumptions on y 0 to be convinced that y ∞ has a chance to have a convex potential. A natural assumption is that y 0 is smooth with a positive jacobian determinant valued in some interval [r, 1/r] with 0 < r < +∞. Indeed, for such an initial condition, the AHT equations have a global smooth solution y (at least in the case when K = -∆, as discussed in the next subsection), with a jacobian determinant that must stay in the same interval [r, 1/r], since y is always a rearrangement of y 0 . So if the convergence to y ∞ , as t → +∞, is strong enough, we expect y ∞ to have a convex potential and, therefore, coincide with the unique rearrangement of y 0 with convex potential y * provided by Theorem 1.2. The results obtained in [AHT] are only partial and leave as an open question this issue.

Wellposedness of the AHT equations

From the PDE viewpoint, it is crucial to check that the AHT system (3,4) is wellposed, which is done by Angenent, Haker and Tannenbaum in [AHT], for a class of dissipative operator K including K = I. Let us briefly discuss the wellposedness issue in the cases K = I and K = -∆.

For K = I, the AHT system is similar to the inviscid Burgers equation,

∂ t y + (y • ∇)y = 0, (7)
since P K behaves like a pseudo-differential operator of order zero. Thus, the local in time existence of smooth solutions for smooth initial conditions can be obtained from rather standard energy estimates. It is a challenging and interesting open question whether the Lipschitz norm, in space, of such solutions may blow up in finite time (as it would be the case of the inviscid Burgers equation). In sharp contrast, in the case K = -∆, smooth solutions are clearly global in time. Indeed, from (3), we immediately get that |y(t, x)| is uniformly bounded by the sup norm of y 0 that we denote by M 0 and suppose, here, to be finite. Thus, because of (4), according to standard elliptic regularity theory, the L ∞ t (W 2,p x ) norm of v(t, x) is controled by M 0 for all finite p. Thus, the same is true for the sup norm of ∇v(t, x). Differentiating (3) in x, we deduce that the sup norm of ∇y(t, x) in x cannot grow, in sup norm, faster than exponentially in t as soon as ∇y 0 has a finite sup norm. So, there is no possible blow up of the Lipschitz norm of both v and y and, therefore, by a standard argument, smooth solutions must be global in time. (Notice that the dissipative operator K = (-∆) 1/2 , with appropriate boundary condition, would be borderline to get such a Lipschitz estimate. In that case v would be a priori only Log-Lipschitz, just like the Yudovich solutions of the 2D Euler equations [MP].) For more details, we refer to the paper by Angenent, Haker and Tannenbaum [AHT], where different kinds of operator K are considered.

Interpretation of the AHT system in terms of Convection Theory

From a Fluid Mechanics viewpoint, the AHT equations look very similar to the Boussinesq equations for convective flows, in particular to their Darcy-Boussinesq version. A classical model for Convection Theory is provided by the Navier-Stokes Boussinesq (NSB) equations that we are now going to review with more details. Using the Boussinesq approximation, the Navier-Stokes equations for an inhomogeneous incompressible fluid subject to gravity along the x d direction read:

ρ 0 (∂ t v + (v • ∇)v) -ν∆v + ∇p = y, ∇ • v = 0, (8)
∂ t y + (v • ∇)y = 0. (9)
Here, v = v(t, x) ∈ R d is the velocity field, p = p(t, x) ∈ R the pressure field, ρ 0 > 0 the average density of the fluid, ν the (constant) viscosity of the fluid, while y has only one component in the "vertical" direction x d , which is -g θ(t, x), where g is the gravity constant and θ(t, x) is the difference between the density of the fluid at (t, x) and the averaged density ρ 0 of the fluid.

(Usually in Convection Theory, a diffusion term is added to the advection equation for y [ID].) We recall that the Boussinesq approximation amounts to consider a variable density incompressible fluid for which the density variations are sufficiently small to be neglected in all terms except the gravity force. This approximation is widely used for ocean and atmosphere modelling [Pe]. (To the best of our knowledge the justification of this approximation is still an open problem in mathematical Fluid Dynamics, mostly because of our rather poor knowledge of the Navier-Stokes equations for inhomogeneous flows, see discussions in [Li, Ma] for instance.) By neglecting the inertia term, or equivalently by setting ρ 0 = 0 in the NSB (Navier-Stokes-Boussinesq) equations, we get the simpler Stokes-Boussinesq (SB) (related to large-Prandtl-number Convection Theory as in [DOR, Wa], for instance).

If, in addition, the diffusion term -ν∆v is replaced by a friction term such as v, we get the Darcy-Boussinesq (DB) model. We immediately see that both the SB and the DB equations are just particular cases of the AHT model (3,4), for which the vector valued function y has only one component along the x d axis. Indeed, the DB and SB models then respectively correspond to the choice K = I and K = -∆ in the second AHT equation ( 4). According to the discussion made in subsection 1.3, we expect, for the AHT model, the y(t, x) to converge, as t → +∞, to a map y ∞ (x) with, hopefully, a convex potential. In the particular case of the convective DB and SB models, y(t, x) has only one component in the x d direction, namely -g θ(t, x). Interestingly enough, the convergence toward a map with convex potential, exactly means, for the DB and SB models, that the density field tends to a density "profile" ρ ∞ (x d ), depending only on the vertical coordinate x d , and monotonically decreasing. This clearly corresponds, in terms of Convection Theory, to a "stable hydrostatic equilibrium". Notice that a similar discussion can be found in Moffatt's paper [Mo] (section 2) as a prelude to his Magnetic Relaxation model that we will consider at the end of the present paper.

2 Generalized Navier-Stokes-Boussinesq equations

The interpretation of the AHT model in terms of Convection Theory suggests the following "GNSB" generalization of the NSB (Navier-Stokes-Boussinesq) equations:

ǫ(∂ t v + (v • ∇)v) + Kv + ∇p = F (x, y), ∇ • v = 0, ( 10
)
∂ t y + (v • ∇)y = G(x, y), (11)
where y = y(t, x) ∈ R m is a vector-valued function (m ≥ 1, in practice m = d or m = 2d for the models discussed below), F and G are given smooth functions with bounded derivatives up to second order, respectively defined on R m and R m × R d , ǫ > 0 is a scaling factor introduced to single out the inertia term, and K is a linear dissipative operator. Depending on the applications in view, only the following cases will be considered: K = 0 (no dissipation), Kv = v (linear friction), Kv = -∆v (viscosity).

Existence theory for the GNSB equations

For simplicity, we consider in this subsection the domain D to be the unit periodic cube T d = R d /Z d , in order to avoid technicalities due to spatial boundary conditions. For the three possible choices of the dissipative operator

K Kv = 0, Kv = v, Kv = -∆v, (12)
the existence and uniqueness of a local in time smooth solution (y, v) of the GNSB equations (10,11), for each smooth initial initial condition (y 0 , v 0 ) given on the torus T d , follow from standard theory on Euler and Navier-Stokes equations (for which we refer to [Li]). We say that (y, v) is a weak solution if: 1) y(t, x) and v(t, x) depends continuously on t with values in L 2 (D, R d ) (with respect to the weak topology of L 2 ); 2) For all smooth time dependent vector fields w(t, x), z(t, x), with ∇•w = 0, we have:

d dt v • w dx = [ǫv • (∂ t + (v • ∇))w -Kv • w + F (x, y) • w]dx, (13)
d dt y • z dx = y • (∂ t + (v • ∇))z + G(x, y) • z)dx, (14)
3) The following energy inequality holds true:

1 2 d dt ( ǫ|v| 2 + |y| 2 )dx + Kv • v dx ≤ [F (x, y) • v + G(x, y) • y] dx. ( 15
)
When K = -∆, the existence of global weak solutions for the GNSB equations follows from standard arguments à la Leray combined with the DiPerna-Lions theory on ODEs [Li, DL]. They are unique in 2 space dimensions. In sharp contrast, as K = I or K = 0, nothing can be said about global weak solutions.

Concerning global smooth solutions, the existence theory is quite challenging, even in 2 space dimensions. Recently, Chae, Hou and Li [Ch, HL] have proven that Navier-Stokes Boussinesq equations (8,9) (just called Boussinesq equations in these papers) have global smooth solutions when d = 2 and K = -∆.

The same result can be readily extended to the GNSB equations (10,11), essentially because we assume the right-hand sides F (x, y) and G(x, y) of each equation to be smooth functions of y and v with bounded derivatives up to order two. (Indeed, these assumptions are enough for a straghtforward adaptation of the proof of Theorem 1.1 in Chae's paper, through estimates (2.1 • • • 17) in [Ch]. Some constants involved in these estimates have just to be modified to take into account the Lipschitz constants of F and G.) So, we can summarize all these results in the following Theorem, which is nothing but a straightforward adaptation of known results:

Theorem 2.1 Assume the dissipative operator K to be of type (12). Then, the generalized Navier-Stokes Boussinesq equations (10,11) admit, for any smooth initial condition, a unique local smooth solution.

If K = -∆, the GNSB equations admit at least a global weak solution (y, v) (in the sense (13,14,15)) for any initial condition (y 0 , v 0 ) in L 2 . If d = 2, these weak solutions are unique. Furthermore, still for d = 2, the solutions are globally smooth for smooth initial conditions.

Zero inertia limit of the GNSB equations

By zero inertia limit of the GNSB, we mean the formal limit obtained by dropping the scaling factor ǫ in front of the inertia terms in (10,11). Namely, in Eulerian coordinates,

∂ t y + (v • ∇)y = G(x, y), (16)
Kv + ∇p = F (x, y), ∇ • v = 0.
We are able to make a rigorous derivation of the zero inertia limit when K is strictly dissipative (K = 0 being so excluded):

Theorem 2.2 Assume the dissipation operator K to be coercive in L 2 , namely K ≥ α for some constant α > 0. Then the zero inertia equations ( 16) admit, for any smooth initial condition, a local smooth solution, which is global if d = 2 and K = -∆. This solution can be obtained as the limit, as ǫ goes to zero, of the weak solutions of the GNSB equations (10,11), with the same initial condition.

For the convergence, we use a simple energy method. Namely, given a weak solution (y ′ , v ′ ) to the GNSB equations (13,14,15) and a solution (y, v) of the HF equations, with same initial conditions (y 0 , v 0 ), we introduce

e(t) = T d ( ǫ|v ′ | 2 + |y -y ′ | 2 ) dx (17)
and try to get an estimate of form:

d dt (e(t) + O(ǫ)) + 1 2 T d K(v -v ′ ) • (v -v ′ ) dx ≤ (e(t) + O(ǫ))c, (18)
where c depends only on the limit (y, v), for any fixed finite time interval [0, T ] on which (y, v) is smooth. From this estimate (18), we immediately get that y -y ′ and v -v ′ are of order O( √ ǫ) in, respectively, L ∞ ([0, T ], L 2 (T d )) and L 2 ([0, T ], L 2 (T d )), using the coercivity of K (K ≥ α for some α > 0). So, we are left with proving (18). Notice first that, from equations ( 15) and ( 16), the following energy balances hold true:

1 2 d dt ( ǫ|v ′ | 2 + |y ′ | 2 )dx + Kv ′ • v ′ dx ≤ [F (x, y ′ ) • v ′ + G(x, y ′ ) • y ′ ] dx 1 2 d dt |y| 2 dx + Kv • v dx = [F (x, y) • v + G(x, y) • y] dx.
Since (y ′ , v ′ ) and (y, v) are respectively supposed to be a weak solution of the GNSB equations and a smooth solution of the zero inertia limit ( 16), we also get

- d dt y • y ′ dx = [((v • ∇)y) • y ′ -((v ′ • ∇)y) • y ′ ]dx -[G(x, y) • y ′ + G(x, y ′ ) • y]dx = ((v • ∇)y -(v ′ • ∇)y) • (y ′ -y) dx + [(G(x, y) -G(x, y ′ )) • (y -y ′ ) -G(x, y ′ ) • y ′ -G(x, y) • y]dx (using that (w • ∇)y) • y) dx = 0 for both w = v and w = v ′ ) ≤ c 1 (|v -v ′ ||y -y ′ | + |y -y ′ | 2 ) dx -[G(x, y ′ ) • y ′ + G(x, y) • y]dx
where c 1 depends on the Lipschitz constants of G (as a function of y and v) and v (as a function of x). Thus, adding up these three equalities, we get by definition ( 17):

d dt e(t) + (Kv ′ • v ′ + Kv • v) dx ≤ c 1 (|v -v ′ ||y -y ′ | + |y -y ′ | 2 ) dx + [F (x, y ′ ) • v ′ + F (x, y) • v] dx.
From ( 13) and ( 16), we also get:

F (x, y) • v ′ dx = (Kv + ∇p) • v ′ dx = Kv • v ′ dx and F (x, y ′ ) • v dx = ǫ d dt v ′ • v dx -ǫ v ′ • (∂ t + v ′ • ∇)v dx + Kv ′ • v dx = √ ǫ( d dt r 1 (t) + r 2 (t)) + e 2 (t) + Kv ′ • v dx, where r 1 (t) = √ ǫ v ′ • v dx, r 2 (t) = - √ ǫ v ′ • ∂ t v dx and e 2 (t) = -ǫ v ′ • (v ′ • ∇)v dx.
Notice that r 2 1 , r 2 2 and |e 2 | are bounded by c 2 e(t) (by definition ( 17)), where c 2 depends on the Lipschitz constant of v as a function of both t and x. So, we have obtained:

d dt (e(t) - √ ǫr 1 (t)) + K(v ′ -v) • (v ′ -v) dx ≤ c 1 (|v -v ′ ||y -y ′ | + |y -y ′ | 2 ) dx + √ ǫr 2 (t) + c 2 e(t).
Using the coercivity of K (K ≥ α > 0) and definition (17), we find c 3 depending on the Lipschitz constants of F , G and v such that:

c 1 (|v -v ′ ||y -y ′ | + |y -y ′ | 2 ) dx ≤ c 3 e(t) + 1 2 K(v ′ -v) • (v ′ -v) dx.
This leads to:

d dt (e(t) - √ ǫr 1 (t)) + 1 2 K(v ′ -v) • (v ′ -v) dx ≤ (c 2 + c 3 )e(t) + √ ǫr 2 (t),
which leads to a differential inequality of the desired type, namely (18), since r 2 1 + r 2 2 ≤ 2c 2 e. Thus, the proof of Theorem 2.2 is now achieved.

A Mechanical interpretation of the GNSB equations

In this subsection, we provide a mechanical interpretation of the GNSB equations and their zero inertia limit. For this purpose, it is worth considering the GNSB equations (10,11) in "Lagrangian coordinates". Assuming the vector field v to be smooth enough, denoting by a ∈ D the position of a fluid parcel at time t = 0, we can recover its position X(t, a) ∈ D at later time t ≥ 0 by solving the ODE ∂ t X(t, a) = v(t, X(t, a)), X(0, a) = a, ∀a ∈ D.

Notice that, for each t, a ∈ D → X(t, a) ∈ D is a measure preserving map as a consequence of the fact that v is a smooth divergence-free vector field parallel to ∂D. Let us also introduce:

Y (t, a) = y(t, X(t, a)) ∈ R m . (20)
Then, the GNSB equations (10,11) read in Lagrangian coordinates:

ǫ ∂ tt X(t, a) + (Kv)(t, X(t, a)) + (∇p)(t, X(t, a)) = F (X(t, a), Y (t, a)), (21)
∂ t Y (t, a) = G(X(t, a), Y (t, a)),
where a ∈ D → X(t, a) ∈ D is Lebesgue measure preserving. Let us now provide a possible mechanical interpretation. We model the atmosphere (or the ocean) as a continuous distribution of infinitesimal rigid balloons floating inside D, each of them having position X(t, a) at time t, with X(0, a) = a, and being attached with probability λ(a) ≥ 0 to an an anchor with an elastic cable. Of course, to be a realistic model with real balloons, λ should be a discrete probability distribution concentrated on a finite collection of points and the corresponding balloons should have a finite extension! Let us rather assume, for mathematical simplicity, that the balloons are just points and that λ is a smooth nonnegative density function on D with unit mass. The cable corresponding to the balloon labelled by a is modelled by a (possibly non Hookean) spring with restoring force -k(ξ, a) = k(-ξ, a) ∈ R d where ξ ∈ R d is the elongation of the spring. Notice that k may depend on a. The location of the anchor attached to the balloon labelled by a is not necessarily fixed and denoted by Y (t, a) ∈ R d . (We may also think of an aircraft, or a boat, or any kind of carrier instead of an anchor.) Notice that we do not require the anchor to be located in D. Neglecting any interaction between the fluid and both the anchors and the springs (which may not be very realistic), we obtain the following dynamical equation for each balloon

ǫ∂ tt X(t, a) + (Kv)(t, X(t, a)) + (∇p)(t, X(t, a)) = (22) -λ(a)k(X(t, a) -Y (t, a), a).
(Observe that as λ(a) = 0, the corresponding carrier Y (t, a) is just fictitious!) Let us consider the special case when the speed of each anchor is constant and given by:

∂ t Y (t, a) = W (a) ∈ R 2 . ( 23
)
Implicitly define a field y

= (ỹ, ŷ) = y(t, x) ∈ R d × D by setting ỹ(t, X(t, a)) = Y (t, a), ŷ(t, X(t, a)) = a (remember that a → X(t, a) is supposed to be a diffeomorphism). Noticing that ((∂ t + v • ∇)y)(t, Y (t, a)) = ∂ t [y(t, X(t, a))] = (∂ t Y (t, a), 0) = (W (a), 0) = (W (ŷ(t, X(t, a))), 0)
and going back to Eulerian coordinates, we recover the GNSB equations (10,11) in the particular case:

F (x, y) = -λ(ŷ)k(x -ỹ, ỹ), G(x, y) = (W (ŷ), 0). ( 24
)
(Notice that assuming λ and the restoring force k to have bounded derivatives up to order 2 is not very realistic! These assumptions are clearly made for mathematical convenience.) We may also consider the following variant of this mechanical model. Instead of prescribing their velocity by ( 23), we may assume that the carriers are driven by a friction-dominated retroaction of type:

η∂ tt Y (t, a) + ∂ t Y (t, a) = -µ(a)k(Y (t, a) -X(t, a), a),
where µ is a given nonnegative function. Dropping the inertia term (η = 0) leads to the following law:

∂ t Y (t, a) = -µ(a)k(Y (t, a) -X(t, a), a), (25)
In that case, we get (keeping unchanged the definitions of Y and y) again the GNSB equations (10,11) with F and G given by:

F (x, y) = -λ(ŷ)k(x -ỹ, ŷ), G(x, y) = (-µ(ŷ)k(ỹ -x, ŷ), 0). ( 26
)
Let us finally consider a second variant where the Coriolis force is added to the model (rotating ocean or atmosphere). Neglecting the vertical extension, so that d = 2, and assuming the rotation vector to be perpendicular to the ocean and of unit length, the Coriolis force Jv = (-v 2 , v 1 ) is completely absorbed by the pressure term and the fluid parcels are not sensitive to it. (Indeed, v = (v 1 , v 2 ) being divergence free, Jv is a gradient and can be removed from the dynamical equation.) However we may think that the carriers are still sensitive to the Coriolis force. Thus, we get for them

η∂ tt Y (t, a) + J∂ t Y (t, a) = -µ(a)k(Y (t, a) -X(t, a), a),
instead of ( 25), neglecting a possible friction term. Neglecting the inertia term (η = 0) leads to the balance equation:

∂ t Y (t, a) = Jµ(a)k(Y (t, a) -X(t, a), a) (27)
(using that J 2 = -J). So we end up with a third version of the GNSB equations (10,11), for which:

F (x, y) = -λ(ŷ)k(x -ỹ, ŷ), G(x, y) = (Jµ(ŷ)k(ỹ -x, ŷ), 0). ( 28
)
To summarize this subsection, let us just say that the GNSB equations (10,11) provide a rather flexible framework to describe the interaction between an incompressible fluid confined in a d-dimensional domain D (each fluid parcel being labelled by its initial position a ∈ D having position X(t, a) at time t), and a set of particles (also labelled by a and of position Y (t, a)) moving in the ambient space R d . The unusual feature of the resulting models is that the interaction is pairwise (X(t, a) interacts only with Y (t, a) for the same label a), except for the mediation by the pressure field p(t, x) wich preserves a → X(t, a) in the class of Lebesgue measure preserving maps of D at each time t.

The Generalized Hydrostatic Boussinesq equations

In this final section, we investigate the most degenerate version of the GNSB equations (10,11), where we neglect not only the inertia terms but also the dissipative operator K. Thus we are left with the strange looking system:

F (x, y) = ∇p, ∇ • v = 0, ( 29
)
∂ t y + (v • ∇)y = G(x, y), (30)
that we call Generalized Hydrostatic Boussinesq (GHB) equations (by seing (29) as a generalization of the hydrostatic balance in Convection Theory).

Let us concentrate on the simpler case when m = d and

F (x, y) = y -x (31)
(which corresponds to cables modelled by Hookean springs, according to the mechanical interpretation of subsection 2.3). Thus, (29) just reads:

y = x + ∇p(t, x). (32)
In Lagrangian coordinates, the Generalized Hydrostatic Boussinesq equations (30,32) become:

Y (t, a) = X(t, a) + (∇p)(t, X(t, a)), (33)
∂ t Y (t, a) = G(X(t, a), Y (t, a)), (34)
where, for all t, a ∈ D → X(t, a) ∈ D is a measure preserving map.

3.1 Formal derivation of some optimal transport models from the GHB equations

We claim that several models involving optimal transport and the Monge-Ampère equation correspond to these GHB equations. In particular, we consider the following generalization of the "semigeostrophic (SG) equations" [Ho, CNP, BB, CGP, Lo]:

∂ t ρ + ∇ • (ρw) = 0, (35)
w(t, x) = B∇ϕ, (36)
det(I + D 2 ϕ) = ρ, (37)
where B is a d × d constant matrix and D 2 ϕ is the "Hessian" matrix, made of all second order derivatives of ϕ(t, x) with respect to x. This system, that we call Generalized Semi-Geostrophic (GSG) equations involves the "fully nonlinear" Monge-Ampère (MA) equation ( 37) which, requires, in order to be of elliptic type, the convexity condition:

I + D 2 ϕ(t, x) > 0, (38)
in the sense of symmetric matrices, for each time t. The 2D SG equations [Ho, CNP, BB, CGP, Lo] just correspond to the special case when d = 2 and B is the rotation matrix of angle π/2. The case when B is just a number can be related to drift-diffusion and Keller-Segel type models [CMPS] for which the MA equation is replaced by the linear Poisson equation

∆ϕ = ρ -1, (39)
which can be seen as a linear approximation of the MA equation (37) as ϕ is small. The drift-diffusion case corresponds to (35,36,39) with B > 0.

The simplified version of the Keller-Segel model [KS] treated by Jäger and Luckhaus [JL] corresponds to B < 0 (with an additional diffusion term for ρ in equation ( 35)).

Let us now show that a solution of the GHB equations (33,34) corresponds to a solution of the GSG equations (35,36,37,38), in a suitable sense. For this purpose, in order to use the Polar Factorization Theorem 1.3, we make the following a priori assumptions for each time t:

A1: The map Y (t, •) is non degenerate, A2: The map x ∈ D → x + ∇p(t, x) ∈ R d has a convex potential.
These assumptions mean that (33) defines the polar factorization of Y (t, •) where X(t, •) is measure preserving and x ∈ D → x + ∇p(t, x) has a convex potential. According to the Polar Factorization Theorem 1.3, the measure preserving factor X(t, •) can be written:

X(t, a) = (∇Φ)(t, Y (t, a)),
where Φ(t, x) is convex and Lipschitz continuous in x ∈ R d , or, equivalently,

X(t, a) = Y (t, a) + (∇ϕ)(t, Y (t, a)), (40)
where ϕ(t, x) = Φ(t, x)-|x| 2 /2. Since Y is supposed to be non degenerate (by assumption A1), there is a nonnegative Lebesgue integrable "density field" ρ(t, x) such that

R d f (x)ρ(t, x)dx = D f (Y (t, a))da, (41)
for all suitable functions f . Thus:

R d f (x + ∇ϕ(t, x))ρ(t, x)dx = R d f (Y (t, a) + (∇ϕ)(t, Y (t, a)))da (by definition (41) of ρ) = D f (X(t, a))da = D f (x)dx
(thanks to (40) and because X(t, •) is Lebesgue measure preserving). So, we have obtained

R d f (x + ∇ϕ(t, x))ρ(t, x)dx = D f (x)dx, (42)
for all compactly supported continuous function f . This, combined with the assumption that x + ∇p(t, x) has a convex potential can be seen as a weak form of the Monge-Ampère equation ( 37) combined with the ellipticity condition (38). Next, using (40), we can write (34) as

∂ t Y (t, a) = w(t, Y (t, a)),
where w is the vector field defined by:

w(t, x) = G(x + ∇ϕ(t, x), x), (43)
which is nothing but a generalization of equation ( 36). Next, we get for all smooth compactly supported function f on R d :

d dt D f (Y (t, a))da = D (∇f )(Y (t, a)) • w(t, Y (t, a))da,
which means, in terms of ρ defined by (41):

d dt R d f (x)ρ(t, x)dx = R d ∇f (x) • w(t, x)ρ(t, x)dx,
that is (35) in a weak sense. So, we have fully recovered the GSG system (35,36,37,38), in a suitable weak form, from the GHB equations (33,34) under Assumption A1,A2. In addition, equation ( 36) can be replaced by the more general relation (43).

A global existence theorem for the GHB equations

We are now going to introduce a suitable concept of solutions of (33,34), by assuming a priori that, in (33), for each time t, the map x ∈ D → x+∇p(t, x) has a convex potential and, therefore, is the unique rearrangement Y * (t, •) of Y (t, •) in the class C of map with a convex potential, as in Theorem 1.2. (This corresponds to assumption A2 in the previous subsection.) Therefore, we can just write (33) as:

Y (t, a) = Y * (t, X(t, a)).
From (34), we get (at least formally) that:

d dt D f (Y (t, a))da = D (∇f )(Y (t, a)) • G(X(t, a), Y (t, a))da.
for all compactly supported C 1 function f . Thus,

d dt D f (Y * (t, X(t, a)))da = = D (∇f )(Y * (t, X(t, a))) • G(X(t, a), Y * (t, X(t, a)))da.
Now, we can factor out X(t, a) and get a set of self-consistent equations for Y * , namely:

d dt D f (Y * (t, a))da = D (∇f )(Y * (t, a)) • G(a, Y * (t, a))da,
without loss of information for Y * . This suggests the following concept of solution for the GHB equations (33,34):

Definition 3.1 Assume that G is continuous and satisfies

sup x,y |G(x, y)| 1 + |y| < ∞. ( 44
)
We say that

Y * ∈ C 0 ([0, T ], L 2 (D, R d ))
is the "convex rearrangement" (CR) solution to the GHB equations (33,34), if: 1) Y * (t, •) belongs to the set C of all maps with convex potenial, for all t ∈ [0, T ],

2) For all compactly supported C 1 function f on R d , we have:

d dt D f (Y * (t, a))da = D (∇f )(Y * (t, a)) • G(a, Y * (t, a))da. ( 45
)
This concept yields the following global existence theorem (without uniqueness) for all initial conditions in L 2 : Theorem 3.2 For each initial condition Y 0 ∈ L 2 (D, R d ), there is at least one CR-solution Y * (t, a), in the sense of Definition 3.1, such that Y * (t = 0, •) = (Y 0 ) * . This solution can be obtained as the limit in C 0 t (L 2 a ) as h → 0 of a time discrete approximation Y h (t, a) defined, first at discrete times t = nh, by:

Y h (nh + h, a) = [Y h (nh, a) + h G(a, Y h (nh, a))] * , n = 0, 1, 2, • • • (46)
(where * denotes the rearrangement operator as in Theorem 1.2) and, then, linearly interpolated in t.

Proof

To get the existence result, it is enough to show the convergence of the time discrete approximation Y h . First, we observe that, Y h (t, •) is valued in C, (the class of maps with convex potential) for all time t. (This is true by definition for discrete times t = nh and preserved by linear interpolation since C is a convex cone.) Next, we deduce from (46) and assumption (44):

D |Y h (nh + h, a)| 2 da ≤ hc + (1 + hc) D |Y h (nh, a)| 2 da, (47)
for some constant c depending only on G and D. We also get, for all compactly supported C 1 function f :

D [f (Y h (nh + h, a)) -f (Y h (nh, a))]da = (48) = h D 1 0 (∇f )(Y h (nh, a) + hθG(a, Y h (nh, a))) • G(a, Y h (nh, a))dθda,
which can be bounded by:

hc sup x |∇f (x)| 1 + |x| D (1 + |Y h (nh, a)| 2 )da, (49)
where c depends only on D and G. So, from (47), we first see that Y h (t, •) is bounded in L 2 uniformly in t ∈ [0, T ] and h ∈]0, T ] by some constant R. Therefore Y h is uniformly valued in C R the set of maps with convex potential and L 2 norm bounded by R, which is a compact subset of L 2 . Next, we deduce from (48,49) that, for each fixed

C 1 function f such that |∇f (x)| 1 + |x| < +∞, t → D f (Y h (t, a))da
is Lipschitz continuous on [0, T ], uniformly in h. Since Theorem 1.2 asserts the continuity of the map µ → y * , we get that Y h (t, •) is uniformly equicontinuous from [0, T ] to L 2 . Then, we deduce from the Ascoli-Arzela theorem that the set of all time-discrete approximations Y h (t, a), for 0 < h ≤ T , is relatively compact in C 0 t (L 2 a ). Thus, there is a sequence of time steps h for which Y h converges to some limit Y * in C 0 t (L 2 a ), which is necessarily valued in C R , and therefore in C. We also easily get (45) by letting h go to zero in (48). So, the proof of Theorem 3.2 is now complete.

Remark: continuous dependence as d = 1

In the very special case d = 1, the rearrangement operator * is well known to be non expansive in L 2 :

D |y * (x) -z * (x)| 2 dx ≤ D |y(x) -z(x)| 2 dx,
for pair (y, z) of L 2 applications from D to R. It follows that two CR solutions y * and z * , obtained as limits of the time-discrete approximations (46), with respective initial condition y 0 and z 0 , must satisfy:

D |y * (t, x) -z * (t, x)| 2 dx ≤ exp(ct) D |y 0 (x) -z 0 (x)| 2 dx,
where c depends only on D and G.

Optimal transport and Magnetic Relaxation

In this section, we discuss a natural "stringy" generalization of the AHT equations (3,4) (discussed in the first part of the paper) and establish a link with the Arnold-Moffat model of Magnetic Relaxation (see [START_REF] Arnold | Topological methods in hydrodynamics[END_REF][START_REF] Moffatt | Magnetostatic equilibria and analogous Euler flows of arbitrarily complex topology. I. Fundamentals[END_REF][START_REF] Moffatt | Relaxation under topological constraints, Topological aspects of the dynamics of fluids and plasmas[END_REF][START_REF] Schonbek | Decay of Solutions to non-oscillating Magneto Hydrodynamics equations[END_REF][START_REF] Vladimirov | On general transformations and variational principles for the magnetohydrodynamics of ideal fluids. IV[END_REF][START_REF] Nishiyama | Magnetohydrodynamic approaches to measure-valued solutions of the two-dimensional stationary Euler equations[END_REF]).

The AHT model as a gradient flow

Using Lagrangian coordinates, we deduce from (3,4), in the case when the dissipation operator is K = 1:

∂ t X(t, a) + (∇p)(t, X(t, a)) = y 0 (a) -X(t, a) (50)
where X(t, •) belongs to the set MP M(D) of all (Borel) Lebesgue measure preserving maps of D. (These equations can be either derived directly from (3,4) or obtained from the GNSB equations written in Lagrangian coordinates (21), by setting ǫ = 0, K = I, F (x, y) = y -x and G(x, y) = 0.) As explained in [AHT], in slightly different words, the AHT equation ( 50) formally corresponds to the gradient flow of the "energy"

X → 1 2 D |X(a) -y 0 (a)| 2 da, (51)
on the "manifold" of all X ∈ MP M(D) for the L 2 metrics. Let us recall that, as seen in the first section, minimizing this energy is equivalent to solve an Optimal Transport problem. The gradient flow structure can be easily understood by considering the standard time discretization of such a gradient flow. Let h > 0 be a time step and let us denote by X h (t, a) the discrete approximation of X(t, a) at discrete time t = nh, n = 0, 1, 2, 3, • • •. At t = 0, we set X h (0, a) = a and, for t = nh, n = 1, 2, 3, • • •, we require X h (t, •) to be a minimizer among all X ∈ MP M(D) of the following functional:

a time discrete approximation X h (t, a, s), by setting X h (0, a, s) = a and asking, for t = nh, n = 1, 2, 3, • • •, X h (t, •, •) to be a minimizer among all curves s ∈ [0, 1] → X(•, s) ∈ MP M(D) of the functional: The formal optimality condition reads:

X h (nh, a, s) -X h (nh -h, a, s) h = = ∂ 2 ss X h (nh, a, s) + (∇p h )(nh, X h (nh, a, s), s), for some scalar function p h . So, we formally obtain, as h → 0, the desired equation ( 56). Equation ( 56) has an interesting interpretation, obtained by assuming a priori that a ∈ D → X(t, a, s) is a smooth orientation and measuring preserving diffeomorphism of D for each (t, s). Then we introduce, for each (t, s), two divergence free vector fields parallel to the boundary ∂D, namely v(t, x, s) ∈ R d and b(t, x, s) ∈ R d , defined by: v(t, X(t, a, s), s) = ∂ t X(t, a, s), b(t, X(t, a, s), s) = ∂ s X(t, a, s).

(58)

Then, we get from (56):

v = ∂ s b + (b • ∇)b + ∇p, (59)
while, from (58), we get the compatibility condition

∂ t b + (v • ∇)b = (b • ∇)v + ∂ s v (60)
(by writting ∂ 2 ts X = ∂ 2 st X), to be added to the divergence free constraints

∇ • v = ∇ • b = 0, v//∂D, b//∂D, (61)
and the boundary conditions at s = 0 and s = 1 induced by (54), namely:

v(t, x, s = 0) = v -(t, x), v(t, x, s = 1) = v + (t, x), (62)
where v + and v -are prescribed. When the fields v and b do not depend on s, we get the Magnetic Relaxation model discussed by Moffatt in [Mo] (see also [START_REF] Arnold | Topological methods in hydrodynamics[END_REF][START_REF] Moffatt | Relaxation under topological constraints, Topological aspects of the dynamics of fluids and plasmas[END_REF][START_REF] Schonbek | Decay of Solutions to non-oscillating Magneto Hydrodynamics equations[END_REF][START_REF] Vladimirov | On general transformations and variational principles for the magnetohydrodynamics of ideal fluids. IV[END_REF][START_REF] Nishiyama | Magnetohydrodynamic approaches to measure-valued solutions of the two-dimensional stationary Euler equations[END_REF]). As t → +∞, we expect, at least for a large class of initial conditions, the solution of equations (59,60,61)

Of course, we are far from being able to provide any rigorous proof of this conjecture.

The "cross-Burgers" equation

In the case when D is the unit ball, The Magnetic Relaxation equations (59,60,61) admit special solutions (b, v, ∇p) which are linear in x: b(t, x, s) = B(t, s)x, v(t, x, s) = V (t, s)x, ∇p(t, x, s) = G(t, s)x, (65)

where B, V are skew-symmetric matrices, while G is a symmetric matrix, all depending only on (t, s). (Notice that the fields b and v are automatically parallel to the boundary ∂D since D is the unit ball.) The resulting equations for B, V and G are:

V = ∂ s B + B 2 + G, (66)
∂ t B + [V, B] = ∂ s V. ( 67
)
Since B 2 is a symmetric matrix, equation (66) reduces to:

V = ∂ s B.
Thus, we get a single equation for B:

∂ t B + [∂ s B, B] = ∂ 2 ss B, (68)
where [A, B] denotes the skew product AB -BA. In the special case d = 3, B can be identified as a 3-vector and [•, •] as the cross product × in R 3 , which leads to:

∂ t B + ∂ s B × B = ∂ 2 ss B. ( 69
)
that we could call the "cross-Burgers" equation. This equation admits interesting special solutions, such as: B(t, s) = (α(t) cos s , α(t) sin s , β(t) -1)

where α ≥ 0 and β are solutions to:

dα dt = -βα, dβ dt = α 2 ,
or, equivalently, d 2 λ dt 2 + exp(2λ) = 0, where λ = log(α).

  The field p can be recovered by solving the Poisson problem: ∆p = ∇.y, inside D with inhomogeneous Neumann condition ∇p • n = y • n along the boundary, where n denotes the outward normal. Then, we get: v = y -∇p.

  -y 0 (a)| 2 da,

  to converge toward an equilibrium, for which v = 0 and b = b(x, s), p = p(x, s) are solutions to the Euler equations [AK, MP] (s acting as the time variable and b as the velocity field): ∂ s b + (b • ∇)b + ∇p = 0, (63) ∇ • b = 0, b//∂D.

Acknowledgment

This work originated while the author was visiting Universität Bonn (Institut für Angewandte Mathematik) and continued at Universität Wien (Wolfgang Pauli Institut). This work is also part of the research made in the LRC CEA-Cadarache/CNRS-Université de Nice and the ANR OTARIE project (ANR BLAN07-2-183172). We thank Adrien Blanchet, Marco Di Francesco, Francis Filbet and François Gallaire for fruitful discussions.

(after rearranging the squares). Thus, assuming a priori that X h (nh-h, •) + y 0 h is non degenerate and using Theorem 1.3, this exactly means:

where y * is a map with a convex potential. If we write this map as

we get

which can be seen just as a finite difference approximation of equation ( 50) as h → 0. This is enough to interpret, at least formally, equation ( 50) as the L 2 gradient flow of energy (51) on the "manifold" MP M(D).

A stringy generalization of the AHT model

This analysis suggests the possibility of more complex models based on similar ideas. A rather natural idea amounts to consider, instead of the "manifold" MP M(D), the "manifold" of "strings" valued in MP M(D): s ∈ [0, 1] → X(•, s) ∈ MP M(D), with fixed end values, say

Then, we may think of the L 2 gradient flow of the following string energy:

We claim that the resulting equation read:

where X(t, •, s) is valued in MP M(D) and end point conditions (54) are enforced. To get this system, as we did before for the AHT model, we define