Holographic Abrikosov Lattices

We study black hole solutions of $D=4$ Einstein-Maxwell theory coupled to a charged scalar field that are holographically dual to a $d=3$ conformal field theory with a non-vanishing chemical potential and constant magnetic field. We numerically construct black hole solutions that are dual to a superfluid phase with a periodic lattice of vortices. For the specific model we investigate, we find that the thermodynamically preferred configuration is given by a triangular lattice and moreover the vortices are associated with the lowest Landau level. We also construct black holes describing a lattice of vortices associated with the next to lowest Landau level and while theses are not thermodynamically preferred they exhibit some interesting features that could be realised for other holographic models.


Introduction
Holography provides a controlled theoretical framework to study strongly coupled quantum field theories. In seeking possible applications to real systems an important development, pioneered by Steven Gubser, was the realisation that holographic matter can exist in a superfluid phase [1]. In the simplest set-up one considers Einstein-Maxwell theory coupled to a charged scalar field with an AdS vacuum solution that is dual to a conformal field theory with a global abelian symmetry. When the CFT is held at finite chemical potential, the unbroken phase at high temperature is described by an electrically charged, planar AdS-RN black hole solution with vanishing charged scalar field. This black hole is unstable below some critical temperature and the system condenses into a superfluid phase which is described by an electrically charged black hole carrying a halo of charged scalar hair [1][2][3].
In this paper we study superfluid phases of holographic matter held at finite chemical potential, with the addition of an external magnetic field. Over the past ten years this topic has been studied from several different points of view. Switching on the magnetic field suppresses the superfluid phase transition, as one might expect. Indeed, and as we will review, the critical temperature at which the superfluid instability sets in decreases as one increases the magnetic field and for large enough magnetic field the instability is longer present [3][4][5]. Below the critical temperature one expects the existence of vortices. The defining feature of a vortex is that the phase of the complex field has non-zero winding as one goes around the vortex and, as a consequence, the complex field vanishes at the core of the vortex. Using certain probe approximations, where there is no back-reaction on the metric, constructions of vortex-like solutions were made in [6][7][8]. A further development, again in a probe approximation, was the construction of a vortex lattice, using vortices in the lowest Landau level [9]. Going beyond the probe approximation the existence of a vortex lattice solution was argued for in [10], again in the lowest Landau level, by considering a perturbative expansion about a purely magnetic AdS 2 × R 2 zero temperature ground state. An approximate back-reacted vortex lattice solution was recently constructed in [11] after imposing by hand a certain circular symmetry on each unit cell.
The purpose of this paper is to report on the first numerical construction of fully back-reacted black hole solutions, without any approximations, that describe a periodic vortex lattice in a superfluid phase. A priori it is not clear what the thermodynamically preferred shape of the vortex lattice will be. If one is just below the critical temperature and one is able to utilise a Landau-Ginzburg description, one finds that a triangular lattice associated with the lowest Landau level is preferred (for a review see [12]). Thus, one might expect that in holography the triangular vortex lattice is also the preferred configuration at least just below the critical temperature, and in the context of the probe approximation some arguments supporting this conclusion were given in [9]. However, it is worth emphasising that in certain holographic situations the Landau-Ginzburg does not effectively capture the properties of the phase transitions near the critical temperature [13] and so this conclusion may not be valid in general. In any event, as the temperature is lowered the Landau-Ginzburg description becomes less useful and it is no longer clear what shape the preferred lattice will take. In fact various different shapes are realised in real superconductors as well as transitions to other phases such as vortex liquids and glasses [12]. It is therefore of significant interest to find out what can happen in the context of holography and this paper is a step in exploring what is possible.
The D = 4 gravitational model that we will consider couples the metric to a Maxwell field and a complex scalar field. The model has an AdS 4 vacuum solution, with vanishing Maxwell and scalar field, which is dual to the underlying d = 3 CFT that we want to study both with non-vanishing magnetic field and at finite chemical with respect to the global U (1) symmetry. It also has another AdS 4 solution with nonvanishing scalar field which describes the IR behaviour of the superfluid phase with vanishing magnetic field, B = 0. When B = 0 the high temperature, unbroken phase is described by the dyonic AdS-RN black hole solution. We review the linearised instabilities of this black hole and, show that the critical temperature at which it becomes unstable just depends on the Landau level of the linearised perturbation, with the lowest Landau level having the highest critical temperature. Within the linearised framework one can then construct vortex lattice solutions, parametrised by the Landau level as well as two additional parameters which determine the shape of the lattice.
For a specific value of the magnetic field, we construct the back-reacted vortex lattice associated with the lowest Landau level. By minimising the free energy of the black hole solutions with respect to the remaining two shape parameters, we show that the triangular lattice is the preferred configuration for the specific temperatures we consider. While a further refinement of our numerics is required in order to construct black holes at very low temperatures, at the end of the paper we discuss a plausible zero temperature ground state solution. We also construct vortex lattice solutions for the second lowest Landau level, which appear at a critical temperature that is lower than those associated with the lowest Landau level. Interestingly, the thermodynamically preferred black holes in this class are associated with infinitely thin and long lattice structures, indicating the existence of an interesting kind of linear vortex defect. For temperatures when both black holes exist we find that the triangular vortex lattice associated with the lowest Landau level is always thermodynamically preferred. It is worth noting, however, that this conclusion certainly depends on the bulk gravitational model, a point we return to in the discussion section.
The holographic black hole solutions that we construct consist of a lattice of vortex tubes that stretch out from the black hole horizon and extend to the asymptotic boundary. Since the proper radius of the vortex tubes grows as one approaches the boundary they have a funnel-type structure. We emphasise that in the boundary theory the Maxwell field is not dynamical and hence the magnetic field is not localised inside the vortices. Instead, the vortices are associated with circulating currents in the boundary theory. It should be noted that this set-up is different from the usual superfluid vortices which carry quantised orbital angular momentum. It would be interesting to know if there are experimental setups where the configurations we discuss can be realised.
The plan of the rest of this paper is as follows. In section 2 we introduce the bottom-up holographic model of interest. In section 3 we discuss various details of the gravitational boundary value problem relevant to the formation of an Abrikosov lattice of vortices. We then discuss the numerical techniques we use to construct the broken phase black hole solutions along with their thermodynamics. In section 4 we discuss the main results of the numerical analysis and we conclude with some discussion in section 5. Appendix A has some details about the asymptotic expansions and one point functions while appendix B contains some comments about our numerical scheme and convergence properties.

Set-up
We will consider a bulk action in D = 4 spacetime dimensions of the form with F = dA and D µ ψ = ∇ µ ψ + i q A µ ψ. We also take V = V (|ψ| 2 ) so the action is invariant under the gauge transformation ψ → e −iqΛ ψ, A → A+dΛ. For simplicity we have fixed Newton's constant so that 16πG = 1. The equations of motion associated with (2.1) are given by We focus on the specific choice of potential given by 3) The equations of motion then admit a unit-radius AdS 4 vacuum solution with A = ψ = 0, which is dual to a d = 3 CFT with an abelian global symmetry. We will choose boundary conditions so that the complex scalar field ψ, of charge q, is dual to an operator O ψ , with scaling dimension ∆ ψ = 2. We want to analyse this CFT at finite temperature T , with constant chemical potential µ, and constant magnetic field B. The high temperature, spatially homogeneous and isotropic phase is described by the planar, dyonic AdS-Reissner-Nordström (AdS-RN) black brane solution. For later convenience it will be useful to write this in the following non-standard form The AdS 4 boundary is located at r → 0 with the metric approaching 6) and note that to get the standard form we can scale r by r + /2. The black hole horizon is located at r = 1 in these coordinates and the temperature of the black hole is given by As T → 0 the black holes approach a dyonic AdS 2 × R 2 solution in the IR with finite entropy density. When B = 0, with ∆ ψ = 2 and any value of q [24], it is well known that below some critical temperature the AdS-RN solution is unstable and the thermodynamically preferred black hole, describing a superfluid phase, has non-vanishing charged scalar hair. Our choice of potential (2.3) is such that the T = 0 limit of these superfluid black hole solutions are domain walls interpolating between the AdS 4 vacuum in the UV and another AdS 4 solution in the IR with |ψ| = 1 and radius squared equal to 12/13. The entropy density of these black holes goes to zero as T → 0, with s ∝ T 2 .
When B = 0, the AdS-RN black hole continues to be unstable below some critical temperature that depends on |qB|/µ 2 , up to some critical value of the magnetic field |qB c |/µ 2 , as we will review shortly. Below this critical temperature it is known that a vortex lattice can form. We will numerically construct such vortex lattice black hole solutions in the sequel.

General considerations
We will construct static black hole solutions with planar horizons that asymptote to AdS 4 in the UV with quasi-periodic boundary conditions for the x and y directions. We will use coordinates (t, r, x, y) which are globally defined outside the black hole horizon. At fixed t, r the x, y coordinates parametrise a two-torus associated with a flat metric dx 2 + dy 2 , with the following identifications (x, y) ∼ (x + L x , y) , (x, y) ∼ x + v L y , y + L y . (3.1) The parameter v ∈ [0, ∞), which can be exchanged for an angle β via cos β = v/(1 + v 2 ) 1/2 , governs the deviation of the shape of the torus from a rectangular torus. Without loss of generality, we can split the gauge field as where B is the external magnetic field and a = a µ (r, x, y)dx µ is a one-form that is globally defined on the spatial torus. Similarly, the metric components g µν are also globally defined on the torus. We now turn our attention to the complex scalar field ψ. It is clear from the gauge field decomposition (3.2) that imposing ψ to be well defined on the torus would fail to satisfy the equations of motion (2.2). Therefore, we take ψ (r, x + L x , y) = e ig 1 (r,x,y) ψ (r, x, y) , ψ r, x + v L y , y + L y = e ig 2 (r,x,y) ψ (r, x, y) , with g 1 and g 2 real functions on the torus and demand that D µ ψ (r, x + L x , y) = e ig 1 (r,x,y) D µ ψ (r, x, y) , D µ ψ r, x + v L y , y + L y = e ig 2 (r,x,y) D µ ψ (r, x, y) .
The general solution for the compatibility of the conditions (3.4) and (3.5) is where c 1 and c 2 are real constants of integration. Thus, the boundary conditions for the complex scalar are given by Furthermore, compatibility between these two gives the quantisation condition where n is an integer. The condition (3.7) implies that the modulus of the complex field is periodic on the torus, while the phase of the complex field has winding number n as we go anti-clockwise around the unit cell. Thus, each unit cell contains n vortices. For a given such solution to the equations of motion, ψ(r, x, y), we can construct another quasi-periodic solution with complex scalar field,ψ(r, x, y), parametrised by two constants x 0 , y 0 , by using the following combination of a translation and a gauge transformationψ and noting that A(x, y) = A(x + x 0 , y + y 0 ) + d(By 0 x). The boundary conditions for ψ(x, y) are as in (3.7) but with c 1 → c 1 − qBL x y 0 and c 2 → c 2 − qBL y (vy 0 − x 0 ). The two constants, c 1 and c 2 can be interpreted as parametrising two Goldstone modes in the boundary theory. These Goldstone modes are associated with translations that are intertwined with the internal abelian global symmetry, due the presence of the magnetic field. An additional Goldstone mode, associated with the abelian global symmetry is obtained by multiplying the solution by a constant phase. All of these Goldstone modes need to be fixed in order to find a solution numerically.

Perturbative zero modes
We now consider zero modes around the dyonic Reissner-Nordström black hole solution given in (2.4). The linearised complex scalar field equation of motion in this background is where Separating variables by writing ψ = Φ(x, y)ρ(r), we deduce that Φ has to satisfy the eigenvalue equation 12) and the radial function will satisfy The task is to solve the eigenvalue equation (3.12) with the boundary conditions (3.7): along with the constraint (3.8). Focussing on the x argument, without loss of generality and consistent with the first boundary condition, we can write Substituting this into (3.12) we discover that W (λ) l (y) satisfy This is essentially the simple harmonic oscillator; the eigenvalues, λ j , are labelled by an integer j = 0, 1, 2, . . . , the Landau level, with and ψ j are the standard harmonic oscillator wave functions and d l are constants. Thus, for a given Landau level j, we have We also need to impose the second boundary condition in (3.14) which leads to a recursion relation on the d l which, using the quantisation condition (3.8), we can solve as where z is a real constant and c is a complex constant. The constant c is just an overall free constant associated with the fact that we are solving a linear equation and so we set c = 1. Now shifting the y coordinate via y → y + c 1 /(qBL x ) combined with the gauge transformation with parameter Λ = c 1 /(qL x )x, as in (3.9), allows us to set c 1 = 0. Having set c 1 = 0, it will be convenient to shift the x coordinate via x → x − L x z/(2πn), again as in (3.9), in order to set z = 0. Thus, in summary we can write the spatial part of our zero mode as which satisfies the boundary conditions given in (3.7) with c 1 = 0 and c 2 = 1 2 vL 2 y qB. We note that we have fixed all of the Goldstone modes. For a given magnetic field, and hence fixed qB, the free parameters specifying these zero modes are the torus parameters (L x , L y , v), constrained via the quantisation condition qBL x L y = 2πn, as well as the Landau level j. The eigenvalue λ in (3.12) only depends on j and thus, the radial equation (3.13), and hence the critical temperature T c at which the Reissner-Nordström solution becomes unstable, also only depend on j.
In figure 1 we have plotted T c as a function of |qB|. When B = 0 we see that the dyonic Reissner-Nordström black hole becomes unstable to forming a superconducting state at T c /µ ∼ 0.090. The instability persists in the range 0 < |qB|/µ 2 0.846. For small enough |qB|/µ 2 there is a finite sequence of instabilities appearing at lower temperatures that are associated with higher Landau levels. In the sequel for q = 2 and |qB|/µ 2 = 0.02 we will construct fully back-reacted black hole solutions for the lowest Landau level with j = 0, which first appear at T c /µ ∼ 0.088, and also the first Landau level with j = 1, which first appear at T c /µ ∼ 0.085. In each class of solutions the moduli space of solutions is parametrised by (L x , L y , v), constrained via the quantisation condition qBL x L y = 2πn. We will determine the thermodynamically preferred configurations both within each class of solutions and also show, for temperatures when they both exist, that the solutions in the lowest Landau level are preferred.  Figure 1: Instabilities of the dyonic Reissner-Nordström black hole. We have plotted the critical temperature, T c , against the magnetic field, B, using dimensionless units, for three different Landau levels: j = 0 (blue), j = 1 (red) and j = 2 (green). We construct fully back-reacted black hole solutions for |qB|/µ 2 = 0.02 with j = 0 and j = 1 that first appear for values of T c /µ ∼ 0.088 and T c /µ ∼ 0.085, respectively.

New coordinates on the torus
For numerical convenience, it is convenient to bring the identifications of the coordinates on the torus into the following form This is achieved by performing the following coordinate transformatioñ In these coordinates the parameters of the torus get encoded in the metric. Indeed, the full boundary metric and gauge field have the form The last term in the gauge field can be removed via a gauge transformation to give Furthermore, if the complex scalar field ψ satisfies the quasi-periodic boundary conditions (3.8) with c 1 = 0 and c 2 = 1 2 vqBL 2 y / 1 + v 2 , precisely as we chose for the zero modes given in (3.20), then we havẽ These are the boundary conditions we will use in the ansatz for the numerical integration. This fixes the two translational Goldstone modes mentioned below (3.9). The Goldstone mode associated with multiplication by an overall phase is fixed by a boundary condition for the complex scalar on the horizon as discussed in the next subsection. Somewhat for historical reasons, in the numerical integration we used the following rescaled quantities 26) in terms of which the torus part of the boundary metric in (3.23) takes the more symmetric form Note that since L x L y =L xLy the form of the quantisation condition is unchanged: For later convenience we also define a quantity k such that which respects the quantisation condition. Note that for the perturbative modes the solution (3.20) is independent of n. The only appearance of n is in the quantisation condition (3.28). Thus, one can obtain different values of n just by rescalingL x andL y . Going beyond the probe approximation, as we discuss in the next section, the solutions to the non-linear equations that we have found, all have n = 1. We do not know if it possible to construct other solutions with n = 1.

Numerical integration
Motivated by the form of the boundary metric (3.23), we consider the following ansatz for the vortex lattice black holes In these expressions we have used the one-forms η t , ηx, ηỹ defined by The functions f (r) and g(r) are precisely the same functions appearing in the Reissner-Nordström solution (2.5) and are incorporated for convenience. The remaining functions defined by F ≡ {Q tt , Q rr , Q, R, W, Q tx , Q tỹ , Q tr , Q rx , Q rỹ , a t , a r , ax, aỹ, φ 1 , φ 2 } are all functions of the radial coordinate r as well as (x,ỹ), where the latter parametrise the torus with the identifications given in (3.21). On this torus, the scalar fields φ 1 and φ 2 satisfy the quasi-periodic boundary conditions (3.25) while all the remaining functions in F are periodic. There is significant redundancy in this ansatz and this will be dealt with momentarily. We demand that the solutions approach an asymptotic AdS 4 boundary, located at r = 0. The boundary conditions that we want to impose are given by Q tt (0,x,ỹ) = Q rr (0,x,ỹ) = Q(0,x,ỹ) = W (0,x,ỹ) = 1 , Notice that these boundary conditions imply that the asymptotic metric approaches (3.27). Furthermore, the asymptotic form of the gauge field is A → µdt−BL xLyỹ dx, also as desired. Due to the presence of g(r) in the ansatz for ψ in (3.30), the boundary conditions on φ 1 and φ 2 are associated with a spontaneous breaking of the global U (1) symmetry and we recall that we have assumed ψ is dual to an operator with ∆ ψ = 2.
We will discuss the one point functions of this operator as well as the stress tensor and U (1) current later.
We also demand that we have a Killing horizon, generated by the Killing vector ∂ t , located at r = 1. This is achieved by demanding that the set of functions F(r,x,ỹ), appearing in (3.30), admit an expansion in powers of (1 − r) of the form The equations of motion impose constraints on the coefficients appearing in (3.33). In particular, for the components Q rr (r,x,ỹ) and Q tt (r,x,ỹ) we obtain the condition for constant surface gravity Q rr (1,x,ỹ) = Q tt (1,x,ỹ), while for all other components we impose the Neumann boundary condition ∂ r F| r=1 = 0. We also impose the following boundary condition on the complex scalar field at the horizon, φ 2 (1, 1/2, 1/2) = 0, which fixes the remaining Goldstone mode, as commented above. At this point, we need to address the fundamental issue that the PDEs obtained after substituting the ansatz (3.30) into (2.2) are weakly elliptic, meaning that they are elliptic only for the physical degrees of freedom, and thus they are unsuitable for numerical solution without gauge fixing. To deal with this issue we use the DeTurck method, following [25,26]. We first modify the Einstein equations given in (2.2), to obtain Einstein-DeTurck equations, by making the replacement Hereḡ denotes a reference metric andΓ is the Christoffel connection ofḡ. We need to choose the reference metric to have a Killing horizon and the same asymptotic behaviour as the solutions we would like to construct; we will take it to be given by the metric in the dyonic AdS-RN black hole solution (2.4). For this choice, we find that ξ µ = 0 on the boundary. After numerical integration we need to check, a posteriori, that we have ξ = 0 everywhere. The condition ξ = 0 corresponds to fixing a set of coordinates. Similarly, we also need to fix the gauge freedom for the gauge-field. As in [22] (see also [21]) we modify the Maxwell equation in (2.2) via whereĀ ν is a reference gauge field. WhileĀ ν is not needed in order to obtain an elliptic system of PDEs, it does allow us to suitably choose A µ −Ā µ for our boundary value problem. Indeed we takeĀ = −BL xLyỹ dx so that A−Ā are periodic functions and, as one can then explicitly check, ϕ then vanishes on the boundary. As in [22], by taking the divergence of the modified Maxwell equation and then integrating over the bulk, one can then deduce that ϕ = 0 everywhere. The vanishing of ϕ fixes the gauge invariance.
To summarise, to obtain our vortex lattice solutions, we will solve the modified equations of motion, given in (3.34),(3.35), which gives a system of elliptic PDEs with a well defined boundary value problem. Having solved them we then check that ξ µ = 0 and since ξ µ is spacelike, this is achieved by checking ξ 2 = 0. Some additional comments on the numerical approach we take, as well as a discussion of the numerical convergence is given in appendix B.
Returning to our ansatz, (3.30), we see that it is left invariant if we make the replacements:L x ↔L y , W → W −1 , Q rx ↔ Q rỹ , Q tx ↔ Q tỹ , ax → aỹ + BL yỹ , aỹ → ax − BL xỹ and interchangex ↔ỹ. This symmetry is then reflected in our solutions. In particular, after recalling (3.29), it implies that any physical quantities that are obtained by integrating overx,ỹ, such as the free energy, for example, will be invariant under the interchange of k → 1/k.

Thermodynamics and one point functions
The thermodynamic properties of the black holes are obtained by analytically continuing the time coordinate via t = −iτ . Demanding regularity of the solution at the black hole horizon we obtain the temperature, T , which has the same form as given in (2.7). We can also read off the area of the event horizon and, since we are working in units with 16πG = 1, we deduce that the entropy density is given by the horizon integral where the integral is over a unit cell.
To calculate the free energy we need to consider the total Euclidean action, I T ot = I + I bdr , where I = −iS, with S as in (2.1), and I bdr given by the following integral on the boundary r → 0: Here K is the trace of the extrinsic curvature of the boundary and γ µν is the induced boundary metric given in (3.27) up to a factor of r 2 + /(4r 2 ) and the dots refer to a term quadratic in the scalar field which does not play a role for the solutions that we consider in this paper. To obtain the free energy density, w, we write the total free energy as T [I T ot ] OS ≡ wvol 2 .
It is similarly straightforward to obtain the expectation values for the boundary stress tensor, T mn , and the abelian current vector, J m , as well as the condensate O ψ (we have presented a few details in appendix A). Since the solutions have a timelike Killing vector ∂ t , we have where the bar refers to a period average over the spatial coordinates e.g.
The free energy density w depends on T, µ, B, the integer n (giving the number of vortices per unit cell), and also on the shape of the lattice, which can be specified by v, k, where we recall k was defined in (3.29): w = w(T, µ, B; v, k; n). To see how the free energy depends on varying B, v and k we follow [27]. A calculation 1 shows that we can write In order to find the thermodynamically preferred black holes for given UV data (T, µ, B) and given n, we need to minimise w over the lattice parameters (v, k). From (3.41) the preferred vortex lattice will therefore satisfyTxx =Tỹỹ, Txỹ = 0. Since the symmetry of the stress tensor implies that 1 we also have Txỹ = 0. Thus, we deduce that where we also used tracelessness of the stress tensor. As in [27] we can also conclude thatQ To see this, it is convenient to work in thex,ỹ coordinates and consider the free parameters to beL x ,L y and v. As we vary these parameters we vary the boundary metric (3.23) and we also notice that the boundary gauge field in (3.24) is given by A ∂ = µdt − 2πn qỹ dx and hence is independent of these parameters. From (2.13) of [27] we deduce and we then use (3.29).
These conditions are a result of the fact that in thermal equilibrium the spatial parts of the heat current, Q i , and the abelian current, J i , are magnetisation currents of the form is the local magnetisation density and the thermal magnetisation density, respectively. A subtlety is that in a superfluid state we can haveJ i = 0 but withQ i = 0, but for the thermodynamically preferred configurations we haveJ i = 0 [28].
Finally, we note that using the above results, the first law for the thermodynamically preferred vortex lattice black holes can be written as In our numerical results, to be discussed next, we have directly verified that (3.42)-(3.45) all hold.

Numerical Results
In the solutions we have numerically constructed which we discuss below, we have fixed 2 q = 2, B/µ 2 = 0.01 . We express physical quantities in terms of the chemical potential µ; this can be achieved by setting µ = 1 and then reinstating µ afterwards using dimensional analysis. For this value of |qB/µ 2 | the critical temperature at which the dyonic AdS-RN black hole become unstable to zero modes associated with the lowest Landau level is T c /µ = 0.088. We constructed the new branch of vortex lattice black holes that exist below T c which are parametrised by three continuous parameters, the temperature, T , and two parameters, v, k, which determine the shape of the lattice (recall k was defined in (3.29)). We also note that the plots of the free energy that we present below are invariant under k → 1/k for reasons discussed at the end of section 3.4. In principle, the solutions are also specified by a discrete parameter n that appears in the flux quantisation condition (3.28), but all of the solutions we discuss below 3 have n = 1.
We have constructed several black hole solutions, associated with the lowest Landau level, for various values of the parameters T, v, k, but we studied in more depth three specific values of the temperature: T /µ = 0.08, T /µ = 0.07 and T /µ = 0.05 with T /T c ∼ 0.905, 0.79 and 0.565, respectively. For each of these temperatures we determine the shape of the thermodynamically preferred configuration by calculating the free energy density w/µ 3 as a function of (v, k). For all these temperatures, within numerical precision, we find that the preferred black holes have k = 1 and v = 1/ √ 3, corresponding to a triangular vortex lattice. This is clearly illustrated in figure 2 for T /T c ∼ 0.79 and T /T c ∼ 0.565. We can also calculate various physical observables in the dual field theory for the preferred triangular vortex lattice black holes. In figure 3 we present the plots for the case T /T c ∼ 0.79. We see that the modulus of the order parameter, |O ψ | , has zeroes at the centre of each vortex, as expected. Furthermore, we have explicitly checked that the phase of the order parameter winds exactly once around each vortex i.e. n = 1. The spatial components of the local current density J i are magnetisation currents (i.e. with no net transport,J i = 0) and circulate around each vortex core. Since we are considering a superfluid these currents do not diminish the background homogeneous magnetic field, which is the constant value of B throughout. The plots also show that the energy density T tt and charge density J t are slightly diminished at cores of the vortices; it would be interesting to know the reason for this phenomenon. Here |O ψ | is the order parameter, J i is the local current density, T tt is the energy density and J t is the charge density. Note that these quantities are plotted in the original, untransformed spatial coordinates (x, y), scaled by L x and L y (see (3.1) and (3.22)).
It is also interesting to examine the structure of the preferred black hole solutions themselves. In particular, the black hole horizons have an inhomogeneous geometry, with spikes that are associated with the position of the vortices. A complementary picture can be illustrated by studying the modulus of the scalar field, |ψ h |, on the horizon, as illustrated in figure 4(a). At the core of the spikes, the value of |ψ h | is zero. Furthermore, the maximum value, |ψ h | max , which is reached at the midpoint between two spikes monotonically increases as the temperature is lowered, as shown in figure  4(b), leading to more pronounced spikes. Our results lead us to conjecture that at zero temperature the horizon will have a singular spiky structure with |ψ| interpolating between |ψ| = 0 at the core of the spikes to a maximum value of |ψ| = 1, the value associated with the symmetry breaking AdS 4 solution, between the vortices.  For T /µ ∼ 0.085 the original dyonic AdS-RN black hole solution becomes unstable to zero modes associated with the first Landau level (i.e. j = 1 in figure 1). Thus below this temperature, there is an additional family of vortex lattice black hole solutions, again parametrised by the temperature T and two shape parameters v, k (and B fixed as in (4.1)). We have constructed some examples of these back-reacted black holes for T /µ = 0.08 and various values of k and v and we find, surprisingly, that the preferred black holes in this branch appear to have v = 0 and k → 0 or k → ∞, corresponding to the unit cell becoming infinitely long and thin (see figure  5). This suggests that the vortex lattice might be trying to form a linear defect. However, we also find that the black holes associated with the first Landau level that we have constructed are never thermodynamically preferred when compared with any of the black holes that are associated to the lowest Landau level that we discussed above and, in particular, the preferred triangular vortex lattice.

Discussion
For a specific holographic model we have constructed fully back-reacted black holes that are dual to a lattice of vortices in a superfluid phase in the presence of a homogeneous, constant magnetic field. We have shown that the thermodynamically preferred black holes describe a triangular lattice of vortices and associated with the lowest Landau level, at least for the temperatures we have considered.
Assuming that this picture persists to lower temperatures, it seems likely that the core of each vortex on the horizon approaches a dyonic AdS 2 × R 2 configuration and these are embedded in an ambient sea of the IR AdS 4 superfluid ground state, in the presence of a uniform magnetic field. We think it is unlikely that this putative horizon configuration will itself be an exact solution to the equations of motion 4 . This is simply because while the dyonic AdS 2 × R 2 configuration is a solution, the AdS 4 superfluid ground state with a uniform magnetic field is not. On the other hand, it may be possible to construct a kind of inside-out solution which describes a bubble of AdS 4 in the superfluid ground state, embedded in an ambient dyonic AdS 2 × R 2 solution.
There are many natural ways to modify the constructions reported in this paper. While we suspect that we will get similar results by simply varying q and B, this may not be the case and other shapes of vortex lattices could in fact be preferred. We think it would be particularly interesting to extend the model to include a neutral scalar a with an aF ∧F coupling in the bulk Lagrangian. On the one hand such a term arises in the top down models considered in [30,31] which also explored a fascinating competition between superfluid phases and phases in which the U (1) symmetry is unbroken but spatial translations are spontaneously broken. In addition given the fact that in field theory a Chern Simons term gives rise to novel vortices carrying both electric and magnetic charges (e.g. see [32]), which are models for anyons, we can expect that the associated gravitational construction 5 may well exhibit novel phenomena too.
The thermodynamically preferred black holes that we constructed in this paper are associated with vortices in the lowest Landau level. However, this is not expected to be the case for other models such as the ones considered in [34,35]. These black holes could have a very interesting structure based on the fact that the preferred black hole vortex lattice in the first Landau level which we constructed here has an interesting linear defect structure.

A Asymptotic Expansion
Using the modified equations of motion as described in section 3.4, the asymptotic boundary expansion for the functions appearing in the ansatz (3.30) take the form Q tt (r,x,ỹ) = 1 + r 3 c tt (x,ỹ) + g 1 (x,ỹ)r (3+ The appearance of the function c tt in the expansion for Q is associated with the fact that the stress tensor is traceless. We have sixteen functions ofx,ỹ which are fixed in the numerical integration. The functions g 1 and g 2 have been seen in similar DeTurck constructions before [22,36], and can be removed by a gauge transformation when the modifications to the equations of motion vanish, ξ µ = ϕ = 0. Similarly, c tr , c rx , c rỹ and c r can be removed by gauge transformations, if desired; they don't appear in the expressions for physical quantities in the boundary theory. In order to calculate the one point functions we need to supplement the bulk action (2.1) with the boundary action given in (3.37) (in the Euclidean frame). The expectation value of the stress tensor, T µ ν , can be obtained from the expression where g refers to the bulk metric, and we find Notice that these expansions imply the Ward identity T t t + Txx + Tỹỹ = 0. Furthermore, the expectation values of the abelian current density are given by and we find The expectation value for the operator, O ψ , dual to the complex scalar, ψ = ψ 1 + iψ 2 , can also be determined 6 . Writing |O ψ | for the modulus of O ψ we find Similarly the sine of the phase of O ψ is given by c 1 / c 2 1 + c 2 2 and we find that the phase winds once around the core of each vortex in the solutions we have constructed, i.e. n = 1.
At fixed µ, B, T , within numerical precision we find that for all of the black holes we have constructed we havē c tx =c tỹ = 0 ,cx =cỹ = 0 , (A. 7) and this corresponds to the fact that in thermal equilibrium the black holes have vanishing current fluxes: Furthermore, for the thermodynamically preferred black holes, identified by varying the free energy density w with respect to v and k (defined in (3.29)), within numerical precision we find for all black hole solutions This corresponds to the conditions

B Numerical Convergence
To solve the system of PDEs discussed in section 3.4, with appropriate boundary conditions, we discretised the computational domain [0, 1]×[0, 1]×[0, 1] using N r , N x , N y points, respectively. We approximated the derivatives of the fields on these points using an interpolation method: since thex coordinate is periodic, we used a Fourier spectral method, while for r andỹ we used the Chebyshev spectral method (recall that we have to implement the quasi-periodic boundary conditions under shifts of y →ỹ + 1 as in (3.25)). The resulting algebraic system of equations was solved iteratively using the Newton-Raphson method. For all plots given in the paper we have taken N r = N y = 50 and N x = 46, unless stated otherwise. We now discuss the convergence of our code with increasing lattice resolution, especially with the number of point N r that we take in the radial direction. One reason this is important is because we have no analytic proof that the only possible solutions of the DeTurck system of equations which we solve actually reduce to Einstein's equations when we impose our boundary conditions. We therefore need to check that ξ µ becomes trivial everywhere in our domain in the continuum limit. This is unlike the case of the auxiliary field ϕ which we introduced in (3.35) in order to fix gauge invariance. In fact, we can show that ϕ must be trivial and it can therefore serve as a check for the convergence of our code. In figures 6 we have plotted ξ 2 max , the maximum value of ξ 2 that we find in our computational domain. Fitting that to a power law we find the behaviour ∼ N −13.7 r . We found it reassuring that in all our solutions we had |ϕ| < 10 −8 . However, this error is close to our numerical precision and so we are not able to check its convergence properties.  Figure 6: Plot of ξ 2 max , the maximum value of ξ 2 , as a function of the radial resolution N r , for the preferred triangular vortex lattice black holes (lowest Landau level) at T /µ = 0.08, B = 0.01. Here N x = 60, N y = 50.
A second reason that the details of the convergence of our code is important concerns the non-analytic terms related to g 1 and g 2 which appear in the near conformal boundary expansion (A.1). In general, such terms turn the exponential convergence of spectral methods to power law with respect 7 to the radial direction (as we saw above for the behaviour of ξ 2 ). For the purposes of our checks, we have found useful to consider the quantities, Using the expansion (A.1), we can easily see that in the continuum we must have c tt = c (1) = c (2) . Notice that this quantity is particularly important since it enters the expression for the free energy of the system via (A.3) and (3.38). From the expansions (A.1), we see that the leading non-analytic power will drop out from c (1) and we therefore expect it to have a better convergence rate than c (2) . In figure 7 we plot c (2) (N r + 5) − c (2) (N r ) as a function of the points in the radial direction, N r . We see that this quantity has power law convergence to zero, with behaviour ∼ N −3.7 r . Our expectation is that c (1) (N r + 5) − c (1) (N r ) would converge with a smaller power. However, we were not able to confirm this since, due the 7 On the one hand this means that the spectral methods which we use in the radial direction become inefficient at high resolutions. On the other hand, however, in terms of memory usage, we have found them to be very efficient in obtaining sufficiently accurate solutions. Indeed it would be significantly more demanding in resources in order to achieve similar accuracy using finite difference methods in the radial direction. It is also worth noting that the precise power law behaviour which we find is for a specific choice of reference metric in (3.34), namely the AdS-RN black hole solution. It would be interesting to explore whether other choices leads to improved convergence. absence of the leading non-analytic terms and faster convergence, it becomes smaller than our numerical precision quite quickly as we increase N r , before it takes the asymptotic form of a power law 8 . It is worth highlighting that these tests for c (1) c (2) are associated with the DeTurck equations we are solving and are not related to the fact that we are not solving Einstein's equations. On the conformal boundary we don't impose the equations of motion but our boundary conditions. We therefore need to check whether the equations of motion are asymptotically satisfied by numerically checking the series expansion (A.1) of our functions.
The main error in the various tests we have discussed arises from different places in the bulk. For figure 7 the main error comes from the near conformal boundary region, as one might expect from the discussion above. On the other hand, the error for the plot in figure 6 arises from a more global test and, in fact, mainly arises from the near horizon region where ξ 2 and |ϕ| take their maximum values.
We have implemented our numerical method in C++ using double precision numbers and class data structures. The code is fully parallelised in shared and distributed memory using a combination of MPI and OpenMP. After fixing a discretisation scheme, the problem reduces to solving a set of non linear algebraic equations for the values of our functions on the grid described above. This is done using the Newton-Raphson method where one starts with an initial guess for the unknown functions at each lattice point and then iteratively corrects it in order to obtain functions that solve the PDEs to a better and better approximation. This boils down to solving a linear system at each step in Newton's method. To do this we use the PETSc library with a block-ILU(0) preconditioner in combination with a GMRES iterative method.