Thermodynamics and transport of holographic nodal line semimetals

We study various thermodynamic and transport properties of a holographic model of a nodal line semimetal (NLSM) at finite temperature, including the quantum phase transition to a topologically trivial phase, with Dirac semimetal-like conductivity. At zero temperature, composite fermion spectral functions obtained from holography are known to exhibit multiple Fermi surfaces. Similarly, for the holographic NLSM we observe multiple nodal lines instead of just one. We show, however, that as the temperature is raised these nodal lines broaden and disappear into the continuum one by one, so there is a finite range of temperatures for which there is only a single nodal line visible in the spectrum. We compute several transport coefficients in the holographic NLSM as a function of temperature, namely the charge and thermal conductivities, and the shear viscosities. By adding a new non-linear coupling to the model we are able to control the low frequency limit of the electrical conductivity in the direction orthogonal to the plane of the nodal line, allowing us to better match the conductivity of real NLSMs. The boundary quantum field theory is anisotropic and therefore has explicitly broken Lorentz invariance, which leads to a stress tensor that is not symmetric. This has important consequences for the energy and momentum transport: the thermal conductivity at vanishing charge density is not simply fixed by a Ward identity, and there are a much larger number of independent shear viscosities than in a Lorentz-invariant system.


Introduction
Nodal line semimetals (NLSMs) are a recently discovered class of materials, in which two electronic bands intersect along a closed curve in momentum space at or near the Fermi energy. As reviewed below, this intersection is protected by the non-trivial topology of the electronic band structure, combined with the discrete symmetries of the system. Examples of materials with band structures containing nodal lines are Ca 3 P 2 [1,2], PbTaSe 2 [3], and ZrSiS [4,5]. As with other recently discovered semimetals, such as the Dirac and Weyl semimetals, these topological materials are of particular interest because of their potential for energy-efficient electronics, (pseudo)spintronics devices, and quantum information processing. This is not only because of the topologically protected and often dissipationless edge states that these materials possess, but also due to the sometimes anomalous bulk properties of these semimetals [6][7][8][9].
A simple model exhibiting a nodal line is a free fermion ψ of mass M coupled to an external antisymmetric tensor field b µν and having the Lagrangian density [10] L = iψ (γ µ ∂ µ − M + ib µν γ µν ) ψ, (1.1) where γ µν = 1 2 [γ µ , γ ν ], and throughout this paper we use units such that = c = 1. Greek indices such as µ and ν label all four spacetime coordinates, while Latin indices such as i and j will correspond only to the three spatial coordinates. Our convention for the Dirac matrices is {γ µ , γ ν } = 2η µν , where η µν is the Minkowski metric in mostly-plus signature, and we take (γ 0 ) † = −γ 0 and (γ i ) † = γ i .
The single-particle Hamiltonian following from equation (1.1) is where k = (k x , k y , k z ) is the momentum of the fermion. To obtain a nodal line we take the only non-zero component of the tensor field to be b xy = −b yx = b for some constant b, explicitly breaking the symmetries of rotations about the x and y axes, Lorentz boosts in the x and y directions, and time reversal, sinceψγ ij ψ → −ψγ ij ψ under time reversal. Without loss of generality, we will take b > 0. With this choice of b µν , the eigenvalues of the Hamiltonian (1.2) are with the two ± signs independent. For |M | < b, two of these eigenvalues meet at ε = 0 along the circle at k 2 x + k 2 y = b 2 − M 2 and k z = 0, as illustrated in figure 1a. This circle is the nodal line.
The degeneracy of eigenvalues at the nodal line is protected by the non-trivial topology of the Berry connection A, along with the discrete symmetries of the Hamiltonian. Recall that for a Hamiltonian depending on some set of parameters k, with corresponding normalised eigenstates |u( k) , the Berry connection is defined as A i = i u( k)|∂ k i |u( k) . If we take |u( k) to be either of the two eigenstates of the Hamiltonian (1.2) with eigenvalues that meet at the nodal line, then the integral of A about any closed contour C in momentum space that winds once around the nodal line is [10] C A = ±π, (1.4) depending on the relative orientation of C and the nodal line. An example contour C is depicted in figure 1b. If the integral of A around any curve C that does not encircle a band-touching vanishes, then a small perturbation to the Hamiltonian that does not break the symmetries cannot destroy the nodal line without also changing the right-hand side of equation (1.4) to zero. Such a large change would violate the assumption that the perturbation is small. The stability of the nodal line therefore rests on the assumption that the left-hand side of equation (1.4) vanishes for any C that does not encircle the nodal line. This is not always the case for a generic Hamiltonian, and typically requires the presence of multiple discrete symmetries. For the Hamiltonian in equation (1.2) with b xy = −b yx = b, the relevant symmetries are parity, rotation about the z axis by angle π, and the composition of time reversal with an x → −x reflection. 1 For further details, see ref. [10].
The non-trivial topology of the band structure becomes clearer if we convert the Hamiltonian (1.2) into a two-band Hamiltonian of the type commonly used in the condensed matter literature. This can be achieved by employing a basis in which the γ matrices are Kronecker products of Pauli matrices, for example γ 0 = iσ 3 ⊗σ 2 , γ 1 = σ 1 ⊗1 2 , γ 2 = σ 2 ⊗1 2 , and γ 3 = σ 3 ⊗ σ 3 [10]. In this basis, the Hamiltonian with b xy = −b yx = b is H( k) = (k y σ 1 − k x σ 2 + M σ 3 + b1 2 ) ⊗ σ 2 − k z 1 2 ⊗ σ 1 , (1.5) 1 The reflection inverts the sign of b, undoing the effect of the time-reversal transformation on the Lagrangian (1.1). and the energy eigenstates take the factorised form |ε = |ε 1 ⊗ |ε 2 , where |ε 1 is an eigenstate of (k y σ 1 − k x σ 2 + M σ 3 + b1 2 ) with eigenvalue ε 1 = b ± k 2 x + k 2 y + M 2 , while |ε 2 is an eigenstate of (ε 1 σ 2 − k z σ 1 ). Hence, if we choose the minus sign in ε 1 , then |ε 2 is an eigenstate of with eigenvalues equal to the two eigenvalues in equation (1.3) that have a minus sign inside the square root. These are the two eigenvalues that meet to form the nodal line, i.e. the inner two eigenvalues in figure 1a, so we may regard H 2 ( k) as an effective Hamiltonian for these two bands. This Hamiltonian has a 'chiral' symmetry as it only depends on two Pauli matrices, arising from the combination of discrete symmetries described above [10]. If we write the two-band Hamiltonian in the form of a Zeeman interaction, H 2 ( k) ≡ − B( k) · σ, the field B( k) therefore has non-zero components only in the x and y directions. If we vary k around a closed loop in momentum space encircling the nodal line, B( k) rotates about the origin of the (x, y) plane an integer number of times. It is this nontrivial winding number that is measured with the integral (1/π) C A. Note also that adding a perturbation proportional to σ 3 to the Hamiltonian will gap out the nodal line, so it is precisely the chiral symmetry that protects it. An important feature of the nodal line band structure is that the density of electron states g(ε) vanishes at the Fermi surface. This leads to weak screening of the Coulomb interaction between electrons, since the screening length diverges as g(ε) −1/2 , see for example ref. [11]. The electrons in an NLSM may therefore become strongly interacting, with significantly different physics from the free model in equation (1.1). Indeed, evidence has recently been found for strong correlations between electrons in the nodal line semimetal ZrSiSe [12]. One tool to model strongly interacting systems is the anti-de Sitter/conformal field theory (AdS/CFT) correspondence, also known as holography, which relates strongly interacting quantum field theories (QFTs) to weakly coupled gravitational theories [13][14][15]. A holographic model of an NLSM was constructed in refs. [16][17][18], by writing down a gravitational theory holographically dual to a strongly interacting QFT with the same operator content and symmetries as the Lagrangian in equation (1.1). 2 The holographic NLSM model includes two couplings, M and b, analogous to the couplings with the same names in the free model (1.1). At zero temperature, the model exhibits a quantum phase transition at a critical value of the dimensionless ratio M/b = (M/b) crit. . Below (M/b) crit. , spectral functions of composite fermionic operators exhibit many circular lines of sharp peaks in momentum space at zero frequency, which are identified as the nodal lines of the system [17]. Above (M/b) crit. no such nodal lines appear. The T = 0 phase structure of the holographic nodal line model therefore appears qualitatively similar to that of the free model, except for the presence of many nodal lines in the low M/b phase, rather than just one.
The many nodal lines in the holographic model pose a challenge to the application of holography to real NLSMs, which have only a finite number of nodal lines. However, we will show in section 3 that at non-zero temperature thermal broadening causes many of the peaks in the fermion spectral function to merge, such that there is a range of temperatures for which only a small number of nodal lines are visible. One might expect holography to provide a more realistic model of nodal line physics in this temperature range.
Another possible method to model a strongly interacting system with only a small number of nodal lines would be to use semi-holography [30][31][32][33], in which an elementary fermion undergoes strong interactions mediated by a holographic conformal field theory. A second advantage of semi-holography is that it yields fermion spectral functions ρ that when integrated over frequency ω obey the sum rule ∞ −∞ dω ρ = 1. In contrast, in the fully holographic approach one only has access to spectral functions of composite fermions, which do not obey this sum rule, and therefore cannot be directly compared to electron spectral functions in real materials, as measured in angle-resolved photoemission spectroscopy (ARPES) experiments.
In this paper we will take the fully holographic approach. We will analyse the thermodynamics and various transport properties of a version of the holographic NLSM model of ref. [16], modified by an additional coupling in the gravitational theory that will allow us to tune some of the infrared (IR) properties of the system. In section 2 we will write down the holographic NLSM model, describe its solutions, and present numerical results for its thermodynamics. As for the model of ref. [16], at zero temperature we find a quantum phase transition between topologically trivial and non-trivial phases. A key result, described in section 2.4, is that the explicitly broken Lorentz invariance of the NLSM implies that the stress tensor is not a symmetric tensor; it has an antisymmetric contribution proportional to the coupling b responsible for the violation of Lorentz invariance.
In section 3 we present numerical results for the spectral function of a composite Dirac fermion operator holographically dual to a pair of bulk probe fermions, in the topologically non-trivial phase. We find that at zero temperature the spectral function exhibits multiple nodal lines, appearing as sharp peaks at zero frequency and non-zero momentum. The nodal lines gradually disappear into the continuum as the temperature is increased.
In section 4 we give results for various transport coefficients of the holographic NLSM, namely the electrical and thermal conductivities, and the shear viscosity. At zero temperature in the topologically non-trivial phase, the DC electrical conductivity in the direction orthogonal to the plane of the nodal line is non-zero, while in other directions it vanishes. All components of the DC conductivity vanish in the topologically trivial phase. The asymmetry of the stress tensor induced by the Lorentz-violating coupling b is crucial to our analysis of energy and momentum transport. It implies that the thermal conductivity at vanishing charge density is not completely fixed by a Ward identity, as it is when b vanishes, and there are multiple shear viscosities, only some of which take the value η = s/4π universal to Lorentz-invariant holographic systems [34][35][36][37][38]].

Action and equations of motion
In this section we describe the holographic model that we will use, based on that of refs. [16][17][18]. The model is a bottom-up construction, built by writing down a gravitational action for a set of fields holographically dual to a set of operators with the same spins and Ward identities as the operators of the free model (1.1). Concretely, the field content on the gravitational side is: • The metric G mn , holographically dual to the operatorT µν . As we will discuss in detail in section 2.4,T µν is closely related to the stress tensor.
• A complex scalar field φ, dual to a complex scalar operator O. Roughly speaking, we think of the real part of O as modelling the fermion mass operatorψψ, while the imaginary part of O modelsψγ 5 ψ.
• A real two-form field B mn , dual to an antisymmetric tensor operator O µν , to model the operatorψγ µν ψ. 3 • A U (1) gauge field A m , dual to a conserved U (1) current, modelling the conserved current iψγ µ ψ in the dual field theory.
One might also choose to include an axial vector field A 5 m to holographically model the axial current iψγ µ γ 5 ψ. However, the axial current will play no role in any of the physics analysed in this paper, so will be neglected.
We take the gravitational action to be where G N is Newton's constant, L is the curvature radius of the asymptotically-AdS spacetime M, F = dA, H = dB, B 2 ≡ B mn B mn and similar for H and F , γ is the induced metric on the conformal boundary ∂M, K is the mean curvature of ∂M, and η, σ, and λ are dimensionless coupling constants. The counterterms S ct are boundary terms necessary to render the on-shell action finite and ensure a well-defined variational principle. The form 3 It has been argued that a more realistic model may use a complex, self-dual two-form to model both ψγ µν ψ andψγ µν γ 5 ψ [18], with the self-duality constraint representing the duality relationψγ µν γ 5 ψ ∝ µν ρσψγ ρσ ψ, where is the Levi-Civita symbol [39]. However, it is extremely unlikely that Bmn is genuinely dual to an operator that may be written asψγ µν ψ, where ψ is an elementary fermion in the dual QFT, and in any case the physics of the model does not appear to depend sensitively on the choice to impose self-duality [18], so for simplicity we will work with real Bmn with no self-duality constraint. When massless, the two-form field was suggested in ref. [40] to model superfluidity, with strings winding the Euclidean time circle in the black hole phase spontaneously breaking the U(1) symmetry associated to the gauge symmetry of the two-form field. The role of this two-form symmetry in hydrodynamics, holography and superfluids was later studied in detail in [41][42][43].
of the counterterms depends on the masses m φ and m B of the fields φ and B respectively, so we will only give them later, once these have been fixed. The action in equation (2.1) is the same as the action used in ref. [16] except for our addition of the λ(B 2 ) 2 self-coupling of the two-form field. As we will show, by tuning the value of λ we can choose some of the IR properties of the model to match expected properties of NLSMs, hopefully making the model more applicable to real-world systems. There are of course many other couplings possible, for example one could add a term proportional to B mn B nk B kl B lm , or terms with higher powers of φ and/or B. A fuller analysis should investigate the effects of such terms.
The equations of motion following from equation (2.1) are where Θ mn is proportional to the bulk stress tensor, explicitly (2.3) To solve these equations of motion we make the black-brane ansatz with all other components of B mn vanishing. Moreover, we take φ to be real. We also take A m = 0 in our ansatz, meaning that the dual field theory will have vanishing net charge density, and will not be subjected to any external electric or magnetic field. Physically, this will later on imply that the Fermi energy lies exactly at the same energy as the nodal line, i.e., our NLSM model is particle-hole symmetric.
In the coordinates used in the ansatz (2.4), the asymptotically-AdS boundary is at r → 0. When the dual field theory is at a non-zero temperature T , there is a horizon at some r = r 0 where f (r) has a double zero while g(r) and h(r) are finite. The temperature T and entropy density s in the field theory are given by the Hawking temperature and the Bekenstein-Hawking entropy density of the horizon, respectively, in units where Boltzmann's constant k B = 1. At T = 0, the spacetime extends to r = ∞. Substituting our ansatz into the equations of motion (2.2), we find four second-order equations of motion and one first-order equation These equations of motion are invariant under the separate, constant rescalings f → Ω f f , g → Ω g g, and (h, B) → Ω h (h, B), reflecting the freedom to rescale the coordinates (t, x, y, z) in the metric ansatz (2.4). As a consequence of these symmetries, the equations of motion imply two radial conservation laws [16]. 4 The first follows directly from a rewriting of equation (2.6a), where on the right-hand side we have used equation (2.5) to evaluate the conserved quantity at the horizon. Note that at T = 0 equation (2.8) implies that f (r) is constant. The second conservation law is Evaluating the factor in square brackets at the horizon, we find that it vanishes at r = r 0 since it is proportional to √ f , and therefore by equation (2.9) it must vanish at all r, Relativistic, (3+1)-dimensional renormalisation group flows obey the a-theorem, guaranteeing the existence of a quantity a that monotonically decreases along the flow [44][45][46][47].
Since non-zero b explicitly breaks Lorentz invariance, the a-theorem does not apply to our system. Nevertheless, following refs. [48][49][50][51] we can use the null energy condition (NEC) to derive two combinations of metric coefficients that must be monotonically decreasing functions of r in our solutions.
The geometric form of the NEC is k m k n R mn ≥ 0, where k m is a null vector field. The first quantity is obtained by choosing the only non-zero components of the null vector to be k t and k r , for which the NEC implies implying that ∂ r C 1 (r) ≤ 0 for any solution of the equations of motion, where Near the boundary at small r, where we can use the asymptotically-AdS boundary conditions f → 1 and g, h → L 2 /r 2 , we find The second quantity is obtained by choosing the only non-zero components of the null vector to be k t and k x , in which case the NEC yields This implies that when T = 0, corresponding to constant f , we must have ∂ r C 2 (r) ≤ 0, where 14) and the second equality is obtained using equation (2.9). At leading order near the boundary we have B ≈ br 2−∆ , where ∆ is the conformal dimension of the operator dual to B, and b is the source for this operator, so that C 2 (r) ≈ (∆ − 2)b 2 r 4−2∆ at small r. The unitarity bound implies ∆ > 2 [52,53], so we find that C 2 (r) diverges as r → 0, perhaps making the monotonicity of C 1 (r) a better candidate for a non-relativistic generalisation of the a-theorem.

Types of solution -zero temperature
At T = 0, equation (2.8) implies that f (r) is constant. The equations of motion (2.6) and (2.7) with constant f (r) admit multiple types of solution, classified by their behaviour in the deep IR r → ∞ [16]. The different types of solution correspond to different phases of the dual field theory.

Topological solutions
The first class of solutions have the large-r behaviour where the equations of motion imply that the constants α, β, γ, and B 0 satisfy This solution is distinguished from the other solutions described below by the fact that the scalar field vanishes as r → ∞, since γ > 0 in equation (2.16), while B/h tends to the non-zero constant value B 0 . The coefficients g 0 and h 0 appearing in equation (2.15) are arbitrary, since the equations of motion are invariant under the rescalings g → Ω g g and (h, B) → Ω h (h, B). This leaves a single parameter family of solutions, parameterised by φ 0 . In this phase, for the values of the couplings used in refs. [16,17] the spectral functions of composite fermionic operators exhibit multiple circular nodal lines. We show in section 3 that the same is true for the values of couplings that we will choose, so we will refer to this phase as topological.
The large-r behaviour of the metric functions written in equation (2.15) imply that the asymptotic metric as r → ∞ has a scaling isometry. Concretely, if we make the coordinate transformation then at large r the metric (2.4) becomes where we have defined L = 2L/α. This metric is manifestly invariant under the combined rescaling (r , t , z ) → Ω(r , t , z ) and (x , y ) → Ω β/α (x, y). Solutions of this type therefore describe an RG flow from a conformally invariant UV fixed point to an IR fixed point with non-relativistic, anisotropic scale invariance with dynamical exponent α/β in the x and y directions and unit dynamical exponent in the z direction. The requirements that B 0 in equation (2.16) is real and that the IR metric in equation (2.18) satisfies the NEC both imply that the dynamical exponent is bounded from below, α/β ≥ 1. The dynamical exponent α/β determines many of the IR properties of this phase. In particular, we will show in section 4 that at small frequency ω, the AC conductivity in the direction orthogonal to the plane of the nodal line scales as σ zz ∝ ω 2β/α−1 . We will choose to engineer a constant σ zz by setting the parameters of our model such that the dynamical exponent is α/β = 2. For example, from equation (2.16) we find that for m 2 B = 1 and λ = 34/13, one has α = 2β = 94/13. The physical motivation for a constant σ zz is that, due to the rotational symmetry around the z axis, we may think of the nodal line semimetal as an infinite collection of graphene-like sheets, with a single sheet for each azimuthal angle, and the conductivity of graphene is known to be finite and nonzero at the charge neutrality point due to particlehole symmetry. While crude, this intuitive picture appears to give correct results for σ zz in real NLSMs [54]. However, the same physical picture would also suggest that σ xx and σ yy would go to the same non-zero constant, which cannot be achieved within the current setup; we always find σ xx = σ yy ∝ ω. A more realistic holographic model should therefore include one or more additional fields to engineer σ xx = σ yy ∝ ω 0 . We leave the exploration of this to future work.

Topologically trivial solutions
This class of solutions has the large-r behaviour where In these solutions, the scalar field φ tends to the constant φ 0 as r → ∞, while B/h → 0 in the same limit. Once more there is a single parameter family of solutions, this time parameterised by B 0 . Setting α = β in equation (2.18), we see that the IR geometry is AdS 5 with curvature radius L = 2L/α. These solutions therefore describe an RG flow between two conformally invariant fixed points. Note that for φ to be dual to a relevant operator we must take m 2 φ < 0, so that from equation (2.20) we see that these solutions only exist for η > 0. This then implies α > 2 and therefore L < L.

Critical solution
Finally there is a single solution with the large-r behaviour 21) where in this case α, β, B 0 , and φ 0 satisfy (2.23) Since g 0 and h 0 are arbitrary, due to the scaling symmetries of the equations of motion, there is only one solution in this class. We will refer to this solution as the critical solution, since it turns out to describe the critical point of a second-order phase transition between the topological and topologically trivial solutions. In the critical solution, both φ and B/h tend to non-zero constants as r → ∞. Like the topological solutions, the critical solution has non-relativistic, anisotropic scale invariance in the IR, but with a different value for the dynamical exponent α/β.

Thermodynamics and one-point functions
We now give expressions for various thermodynamic quantities and one-point functions in the dual field theory. To do so we will need to fix the masses of the fields φ and B mn , as the masses will determine the near-boundary expansions of these fields, and consequently the form of the counterterms. Following ref. [16], we will take m 2 φ = −3 and m 2 B = 1. With these choices, both φ and B mn are dual to operators of conformal dimension ∆ = 3 in the dual field theory, the same as the dimensions as the operatorsψψ andψγ µν ψ in the free field model (1.1), respectively. Of course, interactions will in general renormalise the dimensions of these operators, so our choice of masses has been made purely for concreteness.
For the choice of masses described above, the fields of our ansatz have the nearboundary expansions where the dots denote terms of higher order in the small-r expansion, m, M , b, φ 2 , and b 2 are integration constants, and g 4 and h 4 are given by The coefficients M and b are proportional to the sources for the operators dual to φ and B, while φ 2 and b 2 will determine the corresponding vacuum expectation values. Substituting the expansions (2.24) into equation (2.8), we find that m is related to the temperature and entropy density as With m 2 φ = −3 and m 2 B = 1, the counterterm action appearing in equation (2.1) may be taken to be [16,55,56] where is a small-r cutoff, the field theory (µν) indices are raised and lowered with the induced metric γ µν ,∇ is the covariant derivative calculated with respect to γ, R µν [γ] is the Ricci tensor computed with γ, and R[γ] = γ µν R µν [γ] is the corresponding Ricci scalar. Note that there is a choice of renormalisation scheme implicit in the counterterms (2.27); one could choose to add additional, finite counterterms that would modify the one-point functions of operators by local functions of the sources. 5 The Helmholtz free energy at temperature T is given by F = T I , where I is the Euclidean signature on-shell action. Using the equations of motion, the on-shell action may be reduced to a boundary term. Upon substitution of the near-boundary expansion (2.24), one finds where V is the spatial volume of our system. The one-point function of the scalar operator dual to φ is given by O = δS /δM , where S is the Lorentzian-signature on-shell action. The equations of motion imply that this variation is a boundary term, and using the near-boundary expansions (2.24) we find With our ansatz, the only non-zero components of the one-point function of the antisymmetric tensor operator dual to B mn are O xy = − O yx = 1 2 δS /δb. We find 6 Differentiation of the on-shell action with respect to the boundary metric, yields the one-point function of a symmetric tensor operator [55] T µν ≡ − lim where K µν and K are the extrinsic and mean curvature of the cutoff surface at r = , respectively. The tensorT µν is not quite the stress tensor, since it is not conserved in the presence of a non-zero source b µν . 7 Instead, as we will show in the next subsection, the conserved stress tensor is T µν ≡T µν + 2O µρ b ν ρ . For our black brane ansatz we find T µν = diag(ε, p, p, p), where the energy density ε and pressure p are given by (2.32) The first law of thermodynamics for our system reads dε = T ds − O dM − 2 O xy db. As consistency checks of the expressions in this section, note that p = −F/V , and that ε and p satisfy the expected thermodynamic relation where we obtain the second equality using equation (2.26), and we remind the reader that we are working at vanishing chemical potential. The one-point function of the trace of the stress tensor also satisfies the appropriate Ward identity where β(M ) = −M and β(b) = −b are the beta functions for the sources, and we find the Weyl anomaly to be

Conservation of the stress tensor
In this subsection we show that the operatorT µν defined in equation (2.31) is not the stress tensor, since it does not satisfy the appropriate Ward identity. We will consider a general quantum field theory on a spacetime with metric g µν , which we will later set to the Minkowski metric, in the presence of an external two-form source b µν . Writing the generating functional of connected correlation functions as W [g, b], the one-point function ofT µν may be calculated by functional differentiation of W with respect to g µν , A similar phenomenon with a one-form source is described in ref. [57].
Similarly, the one-point function of the operator O µν sourced by b µν is The method we use to obtain the Ward identity is well known, see for example ref. [58]. We begin by considering the effect of an infinitesimal coordinate transformation generated by a vector ξ µ . The corresponding infinitesimal changes in the metric and source are Demanding that the generating functional is invariant under this coordinate transformation, to linear order in δg and δb we have which implies that the second term on the right-hand side must vanish. Substituting equations (2.36) and (2.37) for the functional derivatives, setting the metric to be flat, and performing an integration by parts, this implies In order for equation (2.40) to hold for any choice of ξ, we find that the factor in square brackets must vanish, yielding the Ward identity The first term on the right-hand side reflects the explicit breaking of translational invariance that arises if b is position-dependent. On the other hand, the second term on the righthand side of equation (2.41) can be non-zero even if translational symmetry is not explicitly broken, in other words even if ∂ µ b ρσ = 0, soT µν is manifestly not the conserved current associated to translational invariance. Instead, the conserved stress tensor is which satisfies the correct Ward identity We show in appendix B that the T µν that we obtain in holography satisfies this Ward identity. Note thatT µν is manifestly symmetric in its indices, since it is calculated from a functional derivative with respect to the symmetric g µν . However, the second term on the right-hand side of equation (2.42) is not symmetric in general, so the stress tensor T µν is not a symmetric tensor. This is an example of the general principle that the stress tensor is asymmetric in any theory with explicitly broken Lorentz invariance [59], see also refs. [57,60].

Numerical results
We now present our numerical results for the thermodynamic quantities and one-point functions of the holographic NLSM. In addition to fixing the masses of the scalar and two-form fields to be m 2 φ = −3 and m 2 B = 1, as discussed in section 2.3, we need to fix the various couplings appearing in the action (2.1). We choose the two-form self-coupling to be λ = 34/13. As discussed in section 2.2 this implies that the dynamical exponent in the topological phase is α/β = 2. With these choices of m 2 φ , m 2 B and λ, we must take the coupling between the scalar and two-form to satisfy σ ≥ 3, so that γ ≥ 0 in equation (2.16), while the scalar self-coupling must obey η > 0 so that φ 0 in equation (2.20) is real. It will be convenient for numerical purposes to choose these couplings so that the various coefficients and exponents in section 2.2 are relatively simple. We will take η = 2 and σ = 8.
At T = 0 we solve the equations of motion numerically by integrating out from the horizon, using the large-r asymptotics in section 2.2 to set boundary conditions. We can then fit the numerical solutions to the near-boundary behaviour (2.24) figure 2 we also plot the free energy and the one-point functions of the scalar and antisymmetric tensor operators as functions of M/b. All three of these quantities are continuous across (M/b) crit. . However, not every physical quantity is continuous across (M/b) crit. . In particular, we find in section 4.1 that the DC electrical conductivity in the z direction has a discontinuous first derivative with respect to M/b at (M/b) crit. , so we identify (M/b) crit. as the location of a second-order quantum phase transition between the two phases dual to the topological and topologically trivial solutions.
At non-zero T we solve the equations of motion by integrating out from the horizon at r = r 0 , imposing the boundary conditions that f has a double zero at r 0 while g, h, φ, and B are finite and non-zero. Since the values of g, h, and the second derivative of f at the horizon may be set to any desired numbers using the scaling symmetries mentioned under equation (2.7), there are only two free parameters in the boundary conditions: the values of φ and B at the horizon. These parameters map to the two dimensionless ratios in the dual field theory, T /b and M/b, through the bulk equations of motion.    To derive these different exponents and values for v 2 s , we construct an approximate near-horizon solution for f (r) valid at low temperatures by solving equation (2.6a) using the zero temperature, large-r results g(r) ≈ g 0 r −α and h(r) ≈ h 0 r −β , where the use of the large-r forms is justified by the fact that sending T → 0 is equivalent to sending r 0 → ∞. The resulting near-horizon solution for f is f (r) ≈ f 0 1 − (r/r 0 ) α+β 2 for some integration constant f 0 . Substituting these solutions into equation (2.5) we can determine how the temperature and entropy density scale with r 0 at leading order at low temperatures, (2.45) From the r 0 dependence of these two results, we see that s ∝ T 1+2β/α at small T , and consequently For our choice of couplings and masses, the topological phase has β/α = 1/2, the critical phase has β/α 0.745, and the topologically trivial phase has β/α = 1. Substituting these ratios into equation (2.46), we reproduce the low temperature results for v 2 s quoted above. The result for v 2 s in equation (2.46) may be understood as arising from the nonrelativistic scaling symmetry of the infrared fixed points. Recall from the discussion in section 2.2 that the asymptotic metric at large r is invariant under the rescaling (t, z) → Ω(t, z) and (x, y) → Ω β/α (x, y), with a different value of β/α in each phase. The corresponding Ward identity is then T t t + T z z + β α ( T x x + T y y ) = 0, which we can solve to find the pressure as a function of the energy density, p = ε/(1 + 2β/α). Substituting this into the formula (2.44) for v 2 s , we recover equation (2.46).

Fermion equations of motion and boundary conditions
In this section we calculate the spectral function of a composite fermionic operator dual to a pair of probe fermions in the gravitational background described in the previous section.
The spectral function is defined as where G R (k) is the retarded two-point function of the fermionic operator, at four-momentum k. Nodal lines should appear in the spectral function as rings of sharp peaks in the (k x , k y ) plane near zero frequency. Since our objective is simply to demonstrate the presence of these peaks in the low M/b phase, for simplicity we will work at M = 0, meaning that we set the scalar field φ of the background to zero. We compute the fermion Green's function holographically following refs. [33,[61][62][63][64]. In order to introduce a Dirac fermion on the boundary we need two bulk Dirac fermions Ψ 1,2 . We denote the five-dimensional Dirac matrices as Γ m , they satisfy {Γ m , Γ n } = 2G mn . Denoting the vielbeins as e m a , where we use underlines to denote the tangent space indices, 8 We also define Γ ab = 1 2 Γ a , Γ b . It will be useful to define projectors P ± = 1 2 (1±Γ r ). On the boundary Γ r is the chirality operator, and these projectors isolate the right-and left-handed components of Dirac fermions. In the bulk, we define Ψ I,± = P ± Ψ I . We choose Γ r to be hermitian.
We take the action for the fermions to be 9 with ω mab the spin connection. The corresponding equations of motion are Our sign conventions for the vielbeins are The minus sign in e r r is convenient for comparison to other works that place the boundary at r = ∞ rather than r = 0. 9 The choice of how to couple the bulk fermions to the two-form field B is not unique. However, it appears to be the only choice that leads to nodal lines in the spectral function [17].
The boundary term in equation (3.2) is evaluated at some small-r cutoff . We leave implicit in our expressions that is to be taken to zero at the end of calculations. The form of the boundary term is chosen such that the variational principle requires the boundary values of Ψ 1,+ and Ψ 2,− to be fixed. Holographically, we therefore interpret the boundary values of Ψ 1,+ and Ψ 2,− as the right-and left-handed components of a source for a composite Dirac fermion operator in the boundary field theory, with the conjugate momenta Ψ 1,− and Ψ 2,+ determining the corresponding one-point functions. It will be convenient to repackage the spinors into i.e. Ψ contains the sources and Π contains the one-point functions.
The two-point functions of the boundary fermion operator are determined from the on-shell action S . Since the bulk part of the action vanishes when the equations of motion are satisfied, only the boundary term in equation (3.2) contributes, yielding To obtain the momentum space Green's function, we Fourier transform in the field theory directions, writing Ψ(x, r) = and similar for Π. The momentum space versions of the equations of motion (3.2) may then be solved for Π, yielding Π(k; r) = −iξ(k; r)Ψ(k; r) for some matrix ξ(k; r). The on-shell action becomes Applying the Minkowski space correlator prescription of refs. [65,66], the matrix of fermion Green's functions is then where the factor of r −2m f L arises because the equations of motion (3.3) imply that Ψ ≈ r 2−m f L near the boundary, while √ −γ ≈ r −4 . In order to actually determine ξ(k; r), and therefore the fermion Green's functions, it will be convenient to rescale the spinors. We define where c is some arbitrary reference point, and F (r) = f 4f + g 2g + h 2h . Since Ψ and Π are rescaled by the same amount, we have ψ(k; r) = −iξ(k; r)χ(k; r). The rescaling effectively eliminates the spin connection from the equations of motion for ψ and χ. Indeed, applying the projectors P ± to the fermion equations of motion (3.3), it is straightforward to show that ψ and χ satisfy where / k = e µ a Γ a k µ . Finally, substituting χ = −iξψ into the second line of equation (3.9), and using the first line to eliminate first derivatives of ψ, we obtain a first-order matrix equation for ξ(k; r), Since we expect to see the nodal lines at k z = 0, using rotational symmetry in the (x, y) plane we can take k µ = (ω, k x , 0, 0). For this choice of k µ , we can solve the first line of equation (3.9) for χ to determine the matrix structure of ξ, finding where the coefficients ξ 0,1,2,3 satisfy four coupled first-order differential equations. These equations may be decoupled by defining the four independent linear combinations where s 1 = ±1 and s 2 = ±1 are uncorrelated signs. In terms of these variables, it is straightforward to show that equation (3.10) becomes, where we have substituted the explicit expressions for the vielbeins. The variables ξ s 1 ,s 2 are holographically dual to the eigenvalues of the fermion Green's function; substituting the expansion (3.12) into equation (3.7), one finds that the four (3.14) The fermion spectral function (3.1) is proportional to the sum of the imaginary parts of these eigenvalues, allowing us to determine various properties of the spectral function by inspection of the equations of motion (3.13). For example, sending Λ → −Λ maps the equations of motion of ξ ±,± to those of ξ ±,∓ , so the spectral function is invariant under this change in the sign of Λ. Similarly, the spectral function is invariant under k x → −k x , since this interchanges the equations of motion for ξ ±,± and ξ ∓,∓ .
Boundary conditions We now determine the appropriate boundary conditions for ξ s 1 ,s 1 , beginning at T = 0. In the deep IR we have For the low M/b phase α > β, so as r → ∞ the terms proportional to ω/ g(r) dominate the equations of motion (3.9), such that they become approximately We can use the first equation to eliminate χ from the second, yielding a second-order equation of motion just for ψ, for some constant spinor ψ 0 . In order to compute the retarded Green's function we impose ingoing boundary conditions on ψ at r → ∞, corresponding to the plus sign in the exponent.
If we then substitute this solution back into the first equation in (3.15), we find χ ≈ Γ 0 ψ. Comparing to χ = −iξψ we then find ξ ≈ iΓ 0 . This means we must impose ξ s 1 ,s 2 = i as r → ∞ on each of the decoupled variables. The same procedure can be used to find the boundary conditions at non-zero temperature. Near the horizon, where f (r) has a double zero, the equations of motion (3.9) are dominated by the terms proportional to ω/ f (r)g(r). Taylor expanding near the horizon and using the expression for the Hawking temperature in equation (2.5), we find that near the horizon the equations of motion become Again using the first equation to eliminate χ from the second, we find which we can solve to obtain the near-horizon behaviour ψ ≈ (r 0 − r) ±iω/2πT ψ 0 for some constant spinor ψ 0 . Ingoing boundary conditions correspond to choosing the minus sign in the exponent. Substituting the ingoing solution into the first equation in (3.18) we then find the near horizon behaviour of χ to be χ ≈ Γ 0 ψ, as for at T = 0. So the appropriate boundary conditions are ξ s 1 ,s 2 = i at the horizon. Suppose we instead chose outgoing boundary conditions for the bulk fermion, meaning we would compute the advanced, rather than retarded, Green's function. Following the same steps as presented above, we find that the appropriate boundary condition on ξ becomes ξ s 1 ,s 2 = −i, either at r → ∞ for T = 0 or at r = r 0 for T = 0. In other words, to compute the advanced Green's function we flip the sign of the boundary condition. Notice that the left-and right-hand sides of equation (3.13) are odd and even functions of ξ s 1 ,s 2 , respectively. Moreover, the sign of the right-hand side may be inverted by simultaneously sending ω → −ω and s 1 → −s 1 . Taken together, these two facts imply that if a given set of functions ξ ±,± solve equation (3.13) at frequency ω, then −ξ ∓,± is also a solution, but at frequency −ω. If ξ ±,± satisfies ingoing boundary conditions, so that its r → 0 limit yields an eigenvalue of the retarded Green's function, then −ξ ∓,∓ satisfies outgoing boundary conditions, and its r → 0 limit yields an eigenvalue of the advanced Green's function. Since the retarded and advanced Green's functions are related to each other by complex conjugation, we therefore find that the eigenvalues of the retarded Green's function satisfy This means that the spectral function, being proportional to the sum of the imaginary parts of the eigenvalues, satisfies ρ(ω, k x ) = ρ(−ω, k x ). (3.21)

Numerical results
We compute the fermion spectral function by numerically solving the equation of motion (3.13), imposing ingoing boundary conditions ξ s 1 ,s 2 = i at the horizon (or at r → ∞, in the case of T = 0). The spectral function is determined by subsituting the resulting numerical solution for ξ into equation (3.7). We expect any nodal lines to appear as sharp peaks in the spectral function at ω = 0 and k x = K, for some K = 0. By rotational symmetry in the (x, y) plane, such a peak would exist everywhere along the circle k 2 x +k 2 y = K 2 , hence forming a nodal line. To obtain numerical results we will need to choose definite values of the fermion mass m f and the coupling Λ to the two-form field. We will choose m f = −1/4L and Λ = 1. Figure 4 shows our results for the fermion spectral function at T = 0. Since the spectral function is invariant under k x → −k x , we only show results for k x ≥ 0. In all of the plots in the figure, the frequency has been given a small imaginary part Im ω/b = 10 −4 . This broadens the peaks in the spectral function, making them easier to resolve numerically. Outside the light cone, we expect these peaks to become delta functions at Im ω = 0. Figure 4a shows the spectral function in units of 1/ √ b as a function of k x /b, at Re ω = 0. We observe multiple nodal lines, visible as the extremely sharp peaks in the figure. Figure 4b shows the same data as figure 4a but with logarithmic axes, making it easier to see most of the peaks. With the logarithmic k x /b axis the peaks appear equally spaced. Indeed, from a fit to our results we find that the locations of the peaks (i.e. the radii of the nodal lines) are very well approximated by the Efimov-like spectrum K n 1.65e −0.782 n b, where the integer n labels the different peaks, starting from the outermost at n = 0. We find a total of seven peaks in the spectral function, 10 however it is plausible that this finite number of peaks arises due to the artificial broadening introduced by the non-zero imaginary part of the frequency, which washes out the very closely-spaced peaks at large n, corresponding to small K n /b, and that the spectrum of peaks may continue to n → ∞. Figure 4c is a density plot of the fermion spectral function in part of the (k x /b, Re ω/b) plane, showing how some of the peaks evolve as we move away from Re ω = 0. Our results for the spectral function are invariant under Re ω → − Re ω, as expected from equation (3.21). We find that each peak splits in two at non-zero Re ω, with one of the daughter peaks moving to small k x as we increase Re ω, while the other moves to larger k x . This splitting qualitatively resembles the momentum dependence of the inner two eigenvalues of the non-interacting toy model, plotted in figure 1a. Also visible in figure 4c is an apparent continuum of states, where the spectral function is non-zero but varies smoothly, rather than exhibiting a peak (this is the region coloured purple in the plot). We expect that this is an artifact of the broadening induced by the non-zero imaginary part of the frequency.
The behavior of the individual four components of the spectral function clearly shows the chiral nature of the excitations around the nodal line. We can straightforwardly use equation (3.14) to compute the individual eigenvalues λ s 1 ,s 2 of the fermion Green's function. In figure 5 we show density plots of the imaginary parts of these eigenvalues over the same range of Re ω and k x as figure 4c. Notice that the imaginary parts of λ ±,± and λ ∓,± are related by Re ω → − Re ω, as expected from equation (3.20). From figure 5 we see that each nodal line is formed by one of the two pairs of eigenvalues: half of the nodal lines are formed by the intersection of peaks in λ ++ and λ −+ at Re ω = 0, while the other half are formed by the intersection of peaks in λ +− and λ −− , with each pair alternating as we increase momentum from k x = 0. Similar behaviour was seen in earlier holographic work [17]. Figure 6 shows our numerical results for the fermion spectral function at non-zero temperature. The plots in the figure show the fermion spectral function in units of 1/ √ b as a function of k x /b at ω/b = 10 −4 i, each at a different value of T /b. We find that as the temperature is increased, the peaks broaden and merge into a continuum, starting at small k x /b and moving outwards. Eventually, for sufficiently large temperatures there are no peaks in the spectral function at all. We find that this occurs for T /b 1.5. Notice that there is a range of temperatures for which there is only a single sharp peak in the fermion spectral function.

Conductivity
In this section we compute the electrical conductivities of our system, using the Kubo formula where AB R (ω, k) denotes the retarded Green's function of operators A and B at frequency ω and momentum k. To compute this two-point function holographically we consider linearised fluctuations of the bulk gauge field A m . We take the fluctuations to be of the form A m (t, r) = dω 2π e −iωt A m (ω; r). In the radial gauge A r = 0, the equations of motion fix A t to be a constant, which we set to be zero by a further gauge transformation. The equations of motion for A x,y,z are then Using the near boundary expansions of f , g, and h written in equation (2.24), one finds that for small r, solutions to the gauge field equations of motion (4.2) take the form i (ω) determined by the boundary conditions. The on-shell action for the gauge field fluctuations is obtained from equation (2.1). Using the equations of motion for the A i , it reduces to a boundary term Applying the Lorentzian correlator prescription of refs. [65,66], we can then read off the expressions for the current-current Green's functions at zero momentum,  Rotational symmetry in the (x, y) plane implies that J x J x R (ω, 0) = J y J y R (ω, 0). We will therefore focus only on the calculation of J x J x R and J z J z R . In order to compute the retarded Green's functions we must impose ingoing boundary conditions on the fluctuations of the gauge field at the horizon. Near the horizon, we find solutions to the equations of motion (4.2) take the form Ingoing boundary conditions correspond to choosing the minus sign in the exponent. At T = 0 there is no horizon. In this case, at large r where g ≈ g 0 r −α and h ≈ h 0 r −β we find that the gauge field equations of motion have the approximate solutions where ingoing boundary conditions correspond to the plus sign in the exponents.

Substituting the expressions for the Green's functions in equation (4.5) into the Kubo formula (4.1) and making use of the near-boundary expansions written in equation (4.3), it is straightforward to show that the DC (zero frequency) limits of the conductivities are
(4.8) The equations of motion (4.2) imply that the combinations r √ f g ∂ r A x (ω; r)/iωA x (ω; r) and r √ f h ∂ r A z (ω; r)/iωA z (ω; r) are independent of r up to O(ω) corrections [67]. We can therefore relax the r → 0 limits in equation (4.8), instead evaluating the right-hand sides at any convenient value of r.
At T = 0 we evaluate the conductivities at r → ∞. Using the large-r behaviour of the gauge fluctuations written in equation (4.7) we find at large r. Substituting this solution into equation (4.8) we find that the DC conductivities in the plane of the nodal line vanish for all phases, since g(r) ≈ g 0 r −α at large r, with α ≥ 0. On the other hand, the DC conductivity normal to the plane of the nodal line depends on the phase of the system through the exponents α and β, In the trivial and critical phases, we have α < 2β, so that σ DC zz = 0. However, recall that by tuning λ = 34/13 we were able to fix α = 2β in the topological phase, which then gives a finite, non-zero value for σ DC zz , , topological phase, 0, other phases.

(4.11)
This behaviour is very special to our particular choice of λ. For λ < 34/13 we find α > 2β, and consequently σ DC zz diverges. Conversely, for λ > 34/13 we find α < 2β, leading to σ DC zz = 0. At non-zero temperature it is convenient to evaluate equation (4.8) at the horizon. Using the near-horizon behaviour written in equation (4.6) we find We can use these expressions to determine the DC conductivities from the numerical solutions presented in section 2.5. We can determine the low-temperature behaviour of the DC conductivities using the T = 0 solutions g(r) ≈ g 0 r −α and h(r) ≈ h 0 r −β in equation (4.12) and making use of  the relationship between r 0 and T , valid at low temperatures, given in equation (2.45). We find σ DC xx ∝ T for all values of M/b, while the behaviour of σ DC zz for a given value of M/b depends on which phase the system is in at T = 0, through the dynamical exponent α/β, σ DC zz ∝ T 2 β α −1 . In the topological phase we have α = 2β, and consequently σ DC zz is finite and non-zero as T → 0, given by the result in equation (4.11). In the topologically trivial phase we have α = β, leading to σ DC zz ∝ T . Finally, in the critical phase we have α/β 1.34, implying σ DC zz ∝ T 0.491 .
To obtain the DC conductivities away from small T we must evaluate equation (4.12) numerically. In figure 7 we plot our numerical results for the DC conductivity. Figure  This also has the advantage that one could compute the derivative without knowing the DC conductivity exactly at T = 0, so it may be more easily experimentally accessible. However, it would be computationally intensive to evaluate ∂σ DC zz /∂T numerically for enough data points to obtain a plot comparable to figure 8, and other than the cusp the qualitative features of such a plot should not be too different, so we have only computed ∆σ DC zz .

AC conductivity
We now compute the AC conductivities of our system. To keep the discussion concise we will restrict to T = 0. We will begin by deriving approximate formulas for the AC 0.9493, suggesting that the DC conductivity may provide a probe of the quantum phase transition at (M/b) crit. . conductivities at small and large frequencies. We will then present numerical results for a wider range of frequencies.
Small frequency Using the near-boundary expansions in equation (4.3), it is straightforward to rewrite the imaginary parts of the Green's functions (4.5) as where and we have used that A i (ω; r) * = A i (−ω; r), as required by reality of A i (t, r). Notice that we have not written r → 0 limits in equation (4.14). This is because the equations of motion (4.2) imply that the right-hand sides of equation (4.15) are independent of r, allowing us to evaluate F x (ω) and F z (ω) at any value of of r. It will be convenient to evaluate them in the limit r → ∞.
We can determine the small-frequency (ω M, b) behaviours of the real parts of the conductivities following ref. [68]. These are determined by the behaviour of the solutions in the deep IR, where g ≈ g 0 r −α and h ≈ h 0 r −β . In this region, the equations of motion (4.2) become (4.16) For ω > 0, the solutions to these equations of motion obeying ingoing boundary conditions at r → ∞ are where H n should be replaced by Hankel functions of the second kind H (2) n . Substituting these solutions into equation (4.15), we find that the r-independent fluxes are F x = 2g 0 α/π and F z = 2h 0 α/π. The real parts of the conductivities are then We now need to determine the near-boundary coefficients A (0) x (ω) and A (0) z (ω). In general these should be determined by matching the IR solutions in equation (4.17) to asymptotic solutions computed at small r. However, in the simple case under consideration the small-r solution is just A x,z (ω; r) = constant, so we can simply obtain A (0) x,z (ω) by taking the r → 0 limit of the solutions in equation (4.17), yielding A z (ω) = −iΓ(β/α)(α √ g 0 /Lω) β/α /π. Substituting these results into equation (4.18), we find that the real part of σ xx is linear in ω for all phases, while the real part of σ zz depends on the phase of the system through the exponents α and β, for ω M, b. In the topologically trivial phase α = β, and so Re σ zz is also linear in ω at small frequency. For our choice of λ = 34/13, the topological phase has α = 2β, and Re σ zz reduces to the DC conductivity written in equation (4.11). Finally, in the critical phase we have α 2.327 and β 1.7346, leading to Re σ zz ∝ ω 0.491 .
Large frequency At large frequencies, ω M, b, the conductivity is determined by the physics of the UV fixed point. Holographically, this means that we can find the large ω limit of the conductivities by taking the metric to be that of AdS 5 , i.e. f (r) = 1 and g(r) = h(r) = L 2 /r 2 . The equations of motion for all three gauge field components are then the same, For ω > 0, the solution obeying ingoing boundary conditions at r → ∞ is [π − i + 2iγ E + 2i log(ωL/2)] + . . . , (4.22) where γ E = 0.577... is the Euler-Mascheroni constant. Reading off the coefficients A (0) i (ω) and A (2) i (ω) and substituting into equation (4.5) we obtain the Green's functions, from which we obtain the conductivities using the Kubo formula (4.1), 11 Numerical results Away from the limits of small or large ω we obtain the conductivities by solving the equations of motion (4.3) numerically, using the ingoing boundary conditions at large r written in equation (4.7). We then perform a fit to the solution at small r to obtain the near-boundary coefficients A (0) i (ω) and A (2) i (ω), which determine the conductivities through equations (4.1) and (4.5).
Since we are most interested in the physics of the topological phase, for simplicity we restrict to M = 0. Logarithmic plots of our results for σ zz and σ xx as functions of ω/b are shown in figure 9. The thick black curves are our numerical results, while the solid grey lines show the small frequency approximation (4.19) and the dashed grey lines show the large frequency approximation (4.23). Both approximations work well within their regimes of validity. An absolute value has been taken in the figures 9b and 9d, as the imaginary parts of the conductivities are not of fixed sign. Both Im σ zz and Im σ xx are negative at small frequencies and positive at large frequencies.

Thermal conductivity
We now calculate the thermal conductivity matrix κ, given by [69] 12 where in T i0 T j0 R (0, k → 0) one should send the frequency to zero first, followed by momentum. Due to the asymmetry of the stress tensor, the order of indices is crucial. Replacing either T i0 with T 0i and/or T j0 with T 0j yields a Kubo formula for a different transport coefficient, special to systems with broken Lorentz invariance, that is discussed in ref. [70].
At b µν = 0, the limits of Green's functions appearing in equation (4.24) are fixed by a Ward identity [71]. This is not the case at non-zero b µν . However, it will still be useful to see what the Ward identities of our system tell us about the thermal conductivity. 11 The factor of L inside the logarithm in equation (4.23) arises from the holographic renormalisation, and can be shifted by the addition of a finite counterterm proportional to Fµν F µν . For comparison to condensed matter systems, one should replace this factor L with some typical ultraviolet length scale. 12 The thermoelectric conductivity αij vanishes for our system, since we have vanishing net charge density.

Ward identities
As we show in appendix C, translational symmetry implies that two-point functions of the stress tensor T µν and antisymmetric tensor O µν satisfy the Ward identities 13 where k µ = (−ω, k) is the four-momentum. At k = 0, the (ν, ρ, σ) = (i, j, 0) component of the Ward identity (4.25a) reads so for non-zero ω we obtain T 0i T j0 R (ω, 0) = δ ij ε. We wish to relate this two-point function to the one appearing in the thermal conductivity (4.24), which involves T i0 rather than T 0i . We can do so using the relation T µν =T µν + 2O µρ b ν ρ from equation (2.42). Since the operatorT µν sourced by the metric is a symmetric tensor, this implies that where we have made use of the fact that b 0 µ = 0 for our system. We then have Now consider the (ν, ρ, σ) = (j, 0, i) component of equation (4.25b) at k = 0, which for ω = 0 implies T 0j O 0k R (ω, 0) = 0. Relating T 0j to T j0 through equation (4.27), we find Finally, for ω = 0 the (ν, ρ, σ) = (0, j, 0) component of the Ward identity is Assuming that the k → 0 limit is smooth, we then have T i0 T j0 R (0, k → 0) = −δ ij p. Substituting this into the Kubo formula (4.24) and using equation (4.29), we obtain where we have used the thermodynamic relation ε + p = T s. For the holographic NLSM, the only non-zero components of the two-form source are b xy = −b yx = b. The form of the thermal conductivity in the direction orthogonal to the plane of the nodal line is therefore unchanged from the well-known b = 0 result, κ zz = is/ω for ω = 0 [74]. Through the Kramers-Kronig relations, the 1/ω pole in the imaginary part of the thermal conductivity implies that there should be a delta-function contribution to the real part at ω = 0, arising from conservation of momentum. The final expression for κ zz , valid at all frequencies, is then On the other hand, the thermal conductivities in the plane of the nodal line depend on the two-point function of O 0y at zero momentum, where κ xx = κ yy due to rotational symmetry in the (x, y) plane, and we have included a delta-function contribution at ω = 0, arising from the Kramers-Kronig relations. We will calculate the two-point function O 0y O 0y R (ω, 0) holographically. We find that this twopoint function is non-zero at T = 0, and hence it provides the dominant contribution to equation (4.33) at low T /b. We therefore find that the thermal conductivities κ xx = κ yy in the plane of the nodal line are significantly enhanced compared to the out-of-plane thermal conductivity κ zz at low temperature. In particular, the Drude weight of κ xx diverges as 1/T as T → 0.
One should be slightly careful when applying the Kramers-Kronig relations to κ xx (ω) in order to obtain the coefficient of δ(ω). The relations only hold for functions which vanish when |ω| → ∞ with Im ω > 0, while the two-point function O 0y O 0y R (ω, 0) may diverge at large ω. One should first subtract these divergences before applying the Kramers-Kronig relations. However, the subtraction does not change the coefficient of the 1/ω pole in the first line of equation (4.33), and therefore does not change the coefficient of δ(ω).

Holographic computation
We wish to use holography to compute the two-point function O 0y O 0y (ω, k = 0) that determines κ zz (ω). We will need to consider linearised fluctuations of the metric δG mn and two-form δB mn . The time dependence of the fluctuations will be written as δG mn (t, r) = dω 2π e −iωt δG mn (ω; r), and similar for δB mn . In a gauge in which δG rm = 0, there are three fluctuation components relevant for computing the two-point function in question: δG tx , δB ty , and δB ry . They couple to each other, and not to any other fluctuations. It will be convenient to define G tx = δG tx /h, B ry = δB tr /L 2 , and B ty = rδB ty /L 2 . The linearised equations of motion for these fluctuations are Note that B ry is fixed algebraically by these equations. In order to solve the equations of motion and determine the two-point functions, it is useful to work with linear combinations of the fluctuations that are invariant under gauge transformations of the background solution [75]. Under an infinitesimal diffeomorphism generated by a vector ξ m , the metric and two-form fluctuations transform as (4.35) We find that B ry is invariant under the set of such diffeomorphisms that preserve the radial gauge condition δG rm = 0, while G tx and B ty transform non-trivially. We can form a gauge-invariant combination Z(ω; r) = G tx (ω; r) − L 2 rB B ty (ω; r). The equations of motion (4.34) imply that Z satisfies a second order ODE, that does not depend on the other fluctuations, The forms of the coefficients c 1 and c 2 appearing in this equation are rather complicated, we give them explicitly in appendix D. Once a solution for Z has been found, the remaining fluctuations are determined by the equations which follow from equation (4.34).
Near the horizon at r = r 0 , we find that solutions to the equation of motion (4.37) take the form Z(ω; r) ∝ (r 0 − r) ±iω/2πT . To determine the retarded Green's functions we impose ingoing boundary conditions at the horizon, corresponding to choosing the minus sign in the exponent. Near the boundary, Z has the small-r expansion with coefficients Z (0) (ω) and Z (2) (ω) determined by the boundary conditions. Crucially, the leading-order coefficient in this expansion is a linear combination of the boundary values of the metric and two-point fluctuations, ty ≡ B ty (ω; r → 0). Expanding the action (2.1) to quadratic order in the fluctuations, and using the equations of motion, we find that the on-shell action reads where the coefficients are Applying the Minkowski space correlator prescription of refs. [65,66], we can then read off expressions for the two-point functions of the operatorsT 0x and O 0y , dual to G tx and B ty , respectively, at zero momentum where we use equation (4.41) to relate Z (0) to the sources G (ω, 0) = 0, as demanded by the Ward identity (4.25b). Further, using equation (2.32) we can confirm as expected from the Ward identity (4.26). 14 We can now compute the thermal conductivity κ xx holographically by solving the equation of motion (4.37) for Z(ω; r) numerically, imposing ingoing boundary conditions at the horizon. Fitting the resulting solution to the near boundary expansion (4.40) at small r, we can determine Z (0) (ω) and Z (2) (ω). These coefficients may then be substituted into equation (4.44) to obtain the two-point function O 0y O 0y R (ω, 0), which then determines κ xx through equation (4.33).
For simplicity, we have only computed the thermal conductivity at M = 0. We show our numerical results for κ xx as a function of frequency in figures 10a and 10b, for different sample values of T /b. In contrast to when b = 0, we find that the real part of κ xx is non-zero for ω = 0, taking a finite, approximately ω-independent value at ω/b 1, and growing as Re κ xx 3.2 b 2 ωL 3 /16πG N T when ω/b 1. In between these two regimes we observe an intermediate region at ω/b of order one, in which Re κ xx grows more rapidly with frequency. The size of this intermediate region decreases with increasing T /b.
At small ω/b 1 we find Im κ xx ∝ 1/ω, as for at b = 0. However, the proportionality coefficient is no longer just the entropy density, instead receiving a contribution from the real part of the two-point function appearing in equation (4.33). At large ω/b, the imaginary part grows as Im κ xx 1.9 b 2 ω log(ω/b)L 3 /16πG N T , with a proportionality coefficient independent of T /b. At low temperatures, we find that the imaginary part of κ xx T is approximately independent of T /b: notice that the solid black and dashed orange curves coincide in figure 10b.
14 The reproduction of the two-point functions fixed by the Ward identities only occurs due to the presence of the counterterm involving∇ λ Bµν∇ λ B µν in equation (2.27). Without this term present, we find T 0x O 0y R(ω, 0) ∝ ω 2 and T 0x T 0x R − ε ∝ ω 2 . In other words, without this counterterm, the two point functions that we compute holographically satisfy Ward identities that differ from those in equation (4.25) by contact terms. We note that T x0 T x0 R(ω, 0), and therefore the thermal conductivity κxx, is independent of the coefficient of this counterterm. The two separate contributions to the Drude weight as functions of T /b. The solid black curve shows the entropy density contribution πs. The dashed orange curve shows the two-point function contribution 4πb 2 O 0y O 0y R (ω → 0, 0)/T . Although the ω → 0 limit of the twopoint function is real, it is not of fixed sign, being positive for small T /b and negative for large T /b. We therefore plot its absolute value. The two-point function makes the dominant contribution to the thermal conductivity at low temperatures, while the entropy density dominates at high temperatures.
The approximate independence of T Im κ xx on temperature at T /b 1 arises because at low temperatures the thermal conductivity is dominated by O 0y O 0y R (ω, 0) in equation (4.33), which is non-zero as T → 0. From equation (4.33), this has the significant consequence that the Drude weight diverges as 1/T in the limit T → 0, with a coefficient determined by the zero frequency limit of the two-point function O 0y O 0y R (ω, 0). We plot the Drude weight as a function of T /b in figure 10c, which clearly shows the expected 1/T divergence as T → 0, as well as a T 3 divergence as T → ∞.
In figure 10d we plot the two separate contributions to the Drude weight as functions of T /b. As temperature goes to zero, the entropy density vanishes as s ∝ T 2 for fixed b, and is therefore much smaller than b 2 O 0y O 0y R (ω → 0, 0)/T , which diverges as 1/T . Conversely, at high temperatures the entropy density grows as s ∝ T 3 , whereas the twopoint function grows more slowly, b 2 O 0y O 0y R (ω → 0, 0)/T ∝ T . The entropy density therefore dominates in this limit, so we have Im κ xx Im κ zz at high temperatures.

Shear viscosity
In this section we compute the shear viscosities of the nodal line system from the Kubo formula The indices should be such that i = j and (k, l) = (i, j) or (j, i).
For rotationally invariant holographic systems, all shear viscosities take the universal value η = s/4π [34][35][36][37][38]. Due to the anisotropy introduced by b this is not the case for the nodal line system. As discussed in section 2.4, a second consequence of non-zero b is that the stress tensor is not symmetric in its indices. Concretely, from equation (2.42) we find that two-point functions of the stress tensor take the form Through the Kubo formula (4.46) this implies that η ij,kl , η ij,lk , and η ji,lk are not necessarily equal. In a hydrodynamic expansion of the stress tensor, the stress tensor's asymmetry implies that there are more possible tensor structures in an anisotropic system compared to a relativistic system. We expect different linear combinations of η ij,kl , η ij,lk , and η ji,lk will couple to different tensor structures, although we leave an analysis of the hydrodynamics of our system to future work. To compute two-point functions of the stress tensor holographically, we consider linearised fluctuations δG mn of the metric [65]. The components of the metric fluctuations relevant for the computation of the shear viscosity are δG xy , δG xz , and δG yz , which decouple from other metric fluctuations due to the subgroup of rotational symmetry preserved by our system [75,76]. These three metric fluctuations are also decoupled from one another. However, at non-zero b, δG xz and δG yz do couple to certain components of linearised fluctuations δB mn of the two-form field: δG xz couples to δB yz , and δG yz couples to δB xz .
Since the Kubo formula (4.46) involves the two-point function at zero momentum, throughout this subsection we will assume our linearised fluctuations have the spacetimedependence δG mn (t, r) = dω 2π e −iωt δG mn (ω; r), and similar for δB mn . It will also be conve-nient to define G ij = δG i j , i.e. G xy = δG xy /h, G xz = δG xz /h, and G yz = δG yz /g. Similarly we define B ij = r −1 δB i j , where the prefactor of r −1 is chosen so that B ij goes to a constant at the boundary.

η xy,xy
We begin by computing η xy,xy , dual to the metric fluctuation G xy . Note that since the only non-zero component of b µν in our system is b xy = −b yx , from equation (2.42) we have T xy = T yx , and therefore η xy,xy = η yx,yx = η xy,yx . We will show that the shear viscosity in this channel takes the universal value for isotropic, holographic systems, i.e. η xy,xy = s/4π.
The metric fluctuation G xy satisfies the equation of motion which is just the equation of motion for a massless scalar field in the black-brane background (2.4). Near the boundary, solutions to this equation of motion take the form with coefficients G (0) xy (ω) and G (4) xy (ω) fixed by the boundary conditions. Expanding the action (2.1) to quadratic order in G xy and using the equation of motion (4.48), we find that the on-shell action for fluctuations in this channel is where the dots denote real contact terms. Applying the Lorentzian correlator prescription of refs. [65,66], the retarded two-point function of T xy at zero momentum is then (4.51) where the dots denote real contact terms descending from the dots in equation (4.50). Since the contact terms are real, they make no contribution to the Kubo formula (4.46) and can be neglected. A useful way to rewrite the two-point function expression is where we have made use of reality of G xy (t, r), which implies that G xy (−ω; r) = G xy (ω; r) * . The imaginary part of the Green's function may then be written as where F(ω) = i r L f gh [G xy (−ω; r)∂ r G xy (ω; r) − G xy (ω; r)∂ r G xy (−ω; r)] . (4.54) The equation of motion (4.50) implies that the right-hand side of this expression is independent of r, so we may evaluate F(ω) at any convenient value of r. We will evaluate it at the horizon at r = r 0 . Near the horizon, the two independent solutions to the equation of motion take the form G xy (ω; r) = c(ω)(r 0 − r) ±iω/2πT . To obtain the retarded two-point function we should choose ingoing boundary conditions, corresponding to the minus sign in the exponent. Then, substituting this solution into equation (4.54) and using equation (2.5) we find F(ω) = −8G N sω|c(ω)| 2 , and therefore (4.55) Substituting this into the Kubo formula we find η xy,xy = (s/4π)|c(0)| 2 /|G We now compute the shear viscosities η xz,xz , η zx,zx , and η xz,zx . 15 These components of the shear viscosity may be computed from the linearised fluctuations G xz and B yz , which satisfy the coupled equations of motion where the coefficient C is given by (4.58) 15 By rotational symmetry in the (x, y) plane, these will be equal to ηyz,yz, ηzy,zy, and ηyz,zy, respectively.
Near the boundary, solutions to the equations of motion take the form xz , B yz (ω), and B (2) yz (ω) determined by the boundary conditions. Expanding the action (2.1) to quadratic order in G xz and B yz , and using the equations of motion (4.57), we find that the on-shell action for fluctuations is 60) where the dots once again denote real contact terms that do not contribute to the shear viscosity. From this expression we read off expressions for the retarded two-point functions ofT xz and O yz using the Lorentzian correlator prescription of refs. [65,66]. Using equation (4.47) we then arrive at expressions for the two-point functions of T xz and T zx , where the dots denote terms descending from the dots in equation (4.60).
To obtain the Green's functions and shear viscosities numerically we use the method of ref. [77]. We begin by constructing two linearly independent sets of solutions to the equations of motion (4.57) that satisfy the ingoing boundary conditions G xz = (r 0 − r) −iω/2πT and B yz = ±(r 0 − r) −iω/2πT (4.62) at the horizon. By taking appropriate linear combinations of these solutions we can construct two new solutions with either G (0) Performing a fit to the nearboundary behaviour (4.59), we can then determine G (4) xz and B (2) yz for these solutions, allowing us to evaluate the derivatives with respect to the sources appearing in equation (4.61).
For simplicity we will restrict to the case M = 0, deep in the nodal line phase. Figure 11 shows our numerical results for the three shear viscosities η xz,xz , η zx,zx , and η xz,zx , each multiplied by 4π/s, at M = 0 and as a function of T /b. The solid black curve shows η xz,xz , which at large T /b is given approximately by the universal result for isotropic holographic systems, η xz,xz s/4π. Moving to small T /b we find that η xz,xz /s decreases monotonically with decreasing T /b, becoming very small as T /b → 0. Unfortunately our numerics for η/s become unstable at small T /b, so it is not possible to determine whether η xz,xz /s → 0 exactly as T /b → 0.
The dot-dashed orange curve shows η zx,zx , which displays the opposite behaviour, increasing with decreasing T /b, becoming very large (and possibly diverging) as we send T /b → 0. Finally, the dashed blue curve shows η xz,zx . We find that this component of the shear viscosity is numerically very close to s/4π for all T /b. The figure shows η xz,zx deviating slightly from this result at very small T /b, but we cannot rule out the possibility that this deviation is caused by the instability of our numerics at low temperatures.
If η xz,zx really does equal s/4π, than one would expect to be able to prove this somehow. To do so, one would presumably need to find a formula for η xz,zx in terms of quantities evaluated at the horizon at r = r 0 . Following ref. [78], one might try to do so by writing the equations of motion (4.57) in the form of a matrix equation for the vector of fluctuations (G xz , B yz ). If one of the eigenvalues of the "mass matrix" M appearing in this equation vanishes, then there is a linear combination of the fluctuations that is independent of r in the limit ω → 0, and we would expect it to be this combination that determined η xz,zx . However, we find det M = 0, so no such combination exists.

Discussion
We have studied a variety of properties of the holographic NLSM model [16][17][18], modified by an additional coupling that gives us greater control over the IR physics. We have paid particular attention to transport phenomena. The two-form coupling b xy , responsible for the presence of nodal lines in the fermion spectral function, breaks rotational symmetry, with important consequences for transport. One obvious consequence of the broken rotational invariance is that the coefficients describing transport in different directions are typically different. There is also the subtler effect that broken rotational invariance implies that the stress tensor is not symmetric in its indices, leading to a larger number of transport coefficients than in a rotationally invariant system. There are many possible directions for future work on holographic NLSMs. We discuss a few examples below.
Throughout this paper we have worked at zero chemical potential, i.e. with a Fermi energy equal to the energy of the nodal line, and considering also non-zero chemical potentials would clearly be of interest. Related to this is the implementation of an energy tilt of the nodal line. Moreover, experiments with nodal line semimetals typically consider quantum oscillations. To discuss these in a holographic setting an external magnetic field should be added to the problem as well. Of course we should realize that the nodal lines in the present model take the form of circles in momentum space, due to the unbroken SO(2) rotational symmetry in the (x, y) plane. Nodal lines in real materials are typically less symmetric, for instance, the nodal line in Ca 3 P 2 has a six-fold discrete rotational symmetry [1]. A more realistic holographic NLSM model should therefore include terms that break the SO(2) to some discrete subgroup, while preserving the discrete symmetries that protect the nodal line.
As discussed in the introduction, in order to obtain more realistic fermion physics, it may be better to work in the framework of semiholography [30][31][32][33], supplementing the fermion action in equation (3.2) with additional boundary terms such that the boundary fermion is elementary, rather than composite. In particular, this means that the fermion spectral function ρ would then satisfy the ARPES sum rule ∞ −∞ dω ρ = 1, satisfied by electrons in real materials. A calculation of the fermion contribution to the electrical conductivity of a semi-holographic Weyl semimetal was performed in ref. [20], and one could attempt to generalise this approach to NLSMs.
By introducing the additional coupling λ into the holographic model (2.1), we were able to obtain a finite, non-zero DC conductivity σ DC zz in the direction orthogonal to the plane of the nodal line at zero temperature, in agreement with the expected physics of NLSMs. However, we always have σ DC xx = σ DC yy = 0, whereas these conductivities are also expected to be non-zero in a real NLSM. It would be interesting to explore ways to obtain finite, non-zero values for all three DC conductivities. Simple dimensional analysis shows that in a theory that has both hyperscaling and an anisotropic scale invariance, such that time scales as t → λt while the spatial directions scale as x i → λ d i x i , the longitudinal conductivity in the x i direction scales as 16 Since the frequency scales inversely to time, ω → λ −1 ω, this dimensional analysis implies that the conductivity depends on frequency as σ ii ∝ ω −∆ i at small ω. For our system, the IR fixed point of the topological phase has just such a scale invariance, with d x = d y = 1/2 while d z = 1, so this dimensional analysis reproduces the low-frequency results σ xx , σ yy ∝ ω and σ zz ∝ ω 0 that we found holographically. To obtain finite, non-zero values of σ DC ii = σ ii | ω=0 for all three directions, we would require ∆ i = 0 for all i, which from equation (5.1) only occurs if d i = 0 for all i. This situation would correspond to fermion self-energies depending only on frequency, which is not so uncommon in the literature of strongly correlated electrons studied by dynamical mean-field theory approximation [79]. In the gravity dual this would arise from an AdS 2 near horizon region that occurs in extremal black-holes, see e.g. ref. [80] and would correspond to a situation where the IR of the field theory is governed by a one-dimensional CFT. An alternative resolution involves hyperscaling violation. Hyperscaling violation modifies the frequency dependence of the conductivity [81][82][83], so may provide a way to avoid having to set d i = 0 for all i.
The bulk fermions discussed in section 3 were treated purely classically. However, important physics arises from quantum effects involving the fermions, see ref. [69] for a review. For NLSMs, possibly the most relevant quantum effect to study would be the fermion contribution to the electrical conductivity, which arises from a loop correction to the propagator of the bulk gauge field A m [84][85][86].
NLSMs with boundaries exhibit surface states [87], as a consequence of the non-trivial topology of the material's band structure, and one could look for similar surface states in the holographic model. A natural approach would be that of ref. [88], which found evidence for surface states in a holographic model of a Weyl semimetal with a boundary in the form of an edge current. Alternatively, in a semiholographic approach the surface states could presumably be directly obtained by considering a half-infinite system and solving the appropriate one-dimensional Schrödinger problem, for instance, with hard-wall boundary conditions.
In section 4 we computed various transport coefficients of the holographic NLSM, using their Kubo formulas. Due to the asymmetry of the stress tensor, arising from the explicitly broken Lorentz invariance, there are a larger number of coefficients governing transport of energy and momentum compared to in a Lorentz invariant system. To better understand the physics of these coefficients, one should work out the full first-order hydrodynamic expansion of the stress tensor, including antisymmetric terms allowed due to non-zero b µν . We leave this for future work. 17 It would also be very interesting to compute entanglement entropy in the holographic NLSM model. For Lorentz invariant renormalisation group (RG) flows, entanglement en- 16 See ref. [31] for a clear discussion of the case when di is the same for all i. The extension to different values of di is straightforward. 17 See refs. [70,89] for related work on hydrodynamics when boosts, but not rotations, are broken.
tropy may be used to define a c-function, a quantity that decreases along the RG flow from the UV to the IR, and therefore provides a measure of the number of degrees of freedom at a given energy scale [90][91][92][93][94]. For the NLSM, Lorentz invariance is explicitly broken by non-zero b. However, there have been proposals, based on evidence from holography, that an entropic c-function also exists in anisotropic systems [95], see also refs. [96][97][98]. One could test this proposal for the holographic NLSM. The entropic c-function has also been proposed as a probe of topological phase transitions, based on the fact that it displays a maximum near the critical point in a holographic model of a Weyl semimetal [99]. One could look for such a maximum in the quantum phase transition exhibited by the holographic NLSM. We thus conclude that many interesting extensions of our present work exist, and we intend to work on some on them in the near future. In particular, we hope to be able to bring the holographic approach for strongly interacting systems closer to the exciting experiments with the recently discovered nodal line semimetals, which indeed appear to show signs of strong correlations due to the reduced screening near the nodal line. However, for a quantitative comparison between theory and experiments these are still early times and much work needs to be done. Nevertheless, we hope that with our present paper we have at least been able to set an additional step towards this ultimate goal.

A Details of holographic renormalisation
In this appendix we give some details on the derivation of the formulas for various physical quantities quoted in section 2.3.
Free energy. The free energy at temperature T is given by F = T I , where I is the Euclidean signature on-shell gravitational action with Euclidean time periodic with period 1/T . Using the equations of motion (2.6) and (2.7), one can show that the bulk part of the Lagrangian in equation (2.1) may be written as a total derivative, f hg , so that the bulk contribution to the on-shell action may be written as a boundary term where in the second line we have substituted the near boundary expansions (2.24), keeping only terms that are non-zero in the limit → 0, and V = d 3 x is the volume of the system. The Gibbons-Hawking term and the counterterms are straightforward to evaluate, yielding Summing these three contributions, we find the free energy F = T (I bulk + I GH + I ct ) is given by where S is the on-shell gravitational action in Lorentzian signature. The change in bulk action under a small change in φ is where EOM denotes an integral over a term proportional to the equations of motion, which vanishes on shell by definition. The variation in the counterterm action is Summming the bulk and counterterm contributions to the variation of the action and inserting the near-boundary expansion (2.24), we find so that the scalar one-point function is Antisymmetric tensor one-point function. The calculation of the one-point function if the antisymmetric tensor operator proceeds similarly. The one-point function is given by The variation of the bulk part of the action under in change of B is while the variation of the counterterms yields Summming the bulk and counterterm contributions to the variation of the action and inserting the near-boundary expansion (2.24), we find so that the non-zero two-form one-point functions are Stress tensor. The naive holographic stress tensor, obtained by differentiation of the on-shell action with respect to the boundary metric, is [55] T µν ≡ − lim As discussed in section 2.4, this tensor is not conserved in the presence of a non-zero source b µν . The conserved stress tensor is instead T µν ≡T µν + 2O µρ b ν ρ . Inserting the near-boundary expansions, we find T µν = diag(ε, p, p, p), where the energy density and pressure are given by respectively.

B Holographic derivation of the Ward identity for translations
In this appendix we demonstrate that our holographic stress tensor satisfies the expected Ward identity for translations (2.43), via an explicit computation of the stress tensor for the case when the antisymmetric tensor operator has a position-dependent source. For simplicity we will neglect the bulk scalar field, since its contribution to the Ward identity is well known [55], and the expressions in this section will be rather lengthy even in its absence. For the latter reason we will also neglect the self-coupling of the two-form field, i.e. we set λ = 0, since this coupling cannot affect the Ward identity. Throughout this appendix we will work in units in which the AdS radius is L = 1. We write the metric of asymptotically locally AdS 5 as where r is the radial coordinate, with the boundary at r → ∞, and x µ are the field theory directions. The radial coordinate is related to the one used in the main text by r = − log r. The matrix γ µν is the induced metric on constant-r slices. For the remainder of this appendix we will drop the prime on r for notational simplicity.
In the metric decomposition (B.1), the bulk Einstein equations may be written as (see for example equation (53) of ref. [100]) where K µν = 1 2 ∂ r γ µν is the extrinsic curvature of the constant-r slices, K = γ µν K µν is the mean curvature, R µν [γ] is the Ricci tensor computed with γ µν , ∇ µ is the covariant derivative with respect to γ µν , Θ mn defined in equation (2.3) is proportional to the bulk stress tensor, and Greek indices have been raised with γ µν ≡ (γ −1 ) µν ..

C Derivation of the Ward identities for two-point functions
In this appendix we derive the Ward identities for two-point functions quoted in equation (4.25), that arise due to the translational symmetry of our system. The discussion follows ref. [72], but with the addition of the antisymmetric tensor operator O µν sourced by b µν . We will first work out the Ward identities for two-point functions in Euclidean signature, and then perform a Wick rotation. As in section 2.4, we will compute the Ward identities by working on a curved spacetime of metric g with general two-form source b, setting g to be flat and b to be constant at the end of calculations. In Euclidean signature the one-point functions are obtained from functional derivatives of W = log Z as T µν = 2 √ g δW δg µν , O µν = 1 √ g δW δb µν . (C.1) The corresponding two-point functions are T µν (x)T ρσ (y) = 4 g(x)g(y) δ 2 W δg µν (x)δg ρσ (y) , T µν (x)O ρσ (y) = 2 g(x)g(y) δ 2 W δg µν (x)δb ρσ (y) , O µν (x)O ρσ (y) = 1 g(x)g(y) δ 2 W δb µν (x)δb ρσ (y) .
Note that the factors of 1/ √ g in equation (C.1) imply that the functional derivative of a one-point function with respect to g ρσ differs from the corresponding two-point function by a contact term, for instance Now consider the one-point function Ward identity (2.43). In a curved spacetime it becomes Taking a functional derivative of this identity with respect to g ρσ , and then setting g µν = δ µν and b µν to be constant, we find ∂ µ T µν (x)T ρσ (0) + T (µρ) δ νσ ∂ µ δ (4) (x) + T (µσ) δ νρ ∂ µ δ (4) (x) − T (ρσ) ∂ ν δ (4) (x) where (·) denotes symmetrisation of indices, with normalisation T (µν) = 1 2 (T µν + T νµ ). The various contact terms in equation (C.5) arise either from equation (C.3) or from the variation of the Christoffel symbols hidden in the covariant derivatives in equation (C.4). We now perform a Fourier transform, defining the momentum-space two-point function where k E = (ω E , k) is the Euclidean four-momentum. From equation (C.5), we find that the momentum space two-point function satisfies where [·] denotes antisymmetrisation of indices, with normalisation T [µν] = 1 2 (T µν − T νµ ). Now we use the relation T µν =T µν + 2O µρ b ν ρ to replaceT with the stress tensor, finding However, it is straightforward to verify that T (µρ) − 2b [µ λ O ρ]λ = T µρ identically, so equation (C.8) simplifies dramatically to A second Ward identity may be obtained by instead differentiating equation (C.4) with respect to b ρσ . After setting g µν = δ µν and taking b to be constant, we find (C.10) Performing a Fourier transformation, we then find that the momentum-space two-point function satisfies We now perform a Wick rotation to Lorentzian signature, assuming that under the Wick rotation AB (k E ) → − AB L (k) for any operators A and B, where k µ = (ω, k) is the Lorentzian four-momentum, and · L denotes some Lorentzian signature two-point function. 18 For the moment we will be agnostic about precisely which Lorentzian signature two-point function (e.g. retarded, advanced, or time-ordered) we obtain through this procedure. Under the Wick rotation, each upper time-like index picks up a factor of i, while each lower time-like index picks up a factor of −i. See for example ref. [102] for an explicit accounting of all of the factors of i. The result is that we find that the Lorentzian signature two-point functions satsify k µ T µν T ρσ L (ω, k) − η νρ T µσ − η νσ T µρ + η µν T ρσ = 0, These are the Ward identities written in equation (4.25), but with · L in place of · R . Now we turn to the question of which Lorentzian correlation functions satisfy the Ward identities (C.12). It turns out that · L cannot be the retarded Green's function, since the components T 0i T ρσ L (ω, k) have the wrong ω → 0 limit [72]. However, one may choose · L to differ from the retarded Green's function only by a constant contact term [72,73], which cancels out in the difference between two-point functions taken in the Kubo formula (4.24). We will denote by · R the set of two-point functions related to the retarded Green's functions only by contact terms. They satisfy the Ward identities (4.25).