Building Bulk Geometry from the Tensor Radon Transform

Using the tensor Radon transform and related numerical methods, we study how bulk geometries can be explicitly reconstructed from boundary entanglement entropies in the specific case of $\mathrm{AdS}_3/\mathrm{CFT}_2$. We find that, given the boundary entanglement entropies of a $2$d CFT, this framework provides a quantitative measure that detects whether the bulk dual is geometric in the perturbative (near AdS) limit. In the case where a well-defined bulk geometry exists, we explicitly reconstruct the unique bulk metric tensor once a gauge choice is made. We then examine the emergent bulk geometries for static and dynamical scenarios in holography and in many-body systems. Apart from the physics results, our work demonstrates that numerical methods are feasible and effective in the study of bulk reconstruction in AdS/CFT.


Introduction
Recent progress [1][2][3][4][5][6][7] in quantum gravity has shown that spacetime geometry can emerge from quantum entanglement. This emergence provides appealing explanations for many intuitive properties of the physical world, including the existence of gravity [8][9][10][11], conditions on the allowed distribution of energy and matter [12,13], and the unitarity of black hole dynamics [14,15]. Most of these developments have taken place in the context of the Anti-de Sitter/Conformal Field Theory (AdS/CFT) correspondence [16,17], although some proposals relating geometry and entanglement also apply to flat or de Sitter geometries [5,7,[18][19][20]. As an example of the holographic principle [21,22], AdS/CFT describes a duality between a (d + 1)-dimensional bulk theory with the presence of gravity in asymptotic AdS spacetime, and a d-dimensional conformal field theory (CFT) without gravity on the boundary. For the purposes of studying quantum gravity, the duality is especially powerful because the CFT is an object whose basic rules we understand. However, much remains to be understood about the relationships between the two sides of the duality.
One such challenge in the AdS/CFT duality is to understand how the boundary degrees of freedom without gravity can reorganize themselves into a higher dimensional bulk configuration with gravity. This is called the problem of bulk reconstruction, and this paper reports two results on this topic. First, we describe a perturbative procedure to reconstruct the bulk geometry given an appropriate set of boundary entanglement data. Second, we show that this reconstruction procedure can detect whether the putative bulk dual is semi-classical in the sense of having only weak fluctuations about an average value.
Our first result builds on a number of works that study bulk metric reconstruction using geodesic lengths [23,24] or entanglement entropy [25]. 1 Early efforts in this area often utilized bulk symmetries to simplify the problem of recovering the bulk metric from minimal geodesic data [33][34][35]. Later work on differential entropy and hole-ography [25,36] furthered our understanding using information theoretic quantities and suggested a method to recover the bulk, although no explicit reconstruction formula was given for generic cases. More recently, it was shown that in certain cases, static geometries [37], or even the full dynamical metric [38], can be fixed non-perturbatively by boundary entanglement data without any prior knowledge of the bulk symmetries. However, barring a few well-known examples with symmetries, there does not exist a reconstruction procedure that directly and explicitly converts entropy data into bulk metrics. Our work addresses this missing element.
Our second result arises from the basic issue that we have a far from complete understanding of what kind of boundary states correspond to semi-classical bulk geometries. Some necessary conditions are expected from holographic entropy inequalities [39] and from consistency relations for any putative metric reconstruction procedure 2 . In general, we do not expect all quantum states from the boundary CFT to correspond to welldefined semi-classical geometries in the bulk. On the contrary, we expect an abundance of non-geometrical states obtained, for instance, by superposing states with macroscopically distinct dual geometries such as the AdS vacuum and an AdS black-hole geometry. Therefore, to better understand holographic duality, we must also address the necessary and sufficient conditions under which a bulk geometry can emerge from the boundary state. Our work also addresses this open question.
In this paper, we study the above issues in the context of AdS 3 /CFT 2 using an approach based on the tensor Radon transform. The method is for metric solutions that are close to AdS 3 , hence it is restricted at present to perturbative problems. To linear order in the perturbation, each quantum state on the boundary corresponds to a constant time slice in the bulk, so we provide a numerical reconstruction algorithm that takes the entanglement entropies of intervals in the boundary state as input, and outputs the best-fit bulk metric tensor of the spatial slice in the linearized regime. This solution is unique up to gauge transformations. The algorithm also provides a quantitative indicator of whether the boundary data admits a bulk geometric description near AdS 3 . This is measured by the quality of the fit, which intuitively quantifies how far the boundary entanglement data of a given state is from being geometric. A poor quality of the fit indicates the boundary data lack consistency with a semi-classical geometry.
As a proof of principle, we explore several reconstructions numerically in holography and with a 1d free fermion CFT. We find that free fermion ground states in the presence of disorder and a mass deformation do not correspond to a well-defined bulk geometries. Likewise, mixtures of states where each is dual to a distinct classical geometry can also fail to have a well-defined bulk geometry. In addition to these static examples, the method is applied to several dynamical scenarios including global and local quenches in the free fermion model and entanglement dynamics in a toy model scrambling system. In the case where the dynamics is scrambling, we find that the bulk description is qualitatively consistent with an in-falling spherical shell of bulk matter experiencing gravitational attraction. These results further demonstrate that it is both feasible and interesting to study AdS/CFT using tensor Radon transform techniques coupled with modest (laptop-scale) computational resources.
The remainder of the paper is organized as follows. In section 2, we briefly review the basic assumptions, especially how entanglement entropy can be tied to metric tensors via the tensor Radon transform. We give a general review on the tensor Radon transform in Appendix A. Its adaptation to a hyperbolic geometry and the gauge fixing procedure are found in Appendix B and Appendix C, respectively. In section 3, we introduce the numerical reconstruction procedure, which we elaborate in detail in Appendix D. In section 4, we apply the reconstruction algorithm to static and dynamical boundary entanglement data. In section 5, we discuss these reconstructions and how geometrical and non-geometrical can be distinguished using the relative reconstruction error. Finally we conclude with some remarks and directions for future work in section 6.

Boundary Rigidity and Bulk Metric Reconstruction
We begin with the vacuum state |0 CFT of a hologaphic CFT. Because this state has conformal symmetry at all scales, it must be dual to empty AdS [16]. When the bulk theory is Einstein gravity coupled to matter, the bulk geometry controls the leading entanglement structure of the CFT state via the Ryu-Takayanagi (RT) formula [2,3]. Given a boundary region A, the RT formula computes the von Neumann entropy S(A) in terms of a minimal area surface, where Area[γ] is the area of γ A , G is Newton's constant, and the minimum runs over bulk surfaces γ A that are homologous to A. In the time symmetric case where RT applies, all of these surfaces can be taken to lie in a time-symmetric spacelike surface Σ. In other words, the Ryu-Takayanagi formula says that the von Neumann entropy of a state on a boundary subregion is given by area of the minimal area bulk surface that subtends the region.
If we have access to a boundary state, in the sense of knowing its von Neumann entropies on all (connected) subregions, then the RT formula translates these entropic quantities into a set of boundary anchored minimal surface areas. Because these minimal surfaces all lie on the spatial slice Σ, recovering the bulk geometry from entanglement reduces to a pure geometry problem where we try to find the interior metric g ij of a Riemannian manifold M while knowing only the areas of minimal surfaces that are anchored to its boundary ∂M. This is precisely the statement of the boundary rigidity problem [40], which is well-studied in the field of integral geometry [41].
For this work, we focus exclusively on the case of AdS 3 /CFT 2 , where minimal surfaces are simply geodesics and spatial slices are 2d Riemannian manifolds. For a class of 2d Riemanian manifolds (called simple manifolds), it is known that the lengths of all boundary-anchored geodesics indeed fixes the bulk metric uniquely up to gauge equivalence [42]. Because the single interval von Neumann entropies in the ground state of a 2d CFT are universally determined by the central charge [43], the RT formula combined with boundary rigidity completely fixes the bulk geometry to be that of hyperbolic space [44,45]. Of course, this is precisely the induced geometry on a time-symmetric slice of AdS, as it had to be based on the grounds of symmetry 3 .
It is natural to ask whether we can exploit the power of the RT formula and the results from boundary rigidity theory to reconstruct dual geometries from boundary states other than the vacuum. For instance, given a generic non-vacuum state |ψ CFT , can we apply the same principles to reconstruct the metric tensor for the bulk geometry from the set of boundary-anchored geodesic lengths? There are several obstacles that prevent us from recovering the bulk metric tensor exactly using the above methods, even if |ψ CFT has a well-defined dual geometry. First, since it may not be possible for minimal surfaces (or extremal surfaces in the dynamical case [46]) to foliate entire space(time) manifold, thus we can only reconstruct regions where there is at least a local foliation with minimal or extremal surfaces. Second, even for the regions whose geometries are fixed by the entanglement data [38], there is no explicit reconstruction formula for the general boundary rigidity problem.
Although such problems are difficult to solve in general, it is typically easier to reconstruct the difference between the dual geometry and a known reference, or background, geometry. This is known as the linearized boundary rigidity problem [40,41]. In this work, instead of a direct reconstruction of the of the dual geometry for |ψ CFT , we reconstruct the differences in the entanglement patterns as linearized metric perturbations around the AdS background.
To do so, we first fix the background geometry to be vacuum AdS 3 . Working with a given constant-time slice, let us suppose that |ψ CFT has a slightly different entanglement structure compared to |0 CFT , and is dual to a bulk geometry with a metric that is close, but not equal, to that of pure hyperbolic space on our time-slice. Then by the RT formula, the boundary-anchored geodesic lengths now differ slightly from those of pure hyperbolic space. For a given boundary subregion A, the change in the geodesic length anchored at A is related to the vacuum subtracted entropy of the state by where L ψ (A) and L 0 (A) denotes the lengths of geodesics anchored at the end points of A for states |ψ CFT and |0 CFT , respectively.
This change in geodesic length corresponds to a change in the bulk metric where g (0) ij is the pure hyperbolic metric and h ij is the perturbation. For the linearized problem, the goal is to find h ij to leading order, with h ij viewed as a rank two symmetric tensor field on the hyperbolic background. Note that geodesics of the background metric remain geodesics of the perturbed metric to first order in h since geodesics satisfy an extremality condition. Changes in entanglement due to variations in the minimal surface itself are of order h 2 [8]. The leading order change in geodesic length can then be written as where γ A is the geodesic of the hyperbolic background anchored at the end points of A andγ i A denotes unit tangent vectors along γ A . The tangent vectors are normalized such that For simplicity, we will use L(γ A ) and L(A) interchangeably, often dropping the explicit dependence on A and simply writing L(γ) when there is no confusion, with the understanding that γ is a boundary-anchored geodesic.
To make concrete progress in this work, we simplify the full problem by taking the linearized bulk result for changes in geodesic length to be equal to the full change in boundary entanglement entropy, This approximation enables calculations; going beyond it may be technically non-trivial, but it is likely not a fundamental obstacle.
With this simplification, the length perturbation becomes the integrated longitudinal projection of the metric perturbation along a geodesic γ. This is precisely the tensor Radon transform R 2 [h] of the metric perturbation h ij [47]. Given a symmetric 2-tensor field h ij , the tensor Radon transform R 2 [h] defines a map from the space of boundary geodesics to the complex numbers given by where γ A is the background geodesic anchored at the boundary of A.
As an aside, we note that there exist several related notions of Radon transform. A standard Radon transform on a Riemannian manifold is defined by integrating some quantity on a minimal co-dimension one surface, whereas an X-ray Radon transform is defined similarly, but for a dimension one surface, i.e., a geodesic. For two spatial dimensions, as is the case we consider here, the two definitions coincide. For more details on the Radon transform, see Appendix A.
Formally, the the bulk metric deformation can be recovered by inverting the tensor Radon transform. Schematically, we can write where R −1 2 denotes some (not yet properly defined) inverse Radon transform, and where ∆L denotes the collection of boundary anchored geodesic length deviations.
Throughout, we work exclusively in the perturbative regime, to leading order in h, which allows us to relate geodesic length deformations to the Radon transform through (4). It also ensures that the resulting geometric solution, when it exists, is uniquely determined by the boundary entropy data [42,45]. 4 We wish to comment here that despite restrictions to the perturbative regime, h ij can still capture highly non-trivial physics. Indeed, standard calculations of gravitational waves and the dynamics of typical stars, planets, and galaxies are all done in the weak-field regime.
Before we can proceed with the inversion process, we must give meaning to the inverse Radon transform. For this purpose, it is important to note that the Radon transform has a non-trivial kernel: given any vector field ξ on M such that ξ| ∂M = 0, we necessarily have Physically, any Radon transform of a pure gauge deformation that reduces to the identity at the boundary is zero. The presence of this kernel is natural because the transform relates geodesic lengths (which are gauge invariant) to metric tensors (which are not), so the Radon transform can only be injective up to gauge.
In the presence of such a kernel, we must fix a gauge prescription in order to recover a metric tensor uniquely. We will use a prescription which we call the holomorphic gauge [48]. In a crude sense, the holomorphic gauge preferentially reconstructs the trace part of the metric at the cost of diminishing non-zero contributions to the off-diagonal. 5 This provides two independent gauge constraints in two spatial dimensions, which would allow us to proceed with the reconstruction. For a more detailed description of the gauge constraints and Radon transform on a hyperbolic background, see Appendix B.
As a final comment, in the above discussion we have explicitly assumed that the boundary state corresponded to a well-defined bulk geometry. Below we will construct algorithmic machinery which can actually carry out the reconstruction in this case. However, we will also see that the algorithm can be applied to a more general class of states, with interesting results.

Numerical Methods for Reconstruction
The inversion formula for the flat-space scalar Radon transform is a well-known classical result in integral geometry [41]. Explicit reconstruction formulas for scalar and vector Radon transforms are also available for constant negative curvature backgrounds [42,49]. However, there are currently no explicit reconstruction formulas available for higher rank tensors on curved backgrounds, although several results in the literature come close to a solution in various regimes [42,[48][49][50][50][51][52]. In the absence of an exact analytic reconstruction formula, we instead draw inspiration from the general principles employed in seismology to study the Earth's interior [53] in developing our numerical method.
In this section, we give a brief overview of the method. The full details of the discretization, gauge fixing, and solutions for the constrained least square problem can be found in Appendix D.

Discretization and Optimization Procedures
The basic idea behind our numerical reconstruction is straightforward. We first discretize the bulk and boundary regions into a finite number of tiles. To each tile T in the bulk, we associate a tensor h ij (T ), and for each interval A on the boundary, we associate a geodesic (of the background metric) γ A anchored on the endpoints of the interval. In two spatial dimensions, a rank-2 symmetric tensor has 3 independent degrees of freedom. Each geodesic anchored at the end points of an interval generates a linear equation via the discretized version of the Radon transform (7), defined by where the tensor W ij contains information about the direction of the tangent vectorsγ A , as well as the arc length ∆s(γ A , T ) of the geodesic segment that passes through each tile T . In equation (10), we sum over all bulk tiles T and over repeated indices. Naturally, one has W ij (T , γ A ) = 0 if the geodesic does not pass through a tile T .
We can abbreviate equation (10) in matrix form as where W and h are vectorized representations of W ij and h ij , and where b denotes the corresponding geodesic length deformations. As a result, given a specific discretization, the discretized forward Radon transform can be written as a linear map W : V B → V γ from the space of tile-wise constant bulk tensor valued functions V B to the space of boundary anchored geodesic lengths V γ . Both spaces are finite dimensional due to the discretization.
Since the forward Radon transform has a non-trivial kernel in the continuum limit, we must impose a gauge fixing condition to recover a unique solution. We give the full detail of the gauge fixing conditions and the accompanying partial differential equations in Appendix C. To ensure the problem is well-posed, we set the holomorphic gauge constraints as discretized partial differential equations, which we formally write as Here, C denotes the constraint matrix representing the partial differential operator associated with the gauge constraint. Reconstruction of the metric perturbation then corresponds to finding the solution to the linear equations (10) above, subject to the linear constraints (12).
In practice, there does not always exist exact solutions h ij which satisfies the constrained system. This can be due to a variety of reasons, such as the presence of discretization errors, or if the boundary data is simply inconsistent with a geometric bulk, i.e., if the boundary entropy function fails to lie within the range of the forward Radon transform [41,42,48]. Instead of trying to look for an exact solution, it is more natural to look for the best-fit solution h * which solves the constrained minimization problem The objective function is linear and we are guaranteed a unique global minimum. Thus we will say that h * is the optimal geometric solution corresponding to boundary data b. We will also write h * (b) when we need to denote the dependence of h * on the initial boundary data.
Even with the existence of an optimal reconstruction h * , we do not expect generic boundary data to correspond to a geometric dual in general. A useful quantity is the relative error of reconstruction, which measures the tension between the best-fit solution and the actual data. We can consider various relative errors. If we know the exact bulk solution, say h 0 , then we can denote the bulk relative error as where h * is the bulk metric tensor reconstructed from the forward transform of h.
More commonly however, we do not have access to an a priori geometric state. Instead, we have a CFT state |ψ CFT from which we can extract discretized boundary data b ψ . In this case, we can likewise consider the boundary relative error, defined by The boundary relative error is simply the normalized distance from b ψ to the subspace of boundary data vectors with geometric duals, which is the same quantity minimized by the constrained least squares problem (13). The boundary relative error therefore serves to quantify the degree to which a state is geometric or non-geometric. We will discuss in greater detail in section 5 how the relative errors can be used to distinguish geometric data from non-geometric data on the boundary.
To ensure the reliability of the reconstruction algorithm for the inverse tensor Radon transform, we also perform the numerical inversion of boundary data b whose bulk tensor field is known. We produce such boundary data by preparing various known bulk tensor valued functions h 0 in the holomorphic gauge then generating its corresponding geodesic data b 0 through a forward tensor Radon transform. Subsequently, we apply the numerical reconstruction to the geodesic data and compare the reconstructed h * (b 0 ) to the original test function h 0 . We find remarkable agreement in our reconstructions. Absent rigorous analytic convergence guarantees, the successful benchmarking of the algorithm on known cases serve to provide confidence in the fidelity of the reconstruction. Details of this benchmarking process are given in Appendix E.

Reconstructed geometries
The forward tensor Radon transform is generally neither injective nor surjective. Therefore, not all boundary data can be interpreted as the Radon transform of a bulk tensor field. However, it may be possible to derive necessary and sufficient characterization of what boundary data corresponds to a bulk tensor field. Such a criterion is known as a range characterization of the tensor Radon transform. Ranges characterizations of various transforms have been discussed extensively for scalar and vector cases on curved backgrounds and tensor cases on flat background [41,42,48]. However, it remains an active topic of research for transforms on curved background for higher rank tensor fields.
Although we lack a rigorous analytic characterization for the tensor Radon transform, our numerical methods can still effectively capture the parts of the entanglement data that do not lie within the range of the tensor Radon transform. Since the Radon transform is linear, any contribution that cannot be fitted to a bulk tensor field in the global best-fit reconstruction effectively captures the non-geometric contribution for the discrete reconstruction. More specifically, the relative boundary error (15) serves as a indicator for the fraction of the boundary entropy data which does not lie within the range, i.e., the fraction which can be considered non-geometric.
In the upcoming sections, we reconstruct geometries from boundary entanglement entropies generated by holographic systems as well as those generated numerically from a 1d free fermion system. We then discuss how their relative boundary errors E bdy can be used as a standard to distinguish states that have a welldefined classical dual geometry from the ones that do not. This provides a direct quantitative condition for whether a state is geometric.
For entanglement entropies generated by holographic systems, we expect a weakly coupled gravity dual, and therefore spacetime geometries that are generically classical at low energies. The same can not be expected for a generic many-body quantum system at criticality [54], although they may still capture interesting features using a dual spacetime prescription. For instance, in a free fermion system, the conformal field theory has neither strong coupling nor a large central charge. Although it is unclear what the dual bulk description would be, it is generally expected that any geometric description must be one where gravity is strongly coupled, the leading order RT formula does not apply, and the spacetime is dominated by quantum effects. Therefore, the reconstruction is likely poor in the absence of a large number of symmetries. We will find that numerical evidence support these basic intuitions.
It is difficult to present certain reconstructions from dynamical systems in a paper. A list of the animated reconstructions for various dynamical processes are linked here. 6

Holographic Reconstructions
As a benchmark for holographic geometry reconstructions, we first reconstruct the metric for the thermal AdS geometry in the linearized regime. The entropy data we use are generated by the minimal geodesic lengths in a BTZ geometry [55] using the Ryu-Takayanagi formula. Since the data corresponds to a bulk geometry by construction, the reconstruction should show good agreement with said geometry at linear order. This includes the correct qualitative behaviour in the bulk, and a positive metric perturbation present deeper into the bulk. It must also afford relatively small reconstruction errors, which can be attributed to factors such as discretization errors, approximations made in linearization, and working on a fixed hyperbolic background despite the changes in bulk geometry.
Indeed, as shown in Figure 1, we find good agreement with our expectations in the reconstruction. By probing deeper into the bulk, the geodesics for a thermal geometry travel through larger distances, resulting in a net positive change in the metric perturbation. Note that at linear order, we can not detect a change in bulk topology, which the entanglement data dual to a BTZ geometry should predict. This is because of the fixed background metric. 7 Motivated by entanglement dynamics in holography, we can also construct heuristics for the entanglement growth of a boundary theory. For instance, under a global quench, we expect qualitative growth of entanglement to be captured by (a) Component h11.
(b) Linearized scalar curvature perturbation. Figure 2: Metric perturbation from volume law entropy growth (16). Plots are ordered from left to right as time increases. We only give h 11 for the sake of clarity, because h 12 ≈ 0 and h 22 ≈ h 11 . The boundary relative error is E bdy ≈ 0.03.
for t ≥ 0, where entanglement for any region will grow linearly in time after the initial stage, until the entropy satisfies a volume law [56,57]. Here s denotes the entropy density, v the speed of entropy growth, and t the time that has elapsed since the quench. The size of the boundary interval is denoted by |A|. The corresponding metric and curvature perturbations are shown in Figure 2.
As the system thermalizes, larger subregions have a volume law entropy. The wavefront of the entanglement spread is reflected in the bulk as a spherically symmetric perturbation moving from the boundary to the center. Assuming Einstein gravity, which holds for holographic CFTs, the curvature perturbation δR also reflects the bulk matter distribution through the linearized Hamiltonian constraint for each instance of time [47]. Therefore, this thermalization process is consistent with the collapse of a spherical shell of matter [56].

Mixture of Thermal States
There are also instances where we do not expect a well-defined geometry to exist. For example, when the state is taken from a theory where the bulk is strongly coupled and/or the state is a macroscopic superposition of certain classical geometries, quantum gravitational effects can dominate, leading to a breakdown of the classical geometric description.
In this section, as an example of a potentially non-geometric holographic state, we look at mixtures of thermal AdS geometries at various temperatures. We will consider states of the form where ρ(T i ) are thermal states of a CFT with distinct temperatures T 1 and T 2 . From [58], the von Neumann entropy of the mixture is estimated as where we write ρ i,A to denote the reduced state of ρ(T i ) on a boundary region A, and where is a Shannon-like term that corresponds to the entropy of mixing [59]. Physically, such a state can be created by superposing two thermofield double states at different temperatures, and then tracing out one side of the wormhole. For CFTs with strong bulk gravity where G N ≈ 1, we find the superposition incurs a large error, indicating non-geometric configurations in the bulk. However, for a CFT with a weakly coupled dual where G N 1, the geometry smoothly interpolates between the two temperatures at the linearized level, consistent with our expectation that the entropy operator is proportional to the area operator at leading order in N . In the language of Radon transform, this is caused by the entropy of mixing H(p) being a non-geometrical contribution. From left to right, we increase the mixing ratio p and the geometry transitions smoothly from one at lower temperature to the other at higher temperature when the gravitational coupling is weak. The reconstruction is dominated by artifacts and have large error if the coupling is strong, leading to a contribution G N H(p) in the leading order.

1D Free Fermion
Given that the reconstruction procedure relies only on the von Neumann entropy of a state, we may also apply it to various quantum many-body systems where we do not necessarily expect the conformal field theory to be dual to a semi-classical bulk theory with weakly coupled gravity. The lack of a semi-classical geometric dual will generally be reflected in the presence of a large boundary reconstruction error. One example where the separation between geometric and non-geometric states is particularly manifest is the contrast between the reconstruction of holographically motivated data and that of a 1d free fermion.
Generic low-energy states for free fermion systems are generally believed to have a highly quantum bulk because the conformal field theory is non-interacting, and has a small central charge. Heuristically, the strong-weak nature of the holographic duality suggests that such systems should result in a strongly coupled bulk geometry, if such geometries even exist. In this section, we give a few examples indicating that generic excited states in the free fermion system are indeed non-geometric. However, we shall also see that certain large-scale geometric properties may nevertheless be present, even if the overall state is ostensibly nongeometric.
For our reconstruction, we consider a 100-site massless free fermion system with periodic boundary conditions, and with HamiltonianĤ where, as usual, the creation and annihilation operators satisfy the canonical anti-commutation relations We shall also consider various quench dynamics and deformations of the free fermion system. The free fermion system is described by a (1 + 1)-d conformal field theory in the thermodynamic limit, with central charge c = 1/2 [60]. Here we study directly the regulated lattice model. Any such model of noninteracting fermionic modes can be studied efficiently numerically with a cost scaling polynomially with the number of fermion modes N . In particular, we can compute all of the single interval entanglement entropies via [61]. The key point is that the state of the system is Gaussian. In particular, given a subset A of the fermions, the reduced density matrix is ρ A ∝ e −K A with K A quadratic inâ andâ † , Here k A is an |A| × |A| Hermitian matrix that determines all correlation functions of fermions in A. It is related to the 2-point function via The entropy of A is then determined by the eigenvalues of k A , but there is also a direct formula in G A : To fix our background geometry, we numerically compute the ground state entanglement of the critical HamiltonianĤ 0 . To consider non-vacuum emergent geometries, we consider states |ψ = |0 , which are not necessarily energy eigenstates ofĤ 0 . For instance, these excited states can be generated by first deforming the Hamiltonian away from criticality through a mass deformation and then finding the ground state |ψ m of the perturbed Hamiltonian. For small m, the new ground state |ψ m is generically an excited state with relatively low energy in the unperturbed system.
We will study reconstructions in both the static and dynamical cases. For the static case, we reconstruct the emergent geometry of the new ground states |ψ m by finding the best-fit geometry using the discrete Radon transform. We do this for a number of distinct values of m. For the dynamical case, we consider the quench dynamics corresponding to the deformationĤ 1 . We start with some fixed deformed ground state |ψ m and then time evolve the state with the free HamiltonianĤ 0 and reconstruct the geometry that corresponds to each time step: |ψ m (t) = exp(−itĤ 0 )|ψ m .

Local Deformations
In this section and the next, we will consider deformations of the form where m is a small positive parameter, and where S is a set of sites where we introduce such deformations. The perturbation serves to deform the ground state wavefunction around the sites near S. We label the sites counter-clockwise from 0 to 99, starting with site 0 aligned with the positive x-axis (i.e., site i sits on the unit circle at angle θ i = 2πi/100). We will consider local perturbations where S consists of one or two sites in this section, and more global perturbations in the next section. Figure 4 shows the best-fit geometry arising from the ground state |ψ m corresponding to a local mass deformation at a single site S = {0}, located along the positive x-axis. We can see that, in contrast to the holographic reconstructions in section 4.1, the geometries shown in Figure 4 are heavily dominated by noise and local artifacts. The corresponding relative error of reconstruction is also larger by more than an order of magnitude as compared to the holographic cases (see Figure 13a in section 5 for a more comprehensive comparison of relative errors). The large error of reconstruction can be seen as a hallmark of the fact that the underlying state is non-geometric. Figure 4: Components of the reconstructed metric tensor perturbation corresponding to the ground state of a mass deformed Hamiltonian, with the deformation located at a single site located along the positive x-axis. The boundary relative error is E bdy ≈ 0.62.
Note that the lack of geometric features in the reconstruction does not necessarily indicate that the reconstruction procedure is flawed. In fact, non-geometric features are expected because not all boundary data should produce geometric reconstructions. Our reconstruction procedure has equipped us with information of telling apart when that will happen. We will discuss this further in section 5.
Nevertheless, there are cases where some large-scale geometric features can be extracted from the plot. This is especially clear in the dynamical case, where we consider the quench dynamics obtained by evolving the deformed ground state |ψ using the free fermion HamiltonianĤ 0 . The reconstruction, performed timeslice by timeslice, is shown in Figure 5 (only the dominant component h 11 is shown). While it may not be entirely clear, Figure 5 reveals a large scale shockwave that originates from the deformation site and then travels across the bulk before being reflected at the left boundary. The large-scale geometry is heavily obscured by non-geometric noise in Figure 5. To better extract the large-scale features that we consider to be relevant, we must filter out the small-scale artifacts what we consider to be "noise".
To more clearly extract the underlying large-scale features of what appears to be waves generated by the perturbations, we perform some pre-processing of the boundary data to smooth out the small-scale noise. This is done by averaging the entanglement entropies of two neighbouring intervals of the same size. More precisely, suppose S(i, j) is the von Neumann entropy on the interval from site i to site j. Then the smoothing procedure we applying is similar to a filter by convolving with a 2-site window function such that This produces a less noisy reconstruction and significantly reduces the amount of non-geometric contributions to boundary relative errors. Intuitively, the smoothing procedure acts as a low-pass filter in the space of boundary intervals. This serves to remove the short distance artifacts normally associated with nongeometricality. While we used a specific filter in this example, one may apply more general filter constructions for other purposes. 8 We wish to emphasize that this filtering procedure is not a part of our reconstruction process. It is merely here to assist us in gaining some intuition regarding the qualitative physical behaviour of this type of dynamical process. One could imagine that the filtering reveals what the bulk physical process might have looked like, had the state actually been geometric.
The corresponding smoothed geometries are shown in Figure 6. It can be seen that smoothing significantly reduces local noise, and clearly reveals large-scale time dynamics associated with the quench that looks like a (shock)wave. 8 It is reasonable to suspect that the two site averaging of entanglement entropy works well as a filter here because UV details in the free fermion Hamiltonian (20) with Fermi momentum ±π/2 contributed to non-geometric noise. We might therefore expect that we can remove such non-geometric contributions through coarse-graining by grouping together adjacent sites in pairs. While this grouping does indeed remove some non-geometric noise, it does not remove it entirely: the reconstruction error after grouping is 0.1 E 0.3, which is in between what is shown in Figure 5 and Figure 6. Given that the smoothed data reconstruction error is still significantly larger than that of completely geometric data, it is not clear if the non-geometric contribution completely reside in the UV. Figure 6: Quenched time dynamics of the reconstructed geometry corresponding to a single mass deformation along the x-axis. The boundary data is pre-processed by smoothing to remove small-scale details (see Figure 5 for the corresponding unsmoothed reconstruction). Only the dominant component h 11 is shown here. The relative error is 0.03 E bdy 0.1.
As another example, the quench dynamics of a state with two distinct deformations at sites S = {0, 30} is shown in Figure 7    In all cases, we see that the free fermion Hamiltonian will generate seemingly non-interacting waves that traverse through the hypothetical bulk spacetime. Since the underlying time dynamics is integrable, the same entanglement feature recurs after the waves traverse the entire system; we show one such iteration in our figures.

Global Deformations
Similar to the previous analysis for local deformations, we can also consider global deformations in the same vein. In this section, we will consider Hamiltonian perturbations of the same form as (26), but now with deformations located at every other site, i.e., S = {0, 2, 4, · · · , 98}.
In Figure 9, we show the geometries reconstructed from the ground state |ψ m of a Hamiltonian with the aforementioned global deformations, plotted across a range of m values. Again, we find the overall geometry to be highly dominated by noise, with a large relative error of reconstruction indicating that the underlying state is non-geometric. The relative error of reconstruction generally becomes worse with increasing values of m (see Figure 13b). Qualitatively, the reconstructions for the different values of m appear similar, the main difference being the overall magnitude of the metric perturbation.  Figure 13b).
Looking at the corresponding quench dynamics (see Figure 10 for the unsmoothed reconstruction and Figure 11 for the smoothed version), we see that the large-scale geometry involves a spread of entanglement that is qualitatively similar to the configuration obtained in the thermalization scenario considered in section 4.1. However, in this case the bulk "matter" is non-interacting. The integrability of free fermion Hamiltonian ensures that the falling shockwave returns to the boundary after some finite time instead of collapsing into a steady configuration.

Random Disorder
Finally, we can take a look at the geometries arising from a random perturbation of the free fermion system. Disorder is introduced by adding a random external field at each site of the spin chain where each parameter w i is a random parameter chosen i.i.d. from a uniform random distribution over the interval [−0.1, 0.1].
In Figure 12, we show a generic sampling of ground states obtained from such random disorder. As would be naively expected, the resulting ground states have generically large relative errors, with no discernible large scale features (quench dynamics also reveal no discernible large-scale patterns) or symmetries. Figure 12: A bulk best-fit reconstruction of a state generated with random disorder. Each row is a reconstruction of a particular instance with random disorder. The boundary relative errors vary across a wide range of values, but typically 0.05 E bdy 1 (see Figure 13a).

Geometry Detection
In this section, we summarize the results of the previous reconstructions and comment on their similarities and differences. We also emphasize that a small relative error of reconstruction is indicative of a classical bulk geometry in which boundary entropies are computed via the RT formula, whereas a large error may indicate some combination of (1) no classical geometry, (2) no classical geometry near the AdS background, or (3) a classical geometry but where entropies have non-trivial contributions from sources other than the area of RT surfaces, e.g. higher order corrections. For the purpose of this discussion, we refer to all three of these negative cases as non-geometric, but it would clearly be desirable to distinguish them further in future work.
From the results of the mass deformation, we confirm that generic low-energy states of a free fermion system do not appear to have a good geometric reconstruction. Comparing the ground states of the mass deformed 1d free fermion with that of the thermal AdS state, we see that the relative errors for the 1d free fermion are significantly larger than those of the holographic data (see Figures 13a,13b,14). Smoothing of the boundary data reduces the relative error of reconstruction by a significant amount (as would be expected due to the reduction of small-scale artifacts). However, even with smoothing, the level of error is clearly distinguishable from the reconstructions of known geometric states such as thermal AdS.
In particular, ground states of the locally deformed 1d free fermion Hamiltonian have a relative error that is of order ∼ 1. Smoothing of the boundary data brings this down to ∼ 10 −1 , which is still an order of magnitude larger in comparison to the thermal AdS reconstructions, which have a relative error on the order of ∼ 10 −2 . The unsmoothed relative error for the global deformation of the 1d free fermion hovers around ∼ 10 −1 . Smoothing brings this down significantly to ∼ 10 −2 . The most likely explanation for the pronounced effect of smoothing here is due to the global symmetry present in the globally deformed state. We also see that the relative errors for reconstructions corresponding to random disorder in the 1d free fermion tends to interpolate between the results for local and global deformations across different random instances, with the best case relative errors being 10 −1 and the worst case being ∼ 1. A summary of these results is plotted in Figure 13a, which shows a histogram of the number of instances for each type of reconstruction, as we vary some of the relevant parameters (i.e., temperature, mass, etc). In Figure 13b, we also plot the relative error of reconstruction as a function of the mass deformation parameter m. As would be expected, we see that a larger mass deformation, and hence a larger deviation from criticality, contributes positively to the relative error of reconstruction regardless of smoothing.
(a) A histogram of relative errors from different boundary data: Thermal AdS (TAdS), global mass deformation (MD), local quench without smoothing (LQ), and system with random disorder (Disorder). Boundary relative error E bdy plotted on the x-axis.
(b) Global mass deformation errors, as a function of the deformation parameter m. Red dots denote the boundary relative error of the original fermion data. Blue squares are relative errors after smoothing.

Figure 13
A summary of the time dynamics of the relative errors is shown in Figure 14. For the dynamical reconstructions, we see that the relative error generally varies with time. This is most noticeable with the global quench dynamics, where the relative error tends to decrease as entanglement spreads. Although one starts with a mass deformed state with large relative error at t = 0, subsequent evolution can effectively wash out non-geometric features the system begins to thermalize. The state at t ≈ 20 captures basic entanglement properties of a thermal state, thus the reconstruction is thermal AdS-like with small relative errors. Nevertheless, the state cannot actually thermalize due to integrability; the system recurs at later times and the relative error rebounds.
At this point, it is not completely clear why the global quench states have such small relative error at intermediate times. A likely guess that the enhanced symmetries for the globally deformed states contribute to a more geometrical reconstruction, with the effect being especially pronounced at intermediate times when the system is maximally thermal, before it is kicked back by integrability. Such behaviour is not seen in the local quench disorder, where the magnitude of the relative error is larger, and remains stable across the time evolution. However, symmetry is likely not the full answer either, since otherwise we would expect the disordered models to generically have larger relative error compared to the global mass deformed state, whereas we observer the disordered model to be peaked at smaller errors. While the full characterization for which states appear more geometric than others is currently unknown, it seems likely that there are many contributing factors to a successful reconstruction, among them properties like symmetry and thermality. Finally, we show the relative errors corresponding to the superposition of thermal states in Figure 15. We find two regimes for the geometricality of the superposition, depending on the magnitude of the mixing term G N H(p). We find the superposition to be non-geometric when the entropy of mixing provides a significant correction to the geodesic lengths. This implies a generally non-geometrical construction when G N H(p) ≈ 1, as would be the case when the gravitational coupling G N is strong, or when the entropy of mixing is made large by superposing a large number of distinct geometries.
Conversely, in the weak coupling limit with few distinct superpositions, the contribution of the mixing term G N H(p) to the geodesic lengths becomes negligible, causing the entropy operator to be approximately linear in this regime. This leads to a smooth, geometric interpolation between the two distinct geometries as we tune the probability of mixing. The clear separation between the two cases is clearly illustrated in Figure 15, where we plot the relative errors for a state constructed as a mixture of two distinct thermal states. In summary, we find that the boundary relative error provides a useful measure for the extent to which a state is geometric or non-geometric. Through the discretized Radon transform, we confirm some existing expectations for the geometric nature of certain holographic states, as well as states arising from many-body systems that are generally believed to be non-geometric. We find that generic low-lying energy states of a free fermion system are indeed dominated by non-geometric contributions as indicated by our algorithm. However, there also exists states in these systems that have small relative error, such as the configurations following a global quench. It is suspected that the enhanced symmetries for the globally deformed states contributes to a more geometrical reconstruction. Moreover, quench dynamics reveal large-scale patterns of geometric evolution, despite the states themselves being non-geometric as characterized by the relative error. Further advances in the tensor Radon transform and its range characterization may help us understand what kind of properties in the boundary data lead to good reconstructions, and the nature of the large-scale behavior revealed by the numerical transform.

Discussion
In this work, we took some first steps in addressing the explicit bulk metric reconstruction problem in holography, as well as the question of whether non-geometrical bulk states can be detected from boundary data alone. Motivated by the tensor Radon transform, we provided -and explicitly implemented -an algorithm that reconstructs the bulk metric tensor without any a priori assumptions on the symmetries or form of the metric tensors, other than the fact that they are perturbatively close to AdS. We applied this reconstruction to entanglement data from holographic systems with large N , as well as to that of a 1d free fermion system. The reconstruction was also applied, time-slice by time-slice, to time-evolving states. We confirmed that, assuming the linearized Einstein's equations, thermalization following a global quench in holographic systems is consistent with a shell of in-falling matter in the bulk that finally settles into a state that appears to be gravitationally bound deep inside the bulk. This kind of behaviour is absent in free fermion systems where the corresponding shockwave of "in-falling matter" repeatedly oscillates between the boundary and the bulk due to the integrability of the system.
We also provided a partial answer to whether a given state in the conformal field theory has a well-defined geometric dual on the gravity side. We find that boundary reconstruction errors provide a quantitative measure that distinguished "geometrical" states, such as the ones we find in large N theories with semiclassical duals, from "non-geometrical" states like the generic excited states of a free fermion spin chain. This is a precise, albeit coarse, way of understanding what portion of the boundary data lies outside the range of the tensor Radon transform, and therefore cannot be interpreted as a tensor field on a hyperbolic background.
In the instances we have examined, this measure was an effective indicator of non-geometricality.
Finally, our initial attempts in developing this algorithm indicate that efficient numerical reconstructions of the bulk metric tensor from entanglement entropy data, at least at the linearized level, are achievable. This is partly because the optimization problem we consider is linear in nature and can be solved in polynomial time with low computational power. As such, it provides an efficiently computable tomographic procedure that translates boundary entropy data, where the underlying spacetime and gravitational dynamics are hidden, to a setting where it is manifest. Although a wealth of previous analytic results are available, the flexibility of numerical studies is also invaluable, as we have found in various areas of physics from quantum manybody systems to numerical general relativity. As we start to move away from the more tractable dynamics of quantum systems that carry simplifying assumptions such as symmetries, the full force of numerical methods will prove to be extremely helpful. Here we provide one such preliminary construction which can hopefully provide a stepping stone for future advances.
Further work is clearly required to improve upon our first efforts. Here we list but a few major directions where progress can be made.
On the front of new results in mathematics, especially related to tensor Radon transform: i To put our geometry detector on a firmer mathematical footing, one needs a rigorous characterization of the range of tensor Radon transforms on curved backgrounds, in particular, the hyperbolic background. In the language of holographers, we are in need of a set of necessary [38,39] and sufficient conditions that characterizes what type of quantum states have semi-classical dual geometries. A complete characterization of the range of Radon transform on hyperbolic space, for example, precisely provides such conditions that checks whether a set of entanglement data is dual to a semi-classical geometry close to AdS.
ii It is also crucial to obtain an explicit reconstruction formula in closed form. Such formulae are known for the Euclidean background, and for the case of scalar or vector Radon transform on curved backgrounds. Development of these results will further enable calculations in the continuum limit, which is relevant for AdS/CFT.
On the front of new results in physics: i We hope to extend this method to higher dimensions, which so far has faced the most obstacles in bulk metric reconstruction. Because area variation can already be cast as a tensor Radon transform in arbitrary dimensions [47], the idea of discretization followed by gauge fixing and linear optimization may simply be extended to minimal surfaces. Similarly, X-ray transforms in higher dimensions that makes use of correlation data can also be viable [23]. This may provide a numerically computable alternative to other methods [38,62,63].
ii Although we examine dynamical cases, we only reconstruct the spatial metrics for each individual time slice. Reconstruction of the full linearized space-time metric perturbation against the AdS background should also be possible by gluing together these spatial slices. However, one has to be careful about gauge fixing and the role played by the time coordinate. More work is needed to clarify these constructions.
iii It is also worthwhile to go beyond linearized level. A known approach is to apply the inverse Radon transform iteratively, such that the background geometry is updated using the reconstructed metric perturbation h ij . Such method is used in the geophysics community [53,64]. However, we have to be careful about the change in extremal surface beyond the linear level, especially in the dynamical cases. This is worthwhile though, since working beyond the linearized level should allow a better understanding of non-trivial topologies, for example. As we have seen in this work, they are invisible at the linearized level despite our expectation that collapse of bulk matter should lead to such changes in geometry.
iv Reconstruction of entanglement data from other quantum systems. Thus far, we have only seen limited application of this reconstruction algorithm, largely limited by the availability of entanglement data. Nevertheless, it may be possible to compute these entanglement approximately for smaller but more interesting quantum many-body systems using tensor networks or other classical techniques. It may also be worthwhile to obtain data for geodesic lengths via other more accessible data, such as correlation functions or mutual information.
v In light of the difficulty in computing von Neumann entropies for quantum systems, we can explore the possibility of using quantum Renyi entropies, which also have a geometric interpretation in holography [65]. A prominent advantage is its computability using numerics and measurability in actual quantum systems. Reconstructions using such data can enable us to directly acquire entropy data from quantum simulations and experimental setups.
vi Similar techniques using Radon transforms to recover the metric tensor of geometry emerged from entanglement are also applicable to near-flat manifolds outside the context of AdS/CFT. In fact, the Radon transform on flat-space is much better understood. Similar reconstructions should be possible for constructions like [18] where, for instance, the relevant quantum states can be low energy states of a gapped local Hamiltonian.
Finally, there is also much progress to be made in our reconstruction algorithm and related numerics.
i One area of improvement is a data-driven modelling. For our current work, we fixed a particular tiling that is easy to implement. However, it has been shown that a tiling depending on the boundary data can lead to improvement for bulk reconstruction [66]. These include Voronoi partitions or improved modelling using Bayesian inference [67].
ii Our discretization, constraint, and interpolation methods are all extremely simplistic, as appropriate for a proof of concept implementation. Numerous improvements can be made to improve the accuracy of the reconstruction algorithm, for example: proper triangular meshes, finite-element methods, better regularization techniques, etc. We expect that the reconstruction procedure here can be drastically improved in fidelity by adapting proper numerical techniques.
Here, let us formally define the geodesic tensor Radon transform. We begin with an introduction to the tensor Radon transform in general, but will quickly specialize to the special case of the 2-tensor Radon transform on the Poincare disk (with finite cutoff).
In short, the m-tensor geodesic Radon transform R m is a map which takes a symmetric m-tensor field defined on a sufficiently well-behaved Riemannian manifold M (with boundary) to the space of geodesics on that manifold.

A.1 General Definitions
Let (M, g) be an n-dimensional Riemmanian manifold with boundary ∂M and metric g. We say that (M, g) is a simple manifold if ∂M is strictly convex 9 and any two points in M are connected by a unique geodesic segment which depends smoothly on the endpoints [68]. Alternatively, a simple manifold is one in which the boundary is strictly convex, and where the exponential map exp p : exp −1 p (M ) → M is a diffeomorphism for every p ∈ M . Fixing some p ∈ M , we may identify a simple manifold with a strictly convex domain Ω of R n .
The simplicity of a manifold is a sufficient condition for the geodesic Radon transform to be well-defined [69]. There are more general conditions available, but simplicity will generally be sufficient for our purposes. From now on, unless otherwise stated, all of our manifolds will be assumed to be simple.
Let SM denote the unit circle bundle of M . The bundle SM is the collection of all pairs (p, v), where p ∈ M and where v ∈ T p M is a unit tangent vector at p. The boundary of the unit circle bundle, denoted ∂SM , consists of all such pairs where p ∈ ∂M . The boundary of the unit circle bundle naturally splits into two components, ∂ + SM consisting of all the inward pointing vectors, and ∂ − SM consisting of all the outward pointing vectors. We will define both components to be closed, i.e., vectors tangent to ∂M will be in both ∂ + SM and ∂ − SM .
Given (p, v) ∈ ∂ + SM , let γ p,v : [0, τ (p, v)] → M denote the unique unit speed geodesic through (p, v), i.e., the unique geodesic such that The parameter τ (p, v) denotes the exit time 10 of γ p,v , i.e., the first non-zero time such that γ p,v (τ ) ∈ ∂M . Note the exit time is well-defined under the assumption that the underlying manifold is simple.
Now, let f i1···im be a smooth, symmetric (covariant) m-tensor field on M . Then the Radon transform of f is defined by Thus, the Radon transform is a map R m : S m (M ) → C ∞ (∂ + SM ) which takes the space S m (M ) of smooth, symmetric (covariant) m-tensors on M to the space C ∞ (∂ + SM ) of smooth functions on the inward pointing boundary unit circle bundle component ∂ + SM . Since we can uniquely identify each (p, v) ∈ ∂ + SM with a corresponding unit speed geodesic γ p,v , it will be often convenient to consider the Radon transform as a map on the space of boundary anchored geodesics.

A.2 s-Injectivity
There are a few natural questions we may ask for the Radon transform, the foremost being the surjectivity and the injectivity of the transform. We will not comment much on the range of the tensor Radon transform, except to note that the tensor Radon transform is generally not surjective. This, of course, corresponds to the well known fact that not all boundary states have a well-defined bulk dual. A useful analytic characterization of the range remains an open problem for the tensor Radon transform on generic manifolds.
Let us now consider the problem of injectivity. The tensor Radon transform has a natural kernel. Let (M, g) be an n-dimensional simple Riemannian manifold. We will let dV denote the canonical volume form, which is locally given by where (x 1 , · · · , x n ) is some oriented chart, and where |g| is the determinant of the metric g ij in that chart. We will let ∇ denote the Levi-Civita connection on M . Now, let S m M denote the space of (covariant) symmetric m-tensors on M . We will define the inner derivative 11 d : where σ denotes complete symmetrization. In local coordinates, we simply have where parentheses indicate the complete symmetrization of the contained indices as usual. We will likewise define the divergence 12 δ : where Tr m,m+1 denotes (Riemannian) contraction between the mth and (m + 1)th arguments. In local coordinates, we have The operators d and −δ are adjoint for compactly supported symmetric tensor fields on M . More generally, for any compact region D ⊆ M , we have where u and v are (sufficiently smooth) symmetric tensors of the appropriate orders, i ν denotes interior multiplication with respect to an outward pointing normal ν, and where dS is the induced volume form on ∂D. The inner product ·, · on S m M is given by complete contraction, i.e., u, v = g i1j1 · · · g imjm u i1,··· ,im v j1,··· ,jm .
The significance of the operators δ and d are as follows. Let H k (S m M ) denote the Sobolev space of msymmetric tensors, i.e., the space of all sections which are k-times (weakly) differentiable, and such that all derivatives are locally square integrable. Each H k (S m M ) can be given the structure of a Hilbert space when M is compact. Then we have the generalized Helmholtz decomposition as follows: Theorem 1 (Generalized Helmholtz Decomposition [41]). Let (M, g) be a compact Riemannian manifold with boundary. Let k ≥ 1 and m ≥ 0 be integers. Given any section f ∈ H k (S m M ), there exists uniquely determined f s ∈ H k (S m M ) and v ∈ H k+1 (S m−1 M ) such that where δf s = 0, and v| ∂M = 0.
The fields f s and dv are called the solenoidal and potential parts of f .
Note that in the case m = 1, Theorem 1 is just the usual Helmholtz decomposition after identifying vectors and covectors using the metric.
The decomposition of a symmetric tensor field into solenoidal and potential parts gives us a natural identification of the kernel for the Radon transform. Indeed, we have the following result: 11 Despite the confusingly similar notation, d is not the exterior derivative d, which acts on forms and not symmetric tensors. Unfortunately, the notation is somewhat well-established in the integral geometry community, and so we will stick to it. We will never need to use the exterior derivative in this paper, but nevertheless we will denote the inner derivative with boldface font d as a reminder that it is not the exterior derivative. 12 Again, δ is not the co-exterior derivative. We will use the same boldface convention.
Therefore, Theorem 2 identifies a natural kernel for the Radon transform, namely the space of all potential tensor fields. Since every sufficiently smooth tensor field can be decomposed uniquely into a potential and a solenoidal part, a natural question is whether the solenoidal part is uniquely recoverable from the Radon transform, i.e., whether the space of potential tensor fields exhausts the kernel of the Radon transform. If this is indeed the case, i.e., if R m [f ] = 0 implies f s = 0, then we say that the Radon transform is s-injective.
The question of the s-inectivity of the Radon transform is a fundamental problem in integral geometry. The general case remains open, although the case for 2-dimensional simple manifolds was settled in the affirmative [70]: Theorem 3 (s-injectivity [70]). Let (M, g) be a simple 2-dimensional Riemannian manifold. Then the tensor Radon transform R m on M is s-injective for all m ≥ 0.
Note that in the case of the scalar transform for m = 0, all scalar functions are automatically solenoidal, so the scalar Radon transform R 0 is injective in the usual sense.
Given the s-injectivity of the Radon transform, we can recover the bulk tensor field by imposing the solenoidal gauge condition The solenoidal gauge is the most commonly employed gauge condition for the tensor Radon transform, but it comes with some inconvenient features. In particular, it does not respect the decomposition of a tensor into its trace and traceless parts. The 2-tensor Radon transform on a purely trace bulk tensor field is identical to the scalar Radon transform on the corresponding trace function, but the inversion of the 2-tensor Radon transform under the solenoidal gauge introduces extraneous gauge degrees of freedom which causes the recovered field to disagree with the original. This is rather undesirable, since the corresponding scalar transform is purely injective and admits a unique recovery. We will instead employ an alternative gauge condition, which we introduce in section C, that makes the tensor Radon transform consistent with the scalar Radon transform for pure trace fields.

B The Radon transform on the Poincare disk
Let us now focus our attention to the case of a single time-slice of AdS 3 . Such a time-slice is isomorphic to the hyperbolic plane, which we will consider in the Poincare disk model. Concretely, the Poincare disk is the Riemannian manifold of constant scalar curvature R = −1, defined within the open unit disk We will cover the Poincare disk using its natural Cartesian coordinates (x, y), or equivalently, using complex coordinates z = x + iy. The metric is then given by where |dz| 2 = dz · dz = dx 2 + dy 2 . We will denote the Poincare disk by H = (D, g). Geometrically, the geodesics of the Poincare disk are circular segments which are orthogonal to the boundary circle S 1 .
Both for physical reasons, and to properly define the Radon transform, we cannot work with the Poincare disk in its entirety. 13 Rather, we must impose a cutoff. We will let κ ∈ (0, 1) be a cutoff radius, and consider the cutoff disk Then we consider the cutoff Poincare disk to be the manifold defined by H κ = (D κ , g).
We will be mainly interested in the Radon transform for symmetric 2-tensors. The Poincare disk (with cutoff) is a simple manifold. As such, the Radon transform is always well-defined on the Poincare disk. To that end, let f ij be a symmetric 2-tensor field on H κ . The Radon transform is then given by Note that since the boundary of the (cutoff) Poincare disk is a circle, we may uniquely identify each geodesic γ p,v with a connected subregion A (i.e., a circular arc) of the boundary. The corresponding geodesic is then the minimal surface with respect to that boundary subregion. In this way, we may regard the image of the Radon transform as a function on the space of boundary subregions.
The metric on the Poincare disk is an isotropic metric. Let us write We may then write a unit tangent vector as v = e −λ (cos θ, sin θ), where θ is the natural angular coordinate on the unit disk. Then we can decompose the Radon transform as where in the last line we've defined the components The components h 0 , h, and h will be called the trace, holomorphic, and anti-holomorphic parts of f , respectively. Note that in this way, we can always identify any symmetric 2-tensor f ij with a functionf ≡ f ijγ iγj of harmonic content (−2, 0, 2).
We can then equivalently write the tensor Radon transform in the h-components as Note that the trace component of the transform simply gives the scalar transform (after absorbing the metric into the definition of h 0 ). Therefore the tensor Radon transform R 2 is essentially equivalent to the scalar Radon transform for tensors which are purely trace.
As previously mentioned, the solenoidal gauge does not respect the decomposition of the tensor into trace and traceless parts. It will therefore be convenient to find a gauge which preserves this decomposition. We do so in Appendix C.

C The Holomorphic Gauge
In this section, we present an alternative to the solenoidal gauge, which we call the holomorphic gauge, which both uniquely fixes a solution space to the inverse Radon transform and preserves the the scalar part of the transform.
To define the holomorphic gauge, we first introduce some background. We employ global isothermal coordinates (x, y, θ) on the unit circle bundle SH of the Poincare disk, where (x, y) are the usual Poincare coordinates on the base manifold, and where θ is a fiber coordinate indicating the angular direction of a tangent vector.
We can define the geodesic flow X, given in global isothermal coordinates on the unit circle bundle by where we write the metric as g = e 2λ (dx 2 + dy 2 ). Writing the vector field in terms of complex coordinates, we can decompose the geodesic flow as X = η + + η − , where where ∂ = 1 2 (∂ x − i∂ y ) and ∂ = 1 2 (∂ x + i∂ y ) are the Wirtinger derivatives. Let Ω k = L 2 (SH) ∩ ker(∂ θ − ik) be the space of square-integrable functions on the unit circle bundle with fixed harmonic content k. Then it can be shown [71] that η ± are smooth elliptic differential operators such that η ± : Ω k → Ω k±1 for any k ∈ Z. In particular, note that ∆ ≡ η + η − = η − η + is also a smooth elliptic partial differential operator.
Note that the geodesic flow X is naturally related to the Radon transform as follows [41]: Let f ij be a symmetric 2-tensor field, and let us define the function where γ denotes the unique unit speed geodesic through (x, y), with initial angle θ, and where τ + (γ) denotes the exit time of the geodesic. Note that we have u| ∂−SM = 0 and by construction. If γ (x0,y0,θ0) is a geodesic through (x 0 , y 0 , θ 0 ), then where θ (x0,y0,θ0) (t) ≡ arg γ (x0,y0,θ0) (t) denotes the angle of the tangent vector to γ (x0,y0,θ0) at time t.
Differentiating with respect to t, we getγ Note that the left-hand side is precisely the expression Xu. 14 Denoting by u f the unique solution to the transport equation (58) with boundary condition u f | ∂−SM = 0, it follows that the Radon transform R 2 [f ] is given by 14 Let us explicitly calculateθ here. We start with the geodesic equation Evaluating the x component, for example, we get d ds (e −λ cos θ) + e −2λ ∂ 1 λ(cos 2 θ − sin 2 θ) + 2e −2λ ∂ 2 λ sin θ cos θ = 0.
The transport equation can therefore be seen as the differential form of the Radon transform.
The key result leading to the holomorphic gauge is then the following theorem. Theorem 4. For any symmetric 2-tensor f ∈ L 2 (SM ), there exists a unique 2-tensor h ∈ L 2 (SM ) such that R 2 f = R 2 h, and such that h is of the form where h 0 ∈ L 2 (M ) ∩ Ω 0 , and where h ±2 ∈ ker η ∓ ∩ Ω ±2 .
Proof. We adapt the proof from Theorem 1 of [48], which covers the case where the underlying metric is the usual flat Euclidean metric. Let f ∈ L 2 (SM ) ∩ Ω k be given. Then consider the differential equation where v ∈ H 1 (SM ) ∩ Ω k−1 . Since ∆ is a smooth elliptic operator, it follows from the standard theory of elliptic differential equations that the above system admits a unique (weak) solution v. Given such a solution, let us define g = f − η + v. It follows that This shows that each f ∈ L 2 (SM ) ∩ Ω k can be written uniquely as where v ∈ H 1 (SM ) ∩ Ω k−1 , v| ∂SM = 0, and where g ∈ L 2 (SM ) such that η − g = 0.
where η − g 2 = 0. Since v 1 | SM = 0, it follows that we have (u + v 1 )| ∂+SM = u| ∂+M = R 2 f . If we define the 2-tensor h by where h 0 = f 0 − η − v 1 and h −2 = f −2 , then h satisfies R 2 h = R 2 f , and is such that η − h 2 = 0. We may repeat this reasoning with the f −2 term (using the complex conjugate of the previous decomposition) to obtain the desired result.
Definition 5. We will define a 2-tensor whose components satisfy the conditions of Theorem 4 to be in the holomorphic gauge.
Importantly, let us note that if a 2-tensor f is purely scalar, i.e., f = f 0 , then it is trivially already in the holomorphic gauge. Thus the holomorphic gauge is a gauge which respects the scalar part of the transform. This is to be contrasted with the solenoidal gauge, which will introduce spurious off-diagonal components even for scalar tensor fields.
Since h ±2 ∈ Ω ±2 , let us write h 2 = he 2iθ and h −2 = he −2iθ , where h, h ∈ L 2 (M ). In this notation, the holomorphic gauge condition reads We have so the holomorphic gauge conditions simplify to the Schrodinger type equations We can solve these equations by introducing an integrating factor of e 2λ . Then the equation for h gives us = ∂(e 2λ h).
This amounts to saying that the solutions are holomorphic functions (up to an exponential factor) on the unit disk D. We can therefore generate any solution as follows: Given any function on the unit circle h| ∂D : ∂D → C, we can extend h into the unit disk using Cauchy's integral formula to get where we write c λ = e 2λ(1) to denote the (constant) value of e 2λ on the unit circle. We will use this convenient reconstruction property of solutions in the holomorphic gauge to benchmark the numerical (inverse) Radon transform in Appendix E.

D The Numerical (Inverse) Radon Transform
In this appendix, we describe the details of the numerical (inverse) Radon transform.
The manifold we work with is the hyperbolic plane H 2 , which is modeled as the Poincare disk, i.e., the unit disk equipped with the canonical hyperbolic metric where (x, y) are global Poincare coordinates. For physical reasons, and to make the Radon transform welldefined, we must impose a cutoff on the Poincare disk. We thus pick a constant κ ∈ (0, 1) and work with the Poincare disk restricted to r ≤ κ. Equivalently, the cutoff Poincare disk can be regarded (after a rescaling of the metric) as a model for a hyperbolic plane with curvature equal to −κ 2 .
To perform the numerical Radon transform, we must first discretize the Poincare disk. A natural first choice would be to use a uniform tiling, a choice which conforms best to the intrinsic symmetries of the Poincare disk. However, the Gauss-Bonnet theorem places limitations on how fine a uniform tiling can be, and the inability to take the tile size to zero is an unwieldy restriction. Instead, we opt for a simple square tessellation which we perform in the Beltrami-Klein model.
The Beltrami-Klein model is related to the Poincare disk model through the change of coordinates Given a tessellation of the Poincare disk, we can discretize any function f by assigning to a given tile T the value of f at the centroid of T . Ordering the tiles arbitrarily, we can regard the discretized functions as vectors, We must also discretize the (ideal) boundary of the Poincare disk for the Radon transform. We do so by placing M equally spaced boundary sites on the unit circle. We will then consider the collection of all geodesics originating from one boundary site and terminating on another. We can order the boundary sites by their angular position on the unit circle, and the geodesics lexicographically by the angular positions of their endpoints.
For each geodesic of the background geometry, the integral (51) can then be discretized by replacing the functions h 0 , h and h with their piece-wise constant discretizations: where W 0 , W , and W contain the information on the geodesic and the remaining parts of the integrand. Explicitly, W 0 and W are defined by W (γ j , T ) = γj e −2λ(γj (s))+2iθ(s) 1 T (γ j (s)) ds, where 1 T (γ(s)) is an indicator function such that Note that W 0 is nothing but the arc lengths of the geodesic segments that intersects a tile T . The function W , on the other hand, also contains a complex weight which captures the directionality of the geodesic in the tile.
It's important to note that W 0 and W depend only on the particular choice of discretization, and can be pre-computed.
We can collect all of the quantities into a matrix equation. Let W be the K × 3N , where K = M 2 , defined by Likewise, let h be the length 3N vector defined by Then the discretized Radon transform, which we denote by ∆L, is given by the matrix equation The discrete inverse Radon transform is then just the inverse problem to the system (97). However, the inverse problem is complicated by the fact that the Radon transform is neither surjective nor injective (recall that the forward transform has a non-trivial kernel that corresponds to the gauge degrees of freedom).
To make the inverse well-posed, we must supplement it with a set of gauge constraints to pick out a unique inverse (when it exists) in the continuum case. In the discrete analog, we expect the kernel to be visible through the presence of zero eigenvalues in the singular value spectrum of W. Due to discretization and numerical errors, the system (97) will either be singular or extremely ill-conditioned. As in the continuum case, we will supplement the system (97) with a set of discretized gauge constraints to make the problem well-posed.
For our implementation, we discretize the holomorphic gauge conditions given by equations (80) and (81). The holomorphic gauge constraints are simple first order partial differential equations. They can be realized discretely by implementing ∂ x and ∂ y are finite-difference operators. We will use a simple three-point stencil for the finite-difference operators, although higher-order variations can be used for increased accuracy. For example, the x partial of a function f at tile i will be approximated by where L(i) and R(i) denote the indices of the tiles to the immediate left and right of tile i. 15 The y-partials are analogous. In this way, we can write discretize the partial derivatives as matrices∆ x and∆ y . 16 Since (99) 16 The entries of∆x are given by if tile i is a non-boundary tile, with the appropriate modifications for the boundaries.
the vector h effectively contains 3 copies stacked on top of each other, we define the operators ∆ y = ∆ y |∆ y |∆ y .
The discretized holomorphic gauge conditions can then be written as where Λ is the matrix defined by whereΛ has entries given byΛ

E Accuracy of the Numerical Reconstruction
In the absence of an analytic reconstruction formula in the continuum case, we need to validate the numerical reconstruction so as to provide confidence that it performs the inverse transform correctly. To do so, we will first benchmark the numerical reconstruction using various test cases obtained by instantiating known rank-2 symmetric tensor fields in the bulk in the holomorhpic gauge. From equation (86), we can see that a bulk solution in the holomorphic gauge can be readily generated by prescribing the boundary values.
As a benchmarking test, we therefore test the reconstruction algorithm as follows: 1. We first pick an arbitrary function h| ∂D : ∂D → C on the unit circle. We also pick another arbitrary scalar function h 0 : D → C. (86), we extend the function h| ∂C into the bulk to get a holomorphic function h : D → D.

Using equation
3. From the functions h 0 and h, we define the bulk tensor where h denotes the complex conjugate of h. The f defined in this way is a real-valued continuum bulk 2-tensor with components fixed in the holomorphic gauge.
4. We run the forward numerical Radon transform to generate boundary data b. From b, we then run the discretized inverse Radon transform with holomorphic gauge constraints to numerically recover a discretized bulk solution h.
5. Finally, we compare the values of h, evaluated at the centroids of the discretized tiling, with the reconstructed value h. Denoting the exact solution at the centroids as h * , we can evaluate the relative error to get a sense of the reconstruction quality.
In the absence of known analytic results, and before we move onto with physically relevant examples, we can run the above test with several choices of boundary functions h| ∂D as a proof of concept that the numerical inverse Radon transform is performing an adequate reconstruction. Figure 16: Top: Original bulk data generated with h 0 (x, y) = xy and h| C (θ) = sin(2θ). Bottom: The bulk data reconstructed after running the forward Radon transform to get boundary data. Visually, the two sets of data are identical. The relative error between the two are E bulk ≈ 0.0227 (relative error 0.004643 without boundary tiles).
Below, we show some sample test cases for the numerical inverse transform. All of the transforms shown below are performed with 1000 bulk tiles and 100 boundary sites (for a total of 4950 geodesics), following the procedure outlined above.
The reconstructions here show good agreement with the original bulk data. On a qualitative level, the plots of the bulk metric profile are essentially indistinguishable between the original and the reconstruction. The bulk relative errors across various test cases range from 0.01 to 0.1, indicating quantitatively successful reconstructions in general. The boundary relative errors are typically an order of magnitude smaller than the bulk relative errors, ranging from 0.001 to 0.01. The magnitudes of the boundary relative errors for these known test cases can serve as an estimate for the general magnitude of numerical errors present in the algorithm. See  It should be noted that the relative errors shown here can actually be slightly misleading. Most sources of error in the reconstruction arise due to large fluctuating values of the boundary tiles. Tiles at or near the boundary are generally underconstrained due to the relatively small number of geodesics which pass through any given boundary region. This allows the boundary tiles to take on arbitrary values in order to minimize the relative boundary error as the algorithm is designed to do, although this comes at the cost of bulk accuracy. We can see that if we exclude the the values of the boundary tiles from the calculation of relative error, that the relative error of the reconstruction is generally an order of magnitude smaller. This suggests that the bulk reconstruction performs very well deep into the bulk, with larger errors towards the boundary. This is important to keep in mind, since it suggests that the qualitative bulk picture provided by the numerical inverse should be generally trustworthy. Figure 17: Top: Original bulk data generated with h 0 (x, y) = 2e −(x 2 +y 2 ) and h| C (θ) = cos(2θ) + sin(3θ). Bottom: The bulk data reconstructed after running the forward Radon transform to get boundary data. Visually, the two sets of data are identical. The relative error between the two are E bulk ≈ 0.0273 (relative error 0.005516 without boundary tiles). Figure 18: Top: Original bulk data generated with h 0 (x, y) = 1/(1+x 2 +y 2 ) and h| C (θ) = cos(5θ). Bottom: The bulk data reconstructed after running the forward Radon transform to get boundary data. Visually, the two sets of data are identical. The relative error between the two are E bulk ≈ 0.0885 (relative error 0.01357 without boundary tiles). Figure 19: Top: Original bulk data generated with h 0 (x, y) = x 2 + y 2 and h| C (θ) = 2 cos(θ) + 3 sin(3θ). Bottom: The bulk data reconstructed after running the forward Radon transform to get boundary data. Visually, the two sets of data are identical. The relative error between the two are E bulk ≈ 0.01886 (relative error 0.003709 without boundary tiles).

E.1 Constrained Optimization
With the discrete Radon transform (97) and the discretized holomorphic constraints (102), we can solve this linear system for h as a constrained optimization problem: subject to Ch = 0, We look for a best-fit solution h * that solves the above system. In general, we do not expect there to exist a solution h * such that Wh * − b = 0, due to either numerical/discretization errors or the boundary data being non-geometric.
Because the objective function is linear, this problem has a unique global solution that can be obtained using standard constrained least squares. We briefly review the method [] below. Theorem 6. Consider the constrained least square problem (107). Assuming the stacked matrix is left-invertible and C is right-invertible, a vector h * uniquely solves the constrained least square problem (107) if and only if there exists some g such that Proof. Suppose that (h * , g) satisfies (109). Clearly h * satisfies the constraint Ch * = 0. Then for any h that satisfies the constraint Ch = 0, we have = ||W(h − h * )|| 2 + ||Wh * − b|| 2 (113) where in the second line we used the definition of the norm v 2 = v † v, in the third line we used the fact that and in the fourth line we used the fact that Therefore, h * is an optimal solution. Furthermore, the optimal constrained solution is obtained if and only if W C (h − h * ) = 0.
By assumption the stacked matrix is left-invertible, therefore the minimum h = h * is also unique.
In the reverse direction, it suffices to show the matrix in (109) is invertible. Suppose on the contrary the matrix is non-invertible. Then there must exist a vector (h, g) = 0 such that Left multiplying both sides by (h, g) † , we have Noting that Ch † = 0, we get Wh 2 = 0, so that Wh = 0. Since the stacked matrix (108) is injective, it follows that we must have h = 0. The system (118) then reduces to We then conclude that g = 0 since C is right-invertible by assumption. This is in contradiction with the assumption that (h, g) = 0. Hence the matrix in equation (109) must be invertible, and a solution for g must exist.

E.2 Interpolation and Regularization
In principle, the optimal solution h * can be obtained through straightforward matrix inversion of equation (109). Suppose h is a column vector with 3N entries, and the constraint matrix is M × 3N , then any well-known polynomial algorithm for matrix inversion is of O (M + 3N ) 3 . The inversion of system (109) is slightly complicated by the fact that the matrix W is generally expected to be extremely ill-conditioned. We can get around this by regularizing the system.
To make the the least squares problem (109) better conditioned, and to smooth out small scale fluctuations in the discretized reconstruction, we employ a derivative type Tikhonov regularization. We replace the matrix W † W in (109) by where γ > 0 is a regularization parameter which controls the strength of the regularization, and where ∆ x and ∆ y are the discretized partial differential operators defined in (101). This replacement effectively changes the least squares problem in (107) to subject to Ch = 0, which both regularizes the system so that the smallest singular values of W are of order √ γ and also takes into account the strength of fluctuations in the resulting solution. Since we expect small scale fluctuations to be mostly due to bulk discretization errors, this choice of regulator serves as a reasonable filter. In this note, all of our reconstructions employ regularization with γ = 10 −8 .
In the case of fixed data but variable number of bulk tiles, the reconstruction can also become ill-conditioned when the number of bulk degrees of freedom exceed the number of boundary constraints in the form of geodesic lengths. Roughly speaking, because we have O(3N ) number of bulk degrees of freedom, one for each tensor component at a particular tile. Suppose there are K sites on the boundary, then the number of linear equations from geodesic lengths is of order O(K 2 ). Hence the reconstruction can become ill-conditioned when 3N > K 2 . While it is possible to decrease the number of bulk tiles, or increase the number of boundary sites, both come at a cost depending on our requirements for reconstruction. As an alternative, we can also interpolate between geodesic data by adding virtual sites and the lengths of geodesics for the extended set of sites.
In the current implementation, we add the virtual sites in between the original lattice sites such that the new lattice scale is half of the original. Let us label all sites sequentially along the counter clockwise direction on the boundary circle from 1 through 2K such that 2K + 1 ≡ 1. To generate the geodesic lengths between a virtual site j and an original site i, we average the lengths of two geodesics that are anchored at (i, j + 1) and (i, j − 1). For a geodesic that ends on two virtual sites, i, j. We take the 4-point average of the original geodesic lengths for the ones anchored at (i, j + 1), (i, j − 1), (i + 1, j), and (i − 1, j).