Instability of Holographic Superfuids in Optical Lattice

The instability of superfluids in optical lattice has been investigated using the holographic model. The static and steady flow solutions are numerically obtained from the static equations of motion and the solutions are described as Bloch waves with different Bloch wave vector $k$. Based on these Bloch waves, the instability is investigated at two levels. At the linear perturbation level, we show that there is a critical $k_{c}$ above which the superflow is unstable. At the fully nonlinear level, the intermediate state and final state of unstable superflow are identified through numerical simulation of the full equations of motion. The results show that during the time evolution, the unstable superflow will undergo a chaotic state, where solitons will appear and disappear. The system will settle down to a stable state with $k


I. INTRODUCTION
As an artificial lattice structure, optical lattice, has a close physical resemblance to the periodic coulomb potential that felt by electrons in solid crystal. While it is not like electrons moving between positive ions in a nanometer dimension, optical lattice can be manipulated with a typical dimension much larger than the conventional crystal-micrometer dimension. Due to this virtue, optical lattice has provided a broad platform for experimental and theoretical physicists to explore much richer and more fundamental properties in condensed matter physics, such as instability [1,2], Landau-Zener tunneling [2][3][4][5][6][7], superfluid-Mott insulator quantum phase transition [8][9][10][11][12][13] and so on. Additionally, by controlling the frequency of the lasers in the course of time, the time-dependent optical lattice can be used for driving the system to explore the floquet dynamics [14][15][16][17][18], which are getting more interesting recently.
Since the first theoretical study about the Landau and dynamical instability for superfluids in optical lattice [1] using the Gross-Pitaevskii (GP) equation, this topic has been studied extensively. [19,20] explore the cases both for repulsive and attractive atomic interactions, [29] adds the influence of three-body interactions. [30] shows that the states in excited bands always have Landau instability and [23] investigates the effect of transverse excitations for higher dimensional systems. [24] reveals that the dynamical instability is the origin of the spatial domain formation in spin-1 atomic condensates. All these research works relied on the mean-field theory, i.e., GP equation, while in addition to this equation the Bose-Hubbard model is also used in [31,32] that is beyond the mean-field theory.
In the above research works, all the superfluid systems have no dissipation. From the instability diagram ( Fig.5 in [1]) we can see that there are two distinct regimes about the instability of superflow. The first is the small optical lattice height regime and in this regime the dynamical instability regime is so narrow that the Landau instability can be studied solely; the second is the large optical lattice height regime where the dynamical instability regime spreads out and it can be detected in experiments, obviously. Experiment [33] shows that the regime where the dissipative process occurs to break down the superfluid phase agrees with the theoretical prediction about the regime of onset of Landau instability for a superfluid in shallow optical lattice; Experiment [34] finds that for deep optical lattice there will be transition from a superfluid into an insulator when the superfluid is under dynamical instability. After that, experiments [35,36] both show that in the presence of thermal component the dissipative process can be related with the occurrence of Landau instability. Since the original GP equation cannot apply to dissipative processes, until [37] based on the modified GP equation with dissipation, the instability of superfluid in optical lattice with thermal component has not been studied theoretically. Its result dramatically shows that when the thermal component is added the dynamical instability will occur throughout the Landau instability regime predicted from the original GP equation, which means the Landau instability and dynamical instability will always appear at the same time when there exists dissipation. So it is no longer suitable to investigate Landau instability or dynamical instability separately.
In this paper, we will use the simplest holographic superfluid model [49,50] to re-explore the Landau and dynamical instability of superfluids in optical lattice simultaneously based on the advantage that such a model contains superfluid and normal fluid (as thermal component with a finite temperature) intrinsically, which introduces the dissipation naturally and consistently [51].
In contrast, the dissipation in the modified GP equation is purely put by hand, without any relation to thermodynamics. In holography, the unstable region can be calculated by quasi-normal mode (QNM) analysis and the results show that there is a critical Bloch wave vector k c above which the Bloch states are unstable due to the Landau instability and dynamical instability. With the help of the calculation of sound speed we confirm that the Landau instability and dynamical instability also appear at the same time in the presence of optical lattice just like the homogenous case in [52,53].
Beyond the linear perturbation analysis, the unstable superflow are also studied by the numerical simulation. In this part we use the unstable superflow state as the initial state for the evolution equation and evolve it for a long time with some small perturbation given at early time. The unstable modes of perturbation will grow exponentially at first and it will lead the whole system into a chaotic state, in which the solitons will appear and disappear (for higher dimension there will be vortices and the system can be considered as under a transient turbulence state [53]). Finally, along with some dissipative processes the system will reduce its current to become stable with the final wave vector k < k c . These results can be tested in experiments. This paper is organized as follows. In Sec. II we introduce the holographic superfluid model that we use and show how the optical lattice is added. Then in Sec. III by solving the corresponding equations of motion we get the static and steady-flow suferfluid states that expressed as Bloch waves with different wave vector k. The QNM analysis and the calculation of sound speed are in Sec. IV. And in Sec. V we present the dynamic processes of the evolution of an unstable superflow with numerical simulation. Finally, in Sec. VI we give a summary.

II. HOLOGRAPHIC MODEL
The simplest holographic model to describe superfluids is given in [49], where a complex scalar field Ψ is coupled to a U (1) gauge field A M in the (3 + 1)D gravity with a cosmological constant related to the AdS radius as Λ = −3/L 2 . The action is where G 4 is the gravitational constant in four dimensional spacetime. The first two terms in the parenthesis are the gravitational part of the Lagrangian with the Ricci scalar R and the AdS radius L, and the third consists of all matter fields: where Since the backreaction of matter fields onto the spacetime geometry is not necessary for our problem, taking the probe limit e → ∞ is a suitable and convenient choice, which means that we fix the background spacetime as the standard Schwarzschild-AdS black brane with metric Here z H is the radius of the black brane horizon. The Hawking temperature of this black brane is T = 3 4πz H , which is also the temperature of the boundary system by the holographic dictionary. Furthermore, due to the existence of the bulk black hole the boundary field theory will intrinsically have dissipation, since there will be energy flow absorbed into the horizon [51,54,55] during dynamic processes. For numerical simplicity, we set L = 1 = z H .
The equations of motion for A µ and Ψ are One can easily see that there is a trivial solution with Ψ = 0. As we increase the chemical potential µ, there will be a critical chemical potential µ c above which a solution with nonzero Ψ appears, indicating the break of U (1) symmetry and the formation of superfluid condensate. Hereafter, we will focus on the case with µ = 4.5 > µ c , i.e., the superfluid phase of the boundary field theory.

A. Equations of motion and optical lattice
Since the considered problem is based on putting superflow onto the optical lattice, the density of the superflow will be periodically modulated by the periodic potential (chemical potential plus external potential). For simplicity and without loss of generality, we take the axial gauge A z = 0 as usual and consider our bulk fields as functions of (t, z, x) with the assumption that the direction of the optical lattice is along x. It turns out that A y can be turned off in this case, so the with d = 3 in our case. The holographic dictionary tells us that with one of Ψ ± being the source the other will be the related response, a t is the total potential, ρ the (conserved) particle number density and a x the source of the particle current density j x in the boundary field theory. The optical lattice is then introduced by including a potential V (x) in a t , i.e. a t = µ + V (x). It is convenient to choose m 2 = −2 to make all calculations easier, in which and we choose Ψ + as source and Ψ − as the response. Here we introduce the bulk field ψ for numerical simplicity.

III. BLOCH WAVES AND BLOCH BAND IN OPTICAL LATTICE
From the Bloch theorem and with the source Ψ + turning off, we can use Bloch waves to describe the static response Here, ψ B (z, x) is the Bloch wave along the x direction, which has the same period as the optical , k is the Bloch wave vector and l the lattice constant (hereafter we choose l = π). Due to the fact that Bloch waves are complex functions, we can separate them into real and imaginary parts: Substituting the above equation into the static equations of motion (10)-(13), we get the following five differential equations: Among them, Eq.
the height (or strength) of the optical lattice. We solve these four static equations of motion (20), (22), (23) and (24) numerically with the Newton-Raphson method. Fig.1 shows the numerical solutions of the bulk fields ψ and A t , respectively, with two different Bloch wave vectors. Since ψ is complex, we plot its amplitude and phase angle separately. The inexistence of node of order parameter means that these Bloch waves sit in the lowest Bloch band.
In Fig. 2a, we plot the Bloch band in our holographic model. To obtain it, we fix the total particle number, defined as and solve the static equations (20) and (22)-(24) with k given different values. The loop structure of the blue line (v = 0.5), at which comes from particle interaction energy is larger than the optical lattice strength [56], is witnessed around the edge k = 1 of the Brillouin zone. And for the red line (v = 4) there is no loop structure. Here, the Bloch band also confirms that the static solution in Fig. 1 are Bloch wave solutions. Actually, the higher value of chemical potential µ, the larger density of total particle ρ in one lattice and hence the larger value of particle interaction energy, which means that when we increase the chemical potential the critical value of the optical lattice strength v c that exceeds the particle interaction energy to prevent the loop structure from appearing will also increase. And we confirm this in Fig. 2b.

IV. INSTABILITY AND SOUND MODES
Now we study the linear instability of the Bloch wave solutions obtained in the previous section via QNM analyses. To calculate the QNM, a better choice is to rewrite the static equations of motions (10)-(13) with the Bloch wave form for ψ in the infalling Eddington-Finkelstein coordinates [57]: where the metric is Here and in the following sections we will not separate ψ into real and imaginary parts. Similarly, there is one constraint equation, which is chosen to be Eq. (27). Next, we give all the bulk field Substitute these perturbed fields into the equations of motion, we will get the linear perturbation In accordance with the background steady solutions, we impose Dirichlet and periodic boundary The perturbations in Eq.(30) will grow exponentially in the linear regime if there exists an ω whose imaginary part is positive. We use ω M to denote the ω with the maximal imaginary part when q and µ are given, so a system is (linearly) unstable if Im(ω M ) > 0. Fig. 3 shows the result for the distributions of ω with two different values for Bloch wave vector k and q. In Fig. 3a and q M means the Bloch wave vector of the most unstable perturbation mode, which has a maximal imaginary value of ω M , since it will increase with time t with the form e ω I t controlled by the linear perturbation equations in early times (the late time behavior will be more interesting and it will be investigated in the next section.). Fig. 4 shows the value of ω M with respect to q ∈ [0, 4], the periodicity and q M can be seen directly.
As we have explained, when q and k are given, ω M can be obtained to determine the stability under the corresponding perturbation. In this spirit, we plot stable and unstable regions as a function of q and k, i.e. the instability diagram, in the half first Brillouin zone in Fig. 5. In this plot, the left region in the figure, which is colored darker grey, is stable (Im(ω M ) ≤ 0), and the right region in the figure, which is colored lighter grey, is unstable (Im(ω M ) > 0). 1 The figure also tells us that in the unstable region the cutoff q increases with k, and when k is large enough (but still in the first Brillouin zone), the system is unstable for all values of q except q = 0.
As a signal of onset of the instability mentioned in [52,53], we also plot the dispersion relation of the sound mode in Fig.6. When there is a non-zero k the sound speed becomes direction-dependent.
There are two directions corresponding to the maximal and minimal values of the sound speed, which are parallel (q > 0) and anti-parallel (q < 0) to the velocity of the superflow, respectively. Re(ω(q)) and right panel shows Im(ω(q)). A universal phenomenon is that Im(ω(q)) > 0 always accompanies Re(ω(q)) < 0.
Hereafter we denote the sound speed as c + for q > 0 and c − for q < 0.
The onset of instability is accompanied by the sign-change of c − . Along the negative axis of q, Re(ω) becomes negative as k is larger than 0.655. As Re(ω) corresponds to the energy of the perturbation, c − < 0 (so Re(ω) < 0) indicates the existence of a perturbed state with lower energy and thus the instability of the system (Landau instability). At the same time, we find that Im(ω) > 0, which implies the dynamical instability of the system. So, our result confirms that the Landau instability and dynamical instability occur at the same time, as we have mentioned.
For the purpose of opening the next section we give a brief discussion about this section. Since it is not always true for a superflow to flow in the optical lattice stably, when k ≥ k c excitations will always appear to slow the superflow down [53], but what kind of intermediate states will this unstable superflow go through and what kind of final states will it end up with? In the next section we will use full nonlinear evolution to find the final state out for the Bloch wave-type superflow and the intermediate states will also be discussed.

V. REAL TIME EVOLUTION FOR STABLE AND UNSTABLE SUPERFLOWS
To simulate the evolution of the holographic superfluid, we adopt the same scheme as in [57], i.e. taking Eqs. (26), (28) and (29)

A. Evolution of particle current density
We introduce perturbations at t = 100, i.e. t 1 = 100, and observe the evolution of the systems with initial Bloch vectors k chosen as 0.3 and 0.75, respectively. Fig. 7 shows the evolution of the particle current density j x , which is defined as From Fig. 7, the evolution of j x of an unstable superflow can be divided into five stages. In the first stage t < 100, no change is observed, which confirms that the static solutions obtained in the Schwarzschild coordinates is correctly transformed to that in the Eddington-Finkelstein coordinates. When 100 < t < 200, the current density changes slightly since the influence of the perturbation is small; when 200 < t < 400, the current density becomes chaotic, which comes from the nonlinear development of the instability. When 400 < t < 1500, the chaotic behavior of the current density disappears but the current is still inhomogeneous; In this stage, the system is tending to a steady state. In the last stage t > 1500, the system becomes steady with the final current density j x f = 0.0931468, which is much smaller than the initial value 0.26048. In contrast, for the stable superflow (k = 0.3 as an example for comparison) the situation is totally different. The influence of perturbations on the system is insignificant and the system will evolve to a state that is same as the initial state in a few moments.

B. Evolution of condensate and intermediate states
As we have seen from the evolution of j x , the intermediate states of the system are chaotic and complicated. However, there are some universal features which can be extracted from the evolution of the condensate |Ψ − | and we plot them in Fig. 8.
We plot the evolution of |Ψ − | from different viewing angles, i.e. Fig. 8a, Fig. 8b and Fig. 8c. We can see that the condensate during 100 < t < 400 is also chaotic and the value of |Ψ − | at later times is larger than its initial value. The planform of |Ψ − | is plotted at Fig. 8c and nodes of the order parameter appear during 100 < t < 400, which indicates the formation of solitons. We select one moment (t = 204) of the soliton formation in Fig. 8d. The formation of solitons is natural here [53] since the unstable superflow is thought of as having a higher free energy than a less unstable superflow with soliton excitations. Along with the dissipative processes, the whole system will go to a stable state with a lower free energy.

C. Final stable state
Since both j x and |Ψ − | become time independent at the end of the evolution, the final state of the system must be able to be described as a Bloch state just like the initial steady flow states that are solved in Sec. III. This is indeed the case. Actually, the final state can also be solved When there is a one-to-one correspondence between j x and the Bloch wave vector k, one can get the additional boundary condition directly by fixing j x = j x f . However, as we can find from Fig. 9a, i.e. the relation between j x and k, a fixed j x normally corresponds to two values of k.
To make the solution unique, we can take account of the average condensate |Ψ − | in one lattice cell. From Fig. 8a we know that the condensate amplitude |Ψ − | of a steady state is not a constant, so the relation between |Ψ − | and k is not clear, while |Ψ − | is a monotonically decreasing function of k, as plotted in Fig. 9b. Therefore, the boundary condition can be uniquely determined by combining the two considerations: j x = j x f and the value of |Ψ − | . From Fig. 8b and Fig. 7c we conclude that the final Bloch state corresponds to a larger |Ψ − | with j x = 0.0931468, thus we The non-monotonicity for j x gives two states with different k, while the monotonicity for |Ψ − | excludes the larger k state.
can solve the equations of motion with the additional boundary condition ∂ z A x | z=0 = 0.0931468, and the solution is marked on Fig. 9b.

VI. SUMMARY
The dynamical process of superfluids that simultaneously experiences Landau and dynamical instability in optical lattice is investigated using the simplest holographic superfluid model. Due to the existence of dissipation the sound mode of an unstable superflow will always have negative energy and grows up exponentially at early times. We give a universal picture of the instability process for the unstable superflow. From real time evolution the unstable superflow is shown to reduce its particle current density to settle down to a steady state with a current density j x f and a Bloch wave vector k f below the critical value k c , where this k f can be solved from the equations of motion by fixing j x = j x f . In the course of the evolution, a chaotic process involving soliton generation is observed, in which the solitons are supposed to play the role of reducing the wave vector k (similar to the case without an optical lattice [53]).
As we have mentioned, when the number of dimensions of the boundary system is greater than one, the process for an unstable state to evolve into a stable one will become much more complicated, which to some extent is described as transient turbulence [53]. So beside the one dimensional boundary superfluid system in optical lattice, it is worthwhile to consider higher dimensional cases. Another direction for future investigation is to include backreaction in this holographic superfluid model, which will enable us to explore the complete interplay between the normal fluid and superfluid components of the superflow at a finite temperature in optical lattice.

ACKNOWLEDGMENTS
We are grateful to Biao Wu  x) ∼ f (t)e ikx with continuous value of k; if there is a discrete space translation symmetry, then we also have F (t, x) ∼ f (t)e ikx but with discrete value of k = 2π L n, where n is integer and L is lattice constant.
But in our case, we do not use the above separation though there is a discrete space translation symmetry but use Bloch wave to describe the x direction of the fields, i.e. F (t, x) = f (t, x)e iqx .
Here F(t,x) and f(t,x) can be complex functions and we can do some deformations, With this deformation we can write (30) and get (31)-(34) straightforward.
Equations (31)- (34) can be written as a matrix form, i.e., with appropriate boundary conditions given in Section IV. These function is indeed the general eigenvalue function and there are many methods in textbook to solve this function to get the eigenvalue ω.
Before we calculate the value of ω we can know that if ω R is real part of one of the eigenvalue ω, then there must be another value of ω whose real part equals to -ω R . To show this we can define a matrix L, which is defined as then Eq.(A2) can be written as Separate ω into real part and imaginary pary and we get that the imaginary part of L only contains ω R , which means From Equations (31)- (34) it is easy to get that when q = 0 the following relation exists, i.e., Since the boundary conditions for u and v are the same, changing the places between u and v do not affect the eigenvalues. While when q = 0, there do not have the relation (A6).