Minimally packed phases in holography

We numerically construct asymptotically AdS black brane solutions of D = 4 Einstein-Maxwell theory coupled to a pseudoscalar. The solutions are holographically dual to d = 3 CFTs at finite chemical potential and in a constant magnetic field, which spontaneously break translation invariance leading to the spontaneous formation of abelian and momentum magnetisation currents flowing around the plaquettes of a periodic Bravais lattice. We analyse the three-dimensional moduli space of lattice solutions, which are generically oblique, and show, for a specific value of the magnetic field, that the free energy is minimised by the triangular lattice, associated with minimal packing of circles in the plane. We show that the average stress tensor for the thermodynamically preferred phase is that of a perfect fluid and that this result applies more generally to spontaneously generated periodic phases. The triangular structure persists at low temperatures indicating the existence of novel crystalline ground states.


Introduction
Spatially modulated phases, in which spatial symmetries are spontaneously broken, are seen in a wide variety of materials in Nature. For example the cuprate superconductors exhibit spin and charge density waves as well as current loop order [1]. Starting with [2,3] it has become clear that a wide variety of spatially modulated phases can be realised in strongly coupled field theories using holographic techniques. These phases are associated with novel black hole solutions in AdS spacetimes which have spatially modulated horizons, 1 with the first example presented in [5]. Exploring the landscape of such black hole solutions and elucidating their possible zero temperature crystalline ground states is interesting and still largely unexplored territory.
In this paper we continue to investigate a class of models in D = 4 spacetime dimensions, which includes the top-down models of [6], that couple a metric to a gauge-field and a pseudoscalar field. The models all admit an AdS vacuum which is dual to a CFT in d = 3 spacetime dimensions with an abelian global symmetry. When the CFT is held at finite chemical potential, µ, with respect to the abelian symmetry, the unbroken spatially homogeneous and isotropic phase in flat space is dual to the electrically charged AdS-RN black hole (brane). Depending on the specific couplings of the model, these black holes can have striped instabilities below a critical temperature T c [7]. In particular, at T = T c a JHEP03(2016)148 static linearised striped mode appears which reveals the existence of new spatially modulated black hole solutions which, furthermore, are associated with abelian and momentum current density waves. Fully back reacted striped black hole solutions, modulated in one spatial direction and with the current density waves running along the orthogonal spatial direction, were subsequently constructed in [8][9][10][11][12] by solving a system of PDEs in two variables.
At the linearised level the static striped modes can be superposed indicating the existence of black holes with spatial order in two spatial directions. By solving a system of PDEs in three variables, black hole solutions dual to checkerboard lattices were constructed in [13], with the currents now running around the unit cell of the checkerboard. It was also shown that the striped black holes are continuously connected to the checkerboard black holes via rectangular lattice black holes. Here we will explain how the checkerboard and rectangular lattice black holes associated with the striped instability of [7], are special cases of a more general three-dimensional moduli space of spatially modulated lattice black holes, defined by a two dimensional Bravais lattice. In general these lattices are oblique and special cases are rectangular, checkerboard, centred rectangular and triangular lattices.
When there is just a chemical potential, it was shown in [13] that the striped phase is thermodynamically preferred over the checkerboard and rectangular lattices (at least for the temperatures and lattices considered in [13]). In order to obtain two dimensional phases that are preferred one can supplement µ with an additional UV deformation. In [13] it was shown that after adding a homogeneous source, φ 0 , for the operator dual to the pseudoscalar, that the checkerboard phase can be preferred over the striped phase. It would therefore be interesting to further investigate whether other lattices are preferred over the checkerboard in this context.
In this paper we will not pursue that line of investigation but instead consider the CFT with µ = 0 and also a uniform uniform magnetic field 2 of strength B. It was shown in [22] that the striped instability of [7] can survive when B = 0. For a specific value of B/µ 2 we will construct and explore the full three-dimensional moduli space of oblique lattice black holes. For all these black holes the dual CFT has abelian and momentum magnetisation currents flowing around the unit cells. The transition from the unbroken phase to the spatially modulated phase is a first order transition and, interestingly, amongst all of the lattices we find that the triangular lattice is the thermodynamically preferred structure. The triangular lattice is a commonly observed lattice since it is associated with minimal packing of circles in the plane. However, why it is preferred in the present context is not so clear given that we have a strongly coupled plasma. We have cooled the black hole solutions down to quite low temperatures and our results indicate the existence of fundamentally new types of holographic crystalline ground state solutions at T = 0.
In section 2 we introduce the class of models that we study as well as discuss the black hole solutions that describe the unbroken phase. Section 3 reviews the striped instability while section 4 discusses a convenient parametrisation of the moduli space of oblique lattices. Section 5 discusses the construction of the lattice black hole solutions and it is shown JHEP03(2016)148 that the thermodynamically preferred solutions are dual to a triangular phase. In carrying out our analysis, extending the results of [23], we show that the average stress tensor is that of a perfect fluid and we argue that this is a general result for spatially modulated phases. Some technical material, including some details on our numerical implementation, are presented in two appendices.

The setup
We consider a bulk theory in D = 4 which couples the metric to a gauge-field, A, and a pseudoscalar field, φ, with action given by yielding the equations of motion where F = dA is the field strength of the gauge field. In the coordinates we choose later we take trxy = √ −g. We have taken 16πG = 1 for convenience. The numerical solutions we construct in this paper will be for specific choices of the functions V, Z ≥ 0, ϑ given in (2.9), below. For the moment, though, we just impose some rather weak conditions on these functions. We demand that φ is a pseudo-scalar field, which is achieved if V, Z are even functions and ϑ is an odd function of φ: We also assume that the equations of motion admit an AdS 4 vacuum with φ = A = 0, which is dual to a d = 3 CFT with an abelian global symmetry and the scalar field is dual to a pseudoscalar operator O φ , with scaling dimension ∆ φ . For convenience we will choose V (0) = −6 so that the AdS 4 vacuum has unit radius.
Our primary focus will be the thermal phase diagram for the CFT with uniform chemical potential µ and uniform (constant) magnetic field B. The F ∧ F coupling to φ in the action implies, generically, that the scalar field will get activated in the dual black hole solutions. In the case that the scalar is dual to a relevant operator, i.e. ∆ φ < 3, we can also consider an additional deformation parameter, φ 0 , associated to sourcing this operator.
At high temperatures, T /µ, T /B 1/2 , T /φ 1, the CFT in flat space will be in an unbroken, homogenous phase which preserves the Euclidean symmetries of the two JHEP03(2016)148 spatial directions. The associated dual black hole solutions are captured by a radial ansatz of the form Notice that this does not fully fix the reparametrisation invariance r → R(r). For our purposes we found it convenient to choose with r + a free constant which sets a scale in the dual field theory. In these coordinates, the conformal boundary is at r = 0. Assuming that f, G → 1 as r → 0, writing r = r + /2 and taking → 0 we find that If we also expand a t = µ + . . . and φ = φ 0 3−∆ φ + . . . , then we identify µ, B, φ 0 as the chemical potential, magnetic field and source for the scalar operator, respectively. The location of the black hole event horizon is taken to be at r = 1 and we impose the following expansions (2.7) After the Wick rotation t → iτ we find that this leads to analytic behaviour provided that we periodically identify τ with period β = 4 π r + c f , associated with a temperature T = β −1 . In the special case that B = φ 0 = 0, the class of theories satisfying (2.3) admits the exact electric Reissner Nordstrom black brane solution given 3 by . We recall that as T → 0 this solution approaches AdS 2 × R 2 in the far IR. For general µ, B, φ 0 , the black holes in the ansatz (2.4) with boundary conditions (2.6), (2.7) need to be constructed numerically, which can be carried out using standard shooting techniques (or the techniques outlined in the remaining of the paper). The specific choices of the functions V, Z, ϑ that we will focus on in this paper are given by where (s, χ) are constants. We will quantise the pseudoscalar field so that ∆ φ = 2. If s = χ = 1 then this is the top-down model discussed in [6] (after rescaling the metric, gauge-field and Newton's constant).

JHEP03(2016)148 3 The striped instability
It was shown in [7] that when B = φ 0 = 0 and µ = 0, the electric AdS-RN solution (2.8) has spatially modulated instabilities below some critical temperature in broad families of theories satisfying (2.3). The striped instability involves both an electric and momentum current density wave in a direction that is orthogonal to the direction of a spatially modulated expectation value for the pseudoscalar operator O φ . A simple way to understand the instability is to analyse perturbations about the T = 0 AdS 2 × R 2 solution in the far IR, seeking modes with non-trivial dependence on the spatial coordinates of R 2 that violate the AdS 2 BF bound. For the class of models given in (2.9) the striped instabilities depend on (s, χ) as displayed 4 in figure 1 of [7] and in particular the top-down model with s = χ = 1 is unstable. When B and φ 0 are non-zero, it was shown in [22] that there can be a more general modulated instability that also involves modulated charge and energy density waves in addition to the modulated currents and pseudoscalar. The argument in [22] for the existence of this instability relied on the fact that the near horizon geometry of the extremal black holes is still AdS 2 × R 2 , which is expected to occur in a regime of small deformations, i.e. where B 1/2 µ and φ /µ is heavily model dependent. In particular, for a given model it is possible that the T = 0 ground state of the unbroken phase black holes is no longer AdS 2 × R 2 and moreover that the new IR geometry is perturbatively stable against such modes. An explicit example of this possibility was realised in [13] where it was also shown that, nevertheless, for small enough deformations there can still be finite temperature spatially modulated instabilities.
We now continue by analysing in more detail the static, perturbative modes about the unbroken phase black hole solutions (2.4) that correspond to the striped instability. Specifically, in a particular gauge we consider where the indices i, j run over the spatial coordinates (x, y) and In addition g tt and g rr are the background metric components in (2.4). By adding this perturbation to the background (2.4) and substituting into the equations of motion for (2.1), at linear order we are led to a system of second order ODEs for {h ⊥⊥ , h t⊥ , h t , h ⊥ , h} and first order ODEs for h tt , h . A solution is therefore specified by twelve integration

JHEP03(2016)148
constants. For the case of the AdS-RN black hole solution (2.8) (with B = φ = 0), we go back to the case studied in [7] with only the set {h t⊥ , h ⊥ , h} being non-trivial. Demanding regularity of the perturbation close to the horizon gives the analytic expansion at r = 1: where {c µν , c ty , c µ , c} are constants of integration to be determined. In the absence of sources, with ∆ φ = 2, the corresponding UV expansion near r = 0 is given by where {d µν , d µ , d} are additional constants, out of which d tt and d xx are fixed in terms of the others. In total we have seven constants from the near horizon expansion and another five from the asymptotic region. Since we are dealing with a linear system of equations, we can scale the unknown functions to set any of the twelve non-zero constants to one. The extra constant needed to uniquely specify a solution for fixed values of µ, k , B and φ 0 is the critical temperature, T c /µ, at which the static mode comes into existence, and below which the black holes are unstable. We mentioned above that for the electric AdS-RN black hole the top down model with s = χ = 1 is unstable. For a given value of s, the critical temperature for the instability is an increasing function of χ. For s = 1, the top-down value χ = 1 is very close to a critical value χ c below which the unstable modes completely disappear. In the rest of the paper we will therefore take s = 1 and χ = 1.5 in order to improve the numerics, but we expect that the results for the top-down values will be similar. Setting the scalar deformation to zero, φ 0 = 0, we have numerically constructed black holes dual to the unbroken phase with µ = 0 and B/µ 2 = 0.05. One finds that the T = 0 limit of these black holes approaches a dyonic AdS 2 × R 2 solution in the IR and hence, using the results of [22], has striped instabilities. We also constructed the critical temperatures at which the striped instability appears and the results are presented in figure 1.

The moduli space of oblique lattices
Before moving on to the discussion of the back reacted solutions that are associated with the striped modes given in (3.1) we first discuss the moduli space we wish to explore. It is clear that at the critical temperature T c , associated with a fixed value of the magnitude of the wave-number, k = k , there is a one parameter family of such striped modes, related by a rotation in the (x, y) plane. At the linearised level we can consider an arbitrary linear superposition of these modes. However, demanding periodicity in the (x, y) plane, so that the variational problem for the back reacted solutions is well-posed, restricts the superposition to be the sum of two wave-number vectors, each with magnitude k and in general non-parallel. Thus, the resulting set of zero modes at T c is parametrised by an angle between the two vectors.
By including the back reaction we therefore expect a family of periodic spatially modulated solutions to appear at T c . For lower temperatures the moduli space of periodic black hole solutions will be parametrised by two wave-numbers k 1 and k 2 , which no longer need to have the same magnitude. In other words we allow the spatial periodicities to change as we lower the temperature. It is always possible to perform a rotation to choose, for example, k 1 to be parallel to the y-direction. This leaves us with a three parameter moduli space of black hole solutions to explore. In the special case that the k i are parallel, with same magnitude, we have stripes, otherwise we will have a two-dimensional Bravais lattice. Generically these are oblique periodic lattices, but there are various special cases: rectangular lattices, checkerboards, centred rectangular lattices and hexagonal lattices. The latter is also known as a triangular lattice and is the lattice associated with minimal packing of circles in the plane. The task is to explore this three-dimensional moduli space of black holes, and identify the configuration that minimises the free energy at each temperature T ≤ T c .
It is clear from this discussion that the back-reacted striped black hole solutions constructed in [8][9][10][11][12] have additional zero modes associated with the larger moduli space of solutions corresponding to the two-dimensional lattices. For the special case of µ = 0, φ 0 = 0 and B = 0, checkerboard lattices were constructed in [13] but were shown to have free energy that is sub-dominant to the striped phase. Some additional rectangular lattices were also constructed in [13] for this case and also shown to be sub-dominant. As an aside, as far as this paper is concerned, we note that some preliminary investigations which we have carried out indicate that the striped phase is also preferred over all oblique lattices. It was also shown in [13] that when µ = 0, φ 0 = 0 and B = 0, the checkerboards can dominate over the stripes and the transition is first order. It remains an interesting open problem to revisit the analysis of [13] and explore the full three-dimensional moduli space of periodic lattices when φ 0 = 0 and B = 0.

JHEP03(2016)148
In the next section we instead construct back reacted black hole solutions with µ = 0, φ 0 = 0 and B = 0 and explore the full moduli space of solutions. Before doing that we briefly discuss a convenient parametrisation of the periodic configurations. Let us now denote the spatial coordinates (x , y ) associated with the asymptotic boundary metric Consider wave numbers with π/2 ≤ α ≤ π. The associated periodic boundary conditions in the (x , y ) directions are given by where n i are integers. The parameters L x , L y and α are the three moduli for the space of solutions that we would like to construct. We also note that in the plots involving the spatial coordinates that we present later we use (x,ŷ) ≡ ( x Lx , y Ly ). An equivalent parametrisation, which we shall use in the next section, is to impose the simple periodic boundary conditions but demand that the conformal metric asymptotes to Indeed the two coordinate systems are related by Instead of (L x , L y , α) we can also parametrise the moduli space of solutions using (k x , k y , R 0 ) with k x ≡ 2π/L x , k y ≡ 2π/L y , R 0 ≡ − cot α and R 0 ∈ [0, ∞). Notice that R 0 = 0 (α = π/2) corresponds to a rectangular lattice and when k 1 = k 2 it becomes a checkerboard. As R 0 → ∞ (α = π) we obtain stripes. The hexagonal, or triangular lattice, is given by k 1 = k 2 and R 0 = 1/ √ 3 (α = 2π/3). All lattices can be obtained by allowing α to lie in the range π/2 ≤ α ≤ π but there is some redundancy. Generically we have an oblique lattice. There are the following special cases with larger symmetry. If k x = k y we have rhombic lattices; if in addition α = 2π/3 it is hexagonal while if α = π/2 it is square. If α = π/2 we have rectangular lattices; if in addition k x = k y we have a square lattice. If k y /k x = −2 cos α we have centred rectangular lattices satisfying |k y | 2 = −2k y ·k x ; if in addition α = 2π/3 it is the hexagonal lattice, while if α = 3π/4 (i.e. k y = √ 2k x ) it is the square lattice rotated by π/4 compared to the previous two. Finally there are also similar centred rectangular lattices with k x /k y = −2 cos α.

JHEP03(2016)148 5 Lattice black holes solutions
We are now in the position to numerically construct the back reacted black hole solutions for the periodic lattices that we described in the previous section. To achieve this, we use the following ansatz which has only a time-like Killing vector: where {Q tt , Q rr , Q, R, M, Q tx , Q ty , Q tr , Q rx , Q ry } and {a t , a r , a x , a y , ϕ} are all functions of r and periodic functions of x = (x, y) with (x, y) ∼ (x + 1, y) ∼ (x, y + 1). The functionŝ f (r) and g(r) that appear in the ansatz (5.1) are fixed by hand for convenience (cf. (2.8)) and their explicit form is taken to bê For convenience we will chooseL Of course, there is a lot of redundancy in this ansatz due to local coordinate and gauge invariance of the theory (2.1) and this will be dealt with below. The boundary conditions near the AdS 4 boundary, located at r = 0, that we want to impose are given by and we will only consider φ 0 = 0, the case of no scalar deformation. Notice that these boundary conditions imply that the asymptotic metric approaches (4.1) and we therefore demand that the coordinates x and y have unit periods, as in (4.4). Also observe that the B dependent part of the field strength for the gauge field takes the formL xLy Bdx ∧ dy = Bdx ∧ dy in the coordinates associated with (4.1) and (4.3). We will demand that there is an analytic Killing horizon generated by the Killing vector ∂ t , located at r = 1. Analyticity demands that our functions will only depend on JHEP03(2016)148 even powers of (1 − r), leading to the Neumann boundary condition ∂ r Φ(1, x, y) = 0 for all functions appearing in the ansatz (5.1) apart from Q tt , for which we impose Q tt (1, x, y) = Q rr (1, x, y).
In order to get a well defined boundary valued problem for an elliptic, quasi-linear system of partial differential equations, as well as dealing with the residual coordinate and gauge invariance, we follow the approach of [24][25][26][27] by suitably modifying the equations of motion. After solving the modified equations, we check that the solutions in fact solve the unmodified equations.
For the Einstein equations we add the term ∇ (µ ξ ν) to the left hand side of Einstein's equation in (2.2). As in [24] we choose ξ µ = g νλ Γ µ νλ −Γ µ νλ where the Christoffel symbol Γ µ νλ is with respect to a reference metricḡ µν which needs to have a Killing horizon and the same asymptotics with the solution we would like to construct. In the case of pure gravity and non-positive cosmological constant, one can actually show [25] that given such a vector, there is no non-trivial solutions to the modified Einstein's equation with ξ µ nontrivial. For the more general kind of theories we are considering, though, a corresponding theorem is not yet available and one has to numerically check that ξ µ = 0 a posteriori. The vanishing of ξ µ fixes the choice of coordinates.
Following a similar logic for the gauge field equation of motion, we will introduce a scalar ψ and, as in [13], add the term Z ∇ ν ψ to the left hand side of Maxwell's equation in (2.2). With the choice of ξ µ above, we choose 5 ψ to be given by By linearising the metric and the gauge-field about a given configuration one can check that this modification leads to an elliptic set of equations. In fact this is true whether or not we include the terms involvingĀ: the utility of these terms is that they allow one to choose A µ −Ā µ conveniently. In the present set up we will chooseĀ =L xLy B 2 (x dy − y dx) so that the components of A µ −Ā µ are periodic functions and so that ψ = 0 on the boundary. We will also choose the reference metricḡ µν to be given by (5.1) with Q tt = Q rr = Q = M = 1, R = R 0 , and vanishing {Q tx , Q ty , Q tr , Q rx , Q ry , }. Notice that ψ is now a periodic function of x which vanishes at the AdS boundary and satisfies Neumann boundary conditions at the black hole horizon.
If we take the covariant derivative of the modified Maxwell equations we deduce that We now multiply this equation by ψ and integrate over r, x. Integrating by parts we find that the boundary contributions vanish and, with Z ≥ 0, we conclude ∇ψ = 0. Since ψ = 0 on the AdS boundary we conclude that ψ = 0 everywhere. Similarly to the metric, vanishing of ψ fixes the gauge invariance.
To summarise, if we solve our modified equations we just need to check that ξ µ = 0 to ensure that we have solved the equations of motion (2.2). Since ξ µ is spacelike, this is JHEP03(2016)148 achieved by checking ξ 2 = 0. Some details about the asymptotic expansion for the modified equations are included in appendix A and some discussion of the numerical implementation and convergence can be found in appendix B.

Currents and free energy
We can extract the current and stress tensor of the dual field theory from the asymptotic behaviour of the metric and gauge-field. As shown in appendix A, and following [28], the spatial components of the abelian and momentum currents are magnetisation currents of the form where M ij = M [ij] is the local magnetisation density and M ij T is the local thermal magnetisation density. This immediately implies that the average current fluxes vanish J i =T ti = 0, where here, and in the following, the bar refers to a period average e.g.
Here we are using the periodic coordinates given in (4.3), √ γ refers to the flat spatial part of the metric, γ ij , in (4.1) and we have used the fact that the volume of the unit cell on the torus is given by sin α/(L x L y ). In all of the black hole lattice solutions that we construct there are both abelian and momentum currents flowing around the unit cells as in figure 4. In order to calculate the thermodynamically preferred black hole solutions we need to determine the free energy. The free energy density is given by where I OS is the total Euclidean on-shell action, including contributions from the Gibbons-Hawking terms as well as standard counter terms. We note that for the black holes of interest we have w = w(T, µ, B; k 1 , k 2 , α). For a fixed type of periodic lattice. i.e. fixing (k 1 , k 2 , α), we have where s is the entropy density (i.e. entropy per unit cell area),J t is the average charge density and m is the magnetisation density given by the bulk integral, as in [22], In appendix A we derive an expression for the local magnetisation density, M ij , which, consistent with (5.7), has an ambiguity of the addition of a constant times ij . This constant can be fixed, thermodynamically, by demanding that m =M x y .

JHEP03(2016)148
If we instead hold fixed the UV data (T, µ, B) and vary (k 1 , k 2 , α), by following the arguments of [23], we find after switching from (x , y ) to (x, y) coordinates whereT is the period average of the stress tensor. In addition, we also have For the thermodynamically preferred black hole solutions, all of the derivatives in (5.12) vanish, and it is straightforward to show that we have where γ ij is the spatial part of the flat metric given in the (x, y) coordinates in (4.1) and the constant p is the average pressure. We mentioned earlier thatT ti = 0. It is interesting to note that in spite of the spatial modulation the average stress tensor for the thermodynamically preferred phase is that of a perfect fluid. 6 We also notice that for the thermodynamically preferred black holes we have sT +J t µ =T tt + (mB + p).

Results
With B/µ 2 = 0.05 we recall that the critical temperature for the onset of the instability is given by T /µ = 0.075 (see figure 1). We have constructed various oblique lattice black holes for temperatures both greater than and less than T /µ = 0.075. Since some oblique lattice black holes exist for T /µ > 0.075 which are not thermodynamically preferred and since they are all thermodynamically preferred for T /µ < 0.075, it is clear that the transition from the homogeneous phase to the spatially modulated phase is a first order transition. Our principal interest is not to investigate the details of the first order transition (for example the critical temperature), but rather to investigate the shape of the thermodynamically preferred lattices and so we just discuss solutions with T /µ < 0.075. In order to explore the moduli space of black hole solutions, we first consider rhombic lattices ("rhombiboards") for which k x = k y (i.e. L x = L y ). We present some results for the black hole solutions that we constructed for the three different temperatures, T /µ = 0.07, 0.05 and 0.04. In the "dactylograms" of figure 2, each of which was made from about 400 individual black hole solutions, we show the contours of the free energy as a function of k x and R 0 . Figure 2. Free energy of the rhombiboard lattices. The contours display the free energy densitŷ w ≡ w/µ 3 as a function of R 0 andk x =k y , wherek x ≡ k x /µ andk y = k y /µ, for temperatures T /µ = 0.07 (top left), T /µ = 0.05 (top right) and T /µ = 0.04 (bottom). Colours with longer wavelength correspond to contours with larger values. The minimum is at R 0 = 1/ √ 3, corresponding to triangular lattices, in all three cases but the period changes with temperature. The minima for the three cases are located atk x = 0.937,k x = 0.828 andk x = 0.777, respectively.

JHEP03(2016)148
In all three cases, quite remarkably, we find that the minimum is given 7 by R 0 = 1/ √ 3 with the value of k x decreasing as the temperature is lowered. In other words at all three temperatures the triangular lattice are the preferred lattices within the class of rhombic 7 The minimum can be obtained in two ways. One can construct the free energy, stress tensor and magnetisation for the individual black holes and then use (5.12) to construct the derivatives on the left hand side of (5.12) as functions of (kx, α) via an interpolation. The vanishing of these derivatives gives R0 = 1/ √ 3 to about one part in 10 7 . An alternative method is simply to use the free energy density for the individual black holes to obtain a function w of (kx, α) via an interpolation. Setting the derivatives of this function to zero gives R0 = 1/ √ 3 to about one part in 10 4 .

JHEP03(2016)148
Rk yk x Figure 3. Surfaces of constant free energy densityŵ ≡ w/µ 3 in the three dimensional module space (k x ,k y , R 0 ), wherek x ≡ k x /µ andk y = k y /µ. The plot was made for the T /µ = 0.04 case of figure 2. The minimum is for the triangular lattice with R = 1/ √ 3 and equal periodŝ k x =k y , = 0.777.
lattices. As an aside we notice in figure 2 that for each case there is a saddle point at R 0 = 0, associated with checkerboards.
To explore whether the triangular lattice is preferred in the full three dimensional moduli space, we also constructed additional black holes, about 600, at the lowest temperature of figure 2, T /µ = 0.04, in a neighbourhood of the triangular lattices. In figure 3 we see the triangular lattice is indeed a local minimum of the free energy and almost certainly the global minimum.
It is interesting to display the spatial variation of various expectation values for the thermodynamically preferred triangular lattices. The general features are the same for all temperatures that we have considered. In figure 4, for the representative case of T /µ = 0.04, we display the spatial variation of the operator O φ , dual to the pseudoscalar field, the charge density, J t , the energy density, T tt , and also the norm and direction of the spatial components of the current vector J i . It is interesting to notice that O φ and J t have a spiky structure and T tt and J i have dimples at the positions of the spikes.
It would be very interesting to determine the ultimate T = 0 ground state of the spatially modulated phase. It seems likely that it continues to be a triangular lattice for all temperatures. Following the argument in section 3.7.5 of [29] we expect that ∂ T k x → 0 as T → 0 and it is plausible this happens at k x = 0, leading to a crystalline ground state. It is numerically challenging to keep reducing the temperature in the systematic manner that we have described so far. In order to get some additional insight into the ground state, we constructed the thermodynamically preferred triangular lattice at T /µ = 0.01, i.e. amongst the triangular lattices we determined the thermodynamically preferred value of the wavenumber, finding k/µ = 0.6456. In figure 5 we have plotted the spatial variation of the pseudoscalar field on the black hole horizon for both this temperature and also for T /µ = 0.04. The appearance of spikes at the horizon as we lower the temperature, combined with the structure of the magnetisation current displayed in figure 4, indicates that the ground state is a fundamentally new type of holographic solution, which will be studied in more detail in the future.

JHEP03(2016)148
xx yŷ φ + φ + Figure 5. A plot of the spatial dependence of the scalar field on the black hole horizon, φ + ≡ φ(r + ), as a function of the spatial coordinates (x,ŷ) ≡ ( x Lx , y Ly ). The left plot is for the thermodynamically preferred triangular lattice at T /µ = 0.04. The right plot is for the preferred triangular lattice at T /µ = 0.01.

Discussion
In this paper we have numerically constructed co-homogeneity three, asymptotically AdS black holes in four spacetime dimensions. The solutions are holographically dual to d = 3 CFTs at finite chemical potential and in a constant magnetic field, which spontaneously form a periodic lattice in two spatial dimensions, with magnetisation currents of the form (5.7) circulating the plaquettes of the lattice. In addition there is also a commensurate modulation of the charge density. These features, which were also observed in the checkerboard lattice of [13], are somewhat reminiscent of the current loop order that is observed in the cuprates [1].
We showed that the black holes come in three-parameter families which are associated with different Bravais lattices. For a specific value of the magnetic field we showed that the triangular lattice is thermodynamically preferred over all other lattices, at least down to very low temperatures. While this is certainly a very natural lattice to appear, being associated with close packing of circles in the plane, it is unclear what aspects of the gravitational model selects this lattice. Perhaps it is possible 8 to illuminate this issue by developing a Landau-Ginzburg type description near the critical temperature. However, such a description will not be valid at low temperatures. It would be interesting to consider other values of the magnetic field and see if the triangular lattice persists. It would be particularly interesting to examine what happens as we reduce the magnetic field to zero since, based on the work of [13], it seems likely that the striped phase is then the preferred configuration. More generally, it is essentially a wide open question which periodic lattices in two spatial dimensions can be preferred in the context of more general holographic constructions. The possibilities in three spatial dimensions, associated with five dimensional gravitational models, are richer still.

JHEP03(2016)148
It would be interesting to examine the transport properties of the new triangular phase that we have found. Powerful techniques have been recently developed that enable one to extract the thermoelectric DC conductivity from the black hole horizon [28,31,32]. However, since the triangular phase spontaneously breaks translations, there will be Goldstone modes and hence the DC conductivity is expected to be infinite. On the other hand the intricate structure of magnetisation currents in our phase could lead to interesting structures in the AC conductivity. It should be noted, however, that while it is conceptually straightforward to calculate the thermoelectric conductivity by perturbing the black holes, it is technically somewhat involved.
The low temperature limit of the black hole solutions that we constructed also revealed some striking features. The spatial modulation of the triangular lattice does not got washed out, but instead persists. Moreover, the spatial modulation at the black hole horizon starts to develop a spiky structure which certainly deserves to be explored in more detail. Indeed it suggests that the zero temperature ground state is a triangular crystal that is being supported by a novel periodic structure on the horizon.

A Asymptotic expansion of fields
Using the modified equations of motion described in section 5, the expansion of the functions appearing in the ansatz (5.1) near the AdS boundary have the form Q tt (r, x, y) = 1 + r 3 c tt (x, y) + g 1 (x, y) r (3+

JHEP03(2016)148
In total we have fifteen functions of x and y which are fixed by integration. Analogues of the functions g 1 and g 2 first appeared in a closely related context in [33]. The condition ξ µ = 0 implies a linear relation between g 1 and g 2 and they can then be removed after doing a coordinate transformation. Similarly, the functions c tr , c xr , c yr and c r can also be removed by coordinate and gauge transformation. In particular, these functions do not appear in the expressions for physical expectation values of the dual field theory. The components of the abelian current vector can be expressed in terms of the above expansion via Similarly, the components of the stress energy tensor are given by In obtaining these results we note that to get the boundary metric in the standard form (2.6) one needs to scale r = r + /2. Note that for the thermodynamically preferred black holes, satisfying (5.14), we should havec which we have checked in our numerics.

A.1 Magnetisation currents
Following [28] we show that the spatial components of the abelian and momentum currents in the dual CFT are magnetisation currents and hence the average current fluxes vanish: and observe that when evaluated at the AdS boundary J a bulk is the conserved current density, √ γJ a , of the dual CFT. We then consider the radial derivative, ∂ r J i bulk , use the

JHEP03(2016)148
gauge equation of motion on the right hand side and then integrate over r from the horizon to the asymptotic boundary. After noticing that J i bulk vanishes at the horizon we deduce that the spatial part of current density of the dual field theory is a magnetisation current: where the local thermal magnetisation density, M ij , is given by In fact this procedure only defines M ij up to a constant times ij . We have fixed this constant in (A.7) by demanding that the zero mode gives rise to the magnetisation density m, given in (5.11), and defined by the first law (5.10). We observe that (A.6) impliesJ i =0. Thus, in the expansion (A.1), by considering (A.3), we must havec tx =c ty = 0 for all black hole solutions, and we have checked this in our numerics. A similar story unfolds for the heat current. We first define where ϕ ≡ i k A and we write i k F ≡ ψ + dθ for a globally defined function θ. We will take θ = −A t , so that ϕ = −θ = A t and ψ ν = ∂ t A ν . We find We now define Q a bulk ≡ √ −gG ar , (A. 10) and observe that when evaluated at the AdS boundary Q a bulk is the heat current density, √ γ(T ta − µJ a ), of the dual CFT. A calculation similar to above implies where the local thermal magnetisation density, M ij T (x), is given by In particular, sinceJ i = 0, we can conclude thatT ti = 0 also. Hence, in the expansion (A.1), by considering (A.3), we deduce thatc tx =c ty = 0 for all black hole solutions, and we have checked this in our numerics.

B Numerical convergence
In order to solve the system of PDEs with the boundary conditions explained in section 5 we discretise the computational domain [0, 1] × [0, 1] × [0, 1] on N r , N x and N y points respectively. As we will describe in the next paragraph, after approximating the partial derivatives on these points by an interpolation method, we are solving the resulting algebraic system of equations by using Newton's method.
Since our boundary conditions are periodic in the x and y directions, we use a Fourier spectral method in order to approximate the partial derivatives in these directions. This amounts to placing our points at equal distances 1/N x and 1/N y in those directions. For the discretisation in the "radial" coordinate r we have tested both Chebyshev spectral method and finite differences.
In figure 6 we present two convergence plots for the maximum value of |ξ| 2 max in our computational grid for a representative lattice (in this case a checkerboard lattice). We have used N x = N y = 60 points in the periodic directions for both cases of spectral and finite difference grids in the radial directions. For both both cases we have found power law convergence, |ξ| 2 max ∼ N −12.6 for spectral and |ξ| 2 max ∼ N −8.01 for fourth order finite differences. We can understand [33] the lack of spectral convergence (i.e. exponential convergence) from the non-analytic field expansion close to boundary (see (A.1)). We have also checked that ψ, defined in (5.5), satisfisfies ψ < 10 −10 for all of the solutions we have constructed.
We have implemented our numerical method in C++. The facility of class templates has been particularly helpful to accommodate the various numerical precisions we have used at low temperatures and in the convergence tests. At certain key points of our code we have specialised our templates to double, long double, and __float128. In order to generate and manipulate our data resulting from the different points of our domain to different computing nodes, we have followed a hybrid approach to parallelisation using a combination of MPI and OpenMP. After obtaining the matrix of variations in Newton's method implementation, we used the library PETSc [34] in order to solve the resulting linear system of equations.

JHEP03(2016)148
In order to achieve that, we used a block-ILU(0) preconditioner in combination with a GMRES iterative method. The typical system we had to solve was a ∼ 10 6 × 10 6 system with ∼ 10 10 non-zero matrix elements for the matrix we needed to invert.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.