Holographic magnetisation density waves

We numerically construct asymptotically AdS black brane solutions of D = 4 Einstein theory coupled to a scalar and two U(1) gauge fields. The solutions are holographically dual to d = 3 CFTs in a constant external magnetic field along one of the U(1)’s. Below a critical temperature the system’s magnetisation density becomes inhomogeneous, leading to spontaneous formation of current density waves. We find that the transition can be of second order and that the solutions which minimise the free energy locally in the parameter space of solutions have averaged stressed tensor of a perfect fluid.


Introduction
Phases that spontaneously break translational invariance are collectively referred to as Spatially Modulated (SM) phases and are widespread in Nature [1]. In condensed matter systems [2] in particular, they manifest themselves in various forms with spin density waves (SDW), and charge density waves (CDW) being the most common ones. Apart from the richness of this class of phases, its potential correlation with the physics of the pseudogap region of the high-T c superconductors' phase diagram [3] and the conjecture that the QCD phase diagram at finite temperature and intermediate density is dominated by a chiraldensity wave state [4] have triggered a lot of interest in understanding these phases.
Holography provides a natural theoretical framework to explore the physics of these phases at strong coupling. They are associated to black hole solutions that asymptote to AdS and have spatially modulated horizons. Many examples of modulated holographic phases have been discussed in the literature, starting with [5] and further explored in [6][7][8][9][10][11] (for earlier related work see [12][13][14]). It should be pointed out that constructing inhomogeneous black objects is not a new concept, dating back to the non-uniform strings. Gregory and Laflamme showed that the uniform black string, described by the product of a Schwarzschild solution and a circle, becomes unstable to modes with wavenumber smaller than a critical value. Since the critical deformation mode is static, they pointed out a new family of black strings without translation invariance: the non-uniform strings [15][16][17]. This new branch of solutions was perturbatively constructed [18] and non-perturbatively in [19]. However, it was shown that the non-uniform solutions have higher mass than the uniform strings and thus they are not thermodynamically preferred in a sufficiently low number of dimensions [20].

JHEP10(2016)038
The electrically charged AdS-RN black branes have been shown to suffer from spatially modulated instabilities. Such instabilities have been investigated for a class of D = 5 gravity theories with a single gauge-field and a Chern-Simons coupling in [11,13,14] and also for D = 4 Einstein-Maxwell theory coupled to a neutral pseudo-scalar field φ and possibly an additional gauge field in [21]. The non-linear striped black branes, connected to the zero modes of [21], have been constructed and studied in a series of papers [6][7][8][9]. It was shown that, depending on the non-linearities of the model, a second order phase transitions occurs at some critical temperature T c . Inhomogeneous solutions that break translation invariance in two directions were studied in [22], where checkerboard black holes with rectangular fundamental domain were constructed. In the same work [22], it was also shown that the striped black holes are continuously connected to the checkerboard black holes via rectangular lattice black holes. In [23], the electric stripes phase was shown to survive the existence of non-vanishing external magnetic, B. The full threedimensional family of solutions corresponding to oblique lattices was discussed in [24] and it was concluded that the free energy is minimised by a triangular lattice for B greater than a critical value B c . It was argued that this could be associated with minimal packing of circles in the plane.
Spatially modulated phases of CFTs placed in magnetic field, in the absence of a chemical potential, have been investigated in [25]. A big class of D dimensional bulk theories coupled to a scalar field, φ and to one or two U(1) gauge fields was studied. It was shown that, for particular choices of Lagrangian parameters, the dual field theory can admit phases that spontaneously break translation invariance via current density waves. In these models, the key feature of the instability is that the relevant mode preserves the internal U(1) symmetries. The zero modes that appear at the onset of the instability have been constructed in [25] and as we later explain, they modulate the magnetisation density of the dual field theory. Furthermore, it was shown [26] that these instabilities exist around magnetic branes solutions of N = 8 gauged supergravity in both D = 4 and D = 5. A different class of instabilities has been considered in e.g. [26][27][28] involving charged degrees of freedom, reminiscent of the instabilities discussed in [29]. In these case, the internal U(1) under which the unstable fields are charged is spontaneously broken and the new phase will necessarily break translation.
Here, we construct the backreacted geometries dual to the magnetisation density waves found in [25]. This involves the numerical solution of a set of coupled partial differential equations in two coordinates, employing the DeTurck trick for dynamical gauge fixing [30] which was first used in a holographic set-up in [31]. For the specific model we examined, we find a branch of spatially modulated solutions that extends to lower temperatures and dominates the thermodynamic ensemble within our assumption of breaking translations in only one direction. For this branch, we construct the complete two-dimensional space of solutions, specified essentially by the temperature T and the periodicity scale k. By minimising the free energy density, we determine the one dimensional branch of preferred solutions k = k(T ). We show that the latter have the stress-energy tensor of a perfect fluid. Down to the temperatures we explored, this branch seems to flow to a spatially modulated ground state. It is interesting to mention that down to the temperatures for which we have JHEP10(2016)038 constructed black holes, the entropy of the new phases remains parametrically large. The existence of spatially modulated ground states with finite entropy density is a possibility which we certainly plan to further explore in the future.
From the field theory point of view, we will examine a particle-hole symmetric medium which is charged under two different U(1)'s. The medium is magnetised with respect to one of them, under the influence of a uniform magnetic field. At temperatures below a critical value T c , we will see that the magnetic densities become modulated giving rise to spontaneous magnetisation current densities. While the connection with holography is not direct, the magnetisation density is a natural variable to discuss SDWs in the context of hydrodynamics [32]. One obvious difference with real systems exhibiting SDWs is the fact that the system we are considering is diamagnetic in its normal phase. Paramagnetic phases of holographic matter are certainly possible [23] and we expect that in a bottom-up approach similar phases can be constructed.
The remaining of this paper is organised as follows. In section 2 we introduce the model of interest and review the spatially modulated instabilities of the magnetically charged AdS-RN black holes of [25]. In section 3 we describe in detail the numerical method used to construct the back-reacted solutions and we report on the findings of this computation. Finally, we conclude in section 4. We have also included two appendices were we discuss in some detail the asymptotic expansion of our solutions and the numerical convergence of the method we used to construct them.

The setup
We consider a four dimensional theory of gravity coupled to a scalar field, φ, and two gauge fields, A and B gauging the corresponding symmetries U(1) A and U(1) B in the bulk. The bulk dynamics will be described by the Lagrangian where F = dA and G = dB. The above form is natural from the point of view of N = 2 supergravity in D = 4 [33,34] which commonly appears in the context of SUSY consistent truncations of D = 11 SUGRA (see e.g. [35][36][37][38]). Note that the wedge product terms between the field strengths that appear in top-down models are not going to be an important structural difference since they would be inactive for the solutions we consider in this paper. The equations of motion deriving from (2.1) are given by

JHEP10(2016)038
The functions V, Z A , Z B and W are chosen to admit the following expansion around φ = 0, where m s , n and s are constants. The equations of motion (2.2) admit AdS 4 of radius R 2 AdS 2 = 1/2 as a solution and will serve as the asymptotics of all the black hole space times we are going to construct.
According to the scenario of the introduction, we now consider the deformation of the boundary theory by a magnetic field β in the U(1) A . At high temperatures the bulk geometry is going to be described by the magnetic AdS-RN black hole, and φ, B trivial. This is a solution allowed by the choice (2.3) which lets us set the scalar φ and the second gauge field B consistently to zero. We have also introduced a length-scale L for later convenience when we fix the period of x. In this coordinate system the conformal boundary is at z = 0 with where we defined ε = z z h . The horizon of the black hole (2.4) is at z = 1 and its Hawking temperature is T = 12−β 2 z 4 h 8π z h . The near horizon limit of the extremal solution with T = 0 reduces to where ds 2 (AdS 2 ) is the metric on a unit radius AdS 2 space and we have scaled the spatial field theory coordinates according to e.g.x = L (β/12) 1/2 x. We now turn our attention to the instabilities of the magnetic AdS-RN black hole (2.4) towards phases with broken translations that were discussed in [25]. These can be understood as near horizon instabilities by examining the perturbation around the background solution (2.6). The functions that appear in (2.7) depend on the coordinates of AdS 2 and k is a constant. The linearised equations of motion (2.2)

JHEP10(2016)038
around (2.6) take the form where v = (δφ, δB) and the Laplacian is with respect to the AdS 2 metric. The mass matrix is given by It can be easily checked that depending on the choices of n and s, the eigenvalues of the above matrix can violate the BF bound signalling an instability. Interestingly, the lightest mode appears generically at a finite value of k.
The finite temperature zero modes related to the above instabilities have been constructed in [25] for m 2 s = −4 and various values of n, s. These were constructed by examining perturbations around the geometry of the black holes (2.4). The operator O φ , dual to the scalar field φ, was chosen to have scaling dimension ∆ = 1 and it was shown that the zero modes appear along a curve T (k), specified by the parameters n and s. At each point on that curve, a new branch of broken phase black hole solutions is expected appear for the corresponding fixed value of k. The maximum of this curve reveals the critical temperature in the case of a continuous transition. In figure 1 we plot T (k) for the case (n, s) = (−1, 2), which corresponds to (T c , k c ) = β 1/2 (0.08, 1.26).
The appearance of the bulk gauge field B in the mode (2.7) and the fact the instability shows up at a finite value of k, suggest that the broken phase will develop an inhomogeneous current density. As we explain in section 3.2, these currents can be seen as magnetisation currents due to the modulation of the field theory magnetisation density. At the same time, the operator dual to the bulk scalar φ will take a modulated VEV with the remaining fields backreacting at second order in perturbation theory close to the critical temperature T c .
In this note, we go beyond perturbation theory constructing the backreacted geometries by solving the non-linear equations of motion (2.2). As expected, the full functional dependence of V, Z A , Z B and W on φ is going to be relevant. The model we are going to consider has (2.10) The choice (2.10) has m 2 s = −4, n = −1 and s = 2. The zero modes appearing in this model correspond to the bell curve shown in figure 1. Furthermore, note that the model (2.10) exhibits a Z 2 symmetry with B → −B and φ → −φ simultaneously. This discrete symmetry becomes important at the non-linear level, relating the two branches of solutions expected to emerge at the critical temperature T c . If the Z 2 is not present, the two branches will be distinct and they will therefore have different thermodynamic properties.

The broken phase black holes
In this section we discuss the construction and properties of the non-linear solutions corresponding to the new branches of black holes proposed to exist in [25] and reviewed in section 2. In section 3.1 we describe the boundary value problem relevant to the present physical situation and the numerical methods we used to solve it. In section 3.2 we discuss the thermodynamics as well as the local magnetisation properties of the field theory states dual to the newly constructed black holes. In section 3.3 we discuss the numerical solutions we constructed along with some of their properties.

The ansatz and the method used
We consider the following ansatz for the back-reacted solutions where Q tt , Q zz , Q xx , Q yy , Q zx , a y , b y and h are functions of the radial coordinate, z, and x. Furthermore, we require that these functions obey periodic boundary conditions in the x direction, with period given by L. The function f (z) is the same with the one in equation (2.4). Note that this ansatz is generic enough to capture both the normal phase so- as well as the static zero modes considered at linearised level in [25]. In terms of the functions appearing in the non-linear ansatz (3.1), this mode takes the form

JHEP10(2016)038
The PDEs obtained when the equations of motions are evaluated on this ansatz are weakly elliptic, meaning that they are elliptic only for the physical degrees of freedom, and thus are unsuitable for numerics without gauge fixing. In this work, we employ the DeTurck method to resolve this issue [30,31]. According to this method, instead of solving the Einstein equations, one solves the Einstein-DeTurck equations which are obtained from (2.2) after making the following shift Hereḡ denotes a reference metric (which is required to have the same asymptotic behaviour as g) andΓ is the Christoffel connection ofḡ; we choose this to be the AdS-RN black hole metric (2.4). The resulting PDEs are then strictly elliptic and, with appropriate boundary conditions, can be solved numerically using a relaxation method.
In more detail, one discretises the coordinates of the PDEs to form a lattice. Our coordinates span z [0, 1] and x [0, 1). To approximate the derivatives of our functions at each grid point, we use a Fourier expansion in the periodic direction x. The treatment of the z direction is a more delicate one since our functions are non-analytic close to the boundary of AdS 4 at z = 0 as we explain in appendix A. We have used both spectral methods on a Chebyshev collocation grid as well as a fourth order finite difference scheme to cross check our results finding essentially the same outcomes within our numerical precision. We have included a convergence test in appendix B showing power law convergence for the spectral method in the z direction with a power compatible with the non-analytic terms in our expansion. A brief discussion on why spectral methods are preferred in this calculation over fintite differences and a quantitative comparison of the results of the two methods can be found in appendix B.
After fixing a discretisation scheme, the problem then reduces to solving a set of nonlinear 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, which presumably does not solve the PDEs. The solution is then iteratively improved using Newton's method in order to obtain functions that solve the PDEs to a better and better approximation. Provided a well posed boundary value problem, this procedure leads to a countable set of (locally unique) solutions to the modified equations of motion. As we will later discuss, in the context of spontaneously broken translations, this involves the fixing a Goldstone mode. A question that naturally arises is whether these solutions are also solutions of the initial equations (2.2). This is true when ξ 2 = 0, 1 i.e. when solutions corresponding to Ricci solitons are discarded. For this reason, after generating the solutions, we have performed convergence tests to verify that ξ 2 smoothly converges to zero within numerical precision at all grid points. This point is further discussed in appendix B.

JHEP10(2016)038
As we explain in more detail in appendix A, a set boundary conditions at z = 0 which are appropriate for AdS 4 asymptotics are On the other end of our computational domain, located at z = 1 we need to impose boundary conditions which guarantee a smooth Killing horizon of temperature T = 12−β 2 z 4 h 8π z h . This boils down to demanding that the functions, F (z, x) that parametrise our ansatz (3.1), admit an analytic expansion of the form The equations of motion impose constraints on the coefficients of the power series (3.5).
In order for the Euclidean signature metric to have a smooth fixed point at z = 1 we must have Q zz (1, x) = Q tt (1, x). By expanding the equations of motion at z = 1 we find another seven relations describing Robin boundary conditions imposed on that surface. Throughout the calculation, the magnetic field will be set to β = 1 while z h will be tuning the temperature.
A final point we need to address is the Goldstone mode associated with the spontaneous nature of the way we break translations. More concretely, if F (z, x) is a solution, then F (z, x + c) is also a solution of the boundary value problem for any constant c. Because of this, we should impose a condition in order to mod out these solutions from the moduli space. It is enough to impose this condition at one point and for the specific solutions we constructed, we imposed b y (1, 0) = 0. In fact, periodicity in the x coordinate guarantees the existence of a discrete number of points at which the partial derivative with respect to x of our functions vanish at any fixed z > 0. The model we are considering has a Z 2 symmetry of the type discussed at the end of section 2. In this case, setting to zero the even Fourier modes (including the zero mode) of b y and h is consistent. In the absence of these modes for b y , we can simply choose to set b y (1, 0) = 0.

Thermodynamics and magnetisation currents
To analyse the thermodynamics of our black hole solutions, we need to regularise our bulk action by adding appropriate surface terms [39,40]. For our bulk theory (2.1), this can be achieved by adding a surface term S ∂ to the bulk action which includes both the Gibbons-Hawking term as well as the necessary counterterms that render the total action finite. More specifically, 2 The specific choice of scalar field counterterms is certainly compatible with the ∆ = 1 choice for the boundary operator. However, when supersymmetry is involved a more careful treatment is required as pointed out in [41].

JHEP10(2016)038
Here K = g µν ∇ µ n ν is the trace of the extrinsic curvature of a z = const surface close to the boundary, with n µ the dual of the outward pointing normal unit vector and g ∞ being the determinant of the projection of the bulk metric on that surface. The ellipsis refers to terms which will not be relevant for the ansatz and boundary conditions that we are considering. Following [39,40], we now compute the expectation value of the boundary stress-energy tensor. The relevant terms are given by where we have included an appropriate power of z h since we are interested in the CFT one-point functions with respect to the metric In order to express this in terms of our asympotic data, we need to plug in the expansion of our fields close to the z = 0 boundary. Using the expansion (A.4) given in appendix A, we obtain and as we discuss in appendix A, the equations of motion imply that T µν is traceless. Defining the unit-norm timelike vector u = 2 −1/2 ∂ t on the boundary, we can express the free energyw, massm and entropyS averaged densities in terms of our numerical datā We now turn our attention to the electric currents of the dual field theory. Varying the on-shell action with respect to the asymptotic value of the gauge fields we obtain their VEVs (3.11) for which the expansion (A.4) gives

JHEP10(2016)038
for the non-zero components. As a comment, the only non-trivial component of the stress tensor continuity equation is satisfied provided that ξ µ = 0 on the background. This consists another check we performed on our numerics. For a fixed period L, the first law of thermodynamics gives whereM A is the thermodynamic magnetisation corresponding to U(1) A . For the backgrounds we construct in this paper, the magnetisation densities of both U(1)'s are inhomogeneous. Our black holes are bulk duals of equilibrium states and one would expect that the currents of the boundary theory should be a total derivative of a periodic antisymmetric rank two magnetisation tensor resulting in zero net transport of charge. To show this from the bulk [42], we integrate the gauge field equations of motion (2.2) over the radial coordinate z. Using regularity on the horizon and the definition (3.11) we can express Note that the definition of the magnetisation through the current density leaves us with a constant unfixed. We fix this constant by demanding that the thermodynamic magnetisation, defined through (3.14), is equal to the average of the local magnetisation. This is achieved by the choice we make when we write equation (3.15).
Following the arguments of [43], we find that the variation of the average free energy density with respect to the period L is given by Even though the solutions we construct numerically break translations in the x direction only, an identical statement is true regarding the variation of the free energy with respect to the period of the y direction. However, since our solutions don't depend on y, for all the solutions we construct herew From the above we conclude that for the black holes holes which locally minimise their free energy as a function of L, we havē and T xy = 0 by construction. For these solutions, we see that the averaged stress tensor corresponds to a perfect fluid.

JHEP10(2016)038
Scaling invariance of the dual CFT constrains further the various thermodynamic quantities we have discussed. In particular, the free energy has to be expressible according tow = β 3/2 f (T / √ β, L √ β) with f a dimensionless function. This symmetry therefore constrains the thermodynamic magnetisation according to Moreover, equations (3.10) and (3.19) imply thatT tt = 4P in agreement with the stronger requirement on the stress tensor being traceless locally. For the normal phase black hole solution (2.4) we can simply evaluate the thermodynamic magnetisation and susceptibility [23] in terms of T and β. Using (3.15) we find where the minus sign in the susceptibility points out that the normal phase describes a diamagnet with respect to U(1) A . It is easy to see from the leading mode (3.2), that close to T c the broken phase black holes will have a modulated U(1) B magnetisation density according to with higher order corrections away from T c . The effects of inhomogeneity on the magntisation of U(1) A will be second order. We show numerical evidence for this behaviour in the next section where we discuss the backreacted solutions.

Numerical solutions
In this section we present the solutions constructed for the model fixed by (2.10). Due to its Z 2 symmetry we expect to find only one distinct branch of solutions. From the numerics we see that this branch extends to lower temperatures, T < T c , covering the two-dimensional region below the bell-curve T (k) shown in figure 1. We present the profiles of some of the functions in figure 2, for a representative solution with T = 0.5 T c and k c = k.
In figure 3, we show a 3D plot of the free energy difference between the normal and the modulated phase, δw, as a function of T and k and we see that through out the moduli, the modulated solutions dominate. To specify the thermodynamically preferred branch branch of black holes one needs to minimise the free energy w with respect to k for fixed T . Doing this numerically yields the red line in figure 3. Various properties of the preferred branch are displayed in figure 4: panel (a) shows the value of the wavenumber k along the preferred branch as a function of the temperature T , and panel (b) and (c) show the free energy and entropy of the preferred solutions. We see that at zero temperature, we approach a ground state with a finite wavenumber k. As we discussed in the introduction, the entropy density remains parametrically large down to temperatures T /β 1/2 ≈ 2 · 10 −2 . Moreover, in figure 5 we display the the free energy and the averaged stress tensor for T = 0.5 T c as JHEP10(2016)038  functions of k. We see that for the preferred branch, at the minimum of the red curve, the green and the blue curves intersect showing that the averaged stress tensor reduces to the one for a isotropic ideal fluid. Moreover, the minimum averaged entropy given by the pink curve is slightly to the right of the preferred k at that temperature. This suggests that the free energy decreaaes faster at that value of k with increasing temperature. This is certainly compatible with the fact that the thermodynamically preferred k changes with T .
In figure 6 we plot the magnetisation densities M A and M B for the broken phase solution corresponding to (T, k) = (0.5 T c , k c ). As we discussed in section 3.1 the leading mode of the instability concerns the magnetisation of U(1) B and figure 6 agrees with expectation that the modulation of M B is a leading order effect. We can use equation (3.20) to find M A at the same temperature, in the normal phase. For T = 0.5 T c = 0.04 β 1/2 we find that M A ≈ −1.787 and a comparison with the top left of figure 6 demonstrates that the modulation in M A is a higher order effect when compared to that of M B .  which certainly conforms with our findings shown in the bottom plot of figure 6. In particular, we find that due to the Z 2 symmetry of our theory, the only non-trivial modes switched on are for n = 1, 3, 5, . . ..

Discussion
In this paper we have constructed numerically a two-parameter family of inhomogeneous AdS black brane solutions by backreacting on the instabilities of [25], for a particular model. These geometries describe a spatially modulated phase of the dual field theory, held at finite temperature and external magnetic field. Below a critical temperature T c the modulated phase becomes thermodynamically preferred over the unbroken one. We have shown substantial numerical evidence that the transition is second order. By exploring the thermodynamics of the entire parameter space of these solutions we find that the preferred solutions has a temperature-dependent periodicity for the modulation, in accordance with [5,9,10]. In contrast with previously constructed examples, we find that the periodicity is a monotonically decreasing function of temperature. Analysing our low temperature solutions, we find indications that the system approaches an inhomogeneous ground state with non-zero entropy. However, this picture could change as soon JHEP10(2016)038 as we achieve even lower temperatures. It would be interesting to try to construct the corresponding ground states directly by considering T = 0.
An obvious question to be addressed in the future is to relax the assumption of translational invariance in the one of the directions and construct the corresponding threeparameter family of solutions. The present work has not excluded the possibility that these solutions could dominate the ensemble. In the context of inhomogeneous phases of field theories at finite temperature and chemical potential, these solutions have been constructed in [22,24] and a competition between phases with modulation in one and two directions was found with the triangular configurations being preferred [24].
A more ambitious direction we plan to pursuit in the future is the fate of the instabilities of top-down models discussed in [26]. We expect a vast number of ground state geometries in these models with the finite temperature phase diagram exhibiting competing orders. 3 Among these ground states, of particular interest are the supersymmetric ones [26,46] and in particular the ones which will be modulated [26].

JHEP10(2016)038
equations then takes the form M · v = 0 where M is an 8 × 8 matrix that depends on δ. Demanding that non-trivial values of v exist implies that det M = 0 and this specifies the possible values of δ. The solutions come in 8 pairs. Apart from the usual modes, corresponding to the scalar, the 2 gauge fields and the metric respectively, we also obtain two modes with scaling dimensions that appear in the variations of the metric functions. These modes appear only in the Einstein-De Truck equations and are manifestations of the dynamical gauge fixing procedure. We then proceed in constructing the actual expansion where the functions Q (4) µµ (x) are fixed in terms of φ 1 , j a and j b . The expansion (A.4) bares a lot of similarities with the expansions that appeared in [24,47] in related topics. As we said above, we choose the operator dual to the scalar field to have scaling dimension ∆ = 1, which gives φ 1 the interpretation of the vacuum expectation value and φ 2 is the source. We want the breaking of translation invariance to be spontaneous and thus, we require the sources of all the operator, except of the background magnetic field, to vanish. For this reason, we impose φ 2 = µ a = µ b = 0. The expansion (A.4), is then specified in terms of nine coefficients {Q yy , Q (2) zx , φ 1 (x), j a (x), j b (x), g 1 (x), g 2 (x)} that will be fixed by solving the PDEs subject to the constraint which reflects the fact that the energy-momentum tensor of the dual field theory is traceless, T µ µ = 0. Demanding that ξ 2 = 0 poses extra constraints, namely that g 2 = −1 2 (3 + √ 33) g 1 and the Ward identity given by Q we will construct in the remaining of this paper the coefficients g 1 , g 2 are non-trivial, but they can be removed by a gauge transformation. However, even though they are pure gauge, these terms do leave an imprint on our numerics as, together with the logarithmic terms, they compromise the convergence rate.

B Numerical tests
In this appendix we discuss two tests we performed in order to check the quality of our numerics. The first one is to check the first law of thermodynamics, while the second one involves the convergence of ξ 2 .
From the differential form of the first law (3.14) we conclude that along a fixed k branch of solutions the quantity W =S +∂ Tw must vanish. In particular, using the values of w for equally spaced values of temperature T to construct an interpolating polynomial of degree 4 and evaluate ∂ Tw . Indeed, in figure 7(a) we see that, for model (2.10), this condition is met along the set of solutions with k = k c , where the temperature step between successive points used for the interpolation was ∆T ∼ 10 −3 β 1/2 . The change in the behaviour of W at T = 0.05 β 1/2 is justified by the fact the we changed our grid resolution at that point from N z = 40 to N z = 60 as we lowered the temperature.
One should also illustrate numerical convergence as the number of grid points is varied. For the DeTurck quantity, ξ, defined in section 3.1, we compute its maximum value on the grid, denoted by |ξ| 2 max . In figure 7(b) we present the convergence of |ξ| 2 max with the number of grid points in the z coordinate, N z , for a representative solution with temperature T = 0.5 T c and momenta k = k c in the model (2.10). The convergence to zero is power law, N −8.4 z , and not exponential as one would naively expect for spectral methods. This is due to the appearance of non-analytic terms with leading power z 8.7 in the UV boundary expansion of ξ 2 .
It is worth commenting on how the fourth order finite differences scheme compares to the spectral methods for the discretisation of the radial direction. The finite differences JHEP10(2016)038 solution at T = 0.5 T c and k = k c with N z = 150 and N x = 50 required 60GB of RAM to be constructed using our sparse solver. The corresponding solution based on spectral methods had N z = 90 and N x = 50 required 20GB of RAM. The finite differences solution yielded |ξ| 2 max ≈ 10 −11 while the solution based on spectral methods had |ξ| 2 max ≈ 10 −20 . Therefore, at these relatively low resolutions, the spectral methods proved to be much more efficient in terms of memory requirement. At these resolutions, the difference of the free energy from the normal phase solution only differed by 0.3% yielding essentially the same result for all practical purposes.
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.