A Non-Equilibrium Approach To Holographic Superconductors Using Gradient Flow

We study a charged scalar field in a bulk 3+1 dimensional anti-deSitter spacetime with a planar black hole background metric. Through the AdS/CFT correspondence this is equivalent to a strongly coupled field theory in 2+1 dimensions describing a superconductor. We use the gradient flow method and solve the flow equations numerically between two fixed points: a vacuum solution and a hairy black hole solution. We study the corresponding flow on the boundary between a normal metal phase and a superconducting phase. We show how the gradient flow moves fields between two fixed points in a way that minimizes the free energy of the system. At the fixed points of the flow the AdS/CFT correspondence provides an equivalence between the Euclidean on-shell action in the bulk and the free energy of the boundary, but it does not tell us about fields away from equilibrium. However, we can formally link static off-shell configurations in the bulk and in the boundary at the same point along the flow. For quasi-static evolution at least, it may be reasonable to think of this link as an extension of the AdS/CFT correspondance.

References 35 The AdS/CFT correspondence provides a possible mechanism for studying strongly interacting quantum field theories in d dimensions via the analysis of dual, weakly coupled gravitational systems in d + 1 dimensions. First proposed by Maldecena [1] in the context of string theory, in recent years the AdS/CFT correspondence has proven useful for studying a variety of condensed matter systems (see [2, 3] for excellent reviews). One of the greatest successes of this holographic approach has been its application to the study of quark gluon plasma [4].
In this paper we are specifically interested in a gravitational theory in 3 + 1 dimensions that gives a holographically dual description of a 2 + 1-dimensional superconductor. The model for a holographic superconductor that we consider was originally developed in [5,6].
The bulk theory consists of an electromagnetic field and charged scalar minimally coupled to Einstein gravity with a negative cosmological constant. Since the gravity theory is weakly coupled the superconducting theory on the boundary is strongly coupled. Useful reviews are found in Refs. [7][8][9].
Here we apply a gradient flow method to study the static solutions and nonequilibrium behavior of the holographic superconductor. Gradient flow is a general analogue of the heat equation that describes the path of steepest descent with respect to the free energy or Euclidean action of a system. It has been used in a variety of mathematical and physical contexts including Perelman's proof [10] of Thurston's geometrization conjecture, string theory [11], superconductors [12] and even image processing [13], [14].
First, we use gradient flow to reproduce a number of known solutions for a charged scalar field in a fixed anti-de Sitter (AdS) black hole background, and to obtain some solutions that are new to the best of our knowledge. Through the AdS/CFT correspondence, the solutions we obtain for the bulk theory allow us to study the properties of the strongly coupled superconductor on the boundary. We determine the critical temperature and magnetic fields below which the material becomes superconducting, and the coherence lengths for the condensate operator and the charge density. In addition to being a numerical tool for finding solutions, the gradient flow itself provides a good test of the nonlocal, as opposed to local, stability of these solutions. We verify this by calculating the (off-shell) Euclidean action (S E ) in the bulk along the flow and show that, as expected, the solutions minimize the free energy.
Secondly, we use the gradient flow to study the approach to thermal equilibrium of a holographic superconductor. Although the gradient flow is not equivalent to, or a replacement for, full time-dependent solutions, it may provide predictions for the approach to thermal equilibrium for holographic superconductors in cases where evolution is quasistatic and dissipation requires the system to be nonisolated from the environment so that energy is not conserved. Our flow calculations are carried out in the bulk, but they create a corresponding flow in the boundary operators. By solving for the flow parameter in terms of the energy, we obtain the path that describes the progression of the boundary operators as a function of the decreasing energy.
The paper is organized as follows: In Sec. II we briefly introduce the holographic superconducting model that we consider and provide a brief general review of gradient flow.
Section III derives the flow equations that we solve, while Sec. IV presents the numerical methods and results. We close in Sec. V with conclusions and prospects for future work.

II. PRELIMINARIES
A. The holographic superconductor Following [5], we start with a planar Schwarzschild anti-de Sitter black hole where L is the AdS radius, M is related to the Hawking temperature of the black hole T = 3M 1/3 /(4πL 4/3 ) and dΩ 2 R 2 is the metric on flat two-dimensional space. The black hole is 3+1 dimensional and is dual to a 2+1-dimensional theory on the boundary. We consider dΩ 2 R 2 in both polar and Cartesian coordinates. For the gravitational part of the action we consider the standard Einstein-Hilbert action with a negative cosmological constant, so that (1) is the unique (uncharged) vacuum solution. We want the boundary theory to describe a superconductor. This requires the boundary theory to contain charged fields whose condensation can lead to superconductivity. A conventional s-wave superconductor has an isotropic condensate described by a charged scalar field. The tensor properties of the corresponding bulk operator must be the same, and therefore we introduce a charged scalar in the bulk theory. We take this field to be minimally coupled to a bulk electromagnetic field. Our bulk action is therefore where ψ is a complex scalar field, and Λ is the cosmological constant.
the black hole horizon, in apparent contradiction to black hole no hair theorems. These theorems are based on the idea that matter outside a black hole wants to either fall into the black hole, or radiate out to infinity in the asymptotically flat case. The proof is formulated in terms of a black hole uniqueness theorem, which says that when gravity is coupled to matter fields, a stationary black hole solution is uniquely characterized by its conserved charges. However, there is no completely general no hair theorem, and counterexamples have been known since the 1990s [15].
In our problem, the formation of scalar hair is possible because we work in AdS space, where the negative cosmological constant acts like a confining box that prevents the charged particles from escaping to infinity. It is easy to see why the vacuum in the theory defined in (2) might be unstable to the formation of scalar hair. The effective mass of the scalar field is m 2 eff = m 2 + g tt q 2 A 2 t + · · · . Since g tt = −f (r) < 0 it is possible that the effective mass becomes sufficiently negative near the horizon to destabilize the scalar field.
We can also see that the formation of the instability could be temperature dependent.
Rewriting f (r) from Eq. (1) in terms of the temperature and horizon radius r 0 , defined from the equation f (r 0 ) = 0, we obtain where r 0 = (M L 2 ) 1/3 = 4πL 2 3 T is the radius of the event horizon and α := 4πT /3 is the temperature. From this expression we see that as the temperature decreases, |g tt | decreases at fixed r and therefore |g tt | increases, which means that the potential instability becomes stronger at low temperature. Following Ref. [16] we define an inverse radial coordinate u = r 0 /r. Using this notation the metric takes the form In these coordinates the black hole horizon is located at u = 1 and the boundary of AdS space is u = 0.
In this paper we look at static field configurations in the bulk theory. One of the variables we use to study their stability is the free energy. We discuss the definition of the free energy below. We start with the definition of the Euclidean action which is obtained from a Wick rotation of the Lorentzian action where we have separated the contributions due to gravity and matter by defining After the Wick rotation we have a new periodic Euclidean time coordinate with period, β = 1/T , inversely proportional to the temperature of the black hole [17]. For static, onshell field configurations, the Euclidean action is proportional to the Gibbs free energy of the system [18]. In the absence of time dependence, we can integrate over the time coordinate and write where the Gibbs free energy, G, involves an integral over the spatial coordinates only. In general, as stated in [7], it is necessary to add extra terms to the action that depend on the boundary conditions imposed on the fields. In this paper we consider a simplified energy functional for the matter fields only which, for the boundary conditions we study in this paper, does not require any additional terms. We have verified this explicitly by checking that our numerical results are independent of the grid spacing and size that is used to do the calculations (see Appendix A).

B. Gradient flow
Gradient flow is a general analogue of the heat equation that describes the evolution (or flow) of fields with respect to a flow parameter, τ [19]. We use a generalized index I to represent a given field, and all internal indices associated with that field. The variable x is taken to represent all spacetime coordinates. For example, using this notation, the QCD gauge field A a µ (t, x 1 · · · x d ) would be denoted Φ I (x), so that the Lorentz index µ, and color index a, and the fact that we refer to a gauge potential (denoted with the letter A) are represented by one index I. A quark field ψ(t, x 1 · · · x d ) would be represented similarly, but with a distinct index, for example as Φ J (x). The idea behind the gradient flow approach is to start with a set of fields Φ I (x) and corresponding generating functional F[Φ], and obtain a set of flow equations whose solutions move arbitrary initial field configurations along lines of steepest decent to an extremum of the generating functional.
Given a functional F[Φ] that depends on a set of fields Φ I , the flow equations are defined where we have assumed the existence of a configuration space metric G IJ , or equivalently, that the field space has a local, invertible and diffeomorphism invariant inner product that preserves any other symmetries of the system, We define an energy functional E in terms of the action by integrating out the time coordinate, similar to how the Gibbs free energy was defined in (8), so that the extrema of E correspond to static solutions of the Euclidean action. We note that for on-shell configurations the energy functional is the Gibbs free energy E = G. We use the energy functional E as the generating functional for our gradient flow equations. Equation (9) shows that, by construction, field configurations that extremize E will be fixed points of the flow and static equilibrium states of the superconductor. These fixed points can be stable, unstable or saddle points. Since the gradient flow finds the paths of steepest descent it is a natural tool for determining the stability of fixed points and studying the potential energy "landscape." The flow parameter is arbitrary up to constant rescalings. Here we are interested in fixed points of the flow and the path that the system takes as it approaches equilibrium neither of which are affected by such rescalings.
The metric G IJ must respect the symmetries of the theory, and is normally read off from the gradient term in the free energy. Specifically, consider a free energy of the form where H IJ depends only on the fields and not their derivatives. In this case it is natural to For our matter fields we take G IJ to be which is the simplest form for G IJ that still ensures the general covariance of the flow equations. In this paper we only consider the flow of the matter fields, so we do not need an inner product for the spacetime metric g µν .
When the system contains different species of fields, as in our case, one can consider the possibility that the different types of fields diffuse at different rates. For example, one can recover the time-dependent Schrodinger equations considered in the standard Ginsburg-Landau model [12] by considering the following configuration space metric with the following block structure: which reflects all the symmetries of the model. In the present case, where we consider bulk equations that come from an underlying supergravity theory, we assume for simplicity that the diffusion rates of the fields are the same, and consider the configuration space metric to come directly from the action.
From the differential flow equations (9) we obtain field configurations for finite values of the flow parameter Φ I (τ, x) starting from specified initial conditions Φ I (τ = 0, x). We use the word solutions to refer to the single parameter family of field configurations Φ I (τ, x).
The fixed points are defined as the configurations that satisfy dΦ I (τ, x)/dτ = 0. In our calculations, the fixed points that minimize the energy are found numerically as the configurations reached at the end of the flow Φ I (τ → ∞, x). The fixed points of the flow therefore satisfy which are the Euler-Lagrange equations of motion.
We note that the flow equations (9) are covariant with respect to field reparametrizations of the form Φ I →Φ I (Φ). However, even if the underlying theory is gauge and/or diffeomorphism invariant, the flow equations (9) are not necessarily gauge or diffeomorphism invariant. Consequently, the gauge choice made for the initial field configuration is not in general preserved by the flow. This problem is addressed by adding a deTurck term to the left-hand side of (9) that compensates the noninvariance of the right side. Suppose the action (and therefore the energy) is invariant under the infinitesmal transformation, for a set of arbitrary parameters ξ α (x) so that The K's generally involve differential operators and functions of the fields. For example, in the case of Yang-Mills theory a gauge transformation of the vector potential is Using the notation in (16) we have so that the index I includes the indices a and µ, and the fact that the covariant derivative involves the gauge field A; the index α becomes c; and the generic parameter ξ is called χ.
The general form of the flow equations, including the required deTurck term, is The role of the deTurck term is to ensure that any change in gauge along the flow can be compensated by a corresponding change in the parameters ξ α .
To see how the deTurck term preserves gauge invariance along the flow we write the rate of change of the energy with respect to the flow parameter τ as and substituting (21) into (22) we obtain where we have used (17) in the last step. The role of the deTurck term is illustrated explicitly in the context of the holographic superconductor in the following section. We also note that Eq. (23) shows that if all the eigenvalues of the metric are positive, the energy is a monotonically decreasing function of the flow paramter τ . In our case, since the eigenvalues of our metric are not all positive due to the term proportional to g tt , the energy may not be monotonic. However, for the solutions we find by starting near the static vacuum fixed point, we find that the energy is monotonically decreasing.

A. Parametrization
In this section we apply the formalism developed in the previous section to the action (2). It is convenient to write the complex scalar field ψ in terms of two real scalar fields, p and ω, We also rescale our vector fieldÃ µ = qA µ , which gives Using this notation the action in Eqs. (5)- (7) becomes

B. Probe limit
In this paper we consider the probe limit. Mathematically we reach the probe limit by taking q → ∞ in Eq. (26), so that the gravity part of the action decouples from the matter part. Physically this means that we neglect backreaction of the gauge and scalar fields on the geometry itself. The basis of the approximation is the assumption that the terms containing the matter fields are negligible in the equations of motion for the metric components, which can then be solved to obtain the metric in (1). The equations of motion for the matter fields are then calculated after we substitute the metric in Eq. (1) into the action, which means physically that we study the dynamics of the matter fields within a background metric.
In this limit we define an energy functional that includes only the matter part of the which in the case of static fields only is where the i, j indices are summed over spatial dimensions only. Since gradient flow moves along lines of steepest descent of E, at any point along the flow, E is stationary in the directions orthogonal to it. In principle these directions can be integrated out to define an off-shell free energy and in a saddle point approximation this free energy is simply equal to the value of the action at that point along the flow [19]. We therefore refer to E as the energy of the system. We are primarily interested in the gradient flow between two fixed points. Since fixed points correspond to on-shell configurations, E is equal to the matter contribution to the Gibbs free energy at the start and end points of the flow. From this point forward we drop the tildes that were introduced in Eqs. (25).
Using (30) as the generating functional, the gradient flow equations for the matter fields are obtained from (21). Before the inclusion of deTurck terms, the gradient flow equations for the matter fields are The right sides of the flow equations are manifestly gauge invariant, but the left-hand side Writing δp = 0, δω = −χ and δA µ = −∂ µ χ Eq. (16) gives K p = 0 for the field p, and K ω = −1, K Aµ = −∂ µ and ξ = χ for the remaining two fields. We see that the deTurck The flow equations (21) become At this stage χ is an arbitrary function.
A choice of the field ω can be interpreted as a gauge choice. We can ensure that the flow in (35) preserves any initial value of ω by choosing Making this choice our new set of flow equations is Since the bulk metric is asymptotically AdS space the asymptotic behavior of the scalar field can be determined from the Klein-Gordon equation in AdS 3+1 . Generally there are two possible falloff rates ∆ ± , but they are not always both normalizable. It is possible to consider tachyonic scalar fields with m 2 L 2 < 0, and Breitenlohner and Freedman (BF) showed that AdS d+1 spacetime is stable if the scalar field mass satisfies m 2 L 2 > −d 2 /4 [20]. Note that this bound is equivalent to the requirement that ∆ ± is real. For masses near the BF bound, −d 2 /4 + 1 > m 2 L 2 > −d 2 /4 both terms in (41) are normalizable [21], but if both coefficients are nonzero the theory is unstable in the asymptotic region [22]. We consider only the case m 2 L 2 = −2, where the asymptotic behavior of p has a simple form We can rewrite this equation in terms of the radial coordinate by recalling that we have defined u = r 0 /r with r 0 = αL 2 . Defining O ± = α ∆ ± c ± , Eq. (43) can be written where O 1 and O 2 are the vacuum expectation values of operators in the boundary theory with dimension 1 and 2 respectively. In this case both terms are normalizable, but we confine our interest to O 2 by taking c 1 = 0.
We also introduce a finite charge density and chemical potential, which are obtained from the scalar potential A t in the boundary theory [7]. The motivation is that an additional scale is necessary to produce a superconducting instability at low temperatures. For u → 0 we write where µ andρ are, respectively, the chemical potential and charge density in the boundary theory (note that we useρ for the charge density because ρ is used as a radial coordinate when studying solutions with axial symmetry). The magnetic field in the boundary theory is obtained from and the current can be obtained from the linear term in the expansion of the gauge potential around u = 0 [see Eq. (66)].

D. Gauge choice and Ansätze
We consider three separate cases, with different symmetries on the boundary.

Spatially independent case
We use coordinates in which dΩ 2 R 2 = dx 2 + dy 2 . We simplify the equations by setting A u = A x = A y = ω = 0. We also take A t and p to be functions of u only. The complete set of conditions we impose is The energy obtained from (30) is The flow equations for the fields in (47) becomeȦ u =Ȧ x =Ȧ y =ω = 0, where the dots denote derivatives with respect to the flow parameter. Thus the conditions in (47) are preserved by the flow. We also note that these conditions give χ = 0, which means that the deTurck term does not contribute. From Eqs. (38) and (40) we obtain the flow equations for the nonzero fieldsȦ by the flow. Equations (50) and (51) give a closed set of equations for two fields that depend on two spatial dimensions, x and u, and the flow parameter.

Translational symmetry
We again consider coordinates in which dΩ 2 R 2 = dx 2 + dy 2 and simplify the equations by setting A u = A x = ω = 0. We now take A t , A y and p to be functions of x and u only, so our problem has a translational symmetry along the y axis. The complete set of conditions we impose is It is easy to see that with these conditions, the flow equations for the fields in (52) becomė A u =Ȧ x =ω = 0, which shows that the conditions in (52) are preserved by the flow.
We also note that these conditions again give χ = 0, so that the deTurck term does not contribute. The energy (30) in this case is From equations (38) and (40) we obtain the flow equations for the nonzero fieldṡ Since the metric depends only on the coordinate u, the right sides of these equations are independent of y and t, which shows that the conditions in (53) are also preserved by the flow. Equations (55) -(57) give a closed set of equations for three fields that depend on two spatial dimensions, x and u, and the flow parameter.
We note that using these coordinates the magnetic field in the boundary theory is obtained from Eq. (46) as

Axial symmetry
We also consider coordinates where dΩ 2 R 2 = dρ 2 +ρ 2 dθ 2 . One motivation is that rotational symmetry allows us to study completely localized solutions that could be created in a lab.
When working with axial symmetry one typically looks for solutions where the phase of the complex scalar field can be written ω = nθ, and n is interpreted as an integer winding number. The value of the winding number is an important property of vortex solutions and leads to flux quantization in superconductors. We once again take ω = 0, which can be thought of as before as a gauge choice, but since we impose axial symmetry on the remaining fields, it also restricts us to solutions with zero winding number. We further take A u = A ρ = 0, and assume that our remaining fields are functions of u and ρ only. The complete set of conditions we use is The energy (30) becomes We note that as before the flow equations giveȦ u =Ȧ ρ =ω = 0 so that the conditions (59) and (60) are preserved by the flow, and we also have again from (37) that χ = 0, which means that the deTurck term does not contribute. The remaining flow equations arė The azimuthal component of the vector potential of the boundary theory can be written where B is the magnetic field and J θ is an azimuthal current.

IV. SOLUTIONS TO THE GRADIENT FLOW
To solve our gradient flow equations we need to proceed numerically, starting from a specified initial configuration for the fields. Since we are particularly interested in studying the flow between two fixed points, we consider the vacuum (hairless black hole) configuration which has a simple analytic form for the translationally symmetric case, and for the axially symmetric case. It is easy to verify that these configurations are fixed points One feature of the gradient flow method is that if our perturbation satisfies δp(u = 0) = 0, thenȦ µ (u = 0) = 0. This is true for all of our Ansätze and can be seen directly from Eqs (50), (55), (56), (63), and (64). We have therefore that in the boundary theory, the chemical potential µ = A t (0) and magnetic field B = ∂ x A y (0, x) or B = 1 ρ ∂ ρ A θ (0, ρ), are specified by our initial configuration and constant along the flow.
A more general statement is that using our method the gauge field on the boundary A µ (0, x), and all its derivatives with respect to boundary coordinates, are fixed by our initial configuration and unchanged along the flow. This means that the boundary theory does not have a dynamical gauge field, which corresponds to a limit where the superconductor is equivalent to a superfluid. It is possible to make the gauge field dynamical by including an additional boundary term in the bulk action and considering a different type of boundary condition on A µ [23]. While a dynamical gauge field is important for many superconductor phenomena such as the Meissner effect, a fixed background is sufficient to study how gradient flow in the bulk creates a corresponding flow in the boundary, and the extension of the flow to more complicated systems is straightforward. We comment that in the boundary theory, while derivatives with respect to boundary coordinates ∂ x A µ (0, x) and ∂ y A µ (0, x) are fixed Using this normalization we find that when the vacuum solution is unstable the quantity ∆E moves from an initial value of 0 into negative values. This type of behavior is typical in systems that exhibit spontaneous symmetry breaking, where a false vacuum decays into a more stable (but less symmetric) configuration.
In numerical calculations we set L = 1, which is equivalent to using the AdS radius as a length scale in all dimensional quantities (including the flow parameter). We also set the chemical potential µ = 1. This is equivalent to defining new ρ (or x, y) coordinates and absorbing µ into the parameter α in the flow equations as follows: In the following we omit the tildes but understand that setting µ = 1 in the solution implies we are in fact referring to the above rescaled quantities in our results. We vary the parameter α = 4πT /3 (or equivalently the temperature) and the magnetic field B.

A. Spatially independent solutions
By varying α = 4πT /3 we can determine the critical value below which the vacuum becomes unstable. We start with an initial perturbation of the scalar field given by δp(u) = 10 −3 × u 2 e −(u−1) 2 .
(73) Figure 1 shows how this initial perturbation evolves for two different values of the temperature. In Fig. 1(a) the bulk scalar field moves away from the (false) vacuum, and in Fig. 1(b) we see that at higher temperature the perturbed system returns to the vacuum configuration.
To compare with the results of [5], we recall O ± = α ∆ ± c ± and plot O 2 versus T . The parameter c 2 is obtained from a second derivative of the scalar field [see Eq. (41)], which is calculated using a finite difference formula on the first three data points. The error bars are obtained by comparing with the value extracted from the first, third and fifth data points.
The result is shown in Fig. 2(a), and agrees well with Fig. 1(b) in Ref. [5]. In Fig. 2(b) we show the same data using different variables: using α = 4πT /3 = r 0 /L 2 we plot c 2 versus α 2 . From Fig. 2(b) we see that the critical temperature is approximately α 2 c ≈ 0.059.
We also want to study the evolution of physical quantities in the boundary theory as a function of the flow parameter. In Fig. 3 we show the evolution of the energy, charge   solitons [24]. We can consider only the portion of the soliton where p > 0 by looking at only x > 0 and enforcing the condition that p(u, 0) = 0. We need to start from a perturbation that satisfies this condition and therefore we choose δp(u, x) = 10 −2 × u 2 e −10(u−1) 2 tanh (5x).
This flow can be interpreted as either a full soliton for −∞ < x < ∞ where p is antisymmetric around x = 0, or as a solution for x > 0 where x = 0 is an interface with a fixed Figure 4: c 2 as a function of −∆E for spatially independent fields at α 2 = 0.03.
vacuum solution for x < 0. In the boundary theory we can interpret the second case as an interface between a superconductor (for x > 0) and a normal material (for x < 0). In the Ginzburg-Landau (GL) theory of superconductivity there is an exact solution for the order where φ ∞ is the value of the order parameter in the pure superconducting phase, and ξ is the coherence length.
Although we would like a solution for all x > 0, to perform the numerical calculation we need to introduce a cutoff. We would like our cutoff to be large enough that our fields are constant for x > x max . From Eq. (75) we find that for We use a cutoff x max = 30 and we find x max /ξ > 5.9 for T ≈ 0.93T c .
Our numerical results indicate that the boundary operator has a similar x dependence, and we can therefore fit to a function of the same form as Eq. (75) to determine the coherence length of the holographic superconductor. The fit for α 2 = 0.03 is shown in Fig. 5(a).
We can see in Fig. 5(b) that, as expected, the coherence length diverges proportional to 1/ 1 − α/α c = 1/ 1 − T /T c as we approach the critical temperature. We note that due to this divergence, the condition x max /ξ > 1 cannot be satisfied very close to the critical temperature and for this reason we have considered only T ≤ 0.93T c . In the GL theory the charge density is proportional to φ 2 , We can therefore find a characteristic length for the charge densityρ(x) by fitting to . Contrary to what we expect from GL theory, we find that the two length scales are different. The difference between ξ and ξ q increases as we move further from the critical temperature, as can be seen in Fig. 6. This result agrees with what was found in [24] (apart from a difference in how coherence length is defined in the GL compared to Gross-Pitaevskii equations).
We can calculate the energy density and the boundary operator for the soliton configuration as functions of x and τ . In Fig. 7 we show the energy density integrated over x as a function of the flow parameter. Since ∆E decreases monotonically with τ , we can look at the evolution of the boundary operator as a function of |∆E| instead of τ . We have shown that our x cutoff is larger than typical values of the coherence lengths [see Fig. 5(b)], which means that field configurations are approximately constant at large x [see for example Fig.   5(a)]. We therefore expect c 2 (x = x max ) ≈ c 2 (x → ∞), and that the soliton solution at large x should be close to the spatially independent solution we considered in Sec. IV A. In Fig.   8 we plot the evolution of the boundary operator for large x, c 2 (x = x max ) , as a function of the energy. In the soliton case we have ∆E Sol ≈ −0.334, and if we calculate the energy of  the spatially independent case on the same interval (0 ≤ x ≤ 30) we find ∆E Ind ≈ −0.395.
These energies depend on the cutoff (x max ), but the difference ∆E Sol − ∆E Ind will be independent of the cutoff and can be thought of as a measure of the effect of the soliton.

C. Droplet solutions
The final case we consider is configurations with nonzero magnetic field. The interesting feature of these solutions is that the magnetic field can limit the formation of the scalar hair/condensate. There is a critical magnetic field B c above which the boundary operator does not condense. For B < B c the condensate forms only in a localized region on the boundary, and these localized solutions are called droplet solutions [25]. For the droplet solutions, the initial perturbation does not need to depend on x, so we use the perturbation δp(u) = 10 −2 × u 2 e −10(u−1) 2 .
(77) Figure 9: c 2 (x) for several values of the magnetic field, and α 2 = 0.03, for a translationally symmetric droplet. As B increases c 2 (x) becomes smaller and more localized.

Translationally symmetric droplets
In Fig. 9 we show how the magnetic field alters the boundary operator c 2 (x). When B → 0 we recover the spatially independent solution [see Fig. 3(b)]. As B approaches the critical value from below, the droplets become narrower and shorter. In Fig. 10 we see how the energy evolves as the flow moves the fields from a small perturbation of the vacuum solution towards a localized droplet of scalar hair. In Fig. 11 we plot the energy versus the value of the operator c 2 at the center of the droplet. We notice in Fig. 10 that when ∆E is small there are noticeable fluctuations due to numerical error, which influence the low ∆E behavior in Fig. 11.

Axially symmetric droplet solutions
The droplet solutions do not require translational symmetry. In this section we start from the same initial perturbation but instead enforce axial symmetry along the flow. This leads to droplet solutions that are entirely localized (as condensate in the boundary theory and as scalar hair on the black hole horizon). Our results are similar to those found in [26], with the important difference that we have a constant magnetic field on the boundary. As discussed at the beginning of Sec. IV, this is a characteristic of the gradient flow method. We study the flow up to a maximum radius ρ = 30, for several different values of the magnetic field.
The final configurations for c 2 are shown in Fig. 12. We can examine several other quantities in the boundary theory, namely the charge densityρ and the azimuthal current J θ . Figure 15 shows the charge density profile of the droplet. Figure 16 shows the axial current. We note that the formation of currents at the edge of the superconducting droplets is expected as the superconductor will attempt to expel any magnetic fields.

V. CONCLUSION
We have demonstrated the versatility of a gradient flow approach to holographic superconductors. In addition to reliably finding solutions to the equations of motion, we are also able to determine the stability of the solutions in a nonperturbative way. We find stable axially symmetric droplets with a constant background magnetic field, whereas the droplets found in [26] would enhance or weaken the magnetic field at the core depending on the temperature. The gradient flow is a much more general approach to finding droplets as we can consider any magnetic field as a fixed background specified by initial conditions.
Although the flow parameter itself does not have a straightforward physical interpretation, we can exploit its connection with the energy to gain insight into how the system can undergo the phase change from hairless black hole (normal phase) to a black hole with scalar hair (superconducting phase) in a quasistatic way.
The AdS/CFT correspondence provides an equivalence between the Euclidean on-shell action in the bulk and the free energy of the boundary. This means that at the fixed points of the flow the free energy of the bulk, E, is equivalent to the free energy of the superconductor. Away from the fixed points the AdS/CFT dictionary does not tell us anything about the relationship between the two energies. However, using the gradient flow method, we can formally link static off-shell configurations in the bulk and in the boundary at the same value of the flow parameter τ . For quasistatic evolution at least, it may be reasonable to think of this link as an extension of the AdS/CFT correspondence.
The extension of the gradient flow method to include the metric flow would allow us to study the system away from the probe limit since the stationary points of the metric flow will be the fully backreacted metric. Although holographic superconductors with metric backreactions have been studied, it is unclear how allowing the metric to change will influence the stability of the scalar hair solutions. To study the system away from the probe limit it is necessary to include additional boundary terms in the action to regulate divergences of the free energy. The evolution of the metric from the vacuum black hole solution is a physically interesting problem in its own right, and this type of metric curvature flow with source terms is an interesting problem in mathematics.
Another application of gradient flow to AdS/CFT that is an active area of research is the Hamilton-Jacobi (HJ) first order equations [27,28]. where the system is never far from equilibrium, it may be possible to study the dynamics of the system using some approximation. The gradient flow method, which takes fields along the path of steepest descent towards the extremum of the action, is a good candidate.
ACKNOWLEDGMENTS This research was supported in part by the Natural Sciences and Engineering Research Council of Canada, and a University of Manitoba graduate fellowship.

Appendix A: Numerical Methods
Since our gradient flow equations are analogous to the heat equation we solve them using a simple explicit finite difference scheme with forward differences in flow time and centered differences in our spatial coordinates. Such methods tend to be stable and convergent as long as the time step is proportionally smaller than the square of the space step, For a simple one-dimensional heat equation it can be shown that C = 0.5, but it is difficult to precisely compute the value of C for the type of nonlinear coupled partial differential equations that we have solved. We typically work with an equal number grid points in both u and x (or ρ). Since u ∈ (0, 1) and x ∈ (0, 30) we have du < dx and therefore we do not need to consider a separate convergence condition of the form (A1) for dx.
We find that dτ obtained from (A1) with C = 0.25 is sufficient for good convergence.
Using this value of C we see that the computer time required to reach τ = 1 is proportional to du −3 for the one-dimensional case and du −4 for two dimensions. To increase computational speed, we therefore want du to be as large as possible without sacrificing accuracy.
Most of the results presented in this paper were calculated with a 300 × 300 grid. The exception is Figs. 3 and 4. where we use 2000 grid points. In this case we found that with fewer grid points there was small but noticeable nonmonotonic behavior in ∆E which suggested that the endpoint of the flow did not minimize the free energy of the system. Figs.
17 and 18 show this behavior with 300 grid points. The nonmonotonicity is more noticeable at lower temperatures since p(u) is larger in those cases, and the fact that it disappears when the grid spacing is reduced (see Figs. 3 and 4) proves that the effect is numerical.
Although we do not see the same behavior in the spatially dependent cases we instead see smaller fluctuations of the energy, in particular when it is small relative to the grid spacing.
These fluctuations are due to the larger error in the spatial derivatives since dx > du. In Fig. 19 we look at the dependence of c 2 , ∆E and p(u = 1) on the number of grid points used (equivalently du −1 ), at the end of the flow for the spatially independent case where the computation time depends least on grid spacing. We see that quantities approach a fixed value as the number of grid points increases. We note that the range for the y axis on these plots is roughly a %10 variation for the energy, and %2 for c 2 and p(u = 1). We expect the larger error in ∆E since it is obtained from a numerically computed integral of fields which introduces additional error. The fields themselves and c 2 depend mainly on the error in first order finite differences which is proportional to du. In Fig. 19(d) we see that p(u = 1) is linear in du, and we can extrapolate to du = 0 to find that p(u = 1) ≈ 2.0034 is within %1 of the value obtained using only 300 grid points.

Boundary conditions
Since we are using centered finite differences we need to take special care of the points at the boundaries. For all cases we treat the u = 0 and u = 1 boundaries in the same way.
At the AdS boundary, u = 0, all terms with derivatives are multiplied by a factor u 2 . We assume that the derivatives of our fields are finite at u = 0, which means that any term with a derivative does not contribute to the flow equations. At the horizon, u = 1, the factor h(u) goes to 0. We require that A t (u = 1) = 0 in such a way that the ratio A 2 t h(u) = 0 at u = 1. Any u derivatives that we need are calculated using a one-sided finite difference.
For the boundaries at x max and ρ max we can usually use one-sided finite differences since typically the fields are approximately constant at this boundary. The one exception is the field A θ for which we enforce the condition ∂ ρ A θ (ρ max ) = Bρ. For a finite superconductor with radius ρ max , this is simply the condition that there is a fixed external magnetic field.
The boundaries at x = 0 and ρ = 0 can be handled by adding an extra grid point at x = −dx and calculating centered finite differences as usual. The value of the fields at this point is determined using the symmetry of the configuration. We take A t to be always symmetric, A y and A ρ are always antisymmetric, and p is symmetric except for the soliton case, where it is antisymmetric.