Jeroen P A Devreese

Michiel Wouters

Jacques Tempere

Competition between the Fulde-Ferrell-Larkin-Ovchinnikov phase and the BCS phase in the presence of an optical potential

published or not. The documents may come

Competition between the FuldeFerrellLarkinOvchinnikov phase and the BardeenCooper Schrieffer phase in the presence of an optical potential

I. INTRODUCTION

In the past decennium, the field of cold atoms has witnessed the experimental realization of various new quantum coherent phenomena in ultracold gases [1]. Amongst the various systems that are being studied, fermionic superfluids occupy a prominent role [2][3][4]. Contrary to bosonic gases, fermionic gases must form pairs between fermions in different spin states in order to form an s-wave superfluid. A fundamental question that has attracted wide attention recently, is what happens when the ratio between these different spin-states is altered, thereby creating a spin-imbalanced or polarized superfluid. This effect is akin to applying a magnetic field to a superconductor, which leads to electrons aligning their spins. In superconductors, and in condensed matter systems in general, it is difficult to control various parameters such as the polarization or the interaction strength. In ultracold gases, however, one has an unprecedented degree of control over these parameters [6,7].

For instance, experimentalists can tune the interaction strength by making use of Feshbach resonances (for a review on this subject see Ref. [8]). Using this technique, the crossover from a Bardeen-Cooper-Schrieffer (BCS) superfluid to a Bose-Einstein condensate (BEC) of bound molecules has been studied [9][10][11][12]. Furthermore, polarization can be created by first loading particles into one hyperfine state and subsequently applying a radio-frequency sweep to send a controlled number of particles to a different hyperfine state. This experimental liberty has led to the realization of many experiments on the imbalanced ultracold Fermi gas [13][14][15][16][17]. Today, the main phases of this system, as a function of temperature, interaction strength and polarization, have been observed experimentally [19]. One of the most fundamental effects that occurs in a polarized Fermi gas is that above a critical imbalance, fermionic superfluidity ceases to exist and the system makes a transition into a normal Fermi gas. This is known as the Clogston-Chandrasekhar limit [18]. One important question that has emerged is whether there exist other pairing-mechanisms by which a gas of interacting fermions can accommodate polarized superfluidity. Recently, there has been an intensive theoretical search for these exotic superfluid phases, such as the breached pair or Sarma phase [20,21], and phase separation [22]. Another phase that has attracted wide attention is the Fulde-Ferrell-Larkin-Ovchinnikov phase (FFLO phase), that was proposed independently by Fulde and Ferrell [23] and by Larkin and Ovchinnikov [24] in 1964. Their idea was that a fermionic system can accommodate spin imbalance, by forming pairs with finite center-of-mass momentum. Up till now, this state has never been observed. A recent paper reports indirect experimental evidence for the FFLO state in one dimension (1D) [25], but in three dimensions (3D), the FFLO state has so far eluded experimental observation. This can be related to theoretical predictions, which state that the FFLO state only occupies a tiny sliver of the BCS-BEC crossover phase diagram [26,27]. It was therefore necessary to find a new method to stabilize the FFLO state in 3D. At the moment, new ideas are emerging to realize this goal. For instance, in two recent papers, it was suggested to stabilize the FFLO state by the use of a 3D cubic lattice [28,29]. As an alternative approach, we proposed to stabilize the FFLO state in an imbalanced 3D Fermi gas by adding a 1D optical potential to this system [30]. This turned out to enhance the stability region of the FFLO state by up to a factor six. In Ref. [30], the proof of principle for this concept was established using the path integral method , but questions remained about the exact effects of the properties of the optical potential on the FFLO state. In this paper, we investigate the role of the optical potential in detail. Our aim here is to present results that relate closely to the experiment. To this end, we construct and discuss the phase diagrams of an imbalanced Fermi gas in 3D, subjected to a 1D optical potential, both at fixed chemical potentials and at fixed densities. Subsequently, we investigate and discuss the effects of the properties of the optical potential on the BCS superfluid state and on the FFLO state. In Sec. II we outline our strategy to construct the phase diagrams at fixed chemical potentials and at fixed densities. Subsequently, we discuss the competition between the BCS phase and the FFLO phase in these phase diagrams, which differ qualitatively from the corresponding phase diagrams for the imbalanced Fermi gas in 3D without the optical potential. In Sec. III, we investigate the effects of the parameters of the optical potential on the FFLO state. Furthermore, we explain the remarkable effect that the maximal polarization at which the FFLO state can exist decreases when the interaction strength is increased above a certain critical value. Finally in section IV we draw conclusions.

To construct the phase diagram of an imbalanced 3D Fermi gas, subjected to a 1D optical potential, at fixed chemical potentials, we start from the saddle-point free energy per unit volume of this system, which is given by [30]:

Ω sp (µ, ζ; ∆, Q) = - 1 (2π) 2 +∞ 0 dk k +Q L -Q L dk z ×   max [|ζ k,Q | , E k ] -ξ k - ∆ 2 2 k 2 + δ 1 -cos πkz Q L   + ∆ 2 8π (1) with          ξ k = k 2 + δ 1 -cos π 2 Q Q L cos πkz Q L -µ E k = ξ 2 k + ∆ 2 ζ k,Q = ζ -δ sin π 2 Q Q L sin πkz Q L , ( 2
)
where k is defined as

k 2 = k 2 x + k 2 y . (3)
To derive expression (1), the 1D optical potential was described by using a modified dispersion relation for the fermionic particles:

ε(k, k z ) = k 2 + δ. 1 -cos πk z Q L , ( 4
)
where δ is a function of the depth of the potential V 0 and of the recoil energy E R , and Q L is the wave vector of the optical potential. Expression (4) consists of a free particle dispersion in the directions perpendicular to the optical potential, and a tight-binding periodic dispersion [32] in the direction parallel to the optical potential. Because we want to describe the FFLO state, we only consider the BCS-side of the BCS-BEC crossover [33,34]. In this region, the scattering length is always negative, hence the following units are used (1):

= 2m = -a s = 1, where a s is the s-wave scattering length. Here it must be mentioned that the derivation leading to the free energy (1) neglected the density modulation of the Fermi gas due to the 1D optical potential. This approximation leads to an underestimation of the interactions between the spin-up and spin-down fermions.

       ∂Ω sp (µ, ζ; ∆, Q) ∂∆ µ,ζ = 0 ∂Ω sp (µ, ζ; ∆, Q) ∂Q µ,ζ = 0 , (5)
which have to be satisfied simultaneously. One must be careful when extremizing the free energy, because only the minima correspond to a stable physical state. We avoided this difficulty by explicitly minimizing the free energy, instead of searching for roots of (5).

Depending on the values of the parameters µ and ζ, three different ground states can emerge:

1) a spin-balanced (S-b) superfluid, with ∆ = 0 and Q = 0 , 2) the FFLO state with ∆ = 0 and Q = 0, and 3) the normal state with ∆ = 0. For the normal state, the value of Q has no physical meaning. Analytically this can be seen by the fact that the gap equation always has a solution ∆ = 0 and that this solution is independent of the value of Q. Mathematically speaking, the value of Q is arbitrary. In the phase diagrams we will make a distinction between the partially polarized normal state (PPN) and the fully polarized normal state (FPN). In the latter, all particles are in the same spin-state, while in the former both spinstates are still present. Following this procedure, the phase diagram as a function of µ and ζ can be constructed. Before doing so, the interval of the chemical potential µ that is relevant to describe the FFLO state must be estimated. As mentioned earlier, we only consider the BCS-side of the BCS-BEC crossover. This means that the variable (k F a s ) -1 , with k F the Fermi wave vector, which characterizes the interaction strength of a two-component Fermi gas in 3D, has to lie in the interval ]-∞, 0]. At unitarity, (k F a s ) -1 = 0, which means that, in our description, the density has to go to infinity since we work in units where a s = -1.

Even at (k f a s ) -1 = 0, the chemical potential is still a fraction of the Fermi energy, and since the Fermi energy goes to infinity it follows that µ = +∞. In practice, it turns out that µ ≈ 6 lies close enough to unitarity to describe the full stability region of the FFLO phase.

The BCS-limit is already attained within good approximation at (k F a s ) -1 = -2. In this limit, µ = E F from which it follows that µ = 1/4.

Figure 1 shows a comparison between the phase diagrams in the (µ, ζ)-plane of an imbalanced Fermi gas in 3D (figure 1 (A)), and an imbalanced Fermi gas in 3D, subjected to a 1D optical potential (figure 1 (B)). In both these phase diagrams the ground state of the system is determined for a discrete set of (µ, ζ)-points and each of these points is displayed. We chose this type of visualization so that later on in this paper, the transformation to the phase diagram at fixed densities will be clarified by the transformation of the set of (µ, ζ)-points to a set of 1 k F as , P -points, where P denotes the polarization. In both phase diagrams in figure 1, the system is always in a spin-balanced superfluid state when ζ = 0. This is not visible in Fig. 1 since in this figure ζ ≥ 0.06. For small values of the total chemical potential µ (µ ≤ 0.5 in Figs. 1 (A) and (B)), the system makes a transition into a partially polarized normal state when ζ becomes of the order of the gap ∆. For larger values of µ, the system goes through all three different states when the imbalance chemical potential ζ is increased from zero. For low values of ζ, the system remains a spin-balanced superfluid. This means that a finite energy is needed in order to convert particles from one spin species to the other.

The minimal amount of energy that is required for this conversion increases when the total chemical potential increases. For high enough values of ζ the Fermi gas becomes polarized.

At this point the system makes a transition into the FFLO state, both in figure 1 (A) and (B). Characteristic of the FFLO state is that it has an oscillating order parameter [23,24]. This order parameter is enhanced by the presence of the 1D potential, which results in a much larger FFLO region in Fig. 1 (B) when compared to the case without the optical potential in Fig. 1 (A). When the imbalance chemical potential ζ increases further, the system will go over into a partially polarized normal (PPN) state and eventually into a fully polarized normal (FPN) state, and superfluidity is destroyed.

In the phase diagram in Fig. 1, we made the following choice for the depth s = V 0 /E R and for the wave length λ of the optical potential: s = 5 and λ = 1000nm, which are both realizable in current experiments. The depth of the optical potential must be chosen carefully, because the dispersion relation ( 4) that is used in our description is only valid in For each point, the ground state of the system is determined. There are four possible ground states: a spin-balanced (S-b) superfluid (red circles), the FFLO phase (green triangles), the partially polarized normal (PPN) phase (light blue squares) and the fully polarized normal (FPN) phase (dark blue diamonds). The presence of the 1D optical potential leads to a significant enlargement of the FFLO region. We used s = V 0 E R = 5 and λ = 1200nm for the depth and wave length of the optical potential, and |a s | = 500nm for the s-wave scattering length.

the tight-binding limit. This will be discussed in more detail in Sec. III A. Both parameters s and λ correspond to specific values of δ and Q L , which are defined in expression (4). Since we use units of scattering length, a fixed value of this quantity has to be chosen, in order to compute the values of δ and Q L . Here we choose |a s | = 500nm, which corresponds to a situation close to a Feshbach resonance.

B. Phase diagram at fixed densities

In the previous section, we chose to construct a phase diagram at fixed chemical potentials µ and ζ. This is not the only possible choice. Another possibility is to work at fixed total density n = n ↑ + n ↓ and fixed polarization P = (n ↑ -n ↓ ) / (n ↑ + n ↓ ), which is equivalent to fixing the density of both spin species n ↑ and n ↓ . To transform between these two descriptions, the following number equations have to be solved:

       n = - ∂Ω sp (µ, ζ; ∆, Q) ∂µ ζ δn = - ∂Ω sp (µ, ζ; ∆, Q) ∂ζ µ . (6)
In principle, we need to use the chain rule in ( 6 The latter can also be written as the polarization P = δn/n. The final step is to find the interaction strength (k F a s ) -1 that corresponds to the density n. Since the scattering length a s is kept constant, this can be done by using the following expression that relates the density to the Fermi wave vector k F :

n =          Q L 2π 2 (k 2 F -δ) (E F ≥ 2δ) Q L 2π 3   (k 2 F -δ) arccos δ -k 2 F δ + δ 1 - δ -k 2 F δ 2   (E F < 2δ) . (7)
This expression (7) is derived by calculating the number of particles at temperature zero that have an energy that lies lower than the Fermi energy, where the particles are assumed to have a dispersion relation given by (4). The procedure outlined above yields the phase the case of zero polarization, the system is in the spin-balanced superfluid state. Figure 2 shows that there exists a portion of the phase diagram at finite polarization, where no points are present. Mathematically, this corresponds to the case when there exists no combination of density and polarization that can satisfy both number equations (6). Physically, this is a state where the system separates into two different phases. To understand which of the possible phases can coexist in this phase separated state, it is instructive to look at the phase diagram at fixed chemical potentials (Fig. 1 (B)). In the case of phase separation, both phases coexist at equal chemical potentials. This means that phase separation corresponds to points in the (µ, ζ)-phase diagram that lie at the boundaries between different phases.

Wherever the FFLO phase exists, this phase is sandwiched between the Sb-SF phase and the normal phase. Hence there exists no boundary between the Sb-SF phase and the normal phase, except at those values of the interaction strength where no FFLO is present. In figure 2, it is visible that the region of phase separation only occurs between the Sb-SF phase and the FFLO phase, and not between the FFLO phase and the normal phase. This region of phase separation thus corresponds to the boundary between Sb-SF and FFLO in figure 1.

This means that all points within this phase separated region represent a phase separated state where a certain amount x of Sb-SF and an amount (1 -x) of FFLO coexist, where

x varies depending on which point is considered. We can thus conclude that the phase separated state only consists of a spin-balanced superfluid phase (S-b SF) and an FFLO phase (except at those values of the interaction strength where no FFLO is present, in which case the phase separated state consists of a S-b SF phase and a normal phase). When polarization increases further, the system makes a transition into the FFLO state. Figure 2 also shows that when the interaction strength increases, both the boundary between the phase separation region and the FFLO state and the boundary between the FFLO state and the normal state, lie at higher values of polarization. Interestingly though, for both boundaries there exists a critical interaction strength at which this polarization becomes maximal. These two maxima are indicated by a red square and a red circle respectively in figure 2. If the interaction strength becomes larger than these two critical values, both the region of phase separation and the FFLO region seem to shrink. This counterintuitive effect will be discussed in the next section. When polarization increases further, the energy cost of forming pairs with finite momentum will eventually become too high and the system will make a transition into a normal Fermi gas. The points with polarization P = 1 correspond to the fully polarized region in Fig. 1. We now have given a general overview of the phase diagrams of an imbalanced 3D Fermi gas subjected to a 1D potential. In the next section we will investigate the effects of the properties of the optical potential on the FFLO state.

III. EFFECTS OF THE 1D POTENTIAL ON THE FFLO STATE

A. Role of the wave length and the depth of the optical potential

There are two parameters that determine the properties of the 1D optical potential: the potential depth s, and the wave length λ. Before discussing the effect of these parameters, it is important to know to which intervals they are limited in our current treatment. For the wave length, we take realistic values corresponding to present day experiments: λ ∈ [600nm, 1200nm]. We have to be careful, however, when choosing values for the depth of the optical potential. There are two physical reasons to set a lower boundary to this depth.

Firstly, since the optical potential is described by a tight-binding dispersion relation (Eq. 4), the depth of the potential must be large enough to justify this approximation. To quantify this, we calculated the energy spectrum of a particle in a periodic potential, and compared it to the tight-binding dispersion given by ( 4). For completeness, we added a concise overview of this calculation in appendix A. This calculation shows that the error made by using expression (4) for the dispersion is less than five percent when s ≥ 4. Secondly, by using the dispersion relation ( 4), we only take into account the lowest Bloch band. This description is justified as long as the system has a Fermi energy which lies below the second Bloch band. A Fermi energy that lies between the top of the first Bloch band and the bottom of the second Bloch band is allowed, since the dispersion in the direction perpendicular to the direction of the optical potential is still the same as in the free 3D Fermi gas case. We have verified that by considering a potential depth s ≥ 4, the system has a Fermi energy that lies below the bottom of the second Bloch band, for all cases considered in this paper. Based on this analysis, we will use the lower boundary s = 4 for the potential depth throughout the rest of this paper. There is also a reason to set an upper limit to the depth of the optical potential. When the potential becomes too deep, the 3D Fermi gas will transition into a series of weakly-coupled two-dimensional (2D) pancakes. Such a 2D Fermi gas must be described using a different renormalized interaction strength [35]. We therefore choose s to lie between 4 and 6, so that both the 3D description and the tight-binding approximation are justified.

We now proceed to investigate the effects of both optical potential parameters λ and s.

Figure 3 shows several phase diagrams at fixed density and polarization, for an imbalanced figures show that the maximal polarization at which the FFLO state can exist increases when the wave length of the optical potential becomes smaller. This effect agrees with our expectations and can best be explained in terms of the Fermi surfaces in momentum space.

When an imbalance exists between spin-up and spin-down particles, an energy gap opens up between the Fermi surfaces of the two spin species. This is akin to a Zeeman-splitting effect.

In order to bridge this gap and to let particles pair with the same energy, the FFLO state forms pairs with a finite momentum Q. The larger the imbalance, the larger the energy gap becomes and the larger Q has to be. In Ref. [30] we demonstrated that the wave vector of the FFLO state cannot become larger than the wave vector of the optical potential. Since a smaller wave length corresponds to a larger wave vector, the maximal value of the wave vector Q that can be attained by the FFLO state increases when the wave length of the optical potential decreases. This means that a larger energy gap can be bridged and thus a larger polarization can be accommodated. Figures 3 (C) and 3 (D) show the phase diagram of the system at the same wave length but at different potential depth. In both figures, the maximal polarization of the FFLO state is equal. We conclude that only the wave length of the optical potential has an influence on the capability of the system to remain in the FFLO state while being polarized, at least within the small range of potential depth s that is considered here.

B. Effect of the edge of the Brillouin zone

One effect that is visible throughout all phase diagrams in figures 2 and 3, is that when the interaction strength exceeds a certain critical value, both the region of phase separation and the FFLO region become less stable with respect to polarization (i.e. the maximal polarization at which they can exist, decreases) . This can be seen by the fact that both the boundary between the phase separation region and the FFLO state and the boundary between the FFLO state and the normal state have a maximum (indicated by a red square and a red circle respectively in figure 2). Moreover, these two maxima seem to correspond to different values of the interaction strength. This effect is counterintuitive, because when the interaction strength increases, the binding energy of the fermion pairs increases likewise, which in general should strengthen superfluidity, not weaken it.

This effect is explained with the help of figure 2. In this figure, the empty blue squares with one of the two maxima can be explained as follows. Consider the case which is shown schematically in Fig. 4. Figure 4 (A) shows the Fermi surface of both spin-up and spin-down particles at low polarization. In order for the system to become superfluid, pairing must occur between particles that lie in a band with thickness ∆ around the Fermi energy (grey bands in figure 4). When these pairing bands of the spin-up and spin-down particles overlap, pairing is possible and the system is in the superfluid state. When polarization is introduced, a gap arises between the two Fermi surfaces, as indicated in figure 4 (A). As long as this gap is smaller than ∆, there is overlap of the pairing bands, and superfluidity is possible in principle. When the polarization becomes too large, the overlap disappears and BCS superfluidity is no longer possible (figure 4 (B)). In order to remain superfluid, the system can transition into the FFLO state, which translates the Fermi surface of the minority spins so that it locally lies on the Fermi surface of the majority spins, making pairing possible again (Fig. 4 (B)). This allows the system to remain superfluid, at the cost of paying some energy to form particles with finite center-of-mass momentum. This is the case for a 3D

Fermi gas, but when a 1D potential is added, an additional effect comes into play. When the density and/or polarization is high enough so that the Fermi surface of the majority spins lies at the edge of the first Brillouin zone (this is true for the area at the right side of the dashed oblique line in Fig. 2), pairing at this Fermi surface becomes more difficult, because for k z > Q L , there are no available states (see Fig. 4 (C)). Hence, the formation of an FFLOtype superfluid phase will be hindered. This now explains why, at equal polarization, the normal phase becomes favored over the FFLO phase when the interaction strength exceeds the value corresponding to the situation where the Fermi energy of the majority spins equals the top of the first Bloch band (indicated by a red circle in Fig. 2). The question can be raised whether the FFLO state can still exist if the fermion pairs acquire a momentum in the x or y direction, since there exists no Brillouin zone in these directions. We did not treat this case in our paper because we assumed the FFLO momentum to be of the form

Q = (0, 0, Q z ).
The question remains why the pairing mechanism for the spin-balanced superfluid phase becomes hindered at a certain value of the interaction strength. The reason for this is the following. As mentioned before, the formation of a spin-balanced superfluid (S-b SF) can only occur when pairing takes place between up-and down-particles lying within a certain energy band ∆ around the respective Fermi energy of both particles. Consider now the case of zero polarization, as shown in figure 4 (D). When the density of the system becomes high enough, the bands around the Fermi surfaces of spin-up and spin-down particles (which fully overlap at polarization zero) will start to overlap with the edge of the Brioullin zone. This hinders the formation of pairs, because there are no states at energy values higher than the Fermi energy to which particles can scatter.

IV. CONCLUSIONS

We have investigated the effects of a 1D optical potential on the FFLO state in an imbalanced Fermi gas in 3D. To allow for a direct connection with the experiments, we constructed the phase diagram of the system, both at fixed chemical potentials and at fixed densities of both spin species. Subsequently, we investigated the effects of the depth and the wave length of the optical potential on the FFLO state. It was shown that this state can exist under a larger polarization, when the wave length of the optical potential becomes smaller.

The depth of the optical potential did not have any appreciable effect, within the range that can be considered in our tight-binding approximation. By studying the phase diagram at fixed density we have discovered an unexpected effect. When the interaction strength exceeds a certain critical value, the maximal polarization at which both the FFLO phase and the phase separation region can exist decreases. The underlying reason was discussed in detail. We have found that when the density of the system becomes so high that the Fermi energy exceeds the top of the lowest Bloch band, superfluid pairing becomes hindered. This is because in that case, particles with the highest energy lie at the edge of the first Brillouin zone. Since no states exist outside of this zone, pairing is not possible and the superfluid state is suppressed. We expect this effect, along with the enhancement of the FFLO state, to be observable in experiments where a 1D optical potential is applied to an imbalanced Fermi gas in 3D.

Re-indexing the second sum over n and identifying the coefficients on both sides in equation (A5) leads to

2 2m (k z + 2nQ L ) 2 + sE R 2 -ε kz c n - sE R 4 (c n+1 + c n-1 ) = 0. (A6)
This expression represents an infinite set of linear equations. In practice, it is enough to consider five equations, since the contributions of the eigenvectors c n decreases exponentially for large n. Solving equation (A6) numerically yields the exact energy spectrum of a particle in a periodic potential, given by (A1). Figure 5 shows a comparison between this exact spectrum and the corresponding dispersion relation in the tight-binding limit, given by (4).

For s ≤ 2 the difference is substantial, but for s ≥ 4, the approximation deviates less than 5 percent from the exact result. This justifies the approximation that we make by considering s ≥ 4. between the two Fermi surfaces is enlarged and eventually, there will be no overlap between the two pairings bands. At this point, BCS superfluidity is no longer possible. However, the system can remain a superfluid by translating the minority Fermi surface (dashed red ellipsoid) so that locally, both Fermi surfaces are re-aligned (indicated by a black rectangle). (C) When the density is increased, the Fermi surface of the majority spins (here spin-up) will lie at the edge of the first Brillouin zone. This will hinder FFLO-like pairing because no states are available at k z > Q L (indicated by a black filled rectangle). (D) Above a critical density, the Fermi surface lies at the edge of the Brillouin zone even for polarization zero. At this point, BCS superfluidity will also be hindered.

FIG. 1 :

 1 FIG. 1: Comparison between the phase diagram in the (µ, ζ)-plane (µ and ζ are in units E 0 = 2 2ma 2 s ) of (A) an imbalanced Fermi gas in 3D, and (B) an imbalanced Fermi gas in 3D, subjected to a 1D optical potential. Both phase diagrams are constructed based on a finite set of points in the (µ, ζ)-plane. For each point, the ground state of the system is determined. There are four possible

  ), since ∆ and Q are separate variables that depend on µ and ζ. But the saddle-point condition (5) will make these terms disappear. If fluctuations around the saddle point are taken into account, these extra terms will remain. The inclusion of fluctuations around the saddle point lies beyond the scope of this paper. Starting from the phase diagram at fixed chemical potentials in figure 1 (B), the phase diagram at fixed densities can now be constructed. Every point in the phase diagram depicted in figure 1 (B) corresponds to a set of values for the parameters (µ, ζ; ∆, Q). Inserting these values into both equations in (6) yields the density n and the density difference δn = n ↑ -n ↓ .

  FIG. 2: Phase diagram of an imbalanced Fermi gas in 3D, subjected to a 1D optical potential, at fixed density and fixed polarization P. The dimensionless parameter (k F a s ) -1 characterizes the interaction strength. Each point in this phase diagram is a transformation of a corresponding point in the phase diagram in figure 1 (B). At polarization zero, the system is in the spin-balanced (S-b) superfluid phase. When polarization becomes finite, there is a region where no points exist. This means that no combination of density and polarization can solve the number equations, and the system separates into two phases: a spin-balanced superfluid component and an FFLO component. Interestingly, both the boundary between the phase separation region and FFLO, and between FFLO and the normal phase have a maximum (indicated by a red square and a red circle respectively). The dashed oblique line indicates the boundary where the Fermi energy of the majority spins becomes equal to the top of the first Bloch band, which is equal to 2δ. The empty blue squares (normal phase) and empty green triangles (FFLO) indicate that the Fermi energy of the majority spins lies higher than the top of the first Bloch band. The vertical dashed line indicates the value of the interaction strength where the Fermi energy at polarization zero becomes equal to the top of the first Bloch band. Both dashed lines intersect with one of the two indicated maxima. (s = 5, λ = 1200nm, and |a s | = 500nm)

FermiFIG. 3 :

 3 FIG. 3: Overview of several phase diagrams of an imbalanced Fermi gas in 3D subjected to a 1D optical potential, for different values of the depth s and the wave length λ of the potential. The maximal polarization P at which the FFLO state can exist (indicated by a dashed horizontal line) increases when the wavelength of the optical potential becomes smaller. At constant wave length, the depth of the optical potential does not have an influence on this maximal polarization, at least within the range 4 ≤ s ≤ 6 that is considered here. (S-b SF = spin-balanced superfluid, PS = phase separation)

(

  normal state) and the empty green triangles (FFLO state) indicate that the Fermi energy of the majority spins lies higher than the top of the first Bloch band, which is equal to 2δ. The boundary at which the Fermi energy of the majority spins becomes equal to the top of the first Bloch band is indicated by the dashed oblique line. The maximum polarization under which the FFLO state can exist (red circle), lies exactly at this line. Furthermore, the vertical dashed line indicates the value of the interaction strength at which the Fermi energy at polarization zero becomes equal to the top of the first Bloch band. This dashed vertical line intersects with the maximum polarization under which the phase separation region can exist (red square). The reason why each of these two dashed lines intersects

FIG. 4 :

 4 FIG. 4: Schematic representation of the pairing mechanisms in an imbalanced 3D Fermi gas subjected to a 1D optical potential. For the sake of clarity, we have omitted the y direction. (A) The Fermi surfaces of the spin-up and the spin-down particles at low polarization, which creates a gap between both Fermi surfaces. Pairing occurs between particles that lie in a band with thickness ∆ around the Fermi surfaces of up-and down spins (grey bands). As long as these two pairing bands overlap, BCS superfluidity is possible in principle. (B) When polarization is increased, the gap

  The free energy (1) depends on two thermodynamic variables µ and ζ, and on two variational parameters ∆ and Q. The variables µ and ζ are the total chemical potential µ = (µ ↑ + µ ↓ ) /2, which fixes the total density, and the imbalance chemical potential ζ = (µ ↑ -µ ↓ ) /2, which fixes the polarization. A description in terms of µ and ζ is equiva-

	lent to a description in terms of the respective chemical potentials of both spin species µ ↑
	and µ ↓ . The variational parameters ∆ and Q denote the binding energy and the center-
	of-mass momentum of the fermion pairs respectively. The parameter Q was introduced to
	include the FFLO state (which is characterized by fermion pairs with finite momentum) in
	our description. To find the ground state of the system, given fixed values of µ and ζ, the
	saddle-point free energy Ω sp (µ, ζ; ∆, Q) has to be minimized with respect to the variational
	parameters ∆ and Q. This defines the saddle-point equations

Acknowledgments

Acknowledgments -The authors would like to thank Gerard 't Hooft, Carlos Sá de Melo, and Sergei Klimin for fruitful discussions. JPAD wishes to thank Nick Van den Broeck for helpful discussions. One of the authors (JPAD) acknowledges a Ph. D. fellowship of the Research Foundation -Flanders (FWO). This work was supported by FWO-V projects G.0356.06, G.0370.09N, G.0180.09N, G.0365.08.

APPENDIX A: RANGE OF OPTICAL POTENTIAL DEPTHS THAT ALLOW THE TIGHT-BINDING APPROXIMATION

The energy spectrum of a particle in a periodic potential can be calculated using Schrödinger's equation. The potential under consideration is given by

where s is the depth of the optical potential in units of the recoil energy E R and Q L is the wave vector of the optical potential given by Q L = π/λ. Given this potential, Schrödinger's equation states that

Since we consider a particle in a periodic potential, Bloch's theorem applies and we can re-write the wave function ψ kz (z) as

where u kz (z) has the same periodicity as the potential. Because u kz (z) is periodic, it can be written in a Fourier series

FIG. 5: Comparison between the exact energy spectrum of a particle in a periodic potential (solid blue lines) and the dispersion relation calculated in the tight-binding limit (long-dashed green lines),

given by expression (4). The horizontal short-dashed red lines indicate the two lowest energy levels of the harmonic oscillator, derived by approximating (4) to second order in the momentum. For s ≥ 4, the difference between the exact energy spectrum and the tight-binding dispersion is less than 5%.


Miguel A Fernández
email: miguel.fernandez@inria.fr

Jimmy Mullaert
email: jimmy.mullaert@inria.fr

Displacement-velocity correction schemes for incompressible fluid-structure interaction

We propose a new class of time-marching schemes for the explicit coupling of an incompressible fluid and an elastic solid (not necessarily thin). We state a general energy-based stability result and illustrate the accuracy of the different variants in a numerical benchmark.

Dans cette note, nous proposons une généralisation (la première, à notre connaissance) des schémas de couplage explicite introduits dans [START_REF] Guidoboni | Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow[END_REF][START_REF] Fernández | Incremental displacement-correction schemes for the explicit coupling of a thin structure with an incompressible fluid[END_REF] dans le cas d'un solide épais. On considère ainsi le problème couplé (1)- [START_REF] Burman | Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility[END_REF], qui modélise l'interaction entre un fluide linéaire newtonien incompressible et une structure élastique linéaire (avec ou sans amortissement), dont la formulation variationnelle est donnée par [START_REF] Guidoboni | Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow[END_REF]. Le détail des schémas de couplage explicite proposés est présenté dans l'Algorithme 1. L'idée principale regroupe les trois ingrédients suivants. Le premier consiste à traiter explicitement l'amortissement du solide (si pris en compte dans (1)) par rapport au fluide. Le deuxième repose sur le calcul découplé de la vitesse intermédiaire du solide (interface et contributions internes, première et deuxième étapes de l'Algorithme 1) par condensation de la de masse dans la structure. Enfin, le dernier consiste à extrapoler à la fois le déplacement et la vitesse du solide (schémas incrémentaux).

Nous présentons également un résultat général de stabilité qui traite à la fois les variantes increméntale et non-incrémentales (voir la Proposition 4.1). En particulier, ce résultat montre que le schéma nonincrémental et le schéma incrémantal avec extrapolation de premier ordre sont inconditionnellement stables. La variante incrémentale avec extrapolation d'ordre deux est stable sous une condition de type CFL (hyperbolique/parabolique). Finalment, la précision supérieure des schémas incrémentaux est illustrée par des expériences numériques (voir la Figure 1).

Introduction

The stability of the numerical approximations of fluid-structure interaction problems involving a viscous incompressible fluid and an elastic structure is very sensitive to the way the interface coupling conditions (kinematic and kinetic continuity) are enforced at the discrete level (see, e.g., [START_REF] Fernández | Algorithms for fluid-structure interaction problems[END_REF] for a review). Stable explicit coupling schemes (i.e., that only involve one, or a few, fluid and solid resolutions per timestep) have only recently been proposed in the literature [START_REF] Burman | Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility[END_REF][START_REF] Guidoboni | Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow[END_REF][START_REF] Fernández | Incremental displacement-correction schemes for the explicit coupling of a thin structure with an incompressible fluid[END_REF]. In [START_REF] Burman | Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility[END_REF], stability is achieved through an appropriate Robin-Robin treatment of the interface coupling and the addition of a weakly consistent interface compressibility term. Yet, defect-correction iterations are required to guarantee accuracy. In the case of thin structures (e.g., string, plate, membrane or shell models), the procedures reported in [START_REF] Guidoboni | Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow[END_REF][START_REF] Fernández | Incremental displacement-correction schemes for the explicit coupling of a thin structure with an incompressible fluid[END_REF] combine the splitting of the solid time-marching with an implicit coupling of the interface fluid stress and the hydrodynamic solid contributions (inertia and damping), fully embedded into the fluid sub-problem through a generalized Robin boundary condition. These coupling procedures are a class of displacement-correction schemes. The non-incremental variant [START_REF] Guidoboni | Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow[END_REF] ignores the displacement in the first step and, hence, may lack accuracy. The incremental schemes [START_REF] Fernández | Incremental displacement-correction schemes for the explicit coupling of a thin structure with an incompressible fluid[END_REF], on the contrary, achieve optimal accuracy by extrapolating the displacement in the first step (corrected accordingly in the second).

In this note, we propose a generalization (the first, to the best of our knowledge) of the explicit coupling schemes reported in [START_REF] Guidoboni | Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow[END_REF][START_REF] Fernández | Incremental displacement-correction schemes for the explicit coupling of a thin structure with an incompressible fluid[END_REF] to the case of thick structures: two-and three-dimensional (possibly damped) linear elasticity. These new schemes are based on the following three main ideas. The first consists of treating explicitly (whenever present) the solid damping with respect to the fluid. The second relies on the uncoupled computation of the intermediate solid velocity (interface and internal contributions) via a mass-lumping approximation of the structure inertia. At last, the third consists of extrapolating both the displacement and the velocity (incremental schemes). We present a general stability result which covers both the non-incremental and the incremental variants. The analysis shows, in particular, that the non-incremental and the first-order extrapolated incremental schemes are unconditionally stable. The incremental variant with second-order extrapolation is stable under a CFL-like (hyperbolic/parabolic) condition. The superior accuracy of the incremental variants is illustrated in a numerical benchmark.

A linear model problem

We consider a low Reynolds regime and assume that the interface undergoes infinitesimal displacements. The fluid is described by the Stokes equations, in a fixed domain Ω f ⊂ R d (d = 2, 3), and the structure by the linear (possibly damped) elasticity equations, in the solid domain Ω s ⊂ R d . We denote by Σ def = ∂Ω s ∩ ∂Ω f the fluid-structure interface and ∂Ω f = Γ ∪ Σ and ∂Ω s = Γ d ∪ Γ n ∪ Σ are given partitions of the external boundary. Our linear coupled problem reads as follows: Find the fluid velocity u :

Ω f ×R + → R d , the fluid pressure p : Ω f × R + → R, and the structure displacement d : Ω s × R + → R d such that      ρ f ∂ t u -divσ f (u, p) = 0 in Ω f , div u = 0 in Ω f , σ f (u, p) n f = -p Γ n f on Γ,            ρ s ∂ t ḋ + c 1 ḋ -divσ s (d, ḋ) + c 0 d = 0 in Ω s , ḋ = ∂ t d in Ω s , d = 0, β ḋ = 0 on Γ d , σ s (d, ḋ) n s = 0 on Γ n , (1)
satisfying the interface coupling conditions

u = ḋ on Σ, σ s (d, ḋ)n s = -σ f (u, p)n f on Σ (2)
and the initial conditions u(0) = u 0 , d(0) = d 0 and ḋ(0) = v 0 . Here, ρ f , ρ s > 0 stand for the fluid and solid densities, σ f (u, p) def = -pI + 2µ (u) for the fluid stress tensor, µ > 0 for the fluid dynamic viscosity, (u) def = (∇u + ∇u T )/2, p Γ for a given pressure on Γ and n f , n s for the exterior unit normal vectors.

The solid stress tensor is given by σ s (d, ḋ)

def = σ(d) + βσ( ḋ), with σ(d) def = 2L 1 (d) + (L 2 divd)I and L 1 , L 2 > 0 the Lamé constants.
The presence of dissipative effects in the structure is described by the term c 1 ḋ -βdivσ( ḋ), the so-called Rayleigh damping (see, e.g., [START_REF] Hughes | The finite element method[END_REF]), with c 1 def = α(1 + βc 0 ) and c 0 , α, β ≥ 0. The zeroth-order term, c 0 d, allows to incorporate the transversal membrane effects that appear in axisymmetric formulations. In what follows, we shall make use of the functional spaces

V f def = [H 1 (Ω f )] d , Q def = L 2 (Ω f ), V s def = {v s ∈ [H 1 (Ω s )] d , v s | Γ d = 0} and V s def = {v s ∈ [H 1 (Ω s )] d , v s | Γ d ∪Σ = 0}. We also set a(u, v f ) def = 2µ (u), (v f ) Ω f , b(p, v f ) def = -(p, divv f ) Ω f , p Γ , v f def = -(p Γ , v f • n f ) Γ , a e (d, v s ) def = σ(d), (v s ) Ω s + c 0 (d, v s ) Ω s , a v ( ḋ, v s ) def = βa e ( ḋ, v s ) + α( ḋ, v s ) Ω s .
The symbol (•, •) ω denotes de inner-product of L 2 (ω) and • 0,ω its associated norm. Problem (1)-( 2) can then be rewritten in variational form as follows: For t > 0, find the fluid velocity and pressure (u(t), p(t)) ∈ V f × L 2 (Ω) and the solid displacement and velocity d(t), ḋ(t) ∈ V s , such that ḋ = ∂ t d and

u| Σ = ḋ| Σ , ρ f ∂ t u, v f Ω f + a(u, v f ) + b(p, v f ) -b(q, u) + ρ s ∂ t ḋ, v s Ω s + a v ( ḋ, v s ) + a e (d, v s ) = p Γ , v f (3) for all (v f , q, v s ) ∈ V f × Q × V s with v f | Σ = v s | Σ .

Displacement-velocity correction explicit coupling schemes

In this section, we address the time and space discretization of the coupled problem (3). The space discretization is based on continuous piecewise-affine approximations, with a lumped-mass approximation in the structure. We denote by

V f h ⊂ V f , Q h ⊂ Q, V s h ⊂ V s and V s h def = V s h ∩ V s the corresponding
finite element spaces, and the bi-linear form s h (•, •) stands for a given symmetric pressure stabilization operator (see, e.g., [START_REF] Burman | Galerkin finite element methods with symmetric pressure stabilization for the transient Stokes equations: stability and convergence analysis[END_REF]). We assume that the fluid and solid discretizations match at the interface, viz.,

X Σ,h def = v f h | Σ , v f h ∈ V f h = v s h | Σ , v s h ∈ V s h .
In what follows, we shall make use of the standard (solid-sided) lifting operator L h : X Σ,h → V s h defined in such a way that L h ξ h | Σ = ξ h and that the nodal values of L h ξ h vanish out of Σ. Note that, for all v s h ∈ V s h , we have the decomposition

v s h = L h v s h + v s h , with v s h def = v s h -L h v s h ∈ V s h .
We denote by (•, •) s,h the lumped-mass approximation of the inner-product (•, •) Ω s (see, e.g., [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF]), and we set a v,h (•, •) def = βa e (•, •) + α(•, •) s,h . At last, we introduce the following discrete Ritz-representation operators with respect to the inner-product (•, •) s,h : For each v s ∈ V s , we define

L e h v s , L v h v s ∈ V s h such that (L e h v s , v s h ) s,h = a e v s , v s h and (L v h v s , v s h ) s,h = a v,h v s , v s h for all v s h ∈ V s h .
The discretization in time is performed with a displacement-velocity correction coupling scheme. Hence, the time-marching of ( 3) is split into three sequential sub-steps: the solid displacement and velocity are treated explicitly (or ignored) in the first two steps and then corrected in the third. The proposed fully discrete schemes are detailed in Algorithm 1, where τ > 0 denotes the time-step size,

t n def = nτ , ∂ τ x n def = (x n -x n-1 )/τ and (d h , ḋ h ) = (0, 0), (d n-1 h , ḋn-1 h ), (2d n-1 h -d n-2 h , 2 ḋn-1 h - ḋn-2 h
) are the, zeroth-, first-and second-order extrapolations, respectively. Algorithm 1 Displacement-velocity correction explicit coupling schemes.

(i) Robin-like fluid step:

Find (u n h , p n h ) ∈ V f h × Q h such that ρ f ∂ τ u n h , v f h Ω f + a(u n h , v f h ) + b(p n h , v f h ) -b(q h , u n h ) + s h (p n h , q h ) + ρ s τ L h u n h , L h v f h s,h = p Γ (t n ), v f h + ρ s τ ḋn-1 h , L h v f h s,h -a e (d h , L h v f h ) -a v,h ( ḋ h , L h v f h ) (4) for all (v f h , q h ) ∈ V f h × Q h . (ii) Intermediate solid velocity: Find w n h ∈ V s h such that ρ s τ w n h , v s h s,h = ρ s τ ḋn-1 h , v s h s,h -a e (d h , v s h ) -a v,h ( ḋ h , v s h ) (5)
for all v s h ∈ V s h . Then set w n h def = L h u n h + w n h ∈ V s h . (iii) Solid displacement-velocity correction: Find d n h ∈ V s h such that ḋn h = ∂ τ d n h and ρ s τ ḋn h , v s h s,h + a e (d n h , v s h ) + a v,h ( ḋn h , v s h ) = ρ s τ w n h , v s h s,h + a e (d h , v s h ) + a v,h ( ḋ h , v s h ) (6) for all v s h ∈ V s h .
Algorithm 1 generalizes the explicit coupling schemes reported in [START_REF] Guidoboni | Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow[END_REF][START_REF] Fernández | Incremental displacement-correction schemes for the explicit coupling of a thin structure with an incompressible fluid[END_REF] to the case of thick structures. The structural damping a v,h (•, •) is explicitly coupled to the fluid. This allows to uncouple the computation of the interface and internal components (L h u n h and w n h ) of the intermediate solid velocity w n h , via a masslumping approximation in the structure. The fluid step (i) has a Robin-like structure, due to the presence of the term ρ s /τ L h u n h , L h v f h s,h (interface lumped-mass contribution of the structure). Step (ii), which only involves the resolution of a diagonal system (internal lumped-mass contribution), is not present in the case of thin structures (see [START_REF] Guidoboni | Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow[END_REF][START_REF] Fernández | Incremental displacement-correction schemes for the explicit coupling of a thin structure with an incompressible fluid[END_REF]). Observe that the addition of (4), ( 5) and ( 6) yields

ρ f ∂ τ u n h , v f h Ω f + a(u n h , v f h ) + b(p n h , v f h ) -b(q h , u n h ) + s h (p n h , q h ) + ρ s ∂ τ ḋn h , v s h s,h + a e (d n h , v s h ) + a v,h ( ḋn h , v s h ) = p Γ (t n ), v f h (7) for all (v f h , q h , v s h ) ∈ V f h × Q h × V s h such that v f h | Σ = v s h | Σ .
On the other hand, from ( 6) and (ii) we get

u n h = ḋn h + τ ρ s L e h d n h -d h + τ ρ s L v h ( ḋn h -ḋ h ) on Σ. (8)
Therefore (as in the case of thin-solid models [START_REF] Fernández | Incremental displacement-correction schemes for the explicit coupling of a thin structure with an incompressible fluid[END_REF]), the explicit coupling schemes reported in Algorithm 1 can be interpreted as implicit time discretizations of (3) 2 (given by ( 7)) that involve a weakly consistent correction of the kinematic constraint (3) 1 (given by ( 8)).

Stability

We define the discrete energy and dissipation of the fluid-structure system, at time t n , as

E n h def = ρ f 2 u n h 2 0,Ω f + ρ s 2 ḋn h 2 0,Ω s + 1 2 a e (d n h , d n h ), D n h def = τ n m=1 2µ (u m h ) 2 0,Ω f + a v ( ḋm h , ḋm h ) .
We have the following stability result, whose proof can be found in [START_REF] Fernández | Diplacement-velocity correction schemes for incompressible fluid-structure interaction: Stability analysis and numerics[END_REF]. The symbol indicates an inequality up to a multiplicative constant (independent of the physical and discretization parameters). Proposition 4.1 Assume that p Γ = 0 (free system).

-Non-incremental scheme (zeroth-order extrapolation): For n ≥ 1, there holds

E n h + D n h E 0 h .
-Incremental scheme, first-order extrapolation: For n ≥ 1, there holds

E n h + D n h E 0 h + τ 2 a e ( ḋ0 h , ḋ0 h ) + τ 2 ρ s L e h d 0 h + L v h ḋ0 h 2 0,Ω s .
-Incremental scheme, second-order extrapolation: Under the conditions βρ s ω 2 e τ + ατ h 2 < ρ s h 2 , τ ω

v ) -1 -τ ) E 1 h + a v ḋ1 h - ḋ0 h , ḋ1 h - ḋ0 h , (10)
where ω e represents a maximum elastic-wave speed in the solid. Some remarks are in order. The above result shows that the non-incremental and the first-order extrapolated incremental schemes are unconditionally stable in the energy norm (the latter under additional regularity on the initial data). Proposition 4.1 also shows that, with second-order extrapolation, the incremental variant is conditionally energy stable. For undamped solids (i.e., α = β = 0), the stability condition (9) becomes a 6/5-CFL constraint and the estimates of Proposition 4.1 are similar to those obtained in [START_REF] Fernández | Incremental displacement-correction schemes for the explicit coupling of a thin structure with an incompressible fluid[END_REF] with a thin-solid model. For β > 0, the condition (9) involves a parabolic-CFL constraint. Note that the estimate (10) involves a constant growing exponentially in time (from Gronwall's lemma).

Numerical experiments

In order to illustrate the accuracy of the proposed schemes, we have performed a series of numerical tests in a two-dimensional pressure-wave propagation benchmark (see [START_REF] Badia | Fluid-structure partitioned procedures based on Robin transmission conditions[END_REF], for instance). The fluid and solid domains are given, respectively, by the rectangles Ω f = [0, 6] × [0, 0.5] and Ω s = [0, 6] × [0.5, 0.6]. The interface is Σ = [0, 6]×{0.5} and we have Γ d = {0, 6}×[0.5, 0.6] and Γ n = [0, 6]×{0.6}. All units are in the CGS system. At x = 0, a sinusoidal pressure of maximal amplitude 2 × 10 4 is imposed during 5 × 10 -3 s, corresponding to half a period. Zero pressure is enforced at x = 6 and a symmetry condition is applied on the lower wall y = 0. The fluid physical parameters are given by ρ f = 1.0 g.cm -3 , µ = 0.035 P. For the solid we have ρ s = 1.1g.cm -3 , L 1 = 1.15×10 6 dyn.cm -2 , L 2 = 1.7×10 6 dyn.cm -2 , c 0 = 4×10 6 dyn.cm -4 , α = β = 0 (i.e., undamped solid, for a numerical study with structural damping we refer to [START_REF] Fernández | Diplacement-velocity correction schemes for incompressible fluid-structure interaction: Stability analysis and numerics[END_REF]). The computations have been performed with FreeFem++ [START_REF] Pironneau | Freefem++[END_REF]. Figure 1 of the solid displacement, in the relative energy-norm at time t = 1.5 × 10 -2 , for the incremental, nonincremental and implicit schemes. The reference solution has been generated with the implicit scheme and a high grid resolution: τ = 2.5 × 10 -7 and h = 3.125 × 10 -3 . We have refined both in time and in space. Figure 1(left) shows that, in the case τ = O(h) (τ = O(h 6/5 ) with second-order extrap.), the incremental and implicit schemes yield an overall O(τ ) optimal accuracy. Note that the second-order extrapolation does not improve the convergence rate, since the overall time-accuracy of the underlying implicit scheme in Algorithm 1 is O(τ ). On the contrary, the non-incremental variant is unable to show a convergent behavior towards the reference solution (saturation of the error). Yet, Figure 1(right) shows that this scheme achieves an overall O(h) accuracy under a parabolic τ = O(h 2 ) time-step restriction. This suggests that the non-incremental scheme introduces a sub-optimal O(τ

2 ) perturbation in the error. The superior accuracy of the incremental variants is also highlighted by Figure 1(right).

6 5 e ≤ c e h 6 5 , βρ s ω 2 e τ 2 + αω 2 e τ h 4 ≤ c 2 v ρ s h 4

 6564 , τ (c 5 e + c 2 v ) < 1, (9) there holds, for n ≥ 2,

Figure 1 .

 1 Figure 1. Solid displacement convergence history. Left: τ = O(h) (τ = O(h 6/5 ) with second-order extrap.). Right: τ = O(h 2 ).

  reports the time-convergence history

		1				implicit
						explicit non-incremental
						explicit incremental (1st order extrap.)
						explicit incremental (2nd order extrap.)
	energy-norm relative error	0.1	implicit explicit non-incremtal	energy-norm relative error	0.1 1	slope 1 slope 2
			explicit incremental (1st order extrap.)
			explicit incremental (2nd order extrap.)
			splope 1
		0.00001	0.0001			0.01	0.1
			time-step size			mesh size


F Lecocq

J Claudon

O Buisson

P Milman

Non-linear coupling between the two oscillation modes of a dc-SQUID

Keywords: numbers: 85, 25, Cp, 03, 67, Lx

We make a detailed theoretical description of the two-dimensional nature of a dc-SQUID, analyzing the coupling between its two orthogonal phase oscillation modes. While it has been shown that the mode defined as "longitudinal" can be initialized, manipulated and measured, so as to encode a quantum bit of information, the mode defined as "transverse" is usually repelled at high frequency and does not interfere in the dynamics. We show that, using typical parameters of existing devices, the transverse mode energy can be made of the order of the longitudinal one. In this regime, we can observe a strong coupling between these modes, described by an Hamiltonian providing a wide range of interesting effects, such as conditional quantum operations and entanglement. This coupling also creates an atomic-like structure for the combined two mode states, with a V-like scheme.

The control of physical systems at the quantum level represents a significant technological challenge. But, this goal is of extreme interest since successful applications would allow us to probe fundamental aspects of quantum physics and to demonstrate the claimed advantages of quantum computation. A number of physical systems have shown, in the last decade, their potentiality as "quantum devices". Examples are atoms in a high finesse cavity [1], trapped ions [2] and solid state quantum circuits [3][4][5][6].

Superconducting circuits (SC) appear as suitable candidates to study controllable macroscopic quantum systems. During the last decade, many experiments on such devices have demonstrated effects previously observed on atomic and quantum optics experiments [3,4]. Furthermore, the extreme flexibility in choosing the characteristic parameters of SC have made possible the realization of new interaction regimes, such as the ultra strong coupling regime [7,8]. In these experiments, SC act as artificial atoms with one degree of freedom, described either as a 2-level system or a multilevel system.

In this work, we show that a dc-SQUID with a large loop inductance is an artificial atom with two degrees of freedom. We make a theoretical study of the two orthogonal phase oscillation modes and analyse the coupling between them. One of the modes, called hereafter the longitudinal mode (LM), was extensively studied and demonstrated a large variety of quantum phenomena [9]. For instance, using the two lowest energy states, all required operations necessary to define a qubit (coherent manipulation, initialization, coupling and measurement) have been performed [10]. The quantum dynamics of the transverse mode (TM), which is perpendicular to LM, have never been studied theoretically or experimentally, although it presents appealing potential advantages, such as longer coherence times and an optimal point for energy relaxation.

Moreover, we also show that the nature and the strength of the coupling terms between the two oscillation modes can be engineered to a large extent. They lead to a wealth of novel effects, that were initially demonstrated or proposed in atomic physics. Examples include quantum non-demolition measurements, the realization of conditional quantum operations, the controlled creation and measurement of entanglement between the two mode. In addition, the coupling term in the dc-SQUID can result in a V-type 3-level structure, that has recently attracted a lot of interest [11]. All these potential experimental achievements involve realistic circuit parameters and should be observable with present day technology.

A dc-SQUID is composed of two Josephson junctions (JJ) with critical current I 0 , capacitance C 0 and phase difference φ 1 and φ 2 , embedded in a loop of inductance L (Fig. 1a). The phase dynamics of this system can be mapped on the evolution of a fictitious particle of mass m = 2C 0 (Φ 0 /2π), with two degrees of freedom x = (φ 1 + φ 2 )/2 and y = (φ 1φ 2 )/2, which is experiencing a potential [12,13]:

U (x, y) = U 0 -cos x cos y - I b 2I 0 x + b y - πΦ b Φ 0 2 , (1)
where Φ 0 = h/2e is the superconducting flux quantum, U 0 = 2I 0 (Φ 0 /2π), I b and Φ b are respectively the current and flux bias and b = Φ 0 /2πLI 0 the ratio between the Josephson inductance and the loop inductance. For I b < 2I 0 , this potential presents a periodic series of equivalent local minima separated by potential barriers (Fig. 1b). For simplicity, we consider a perfectly symmetric dc-SQUID, with symmetric inductances, JJs critical currents and capacitances. To describe the 2D dynamics of the particle trapped in a potential minimum, we introduce the basis (x , x ⊥ ) corresponding to the motion parallel ( ) or transverse (⊥) to the trajectory of minimum energy [14]. It is obtained from a rotation with an angle θ(I b , Φ b ) of the basis (x, y). The system dynamics is then governed by the Hamiltonian H = H + C + H ⊥ which describes two anharmonic oscillators whose interaction is mediated by a coupling term C (Fig. 1c). Up to order 4,

H α = ω α × [ 1 2 (p 2 α + x2 α ) -σ α x3 α + δ α x4 α ] with α = or ⊥.
Here, xα and pα are the reduced (dimensionless) position and momentum operators in the direction α; they satisfy the commutation relation [x α , pα ] = i 1; the ω α are the characteristic frequencies of each mode; in general, they are different with ω < ω ⊥ . For relevant bias points, the oscillators are weakly anharmonic: δ α ≪ 1 and σ α ≪ 1. Thus the energy E nα associated with the eigenstate |n α of H α can be obtained from a simple perturbative approach. Specifically,

E nα -E nα-1 = (1 -Λ α n α ) ω α with Λ α = 15 2 σ α + 3δ α [10]. The coupling Hamiltonian reads C = g 21 x2 x⊥ + g 12 x x2 ⊥ + g 22 x2 x2
⊥ . Note that here the linear coupling term x x⊥ is cancelled by the (x , x ⊥ ) basis choice. The fourth order coupling terms x x3

⊥ and x3 x⊥ are negligible at finite bias current and go to zero at zero bias current. We have also the frequency hierarchy g 21 , g 12 , g 22 ≪ ω , ω ⊥ .

The LM alone was previously studied [10,[15][16][17]. Quantum state detection was realized using a nanosecond pulse to produce a state dependent escape rate, through the potential barrier. At zero-current bias, the LM mode is insensitive to current noise [10]. Nevertheless, away from the zero-flux bias point, it is still flux noise sensitive as other phase qubits [18]. Until now, the TM has not eff . The area of maximum entanglement between these two states is found inside the dotted gray rectangle.

been studied and the resonant dynamics of the dc-SQUID have been restricted to one dimension only. Indeed the plasma frequency of the TM strongly increases with b. In previous dc-SQUID studies, the loop inductance was significantly smaller than the Josephson inductance, i.e b > 1, leading to a high TM plasma frequency [10,15]. By making the loop inductance larger than the Josephson inductance one can decrease the TM frequency and make it closer to the LM one, as shown in Fig. 2. The TM presents several interesting theoretical features that demand to be experimentally investigated and exploited. For instance, current fluctuations induced by the environment lead to a TM relaxation rate Γ r = g I (θ)S I (ω ⊥ ) where S I (ω ⊥ ) is the spectral density of the current noise and g I (θ) the coupling strength between current and the TM. For a bias current I b = 0, this coupling goes to zero for the TM. Thus, a longer relaxation time is expected. Moreover, the dephasing time depends on ∂ω ⊥ ∂I b , which goes to zero at I b = 0, and on ∂ω ⊥ ∂Φ b which is one order of magnitude smaller than for the LM (see Fig. 2). Therefore, the TM is expected to suffer from a smaller dephasing than the LM.

In addition, the coupling between the TM and the LM leads to a number of physical effects, expanding the range of applicability of SC devices in quantum information field. Hereafter we will concentrate our study on these coupling effects. The Hamiltonian of the uncoupled degrees of freedom is given by

H 0 = H + H ⊥ .
The energies E n ,n ⊥ = E n + E n ⊥ of its eigenstates |n , n ⊥ = |n ⊗ |n ⊥ are shown in Fig. 2 as dashed lines, for a particular set of circuit parameters. Since g 22 , g 12 , g 21 ≪ ω ⊥ , ω , the interaction term C is treated as a perturbation. To perform a detailed analysis, we introduce the annihilation operators â and b of the LM and TM, and the occupation number operators n = â † â and n⊥ = b †b . One then obtains

C = 1 2 √ 2 g 21 (â 2 + â †2 + 2n + 1) × ( b + b † ) + 1 2 √ 2 g 12 ( b2 + b †2 + 2n ⊥ + 1) × (â + â † ) -1 4 g 22 (â 2 + â †2 + 2n + 1) × ( b2 + b †2 + 2n ⊥ + 1
). In the interaction picture, the interaction Hamiltonian is H int = UCU † , with U = e -i(H +H ⊥ )t . The effect of the coupling is captured by the matrix

elements n ′ , n ′ ⊥ |H int |n , n ⊥ = e -iǫt n ′ , n ′ ⊥ |C|n , n ⊥ , where ǫ = E n ′ -E n + E n ′ ⊥ -E n ⊥ .
) eff = 1 2 √ 2 g 21 ( bâ †2 + b † â2 ).
The first term is present independently of the values of n and n ⊥ . It induces a conditional energy shift for both modes. The second term leads to a single desexcitation of the TM and the simultaneous double excitation of the LM (and the hermitian conjugated process). It is relevant only when a quasi-resonant condition is verified between states |n , n ⊥ and |n ∓ 2, n ⊥ ± 1 . The deviation from the resonant condition between these two states is quantified by

ǫ n ,n ⊥ (Φ b ) = (E n + E n ⊥ -E n -2 -E n ⊥ +1 )/ ≈ 2ω -ω ⊥ -(2n + 1)Λ ω + n ⊥ Λ ⊥ ω ⊥ .
One can choose a set of circuit parameters so that 2ω ≈ ω ⊥ . As shown in Fig. 2, tuning ω and ω ⊥ with the bias flux, we either satisfy the exact resonant condition (ǫ n ,n ⊥ (Φ b ) = 0), or ensure an off-resonant condition (|ǫ n ,n ⊥ (Φ b )| ≫ g 21 ).

We now study the rich physics associated with H (1) eff and H

(2) eff . Depending on the experimental parameters, each of these two terms can be chosen to dominate the coupling Hamiltonian.

We start with H

(1) eff , setting |ǫ n ,n ⊥ (Φ b )| ≫ g 21 , so that H (2)
eff can be neglected. H (1) eff is very similar to the Kerr effect Hamiltonian [START_REF] Butcher | The Elements of Non-linear Optics[END_REF], with typical values of g 22 /2π ranging from 50MHz to 500MHz. It allows i) the realization of quantum non-demolition (QND) measurement of the LM or TM [START_REF] Imoto | [END_REF], ii) conditional quantum logic (CQL) gates [21] and iii) the realization of a V-type 3-level atom, introduced initially in atomic-CQED. For example, QND measurement of the LM can be realized by TM spectroscopy:

H (1)
eff shifts in frequency the transition between the H 0 eigenstates |n , 0 ⊥ and |n , 1 ⊥ proportionally to (2n + 1). Using the two lowest energy states of both modes as qubits, this frequency shift can also be used to implement CQL gates. A microwave π-pulse resonant . The blue continuous line shows the population of state |0 , 1 ⊥ , set equal to 1 at t = 0, and the pink dashed line is the population of state |2 , 0 ⊥ , initially equal to 0. The inset shows the absolute value of the negativity as a function of time, providing evidence that maximally entangled states between the two modes can be created.

with the |0 , 1 ⊥ ↔ |1 , 1 ⊥ longitudinal transition can switch both states without affecting |0 , 0 ⊥ and |1 , 0 ⊥ . Therefore the qubit defined by the TM acts as a control qubit in a CNOT gate. Such gate have been successfully realized with flux qubit circuit [22], but away from the optimal point. In the proposed dc-SQUID, such operations can be realized at the optimal point (I b = 0,Φ b = 0), insensitive to flux and current noise in first order. Finally, the three states |0 , 0 ⊥ , |1 , 0 ⊥ and |0 , 1 ⊥ realize a V-like level structure, that was also proposed recently in charge qubit [26]. As proposed in Ref. [11], inserting such a V-type atom in a microwave resonator opens avenues for a wealth of effects discussed up to now in the context of atomic physics, such as quantum beats [27], the realization of dark states and the implementation of non-linear quantum field Hamiltonians [28].

We now discuss H

eff in the case g 22 ≪ g 21 , so that H (1) eff can be neglected. Experimentally, this condition is satisfied in a dc-SQUID with the parameters presented in Fig. 2. In this case, H

eff couples states |n , n ⊥ to states |n ∓ 2, n ⊥ ± 1 when |ǫ n ,n ⊥ (Φ b )| ≪ g 21 . In the following, we show that H (2) (2)
eff can be used to i) create the Fock state |2 in the longitudinal mode, starting from a single excitation quanta in the transverse mode, and ii) entangle the states of the two oscillation modes.

To illustrate the first point, we consider the subspace spanned by the pair of states |2 , 0 ⊥ and |0 , 1 ⊥ . The dependences of the associated energies on the bias flus are shown in Fig. 2 as solid blue curves. These two states satisfy the resonance condition at a bias flux Φ r = 0.54Φ 0 , and H

(2) eff then results in an avoided crossing (inside the gray rectangle). Away from the resonance, the coupling is negligible and the system eigenstates are those of the uncoupled Hamiltonian H 0 . Initially, Φ b is set to 0 and a π-pulse drives the TM from |0 , 0 ⊥ to |0 , 1 ⊥ (this state has a lower energy than |2 , 0 ⊥ ). Increasing Φ b slowly enough, it is then possible to realize an adiabatic passage through the avoided crossing. This results in the conversion of the state |0 , 1 ⊥ into the state |2 , 0 ⊥ . To precise the adiabaticity condition, we suppose that in the vicinity of the resonance, the frequency difference between the states varies linearly in time: ǫ n ,n ⊥ (Φ b ) = λt. The adiabatic regime is obtained for a frequency sweep parameter λ ≪ (π/2)g 2 21 . This process is illustrated in Fig. 3 using realistic experimental parameters. The initial frequency splitting between the two involved states is ǫ 2 ,0 ⊥ (0)/2π ≈ 5 GHz and g 21 /2π ≈ 600 MHz close to the level anti-crossing. With λ = 0.69g 2 21 , corresponding to a variation of 1.8GHz per nanosecond, one can perfectly transfer the population from |0 , 1 ⊥ to |2 , 0 ⊥ . Experimentally, this frequency sweeping speed can be realized using a pulse generator with a typical rise-time of 5 ns. This time should be compared to the coherence damping time T 1 . For a single mode dc-SQUID in the state |1 , T 1 ≈ 100 ns [10]. Since we expect the TM to present longer values of T 1 , these relaxation times are more than enough to perform the state transfer protocol.

Here, we must stress that state transfer is made possible by the adiabatic evolution of Φ b . In the instantaneous limit (λ ≫ (π/2)g 2 21 ), if the system is initially prepared in |0 , 1 ⊥ , for instance, it does not have the time to couple to state |2 , 0 ⊥ during the fast evolution. Consequently, coupling can be disregarded during the flux time variation. Such a condition is met using a pulse generator with a typical 500 ps rise-time.

Inspecting again Fig. 3 around the resonance, one observes strong oscillations in the population of |0 , 1 ⊥ . As we approach resonance, |0 , 1 ⊥ is no longer an eigenstate of H 0 . Interestingly enough, the system eigenstates are in fact entangled states of the two oscillation modes. The inset of Fig. 3 shows the time evolution of entanglement, quantified by the negativity N [23]. During the adiabatic evolution, the system passes through a maximally entangled state, evidenced by N = 1/2. This consideration, together with the time evolution regimes explored in the previous paragraphs, suggest a procedure to create and detect controlled entanglement between the longitudinal and transverse modes. To prepare a maximally entangled state, for instance, one can stop the adiabatic evolution at the time corresponding to N = 1/2. By 'instantaneously' restoring the system to the non-resonant regime, it remains a maximally entangled state. One can then measure this state and its entanglement by measuring the population and coherences of the TM and the LM. The same type of analysis can be extended to other, more excited, TM-LM states, where coupling between modes can also occur depending on the value of ǫ n ,n ⊥ (Φ b ).

The double excitation of one mode and the simultane-ous single desexcitation of another mode (and the hermitian process), has been experimentally observed in cavity QED systems [24]. In this system, the two modes of a high quality microwave cavity (formally equivalent to two harmonic oscillators) interact, mediated by an atom, via a third order process in perturbation theory. In trapped ion systems, such Hamiltonians can also be produced by tuning a laser to specific sideband frequencies corresponding to the excitation of the two-dimensional vibrational motion [25]. Usually, the coupling constant is proportional to the fourth power of the Lamb-Dicke parameter, so the coupling strength is significantly reduced. In the present system the non linear coupling is observed in the strong coupling limit, i.e. g 21 ≫ 2πT -1 1,2 . In conclusion, we have shown that a dc-SQUID is an artificial atom with two coupled degrees of freedom, the LM and TM. This coupling is fundamentally interesting and can be described by a number of effective Hamiltonians whose strengths vary according to experimentally tunable parameters. Each Hamiltonian gives rise to a number of applications, such as QND measurement, quantum logic, entanglement creation and the introduction of a more complex level structure to a superconducting artificial atom. Further applications arise when we consider coupling this system to a microwave cavity, opening the way to generalizing a larger collection of results already existing in cavity QED.

FIG. 1 :

 1 FIG. 1: (a) Electrical scheme of an inductive dc-SQUID (b) 2D potential for a dc-SQUID with b = 3, I b = 0 and Φ b = 0.51Φ0; M and M ′ indicate two different local minima separated by a saddle point C. The dashed black line represents the trajectory of minimum energy. (c) potential of the LM (left side) and TM (right side), related to the motion of the fictitious particle trapped at the local minimum M .

FIG. 2 :

 2 FIG. 2: Energy spectrum. Transition frequencies between ground state |0 , 0 ⊥ and excited states |n , n ⊥ versus Φ b at zero current bias of the artificial atom described by two oscillators without coupling (dashed lines) and taking into account the coupling (continuous lines) with I0 = 710 nA, C0 = 510 fF, b = 0.73. The two solid blue curves highlight the states |0 , 1 ⊥ and |2 , 0 ⊥ , chosen to illustrate the effect of the coupling term H (2)

21 NegativityFIG. 3 :

 213 FIG.3:(color online) Population transfer between states of the TM and the LM as a function of time for λ = 0.69g2 21 . Time is in units of g -121 . The blue continuous line shows the population of state |0 , 1 ⊥ , set equal to 1 at t = 0, and the pink dashed line is the population of state |2 , 0 ⊥ , initially equal to 0. The inset shows the absolute value of the negativity as a function of time, providing evidence that maximally entangled states between the two modes can be created.

  In order to couple both modes in a non negligible way, n ′ , n ′ ⊥ |H int |n , n ⊥ must oscillate at a frequency ǫ much slower than the coupling strength between the two states. Only two terms can potentially verify this condition in C: H

	and H	(2	(1) eff = 1 4 g 22 2n + 1 × (2n ⊥ + 1)

We thank W. Guichard and F. W. J. Hekking for fruitful discussions. This work was supported by the european EuroSQIP and SOLID projects and by the french ANR "QUANTJO".