Multiple period states of the superfluid Fermi gas in an optical lattice

We study multiple period states of a two-component unpolarized superfluid Fermi gas in an optical lattice along the Bardeen-Cooper-Schrieffer (BCS) to Bose-Einstein condensate (BEC) crossover. The existence of states whose period is a multiple of the lattice spacing is a direct consequence of the non-linear behavior of the gas, which is due to the presence of the order parameter associated with superfluidity. By solving Bogoliubov-de Gennes equations for a superfluid flow with finite quasimomentum, we find that, in the BCS side of the crossover, the multiple period states can be energetically favorable compared to the normal Bloch states and their survival time against dynamical instability drastically increases, suggesting that these states can be accessible in current experiments, in sharp contrast to the situation in BECs.

Density structures and patterns caused by the interplay of competing effects are ubiquitous in nature.The competition between dispersion and non-linearity yields solitons [1] and the interplay of the original periodicity of the crystalline lattice and the Peierls instability of conduction electrons results in charge and spin density waves [2].A layered density pattern in supercritical flow of 4 He in a narrow channel is predicted to appear due to the interplay of effects of the Landau roton instability and the interaction between rotons [3].Nuclear halo [4] is due to the interplay between the effects of the quantum mechanical spreading of the wavefunction and the attractive nature of the nuclear force.Even in astrophysical environments such as in neutron stars and supernovae, the competition between the Coulomb energy and the surface energy of nuclei causes the so-called "pasta" phases, which consist of crystalline lattice of rod-like and slablike nuclei ( [5] and references therein).
In the case of superfluids in a periodic lattice, nonlinearity due to the presence of the order parameter of the superfluid phase can give rise to non-trivial phenomena such as a loop structure called "swallowtail" in the energy band [6][7][8][9][10][11][12].Superfluidity and periodic lattices are basic elements of condensed matter physics.Density patterns caused by the subtle interplay of these two elements are hence an intriguing problem.Ultracold atomic gases in optical lattices are suitable to study these effects because of the high controllability of both the lattice geometry and the interatomic interaction [13,17].In addition, by using a Feshbach resonance to increase the interatomic attraction, one can continuously go from the Bardeen-Cooper-Schrieffer (BCS) to the Bose-Einstein condensate (BEC) regime [14,15], thus allowing one to study Bose and Fermi superfluids from a unified perspective.
For BECs in a periodic potential, it was found that non-linearity of the interaction can cause the appearance of stationary states whose period does not coincide with that of the lattice; instead, a multiple of it [16][17][18].Such states are called multiple (or n-tuple) period states.In BECs without long-range interaction, however, these states are energetically unfavorable compared to the normal Bloch states whose period is equal to the lattice constant, and the lowest multiple period states are dynamically unstable [16].
Period-multiplying phenomena are particularly interesting in Fermi superfluids because of possible wide implications in condensed matter physics and nuclear physics, such as for superconducting electrons in solids and superfluid neutrons in neutron stars [19,20].However, unlike the case of Bose gases, little has been studied about multiple period states in superfluid Fermi gases and their existence itself along the BCS-BEC crossover is still an open question.In this work, we predict that multiple period states indeed exist in Fermi superfluids and we determine their energetics and dynamical stability along the BCS-BEC crossover, pointing out peculiar features of fermions which cannot be seen in the case of bosons.
We consider a two-component superfluid Fermi gas in the BCS-BEC crossover flowing through a onedimensional (1D) optical lattice, where q B ≡ π/d is the Bragg wave vector, d is the lattice constant of the primitive cell of the periodic potential, E R ≡ 2 q 2 B /2m is the recoil energy, and m is the atom mass.The gas is uniform in the transverse directions.We look for stationary states of the system by numerically solving the Bogoliubov-de Gennes (BdG) equations: where H(r) = − 2 ∇ 2 /2m + V ext (r) − µ, u i (r) and v i (r) are quasiparticle amplitudes, and ǫ i is the corresponding quasiparticle energy.The chemical potential µ is determined from the constraint on the particle number N = 2 i |v i (r)| 2 dr and the pairing field ∆(r) should satisfy a self-consistency condition ∆(r) = −g i u i (r)v * i (r), where g is the coupling constant for the s-wave contact interaction which needs to be renormalized.
In the presence of a supercurrent with quasimomentum P (Q ≡ P/ ) per atom moving in the z-direction, one can write the quasiparticle amplitudes in the Bloch form as u i (r) = ũi (z)e iQz e ik•r and v i (r) = ṽi (z)e −iQz e ik•r leading to the pairing field as ∆(r) = e i2Qz ∆(z).Here ∆(z), ũi (z), and ṽi (z) are complex functions with period ν times d, with ν ∈ {1, 2, 3, ...}, and the wave vector k z lies in the first Brillouin zone for a supercell (a cell containing several primitive cells) with period ν, i.e., |k z | ≤ q B /ν.This Bloch decomposition transforms (2) into the following BdG equations for ũi (z) and ṽi (z): where Here, k 2 ⊥ ≡ k 2 x + k 2 y and the label i represents the wave vector k as well as the band index.We solve this BdG equations (3) for a supercell with period ν (−νd/2 ≤ z ≤ νd/2) to obtain the period-ν states.
In the following, we mainly present the results for s = V 0 /E R = 1 and 2 with E F /E R = 0.25 as examples, where E F = 2 k 2 F /(2m) and k F = (3π 2 n 0 ) 1/3 are the Fermi energy and momentum of a uniform free Fermi gas of density n 0 .These values fall in the range of parameters of feasible experiments [21].We denote by P edge the quasimomentum P at the edge of the Brillouin zone for the normal Bloch states with period 1 (P edge ≡ q B /2 for superfluids of fermionic atoms).
In Figs.1(a) and 1(b), we show the profiles of the pairing field |∆(z)| and the number density n(z) of the lowest period-doubled states (period-2 states), respectively.Here we set P = P edge /2 = q B /4 at the Brillouin zone edge of the period-2 states.The nature of the period doubling and the difference between the regions of −1 < z/d ≤ 0 and 0 < z/d ≤ 1 can be clearly seen in |∆(z)| at any value of 1/k F a s .At P = P edge /2, |∆(z)| of the period-2 states has a node [see z/d = 0.5 in Fig. 1(a)] and consequently the supercurrent is zero (∂ P E = 0; see Fig. 2).On the other hand, especially in the deeper BCS side (1/k F a s = −1), the difference in n(z) between the regions of −1 < z/d ≤ 0 and 0 < z/d ≤ 1 is small [see the red line in Fig. 1 imbalance between the two regions becomes larger with increasing 1/k F a s toward the BEC regime.
It is instructive to consider the deep BCS limit of 1/k F a s → −∞.There, ∆(z) vanishes and n(z) of the neighboring sites becomes identical so that the nature of the period doubling disappears in this limit.We observe that by going to the deep BCS regime, where ∆(z) and the supercurrent are infinitesimally small, the energy difference E(P ) − E(0) for period-2 states decreases [i.e., E(P ) becomes more flat] and period-2 states at P = P edge /2 approach the normal Bloch state at P = 0.These observations are consistent with the fact that, if ∆(z) = 0, our model reduces to the Schrödinger equation, whose solutions have the periodicity of the lattice due to the Bloch theorem.Multiple period states are hence possible only in the superfluid phase.
In Fig. 2, we plot the energy per particle of the lowest band as a function of the quasimomentum P for the normal Bloch states (blue dotted line) and the multiple period states.We show the results at 1/k F a s = −1 in the BCS side.In the region of small P , all the lines of the multiple period states coincide with the line of the normal Bloch states.This is due to the fact that the multiple period states are equivalent to the normal Bloch states in this region: The states with period 1 are a subset of all the multiple period states.It is remarkable that the multiple period states appear energetically more stable than the normal Bloch states near the first Brillouin zone edge of each multiple period state, i.e., P P edge /ν for period-ν states.Particularly, the period-doubled states are energetically the most stable in a wide range of P in this case [22].This is in contrast to multiple period states in BECs and in the BEC regime of the BCS-BEC crossover (see later), where the energy of these states are higher than that of the lowest band of the normal Bloch states [16].While, in those cases, the lowest multiple period states appear as an upper branch of the swallowtail band structure (thus, ∂ 2 P E > 0) around the Brillouin zone edge of the respective multiple period states [16], , the green dashed line with × (period 3), the cyan dashed-dotted line with (period 4), and the magenta solid line with (period 5).The dotted lines at P/P edge = 1/3, 0.25, and 0.2 show the first Brillouin zone edge for the period-3, 4, and 5 states, respectively.The inset shows that the lowest band of the period-5 states continues beyond the first Brillouin zone edge and the swallowtail appears.The critical quasimomentum Pc for the pair-breaking instability of the normal Bloch state is 0.147P edge in this case.
Fig. 2 shows that the lowest band of the multiple period states in the BCS regime yields ∂ 2 P E < 0 around the Brillouin zone edge.
At first sight, this seems to imply a pathological situation for ν-period states with large ν, whose energy would be smaller than the energy of the normal Bloch states even in the limit of P → 0. However, this pathological situation is saved by the emergence of the swallowtail: The multiple period states with large ν continue being almost identical to the normal Bloch states, and keep their nearly quadratic dispersion around P ∼ 0 even beyond their first Brillouin zone edge at P edge /ν, which results in the swallowtail band structure for the period-ν system.In the case of Fig. 2, the swallowtail starts to appear at ν = 5 (see the inset of Fig. 2).
In Fig. 3, we show the total energy difference ∆E ≡ E 2 − E 1 between the normal Bloch states (E 1 ) and the period-2 states (E 2 ) at P = P edge /2 = q B /4 along the BCS-BEC crossover.As we have seen in Fig. 2, the period-doubled states are energetically more stable (i.e., ∆E < 0) in the BCS regime.With 1/k F a s increasing from the deep BCS limit, ∆E increases from a negative value and finally period-doubled states become higher in energy than normal Bloch states (i.e., ∆E > 0) in the BEC side.We also show the results in the BEC side obtained by solving the Gross-Pitaevskii (GP) equation for parameters corresponding to s = 1 and E F /E R = 0.25 (the green dashed line) [23].We note that ∆E of the GP results approaches zero from positive values with increasing 1/k F a s .This suggests that, before the BdG results (the red and blue solid lines) converge to the GP ones in the deep BEC regime, ∆E a maximum value.
Why the multiple period states can be energetically more stable than the normal Bloch states in the BCS regime can be physically understood as follows.For clarity, we consider a period-doubled state and a normal Bloch state at P = P edge /2.The key point is that ∆(z) and n(z) can behave in a different way.Suppose we have a normal Bloch state at P = P edge /2.Since |∆(z)| is exponentially small in the BCS regime, we can distort the order parameter ∆(z) to have a node, such as the one in the period-doubled state at P = P edge /2, with a small energy cost (per particle) up to the condensation energy On the other hand, making a node in ∆(z) kills the supercurrent j = V −1 ∂ P E, which yields a large gain of the kinetic energy (per particle) of the superfluid flow of order ∼ P 2 edge /m ∼ E R .Even if ∆(z) is distorted substantially to have a node, the original density distribution of the normal Bloch state is almost intact so that the increase of the kinetic energy and the potential energy due to the density variation is small.Therefore, the period-doubled state is energetically more stable than the normal Bloch state in the BCS regime [24].On the other hand, in the BEC limit, the density is directly connected to the order parameter as n(z) ∝ |∆(z)| 2 , and distorting the order parameter accompanies an increase of the kinetic and potential energies due to a large density variation.
Finally, we examine the dynamical stability of the  period-doubled states using numerical simulations based on the time-dependent BdG (TD-BdG) equations.Our TD-BdG simulations are performed for systems with the periodic boundary conditions, which are different from the ones with the Bloch-wave boundary conditions used for the above calculations of the stationary BdG equations.To mimic the setup for the stationary BdG solutions, we use large systems with L z = 32d to 64d in the z-direction [25].We construct the initial condition for the TD-BdG simulations using the stationary solution of the BdG equations by taking its quasiparticle amplitudes u i and v i whose quasimomentum k z is equal to a multiple of ∆k z .We integrate the TD-BdG equations using a 4-th order predictor-corrector method [26].
The right and middle panels of Fig. 4 show the time evolution of |∆(z)| of the period-doubled state at P = P edge /2 in the BEC side (1/k F a s = 0.5) and at unitarity, respectively.We see that |∆(z)| [and n(z)] does not keep its initial profile and large-amplitude oscillations triggered by the dynamical instability set on at t ≃ 55 /E R in the right panel and t ≃ 130 /E R in the middle panel.
We also notice that the TD-BdG simulations allow us to identify the spontaneously growing excitations which trigger the instability [27].
It is remarkable that the survival time τ surv of the period-doubled states until they are destroyed by the large-amplitude oscillations drastically increases as going toward the BCS side.In the left panel of Fig. 4, we show the time evolution of |∆(z)| of the period-doubled state at P = P edge /2 for 1/k F a s = −1.In this realization, the period-doubled state almost keeps its initial profile of |∆(z)| until t ≃ 900 /E R , and even longer for n(z) because only a small fraction of particles participate in the pairing in the BCS regime.
To further analyze the time scale of the deviation |∆(z, t)| − |∆(z, 0)| from the initial state, we take its spatial Fourier transform and look for modes with exponentially growing amplitudes |η(t)| = |η(0)| e γt .From a fit we extract the growth rate γ of the fastest growing mode [28].The result is shown by the black solid line in Fig. 5, which clearly shows the suppression of γ with decreasing 1/k F a s .In practice, the survival time τ surv of the period-doubled states depends on the accuracy of their initial preparation.We estimate τ surv with η(0)e γt ∼ 1, where η(0) is the relative amplitude of the initial perturbation with respect to |∆(t = 0)|.In Fig. 5, we show τ surv for three values of η(0).This result suggests that if the initial stationary state is prepared within an accuracy of order 1% or 0.1%, this state safely sustains for O(100) /E R to ∼ 1000 /E R in the BCS side corresponding to O(1)msec < τ surv 10msec for typical experimental parameters [21,30].
In summary, we have studied multiple period states, especially period-doubled states, of superfluid Fermi gases in an optical lattice and have found that they can be energetically more stable than the normal Bloch states and their survival time can be drastically enhanced in the BCS side.There, the multiple period nature distinctly appears in the pairing field, which could be observed by the fast magnetic field sweep technique [31,32].We hope our work will stimulate future experimental studies.

Momentum distribution
The period n-tupling studied in the present work is different from the charge density wave due to the nesting of the Fermi surface.In the following Fig. 7  plateau.Note that even though the system is in a periodic potential with nonzero s, the Fermi surface is almost spherical.Figure 8 in this Supplemental Material is the same as Fig. 7(a) there, but for a stronger periodic potential with s = 2.The plateau region of the momentum distribution n(k) is significantly smaller compared to that for s = 1 shown in Fig. 7(a).This is due to the formation of bound bosonic dimers induced by the stronger peri-odic potential [1][2][3].However, also in this case, n(k) is almost isotropic although it is more compressed in the k ⊥ -directions compared to the case of s = 1 [Fig.7(a)].

Comparison of the energy between period-2 and -3 states
In Fig. 9 in this Supplemental Material, we show the energy difference E 2 − E 3 between the period-2 and -3 states at the Brillouin zone edge of the period-3 states, P = P edge /3.For the parameter of Fig. 2 in the paper, E F /E R = 0.25, we see that E 2 − E 3 < 0 and thus the period-2 states are energetically lower than the period-3 states in the whole Brillouin zone of the latter (see Fig. 2 in the paper).This figure shows that the period-2 states are always energetically lower than the period-3 states in the region of 0.245 E F /E R 0.4 for s = 1 and 1/k F a s = −1.

Further comments about multiple periodicity
In principle, multiple period states with rational number periods are covered by our calculations with a supercell.Specifically, using a supercell with period ν, we can describe period-(ν/η) states, where ν and η are natural numbers.On the other hand, states with irrational number periods are excluded, which are beyond the scope of the present study.In our numerical calculations, which cover multiple period states with rational number periods, neither states with ν < η nor states whose period is incommensurate with the lattice constant appear as the lowest energy state.Therefore, it is probable that the multiple period states with irrational number periods could not be the energetically minimum states.

2 FIG. 2 :
FIG.2:(Color online) Energy E per particle in units of ER as a function of the quasimomentum P .Here we set s = 1, EF /ER = 0.25, and 1/kF as = −1.The normal Bloch states with period d are shown by the blue dotted line with • symbols, and the multiple period states are shown by the red solid line with + (period 2), the green dashed line with × (period 3), the cyan dashed-dotted line with (period 4), and the magenta solid line with (period 5).The dotted lines at P/P edge = 1/3, 0.25, and 0.2 show the first Brillouin zone edge for the period-3, 4, and 5 states, respectively.The inset shows that the lowest band of the period-5 states continues beyond the first Brillouin zone edge and the swallowtail appears.The critical quasimomentum Pc for the pair-breaking instability of the normal Bloch state is 0.147P edge in this case.

FIG. 3 :
FIG. 3: (Color online) Difference ∆E ≡ E2 − E1 of the total energy per particle in units of ER between the normal Bloch states (E1) and period-doubled states (E2) at P = P edge /2.The parameters we have used are s = 1, 2 and EF /ER = 0.25.The red solid line with + (s = 1) and blue solid line with × (s = 2) show the results obtained by solving the BdG equations and the green dashed line shows the results by the Gross-Pitaevskii equations for corresponding parameters.

FIG. 4 :
FIG. 4: (Color online) Time evolution of the magnitude of the pairing field |∆(z, t)|/|∆(z = 0, t = 0)| of the period-doubled state at P = P edge /2 for 1/kF as = −1 (left panel), 0 (middle panel), and 0.5 (right panel).Here, s = 1 and EF /ER = 0.25.Actual calculation has been done for Lz = 32d in the cases of 1/kF as = 0.5 and 0 and for Lz = 48d in the case of 1/kF as = −1; a part of the system is shown in the figure.

FIG. 6 :
FIG. 6: Profiles of (a) the pairing field |∆(z)| and (b) the density n(z) of the normal Bloch states (magenta dashed line for 1/kF as = −1 and cyan dashed-dotted line for 1/kF as = 0) and period-doubled states (red solid line for 1/kF as = −1 and blue dotted line for 1/kF as = 0) at P = P edge /2.Here we set s = 1 and EF /ER = 0.25 as in Fig. 1 in the paper.