Dynamical stability of the Holographic System with Two Competing Orders

We investigate the dynamical stability of the holographic system with two order parameters, which exhibits competition and coexistence of condensations. In the linear regime, we have developed the gauge dependent formalism to calculate the quasi-normal modes by gauge fixing, which turns out be considerably convenient. Furthermore, by giving different Gaussian wave packets as perturbations at the initial time, we numerically evolve the full nonlinear system until it arrives at the final equilibrium state. Our results show that the dynamical stability is consistent with the thermodynamical stability. Interestingly, the dynamical evolution, as well as the quasi-normal modes, shows that the relaxation time of this model is generically much longer than the simplest holographic system. We also find that the late time behavior can be well captured by the lowest lying quasi-normal modes except for the non-vanishing order towards the single ordered phase. To our knowledge, this exception is the first counter example to the general belief that the late time behavior towards a final stable state can be captured by the lowest lying quasi-normal modes. In particular, a double relation is found for this exception in certain cases.


Introduction
The correlated stability, namely, the equivalence between thermodynamical and dynamical stability in gravitational physics, is a long-standing problem and has recently received renewed attention due to AdS/CFT correspondence. AdS/CFT correspondence states that a gravitational system in the AdS bulk is dual to an ordinary laboratory type system without gravity on the boundary. In particular, the bulk black hole is dual to the boundary system at finite temperature. So it is tempting to conjecture that the correlated stability should hold for such holographic gravitational systems [1,2]. Actually there are plenty of holographic models available for such a test, where a lot of works have been done [3][4][5][6][7][8][9].
Of special interests, the authors of [10] investigated a holographic model with two charged scalar fields not coupled directly to each other, which are correlated to two order parameters in the dual field theory, in the four dimensional bulk AdS spacetime under the probe limit, which we call the BHMRS model. More than one order parameters may describe a more complicated system, which has been extensively studied in condensed matter physics (e.g., in multicomponent superconductors and superfluids) [11][12][13].
In [10], a new superconducting phase with two coexisting orders was found in a certain range of chemical potential in the dual theory, besides the ordinary phase and the phase with only one superconducting order. The coexisting phase is competitive between the two order parameters, or in other words, when one of the two scalar fields' condensation increases, the condensation of the other would decrease. Taking account of the backreaction, the BHMRS model was studied in detail in [14] and a much richer phase structure was found, which is determined by calculating and minimizing the free energy (grand potential) of each phase when there are more than one possible phases for a fixed set of thermodynamic parameters.
In this paper, we focus on the dynamical (in)stability of the BHMRS model in the probe limit. In order for comparison with the thermodynamical (in)stability, we first calculate the free energy of the system in the rigorous probe limit. Then, we investigate the quasinormal modes (QNMs) of the system. Since the QNM investigation only involves linearized equations of motion, it is a quick and reliable way to find out the phase structure of a model.
Actually, we propose a new method to calculate QNMs without using the gauge invariant forms in the coexisting phase. Our new method looks more convenient than other widely used methods, e.g., the determinant method proposed by [16] for which one should first find a maximal set of linearly independent numerical solutions of the linearized equations of motion. In particular, compared with the gauge invariant formalism, the boundary conditions and gauge fixing of the linearized equations of motion are treated in a very natural way in our method, which is remarkably adapted for comparison with the late time behaviors of the real time evolution.
In addition to describing the late time behaviors of the evolution, QNM can also be used to investigate the phase transitions. Although we do not show the figures in this paper, we have made sure that there will be the lowest lying quasi-normal mode passing through the real axis around the critical values where a stable phase converts to be unstable.
We investigate the real time dynamical (in)stability of the system by giving the system various strength of Gaussian wave packets as perturbations to different equilibrium states in the initial time and then numerically evolving the system to find out which phase is dynamically stable when there are more than one possible phases. We find that the phase whose free energy is the lowest in certain area of charge density is also dynamically stable, i.e., the dynamical stability is consistent with the thermodynamic stability in this holographic system.
Interestingly, we find that the relaxation time of the system is generically much longer than the simplest holographic system, both in the case of far from criticality. Actually, we first observe this phenomenon in the dynamical evolution, and later confirm it through careful identification of the corresponding QNM, which would otherwise be confused with the mode w = 0. It needs further investigation to see if there is some physical mechanism behind this phenomenon.
We also investigate the late time behaviors of the above evolution, in comparison with the QNMs as linear perturbations over the final stable configuration, and find the late time behaviors of the coexisting phase and the order parameter that will decay to zero in noncoexisting phases can be described well by the QNMs. Unexpectedly, the late time behaviors of the order parameters that will condense in non-coexisting phases are not consistent with their own QNMs, which is believed to be due to the nonlinear effects, which we will discuss in some detail. This paper is organized as follows. In the next section, we will review the aforementioned holographic model for our two band superconductor. Then in Section 3 we confirm the phase diagram presented in [10] by comparing the free energy for all possible static solutions in the canonical ensemble. As a warm-up, in the subsequent section we develop the gauge dependent formalism with the gauge fixed to work out the quasi-normal modes, which indicates the thermodynamically stable phases are also dynamically stable at the linear level. In Section 5, by evolving our system with different Gaussian wave packets as perturbations at the initial time we find the thermodynamically stable phases are dynamically stable at the fully non-linear level. Moreover, it is found that the late time behavior can be well captured by the lowest lying quasi-normal modes except for the non-vanishing order parameters towards the single ordered phases. The reason for this exception comes essentially from the non-linear effect, or put it in a physical way, the competition between the two bands such that the relaxation time for the non-vanishing order is generically elongated by the dying-out order. We conclude our paper in the end.

The Holographic Model
We will consider a holographic superconductor model in (3+1) dimensional AdS. The action is given by [10] Here L is the radius of AdS, κ 2 = 8πG with G the gravitational constant in the bulk, and the matter field Lagrangian is given by where D 1 = ∇ − i e 1 e 2 A, D 2 = ∇ − iA with e i and m i (i = 1, 2) the charge and mass of the scalar field ψ i (i = 1, 2), respectively.
In this paper we shall work with the probe limit in which the back reaction of matter fields to the metric can be ignored by taking e 2 → ∞ but keeping the ratio e = e 1 /e 2 finite. This allows us to start with the planar black hole as follows where f (z) = 1 − ( z z h ) 3 with z = 0 the AdS conformal boundary and z = z h the position of horizon. The temperature of the dual boundary system is given by the Hawking temperature Then on top of this background geometry the probe matter fields are governed by the following equations of motion, i.e., For simplicity but without loss of generality, below we will focus exclusively on the case of m 2 1 L 2 = 0, m 2 2 L 2 = −2. Moreover, to make our life easier, we shall turn off A x , A y , and require all the other quantities depend only on coordinates t and z. With this, the asymptotic behavior of the solutions near the AdS boundary can be written as where we have been working in the axial gauge A z = 0. According to the usual dictionary in the AdS/CFT correspondence, we know that the constants µ and ρ can be regarded as the chemical potential and charge density of the dual boundary theory, respectively. In addition, when the sources ψ −,1 and ψ −,2 are switched off, ψ +,1 and ψ +,2 correspond to the spontaneous broken order parameters of our two band holographic superconductor with dimension three and dimension two, respectively.
In what follows, we will work in the units where L = 1 and z h = 1, making use of the scaling symmetry of AdS. Then the evolution equations can be written as with the constraint equation

Phase Diagram and Thermodynamical Stability
In this section, we shall reproduce the phase diagram for the homogeneous and isotropic two band holographic superconductor. The corresponding bulk static equations of motion are There is a trivial solution ψ 1 = χ = 0 and A t = ρ(1 − z), which corresponds to the normal phase of our boundary system by holography. But generically it is hard to find out the analytic expression for the static bulk solution dual to the boundary superconducting phase. So we would like to obtain the corresponding numerical solution by pseudo-spectral method. To achieve this, we first set ψ 1 = |ψ 1 |e ieθ 1 = φ 1 e ieθ 1 and χ = |χ|e iθ 2 = φ 2 e iθ 2 , and then rewrite the above static equations to get the following five independent ones (see Appendix for details): Therefore we shall solve the above equations of motion to obtain our numerical solution. To this end, we are required to specify the appropriate boundary conditions. θ 1 = θ 2 = 0 can be achieved on the horizon by gauge fixing. In addition, the above equations of motion give rise to the following natural boundary conditions on the horizon. Finally since we want to find out the spontaneously broken phase for our boundary system at a fixed charge density, we impose the boundary conditions on the AdS boundary as follows: When there is only ψ 1 under the condensation, we call the phase as Phase-I. Phase-II represents the phase in which only ψ 2 is under the condensation. By our numerics we obtain the phase diagram in Figure 1 for the one band holographic superconductor, and the phase diagram in Figure 2 for our two band holographic superconductor. Here the condensations of the two scalar fields can coexist in a certain range of ρ, which we call as Phase-III.  In Figure 1, we denote the critical charge density for the condensations of scalar fields ψ 1 and ψ 2 as ρ 1 and ρ 2 , respectively. In Figure 2, the occurrence of the attenuation of the condensation of ψ 1 is also the beginning of the emergence of condensation of ψ 2 , and the corresponding critical charge density is denoted as ρ c1 . We use ρ c2 to denote the critical charge density where the condensation of ψ 1 vanishes.
To make sure Figure 2 represents the genuine phase diagram for our two band holographic superconductor, we are required to check whether the corresponding free energy density is the lowest among all possible phases 1 . Upon taking into account the back reaction, [14] shows that the grand potential results in the case of weak back reaction are consistent with the claim in [10], but direct confirmation in the probe limit is still lacking. Now we will calculate the free energy in the real probe limit. By holography, the free energy density can be obtained from the on shell Lagrangian of matter fields as follows: In the first line we have added the first boundary term to guarantee that we are working in the canonical ensemble rather than the grand canonical ensemble, and the last two boundary terms as the counter terms in the holographic renormalization to make the bulk on shell action finite, where h is the determinant of the induced metric on the AdS boundary. In the second line, we have made use of the equations of motion and the asymptotic AdS boundary behavior for the static solutions of matter fields. As shown in the left panel of Figure 3, compared to the normal phase the free energy density of three superconducting phases is much lower. Furthermore, as demonstrated in Figure 4 and the right panel of Figure 3, by scrutinizing the free energy density difference among these three superconducting phases, we lead to the above desired phase diagram in the canonical ensemble.  Figure 3. The left panel shows the free energy density difference of the three superconducting phases from the normal phase. The right panel shows the free energy density difference of Phase-I from Phase-II. Whence Phase-I is thermodynamically stable when the charge density is between ρ 1 and ρ c1 , while Phase-II is thermodynamically stable when the charge density is larger than ρ c2 .  Figure 4. The two panels show the free energy density difference of Phase-III from Phase-I and Phase-II between ρ c1 and ρ c2 , respectively. Whence Phase-III is the thermodynamically stable phase in this region.

Quasi-normal Modes
Before we move onto the fully non-linear dynamical stability by the real time evolution of our bulk fields, we would like to pause to play a little bit with the quasi-normal modes for our holographic model as a warm-up because the spectrum of quasi-normal modes can be regarded as a diagnosis of dynamical stability in the linear regime. If the phase in consideration is dynamically stable, then all the quasi-normal modes sit in the lower complex frequency plane. Actually as we shall show later on, all the thermodynamically stable phases in consideration is dynamically stable at such a linear level. In addition, he onset of Goldstone mode at the origin generically signals a transition from one phase to another. Moreover, among others, the lowest lying QNM is believed to capture the late time behavior during the real time fully non-linear evolution towards the desired final state as the late time perturbation is supposed to be small enough to validate the linear perturbation theory.
To our knowledge, so far there have been two methods developed to calculate the quasinormal modes for the gauged systems. One is the gauge invariant formalism [15], and the other is the gauge dependent formalism, where the gauge degree of freedom is freed [16]. We would like to take a third route by the gauge dependent formalism but with the fixed gauge. As seen later on, this alternative method is highly efficient and particularly suited for the comparison with the late time behavior of fully non-linear dynamical evolution of bulk fields.
To begin with, we first write down the linear perturbation equations of (2.7), (2.8) and (2.10) as Here we simply ignore the linear perturbation equation of (2.9) because it will be satisfied automatically in the whole bulk once it holds at the AdS boundary by our later boundary condition there. To calculate the quasi-normal modes, we further make the following consistent ansatz for the perturbed fields, i.e., with p,p, q,q and a complex functions of z. It is noteworthy thatp is independent of p andq is independent of q. Then substituting (4.4) into the above perturbation equations (4.1)-(4.3), we eventually end up with the following equations with p =p * = q =q * = a = 0 (4.10) on the AdS boundary. Note that for a small parameter λ and arbitrary w the gauge transformation always induces a spurious perturbation solution as follows a = −iwλ, p = eλψ 1 ,p = −eλ * ψ 1 , q = λχ,q = −λ * χ. This spurious perturbation solution can be removed by setting λ = 0, which can be further implemented by requiring a = 0 on the horizon 2 . The boundary conditions (4.10) together with a = 0 are naturally proposed with clear physical meaning. Whereas in the gauge invariant formalism, it is difficult to propose proper boundary conditions that exactly correspond to the gauge free boundary conditions that we need in the equations of dynamical evolution in the coexisting phase. This is one of the advantages of our method.
Finally, we cast the above linear perturbation equations and boundary conditions into the form L(w)v = 0 with v the perturbation fields evaluated at the grid points by the pseudo-spectral method. The quasi-normal frequencies are then obtained by the condition det[L(w)] = 0, which can be further identified by the density plot of |det[L(ω)] /det[L(ω)]| with the prime the derivative with respect to ω. The relevant results are summarized as follows. When ρ < ρ 1 , the system is in the normal phase and we thus have two sets of decoupled quasi-normal modes from δψ 1 and δχ, respectively. Both of the lowest lying modes are sitting symmetrically with respect to the imaginary axis. They migrate towards the origin with the increase of the charge density and the lowest lying modes from ψ 1 meet at the origin when ρ 1 is reached, which indicates the occurrence of phase transition to Phase-I.
When ρ 1 < ρ < ρ c1 , the system is in Phase-I and we still have two sets of decoupled quasi-normal modes, where one is from δψ 1 coupled to δA t , and the other is from δχ. The lowest lying modes from δχ keep migrating with the increase of the charge density and meet at the origin when ρ c1 is reached, signaling the phase transition to Phase-III. Among the quasi-normal modes from δψ 1 and δA t , one mode is pinned at the origin as Goldstone mode, and the other travels down the imaginary axis as Higgs mode with the increase of the charge density, eventually exceeded by the aforementioned lowest lying mode from δχ at a certain charge density ρ d as plotted in Figure 5.
When ρ c1 < ρ < ρ c2 , the system is in the coexistent phase and we only have one set of quasi-normal models from the coupled δψ 1 , δA t , and δχ. The lowest lying mode sits right at the imaginary axis. With the increase of the charge density, this mode travels down the imaginary axis from the origin till a certain charge density ρ m and then travels up to the origin as plotted in Figure 5.
When ρ > ρ c2 , the system is in Phase-II and we again have two sets of decoupled quasinormal modes, where one is from δψ 1 and the other is from δχ coupled to δA t . The lowest lying modes from δψ 1 migrate symmetrically away from the origin with the increase of the charge density. While the Higgs mode from δχ and δA t keeps traveling down the imaginary axis, eventually exceeded by the climbing-up modes, which sit symmetrically with respect to the imaginary axis.
In addition, we ascertain the points of phase transition by seeing that the lowest lying QNM crosses the real axis, and confirm the results that have been accurately obtained from the nonlinear static solutions. In some sense, the QNM provides a quicker and more reliable method to fix the points of phase transition, since it only involves solving the linear equations of perturbation, while the Newton-Raphson method used to solve the nonlinear static equations is inevitably sensitive to initial approximations.
As alluded to in the very beginning, suppose that the fully non-linear dynamical evolution ends up with one of the above phases, then the late time behavior is believed to be captured by the lowest lying quasi-normal modes on top of the final equilibrium state. To be more specific, if the final equilibrium state has a non-vanishing condensate | O f | = 0, then the late time behavior is supposed to be approximated by (4.14) On the other hand, if the final condensate | O f | = 0, then the late time behavior can be captured by In the next section, we shall check whether these thermodynamically stable phases are also dynamically stable at the fully non-linear level by the real time evolution and see whether the late time behavior can be captured by our lowest lying quasi-normal modes(QNMs).

Initial data and evolution scheme
Although the parameter space of ρ is divided into five parts by points ρ 1 , ρ 2 , ρ c1 , and ρ c2 , in what follows we shall focus only onto the regime where there can exist at least two superconducting phases, namely ρ 2 ≤ ρ ≤ ρ c1 , ρ c1 ≤ ρ ≤ ρ c2 and ρ ≥ ρ c2 .
Because Phase-I and Phase-II are always the possible phases in the above three regions, we will take two different kinds of initial data in each region. The first case is that we only let ψ 1 under the superconducting phase while give χ a perturbation of Gaussian wave packets, and we call this Case-I. The second case is that we only let χ under condensation while give ψ 1 a perturbation of Gaussian wave packets, and we call this Case-II.
In each case, the perturbation of Gaussian wave packets takes the form as s * e −30(z−0.5) 2 with s = 1, 0.1 or 0.01, which satisfies the source free boundary condition on the AdS boundary within our numeric error. In order to decrease the error during the long-time evolution, we do not evolve the dynamical equation (2.9) of A t , but obtain A t by solving the constraint equation (2.10), combining the boundary condition A t = −ρ on the AdS boundary and A t = 0 on the horizon.
With the above fixed charge density, we resort to the pseudo-spectral method to discretize the spatial direction, while in the time direction we take the forth order Runge-Kutta method to perform the evolution.The relevant numerical results will be presented in the following sections.

ρ 2 ≤ ρ ≤ ρ c1
Phase-I is thermodynamically stable in this region of ρ and we take ρ = 7.5 as an example. As shown in Figure 6, the larger the perturbation strength is, the more time it takes for the system to approach Phase-I in Case-I, while for Case-II, the system exhibits an inverse behavior as the perturbation strengthens, which is demonstrated in Figure 7.

ρ c1 ≤ ρ ≤ ρ c2
Phase-III is thermodynamically stable in this region of ρ and we take ρ = 8.5 as an example. As shown in Figure 8 for Case-I and Figure 9 for Case-II, the larger the perturbation strength is, the less time it takes for the system to approach Phase-III for both cases.

ρ ≥ ρ c2
Phase-II is thermodynamically stable in this region of ρ and we take ρ = 9.5 as an example. As shown in Figure 10 for Case-I, the larger the perturbation strength, the less time it takes for the system to approach Phase-II . For Case-II, the system exhibits an inverse behavior, which is demonstrated in Figure 11.

Late time behavior towards the equilibrium state
We see that the dynamical stability of the system is consistent with the thermodynamical stability in all the regions we are concerned with. But now we would like to turn to the detailed analysis of late time behavior of such a relaxation. One of the remarkable characteristics of the evolution is that the relaxation time of the system is generically much longer than that of the simplest holographic system with only one scalar field. 3 The relaxation time in such simplest system at 6.5 < ρ < 7.5 is less than one hundred time units. However, it is several thousand time units here at 6.5 < ρ < 7.5 (in the non-coexisting phase), so the difference is two orders of magnitude. This characteristic in return indicates us that the lowest lying QNM must be very close to the real axis (see Figure 12 for example), which would otherwise be confused with the mode w = 0. It is interesting to investigate the possible physical mechanism behind such long relaxation time in this model. Without loss of generality, we focus only onto the case of s = 1. As fitted in Figure 13 for χ at ρ = 7.5, there is an exponential decay with the decay frequency consistent with the lowest lying QNMs for χ, which is obtained by the aforementioned density plot in the right panel of Figure 14. so the decay frequency is also −0.0010934 √ ρ. This decay frequency is consistent with the lowest lying QNMs in the right panel of Figure 14. We find the late time behavior of ψ 1 is also controlled by the lowest lying QNMs at ρ = 9.5, where the decay frequency is obtained as −0.000718787 √ ρ. The late time behavior towards the coexisting phase is controlled by an exponential decay as well, with decay frequency −0.00025727 √ ρ for ρ = 8.5, which is consistent with its lowest lying QNM.
As we see, so far the late time behavior can be well captured by the lowest lying QNMs from the linear perturbation theory. But we find the non-linear effect starts to play a role in the late time behavior for ψ 1 in Phase-I and χ in Phase-II.
To be concrete, we denote the ratio of the absolute value of the imaginary part of the lowest lying QNM from δψ 1 coupled to δA t over that from δχ as r 1 for Phase-I, and the ratio of the absolute value of the imaginary part of the lowest lying QNM from δχ coupled to δA t over that from δψ 1 as r 2 for Phase-II. In the case of ρ = 7.5, we can see from Figure 14 that the lowest lying QNM from δψ 1 coupled to δA t sits at the imaginary axis, and r 1 is around 200. But as fitted in Figure 15, the real decay frequency for ψ 1 turns out to be the double of that for χ rather than its own lowest lying QNM frequency. In the case of ρ = 9.5, the decay frequency for χ is −0.00144 √ ρ which is also the double of that for ψ 1 mentioned before although the lowest lying QNM frequency for χ has the ratio r 2 = 500 or so. Such a double relation is believed to come from the very non-linear character of our full equations of motion. To be more specific, due to the above large ratio r 1 or r 2 , the lowest lying QNM for the non-vanishing order has died away long before the late time evolution in consideration. Then it follows from our non-linear equations of motion that such a double relation is generated by taking the lowest lying QNMs for the dying-out order as a source. This can also be regarded as the competition between the two orders, namely the dying-out order always tends to keep the other from relaxing towards the equilibrium state. Actually further numerical investigation indicates that this double relation occurs when r 1 or r 2 is bigger than 6. Beyond this regime, such a double relation becomes obscure and the situation becomes more complicated due to the non-linear effect.

Conclusion
In this paper, we investigate the dynamical stability of the two band holographic superconductor (the BHMRS model). In passing, we find it more convenient (if not inevitable) for one to calculate the QNMs by gauge fixing the gauge dependent formalism. Actually, this method is especially adapted for models with more than one scalar fields and for direct comparison with the late time behavior of the dynamical evolution as we performed later in this paper. As a result, the thermodynamically stable phase is also dynamically stable at the linear level. Furthermore, by evolving the system with different Gaussian wave packets as perturbations at the initial time, we find the dynamical stability of the system is also consistent with the thermodynamical stability at the fully non-linear level. We also analyze the late time behavior towards the final equilibrium state, and find it consistent with the lowest lying QNM in the coexisting phase as well as the lowest lying QNM for the field whose order dies out in the single ordered phases, while the late time behavior for the field whose order parameter does not vanish in the single ordered phases obeys a double relation or behaves in a much complicated way because of the non-linear effect induced competition between the two orders. In this sense, not only does the competition exhibit its existence in the phase diagram, but also does leave a footprint in the dynamical evolution. In addition, to our knowledge, this footprint is supposed to be the first counter example to the general belief that the late time behavior towards a final stable state can be captured by the lowest lying QNM obtained from the linear perturbation theory.
In comparison to the simplest holographic model with only one scalar field, the long relaxation time of the BHMRS model far from criticality is unexpected and interesting. This phenomenon may need further investigation, in particular for various other models with more than one order parameter.
Multiplying the equation (6.1) by φ 1 , we obtain (6.4) which gives rise to Note that we have f = 0 at the horizon z = 1 as well as A t = 0 there by (6.3), so we obtain C = 0, which further implies A t + f θ 1 = 0. (6.6)