Holographic Checkerboards

We construct cohomogeneity-three, finite temperature stationary black brane solutions dual to a field theory exhibiting checkerboard order. The checkerboards form a backreacted part of the bulk solution, and are obtained numerically from the coupled Einstein-Maxwell-scalar PDE system. They arise spontaneously and without the inclusion of an explicit lattice. The phase exhibits both charge and global U(1)-current modulation, which are periodic in two spatial directions. The current circulates within each checkerboard plaquette. We explore the competition with striped phases, finding first-order checkerboard to stripe phase transitions. We also detail spatially modulated instabilities of asymptotically AdS black brane backgrounds with neutral scalar profiles, including those with an hyperscaling violating IR geometry at zero temperature.


Introduction
The use of holography to describe strongly coupled gauge theories, including models for metallic, insulating and superconducting phases has received a great deal of attention recently. A variety of metals exhibit phases in which translational invariance or the underlying lattice symmetries are spontaneously broken, see [1] for a review. One possible outcome is striped order where translations are only broken in one direction resulting in a periodic configuration. Another possibility is checkerboard order, where no continuous translational symmetry remains and the phase becomes periodic in two spatial directions.
The spontaneous breaking of spatial symmetries can be modelled using holography, using a variety of holographic models [2,3,4,5,6,7,8,9]. Charge density may be introduced by turning on a chemical potential which sources a U (1) gauge field in the bulk, and a much used gravitational model which provides these minimal ingredients is Einstein-Maxwell theory, admitting asymptotically AdS Reissner-Nordstrom black brane (RN) solutions describing the normal phase. In 4d the spatial symmetries of the R 2 brane directions of RN may be spontaneously broken.
A symmetry breaking scenario of this type occurs in a Einstein-Maxwell-(pseudo)scalar model, where we set 16πG = 1. In [2] a perturbative analysis about the RN solution revealed marginally unstable modes breaking translational, as well as P and T symmetries. 1 The perturbations contain charge density waves and, due to the parity violating term in (1.1), also include a spatially modulated current. These marginal modes indicate the existence of a corresponding branch of black brane solutions which break translations.
One possibility is that the branch is spatially modulated in only one direction, referred to as striped black branes, dual to the striped phases of the field theory. Solutions of this type have been constructed [11,12,13,14,15], and are obtained numerically, solving a set of nonlinear elliptic 2d PDEs resulting from the bulk Einstein-matter equations.
In [14] the space of striped solutions was constructed with the dominant stripe's wavelength dependent on the temperature. It was also shown numerically that the line of dominant stripe solutions obeyed a particular averaged thermodynamic relation, which was subsequently derived as the result of a first-law relation involving the periodicity of the solution [16].
Another possibility is that there are black brane solutions which break translations in both boundary spatial directions. At least, at the linear level the marginal modes may simply be rotated and superposed. It is also possible that the resulting inhomogeneous black brane solutions are thermodynamically dominant to the striped solutions. The purpose of this paper is to seek these cohomogeneity-three black brane solutions with no surviving continuous translational invariance, which we shall refer to as checkerboard phases. We explicitly construct backreacted checkerboard configurations and explore their competition with striped phases. We find normal-to-checkerboard and checkerboard-to-stripe first order phase transitions as the temperature is lowered. 2 Finally, returning to the instability analysis [2], the linear marginally unstable modes of the RN solution at wavenumber k involve a spatially modulated current for the global 1 Instabilities which do not break P and T can also exist within this general class of models, without the ϑ term. This has been shown by a near-horizon, AdS 2 × R 2 analysis at T = 0 in [10]. 2 There are further possibilities of course, such as a triangular lattice. We have not investigated such possibilities in this work, though they would be a interesting future direction. U (1), J y ∝ cos kx. However, at linear order they do not involve a modulation of the charge density, J t . A charge density wave does occur at higher order in perturbation theory and with half the period, i.e. J t −J t ∝ cos 2kx whereJ t denotes a homogeneous component. In this paper we extend these instabilities to include deformations by the operator dual to φ. This is achieved by setting a constant source in the near-boundary expansion for φ. This changes the nature of the instability; whilst it still breaks P and T , the charge density is modulated at leading order in a perturbative instability analysis and with wavenumber k.
The paper is arranged as follows. In section 2 we outline the specific model choice, discuss the numerical method, boundary conditions and associated renormalised onepoint functions of the dual theory. In section 3 we present a checkerboard solution for the model investigated in [12] and investigate the space of rectangular checkerboards at fixed temperature. In section 4 we discuss spatially modulated instabilities with O φ deformations. In section 5 we present our main results for the deformed model -the existence of a dominant line of checkerboard solutions, reached by first order transition from the normal phase, and in some cases a first order transition to the striped phase at lower temperatures. We conclude in section 6.

Setup
We make the following choices for the model (1.1), introducing the parameters (n, c 1 ) corresponding to those used in the linear instability analysis of [2]. In all cases φ is dual to an operator of dimension ∆ = 2 with a ∆ = 1 source, φ (1) which appears as the leading piece in a near boundary expansion: We shall consider checkerboards in the presence of homogeneous deformations of this type and, as we shall see, it facilitates the competition of checkerboard phases with striped phases for the models considered.
The case n = 36, c 1 = 6 √ 2 arises in a consistent truncation from M-theory on an arbitrary SE 7 space [17], and motivates the functional form of (2.1). This particular choice has been well-studied, utilised in the study of holographic superconductors [18] and striped instabilities [2], including competition with superconductors under the influence of a magnetic field [19]. In [12] this model was used to construct striped black branes.
For these parameters the critical temperature for the striped instability is particularly low, consequentially a small adjustment to the model was made, setting n = 36, c 1 = 9.9 to improve the critical temperature. In section 3 we construct checkerboard black branes using this model; we find that only striped phases dominate, though we have not performed an exhaustive analysis. In section 5 we study the case n = 0, c 1 = 9.9, where we find that the checkerboard phases dominate with first order phase transitions once deformed by turning on φ (1) .

Gauge fixing and numerical method
We will directly construct stationary solutions of the coupled Einstein-Maxwell-scalar equations. Without modification or gauge fixing this is not an elliptic pde system as we require for this boundary value problem. An elliptic system can be reached via a suitable modification of the equations using the DeTurck method [20,21,22]. For the Einstein equations the Ricci tensor is replaced by where we have introduced the DeTurck vector ξ M = g N P (Γ M N P −Γ M N P ) and whereΓ is the Christoffel symbol for a reference metricg ab . We demand that ξ vanish on a solution of the modified equations so that we also have a solution to the unmodified system. In addition we find it necessary to remove longitudinal gauge modes of the U (1) gauge field A. A sufficient modification for this purpose is, where we have defined the scalar ψ =g M N∇ M (A N −Ã N ) with reference gauge fieldÃ M . With this additional term the principle part of gauge field equation no longer vanishes for longitudinal/gauge modes. 3 In addition to ξ M , we now require that the scalar ψ vanishes on any solutions.
The resulting set of Einstein-Maxwell-scalar equations with the additional terms are then solved numerically using a Newton-Raphson method. We use a grid of size N 3 = N z × N x × N y , represented spectrally with Fourier modes in x, y, the spatial boundary directions and Chebyschev polynomials in z, the radial holographic direction. Solutions are periodic in the x and y directions with wavenumbers k x and k y respectively. For lower grid sizes we find it convenient to take N z = 2N x = 2N y , primarily as it assists in the accurate extraction of the free energy from the numerical data. We will provide convergence checks of ξ and ψ and the free energy with N .
Finally, some comments on implementation. With the inclusion of currents resulting from P and T breaking through the ϑF ∧ F -term, the system studied has 15 equations of motion (for 15 fields: 10 metric, 4 gauge field and 1 scalar), as described in the next section, 2.2. We find it convenient to generate these only numerically, at runtime. Then, each step in the Newton-Raphson iterative method consists of two distinct computational stages: first constructing the Jacobian matrix of first derivatives of the equations (including boundary conditions) with respect to the fields. This may be done with a numerical finite differencing and may be efficiently parallelised. Because it is a 3d problem containing at most second derivatives, the resulting matrix is sparse despite being spectrally represented. The second Newton-Raphson stage is solving the resulting sparse linear system, for which we use a bi-conjugate gradient method. Good convergence occurs within only a handful of Newton-Raphson steps.

Ansatz and boundary conditions
To start we can consider a metric ansatz for a coordinate system adapted to the timelike Killing vector field T = ∂ t , where a = 1, 2, 3 with x a = (z, x, y) a . The horizon position dictated by g tt (x h ) = 0 where T is null. We do not restrict the base manifold with metric h to have any particular symmetries. Similarly, A a , will be non-zero in general. In the absence of gauge fixing this ansatz is then just a way of packaging the full quota of 10 metric functions, each of which is t-independent but depend on the three remaining x a coordinates. In practise we repackage the above ansatz for convenience: where i = 1, 2 labels boundary spatial directions, x i = (x, y) i . Each of G M N are functions of (z, x, y) and f (z) = 1 4 (4+µ 2 z 4 −z 3 (4+µ 2 )), with f (z h ) = 0 giving the horizon location, chosen to lie at z = z h = 1.
we recover the RN solution. The reference metricg and reference gauge fieldÃ are chosen to take on these values.

UV
In the UV near z = 0 for solutions where ξ = ψ = 0 we have the following expansions for the matter fields and for the metric functions, with additional relations constraining subleading data such that they satisfy U(1)-current conservation, stress tensor conservation and a conformal Ward identity, derived in section 2.3. From these expansions we can simply read off the Dirichlet boundary conditions required for the numerics. In the case of φ we work with φ c (z, x i ) ≡ φ(z, x i )/z and set a (constant) Dirichlet condition for it in the UV.

Horizon
At z = 1 we seek a regular, non-extremal horizon. This entails a near-z = 1 expansion for the matter fields, where the behaviour of A z is determined by the near horizon expansion of the ψ = 0 condition. For the metric functions, with a Dirichlet boundary condition G + zz = G + tt . Subleading orders in this expansion are determined by these leading data. We use a second Dirichlet condition for A t (z, x i ) = 0.
For the remaining fields we impose the leading form of the expansion. In particular, As a check of this boundary condition, of the near-horizon behaviour and as a more general check of the solutions, we confirm that the analytic relations between the coefficients of the O(1 − z) terms and leading order data hold numerically on the solutions constructed. Similarly the leading part of (2.17) holds with the above boundary condition, or alternatively may be fixed using a Dirichlet condition giving equivalent results.
Finally, the temperature is given by 22) and the entropy density by (2.23)

One-point functions
To perform renormalisation of the action (1.1) and extract the one-point functions we first transform the metric to Fefferman-Graham form near the AdS boundary. Details are provided in appendix A. Here we will be brief and just quote the key results. The renormalised action including the appropriate Gibbons-Hawking term is given by In total the one point functions are summarised by the variation, x, y) µ , γ is the induced metric on the boundary, and sub/superscripts in brackets label coefficients in the small-z near boundary expansion. For instance, φ (1) is as defined in (2.2). The one point functions are given by, In addition invariance under the Weyl transformations δγ (1) gives the conformal Ward identity, Inserting the expressions for the one-point functions above we verify this holds using the equations of motion. We can also use this identity as a check of the numerical extraction of the relevant stress tensor components.
An expression for the free energy density is given by, where T and s(x i ) are defined in (2.22) and (2.23). Throughout we will adopt thermodynamical quantities averaged over the spatial boundary directions of the computational domain, these will be denoted with a bar, for instance, Finally note we work in the grand canonical ensemble and appropriately construct dimensionless quantities using powers of the chemical potential, µ. This allows us to set µ = 1 where it streamlines presentation, specifically we shall do this where we present numerical results.
First consider the linear instability of RN, details of which can be found in [2,14]. This is formed from the following set of consistent perturbations with momentum k, δφ = λ(z) cos kx, δA y = a y (z) sin kx, δg ty = h ty (z) sin kx (3.1) subject to horizon regularity and normalisability at the boundary. This results in a 'bellcurve' of k-dependent critical temperatures. At higher orders the charge density becomes modulated with period 2k. In what follows we study the nonlinear branch of black brane solutions which emerges from a linear combination of two such modes, one modulated in the x-direction as above, and one modulated in the y-direction.
The solution is a checkerboard, periodic in x and y with momenta k x and k y respectively. The highest temperature at which the checkerboards connect with RN is the same as for the striped solutions and is given by T = T * = 0.0236 for a square config- The scalar vev O φ and charge density distribution is presented in Figure 1, together with integral curves of the (divergence free) boundary current, J i . The current is seen to circulate with the sense alternating between adjacent plaquettes of the checkerboard.
This should not be too surprising; this would be qualitatively the picture formed for the current near T * simply by superposing two of the linear modes (3.1). We emphasise however that this solution is not near the critical temperature, and is backreacted in the nonlinear regime. The convergence of ξ and ψ with the number of grid points for this solution shown in Appendix B, exhibiting exponentially fast convergence.

Varying k
At fixed T there is a 2-parameter family of solutions parameterised by k x , k y . Both checkerboards and stripes exist within this family and are continuously connected, as we shall demonstrate in this section. We wish to minimise the free energy in this family. Defining, where ω normal denotesω in the normal phase. ∆ω for striped and 'square' checkerboard solutions with k ≡ k x = k y at fixed T = 0.55T * and varying k are shown in Figure 2.
For this model we can see that the striped solutions are thermodynamically dominant, at least to the classes of checkerboards we have investigated in this section, where the dominant stripe has ∆ω −7.5 × 10 −4 at k 0.73 (more details can be found in [14]).
In the section 5 we will present a different set of model parameters where checkerboards are the dominant phase. Following [16] (see also [15]) we may compute the variation of theω with respect to the momenta. We find, and similarly for y. Note that if we look at the 'square' solutions where k x = k y we have by symmetry k x ∂ω ∂kx = k y ∂ω ∂ky . If we then use (3.3) in an iterative scheme such as Newton-Raphson in order to find the minimumω starting with a k x = k y checkerboard, we will stay within the k x = k y family. In practise we do not implement k x ∂ω ∂kx = 0 numerically e.g. as a boundary condition, but we do utilise the relation (3.3) both in order to check our results and to assist in the determination of the minimumω. We illustrate the agreement of the relation (3.3) with the solutions constructed for the case of k x = k y in Figure 2; the gradients computed in (3.3) shown as tangents to the data, match the gradients of the data itself.
We would also like to understand how the checkerboard solutions presented in Fig The normal phase may be constructed by numerically solving a set of ODEs. We seek solutions of the form, where f (z) is defined in section 2.2. This results in a system of second order ODEs for matter fields a, Φ with first order equations for the metric functions T, Z. The construction proceeds via a standard shooting problem to enforce horizon regularity and boundary normalisability, see for example [18,19]. Counting the number of undetermined coefficients in the near horizon expansion (there are 3 with the horizon position fixed at z = 1) and near boundary expansion (5) reveals solutions will exist in two-parameter families, as expected. We may take convenient parameters to be the source for O φ (this is φ 1 /µ for the m 2 = −2 case) and the temperature T /µ. The zero temperature solutions will be discussed for specific examples later in this section.
The spatially modulated instability of the RN solutions involving the fluctuations (3.1) continues to the φ = 0 branes. However, because the normal phase has Φ(z) = 0 (3.1) no longer forms a consistent set of perturbations. We find it convenient to work in a particular gauge with δg tt = 0 and δg zx = 0, which fixes a set of gauge modes arising from first-order coordinate transformations of the background solution. Here a consistent set of perturbations for any τ, V, ϑ consists of second order ODEs for (3.1) together with We may write the equations as second order for a t (z), and first order for h xx (z), h yy (z) and h zz (z).
Note that due to a t , the charge density may become modulated at leading order with a momentum k rather than at sub-leading order with momentum 2k as for the instabilities of the normal phase with φ = 0. This leading order modulation can be seen explicitly in the nonlinear checkerboard and stripe solutions constructed in section 5, with the charge density modulated with the same period as the scalar field φ, both near the critical temperature and also extending to solutions with significant backreaction.
The fluctuations comprise a coupled system of ODEs with total differential order 11.
At the horizon the fields take the form, where we have shown any undetermined coefficients in the near-horizon expansion. In the UV we enforce normalisability, in particular note that none of the boundary metric components, µ or the source φ (1) are affected by this mode, and remain constant as required for spontaneous modulation.
There are 11 undetermined coefficients at linear order, one of which may be fixed by linearity. Coupled together with the equations for the normal phase, the total differential

n = 0
When n = 0 it is known that the AdS 2 ×R 2 near-horizon geometry RN is linearly unstable to k = 0 modes involving the scalar [23,2]. Thus we may expect that turning on φ (1) can drive the IR away from the φ = 0 AdS 2 × R 2 to something else, such as a hyperscaling violating (HSV) geometry with a running scalar. Indeed, the theory admits the HSV solutions [24,25], This solution has a hyperscaling violation exponent θ = 4 and dynamical critical exponent z = −9. Based on this scaling property we can infer the low temperature scaling of the entropy density, s ∝ T 2/9 . We confirm this expectation by numerically constructing the finite temperature branch, for which the entropy density is shown in Figure 4. Hence we see that dialling the parameter φ (1) results in at least one quantum phase transition in the normal phase, moving from the φ = 0 AdS 2 × R 2 to the HSV geometry (4.18)- (4.20). It is interesting to note the effect that this change of IR has on the spatially modulated instabilities (4.3)-(4.5). Since the AdS 2 × R 2 case is k = 0 unstable [2] by continuity we expect that the instability survives for small φ (1) = 0, at least for finite temperature. Indeed, the critical temperature bell-curves do continue to φ (1) = 0 as shown in Figure 5. The change is most striking at lower temperatures where the linear instability at T = 0 seems to disappear entirely, at least for large enough φ (1) . It would be interesting to investigate this feature directly, and any connection with k = 0 stability properties of the HSV infrared geometry.

Checkerboard transitions
The checkerboards of section 3 are thermodynamically subdominant to the striped solutions. In this section we explore the role of introducing a homogeneous source φ (1) for the operator dual to the scalar φ. Linear instabilities in this case were constructed in section 4. Here we consider the case n = 0 with associated instabilities mapped out in section 4.1. The key results of this section are checkerboards which are thermodynamically dominant to the stripes, including first order normal-to-checkerboard first order transitions followed by checkerboard-to-stripe first order transitions at lower temperatures.
An example of a checkerboard in the φ (1) -deformed model is shown in Figure 6, at T = 0.57T c where T c 0.0707 is the critical temperature of the normal-to-checkerboard phase transition demonstrated shortly. We show convergence tests of the numerics for this solution, presented in Appendix B.
At each temperature we construct the thermodynamically dominant solution of the square checkerboard family by minimising ∆ω with respect to k x = k y , as described in section 3.1. Similarly we construct the dominant stripe solutions. The temperature dependence for the φ (1) = 0 case is shown in Figure 7. Clearly the striped solutions are thermodynamically dominant, as in the n = 36 case discussed in section 3.
However the crucial behaviour of this model lies in the φ (1) deformations. In Figure 8 we show the temperature dependence of ∆ω for the case φ (1) = 0.13. The introduction of φ (1) has opened a swallow-tail structure near the instability threshold, pushing the checkerboard transition temperature higher than that of the stripes. This results in a first order phase transition from the normal phase to the checkerboard phase, followed  Considering several other fixed-φ (1) slices, we find that the temperature of the checkerboard to stripe first-order phase transition decreases as φ (1) is increased. At sufficiently small, positive φ (1) there is no checkerboard transition at all. At higher values of φ (1) the checkerboard phase is dominant for all temperatures considered, although note we have no data at very low temperatures. Based on these observations we anticipate the qualitative structure of the phase diagram illustrated in Figure 9. Note that the tri-

Discussion
We have constructed cohomogeneity-three, finite temperature stationary black brane solutions dual to a CFT exhibiting checkerboard order. Due to the parity violating term (1.1) the phase breaks P and T resulting in J i = 0 which circulates within each checkerboard plaquette. The phase appears spontaneously, without any UV lattice or other inhomogeneities introduced by hand.
We found qualitatively different behaviour depending on the model parameters, (n, c 1 ) as well as the constant deformation governed by the source φ (1) . The checkerboards constructed are thermodynamically preferred over the normal phase below a critical temperature, with a first order transition depending on the model. We also explored competition with the striped phases. At n = 0, c 1 = 9.9 and φ (1) = 0 we found that there can be up to two first-order phase transitions, from the normal phase to the checkerboard phase and from the checkerboard phase to the striped phase at lower temperatures.
We constructed linear k = 0 instabilities of normal phase solutions with φ = 0. Here, the charge density becomes modulated at leading order in the instability at wavenumber k. We investigated one model where the IR becomes hyperscaling-violating as the deformation parameter φ (1) is increased. It would be interesting to extend this analysis and investigate T = 0 k = 0 stability directly in the absence of an AdS 2 × R 2 IR solution.
A remaining open question is the nature of the T = 0 ground states of both the striped and checkerboard phases. It would be particularly interesting to see the explicit T = 0 manifestation of the checkerboard to stripe phase transition in the phase diagram of Figure 9. Our numerical analysis has restricted to the case where the checkerboards are rectangular. There may be a larger space of solutions and it would be interesting to investigate the dominant shape of these phases. 4 Finally we comment on recent progress in the construction of explicitly modulated phases utilising spatially dependent field theory sources and bulk fields in such a way that only ODEs need to be solved in the bulk [27,28] 5 . Note that the approach taken in [28] allows for analytic construction of the background geometry 6 though note it may be generalised to other contexts by adding a dilaton, as explored in [34]. These are technically similar in spirit to earlier work in 5d where helical symmetry is exploited to reduce a problem to ODEs [3,35,36]. It is conceivable that the mechanisms of [27,28] may be used in the study of spontaneous modulation. 7 One concern in the present context is that, since they are restricted to special scalar configurations which allow the trick to work, they may not be capable of capturing the dominant phase.

Acknowledgements
I would like to thank Toby Wiseman for valuable discussions on numerical issues. I would also like to thank Aristomenis Donos and Jerome Gauntlett for early discussions.
I acknowledge the use of the IRIDIS High Performance Computing Facility and associated support services at the University of Southampton, in the completion of this work.

A Holographic renormalisation
First consider the components of the metric organised according to z-directions and all others, The near boundary expansion to the order of interest takes the following form, We now consider the computation of the one-point functions, performing the holographic renormalisation procedure in Fefferman-Graham coordinates near the boundary. In particular, near the boundary we expand in the following coordinate system, where r = 0 gives the boundary.
The relation between these two coordinate systems in the vicinity of the boundary may be expressed as, where the transformation is given by the following functions, The asymptotic form of the metric γ can now be written in terms of the components of the original metric using the transformation above, The near-boundary expansion of the scalar field becomes, Similarly, we may re-write the asymptotic expansions for the gauge field, where the cross terms appear too high in order z, r to affect the expressions. We may also perform a gauge transformation to eliminate the A r component near the boundary,

A.1 One point functions
In this section we use the coordinates (r,x), defined in (A.6),(A.7). First we define the unit one-form This is normal to the boundary at r = 0. We can then define the orthogonal projector, In particular, p rA = 0, whereas, similarly, K rA = 0 and, Expanding for small r, we find, The renormalised action is given by, with associated stress tensor, A standard analysis [37,38] reveals the counterterms, Here we focus on flat boundary metric and we have omitted curvature counterterms. we arrive at, from which we identify the expectation value of the field theory stress tensor, Similarly we may compute the one point function for the scalar, 30) and the current, In Figure 10 we present the convergence with the number of grid points for the solution presented in section 3. The values of |ξ| max and |ψ| max converge towards zero exponentially. The accuracy to which we can extract ω N (we use three numerical derivatives in this case) saturates, but at that point ω N is very small.
In Figure 11 we present the convergence with the number of grid points for the data of section 5 near the checkerboard-stripe first order transition, at a temperature of T = 0.04 and momenta k x = k y = 0.6. Comparing with the other model in Figure 10 the values of |ξ| max and |ψ| max are not quite as small. However, they still converge towards zero exponentially fast with N as required. Moreover, the free energy has converged sufficiently at the values of N shown. The situation improves at higher temperatures, e.g. near the swallowtail where the values once more mirror those of Figure 10.  Figure 11: Convergence with the number of grid points (for definitions see Figure 10) for the data of section 5 near the checkerboard-stripe first order transition, at a temperature of T = 0.04 and momenta k x = k y = 0.6.