Bit threads, Einstein’s equations and bulk locality

In the context of holography, entanglement entropy can be studied either by i) extremal surfaces or ii) bit threads, i.e., divergenceless vector fields with a norm bound set by the Planck length. In this paper we develop a new method for metric reconstruction based on the latter approach and show the advantages over existing ones. We start by studying general linear perturbations around the vacuum state. Generic thread configurations turn out to encode the information about the metric in a highly nonlocal way, however, we show that for boundary regions with a local modular Hamiltonian there is always a canonical choice for the perturbed thread configurations that exploits bulk locality. To do so, we express the bit thread formalism in terms of differential forms so that it becomes manifestly background independent. We show that the Iyer-Wald formalism provides a natural candidate for a canonical local perturbation, which can be used to recast the problem of metric reconstruction in terms of the inversion of a particular linear differential operator. We examine in detail the inversion problem for the case of spherical regions and give explicit expressions for the inverse operator in this case. Going beyond linear order, we argue that the operator that must be inverted naturally increases in order. However, the inversion can be done recursively at different orders in the perturbation. Finally, we comment on an alternative way of reconstructing the metric non-perturbatively by phrasing the inversion problem as a particular optimization problem.


Motivation
Recent progress in the joint program on quantum information and holography has uncovered striking connections between entanglement and spacetime. Arguably, the most exciting discovery in this context, and the one which ignited most of the research in this field, was the proposal of Ryu and Takayanagi that relates the entanglement entropy of a region A in the boundary to the area of a minimal codimension-two bulk surface γ A [1], (1.1) This formula was further generalized to a fully covariant setting in [2] and proved formally in [3,4] using the known entries of the holographic dictionary. The RT prescription (1.1) and its covariant version generalize in an elegant way the well-known Bekenstein-Hawking formula for black hole entropy and provide a natural way to interpret it directly in terms of a microscopic CFT description. Given its elegance and simplicity, entanglement entropy became a robust tool to investigate fundamental aspects in holography, ranging from the problem of bulk reconstruction [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21], to the emergence and dynamics of spacetime [22][23][24][25][26][27][28][29].
Recently, Freedman and Headrick proposed an alternative way to compute entanglement entropy that does not rely on bulk surfaces, but instead, is phrased in terms of a specific flow maximization problem [30]. More specifically, the new prescription states that and can be shown to be equivalent to the RT formula through the continuous version of the max flow-min cut theorem of network theory. The maximization above is an example of a convex optimization program and, hence, the equivalence between (1.1) and (1.2) can also be proved using techniques borrowed from convex optimization [31]. Soon after this paper appeared, it was realized that various geometric problems could likewise be translated to the realm of convex optimization leading to interesting new results [32,33]. The connection with convex optimization has also helped uncover various properties of entanglement entropy from the bit thread perspective [34], as well as some generalizations and applications to other entanglement related quantities [35][36][37][38][39][40][41]. A complementary approach that departs from the realm of convex optimization was put forward in [42,43] and studies aspects of bit threads and entanglement by considering explicit constructions of max flows. This is the line of work that we will mostly follow in this paper. There is one crucial distinction between the two prescriptions to compute entanglement entropy that we believe deserves further investigation: while the minimal surface γ A is in most cases unique, the solution to the max flow problem v is highly degenerate. More specifically, it can be shown that v is uniquely determined only at the bulk bottle-neck γ A , but is highly non-unique away from it. This non-uniqueness raises the question: Out of the infinitely many thread configurations that could be associated with a boundary region, is there any meaningful separation or classification that could be associated with states of special "entanglement classes" in the dual field theory?

JHEP01(2021)193
Intuitively, it would seem that this large degeneracy could indeed be associated to a choice of microstate (or a particular class of microstates) that give rise to the same amount of entanglement between the region A and its complement, 1 however, a precise version of this statement is not settled yet. On the other hand, one can try to exploit this nonuniqueness to gain new insights on the gravity side. The utility of the non-uniqueness property stems from the observation that, if a version of this statement is true (even if we do not know it yet), then a particular solution to the max flow problem v could potentially carry more information than the minimal surface γ A itself: it could encode in detail how the local correlations between the degrees of freedom in the region A and in its complement are distributed for a particular choice of microstate. If so, then, one could imagine that specific questions related to bulk reconstruction and the emergence of spacetime could be answered in a more efficient way by properly selecting a class of configurations/states adapted to the specific problem at hand.
In this paper we will give some steps in this direction. Specifically, our main objective is to understand how the program of gravitation from entanglement [22][23][24][25][26][27][28][29] unfolds in the language of bit threads and to explore an alternative way of metric reconstruction based on this framework. The particular questions that we want to address are the following: • How are the metric and Einstein's equations encoded in generic thread configurations?
• Can bulk locality be manifest in particular constructions?
• Is it possible to reconstruct the bulk geometry from a max flow solution?
• If so, how does the method compare to the ones based on RT surfaces?
Following [22][23][24][25][26][27][28][29], we begin by considering these questions in a perturbative setting in which we study small deformations continuously connected to a reference state. An important motivation of such continuous construction comes from the study of the phase transition of RT surfaces that happens for disjoint regions as one varies their separation. It is known that, close to the phase transition, the RT surface can change from a connected to a disconnected configuration. Such jumps posit a puzzle to a potential quantum information interpretation of the RT surfaces from the bulk perspective, which is solved in the language of bit threads by imposing the additional property of being continuous across phase transitions [30]. Continuity is, then, a desirable feature of bit threads under continuous deformations.
Before we study the above questions, let us review some of the features of the standard methods of metric reconstruction using RT surfaces [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21], and explain potential advantages of studying this problem with bit threads. While there are other methods for bulk reconstruction, e.g. [45][46][47], our comments and comparisons refer only to approaches that make explicit use of RT surfaces. Quite generally, if one hopes for a reconstruction of JHEP01(2021)193 the metric everywhere in the bulk one must start with a sufficiently dense set of extremal surfaces that probe the full manifold M . This is in fact possible in some simple cases, at least for the subset of M that can be foliated by boundary anchored extremal surfaces. For static (2 + 1)−dimensional bulk geometries this was achieved in [9] starting from the full set of extremal surfaces associated with all CFT intervals, and using ideas from holeography [5][6][7][8]. More recently, it was shown that the same ideas could be extended to the time-dependent case in [18] and to higher dimensions [19], here by focusing on the subset of extremal surfaces associated with spherical regions (or topologically equivalent, in the approach of [20]).
The problem of metric reconstruction using bit threads has a major advantage over the ones described above: it does not rely on the ability of the manifold M of admitting foliations by boundary anchored extremal surfaces. In fact, threads can probe regions in the bulk that extremal surfaces cannot, such as entanglement shadows near the vicinity of (spherical) black hole horizons [43]. It is important to point out that bulk shadows do not appear exclusively in cases where gravity is strong; one simple counterexample is the metric of a conical deficit geometry, which arises by the backreaction of a point particle in AdS [48]. Consequently, formulating the problem of metric reconstruction in the language of bit threads, even for the simpler case of perturbative states, is interesting on its own right. In particular, it will shed new light on the issue of emergence of spacetime from entanglement entropy [49,50], without resorting to other measures of entanglement such as entwinement [48].
Another important difference with respect to the problem of metric reconstruction using extremal surfaces is that the latter requires as a starting point the knowledge of a dense set of surfaces that probe the bulk geometry. While we can do the same in the language of bit threads, i.e., start from a dense set of thread configurations, the fact that one single solution to the max flow problem already probes the full bulk geometry presents us with an interesting possibility: we can start from a finite set of thread configurations, containing one, or possibly only a few solutions of the max flow problem. We will consider both approaches in this paper, and show that the explicit reconstruction is possible in both cases. In the remaining part of the introduction we will provide a quick guide to help navigate our paper and enumerate the most important findings of each section.

Road map and summary
We begin in section 2 with a short discussion of various topics that we constantly refer to in our paper. Most of this material is a review of previous works, covering known results about perturbations around AdS and the calculation of entanglement entropy, both in the language of extremal surfaces and bit threads. We also include a short analysis of bit threads in perturbative excited states in subsection 2.3.1 which is new. The main message of this analysis is that, to leading order in the perturbation, it is consistent to use the prescription (1.2) on a constant-t slice, even if the perturbation includes time dependence.
In section 3 we study simple explicit realizations of max flows for bulk geometries that are perturbatively close to pure AdS. We begin in section 3.1 by discussing some general properties about these max flows: the boundary condition at the minimal surface and JHEP01(2021)193 how this condition is sufficient to encode the first law of entanglement entropy. We then proceed to study two particular constructions in subsections 3.2 and 3.3, respectively. The first method that we consider is a generalization of the geodesic method developed in [43]. This method assumes a particular set of integral curves as a starting point, which we take to be the family of space-like geodesics that intersect normally the minimal surface. Given this assumption, one then determines the norm by imposing the divergenceless condition, implemented through the Gauss's law. We show that this construction works both for geodesics of the unperturbed and perturbed geometries under some mild assumptions. In subsection 3.3 we study a slightly more general method. Here, our starting point is to propose a family of level set functions for the flow and then determine its norm based on the divergenceless condition, now implemented directly by solving a differential equation. The flows constructed via this method are a generalization of the maximally packed flows presented in [42,43], where the level set functions are now arbitrary (not necessarily a nested set of minimal surfaces). This method is therefore fully non-perturbative and easily adapted to any boundary entangling region.
Importantly, both constructions presented in section 3 assume as an input a solution to the Einstein's equations in the bulk. Given an explicit metric one can determine the norm of the vector field from the divergenceless condition, which requires an integration from the minimal surface (where the norm is known) to the points of interest. Such integration generically introduces a nonlocal dependence on the background metric which renders these methods non-suitable for addressing questions of bulk reconstruction. However, this also suggests a way around it. More specifically, since the nonlocality is introduced in both cases through the implementation of the divergenceless constrain, it suggests that a construction that implements this condition in a background independent way would be absent of such nonlocalities, which is possible if pose the question in the language of differential forms.
Motivated by the above observation we start subsection 4.1 by rewriting the bit thread framework in the language of differential forms. We study in detail the case of perturbative states and show, in subsections 4.2 and 4.3, that the Iyer-Wald formalism provides a candidate for a thread perturbation which is explicitly local in the metric and furthermore connects the closedness condition with the linearized Einstein's equations. Further, we explore the problem of metric reconstruction in subsection 4.4 and show that it can be cast in terms of the inversion of a particular differential operator. We provide explicit inversion formulas for the case of spherical entangling regions in two distinct scenarios: i) assuming knowledge of a dense set of forms parametrized by their radii and centers and ii) assuming knowledge of a finite set with a minimal number of forms. The second approach turns out to be very powerful; for instance, it suffices to have a single form to provide a full solution for the bulk metric in asymptotically AdS 3 and AdS 4 spaces, which we construct explicitly. We also show that the problem is well-posed in higher dimensions, starting with a carefully selected finite set. We end the section with a detailed analysis of how to recover the time components of the metric via boosts and translations of the space-like hypersurface on which the threads are defined, and a thorough discussion on how to generalize the reconstruction problem to higher orders in the perturbation.

JHEP01(2021)193 2 Preliminaries
In this section, we will start with a brief discussion of a number of topics that we will be essential throughout the paper. We include this discussion for completeness. However, since most of this material is a review, it can be safely skipped by the cognoscenti.

Linear perturbations around AdS
Let us start by reviewing basic properties of linear perturbations around empty AdS. In Fefferman-Graham coordinates, any asymptotically AdS metric can be written as 1) where x σ are boundary coordinates and z is the holographic coordinate. For concreteness we have assumed a Minkowski boundary geometry. With this parametrization, one can extract the expectation value of the stress-energy tensor from the asymptotic form of the perturbation, Plugging the above ansatz into the vacuum Einstein equations, we obtain the following expressions for the zz, zµ and µν components [22]: respectively, where the box operator is the standard Laplace operator in Minkowski space, i.e., ≡ ∂ µ ∂ µ . Alternatively, one can write down the perturbation as follows: where G(x, z) is the Green's function of the graviton in empty AdS, (2.6) A somewhat useful expression can be obtained by expanding δg µν in powers of z [51]: The strategy is to use the linear Einstein equations order by order in z to determine T (n) µν for n > 0 in terms of the expectation value of the stress-energy tensor T µν . A simple calculation shows that the zz and zµ equations imply µν is traceless and conserved for all n. Finally, the µν equations imply that It is convenient to go to momentum space, where ∂ z C µν = 0. Moreover, since H µν must be finite at z = 0, the only possibility is that C µν = 0. This means that only the n = 0 term in (2.7) survives, while all other higher order terms vanish. This can also be seen from the recursive formula (2.9): for d = 2 we have that T (0) µν = 0, therefore all n ≥ 1 terms vanish! The above analysis implies that, to linear order in the perturbation, the general solution for the metric in d = 2 is given by: Since the stress tensor should be traceless and conserved the general form it can take is the sum of right-moving and left-moving waves, Specific examples can be obtained by specifying the profiles of f (t − x) and g(t + x). In appendix A we will explore in detail the case corresponding to a local quench state.

Linear corrections to entanglement entropy
Entanglement entropy can be computed via the RT formula [1],

JHEP01(2021)193
or its covariant HRT version [2], where the minimality condition is replaced by extremality, (2. 16) We are interested in computing the leading correction to entanglement entropy, assuming that the geometry is a small perturbation over AdS. At linear order in the expansion parameter, λ, entanglement entropy can in principle receive two types of contributions. To see this we can expand the area functional L and embedding functions φ(ξ i ) parametrizing the codimension-two surface γ A as follows, (2. 18) This means that at this order, the embedding φ(ξ i ) can be taken to be the (unperturbed) embedding in pure AdS. This is a useful property, because there are many exact solutions for the embending functions of various regions in empty AdS. For our purposes, it will suffice to recall the explicit embedding for spheres in empty AdS in Poincaré coordinates, We will make use of this expression in later sections when we discuss concrete realizations of perturbative bit threads.

Bit threads in dynamical scenarios
The original formulation of bit threads [30] is equivalent to the (non-covariant) RT formula [1], equation (2.15), so it only applies to situations with time reflection symmetry (e.g. spatial regions in static spacetimes). In this section we will explain one way to extend this prescription to fully dynamical cases and show that the formulation of [30] extends straightforwardly to the case of perturbative excited states. One way to include time dependence is by using the maximin reformulation of HRT [52]. To do so, we pick a particular Cauchy surface Σ that contains the boundary of the region, ∂A, perform the area minimization on it, and then maximize over all possible Σ. We can then use the standard bit thread prescription for each Cauchy surface Σ by maximizing the flux through the boundary region Σ ∩ D[A], 2 and then maximizing over all Cauchy JHEP01(2021)193 surfaces: Here, X(Σ) is the space of vector fields on Σ. We note that this formula was recently studied in the context of the membrane theory [40]. There also exists a fully covariant bit thread version of the correspondence [53], but we will not use it in this paper. A solution to the maximin prescription given by the left-hand side of (2.20) consists of a codimension-two surface γ A that solves the two optimization steps. Such a solution would naturally be accompanied by a specific choice of a codimension-one hypersurface Σ on which γ A is a minimal surface. However, in [52] it was shown that such Σ is highly non-unique away from the maximin surface γ A . This fact was used in [52] to argue that one could pick a particular Σ that simultaneously contains the maximin surfaces of various disjoint boundary regions required to prove the strong subadditivity property of holographic entanglement entropy. Below, we will use this freedom to argue that to first order in a general time-dependent perturbation of a static metric, one can always choose Σ to be the constant-t hypersurface associated with the unperturbed metric Σ 0 , or more in general, any space-like surface that is perturbatively closed to it and passes through γ A , Σ λ .

The case of perturbative excited states
Even though the choice of Σ is highly non-unique, it can be shown that not any slice that passes through γ A is a good one. The reason is that γ A is not necesarily minimal on any of such slices Σ. To see this, consider a null congruence shot out from γ A . The surface γ A is extremal, hence, its expansion vanishes: θ = 0. However, the Raychaudhuri equation implies that dθ/dλ < 0 [52]. This means that in this case γ A is a local maximum of area rather than a minimum and, by continuity, the same should hold for space-like surfaces Σ that are close enough to the null congruence. In the left panel of figure 1 we give an example to illustarate this fact. Notice that in one of these slices Σ the minimal area surfaceγ A is not the same as extremal surface γ A . Therefore, finding a max flow in Σ is not equivalent to computing the entanglement entropy of region A.
For the case of perturbative excited states, a natural candidate for a good Cauchy slice would be a slice Σ λ that is perturbatively close to the t = t 0 hypersurface associated with the unperturbed metric, Σ 0 . We can parametrize such a slice as t = t 0 + λ δt(z, x), with the constraint that δt(z, x) must vanish at γ A . The question here is if we can find surfaces on Σ λ that are homologous to A but have smaller area than γ A at order λ. Supposing there are such surfaces, we denoteγ λ A as the one with the minimal area. However, we know that γ A is a minimal area surface in the unperturbed background, therefore, by continuity we know thatγ λ A → γ A as λ → 0. Without loss of generality we can then parametrize such a surface with embedding functions as in (2.17). On the other hand, the calculation in (2.18) shows that corrections to the embedding do not affect the area at linear order. This means that area(γ λ A ) = area(γ A ) + O(λ 2 ), so we can conclude that γ A is a minimal area surface on Figure 1. A solution to the maximin problem γ A is naturally accompanied by a specific choice of a codimension-one slice Σ on which γ A is a minimal area surface. Such a slice is highly non-unique, however, not all slices that pass through γ A are allowed. Left: in this example Σ is perturbatively close to the null congruence shot out from γ A . In this case the minimal area surfaceγ A (orange curve) does not coincide with γ A (red curve). Right: for perturbative excited states, it can be shown that γ A is a minimal area surface on any slice Σ λ that is perturbatively close to Σ 0 . This means that we can pick any of these surfaces, and in particular Σ 0 , to construct relevant bit thread configurations.
any Σ λ perturbatively close to Σ 0 . We illustrate this result in the right panel of figure 1. This also implies that on any of these surfaces, and in particular on Σ 0 , the solution to the max flow ploblem computes the entanglement entropy of region A, and hence all of them are equally good for the construction of bit thread configurations.

Simple realizations of perturbative bit threads
Given the enormous simplification that happens at O(λ) from the point of view of the HRT prescription, we would like to study and understand the general properties of perturbative thread configurations based on the constructions developed in [43]. We will start by stating simple constraints that the O(λ) HRT surfaces induce on general bit threads, and then proceed with the specific constructions. We will show that these methods lead to thread configurations that successfully encode general properties of the CFT state and the bulk geometry, such as the first law of entanglement entropy and its relation to the (linearized) Einstein's equations, albeit in a highly nonlocal form. Along the way, we will state the precise problem of metric reconstruction that we look to solve and enumerate the challenges that these simple constructions face, leading to a quest for a new method that exploits bulk locality in a more explicit way.

Generalities
Let us begin by considering empty AdS d+1 in spherical coordinates. The geometry of a constant-t slice Σ is given by

JHEP01(2021)193
The minimal surface γ A for a ball of radius R is given implicitly by and its outward-pointing unit normal vectorn at a point (r, z) on the minimal surface is: For simplicity we have omitted the angular coordinates, since both the minimal surface and the state are invariant under rotations. A simple realization of a vector field/thread configuration, v = |v|τ , based on geodesics is given by (see [43] for details) As a check, notice that this vector field (i) satisfies the divergenceless condition ∇ · v = 0 and (ii) is equal ton at the location of the minimal surface v| γ A =n. Combining these two, it immediately follows that the flux along any bulk surface Γ A homologous to A (not necessarily the minimal surface γ A ) yields the entanglement entropy of the ball (in units of 4G N ), We emphasize that while the minimal surface γ A is in most cases unique, the choice of vector field v is highly non-unique; it is uniquely determined only at the bottle-neck γ A . Next, we would like to find the perturbed vector field in a perturbatively excited state, i.e., a state with bulk metric g λ µν = g µν + λδg µν + O(λ 2 ) (satisfying Einstein's equations): at linear order in λ. While the perturbation in the vector field δv is on its own highly non-unique, any consistent realization must satisfy some nontrivial properties, including the first law of entanglement entropy in the CFT and the linearized Einstein's equations in the bulk. The problem that we want to address is the following: Given a consistent thread configuration for an excited state v λ , is it possible to reconstruct locally the bulk geometry at the same order in the perturbation?

JHEP01(2021)193
A couple of comments are in order. First, note that we are focusing on excited states. While it is true that the same question would make sense even in the vacuum state, we recall that the bulk metric in this case is fixed by symmetries, rendering the problem exceptionally simple. Second, the non-uniqueness of v λ for a given metric indicates that the correspondence is not one-to-one. Even if we isolate a family of thread configurations that follow from the same bulk metric, the way they encode this information may be nonunique and, generically, highly nonlocal. In the following, we will identify basic constraints that generic realizations of v λ must satisfy and then, study how the particular constructions of [43] encode the information about the bulk metric.

Boundary conditions for the perturbed threads
In order to find a solution v = |v|τ for a thread configuration, we need to solve for the divergenceless condition ∇ · v = 0 subject to the norm bound |v| ≤ 1. One way to proceed is to use the fact that the norm bound is saturated |v| = 1 at the bottle-neck γ A . In other words, we need to impose that at the minimal surface γ A , v is equal to its unit normal, v a λ | γ A =n a . (3.9) Notice that this does not uniquely determine the vector field everywhere in the bulk; intuitively, the ambiguity of the thread configuration away from γ A corresponds to a choice of microstate in the dual CFT, such that all the macroscopic properties of the system are satisfied, including the entanglement entropy S A . Let us now determine how (3.9) looks like in the perturbed geometry. Fortunately, at the linear order in the perturbation the RT surfaces are unchanged and we can use this to our advantage. This implies that at this order, the change in the normal vector is only induced by the change in the geometry. To see this, consider the metric on a constant-t slice Σ of the perturbed geometry 3 where We will keep the λ's explicitly throughout our calculations as a bookkeeping device (to count the order of the perturbations), but at the end we will set it to unity. Also, for future reference, we give an explicit expression for the inverse metric at linear order in λ, where:

JHEP01(2021)193
As explained in the previous section, the embedding function (3.2) is not corrected at this order. Therefore its normal covectorn a remains the same, up to an overall constant N , (3.14) Ensuring thatn is properly normalized to one, we find that at linear order in λ: Finally, raising the index with the inverse metric we find that For example, in d = 2, we find that: For d ≥ 3 we can obtain similar but more longwinded expressions but for the sake of simplicity we will not transcribe them here. Finally, from (3.9) we find that at linear order in λ, our boundary condition at the bottle-neck γ A is: We emphasize that this condition does not uniquely determine v λ in the bulk, specially in regions far away from γ A where v λ is highly non-unique.

First law of entanglement entropy
Since v λ is divergenceless, the flux across any bulk surface homologous to A is constant. Hence, the boundary condition (3.18) should be enough to demonstrate the first law of entanglement, provided we pick γ A itself as our homologous region.
To illustrate this, we can perform a simple analysis in d = 2 dimensions. The area element dS γ A in this case is given by The order O(λ) term gives the change in entanglement entropy,

JHEP01(2021)193
Finally, according to (2.7) we can expand H xx (t, x, z(x)) as However, as emphasized in the previous section, for d = 2 only the n = 0 survives. By the traceless condition we know that T x), so we arrive to the first law of entanglement entropy with the right modular Hamiltonian in 2d [23] For d > 2 the proof is slightly more complicated, but it can be shown by working out the above expansions in momentum space, and resuming the resulting series. We refer the reader to [51] for a detailed analysis in these higher dimensional cases. The crucial insight here is that any divergenceless vector field satisfying (3.18) will automatically encode the first law of entanglement entropy, which for arbitrary dimensions takes the form Since the first law of entanglement entropy has been shown to be equivalent to the bulk Einstein's equations at the linear level [23], then all consistent thread configurations should also encode them in some form. It remains to be seen how are the Einstein's equations encoded in the specific thread configurations, and how easy would be to recover the metric from particular constructions.

Method 1: geodesic bit threads
Following [43], we will now present simple methods to construct explicit thread configurations satisfying the boundary condition (3.18) for perturbative excited states. The first method consists on picking a family of integral curves with good properties, and then fixing the norm by ensuring that Gauss's law is satisfied everywhere. In the following, we will describe this construction in some detail and study how the information of the bulk metric is encoded in the resulting thread configuration.

Integral curves
A good family of integral curves must satisfy the following properties: 1. They must be orthogonal to the minimal surface γ A .
2. They must be continuous and not self-intersecting.
3. They must start and end at the boundary, or possibly at a bulk horizon.

JHEP01(2021)193
Given a family with these properties, it is then straightforward to construct a divergenceless vector field with the desired boundary condition. There is a small caveat here, however: one can only check if the norm bound is satisfied |v| ≤ 1 a posteriori.
One crucial result of [43] is that a thread construction based on space-like geodesics automatically satisfy the norm bound, provided that the metric background satisfies some simple geometric properties. This conclusion followed from a systematic analysis of geodesic foliations of an arbitrary Riemannian geometry, so it must also hold true for the case in consideration, i.e., for geometries dual to perturbative excited states. Therefore, our first candidate for the family of integral curves will be the space-like geodesics of the perturbed background.
Corrected geodesics: let us consider the d = 2 and d > 2 cases separately. In [43] it was shown that space-like geodesics in an arbitrary (2 + 1)-dimensional (d = 2) background lead to a vector field satisfying the norm bound |v| ≤ 1, provided that the Ricci scalar on a constant-t slice (a Riemannian submanifold) is negative everywhere, i.e.
We can check that this condition is indeed satisfied for the perturbative states that we are considering. Working in coordinates adapted to the geodesics, and using the same notation of [43], we will write the bulk metric as follows: where x labels different points along the minimal surface and λ is an affine parameter that runs along geodesics orthogonal to it. 4 The above metric is a solution of Einstein's equations: 5 where T µν is the bulk energy momentum tensor. A quick calculation shows that the induced Ricci on a constant-t slice is: hence, for negative cosmological constant Λ < 0, we have that R < 0 if and only if the local energy density is bounded from above: Since the kind of perturbations that we are considering are all vacuum solutions, i.e. we have T µν = 0, then we conclude that the corrected geodesics can indeed be taken as a good family of integral curves. 4 This coordinate system does not need to foliate the full manifold; points that are not covered by these coordinates have by definition a vanishing vector field v = 0. 5 We have set 8πGN = 1 for simplicity.

JHEP01(2021)193
For spheres in higher dimensional spaces (d > 2) the situation is a bit more complicated. Assuming that the state is invariant under rotations, we can pick a plane that intersects the origin and find the geodesics within this plane. Then we foliate the full spacetime by surfaces of revolution generated by rotating such geodesics along all possible angles. With this construction, the bulk metric can be written as After some algebra, one finds that the criterion (3.25) generalizes to [43] where R 2 is the induced Ricci on the auxiliary 2-dimensional metric defined in (3.30). On a pure AdS background, one finds that R 2 = −2, while the terms on the right hand side of (3.31) are strictly positive. This means that there is a finite gap, or in other words, that the bound is O(1) far from saturation. On the other hand, linear perturbations of the metric would lead to corrections on both sides of the equation but these corrections can only be of order O(λ). This means that for sufficiently small λ, the condition (3.31) will still hold true, regardless of the fluctuations. Similar arguments could be made for metrics that are perturbatively close to AdS but are not rotationally invariant, however, the analysis would be certainly more complicated. In these situations one would need to find corrected geodesics within infinitely many planes intersecting the origin and repeat the above steps. But, again, since the pure AdS case is far from saturating (3.31), the analysis at linear order would only lead to corrections of order O(λ), meaning that the bound would always be satisfied for sufficiently small λ.
The above arguments show that the O(λ) geodesics are good candidates for integral curves for any number of dimensions. There is a slight technical problem, however: it is practically impossible to obtain closed expressions for the corrected geodesics in a generic perturbed background. In practice, rather than working with the corrected geodesics, it is more convenient to propose an alternative family of integral curves. In the following we will explore this possibility in more detail.
Uncorrected geodesics: the corrected geodesics are far from saturating the bound (3.25) in d = 2 or, more generally, (3.31) in higher dimensions. Therefore, it is clear that a continuous family of curves that are perturbatively close to them will similarly do the job. The most natural and simplest candidate for this are the uncorrected space-like geodesics.
To illustrate this point we will consider the d = 2 case, where we can make a precise analytic statement. In this case, the minimal surface (3.2) is given implicitly by We have added subindexes 'm' to point out that these coordinate points are on γ A . The geodesics in pure AdS are given by semicircles anchored at the boundary. These semicircles form a two-parameter family of curves and are defined implicitly by

JHEP01(2021)193
where x s is the center of the circle and R s its radius. The tangent vector with unit norm at an arbitrary point is given bŷ where H(t, x) ≡ H xx (t, x). As expected, the tangent vector still points in the same direction but its normalization is corrected at leading order in the perturbation. Since the integral curves must be orthogonal to the minimal surface, we must enforce thatτ where the latter is given in (3.18). At order O(λ), this requirement leads to 6 In order to arrive to these expressions we have made use of the equation ( Next, we need to check if the proposed integral curves are properly nested [43]. In order to check this, we find the point x a at which they intersect A, 7 37) and the dual point xā at which the curves intersectĀ, One can check that self-intersection is avoided if and only if dx a /dx m > 0 and dxā/dx m < 0. A quick calculation leads to 6 With these definitions Rs can take negative values. We can take an absolute value of xm in the denominator of (3.35) to make Rs positive. However, allowing Rs to take any value will be useful below, in the definitions of xa and xā. 7 If we insist that Rs ≥ 0, these definitions for xa and xā would only be valid for xs ≥ 0, while for xs ≤ 0 one should interchange the two.

JHEP01(2021)193
One can check that at order O(1) both conditions are satisfied, i.e., dx a /dx m > 0 and dxā/dx m < 0. At linear order in the perturbation, we get a term that does not have a definite sign (the last term in the square brackets), but one can always choose a small enough λ such that these inequalities are still satisfied. As an example, let us consider a plane wave, 8 The last term in the square brackets can become order O(1) if the frequency ω is large enough. To prevent this to happen, one must take λ 1/( ωR 3 ). If the background is decomposed in Fourier modes, then the maximum frequency will be the relevant one, and the above condition is replaced by This means that for smooth functions we can always find a small λ that satisfies the conditions. For sharply peaked functions this might not be the case, since the Fourier spectrum could contain arbitrarily high frequency modes. We will therefore restrict our attention to states with smooth stress energy tensor. Notice that this is not an important restriction. In CFT language, a state with a sharply peaked stress energy tensor will not be perturbatively close to the vacuum, and hence, the gravity dual would have important higher order contributions that we have ignored in the approximation of linearized gravity.

Magnitude
Given a set of integral curves the next step is to find the appropriate norm of the vector field |v λ |. We will denote X(x m , ξ) the proposed family of curves; x m labels points on the minimal surface and ξ is a parameter that runs along the curve. As explained above, the curves X(x m , ξ) can be the uncorrected geodesics. The parameter ξ can be taken as the proper length from the given point to the minimal surface. Following [43], we now fix the norm by implementing a version of Gauss's law for an infinitesimal cylinder enclosing each curve. 9 More specifically, we impose that the flux through an infinitesimal area element δA transverse to one of the threads is constant, where h ab = g ab −τ aτb is the projection of the metric on a plane orthogonal to the integral curve. Using the fact that at the minimal surface |v λ (x m , ξ m )| = 1, and letting δA → 0, we arrive to the following expression for the norm In this example even the second term in the square brackets can have a negative sign, but this problem goes away when one impose energy conditions. The last term, however, will still be indefinite after imposing energy conditions. 9 Alternatively, we could fix the norm by solving the first order differential equation for |v λ | resulting from the divergenceless condition, subject to the appropriate boundary condition at γA. This would be completely equivalent to the Gauss's law method described here, since the latter condition is the differential form of Gauss's law. However, since we have the explicit form of the integral lines, the Gauss's law turns out to be more convenient in this case, providing a final answer in closed form, as shown below in equation (3.43).

JHEP01(2021)193
where ξ m is the parameter at which the curve intersects the minimal surface γ A . Notice that we do not need to verify whether the norm bound |v λ | ≤ 1 is satisfied everywhere. This is already guaranteed given our choice of integral curves and the argument based on the negativity of the scalar curvature presented in section 3.2.1. Reference [43] provides various explicit examples of geodesic flows constructed with this method, including the case of spherical regions in empty AdS, given in equation (3.4). In appendix A we complement this study by constructing a new explicit example, now for the case of the specific perturbative excited state corresponding to a local quench.
We can now inquire about how the bulk metric and the Einstein's equations are encoded in this particular construction. Unfortunately, at this level we can already see that such information is encoded in the vector field v λ in a highly nonlocal fashion. On one hand, one needs to solve for the geodesics in the unperturbed background subject to a boundary condition that depends on a particular metric perturbation. And, on the other hand, the magnitude of the vector is found by transporting the boundary condition along the geodesic, ensuring that the vector field is divergenceless. This process is inherently nonlocal; in particular, the final result for |v λ | exhibits an explicit bilocal dependence on the metric perturbation, since it must be evaluated at the points labeled by ξ m and ξ. The latter parameter, in particular, encodes the proper distance between the point in consideration and the minimal surface γ A , which is nonlocal information on its own. These observations imply that it would be rather difficult to invert the problem and recover the metric from the resulting thread configuration. Similarly, the same remarks apply for the Einstein's equations: even though they are assumed as a starting point for this construction (the perturbations we consider are on-shell), they are ultimately encoded nonlocally in the resulting thread configuration.

Method 2: level set construction
The second method of constructing thread configurations consists on starting with a specific family of level set hypersurfaces and then building up a vector field that is orthogonal to them and, of course, divergenceless. This is a slight generalization of a method initially proposed in [43], as we will see below. In the following, we will spell out the details of the general construction for arbitrary metrics, and then specialize to the case of perturbative excited states, where the construction simplifies drastically.

General metrics
We begin by proposing a family of level set surfaces with the following properties: 1. They must contain the minimal surface γ A as one of its members.
2. They must be continuous and not self-intersecting.
3. They must not include closed bulk surfaces.
Given a family with these properties, it is then straightforward to construct a divergenceless vector field with the desired boundary condition. We can understand this as follows: given JHEP01(2021)193 a family of level set hypersurfaces, one can first generate the corresponding integral lines by imposing that they must be orthogonal to each member of the family. Having the integral lines, then, the problem reduces to that of section (3.2) so we could follow the steps outlined there. This means that, in general, we can only check if norm bound is satisfied a posteriori. There is however, one clever exception to the rule. We can ensure that |v| ≤ 1 is satisfied everywhere by construction if we impose the following extra condition on the level set surfaces: 4. They must be homologous to A.
If this condition is satisfied, then, the max flow-min cut theorem guarantees that |v| will be maximal at the bottle neck γ A . Since |v| γ A = 1 then, this implies that |v| ≤ 1 at any other member of the family. Notice that condition 4 is not a strict requirement, but a useful one. In fact, simple examples of vector fields generated by level sets that are not homologous to A are the maximally packed flows constructed in [43]. In that construction the level set surfaces were picked as a family of nested minimal surfaces, containing γ A as one of its members. The motivation there was to find a flow with maximal norm |v| = 1 in a given codimension-one region of the bulk, which was possible due to the nesting property of bit threads [30,31]. 10 For the purposes of this paper, however, we are not interested in the above requierement, so we can explore other possibilities. In the remaining part of this section we will in fact assume that the condition 4 is satisfied, so we do not have to deal with the norm bound.
Let us now describe in detail the construction from level sets. To begin with, we need an efficient way to specify our level set hypersurfaces. In practice, we can do so by picking an appropriate scalar function ϕ(x i ) such that the ϕ = constant surfaces give us our desired level sets. We can then write the following equation At first glance, (3.44) seems more general than a gradient flow, but in fact it is not. In principle one could always redefine the scalar function ϕ →φ = ϕ Υ(ψ)dψ and therefore simply write v = ∇φ. However, the functionφ would not only encode information about the level sets, but also about the norm, so it would be extremely difficult to guess a good function that gives us our desired level sets and that also satisfies the divergenceless condition ∇ 2φ = 0. In practice, then, it is much easier to start with (3.44) and determine Υ(x i ) through the divergenceless condition. We emphasize that, in this scenario, the specific values of ϕ do not have a particular meaning and are in particular not related to the norm of v. The field ϕ here only determines the unit vector in τ = v/|v|, through One crucial observation that follows from the definition (3.44) is that the covector v a (i.e. v with lower index) only depends on the metric g ab through Υ(x i ). To make this point

JHEP01(2021)193
self-evident, and in a form which is partially "independent" of the metric, we can write v a = Υ(ϕ, g)∂ a ϕ . (3.46) We will exploit this observation below, for the case of perturbative states. For now, let us notice that the boundary condition at the minimal surface implies: 47) or, equivalently, All we have left is to determine Υ away from the minimal surface, which can be done by imposing the divergenceless condition. Here we have two options: i) we can use Gauss's law as we did in section 3.2 or ii) we can directly attempt to solve ∇ · v = 0, which should give us a first order differential equation for Υ. As mentioned in section 3.2, the two methods are completely equivalent, since Gauss's law is the integral form of the divergenceless condition. However, since we do not have explicit expressions for the integral lines, then, the first option turns out to be more complicated in this case. 11 We therefore proceed by deriving a differential equation for Υ, which can be derived from ∇ · v = 0. Plugging (3.46) into this condition and massaging the equation leads to: As advertised, this is a first order differential equation for Υ in terms of the scalar field ϕ and the background metric g. Solving this equation subject to the boundary condition (3.48) would then give a unique solution for the vector field v.

Perturbative excited states
The above construction simplifies drastically for the case of perturbative excited states. In the following we will specialize to this situation and study in detail how the information about the bulk perturbation is encoded in the resulting thread configuration. For a metric of the form g λ ab = g ab + λδg ab we are only interested in obtaining the vector field v λ (3.8) to linear order in the perturbation around the zeroth order solution. Since the minimal surface γ A does not change at linear order in λ, a simple choice for the level sets consistent with all requirements would be to pick the same surfaces as for the unperturbed geometry. In this case we have that: In other words, with this choice of level sets, only the function Υ(x i ) gets corrected at linear order in λ, so the first correction of the vector field δv a turns out to be proportional to the zeroth order solution,

JHEP01(2021)193
The function Ψ is determined at the minimal surface by the boundary condition |v λ | = 1.
Expanding at linear order, we obtain Since the zeroth order term is already normalized to one, the terms inside the parenthesis must vanish. Using (3.51) we arrive to: In the last equality we have used the fact that δ(δ a b ) = δ g ab g bc = δg ab g bc + g ab δg bc = 0. We did this, because it would be particularly convenient to have an expression for the boundary condition of Ψ(ϕ, g λ ) in terms of the background v with upper index.
Next we would like to determine the function Ψ away from the minimal surface, which can be done by imposing the divergenceless condition. Again, we proceed by deriving a differential equation for Ψ akin to (3.49). In order to do so, first notice that Taking the divergence of v λ and using the fact that ∇ · v = 0 (at zeroth order), we obtain: Taking the explicit variation of √ g λ and using (3.55) we obtain: where δg ≡ g ab δg ab . In summary, given a background metric g ab and a solution to the max flow problem v a , one can always solve the problem of maximizing the flux in a state where the metric g λ ab is perturbatively closed to the original one. Assuming that the level set surfaces remain the same in the perturbed geometry, the solution for the perturbation of v is given by equation (3.55), which is determined in terms of a scalar function Ψ and the metric perturbation δg ab . This function can be obtained by solving the first order differential equation (3.58) subject to the boundary condition (3.53).
In retrospective, the only non-trivial input required for this kind of construction is the choice of background vector field v, which is in turn used as a seed for the perturbed solution v λ . Specializing to spherical regions, one simple choice would be to pick v as a JHEP01(2021)193 Figure 2. Contour plot for the magnitude |v| of the geodesic flow given in (3.4), in d = 2 dimensions (i.e., empty AdS 3 ). The contours correspond to the level set surfaces of v, which are all homologous to A and, in particular, include γ A as one of its members. This implies that this vector field v can indeed be used as a seed to generate a good solution v λ in a perturbative excited state. geodesic flow, which is known in closed form if the background metric is empty AdS. This background v is given explicitly in equation (3.4). It is easy to check that the level sets of this vector field are all homologous to A, as is shown in figure 2. Since this construction assumes that the level sets are kept fixed, this implies that any perturbative solution v λ build up from this background field v will automatically respect the norm bound |v λ | ≤ 1.
In appendix A we present an explicit example of such perturbative solutions, for the case of a local quench.
Finally, we can comment on how the metric perturbation and the Einstein's equations are encoded in this particular construction. Although the explicit use of metric is reduced in comparison to the construction via integral curves, the last step in the level sets method introduces the same level of nonlocality. In particular, the way we fixed the scalar field Ψ was by solving the divergenceless condition (3.56). Even though this equation is local, the nontrivial boundary condition (3.53) introduces nonlocalities in the solutions, because the equation effectively transports information from γ A to other regions in the bulk. From the Gauss's law perspective the situation is perhaps easier to understand. In that case, the final answer for v λ exhibits an explicit bilocal dependence with respect to the metric perturbation, through its magnitude (3.43). The way we solve for Ψ in this formalism is completely equivalent to that case, because Gauss's law is nothing but the integral form of the divergenceless condition. Hence, even though this construction seems particularly efficient for building up perturbative solutions v λ , it ultimately contains the same kind of nonlocalities than the construction via integral curves. Hence, the inversion problem to recover the bulk metric and the Einstein's equations is equally difficult in both constructions.

Bit threads and bulk locality
The simple perturbative realizations of bit threads of the previous section highlight the need of a bit thread construction that does not make explicit use of the metric. Fortu-

JHEP01(2021)193
nately, we know how to reformulate this formalism in a framework that makes background independence explicit: using the language of differential forms. The equivalence between divergenceless vector fields v and closed (d − 1)−forms w was already emphasized in [30] and was used in [31] to efficiently deal with some subtleties of the max flow problem for null intervals. 13 In this section we will first break down this equivalence in detail, giving explicit formulas that translate various relevant expressions between the two languages. We then argue that the Iyer-Wald formalism provides us with a particular realization of the perturbed thread configuration δw that makes explicit use of bulk locality. In particular, we show that the linearized Einstein's equations are explicitly encoded in this construction through the closedness condition, i.e., dδw = 0. We exploit this unique property of the Iyer-Wald construction to tackle the question of metric reconstruction and show that this problem can be phrased in terms of the inversion of a particular differential operator. Finally, we carry out the explicit inversion at linear order and discuss how to generalize our results to higher orders in the perturbation.

Bit threads in the language of differential forms
In the presence of a metric g ab , the explicit map between flows, i.e., divergenceless vector fields v and closed (d − 1)−forms w, is given by where w represents the Hodge star dual of w, defined via In the above formula ε a 1 ...a d represents the totally antisymmetric Levi-Civita symbol, with sign convention ε i 1 ...i d−1 z = 1. Furthermore, the indices of w are raised with the Riemannian metric g ab , and its determinant is denoted by g. At this point we can already notice an important difference between the two objects, namely that, while the notion of a flow requires a background metric, w can be defined independently of g ab . This will play a crucial role below, specifically, when we address the problem of metric reconstruction. Let us carry on with our analysis. The inverse of the map (4.1) can be stated in terms of the natural volume form , given by where a 1 ...a d is proportional to ε a 1 ...a d and normalized such that i 1 ...i d−1 z = √ g. In terms of , the (d − 1)−form w is given by 4) or in components,

JHEP01(2021)193
Following standard manipulations one can relate the divergence of v a with the exterior derivative of w. 14  This formula shows explicitly the anticipated fact that divergenceless vector fields, or "flows", are mapped to closed (d − 1)−forms. The precise relation between the two is given by (4.4). Now, it is well known that k−forms have well defined integrals over k−dimensional hypersurfaces. Therefore it is convenient to write down an explicit formula for the restriction of w on a codimension-one surface Γ in terms of intrinsic geometric quantities of that surface. Such a formula can be derived using the fact that the volume d−form induces a volume (d − 1)−form˜ on Γ via where n is the unit normal to the surface. Contracting the last index of (4.7) with v and using (4.5) leads to an explicit expression for the form w evaluated at an arbitrary codimension-one surface Γ, with local unit normal n, in terms of the (d − 1)−form˜ This result is equivalently derived in the language of forms, using Stoke's theorem: This leads to which is equivalent to (4.10), given (4.8).
With the ingredients described above, we are now in a position to translate the max flow-min cut theorem to the language of differential forms. First, we have (4.13)

JHEP01(2021)193
The inequality here comes from the standard norm bound, |v| ≤ 1, which in terms of forms can be rewritten as (4.14) In short, equation (4.13) implies that, locally, the form w evaluated on any codimensionone hypersurface is bounded by the natural volume form defined on it˜ . The max flow-min cut theorem then implies that where W is the set of closed forms obeying the bound (4.14). This means that at the bottle-neck γ A , an optimal bit thread form w * should be equal to the volume form˜ , i.e., Finally, combining with the RT formula for entanglement entropy, (4.15) becomes which is the differential form version of the max-flow formula (1.2). There are many situations in which one might want to define the threads in terms of forms w instead of vector fields v. In particular, this reformulation will prove extremely useful for the problem at hand, namely, for the study of perturbations around a given background and the corresponding solutions to the flow maximization problem.

The case of linear perturbations
Having understood how the bit threads formalism translate to in the language of differential forms, it is now time to go back to our original problem. We will assume that the following data is given: a background metric g ab on a manifold M with boundary ∂M , and an optimal flow v that maximizes the flux through a boundary region A. Using (4.5), then, this would also imply the knowledge of an optimal closed form w. In the following, we will consider the max flux problem in geometries that are perturbatively close to g ab , i.e., g λ ab = g ab + λδg ab . We will denote a solution to the problem as w λ , where w λ = w + λδw. First, notice that the closedness condition implies We can also use the fact that the minimal surface γ A does not change at first order in the perturbation, so γ λ A = γ A . Since this is a bottle-neck for the flow, both v and w are fixed at its location. In particular, from (4.16) it follows that Then, given a max flow w for the unperturbed geometry, we are set to find a closed (d − 1)−form δw that satisfies the boundary condition (4.19) and it is such that the norm

JHEP01(2021)193
bound constraint (4.14) holds everywhere in the bulk for the sum w λ = w + λδw. For simplicity, let us introduce the following notation for the inner product 20) and for its first order variation with respect to the metric With these notations, the norm bound (4.14) at first order in λ is given by which looks more difficult to implement than its vector field counterpart. From (4.22) it is clear that the norm bound will typically depend on w so a priori it seems unlikely that a generic δw obeying (4.18) and (4.19) could satisfy (4.22) independent of w. The task becomes even more untractable if one requires δw to be given in terms of a linear local functional of δg ab and its covariant derivatives ∇ (a 1 · · · ∇ an) δg ab (see however [56]). In the remaining part of this section we will show that, despite the above remarks, the Iyer-Wald formalism provides a concrete realization of such perturbed form.

Iyer-Wald formalism and Einstein's equations
One of the crucial breakthroughs in the joint program of holography and quantum information is that the first law of the entanglement entropy, together with the Ryu-Takayanagi formula, imply the linearized Einstein's equations in the bulk. This was originally proven using Hamiltonian perturbation theory [22]. In a beautiful paper [23], it was further shown that it is possible to make this connection more explicit by the proper implementation of the Noether's charge formalism in the bulk, also known as the Iyer-Wald formalism. In this new language, the problem of linearized perturbations is cast in terms of differential forms, a more natural and elegant approach that bridge the CFT and bulk quantities in an efficient way. In this section we will briefly review the basic ingredients of [23], making the connection between entanglement entropy and Einstein's equations manifest. Later in section 4.3 we will show that the Iyer-Wald formalism provides us with a canonical choice for the differential form δw that solves the max flux problem in a perturbed geometry. As a byproduct, we will show that such a canonical form will automatically encode (locally) the linearized Einstein's equations in the bulk which, in turn, will prove useful for the problem of metric reconstruction. Let us first state the problem that [22] sought to solve and then discuss the approach of [23]. In general quantum field theories (holographic or not), for small perturbations over a reference state, ρ = ρ (0) + λδρ, entanglement entropy satisfies the first law For generic CFTs, this setup can be conformally mapped to the case where A is a ball of radius R, centered at an arbitrary point x = x c , in which case [59,60] On the other hand, the left-hand side of (4.23) is computed via the Ryu-Takayanagi formula in the bulk. For ball shaped regions in pure AdS, or small perturbations around it, the RT surface γ A is given by a half hemisphere of radius R extended on the extra dimension, centered at x = x c and z = 0. The Ryu-Takayanagi formula then adopts the form where δ˜ is the variation of the natural volume form on the surface γ A . A further ingredient is the relation between the expectation value of the boundary stress tensor, appearing in the right-hand side of (4.23), and the fluctuations of the bulk metric δg µν . In the Fefferman-Graham gauge, where the latter is given by (2.1), the former can be identified as the first subleading (normalizable) mode in a near boundary expansion (2.2). Taking into account that the boundary field theory is conformal and that the stress tensor conserved, then, this identification imposes non-trivial boundary condition for the metric fluctuations H µν , (4.28) Using the above, the right-hand side of (4.23) becomes Similarly, evaluating the left-hand side of (4.23) using (4.27) yields

JHEP01(2021)193
The first law (4.23), together with (4.29) and (4.30), then establishes a relation between integral functionals of H µν on A and on γ A . It turns out that this functional dependence in turn implies the Einstein's equations for H µν , linearized around pure AdS. This was shown originally in [22] by a direct comparison between the two sides of (4.23). From the gravitational perspective, (4.23) was then proven to be equivalent to the generalized first law of black hole thermodynamics applied to the bifurcate killing horizon of Rindler AdS [23]. This was made explicit by a clever implementation of the Noether's theorem in the bulk, using a formalism developed a couple of decades back by Iyer and Wald in [61]. In order to apply this formalism to the problem at hand, the key observation was that the RT surface for ball-shaped regions in pure AdS, or perturbations around it, coincides with the bifurcate horizon of the time-like killing vector with respect to which a notion of energy and entropy are possible. In fact, a specific conformal transformation (known as the CHM map [60]) maps the interior of the Rindler wedge associated with A to the exterior of an hyperbolic black hole in AdS, where the killing vector ξ coincides with the generator of time translations. Following Iyer and Wald [56,61,62] one then investigates the Noether's theorem for the Killing symmetry generated by ξ. This leads to the definition of a (d − 1)−form with zti 1 ···i d−1 = √ −G. As noted in [23], the form χ satisfies the following properties which can be more easily verified by evaluating (4.32) on a Cauchy hypersurface Σ, containing both γ A and A. For instance, taking Σ to be the t = t 0 slice, one obtains the (d − 1)−form from which both equations in (4.34) follow trivially. The key point here is that χ is closed provided that the bulk equations of motion are satisfied. For instance, in the constant-t Cauchy slice used above one finds

JHEP01(2021)193
where δE g tt is the tt component of the linearized Einstein's equations, and t is the induced volume form on Σ. Similarly, other components of the Einstein's equations are obtained by specializing to different Cauchy slices. Thus, provided that the metric perturbations satisfy these equations, the form χ is closed and the Stokes theorem implies the equality between the left and right equations (4.34). This equivalence also applies in the converse way, given the nonlocal form of (4.34) and the arbitrariness of R and x 0 . This concludes the proof of the statement that they were after, namely, that for theories where the Ryu-Takayanagi formula computes entanglement entropy, the first law of entanglement entropy in the CFT is equivalent to the Einstein's equations in the bulk, linearized over empty AdS.

Method 3: canonical bit threads from Iyer-Wald
Taking into account the nice properties of χ defined via the Iyer-Wald formalism, here, we propose that specializing this form to a spacelike hypersurface Σ containing both γ A and A can be taken as a canonical candidate for the perturbed thread (d − 1)−form 15 Given the integral properties of this form, it is straightforward to check that the flux through any surface homologous to A yields the change of the entanglement entropy in the perturbed state, δS A , as expected. Furthermore, this construction fully exploits the property of bulk locality, in particular, connecting the required closedness of δw with the linearized Einstein's equations via (4.36). We will see below that this property will play a very important role for the problem of bulk reconstruction. It only remains to be checked whether the norm bound constraint (4.14) is satisfied at the desired order (4.22) everywhere in the bulk. This condition will depend on the background form w and then might not hold in general. However, for our purposes it will suffice to find one w such that the combination w λ = w + λδw respects the bound for any perturbation. We will devote the remaining part of this section to check that this is indeed possible.
To begin with, we note that the norm constraint in the form (4.22) is slightly more complicated than its equivalent in terms of vectors. Hence, we will first translate the form (4.37) into the language of flows and then check the condition in terms of the latter. For this purpose, we will need an explicit expression relating δv and δw in the presence of a perturbed metric g λ ab = g ab + λδg ab . In terms of the Levi-Civita symbol ε, the variation of (4.5) reads It is convenient to define a new vector field δv a Φ , given by

JHEP01(2021)193
which is divergenceless with respect to the unperturbed metric g ab , i.e., 40) and is related to δw via its Hodge dual (again, with metric g ab ), The subindex Φ here highlights the fact that the flux of this vector field across any bulk surface homologous to A, computed with the original metric, equals the change in the entanglement entropy δS A . We emphasize that this vector field should not be thought as the variation of the flow v, but just as an auxiliary object. However, given a δv a Φ obtained from (4.41), we can easily recover the true variation of the flow δv a from (4.39). In the Fefferman-Graham gauge (2.1), the metric perturbation takes the form δg ij = z d−2 H ij (with δg zz = δg zi = 0) and δv a reads where H i i = δ ij H ij . Thus, δv a depends not only on δw but also on the background flow v a . In fact, the extra piece in (4.42) is precisely what is needed such that the divergence of v a taken with the full metric g λ ab vanishes at the desired order, Next, we need to make a choice for the background flow v in order to get an explicit δv a and test the norm bound |v λ | ≤ 1. Since the background v should already respect the bound |v| ≤ 1, it is clear that v λ can only exceed this bound by an amount of order O(λ). This can indeed be the case for bulk points that saturate the bound at leading order |v| = 1 (e.g., at the bottle-neck γ A ), or in their vicinity. On the other hand, points that are parametrically far from saturating the bound at leading order are safe, in the sense that we can always take λ to be arbitrarily small such that |v λ | = |v| + O(λ) ≤ 1.
Given the above discussion, then, we should ideally pick a background flow v such that its magnitude decays rapidly away from the minimal surface γ A . Fortunately, we already know good examples of flows respecting this property, e.g., the so-called "geodesic flows" [43], which for spheres in empty AdS take the form (3.4). In the following we will take these background solutions and verify that the norm bound is satisfied at the desired order in the perturbation. First, notice that from (3.4), and using (4.41)-(4.42), it immediately

JHEP01(2021)193
follows that By construction, v λ = v+λδv saturates the norm bound on γ A since the form w λ = w+λδw from which it is derived obeys the appropriate boundary condition for a max flow (4.16). We can check this explicitly: at γ A we have that ξ t = 0, so 16 This leads to the expected saturation at first order in λ, Away from γ A the norm bound is not guaranteed to hold, but since |v| decays as a power law, it would suffice to study |v λ | in a neighborhood of γ A . In order to see this in detail, we note that the level sets of the background flow v (depicted in figure 2) have the form where ∆ ∈ R, i.e., spheres with radius √ R 2 + ∆ 2 , centered at ( x c = x 0 , z c = −∆). It can be checked that, in the vicinity of the minimal surface (∆ 2 R 2 ) Since d > 1, then |v| < 1 for any ∆ = 0 as expected. Now, we want to check whether the norm bound is still satisfied for v λ at the leading order in the perturbation. More precisely, what we really want is to check that for a fixed λ, |v λ | ≤ 1 at linear order in λ for an arbitrary ∆. A short calculation shows that 16 A brief comment is in order. The expression for δva|γ A in (4.46) does not agree with the expected boundary condition at the bulk bottle-neck (3.18). The explanation of this mismatch is simple: the difference between the two vector fields is proportional to a vector that is tangential to the minimal surface δv a T = (δ a b − v a v b ) δv b so its first order contribution to the norm constraint vanishes, g ab v a δv b T = 0, because δvT is orthogonal to v. Therefore, δva|γ A in (4.46) is equally good as (3.18) to our order of approximation.

JHEP01(2021)193
which after some algebra can be put in the following form: The parenthesis in the left-hand side of (4.51), or (4.50), can in fact be non-negative, which implies that the above inequality will not hold for arbitrarily small ∆. Nevertheless, it is interesting to estimate the order of magnitude of the potential violation of the norm bound constraint as it could still be consistent with our order of approximation. From (4.51) it follows that the norm bound can be violated provided that Plugging (4.52) back into (4.50), we observe that this would only lead to a violation of order O(λ 2 ). Since all of our analysis is at linear order in λ, we can safely ignore this issue. In other words, up to our order of approximation the norm bound is not violated and this means that the "canonical" thread configuration constructed from the Iyer-Wald formalism satisfies all the defining properties for a max flow. We relegate to appendix A the study of a explicit example of these canonical thread configurations. Finally, we can comment on how the information of the metric perturbation and the Einstein's equations are encoded in this particular thread construction. First, notice that the variation of the flow δv constructed from Iyer-Wald fully exploits the property of bulk locality. This is evident, since for this particular construction the divergenceless condition, ∇ λ · v λ , maps directly to the Einstein's equations, which are defined locally in the bulk. Moreover, the fact that the δv constructed here (4.45) can be written in terms of a linear local functional of δg ab and its derivatives, present us with an interesting possibility: we can use the information of this canonical solution to invert the problem and recover the bulk metric from it! We will explore this problem in more detail in the next subsection, and comment on the implications and possible generalizations to the full non-linear regime.

Explicit reconstruction at linear order
Our bit thread construction based on differential forms makes explicit use of the property of bulk locality, hence, it should be possible to invert the problem and recover the metric for a generic linear excitation of the boundary quantum state. In this section we will study this problem in detail. More specifically, we will consider a manifold M with boundary ∂M and a set of forms δw that encode the local pattern of entanglement of boundary regions. We will assume the knowledge of the zeroth order -pure AdS-metric g µν , which is otherwise fixed by conformal symmetry (i.e. kinematics), and set up the problem of how to reconstruct the metric perturbations δg µν from the above data.
Our starting point is the knowledge of the change in the entanglement structure of the CFT, which in this case is encoded in the set of (d − 1)−forms δw. We emphasize that these canonical forms can be uniquely specified solely from CFT data. Given a perturbative excited state in the CFT, one can first evaluate the expectation value of the stress-energy

JHEP01(2021)193
tensor T µν and thus the modular Hamiltonian H A associated with a spherical region A. This information can then be used as a boundary condition for δw on A ⊂ ∂M . For instance, specializing to a constant-t slice Σ, this yields 17 where¯ is the natural volume form in the boundary CFT. In fact, we can analytically continue this form to the whole boundary ∂M , so that 18 One way to see that this is consistent would be to conformally map the interior of the sphere to the exterior. Upon implementing this transformation one finds the same functional form for the modular Hamiltonian but integrated along x ∈ A c , hence, providing a boundary condition also at A c = ∂M \A. With the above boundary condition, the full (d−1)−form in the interior of the manifold M is then uniquely determined if we assume bulk locality [56]. To see this, notice that the Iyer-Wald construction yields a form δw such that dδw = 0 on-shell, which is a local condition. If we want to maintain this condition, then, the only ambiguity in δw would be the addition of a term δw → δw + dC where C is a (d − 2)−form such that dC vanishes on ∂M . This is of course a gauge redundancy, which we fix by working in Fefferman-Graham coordinates. Therefore, the boundary condition (4.54) together with the condition of bulk locality are enough to uniquely specify the full (d − 1)−form on M . Before proceeding with the specifics of this analysis, let us first quickly review how the problem of metric reconstruction is normally addressed. In the usual HRT story, given background metric g µν , the change in the entanglement entropy of a region A at first order in the perturbation δg µν is given by 19 This means that δS A encodes information about the first order change in the trace of the induced metric h ij over the extremal surface γ A . On the other hand, the induced 17 Notice that a choice of boundary condition on A is equivalent to picking a specific entanglement contour in the dual CFT. We emphasize that this corresponds to focusing on a particular class of microstates with a given local entanglement pattern. Although this boundary condition is in general non-unique, (4.54) is the boundary condition singled out by the Iyer-Wald construction. 18 More covariantly, on a general Cauchy slice Σ containing ∂A, the boundary condition would be where N µ is a future pointing unit normal vector, and ζA is the conformal killing vector that generates D[A], Upon conformally mapping the causal development of the region D[A] to the hyperbolic cylinder H d−1 ×Rτ , it can be shown that ζA coincides with the time-like Killing vector 2πR ∂τ . 19 The analysis in terms of extremal surfaces can be done for general states, not necessarily perturbative.
Here we are discussing only this simpler case to highlight an important difference with our approach.

JHEP01(2021)193
metric on γ A depends on the bulk metric as well as the explicit embedding of γ A on the geometry. Therefore, by cleverly considering different boundary regions with extremal surfaces intersecting at a bulk point, one could access to the various components of the bulk metric at the given bulk point and hence derive an inversion formula for the metric perturbation δg µν .
As is evident from the previous paragraph, the problem of metric reconstruction by extremal surfaces heavily relies on the possibility of foliating the full manifold M with boundary anchored extremal surfaces. In particular, we would necessarily need to start from a dense family of surfaces that pass through all (reachable) bulk points multiple times. While we can do this in the language of bit threads, i.e., start from a dense set of thread configurations, the fact that one single solution to the max flow problem already probes the full bulk geometry presents us with an interesting possibility: we can start from a finite set of thread configurations, containing one, or possibly only a few solutions of the max flow problem. The minimal number of thread configurations needed in such a set will generally depend on symmetry considerations as well as the number of dimensions, as will be discussed below. For the time being, let us summarize the two approaches to metric reconstruction that we can explore. For simplicity, we will frame the discussion by focusing on a constant-t slice Σ, so that we will aim to recover the spatial components of the metric δg ij . The δg tt and δg ti components can be recovered in a similar way, by choosing appropriate boosted slices Σ , as we will explain at the end of the section. The two methods that we will explore are: • Reconstruction from a dense set of thread configurations. Here we will assume knowledge of δw for all spheres in the CFT, with arbitrary radius R and center point x 0 .
• Reconstruction from a minimal set of thread configurations. Here we will assume knowledge of δw for a few spheres with radius R (i) , and center point x (i) 0 , with i = 1, . . . , n. The precise value of n will be fixed so that the inversion problem is well defined.
We will now discuss these two methods in detail.

Reconstruction from a dense set of thread configurations.
Given the infinite set of (d − 1)−forms δw encoding the canonical entanglement pattern of spheres of arbitrary radius R and center point at x 0 , on a constant-t slice Σ, our goal is to extract the components of the perturbed metric δg ij . We recall that, in the presence of a metric, the set of (d − 1)−forms δw define a set of covector fields δw a (R, x 0 , z, x), instead of the numbers δS(R, x 0 ), so it is clear that in this framework we have infinitely more information in comparison to the standard setup using extremal surfaces. Hence, we can expect to be able to reconstruct the metric in a more straightforward way.
If the full metric is given in the Fefferman-Graham gauge (2.1), the components of the metric perturbation take the form δg ij = z d−2 H ij (with δg zz = δg zi = 0). In this gauge,

JHEP01(2021)193
the components of the covectors δw a can be related to the metric perturbations as follows: These equations can be easily inverted using the dependence on R and x 0 of (4.56) and (4.57). In fact, there are infinitely many ways to invert these equations. The simplest way is to get rid of the derivative terms, so that we obtain a system of algebraic equations. However, we have several ways to accomplish this. Below we will discuss two different options. The first option is by evaluating both sides of (4.56) and (4.57) on the set of parameters (R, x 0 ) that satisfy ξ t (R, x 0 ) = 0, i.e., (4.59) Notice that the requirement given by (4.59) means that our reconstruction is limited to the points that are accessible via extremal surfaces. This means that this option is, in a sense, analogous to the metric reconstruction via the HRT prescription and does not exploit the full reach of the bit threads. We will continue for now and then explain an alternative that does not impose this limitation. Let us first analyze (4.56). From this equation, we can immediately find an algebraic expression that gives the perturbed trace H i i (z, x), (4.60) In fact, we can extract the trace H i i at a point (z, x) from (4.60) using any single covector with parameters (R, x 0 ) such that (4.59) is satisfied. This is an example of the nonuniquess of the inversion formulas. Notice that equation (4.60) provides the solution to the full inversion problem for d = 2, in which case (4.57) is identically zero. In fact, for d = 2 the only component of the metric perturbation that we need to solve for corresponds to H xx (z, x) which equals the trace (4.60). In higher dimensions, we can use (4.57) in addition to (4.56), and proceed in a similar way to extract the information of the individual components of the perturbed metric H i j (z, x). In order to do that, first replace the solution for the trace (4.60) in equation (4.57), so that the latter equation becomes: Further, for a given j and within the set of allowed parameters, we can take x j 0 = x j and x k 0 = x k for k = j. This leads to the following solutions for the diagonal and non-diagonal

JHEP01(2021)193
components of the perturbation: which provides the solution to the full inversion problem for d > 2.
As mentioned above, the method outlined in the previous paragraph is limited to the points that are accessible via extremal surfaces. An alternative way of inverting the system of equations in a less restrictive setting is the following. First notice the relations Using the above, one finds that from which one can easily invert the diagonal and off diagonal components of the perturbation for d > 2. 20 More explicitly, we obtain that which gives the solution to the full inversion problem for d > 2. On the other hand, for d = 2 we can simply use equation (4.60) or, alternatively, the second reconstruction method that we explain below. Notice that summing over j in equation (4.66) leads to a different expression for the trace of the perturbation as (4.60). This represents another example of the non-uniqueness in the inversion equations.
Reconstruction from a minimal set of thread configurations. One may wonder whether the extended nature of the entanglement pattern information present in a single form δw(R, x 0 ; z, x) for fixed (R, x 0 ), or a finite number of forms n, could suffice to recover the components of the metric perturbations δg ij . Naively, the number of unknown variables that we need to solve for is d(d − 1)/2, which is the number of symmetric components of H ij ({i, j} run over d − 1 values). On the other hand, the number of equations that we would have at our disposal is given by nd, i.e., d components of a single covector δw a , times

JHEP01(2021)193
the total number of forms n ∈ N. Assuming they are all linearly independent of each other, we have that for a fixed d, the minimum number of formsn that we would need is where • represents the ceiling function (i.e., the smallest integer greater or equal to its argument). For d = 2 and d = 3 (3d and 4d gravity respectively) we obtainn = 1, which means that we can hope to recover the full metric from a single form. In higher dimensions the problem would be underdetermined so we would need to increase the number of forms, although, to a finite number set by d. In the following, we will focus on the cases for which n = 1 but we will come back to the higher dimensional cases at the end of the section. We start with equation (4.56) and notice that its right-hand side can be rewritten as This equation can be in principle easily inverted for the trace H i i . However, we must proceed with some care because ξ t can attain a zero value in multiple points in the bulk. We refer the reader to appendix B for a detailed analysis of this equation, and we will simply state the final result here. The analysis is naturally split for bulk points x ∈ A c and x ∈ A (∀z). In the former case, ξ t never vanishes in the bulk, so the analysis is particularly simple. In this case we obtain This equation is also valid for x ∈ A provided that z 2 < z 2 * . For z ≥ z * the integrand in (4.70) has a double pole at λ * = z * /z (with λ * ∈ [0, 1]). This divergence can be removed by an appropriate regularization, e.g., using the principle value prescription. However, given the simplicity of our problem we can find the answer directly from (4.69) by taking one of the end points of the integral to be arbitrarily close to the zero locus of ξ t . After a series of careful manipulations explained in appendix B we arrive at a formula that is valid in the region x ∈ A and for all z: (4.71) This new integral seems to still have a single pole at λ = λ * . However, a close inspection shows that this term is proportional to which can be checked from (4.56). Therefore, the integral is manifestly finite. We note that, indeed, the naive principle value regularization of (4.70) results in (4.71), and likewise

JHEP01(2021)193
there are many ways to derive (4.71) from (4.70). Perhaps the simplest way to do it is by slightly changing the integration contour: where ∈ R, and then letting → 0. It can be shown that this prescription is consistent with both (4.70) and (4.71) so it is valid everywhere in the bulk. Notice that in the above formulas we have not specified the number of dimensions d. This means that given the knowledge of a single δw(z, x) one can always find an inversion formula for the trace of the perturbation, given explicitly by (4.70)-(4.71) or its equivalent (4.73). We will now recover the other components of the metric for the cases d = 2 and d = 3.
The d = 2 case: the inversion problem in d = 2 is exceptionally simple because in this case there is only one metric component to solve for. Therefore, the trace of the perturbation provides the full solution to the problem, i.e., H i i (z, x) = H xx (z, x). Nevertheless, we will analyze this case in some detail and check that our formulas for the trace ar consistent with the expected results.
First, notice that in this case there are further simplifications that considerably reduce our problem: equation (4.57) vanishes exactly so δw x = 0 everywhere in the bulk. The closedness relation dδw = 0 then implies ∂ z δw z = 0, so δw z (z, x) = δw z (0, x). Using this fact, and applying the formula (4.73) which is valid everywhere in the bulk, we obtain This is precisely what is expected from the analysis of section 2.1, specifically from equation (2.13). For x ∈ A c , (4.70) coincides explicitly with (4.73) so the integral above is the same. For x ∈ A we can alternatively use (4.71). Notice that for d = 2 the integrand in (4.71) is identically zero since δw z (z, x) is constant so, in this case, the full result is given by the first term in (4.71). Indeed, this term coincides with (4.74), as expected. Again, since for d = 2 this is the only component of the perturbed metric that we need to solve for, then, equation (4.74) completes the inversion problem.
The d = 3 case: to solve the inversion problem in d = 3 we need equations (4.57), in addition to (4.56). These equations involve derivatives with respect to the spatial coordinates ∂ i so, as a system of first order differential equations, we would need information about H ij at some fixed x i in order to have a well defined boundary value problem. We will deal with the choice of such boundary conditions below, but for now, let us explain how to setup the inversion problem.
First, since equation (4.73) already provides the solution for the trace part of the metric, we can already replace this back into (4.57) and solve for the remaining metric components. We find it convenient to write the fluctuations as (4.75) so that h = H i i gives the trace while φ and χ are the two fields that we still need to solve for. Defining x 1 ≡ x and x 2 ≡ y, equations (4.57) then take the form: Thus, we have a system of two coupled partial differential equations of first order. It is possible to decouple these equations by taking one further derivative and combining appropriately the two equations. To do so, we first rewrite (4.76)-(4.77) as: (4.79) Equations (4.78) and (4.79) can now be combined as a pair of Poisson's equations: where ρ φ and ρ χ act as sources for φ and χ, respectively. These equations can be solved using standard Green's function methods. Assuming knowledge of the solutions at a closed surface ∂V, one can formally write the solution in the interior of the surface V as follows: where G( x, x ) is the Green's function for the Laplace operator in 2d, given by Next, we need to impose appropriate boundary conditions. To do so, let us start discussing the region z > R for which ξ t < 0 so that there will be no subtleties with poles in the integrals. We will consider a closed surface ∂V at infinity, so that dS =r r dθ , with r = | x | → ∞. We note that, normally, the integral over ∂V in 2d would give a finite contribution, assuming that the fields that we solve for are finite at infinity. This is because ∇ G( x, x ) · dS → constant at large r . However, in our particular problem, we JHEP01(2021)193 need to solve for the combination φ/ξ t , or χ/ξ t , which decay (at least) as 1/r 2 at large r . 21 This means that the surface integrals in (4.82) and (4.83) vanish, regardless of the specific values of φ and χ at r → ∞, and therefore, where the boundary conditions once again force c 0 = 0. 22 For the region z < R, we note that ξ t can become zero at multiple points in the bulk, so we would need to deal with potential divergences of the integrands. Following the calculation of the trace (explained in detail in appendix B), we can expect to get rid of these non-physical divergences by a simple regularization procedure. One easy way to do this is by adding a small imaginary part to z, (4.88) and at the end of the calculation let → 0. The integration region should now be free of singularities, so these formulas apply for all values of z. Hence, together with the trace formula (4.73), they provide a complete solution of the reconstruction problem in d = 3.
Higher dimensional cases: as discussed above, by comparing the number of independent components of the perturbed metric H ij against the number of equations that we obtain from a single form δw, one can conclude that the minimum number of formsn that are in principle required to invert the problem is given by (4.68). However, this does not imply that any set ofn forms should lead to a well defined inversion problem since, depending on the choice, one could end up with a set of equations that are not completely linearly independent. In the following we will spell out the precise conditions that we must impose on a minimal set of forms and give a concrete example of how these conditions can be satisfied. First, notice that for a spherical region, a single form δw is parametrized by d real numbers P = {R, x 1 0 , · · · , x d−1 0 }. Moreover, one can easily check that for each choice of such numbers, the d components of the form δw a , with a ∈ {z, 1, . . . , d − 1}, are all linearly independent since each component involves different sets of metric components H ij . Therefore, the task at hand reduces to finding a convenient set of parameters P k , with 21 To see this, notice that both φ and χ are components of the perturbation Hij, which can generally be written as (2.7), i.e., Hij ∼ z 2n T (n) ij . The leading order term gives the CFT stress-energy tensor, which scales as T ij , so these terms they decay faster at infinity. Hence, Hij scales like T (0) ij and φ/ξ t ∼ 1/r 2 and χ/ξ t ∼ 1/r 2 at large r (in the worst case scenario). 22 Notice that if c0 = 0, both φ ∼ r 2 → ∞ and χ ∼ r 2 → ∞ at large r, which would be unphysical.

JHEP01(2021)193
k ∈ {1 , . . . ,n}, such that the individual components of each form are linearly independent across the set. More concretely, if we label then forms as δw k , then we would need to impose that for a fixed a the set {δw 1 a , . . . , δwn a } must be linearly independent. In order to visualize explicitly the dependence of each component δw a on the parameters {R, x 1 0 , · · · , x d−1 0 }, we rewrite equations (4.56) and (4.57) as follows: (4.90) In this way, we can express each component of the form as a linear combination of (d + 1) linearly independent functions with coefficients α l given by Note, however, that by choosing a set of parameters P we can only specify up to d of the above coefficients, while one of them will necessarily be determined in terms of the others. This is in fact not a problem. If we consider the set of forms δw k and repeat the above analysis, we find now that where k ∈ {1, . . . ,n}. We note that the number of coefficients α k l that we can fix for each k is larger than the total number forms that we have at our disposal, i.e., d > d− 1 2 . Therefore, we still have a lot of freedom on the choice of parameters P k to be able to make the set {δw 1 a , . . . , δwn a } linearly independent. One way to achieve this is by focusing only on a subspace of forms obtained by varyingn out of the d free coefficients of each δw k a , which we denote as β k l , while keeping the rest fixed. More explicitly, we can split the sum in (4.93) as where β k l is now an ×n matrix, with the coefficients that we vary, andβ k l denote the ones that we keep fixed. The condition for the above linear independence is then given by the non-vanishing of the determinant of the matrix β k l . For instance, we could take

JHEP01(2021)193
and vary the parameters {R k , (x i 0 ) k } for i ∈ {1, . . . ,n − 1} such that det(β k l ) = 0 while keeping (x i 0 ) k fixed for i ∈ {n, . . . , d − 1}. More generally, notice that there can be many other possible choices. One one hand, the choice of the subset β k l ⊂ α k l can be arbitrary and, on the other hand, the remaining free parameters can be random; they must not necesarily be kept fixed.
General time slices: recovering the full bulk metric. Let us now go back to the problem of recovering the full bulk metric δg µν . First, notice that when we picked a constant-t slice Σ, we were able to recover the components of the metric tangent to it, namely, δg ij . This means that we still need to find the time components, δg tt and δg ti , in order to complete the reconstruction problem. Below, we present a simple algorithm to recover these extra metric components.
From the bulk point of view it is easy to see that δg tt and δg ti are, in fact, constrained from the equations of motion (2.4). Specifically, they can be determined from the zz and zµ components of Einstein's equations, These equations imply that δg tt must equal the spatial trace δg tt = δg i i , while δg ti can be determined from a first order equation ∂ t δg ti = ∂ j δg ij . This can be easily implemented in practice. However, the problem is that these particular bulk equations of motion are not known from the CFT perspective, so their use cannot be justified. Indeed, once the surface Σ is chosen to be a constant-t slice, the closedness condition dδw = 0 only encodes the tt component of the Einstein's equations, as explained at the end of section 4.2.
One simple solution to this problem, is to pick a more general time slice Σ and repeat the reconstruction analysis outlined above. For simplicity, we will pick here a boosted slice parametrized by coordinates (t , x i , z ), but a similar analysis can be implemented from more general choices of Σ . We will denote the boosted coordinate x i = x and label the coordinates orthogonal to it as x j = y j for j = i. With this notation, a boost with rapidity β (which we can take to be arbitray) is given by the bulk transformations t = t cosh β + x sinh β , (4.97) x = x cosh β + t sinh β , (4.98) We can now perform the reconstruction analysis in this new slice, and recover the components δg i j on Σ . Indeed, a quick calculation shows that the components that we can recover are δg x x | β = sinh 2 β δg tt + 2 sinh β cosh β δg tx + cosh 2 β δg xx , (4.101) δg x y i | β = sinh β δg ty i + cosh β δg xy i , (4.102)

JHEP01(2021)193
If this analysis is done for at least a few values of the rapidity β 1 = β 2 = 0, then it is clear that we will have enough information to recover the extra components δg tt and δg tx from linear combinations of δg i j | β i . This completes the reconstruction problem for the locus of all points in the intersection between the original slice Σ and the new slices Σ , as shown in left panel of figure 3. In particular, from the transformations (4.97)-(4.100), it follows that the constant-t slices parametrized by β intersect the surface Σ at the line 105) x = x cosh β + t sinh β + σ , (4.106) The new slices Σ , depicted in the right panel of figure 3, intersect the original slice Σ at so it is clear that, if we perform the reconstruction problem for a general slice Σ with arbitrary β and σ we would have enough information to reconstruct the full metric in Σ. Finally, note that the above algorithm has a trivial extension to generic choices of Σ. Thus, picking a family of slices Σ that foliates the full manifold M , and repeating the same analysis for each Σ gives a complete solution to the reconstruction problem in M .

Going beyond linear order
In the past section we have shown that the problem of metric reconstruction can be carried out explicitly at linear order for the case of perturbative excited states. This was accomplished using the canonical bit thread construction based on differential forms. But, can this methodology be generalized to the non-linear regime?
To answer this question, we will start with a brief review of our findings and then discuss how the different aspects of our proposal can be generalized. To begin with the reconstruction problem, we assumed knowledge of a set of canonical forms that codify the entanglement pattern of subregions in the dual CFT. We recall that these canonical forms δw can be uniquely specified from CFT data. Specifically, they are completely fixed by the boundary condition on ∂M (4.54), which is given in terms of the one-point function of the CFT stress-energy tensor T µν , and the requirement of bulk locality. In the presence of a metric, we found that one of these forms specifies a covector field that can be schematically written as ( δw) a = F bc a δg bc , (4.110) where F bc a is a specific linear differential operator. In low-dimensional cases, the equations resulting from a single form provide enough data to invert the problem, so that This kind of inversion works for d = 2 and d = 3, i.e., in AdS 3 and AdS 4 , respectively. For higher dimensional cases the problem becomes underdetermined, but it can be easily generalized by starting from a set of differential forms δW = {δw 1 , . . . , δw n }, that encode the entanglement pattern of a family of subregions in the CFT. In this case, we find that so, for large enough n one can always invert the system as The optimal value of n depends on the number of dimensions, and is given by (4.68). Of course, the larger the value of n, the more information that we have at our disposal, and the easier the inversion problem becomes. In fact, in the limit of n → ∞, or when the set of differential forms is dense enough, there is even enough maneuver to turn the inversion problem into a simple algebraic system of equations as shown explicitly in section (4.4.1).
Let us now comment on what to expect in the non-linear regime. To start with, we can treat the problem perturbatively but extending the above results to higher orders in the JHEP01(2021)193 perturbation. In this case, the perturbation of the metric can be given as an expansion 23 g λ µν = g µν + λδg λ µν , δg λ µν ≡ δ (1) g µν + λδ (2) g µν + λ 2 δ (3) g µν + · · · , (4.114) and, similarly, any solution to the max-flow problem (4.115) On general grounds we expect that, in the presence of the metric (4.114), the change in the form δw λ should follow an equation similar to (4.110), i.e., ( δw λ ) a =F bc a δg λ bc , (4.116) but now withF bc a being a higher order differential operator. For example, at second order in λ we expect a second order differential operator, and similarly for higher orders. This introduces the standard non-uniqueness problem for the inversion of a non-linear operator. However, this issue can be circumvented by solving the reconstruction problem recursively in λ. To see this, notice that the different terms in (4.115) should depend on the different metric perturbations and their derivatives as follows: for j ∈ {1, . . . , i} and k ∈ {1, . . . , 1 + i − j}. In other words, for a given value of i, in the right-hand side of (4.117) we expect up to i th derivatives of δ (1) g µν but only 1 st derivatives of δ (i) g µν . Thus, if we solve for the metric at first order in λ, then we can reformulate the problem of bulk reconstruction at second order as a linear problem above the g µν +λδ (1) g µν solution. This also generalizes to higher orders in λ, so that the inversion problem at a given order can also be cast as a linear problem above one lower order. We should also comment on how the boundary condition (4.54) generalizes to higher orders in λ. At linear level, we saw that it is fully specified by the one-point function of the CFT stress-energy tensor T µν . However, at higher orders we would need to specify further data, e.g. [66]. A useful case study would be to consider the reconstruction problem at second order in λ, which was already worked out in [27] in the framework of extremal surfaces. This paper generalized the Iyer-Wald construction, and thus [23], to include second order variations in the metric, hence their construction can be used to obtain canonical thread configurations at second order in λ. More specifically, [27] focused on a class of CFT states which are expected to have a classical gravity description, and are defined by adding sources for primary operators to the Euclidean path integral defining the vacuum state [67][68][69], (4.118) 23 As shown in [64], when going to higher orders in λ it is convenient to work in the so-called Hollands-Wald gauge [65] so that the coordinate location of γA is fixed and L ξ A (g λ µν )|γ A = 0. In this gauge the argument presented in section 2.3.1 about the choice of Σ can be extended to higher orders in λ.

JHEP01(2021)193
Among other things, their calculation led to a closed expression for the entanglement entropy in these CFT states at second order in the sources: 24 A (x 1 , x 2 ) andK (2) A (x 1 , x 2 ) are some specific kernels. A few comments are in order. First, notice that at this order the entanglement entropy depends on the specific CFT only through the central charges a * and C T . In fact, for theories that are dual to Einstein gravity, they must be proportional to each other, i.e., a * ∝ C T , so the answer depends only on one CFT parameter. Second, (4.119) correctly encodes the first-law of entanglement δS A = H A at linear order in λ but, in addition, it also contains information about the relative entropy in the excited state δ (2) A ), to second order in λ. Finally, note that equation (4.119) can now be used to specify a canonical boundary condition for δw λ in ∂M , and so, it should suffice to uniquely specify the full (d − 1)−form on M by further requiring bulk locality. This means that at second order in λ, we need also the one-point functions of all primary operators O α , in addition to stress-energy tensor. From the bulk perspective, [27] showed that the closedness of the Iyer-Wald form χ, related to δw λ through (4.37), encodes now the following data: at first order one recovers (4.36), which specialized to an arbitrary slice Σ leads to the linearized Einstein's equations E (1) ab = 0. At second order, one obtains where (E (2) ab ) grav is the second order Einstein tensor and T (2) ab is the second order stressenergy tensor associated with the bulk fields φ α dual to the CFT operators O α . Thus, for theories with a * = C T , a CFT state of the form (4.118) secretly encodes Einstein's equations (at least up to second order in the perturbation) with matter determined by the CFT one-point functions. Therefore, on general grounds we can expect that the inversion problem using the canonical bit thread prescription to be well defined at second order in λ.
One can also go to higher orders in perturbation theory, however, it is clear that one would ultimately need infinitely more CFT data in the full non-linear regime, rendering the problem untractable. To see this, notice that to the i th order in the perturbation, the JHEP01(2021)193 relative entropy δ (i) S(ρ A ||ρ (0) A ) will generically involve i-point functions of the primary operators O α , so at the full non-linear level we would need to specify all correlation functions of the theory. Here we have two options to proceed. One could try to pursue the reconstruction problem in a theory with a know gravity dual, e.g., N = 4 SYM, and with some hard work, one should be able to recover the metric order by order in the perturbation. Alternatively, one could start with a generic CFT and try to constrain the structure of the CFT i-point functions such that the calculations match with the gravity side. Indeed, such constraints must start appearing for i ≥ 4, since these correlation functions are not fixed by conformal invariance.
Another possibility would be to consider the full reconstruction problem without resorting to perturbation theory. 25 For example, given a CFT with a known holographic dual, one can start by computing entanglement contours [70,71] for several regions A i in a state of the form (4.118). These contours can then be used as boundary conditions in ∂M for a set of closed bulk forms w i that encode the entanglement pattern of the individual regions. Via the closedness condition, dw i = 0, and assuming further input such as bulk locality, one should then be able to reconstruct particular realizations of the set of forms w i on M that solve the max flow problem in the bulk. If this set is sufficiently dense, then, one could set up the problem of metric reconstruction as a particular optimization problem. More specifically, given a set of such (d − 1)− forms W = {w i }, the metric g ab should emerge as the minimal positive definite symmetric (0, 2) tensor for which the norm bound constraint holds for all the elements w i of the fundamental set W . It would be interesting to develop a more precise algorithm based on this optimization problem and understand how the Einstein's equations would emerge upon its implementation.

Conclusions and outlook
In this paper we developed a new framework for metric reconstruction based on the bit threads reformulation of entanglement entropy. Our work can be divided roughly into two main parts. In section 3 we explored simple constructions of perturbative thread configurations based on the general methods originally developed in [43] but expanded in this paper in various ways. We explored in detail two particular constructions, one that starts by specifying the class of integral curves and a second one that assumes a specific family of level set surfaces. We showed that both methods are efficient and can be easily implemented for the case of perturbative excited states, as we discuss in detail in the concrete example of a local quench in appendix A. However, we realized that both constructions encode the information about the bulk metric in a highly nonlocal way. This implies that these realizations are not particularly useful to tackle the question of metric reconstruction and highlights the necessity of reformulating bit threads in a language that makes background independence manifest.

JHEP01(2021)193
Motivated by the above results, we started section 4 by reformulating the bit threads framework in terms of differential forms. We gave general formulas that translate the relevant equations of the standard description in terms of flows into this new language and studied in detail the case of perturbative excited states. We pointed out that the Iyer-Wald formalism provides us with a canonical choice for the perturbed thread configuration that makes explicit use of bulk locality. More explicitly, we showed that the Iyer-Wald construction yields a particular solution to the max flow problem in the bulk that can be uniquely determined from CFT data, and encodes the Einstein's equations in the bulk through its closedness condition. Assuming that a set of such forms is given, we then showed that the problem of metric reconstruction is equivalent to the inversion of a particular differential operator. We gave explicit inversion formulas for the case of 2d and 3d CFTs, and argued that the problem is also well possed in higher dimensions. Finally, we discussed the generalization of our results to higher orders in the perturbation and its relation to the full Einstein's equations.
There are some open questions related to our work which we think are worth exploring: • Explicit inversion at higher orders: in section 4.4.2 we have sketched out how to generalize the metric inversion problem to higher orders in λ. We believe that this would be fairly straightforward to second order in the perturbation if one uses [27] as a starting point, while it would be illuminating to have explicit inversion formulas for the differential operator at this order. Generalizing this story to higher orders should be possible but may require some extra work. To some extent, this study would be even more rewarding since it could yield non-trivial constraints on the space of theories with classical gravity duals, specifically on the structure of their correlation functions. It would also be interesting to work out a more precise algorithm for metric reconstruction at the full non-linear level following the discussion at the end of section 4.4.2 and, in particular, understand how the Einstein's equations would emerge in this context.
• More general states and entangling regions: in our work we have considered perturbations around the vacuum state and focused on spherical regions in the boundary theory. While it is true that in this setting one could in principle recover Einstein's equations systematically (order by order) and hence the bulk geometry, one could relax these two points, with the latter one being arguably the easiest of the two (at least for regions with local modular Hamiltonians). We note that substantial progress on generalizing the former using the Iyer-Wald formalism appeared in [72]. This paper also argued that doing the linear analysis for arbitrary states and shapes of the entangling region is sufficient to capture the full non-linear Einstein's equations in the bulk. It would be interesting to relax these conditions in our method of bulk reconstruction using bit threads, and try to make these statements a bit more precise in our context. It would also be interesting to explicitly start with a state for which the bulk geometry has entanglement shadows, and understand how our approach encodes information about these regions.

JHEP01(2021)193
• Covariant bulk reconstruction: in this work we have incorporated time-dependence by combining the maximin prescription introduced in [52] and the non-covariant formulation of bit threads [30]. As explained in section 2.3.1, this was possible in virtue of various crucial simplifications of the perturbative setting that we consider. However, it should be possible to pose the same question in the fully covariant formulation of bit threads [53]. We believe that in this case, the full modular flow in the bulk should play a role, and could serve as a guide for constructing canonical bit thread configurations in other special cases. Related to this, it would be interesting to ask the question about bulk reconstruction in time-dependent situations that are not easily accessible to the perturbative setting we consider here, e.g., completely far-of-equilibrium states that could lead to black hole formation in the bulk.
• Higher derivative theories: finally, we believe it would be worthwhile to explicitly generalize our method of bulk reconstruction to the case of higher derivative theories in the bulk. We point out that the program of gravitation from entanglement using extremal surfaces has been worked out in detail for these theories in [28], using the Iyer-Wald formalism both to first and second order in the state deformations around the vacuum. Likewise, the bit thread reformulation of entanglement entropy has already been generalized to the case of higher derivative gravities in [36], incorporating corrections to the local norm bound that depend on the specific theory. It would be interesting to see how our formulas are corrected if we turn on these extra gravity couplings.
We hope to come back to some of these points in the near future.

JHEP01(2021)193
In a general dimension, the metric perturbation can be expanded as H µν = 16πG N d z 2n T (n) µν . However, in d = 2 all n ≥ 1 terms vanish and we have, The stress tensor should be traceless and conserved, therefore, its general form is For general linear perturbations we can perform a Fourier decomposition, and take f (t + x) and g(t + x) to be an appropriate superposition of plane waves. However, for concreteness we will consider an example that is physically motivated: a local quench that arises by the insertion of a local primary operator. The stress tensor of a local quench is given by the sum of two shock waves [73,74], where the (small) parameter λ gives the total energy inserted and α acts as a UV regulator.
As discussed in section 3.2 the unperturbed geodesics provide a family of integral curves that satisfy the criteria given in [43]. In pure AdS, the geodesics are semicircles anchored at the boundary. These semicircles form a two-parameter family of curves, defined by where x s is the center of the circle on the x-axis and R s is its radius. If we denote (x m , z m ) a point on the minimal surface γ A , we showed in section 3.2.1 that where H(t, x) ≡ H xx (t, x). Plugging A.6 and A.7 in A.5 we obtain a family of geodesics orthogonal to γ A , parametrized by the point x m ∈ [−R, R] on the minimal surface. Likewise, we can use the formula (3.43) to obtain the magnitude of the vector field. This calculation can be done following the examples of [43]. The final result for the integral curves and magnitude are plotted in figure 4. We can study the same example using the level set construction discussed in section 3.3. In this approach, given a solution to the max flow problem in the unperturbed geometry v, the solution for the perturbation of δv is  where Ψ is a scalar function that is determined by solving the first order differential equation v · ∇Ψ + ∇ a (δg ab v b ) + 1 2 v · ∇(δg) = 0 , (A.9) with boundary condition In figure 5 we present the results obtained using this method. The third method explored in this paper, discussed in section 4.3, relies on the Iyer-Wald formalism to define a canonical bit thread perturbation. This approach uses the language of differential forms and relates a max flow vector field v, to an optimal closed form w. In a background that is perturbatively close to a given geometry, i.e. g → g + δg, the optimal closed form w → w + δw. Since knowing w + δw determines the max flow v + δv, the problem now amounts to finding δw. The Iyer-Wald formalism provides a form, χ, defined in (4.32), that can be taken as δw, δw = 4G N χ. The form χ is closed when the equation of motions are satisfied. Having δw it is straight forward to determine δv. However, the configuration should satisfy the norm bound (4.22) and that is not guaranteed, in general, in this construction since this condition depends explicitly on the metric. Nevertheless, the more detailed analysis performed in section 4.3, reveals that up to our order of approximation, the norm bound is indeed satisfied. To illustrate this point we plot the norm for the same perturbation studied with the previous approaches. The resulting v and its norm is plotted in figure 6 shows that for a perturbative small λ the norm bound is indeed satisfied.
In figures 4-6 we have presented the perturbed vector field, v λ = v + λ δv, and its magnitude obtained using the three different methods developed in this paper. It is also interesting to look just at the perturbation δv to gain insight into the time-dependence of the local pattern of entanglement induced by the quench. For concreteness we only show results using the level-set method, which are presented in figure 7. In [73] it was conjectured that the quench insertion generates an entangled pair of wave packets that move in opposite directions (see figure 4 of [73] for a pictorial representation). However, it was recently shown in [74] that this intuition is only true if one includes the leading JHEP01(2021)193 1/N corrections coming from the entanglement of bulk fields. More specifically, in this paper it was argued that at the leading order in G N the two wave packets are effectively unentangled, and only the quantum correlations between the degrees of freedom in each individual packet contribute to the total entanglement entropy. Remarkably, we can reach the same conclusion from the plots in figure 7, which exhibit the following features: i) two wave packets moving together with the shocks, i.e., in opposite directions at the speed of light and ii) threads around each wave packet connecting degrees of freedom in their fronts with those in their tails. The fact that we do not see threads connecting the two wave packets implies that they are effectively unentangled at the leading order in G N , in agreement with the result of [74]. Moreover, the precise pattern of the threads explains why S A peaks at t = R (see e.g. figure 7 of [74]): at this time most of the threads connect the degrees of freedom of A with those in its complement (recall that threads connecting points within A do not contribute to the entanglement entropy of the region). It will be very interesting to repeat this analysis for the case of a global quench, and understand how the local pattern of entanglement evolves in time for cases that admit an entanglement tsunami interpretation [75,76] (large regions) and cases that do not [77,78] (small regions).