Geometric Aspects of Holographic Bit Threads

We revisit the recent reformulation of the holographic prescription to compute entanglement entropy in terms of a convex optimization problem, introduced by Freedman and Headrick. According to it, the holographic entanglement entropy associated to a boundary region is given by the maximum flux of a bounded, divergenceless vector field, through the corresponding region. Our work leads to two main results: (i) We present a general algorithm that allows the construction of explicit thread configurations in cases where the minimal surface is known. We illustrate the method with simple examples: spheres and strips in vacuum AdS, and strips in a black brane geometry. Studying more generic bulk metrics, we uncover a sufficient set of conditions on the geometry and matter fields that must hold to be able to use our prescription. (ii) Based on the nesting property of holographic entanglement entropy, we develop a method to construct bit threads that maximize the flux through a given bulk region. As a byproduct, we are able to construct more general thread configurations by combining (i) and (ii) in multiple patches. We apply our methods to study bit threads which simultaneously compute the entanglement entropy and the entanglement of purification of mixed states and comment on their interpretation in terms of entanglement distillation. We also consider the case of disjoint regions for which we can explicitly construct the so-called multi-commodity flows and show that the monogamy property of mutual information can be easily illustrated from our constructions.


Introduction
An important lesson from the past decade's research program in holography is that quantum information theory provides a powerful tool to structure our thinking about quantum gravity and quantum field theories. By now, significant progress has been made in understanding how quantum entanglement is encoded in the bulk [1][2][3][4][5] and, in particular, how Einstein's equations emerge from its dynamics [6][7][8]. See [9][10][11][12] for useful reviews. According to the RT prescription [1], entanglement entropy of boundary regions is encoded by the area of certain codimension two surfaces in the bulk, in a manner reminiscent of the well-known Bekenstein-Hawking formula of black hole entropy. Due to its elegance JHEP05(2019)075 and simplicity, entanglement entropy has become a powerful tool to investigate some of the most fundamental aspects of the holographic correspondence, from bulk reconstruction to the emergence of spacetime. Accordingly, an enormous amount of work has been devoted to finding both analytic and numerical solutions of minimal surfaces in different gravitational settings.
Very recently, Freedman and Headrick suggested that the computation of entanglement entropy can be reinterpreted in terms of a specific convex optimization problem [13]. More specifically, they showed that the problem can be cast in terms of finding a bulk vector field V with maximum flux through the corresponding boundary region, divergenceless ∇·V = 0, and with an upper bound set by the Planck scale |V | ≤ 1/4G N . 1 Besides the enormous conceptual advantage that this new picture brings into the bulk interpretation of different quantum information properties, it also presents an alternative method with potential utility for the computation of holographic entanglement entropy of boundary regions.
From a more practical standpoint, one can use the tools borrowed from Riemannian geometry and apply them to the construction of the auxiliary objects -referred to as flows or bit threads -of the max-flow min-cut proposal. Convex optimization techniques have been shown to be useful in the proof of statements about such objects [16]. However, explicit examples of such vector fields do not make a significant appearance in the entanglement entropy literature up to now. An obvious reason for this is the fact that such vector field configurations are not unique and therefore have no physical meaning on their own. Nevertheless, it is illustrative to work out examples of such objects (representatives of the homology class of vector fields that satisfy the previously described conditions) to gain intuition about their behavior and hopefully be able to make more general claims that are independent of the specific class of vector fields used in the constructions.
In this work we will start filling this gap by constructing explicit examples of such vector fields configurations, starting from situations in which the minimal surface is known. Even though this is a crucial restriction, the main goal of this work is to initiate a broader program that aims at the understanding of more general bit threads configurations. Indeed, in all of our constructions, we are able to see a large amount of non-uniqueness as it is expected from the previous arguments. By focusing on particular examples and studying the freedom allowed in the constructions, we then are able to uncover a large class of new possibilities for consistent bit thread configurations.
In summary, the kind of constructions that we study throughout this paper fall in one of the following three categories: • Symmetric or geodesic flows: these are obtained from general geodesic foliations of a suitable (possibly unphysical) bulk geometry. The geodesics here are identified as the integral lines of the corresponding vector field, while its norm is chosen such that the divergenceless condition is satisfied locally.
• Maximally packed flows: these are obtained based on the entanglement wedge nesting property of holographic entanglement entropy. The integral lines of the vector field JHEP05(2019)075 are taken to be orthogonal to the geodesics (or minimal surfaces) associated to a set of nested subregions. By construction, the norm of the vector saturates the bound |V | = 1.
• Mixed flows: these are obtained by consistently cutting and pasting symmetric and/or maximally packed threads in different patches of the bulk geometry.
The last category is the most generic case that we consider, however, we emphasize that there are infinitely many more configurations that are not covered by our constructions.
In the future, it would be interesting to study more in detail this degree of ambiguity, e.g, by determining the most general kind of transformations that can be implemented to a particular set of threads and determine how would they modify our specific vector field configurations. 2 This paper is organized as follows. In section 2 we propose a general algorithm to construct flows in cases where the minimal surface is known. Specializing to the case where the integral curves are taken as geodesics (or minimal surfaces) of the underlying geometry, we then present several explicit examples of the symmetric flows advertised above, including spheres and strips in vacuum AdS, and strips in a black brane geometry. In section 3 we study general constraints on curvature and matter fields that must be fulfilled in a generic gravitational setting such that a similar construction based on geodesic foliations leads to a valid flow. Based on this study, we revisit the explicit examples constructed in section 2 and explain why they work based solely on the geometric properties of the underlying geometry. In section 4 we describe the construction of maximally packed flows in a general setting and explain the reasoning behind these constructions. Finally, in section 5 we present our two main applications, which are based on mixed flows: specific flow configurations that illustrate the computation of entanglement of purification and the monogamy property of mutual information. We close in section 6 with a summary of our results and a discussion about potential topics for future research.

Symmetric flows: a general algorithm
In [13] it was proposed that the entanglement entropy of a holographic theory associated to a spatial boundary region A is given by the maximum flux of a bounded, divergenceless, bulk vector field through the region in question. As noted there, there are infinite many vector fields that maximize that flux and respect the defining properties. Indeed such non-uniqueness play an important role in its tempting interpretation as a holographic representation of the pattern of spatial entanglement of the field theory degrees of freedom.
In this section, we would like to consider a particular construction of a representative of the homology class of vector fields that solve the above problem. We will start by proposing a general method before specializing to the AdS case. For this construction we JHEP05(2019)075 assume that the following data is given: a space-time metric g ab with asymptotic boundary at z → 0, i.e. g ab → z −2 η ab , a connected boundary region A and its corresponding minimal area bulk surface m(A) (such that m(A)| z→0 ∼ ∂A). The method has two parts. First, we propose a family of integral curves (also referred to as flow lines or threads) with certain properties (discussed in detail below), from which we extract their tangent vectorτ . And second, we find the magnitude |V | by imposing the divergenceless condition. With these two ingredients, we can then construct the desired vector field V = |V |τ .
Let us explain in more detail the method. The first step is to propose a family of integral curves that satisfy the following properties: 1. They must be continuous and not self-intersecting so that the tangent vectorτ is unique and well defined everywhere. Notice that we do not need to impose that the integral curves foliate the full manifold. Points that are not reached by any of the curves will have by definition a vanishing tangent vectorτ = 0.
2. At the minimal surface m(A), the tangent vectorτ must be equal to the unit normal of the minimal surfacen m . This is required since the bound on the magnitude of the vector field must be saturated there, a key fact in the equivalence between the minimal surface and the maximum flux prescriptions for the holographic entanglement entropy.
3. And finally, since the vector field must be divergenceless, the integral curves should start and end at the boundary.
These three conditions are completely generic, however, since the choice of integral curves is highly non-unique we will impose a further constraint to reduce the arbitrariness in the possible vector fields, without eliminating it completely: 4. We will consider curves that leave from region A and end in regionĀ. Physically, this condition is equivalent to restrict our attention to microstates without entanglement between degrees of freedom contained in region A orĀ.
Once a family of integral curves with the above conditions has been found, the second step is to find the appropriate norm |V |. A natural way to do it is by enforcing that the flux through an infinitesimal area element is constant throughout the threads. This construction automatically ensures that the divergenceless condition is satisfied. In practice, we proceed in the following way: first, we parametrize the curves X(x m , λ) by the point at which they intersect the minimal surface, x m , and a parameter λ that runs along each curve. Then, we follow the set of integral curves that leave from an infinitesimal region δA(x m ) around the point x m at m(A). This area will propagate along the threads and define a co-dimension one bulk region with the topology of a cylinder, which we denote by T (δA(x m )). Since the surface of T (δA(x m )) is made out of integral curves, this means that there is no flux leaving or entering that surface. Therefore, the divergenceless condition can be imposed by choosing |V | such that the flux through any transverse section of T (δA(x m )) is constant -see figure 1 for a schematic representation. Mathematically, this condition implies where h| λ is the induced metric on the plane orthogonal to the flow line at the point X(x m , λ). Using the fact that at the minimal surface |V (x m , λ m )| = 1, and our knowledge of the exact location of that surface, one can write an explicit expression for |V | along the integral curves, this is where λ m is the parameter at which the flow line intersects the minimal surface m(A).
The final non-trivial check is to verify that the magnitude of the vector field respects the upper bound away from the minimal surface, |V | ≤ 1. If this is not the case then, the choice of integral curves should be modified and the process should be repeated until the bound is achieved. In that sense, the critical/essential step of our construction relies on the choice of what we call "good" integral curves.

Explicit constructions in pure AdS
For concreteness, we would like to start by illustrating the method for the case of pure AdS d+2 in Poincaré coordinates and for highly symmetric cases: spherical regions (2.1.1) and strips (2.1.2). In cases like these the natural choice is to consider curves that respect the symmetry, so the problem simplifies and becomes effectively lower dimensional.

The sphere
The sphere is the most symmetric region on can think of. The minimal surface as well as the explicit bulk and boundary modular Hamiltonians are known explicitly and have been widely exploited leading to many deep results connecting properties of entanglement entropy in field theory with the dynamics of the emergent bulk spacetime.

JHEP05(2019)075
Minimal surface: consider the minimal surface associated to a spherical region A of radius R located at the boundary of AdS d+2 . In Poincare coordinates (a constant-t slice of) the metric is given by 3 The minimal surface is given by a hemisphere of radius R that extends into the extra dimension z ∈ R + . The minimal surface can be defined implicitly by the collection of points (r m , z m ) that satisfy where r m = i (x i m ) 2 and x i m are the standard cartesian coordinates. Notice that, without loss of generality, we have chosen to locate our minimal surface centered at the origin of the Poincaré coordinates.
For later convenience, we give here the outward-pointing unit normal vectorn m at a point (r m , z m ), which in the above coordinate system is given bŷ For simplicity we have chosen to suppress all angular coordinates, so the index a runs over the coordinates (r, z).
Integral curves: since the problem has spherical symmetry, we restrict our attention to one plane, i.e., we set x 1 = r and x i = 0 ∀ i > 1. The angular dependence can be easily restored at the end by implementing a rotation on the resulting integral curves. 4 Now, let us consider the space of geodesics that lie on the (r, z) plane and find the ones that intersect m(A) at (r m , z m ) and whose tangent vectorτ is parallel to the normal n m at that point. We will argue below that such geodesics are a natural candidate for the integral curves of our vector field V .
The two-dimensional effective metric is the following: 6) and the set of geodesics that lie in this plane is given by the two-parameter family of circumferences defined implicitly by where r s is the center of the circle on the r-axis and R s is its radius. The tangent vector with unit norm at an arbitrary point is given bŷ (2.8) 3 We will set the AdS radius to unity throughout this paper, L = 1, but it can be easily restored at any point by dimensional analysis. 4 In addition, notice that since there is a reflection symmetry in this plane, x1 → −x1, we can further restrict our attention to x1 = r ∈ R + .

JHEP05(2019)075
Enforcing thatτ =n m at a point (r m , z m ) on the minimal surface leads to Finally, plugging (2.9) into (2.7) we obtain an implicit expression for the family of geodesics orthogonal to m(A), parametrized by the point r m on the minimal surface. Before proceeding further, we need to check if the proposed flow lines satisfy our criteria. Based on our construction, the first condition (namely, that the curves do not intersect with each other) is the only non-trivial requirement. In order to check this, we parametrize the curves with respect to the point r a at which they intersect A, Similarly, we can find the dual point rā at which the curves intersectĀ, Self-intersection is avoided if and only if the curves parametrized by r m are nested. Since the proposed curves are geodesics, this condition is guaranteed provided that dr a /dr m > 0 and drā/dr m < 0. Indeed, a quick calculation leads to which satisfy the above conditions. Therefore our choice of integral curves is validated.
Vector field: we can now proceed to find the appropriate norm of the vector field |V |.
In order to compute the right hand side of (2.2) we consider the family of integral curves parametrized by the point r m at which the curves cross the minimal surface m(A). These curves are given implicitly by equation (2.7), with R s (r m ) and r s (r m ) given in (2.9). Then, we compute the orthogonal metric at different points along the curve, whereτ is the unit tangent vector, given in (2.8). A brief calculation leads to is convenient to eliminate z (in order to avoid multi-valuedness), and choose r = λ as the affine parameter. We therefore solve for z(r) in (2.7), which can be done analytically: On the other hand, from implicit differentiation of (2.7) we obtain where in the last equality we have used the explicit expressions for R s (r m ) and r s (r m ) given in (2.9). Putting all together, and restoring the angular dependence, we find that The magnitude of the vector field follows from equation (2.2), which leads to where z m = z(r m ). Finally, we would like to express our vector field as a function of (r, z) without reference to the minimal surface. This can be achieved by solving for r m = r m (r, z) from the geodesic (2.16) and plugging it back into (2.19). A short calculation leads to (2.20) It is easy to check that |V | ≤ 1 everywhere, an is only saturated at the position of the minimal surface m(A). This confirms that our construction was successful. Finally, we can plug the function r m (r, z) into the expression forτ (2.8), to obtain τ a = 2Rz The full vector field V = |V |τ is then given by In figure 2 we show plots of the vector lines as well as the magnitude of V for d = 1, d = 2 and d = 3 spatial dimensions.

The strip
Given the success achieved in the spherical case, it is natural to conjecture that geodesics would always be good candidates for the integral curves of an arbitrary entangling region. However, it is easy to check that geodesics do not always satisfy the minimal requirements JHEP05(2019)075 outlined above. Moving beyond spherical symmetry, the simplest class of subsystems that can be studied analytically are strips. Indeed, it is possible to show that for strip geometries it is not possible to build a family of non-intersecting geodesics which are perpendicular to the minimal surface for dimensions other than d = 1 and d = 2. A proof by explicit construction is given in appendix A. The goal of this section is to propose a family of integral curves that satisfy all the condition for strips in arbitrary dimensions.
Minimal surface: the strip is defined by the collection of points with − /2 ≤ x 1 ≤ /2 The corresponding minimal surface can be obtained analytically by exploiting the translational symmetry of the problem, see e.g. [17], and it is given by where z * is given in terms of the strip length by . (2.24) The two signs ± in (2.23) correspond to two different branches of x m (z m ); the '+' corresponds to the x > 0 part of the surface, while the '−' corresponds to x < 0. However, since there is a reflection symmetry around the origin, x → −x, we will restrict our attention to the right branch, with '+' sign and x > 0. We note that the minimal surface for the strip (2.23) is more complicated than the minimal surface for the sphere (2.4). In particular for general d we cannot invert (2.23) to obtain z m = z m (x m ). We will see below that this represents an obstacle for finding a closed expression for vector field V (x, z). Finally, we give the outward-pointing unit normal vectorn m at a point (x m , z m ) on m(A), (2.25) As a check, notice that if we set d = 1 and z * = R we recover (2.5).

JHEP05(2019)075
Integral curves: the first step in our construction is to propose a family of integral curves that satisfy the consistency requirements outlined at the begging of section 2. As mentioned above the geodesics are not a good choice for strips in general dimensions. A second natural proposal is to consider curves in the (x, z) plane based on minimal surfaces for strips. This new choice is equally intuitive since these surfaces satisfy all the symmetry requirements and can be defined covariantly from the given bulk geometry. More specifically, our proposed family of integral curves will be given by the curves of the form where x m (z) is a minimal surface of the form (2.23), centered at x = 0 and with a maximum depth z * . It is clear x s corresponds to a shift along the x-direction, while α rescales the depth, The choice of sign ± corresponds to the two different branches of x m (z); the '+' corresponds to the x > x s part of the surface, while the '−' corresponds to x < x s . It is important to keep both signs in this case in order to have a full coverage of the region x > 0. It is easy to check that the above curves are geodesics in the following two-dimensional effective metric: 5 where the extra d in the denominator and the constant V ⊥ ≡ d−1 ⊥ arise from integrating out the transverse coordinates.
In section 3.2 we will justify the above choice of integral curves in more detail. In particular, we will show that a foliation with minimal surfaces generally leads to a good family of integral curves for a strip entangling region in any translational invariant bulk geometry. For the time being, we will proceed with the explicit construction of the vector V in the particular case in consideration, i.e., empty AdS. The next step is to impose the orthogonality of the integral curves with respect to m(A). In order to do so we first compute the unit tangent vector along the curves, Notice that both (2.29) and (2.25) are normalized with respect to the original AdS metric. The '−' here correspond to the x > x s branch, while the '+' corresponds to x < x s . It is clear that we need the latter sign to enforce orthogonality, since we are focusing on the x > 0 portion of m(A) and the outward-pointing normal vectorn m for this branch (2.25) has this sign. Enforcing thatτ =n m at a point (x m , z m ) on the minimal surface leads tõ This two-dimensional effective metric has the form of a hyperscaling violating metric. Indeed, hyperscaling violation is known to arise generically from the dimensional reduction of higher dimensional AdS spaces [18].

JHEP05(2019)075
Notice that if we set d = 1, z * = R,z * = R s and x s = r s we recover the expressions in (2.9). Finally, plugging (2.30) into (2.26) we can obtain an explicit expression for the family of geodesics orthogonal to m(A), parametrized by the point z m on the minimal surface: where x m (z) is the function given in (2.23) (with '+' sign) and z ∈ (0,z * ). Again, we have two branches depending on the sign of the last term in (2.31). The '−' branch intersects A at while the '+' branch intersectsĀ at With these expressions we can now check if the curves parametrized by z m are properly nested. Since the curves are minimal surfaces, this condition is guaranteed provided that dx a /dz m < 0 and dxā/dz m > 0. Indeed, after some algebra we find that and which justify the choice of minimal surfaces as good integral curves for the case of the strip.
Vector field: we can now proceed to find the appropriate norm of the vector field |V |, through equation (2.2). In order to do so, we first compute the orthogonal metric at different points along the curve parametrized by z m , JHEP05(2019)075 Putting all together, and restoring the transverse coordinates, we find that The magnitude of the vector field follows from equation (2.2), which leads to 6 The subscript in x(z m ; z) refers to the choice of sign in (2.31). In particular, it is important that x − (z m ; z) appears in the numerator of (2.40) since this is the branch that intersects the minimal surface m(A). The two signs of x ± (z m ; z) in the denominator of (2.40) cover the regions x > x s and x < x s , respectively, and together they span the whole axis x ∈ R + . The explicit expressions for |V | as a function of z and z m are straightforward to obtain but are lengthy and not particularly illuminating, so we will refrain from writing out these results here. The final step would be to solve explicitly for z m (x, z) from the geodesic equation (2.31) and replace the result in equation (2.40) to obtain |V | purely in terms of x and z. As mentioned above, this is not possible to do analytically for d > 1. However, we can perform a parametric plot by varying z and z m . The final result is plotted in figure 3 for strips in various dimensions. From these plots we can easily check that the magnitude of the vector is indeed bounded, |V | ≤ 1, and is saturated precisely at the location of the minimal surface. 6 We note that the same result can be obtained by performing the same analysis starting directly with the effective metric (2.28). This includes: normalizing the vectorsn andτ with respect to this metric and using it as a background metric in (2.36). We will come back to this point in section 3.2.

The strip
In terms of minimal area surfaces, going beyond pure AdS is challenging and often require the use of numerical techniques. Of particular physical relevance is the case of minimal area surfaces in an AdS black hole geometry, from which one can compute entanglement entropy in an excited state at finite temperature. Perhaps the only non-trivial example that can be treated analytically is the case of strips in a black brane geometry which was recently explored in [19]. In d = 1 the problem simplifies drastically. On one hand, the codimension-two surfaces required to compute entanglement entropies are space-like geodesics. On the other hand, Schwarzschild-AdS 3 (the BTZ black hole) is diffeomorphic to AdS 3 so the geodesics can be obtained in a compact form. In the following we will specialize for simplicity to the d = 1 case. However, the generalization to higher dimensional cases is straightforward, since the minimal surfaces are known analytically for any d [19]. We expect that the general method and the qualitative results will apply also to these cases.
Minimal surface: consider the minimal surface associated to an interval of length in (a time slice of) a BTZ black hole, The metric is invariant under translations in x, so without loss of generality we focus on an interval centered at the origin. The minimal surface is given by the collection of points (x m , z m ) that satisfy 42) or, equivalently, In either of these expressions, z * indicates the maximum depth of the geodesic. This parameter is related to the length of the interval through any of the following equivalent relations: The outward-pointing normal unit vectorn m at the point (x m , z m ) is given bŷ

JHEP05(2019)075
For future reference we also give an explicit expression for the slope s m = n z m /n x m of the normal vector, at a given point of the minimal surface Integral curves: now lets consider the space of geodesics that lie on the (x, z) plane that intersect m(A) at (x m , z m ) and whose tangent vector is parallel to the normaln m at that point. These geodesic will represent the integral curves of our vector field V . As we will see below, for the case of the BTZ black hole, some integral curves will necessarily end up at the horizon, while others will go back to the boundary. For this reason it will be convenient to express the geodesics in two different parametrizations, as z(x) and as x(z). The former parametrization will be useful to describe threads that start and end at the boundary, while the latter will be more convenient for threads with one end at the horizon.
• Parametrization I: a general geodesic that goes through the point (x 0 , z 0 ) with a slope dz/dx = s 0 can be written as follows: where the constants C i are given by As a simple check, notice that if we set x 0 = 0, z 0 = z * and s 0 = 0, we recover the minimal surface (2.42). Enforcing that the tangent to the geodesic is equal to the normal of the minimal surface (2.45) at the point (x m , z m ) we arrive to the corresponding family of integral curves. These can be obtained simply by setting (2.46). Furthermore, using the equation for minimal surface one can eliminate either x m or z m and write the final result as z(x m ; x) or z(z m ; x).
• Parametrization II: in the second parametrization, a general geodesic that goes through the point (x 0 , z 0 ) with a slope dx/dz = 1/s 0 can be written as follows: where D and σ are the following constants:

JHEP05(2019)075
The solution has two branches, depending on the sign of σ. For s 0 > 0 the correct solution has σ = −1, while for s 0 < 0 the right choice is σ = 1. As a check, notice that if we set x 0 = 0, z 0 = z * and s 0 = 0, we recover the minimal surface (2.43). Enforcing that the tangent to the geodesic is equal to the normal of the minimal surface (2.45) at the point (x m , z m ) we arrive to the corresponding family of integral curves. These can be obtained simply by setting A comment about the parametrizations is in order. As mentioned above, the parametrization II is more convenient for geodesics that reach the horizon. However, if the geodesic has the two endpoints at the boundary one can still use the two branches in equation (2.48) to fully describe them. Notice that the fact that we need two solutions is to be expected because x(z) is double valued for this type of geodesics. The extra branch in each case describes the analytic continuation of the geodesic after the maximum depth is reached which, by symmetry, is located at the center of the interval. Moreover, at this point one has that dz/dx = 0 (or, equivalently, dx/dz → ∞). Similarly, for geodesics that reach the horizon one can alternatively use equation (2.47) from the parametrization I to describe them, but in that case one would need to truncate the solution at z = z h . Beyond this point the solution would be unphysical. Before proceeding further, let us investigate the exact conditions that must be satisfied in order to have threads connecting the boundary and the horizon. Assuming that s 0 = 0, from our solution (2.48) we obtain that the point at which dx/dz → ∞ is given by: For z max < z h the curve has its two endpoints at the boundary, but for z max ≥ z h one of its endpoints inevitably reaches the horizon. The transition takes place exactly when z max = z h , which leads to a simple relation between the maximum slope and the radial depth z 0 , For |s 0 | < s max the curve has its two endpoints at the boundary, however, for steep enough slopes |s 0 | > s max the curve will necessarily have one of its endpoints at the horizon. The above bounds hold for any constant-t geodesic in the BTZ background (2.41). On the other hand, using explicitly the information about the minimal surface (2.42) and the equation for the normal slope (2.46) we can arrive to the following bounds for threads of our vector field V :

JHEP05(2019)075
which translates into For points of the minimal surface m(A) close to the center |x| < x max (z > z max ) the normal vector is steep enough so that |s m | > s max , hence the corresponding threads reach the horizon. Conversely, points of m(A) near the edges, with 2 > |x| > x max (z < z max ), have threads that connect the region A with its complementĀ. This phenomenon is illustrated in figure 4. Finally, in order to see if these geodesics are indeed good candidates for integral curves we need to check that they are properly nested. There is a subtlety, however, since some of the threads end at the horizon. Without loss of generality, we focus on the region x m > 0 and split the analysis in two, 0 < z m < z max and z max < z m < z * : • 0 < z m < z max : in this range the threads connect a point of region A, with a point in the complementary regionĀ, The subscript in x(z m ; z) in these expressions refers to the choice of σ or, equivalently, the branch of the geodesic. By the property of entanglement wedge nesting, we know that these integral curves are properly nested if and only if dx a (z m )/dz m < 0 and dxā(z m )/dz m > 0. Indeed, a brief calculation leads to: and The only term that is not explicitly positive definite is the denominator of the first term in the square brackets, however, it is easy to check from The above two inequalities then show that the corresponding integral curves are nested.
• z max < z m < z * : in this range the threads connect a point of region A, x a (z m ) given in (2.54), with a point at the horizon,

JHEP05(2019)075
Besides the inequality (2.56) it is possible to show that All terms in this expression are strictly positive, except for the piece of the denominator in the square brackets. However from (2.52) one can check that if z m > z max then z 2 In order to prove the nesting of the integral curves we point out that the excised geometry outside the bulk horizon has been conjectured to be dual to a pure state with extra degrees of freedom living at the (stretched) horizon [20,21]. 7 In this setup, geodesics that end at the horizon can indeed be interpreted as entanglement entropies in the purified state [22], so from the property of entanglement wedge nesting, we can conclude that (2.56) and (2.59) are indeed sufficient to prove the nesting of the corresponding integral curves.
Vector field: we can now proceed to find the appropriate norm of the vector field |V |, through equation (2.2). To do so, we find it convenient to work with the second parametrization, so we label the geodesics as x(z m ; z). The orthogonal metric at different points along one of these geodesics is whereτ is the corresponding unit tangent vector. In this parametrization, this is given bŷ (2.61) Sinceτ z = 0 precisely at z =z * , this means thatz * gives the maximum depth of the geodesic labeled by z m . Plugging (2.61) into (2.60) we arrive to This line element can be rewritten in terms of the transverse element dz m by implicitly differentiating the geodesic x(z m ; z), which yields The magnitude of the vector field can be read from equation (2.2), which leads to Again, the subscript in x(z m ; z) refers to the choice of σ or the branch of the geodesic. In particular, it is important that x − (z m ; z) appears in the numerator of (2.64) since this is the branch that intersects the minimal surface m(A). The two signs of x ± (z m ; z) in the denominator of (2.64) cover the full axis x ∈ R + . The explicit expression for |V | as a function of z m and z is straightforward to obtain but is lengthy and not particularly illuminating, so we will not write the final result here. The final step would be to invert the geodesic equation x(z m ; z) to obtain z m (x, z), and replace the result in equation (2.64) to obtain |V | purely in terms of x and z. However, this is not possible to do for the case in consideration. Here we proceed as we did for the strip geometry in higher dimensions, i.e., by varying z and z m and performing a parametric plot. The final result is shown in figure 4. From these plots we can easily check that the magnitude of the vector is indeed bounded, |V | ≤ 1, and is saturated precisely at the location of the minimal surface.

Constraints on geometry and matter
In the previous section, we provided a general algorithm to construct bounded, divergenceless vector fields, with maximum flux through a given boundary region, for a given bulk geometry. With this algorithm studied specific examples of vector fields V representing holographic bit threads, i.e., satisfying ∇ µ V µ = 0 and |V | ≤ 1, for some particular cases of interest: spheres and strips on empty AdS in arbitrary dimensions and strips on an AdS black brane geometry. The purpose of this section is to relax some of the conditions used there and investigate the possibility of implementing a similar construction in more general bulk geometries. Let us first summarize the main ingredients of our constructions: • Perhaps the most essential feature of our examples was the high degree of symmetry in both, the geometry of the problem (including the bulk geometry and the shape of the entangling region A) and the vector field V . Indeed, assuming that the symmetries of the problem are inherited by the vector field V , our analysis became effectively two dimensional. For spheres, we studied vector field configurations by focusing on JHEP05(2019)075 a given plane and assuming that V had no angular components nor angular dependence. Similarly, for strips we assumed that V had no components along the tangent directions so we were also able to reduce our analysis to a fixed plane.
• Another important ingredient was our canonical choice of integral curves. In all cases, the integral curves were taken to be geodesics of a suitable auxiliary geometry. For spheres such curves were given by the geodesics of the global metric, while for the case of strips they were given by geodesics of a 2-dimensional effective metric. Once the choice of integral curves was made, the final step was to obtain the norm |V | from the continuity equation (2.2). From this equation it was clear that |V | ≤ 1, if and only if the area transverse satisfies h(x, λ) ≥ h(x, λ m ) everywhere away from the minimal surface. We emphasize that this condition was verified only a posteriori, but was not explicitly imposed in our constructions.
In this section we will investigate more generally when the bound |V | ≤ 1 is satisfied for this kind of constructions. More specifically, we will consider the same class of foliations in a generic bulk geometry and find constraints on the curvature in order to satisfy the aforementioned bound. We do so by studying the transverse area of the hypersurfaces orthogonal to the given foliations and finding the requirements that must be satisfied such that the transverse area element increases as one moves away from the minimal surface. In addition, we will use the full Einstein equations to translate these constraints in terms of the stress-energy tensor of the bulk matter fields. On general grounds, we expect that such constraints must hold for (semi-)classical gravity in the large N limit (G N → 0), but should be modified once quantum corrections in the bulk are taken into account.

General geodesic foliations
Lets consider a manifold Σ with boundary ∂Σ endowed with a Riemannian metric g ab , having the property that given a minimal co-dimension one bulk surface m(A) anchored at the boundary ∂Σ, this is m(A) ∂Σ = ∂A such that A ⊂ ∂Σ is a connected boundary region, there exists a global or partial foliation of the manifold by geodesics that intersect orthogonally the surface m(A). 8 If such a foliation exists, then, it is possible to prove a theorem for the area of the family of hypersurfaces that are orthogonal to the foliation.
For the class of foliations that we consider, the most general form of the metric is the following: generality, we take λ = 0 to be the location of the minimal surface, so that |λ| measures the proper length from m(A) to an arbitrary slice of the foliation. We further assume that λ > 0. This means that we generically require two different patches to cover the regions inside and outside of m(A). Let us write our vector field V as V = |V |ξ, and assume that ξ is affinely parametrized. Thus, in the above coordinate system ξ takes a simple form: We would like to find a coordinate invariant criterion for the transverse area element induced by ξ. Namely, we would like to find an equation for where γ = det γ ij . Indeed, since we know that at the minimal surface ∂ λ ln √ γ = 0, what we really want is to argue that ∂ λ ln √ γ > 0 away from the minimal surface so that |V | decreases. We note that this condition may be too strong in certain situations. Namely, there can be cases where ∂ λ ln √ γ < 0 locally, but |V | < 1 globally nevertheless. For simplicity, we will assume that |V | decreases monotonically along our geodesic foliations. We notice that our problem can be recast as an electrostatics problem with a constant surface charge density at the location of the minimal surface. In flat space, the magnitude of the electric field indeed decreases monotonically away from the sources, given that i) the transverse area increases with the distance and ii) by Gauss's Law, the magnitude of the electric field is inversely proportional to the transverse area. Here, we want to find an equivalent statement for the transverse area in a general Riemannian metric g ab . In order to do so, we will impose the condition that ∂ 2 λ ln √ γ > 0 everywhere so that the magnitude |V | decreases monotonically. 9 We begin by studying how the curvature of the ambient space g ab varies along the geodesic flow. This is, we compute the Ricci tensor R ab and study its λλ component, Using the explicit expressions for Γ a bc it is easy to check that: leading to the following relations:

JHEP05(2019)075
where in the last line we have used that γ ik ∂ λ γ il = −γ il ∂ λ γ ik . Equation (3.8) can be rewritten in the following way, which is the Riemannian version of the Raychauduri equation. Alternatively, this equation can be put in a coordinate invariant way by noticing that the extrinsic curvature tensor and its trace are given by K ij = 1 2 ∂ λ γ ij and K = γ ij K ij = ∂ λ ln √ γ, respectively, which lead to the more familiar form, We can now study the geometry of the constant-λ slices. The induced metric on these hypersurfaces is the following: (3.11) It can be shown that Ricci tensorR ij associated to this metric satisfies the following relation: which relates its components with those of the global Ricci tensor and the extrinsic curvature K ij . To arrive at this equation one can start from the definition of R ab in the coordinates (λ, x i ) and then separate the components involving only γ ij , ∂ k γ ij which will form theR ij . The remaining terms will involve functions of the extrinsic curvature K ij .
With this relation at hand we can rewrite the Raychauduri equation (3.10) as Finally, adding (3.10) and (3.13) we obtain (3.14) whereR = γ ijR ij is the Ricci scalar associated to the induced metric (3.11). This equation has the advantage that involves only geometric quantities.
In the following, we will use equations (3.10), (3.13) and (3.14) to argue for the monotonicity of the transverse area in various cases of interest.

Strips in a general translationally invariant background
For strip entangling surfaces in a general translationally invariant background the problem can be formulated in terms of an effective two-dimensional bulk geometry, which can be obtained by dimensional reduction of the original metric. We assume that this effective geometry is foliated by geodesics so that

JHEP05(2019)075
where γ is the determinant of the induced geometry for the constant-λ slices. A geometry of this type can be obtained by integrating out the transverse coordinates, and then foliating the effective geometry with geodesics, which are nothing but codimension-one minimal surfaces in the ambient space. Now, since the constant-λ hypersurfaces are one dimensional in this case, a number of simplifications take place: i) the extrinsic curvature has only one component K ab = K xx and ii) the induced Ricci tensor vanishes identicallyR ij = 0. 10 Therefore, we can use the general formula (3.14) in our analysis, which in this case reduces to Rewriting this equation purely in terms of γ we obtain, Therefore, if R < 0 we obtain that ∂ 2 λ √ γ > 0. Since ∂ λ √ γ = 0 at the minimal surface this implies that ∂ λ √ γ > 0 everywhere so that √ γ is a monotonically increasing function of λ, i.e. it increases as one moves away from the minimal surface.
In the last part of this section, we will revisit the vacuum solutions studied in section 2 and then study more general backgrounds supported by matter fields. In the latter case, we will be able to translate the curvature condition found above in terms of an upper bound for the energy density of the fields.
Vacuum solutions: let us start by revisiting the known cases, strips in pure AdS and in a black brane background. As it was shown in the previous section, the two dimensional effective geometry for the former case is the following: The factor of d in the denominator and the constant V ⊥ ≡ d−1 ⊥ arise from integrating out the transverse coordinates. One can easily repeat the same exercise for a strip in an AdS black brane, which yields the following effective metric: It is easy to see that (3.18) is a special case of (3.19), where one sends z h → ∞. So without loss of generality we can focus only on the second case. A short calculation shows that the Ricci scalar associated to the metric (3.19) is:

JHEP05(2019)075
Since it is negative definite, then it follows that |V | decreases monotonically away from the minimal surface. Therefore |V | ≤ 1 everywhere for strips in pure AdS or in an AdS black brane in arbitrary dimensions.
Non-vacuum solutions: to begin with, consider a generic static gravitational system in 3-dimensions. Let us assume that the full bulk metric is a solution of Einstein's equations with an appropriate stress-tensor: 11 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: Notice that the right hand side is parametrically large in the G N → 0 limit. Indeed, in this limit the negative curvature is fully supported by the cosmological constant term, so a flow based on geodesic foliations would be allowed for arbitrary matter content. On the other hand, if one takes G N to be small but finite, equation (3.24) would indeed provide a sharp upper bound for the energy density. Once this bound is violated, geodesics start to focus and the transverse area start to decrease instead of increase. It is possible to generalize this inequality to higher dimensions. To show this, we dimensionally reduce a translationally invariant metric in 3 + k dimensions down to 3 dimensions. We write the full metric as follows: where ds 2 3 is of the form (3.21). The Einstein-Hilbert term in the action transforms as follows: where κ 2 3 = κ 2 /V ⊥ . Assuming that the 3 + k dimensional metric satisfies Einstein's equations (3.22), after some algebra we arrive to the following relation: 27) 11 We have set 8πGN = 1 for simplicity.

JHEP05(2019)075
where R is the induced Ricci on a constant-t slice of the effective 3-dimensional geometry. Hence, we have that R < 0 if and only if the local energy density is bounded from above: The last term in this equation is contracted using the full metric (3.25) so the criterium is fully covariant. This term is strictly negative, so in higher dimensions there is an interplay between the cosmological constant term and the kinetic term of the scalar field that arises from the dimensional reduction. In particular, we notice that at finite G N the bound becomes more rigid as one increases k.

Spheres in a general rotationally invariant background
For spherical entangling surfaces one can try to repeat the same steps as for the case of strips. However, since the construction here is directly based on geodesics of the full geometry (instead of minimal surfaces), it is easy to see that it is not possible to recast the problem solely in terms of a two dimensional effective geometry. Nevertheless, we will see below that the spherical symmetry is powerful enough to allow some simplification.
To begin the analysis, we start from a generic metric of the form: We emphasize that the two-dimensional metric ds 2 2 cannot be obtained by a dimensional reduction of the full geometry, as was done in the case of strips, but instead should be thought of as an auxiliary object that is part of the full metric (3.29). With the above ansatz, one can compute the quantities of interest, namely the curvature tensors R ij ,R ij and the extrinsic curvature tensor K ij . After some algebra, we can rewrite (3.14) as follows: In this equation γ ij involves the metric functions of the r-part of the metric as well as the angular components; moreover, R 2 refers to the induced metric of the auxiliary 2dimensional metric defined in (3.29). We must require the sum of these terms to be positive, in order to have a transverse area that grows away from the minimal surface. This implies that for spheres, the condition in terms of curvature generalizes to: As expected, for k = 0 we recover the 2-dimensional condition (3.17) which implies negative curvature. One way to realize the above condition for general k is to assume that R 2 < 0 and simultaneously have all terms in the right hand side of (3.31) to be positive. As we will see below, this is indeed the case for spheres in empty AdS. We leave a more general analysis of (3.31) and its geometrical and physical interpretation to future work. In the last part of this section we will revisit the case of spheres in empty AdS and offer some comments on more generic cases.

JHEP05(2019)075
Vacuum solutions: let us revisit the examples studied in section 2, namely, spheres in empty AdS for an arbitrary number of dimensions. In this case, it can be shown that the spatial part of the metric can be put into the more symmetric form so that e 2τ (λ,r) = r 2 γ(λ, r) . (3.33) To prove this we start with the spatial part of the metric of an hyperbolic AdS black hole: Indeed, this metric covers the inside of the entanglement wedge associated to a spherical region in empty AdS. Since the 2-dimensional metric (at constant angles) is invariant under translations in u, it is clear that the vertical geodesics with u = constant correspond to the integral curves of V in this coordinate system. Finally, the coordinate transformation brings the metric (3.34) into the form (3.32). Since and then, it follows immediately that the right hand of (3.30) is positive. Hence, the transverse area grows monotonically away from the minimal surface. The other example in the vacuum that one can consider is the case of an AdS black brane. However, we do not know the minimal surface associated to a spherical region in this case, nor the explicit spacelike geodesics. It is also not clear that the foliation admits a simplification of the form (3.32). It would be interesting to attempt such a construction numerically.
Non-vacuum solutions: adding time to the metric (3.29) we get One can assume that this metric satisfies Einstein equations (3.22), and proceed in the same way as we did for the case of strips. By doing so, one obtains an upper bound on the energy density ε(λ, r) in terms of the cosmological constant and derivatives of the metric functions, γ(λ, r) and τ (λ, r). However, the final result is not fully covariant and does not seem to be particularly illuminating. Hence, we will refrain from writing out these results here.

JHEP05(2019)075
4 Nesting property and maximally packed flows So far we have considered only one kind of thread configurations, which are based on specific foliations of the spacetime by geodesics or minimal surfaces. This particular class of constructions ensures that the flux across the minimal surface is maximal, but they generically lead to a flow that decreases in magnitude as one moves away from the minimal surface. In certain cases, the entanglement wedge itself might contain one or multiple bottle necks that might admit more flux than the one supplied by the geodesic flow. The idea here is to propose a strategy to simultaneously maximize the flux across the minimal surface and another minimal surface of interest. In order to do so, we will use a general property of holographic entanglement entropy, which is known as entanglement wedge nesting. Let us begin for simplicity with an arbitrary two dimensional Riemannian geometry with a boundary; this can be, for instance, a constant-t slice of pure AdS 3 or a BTZ black brane. Consider a family of intervals of lengths l n , such that their left boundary it located at x L = 0 while the right boundary is at x R = l n , such that l n+1 > l n with n ∈ {0, 1, · · · , N }. This family will have minimal surfaces m(A n ) such that their entanglement wedges r(A n ) are nested, i.e., r(A n+1 ) ⊃ r(A n ). Now, the nesting property for bit threads implies that one can find a flow that maximizes the flux through A N , while simultaneously maximizing the flux through all A n with n < N . If one takes l n to vary continuously then the above construction would then generate a thread configuration with magnitude |V | = 1 everywhere in r(A N ) \ r(A 0 ), which smoothly interpolates between the unit normal on m(A 0 ),n 0 , and the one on m(A N ),n N . Since the norm of vector saturates the inequality |V | ≤ 1 in this region, we call this class of constructions maximally packed flows. 12 Let us make the above example more precise. Assuming that the geometry is given by a time slice of AdS 3 and taking the intervals to have endpoints at: then, the corresponding family of minimal surfaces will be given by: i.e., a family of semicircles with both, radii and centers given by l n . The continuum limit of this family can be easily obtained by taking the limit N → ∞ and defining a parameter λ ≡ n/N ∈ [0, 1], such that: Taking this limit leads to the continuous foliation by semicircles: For this foliation, the outward-pointing unit normal vectorn(λ) at a point (x, z) of a constant-λ slice is given by: Since V is equal to the normal vectorn(λ) in the region r(A N )\r(A 0 ), it would be desirable to eliminate λ from from this equation and obtain an expression that does not refer to a specific slice of the foliation. We can do so by explicitly inverting l(λ) from equation (4.4). A brief calculation yields Plugging (4.6) into (4.5) we then arrive to the following formula for the vector field, V a ≡n a = z x 2 + z 2 x 2 − z 2 , 2xz . (4.7) Notice that since V has been constructed from the unit normal vectorn n associated to a set of minimal surfaces, then, the divergenceless condition of V follows trivially. The argument is the following: at every minimal surface m(A n ) one can define its extrinsic curvature K n = ∇ an a n . However, the minimality condition implies that the extrinsic curvature vanishes identically, therefore ∇·V = 0. Since V is smooth in the interpolating region, divergenceless and satisfies the bound |V | ≤ 1, then it represents a valid thread configuration.
By construction, then, we have that the portion of the vector field V that goes from m(A 0 ) to m(A N ) is equal to the unit normal vectorn n at the intermediate minimal surfaces. We illustrate this example graphically in figure 5. A few important observations are in order:

JHEP05(2019)075
• The integral curves can be easily obtained from the vector field itself (4.7). For the above example they satisfy the following equation: which has solutions of the form: These curves describe circles centered at x c = 0, z c = c with radius c.
• Since the volume of space is infinite near the boundary, the flux across any of the regions A n is also infinite. Therefore, we have to be careful with the regularization. A particular choice of cutoff can be obtained if one imposes that the flux across all surfaces are the same. In this case the cutoff surface should follow one of the integral curves defined by (4.9), i.e., it must be orthogonal to all slices of the foliation. If we assume that the cutoff at the A 0 slice is , then a brief calculation shows that In other words, since the area of nested intervals grows monotonically with the length, it means that the cutoff should also increase so that the flux is preserved. However, other choices of cutoff can be made by extending this construction, but will generally lead to different fluxes across the surfaces. For example, for the standard choice (x) = constant, the flux increases as the region increases.
• Notice that so far we did not say anything about what V should be in the region inside m(A 0 ) or the region outside m(A N ). One simple option is to glue part of a geodesic flow constructed in the previous sections. We do this explicitly in figure 5.
Since the magnitude and tangent vectors across the minimal surfaces are continuous, the vector field V is continuous and once differentiable at the gluing surfaces.
• We can generalize the above construction to arbitrarily nested intervals. For a two dimensional geometry, the most general case has endpoints at x L = l L n and x R = l R n such that l L n+1 < l L n and l R n+1 > l L n . These two conditions guarantee that the corresponding intervals are properly nested. Notice that this simple generalization already lead to infinite many new possibilities, which are inequivalent in the continuum limit. Fixing one endpoint uniquely determines the flow, because all choices of l(λ) would be equivalent to each other by considering reparametrizations of λ. Conversely, having two arbitrary functions l L (λ) and l R (λ) make the choices inequivalent, up to a function that continuously maps points in the range of l L (λ) to points in the range of l R (λ). As an example, consider a time slice of AdS 3 . Repeating the steps of the previous example leads to a vector field of the form: Having two functions l L (λ) and l R (λ) related through only one equation, namely, the minimal surface equation, it is clear that generically it is not possible to explicitly solve for V (x, z). This is clear from the previous argument, since we generically expect that different choices of functions would lead to different vector fields V . For particular choices of l L (λ) and l R (λ), however, it is still possible to invert (4.12) to obtain λ(x, z). Having done that, then one can simply replace the result in (4.11) to obtain an explicit equation for V (x, z). Figure 6 illustrates one of these generalized maximally packed flows.
Needless to say, the above constructions can be easily generalized to arbitrary 2dimensional geometries, and to higher dimensional cases for any family of nested subregions. In the following section, we will use both, maximally packed flows and geodesic flows, to study more general thread configurations for particular cases of interest.

Entanglement of purification
An interesting byproduct of the bit threads picture of holographic entanglement entropy was the discovery of the holographic dual to the quantum information quantity known as entanglement of purification [27]. Already in [13], the authors studied the minimal cross section for disjoint regions and interpreted it as the maximum flux among the thread configurations that live completely inside the corresponding entanglement wedge, connecting the two disconnected regions. Using intuition from tensor networks, this minimal cross section was later identified with the concept of entanglement of purification in [20,21].

JHEP05(2019)075
Interestingly, the minimal cross-section of the entanglement wedge is an example of a bottle neck, discussed in the previous section. Hence, we expect that a particular construction with the so-called maximally packed flows can be designed to specifically compute this quantity. We will devote this subsection to construct and study this kind of flows.
The definition of entanglement of purification is the following. Consider a quantum system Q bipartitioned into sets of degrees of freedom A and B, in a state described by a density matrix ρ AB . If the state is mixed, the 'entanglement entropies' associated to A and B differ from each other, S A = S B . More importantly, these entropies quantify both quantum and classical correlations between A and B. In order to quantify the amount of correlations that are purely quantum, one option is to purify the system, i.e., to consider a set Q of additional degrees of freedom, and a choice of pure state |ψ for the overall system QQ , such that Tr Q |ψ ψ| = ρ AB . If we further partition the auxiliary system Q into A and B , we can compute the entanglement entropies S AA = S BB , which should indeed arise from purely quantum correlations. By optimizing among all possible purifications and all possible partitions A and B , the entanglement of purification between A and B is then defined as According to the proposal of [20,21], the holographic dual of this quantity is given by the minimal cross section Σ of the entanglement wedge of AB, namely Hence, the resulting density matrix ρ AB describes a mixed state. It is well known that, the entanglement wedge of the system AB undergoes a first order transition. When A and B are close enough r(AB) is connected, however, when A and B are sufficiently far r(AB) is disconnected. More specifically, the transition depends on the cross ratio and it is only for z ≥ 1 that the corresponding entanglement wedge is connected. An example of a connected r(AB) is shown in figure 7. Since the intervals are symmetric, the minimal cross-section Σ AB in this example is simply given by the green vertical line centered at the origin. For the disconnected case, the minimal cross-section would be identically zero.
The second example shown in figure 7 consists of a single interval in a black brane background A = [a 1 , a 2 ] and its complement B =Ā. In this case, one starts directly with a mixed state and, in particular, ρ AB is given by a thermal density matrix. There are two bottle necks in this case, the two vertical lines at the edges of A, Σ  Figure 7. Holographic entanglement of purification for two particular cases of interest: two disconnected intervals in pure AdS and an interval and its complement in a black brane background. The corresponding entanglement wedges are shaded in yellow, while the candidates for minimal cross-section Σ AB are shown in green. In the first case, there is only one bottle neck in r(AB), for the case connected case and hence one candidate for Σ AB . For large enough separations, the entanglement wedge undergoes a first order transition, and becomes disconnected. In this case the entanglement of purification vanishes. In the second case there are two bottle necks in r(AB) and hence two candidates for the minimal cross-section Σ AB . The entanglement of purification is this case is the minimum of the two, which depends on the ratio between the length of the interval and inverse temperature. minimal surface that computes entanglement entropy of region A, Σ AB . The entanglement of purification is then the minimum of these two, It can be shown that the two above areas exchange dominance at a particular value of the interval length , in units of the inverse temperature β, so that for large enough /β the disconnected solution is favored. In that case the thread configuration that computes the entanglement of purification has an interesting interpretation, when compared with the one that computes the entanglement entropy of region A. If one recalls the interpretation of the threads as quantum bits of holographic entanglement, then it means that the threads that go through the surfaces Σ (1) AB compute the maximum amount of entanglement entropy present in S(A) whose source is purely quantum mechanical. In other words those threads could be related to the maximum number of Bell pairs which can be distilled fromρ AB using only local operations and classical communication (LOCC). The threads that go into the horizon, on the other hand, could be interpreted as the minimum amount of correlations present in S(A) that are thermal or classical. The above separation makes sense since, as we will see below, one can explicitly construct flows which simultaneously compute the entanglement entropy S(A) (as the flux through the boundary region A) and the entanglement of purification E(A : B) (as the part of the flux that goes through Σ (1) AB ). We will now proceed to construct the vector fields that computes the entanglement of purification in these two examples. We will start with one interval in a BTZ black hole, since most of the formulas needed for this case were already worked out in section 2.2. Then, we will consider the case of two disjoint intervals in AdS. The constructions of this example will in turn set the grounds for the topic that we will discuss next, namely, the monogamy property of mutual information.

One interval in a BTZ black brane
Consider an interval of length centered at the origin, and its complement, in a constant-t slice of the BTZ geometry (2.41). This example was studied in detail in section 2.2 so we will borrow some of the formulas presented in that section. There are two bottle necks for this configuration, as shown in figure 7. The first one consists of two straight surfaces located at x = /2 and x = − /2, respectively. They have area: where is a UV regulator and β is the inverse temperature β = 2πz h . The second bottle neck corresponds to the minimal area surface that computes entanglement entropy for region A. This minimal surface can be written either as (2.42) or (2.43), depending on the parametrization, and is characterized by a single constant z * that indicates the maximum depth of the surface. This constant is related to the length of the interval through any of the equivalent expressions presented in equation (2.44). The area of this surface can be obtained by direct integration and can be written in terms of the interval length as follows Area Σ (2) Comparing the two areas, we arrive to the following expression: where c = 3/2G is the central charge of the 2-dimensional CFT and c ≡ β π log( √ 2 + 1) ≈ 0.28β (this yields a critical radial depth z c * = z h / √ 2 ≈ 0.71z h ). Thus, for large enough regions, the surface Σ (1) AB gives the minimal cross-section and hence the entanglement of purification, 13 while for small regions Σ (2) AB computes both entanglement entropy of region A and its entanglement of purification. In this latter case, then, there should be particular microstates where S(A) can be interpreted entirely as quantum correlations between Bell pairs in region A and its complement B. In the former case, on the other hand, one finds that there is a maximum number of Bell pairs that can be distilled from the mixed state, so at least part of S(A) must be thermal in nature.
In the following, we will construct explicitly flows that maximize the fluxes through A and Σ (1) AB simultaneously, for > c . In other words, we will find thread configurations that most efficiently avoid the horizon (since those would necessarily cross Σ (1) AB ) when maximizing the flux on A and whose magnitude saturates the bound |V | = 1 at Σ (1) AB . Interestingly, one can see that the maximally packed flows studied in the previous section provide a clean construction of such configurations.

JHEP05(2019)075
Quantum and classical entanglement in S(A): we would like to construct a vector field that computes E(A : B) for an interval of length > c . Without loss of generality, we will discuss the construction of the x > 0 portion of the vector field, given that the configuration has reflection symmetry around the origin.
We start by applying the same construction of section 4 based on the nesting property of entanglement entropy. In order to do so, we take a family of geodesics that interpolate between Σ (1) AB and m(A) and are nested in the appropriate sense. Here, we recall that the geometry outside the bulk horizon has been conjectured to be dual to a pure state with extra degrees of freedom living at the (stretched) horizon [20,21]. In this setup, geodesics that end at the horizon can indeed be interpreted as entanglement entropies in the purified state [22], Moreover, the nesting property applies in the usual sense, but one must in turn consider the horizon as part of the boundary of spacetime. With these points in mind, then, we can consider the family of geodesics with the right endpoint fixed at x R = /2 and either (i) a second endpoint at the horizon with z h ∈ [ /2, −∞] or (ii) a left endpoint at x L ∈ [−∞, − /2]. In order to continuously parametrize this family of geodesics, we will use the notation of (2.48). We set z 0 = z * so that we can vary s 0 ∈ [−∞, 0] continuously, and use the branch with σ = 1. Furthermore, we set JHEP05(2019)075 AB ). The extra threads that we have added in the central bundle cross the minimal surface at points x < x e and eventually reach the horizon. Hence, these threads have a purely thermal interpretation.
Before finishing this example, let us offer more insights on the interpretation. Since we have obtained both S(A) and E(A : B) from the same vector field configuration, then it is reasonable to subtract these two quantities. The difference corresponds to the smallest amount of flux that goes thorugh the horizon and therefore has a natural interpretation as the minimal part of S(A) whose nature is undoubtedly thermal, this is

Two intervals in pure AdS
The second example illustrated in figure 7, consists of two strips in a pure AdS geometry. As for the BTZ case, we have a situation where there is a bottle neck in the entanglement wedge associated to ρ AB , and therefore we can use the nesting property of entanglement entropy to construct a thread configuration that simultaneously computes entanglement entropy of one of the regions, say S(B), and the entanglement of purification E (A : B). For the sake of simplicity we will focus on AdS 3 , in which case the strips are just 1-dimensional intervals. A similar construction would also work for higher dimensions and we expect the results to be qualitatively the same.
We will begin by considering two disjoint and sufficiently close intervals, so that the corresponding entanglement wedge is connected. As we will see below, we will need a slight generalization of the method used in the BTZ example. Namely, the set of nested intervals needed to construct the maximally packed flow will need both endpoints to vary in a delicate synchronized fashion, as was considered at the end of section 4 (see figure 6). We will show, however, that the geodesic flows constructed in section 2 provide an appropriate family of nested boundary regions needed for this example. We will conclude by considering the case of two adjacent intervals, which can be easily obtained by a limiting case of the disconnected configuration.
Disjoint intervals: for clarity of exposition we will discuss only the symmetric case, where the two intervals are of equal size 1 = 2 = , and are separated a distance ∆ from center to center. 14 When the separation is small enough, the entanglement wedge is a connected surface whose boundary consists of two concentric semi-circles of radii R 1 , R 2 with R 2 > R 1 , such that = R 2 − R 1 and ∆ = R 2 + R 1 . The condition to obtain a connected minimal surface in terms of these parameters is In the following, we will show how to construct the vector field configuration whose flux across Σ AB computes the entanglement of purification for this system. Since the problem has reflection symmetry around x = 0, we will focus on the region x > 0 for simplicity. The first step is to consider a family of nested intervals that interpolate between the minimal surface m(B) and the vertical surface Σ AB . There are infinitely many ways to pick a family, since the two endpoints have variable locations. In order to choose one, we consider an auxiliary minimal surface m(Ã) for an intervalÃ centered at the origin.

JHEP05(2019)075
We impose that m(Ã) intersects normally the surface m(B) (and by symmetry, the surface m(A)), which implies that m(Ã) is a semi-circle with radius (5.20) Let us now takeṼ to be a geodesic flow (associated to regionÃ) of the kind constructed in section 2. The full set of integral curves ofṼ between m(B) and Σ AB gives us a natural family of nested intervals, which are specified by their endpoints. Therefore, using the nesting property one can construct a maximally packed flow which simultaneously maximizes the flux through B and all the intervals of the family. The associated integral curves of this maximally packed flow are simply given by the level set curves, with |Ṽ | = constant. We illustrate this explicitly in figure 9. More specifically, the level set curves are naturally parametrized by and angular parameter χ in terms of which |Ṽ | = sin χ, where This is a good parametrization for region inside m(Ã), namely, for x 2 + z 2 <R 2 . Alternatively, one can chose to parametrize the level set curves with a parameter ψ related to χ by in terms of which the |Ṽ | = constant surfaces adopt the simple form This is, they correspond to circles centered at (r c , z c ) = (0, ψR) and radius R c = 1 + ψ 2R . If one focus only on those integral curves lying entirely within the entanglement wedge associated to AB, r(AB), then, their flux will compute the entanglement of purification E(A : B). Two final comments are in order: (i) Notice that all integral lines of this maximally packed flow end at a single point of region A and B, respectively. This is perfectly valid. However, we can alternatively delete the portion of the threads inside m(A) and m(B) and put there instead a part of a geodesic flow. This is in complete analogy to the construction done in the BTZ case. And (ii) If desired, one can add an extra thread bundles that connect AB and AB and cross the regions of m(A) and m(B) that are not covered by the maximally packed flow. Such construction would give us a vector field that simultaneously computes S(A) (or S(B)) and the entanglement of purification E(A : B).
Adjacent intervals: when the two intervals share a boundary, the family of nested intervals required to construct our maximally packed flow is fully determined, since one of the endpoints is fixed. Therefore, one can easily obtain the corresponding flow by following the BTZ construction. However, since we already have a general construction for the disconnected case, one should be able to recover the case of adjacent intervals by taking an specific limit of the above. Indeed, by considering the limit in whichR → 0 while keeping ψR finite and arbitrary, it is possible to show that we can recover the correct flow for adjacent intervals. In this limit equation (5.23) becomes: The integral lines in this case correspond to circles centered at (r c , z c ) = (0, ψR) and radius R c = ψR. We illustrate this limit in figure 10. We notice again that all integral lines of this maximally packed flow end at a single point, although, in this case at a coincident point (x, z) = (0, 0) that lies in the boundary of A and B. This seems problematic at first glance, but all issues are alleviated by the presence of the UV cutoff . Alternatively one can delete the portion of the threads inside m(A) and m(B) and continue them with part of a geodesic flow, as discussed in the previous example.

Monogamy of mutual information
Very recently, the authors of [28] introduced the notion of multicomodity flows or multiflows: a set of flows that can coexist simultaneously in a given Riemannian geometry while satisfying some non-trivial properties that define them. Furthermore, they showed a set of theorems about multiflows that were used to prove the monogamy property of holographic mutual information. In this section we will give a quick review of the multiflow proposal and briefly summarize the connection with the monogamy property. We refer the reader

JHEP05(2019)075
A B Figure 10. A graphic illustration of the construction of a maximally packed flow for two adjacent intervals in AdS. The flow is obtained as a specific limit of the disjoint case, presented in figure 9. All conventions of figure 9 apply to the present one as well.
to [28] for a complete presentation, a detailed proof of the theorems as well as some conjectures about the entanglement structure of holographic states suggested by the bit thread picture of holographic entanglement entropy. Later in the same section, we will also show that our previous constructions can be used as individual components of such multiflows. As an application, we will show how the geodesic and maximally packed flows presented in the previous sections can be used to illustrate the monogamy property of holographic mutual information in concrete settings, complementing the alternative constructive proof presented in [26]. We do so by explicit construction in the case of two disjoint intervals in AdS 3 . However, there is in principle no limitation in considering more regions, other geometries or strips in higher dimensions.

Quick review of the multiflow proposal
Consider a Riemannian manifold M with boundary ∂M and take a set of non-overlapping boundary regions {A 1 , · · · A n } such that they cover the full boundary ∂M (∪ i A i = ∂M). A multiflow is a set of vector fields V ij on M that satisfy: Each V ij is by itself a flow which has non-zero flux only in the regions A i and A j and therefore, from (5.25), it follows that the flux through each of these two regions is equal in magnitude but opposite in sign, Furthermore, the condition (5.26) implies that the set of flows V i defined by

JHEP05(2019)075
are also good flows and therefore satisfy the bound A non-trivial property of these multiflows, proved in [28], is the existence of what is called a max multiflow.
Max multiflow: is a multiflow {V ij } such that for each i, the flow {V i } given in (5.28) is a max flow for A i , i.e., Based only on the existence of this max multiflow, the MMI follows almost trivially. In order to simplify the presentation we will use a slight generalization of the bit thread picture introduced in [28]. According to them, the threads are non-oriented one dimensional objects with a density ρ bounded by the Planck scale as 15 ρ ≡ length of threads inside a small ball volume of the ball The individual threads are allowed to intersect each other as long as the density bound is not violated. In turn, the correspondence with a smooth vector field is hence not one to one. Denoting N A:Ā as the number of threads connecting A and its complementĀ, then This is, the entanglement entropy of A is given by the maximum N A:Ā subject to the density constrain (5.31). The max multiflow (or multi-thread) theorem states that there exists a multiflow configuration N A i :A j obeying the density bound (5.31) such that 33) where N A i :A j is the number of threads connecting the region A i with A j , not necessarily the maximal one. which is known as the monogamy of mutual information.

Max multiflow and MMI for disjoint intervals
In this section we will present an explicit realization of the max multiflow for two disjoint intervals in pure AdS starting from the constructions of section 5.1.2. Let us consider two disjoint intervals in AdS 3 . Let A be the interval on the left and B the one on the right. Between A and B there is an interval, which we denote by C. Finally, we call D to the complement of ABC, namely, D = ABC. We are interested in the minimal surface and entanglement wedge associated to the disjoint region AB. As discussed earlier in section 5.1, depending on the relative separation of the intervals there are two qualitatively different scenarios or phases, which yield a connected or disconnected entanglement wedge. These two scenarios are related by purity: if m(AB) is connected then m(CD) is disconnected and vice-versa. This means that the thread configuration that describes the max flow of m(AB) in the connected case would correspond to the thread configuration that describes the max flow of m(CD) in their disconnected configuration. For concreteness, then, we will concentrate on the connected phase only. Additionally, we will consider the simplest case of two intervals of equal sizes since the most generic case with different sizes can be obtained from the former by a specific conformal transformation.
Consider the set of intervals X = {A, AB, ABC}. Given the nesting property of entanglement entropy, there should be a flow that simultaneously maximizes the flux through A, AB and ABC. Such a flow is fixed (up to a sign) at the location of the minimal surfaces m(A), m(AB) and m(ABC), this is V | m(X i ) = ±n| m(X i ) wheren| m(X i ) is the outward-pointing unit vector normal to m(X i ). The existence of such minimal surfaces m(X i ) suggests a natural separation of the bulk for a given boundary region X i into the co-dimension one bulk region 16 r(X i ) bounded by X i ∪m(X) and its complement M\r(X i ). The causal domain of dependence of r(X i ) is the so called entanglement wedge and plays a crucial role in the subregion duality of the AdS/CFT correspondence. However, r(X i ) is also informally called entanglement wedge so we do so throughout this paper. The idea now is to construct a flow maximizing the flux through A, AB and ABC by considering the various patches of the bulk that are naturally separated, this is r(A), r(AB)\r(A), r(ABC)\r(AB) and M\r(ABC). We do so by imposing the appropriate boundary conditions at m(A), m(AB) and m(ABC) and constructing flows for the individual patches.
First, let us consider one single interval X i , as given in section 2.1.2 for d = 1. In this case, the bulk is divided into two parts: M = r(X i ) ∪ (M\r(X i )). We will call the 16 Or co-dimension zero, if we are talking about a constant-t slice of the bulk geometry.

JHEP05(2019)075
part of V inside r(X i ) as the interior flow associated to X i and the part of V inside M\r(X i ) as the exterior flow. Now, consider the system of two disjoint regions A and B discussed above. Given the previous arguments, we can then construct flows in the regions inside r(A) and outside r(ABC) by taking the inside and outside parts of a geodesic flow associated to a single interval, respectively. Next, since we are considering the case with connected r(AB), i.e., m(AB) = m(ABC) ∪ m(C), it is easy to see that one must pick V | m(C) = −n| m(C) . This implies that the flow inside of r(C) can be chosen to be the negative of a geodesic flow of an interval of size C. 17 In summary, using the geodesic flow constructions of section 2.1.2 we have found the part of the flow in the bulk regions r(A), r(ABC)\r(AB) = r(C) and M\r(ABC). The only missing piece now is to find the part of the flow inside region r(AB)\r(A), which we proceed to do next.
Before describing how to construct the flow inside region r(AB)\r(A) we can use simple symmetry arguments put some constraints. First, notice that since the two intervals are of equal size then there is a reflection symmetry around the origin. Strictly speaking, this symmetry is broken by our choice of maximizing the flux through A instead of B. However, if one ignores the orientation of the threads then the flow can indeed be taken to respect this symmetry. A number of consequences follow from it. For example, since the number of threads connecting the region A andĀ is maximal the above symmetry implies that the number of threads connecting B andB is also maximal. This in turn implies that flow inside m(B) can also be taken to be the interior flow of a single interval. Another consequence is that the threads connecting D and A, and D and B, through m(ABC) are divided precisely by the symmetry preserving surface. The same also applies for the threads connecting C and A, and C and B, through m(C). Indeed, the part of the symmetry preserving surface that lies inside r(AB) is precisely equal to the crosssection whose area computes the entanglement of purification E(A : B). We therefore see that the entanglement of purification plays an interesting role in separating the multiflow components of a max multiflow. It would be very interesting to explore the extent to which this statement holds in more general situations, including general multiple regions, higher dimensions, and arbitrary geometries.
Finally, let us explicitly construct the advertised flow inside region r(AB)/r(A), focussing on the x < 0 part of the geometry. Recall that in section 5.1.2 we constructed maximally packed flows connecting minimal surfaces of intervals that either share a common boundary or are separated by some distance -see equations (5.23) and (5.24). In fact, these are precisely the kind of flows that we need to connect the different minimal surfaces in r(AB)/r(A). First, let us start from the threads connecting m(A) to m(ABC). These two surfaces share the left boundary, so it is easy to find a family of nested intervals whose minimal surfaces interpolate between the two. The threads that correspond to such maximally packed flow are called N A:D . Similarly, m(A) and m(C) share a common boundary so we can repeat the same process to find a maximally packed flow that connects them. We refer to these threads as and the symmetry preserving surface at x = 0 (or, equivalently, threads connecting m(A) and m(B)). These surfaces are separated by some distance, however, using the method explained in 5.1.2 we can propose a family of nested intervals whose minimal surfaces interpolate between the two, and hence construct the corresponding maximally packed flow. We use this construction for the portion of m(A) that has not been covered so far. We call these new threads N A:B . The parts of r(AB)/r(A) that are not covered by one of these thread bundles are assumed to have a vanishing flow. In figure 11 we represent graphically this construction. The flow constructed above is an explicit realization of the max multiflow theorem introduced in [28], for two disjoint intervals A and B. As discussed in section 5.2.1, the monogamy property of mutual information (5.39) follows immediately from the mere existence of such a flow. Indeed, one can illustrate the inequality by representing graphically the mutual information of the relevant pairs of intervals, as is done in figure 12. As it is shown in this figure, the flow that computes the mutual information between A and BC, I(A : BC), contains mutually disjoint flows that compute I(A : B) and I(A : C) separately. Therefore I(A : BC) is manifestly larger or equal to I(A : B) + I(A : C). This illustration has the flavor of the separation of flows that came from the commodity property in a slightly more detailed fashion.

Summary and discussion
In this paper, we have given explicit examples of vector field configurations that compute entanglement entropy in a variety of physical scenarios. As advertised in the introduction, our constructions can be categorized in three types: (i) symmetric or geodesic flows based on specific foliations of the bulk geometry, (ii) maximally packed flows, based on the nesting JHEP05(2019)075 Figure 12. In this figure we plot the minimal surfaces of all the individual regions as well as the different pairwise combinations. In particular we are interested in AB, AC, and ABC. We represent the different flows whose flux compute the mutual information between A and B, A and C, and A and BC in blue, red and yellow, respectively.
property of entanglement entropy, and (iii) mixed flows that are consistently composed of the two kinds of flows above.
The symmetric flows were used to explicitly construct flows for strips and spheres in empty AdS, in an arbitrary number of dimensions, and an interval in a BTZ background. In addition, we showed that a similar construction can be implemented in a more general background, provided that the metric satisfies certain curvature condition. For systems that can be dimensionally reduced to a 3-dimensional gravitational setup, e.g., strips in a general translationally invariant background, we showed that the condition simplifies drastically and reduces to having negative spatial curvature in the effective lower dimensional metric. Intuitively, this result can be explained by the fact that minimal length curves in a space of negative curvature tend to deviate from each other, as is the case for hyperbolic space. This is in contrast to the more familiar case of positive curvature (e.g. a sphere), where the curves tend to focus. For backgrounds supported by matter fields, we were also able to translate this statement in terms of the local energy density, for which we found a sharp upper bound.
We also presented two concrete applications of mixed flows to particular situations of interest: (i) a strip in a BTZ background and (ii) two disjoint intervals in empty AdS. The first case was already considered using geodesic flows; however, we showed that combining geodesic and maximally packed flows, it is possible to construct a configuration that simultaneously computes entanglement entropy and entanglement of purification of the system. Since the specific flow is found by maximizing the number of threads in the smallest cross-section of the entanglement wedge, we showed that the latter quantity can be naturally interpreted as the maximum number of Bell pairs that can be distilled from the mixed state. Finally, the case of two disjoint intervals was used to illustrate the monogamy property of mutual information. In order to do so, we provided an explicit example of the max multiflow theorem recently proposed in [28], and we showed that it can be graphically illustrated in an elegant way.

JHEP05(2019)075
There are some open questions related to our work that are worth exploring: • Reduced symmetry. It would be interesting to relax some of the symmetries that we have assumed in our constructions. For example, what can be said about situations in which the minimal surface is not known? This applies in particular to higher dimensional regions that cannot be dimensionally reduced. Are there special foliations that are useful in these cases? Perhaps by restricting the attention to families of lower dimensional cross-sections?
• Finite volume. All of our explicit constructions were carried out in the Poincaré patch of AdS or in a black brane geometry, hence, the dual field theory is defined on R 1,d+1 . It would be interesting to study situations where the field theory lives on a compact space, such as R × S d+1 . This would amount to study the problem in global AdS or in a black hole geometry. Both constructions, the one based on geodesics and the one based on the nesting property, are still valid in these scenarios. However, finite volume introduces new physical features such as the so-called entanglement plateaux [29], entanglement shadows [30] and the existence of winding/entwinement geodesics [31]. It would be interesting to study the interplay between all these new phenomena, their role in the construction of bit threads configurations and their interpretation in this language.
• Covariant generalization. In this paper we have focussed on static situations, where the entanglement entropy is computed by the standard RT prescription [1,2]. In dynamical situations, however, one needs to upgrade the recipe and use the covariant HRRT formula instead [3]. A covariant proposal of bit threads is currently under development [32]. It would be very interesting to explicitly construct examples of thread configurations in relevant dynamical situations, e.g., a global quantum quench. The bit thread picture in this scenario would yield important insights regarding the entanglement tsunami proposal of entanglement propagation [33,34] and its breakdown [35].
• Dynamics. in perturbative excited states, entanglement entropy satisfies the so-called first law of entanglement, in which the change in entanglement entropy is given by the change in modular energy. In the bulk, this statement is mapped into the linearized Einstein equations around empty AdS [6,7]. It would be very interesting to see how this translates into the bit thread picture. It would also be interested to go beyond linear level and investigate the imprints on quantities such as relative entropy and quantum Fisher information [36,37]. Some work in this direction is already underway [38].
• Multipartite entanglement. Our explicit construction for two disjoint intervals can be straightforwardly generalized to an arbitrary number of parties. It would be very interesting to explicitly do so, at least for 3 and 4 intervals, and find how the holographic inequalities that must be satisfied in each case [39,40] are reflected in terms of bit threads. In particular, it would be insightful to represent visually these JHEP05(2019)075 inequalities in terms of thread bundles connecting the various regions, as we did for the monogamy of mutual information. We expect the construction of flows in these scenarios to be intrinsically related to the recent generalizations of entanglement of purification for multipartite systems [41][42][43].
• Bulk reconstruction. The program of hole-ography [44][45][46][47] aims to reconstruct arbitrary bulk surfaces by cleverly adding and subtracting the entanglement entropies of a family of intervals associated to the surface. Interestingly, after taking a continuum limit and shrinking the surface to a point it is possible to recover an intrinsic definition of a bulk point in terms of a CFT quantity dubbed as differential entropy [48]. It would be interesting to explore similar questions in the language of bit threads and shed light on the issue of bulk reconstruction and bulk locality [49].
• Higher derivatives. In higher curvature gravity, the holographic entanglement entropy is computed by the Dong-Camps prescription [50,51] which includes Wald's formula for black hole entropy, as well as corrections involving the extrinsic curvature. The bit thread proposal was recently generalized to arbitrary higher curvature gravity in [52].
In would be interesting to construct explicit examples of flows using the generalized prescription to gain insight into the stringy or α corrections to entanglement entropy.
• Quantum corrections. The leading quantum or 1/N corrections to holographic entanglement entropy are given by bulk entanglement entropy across the RT surface [53]. One immediate consequence is that threads can now connect points in the bulk that are not necessarily at the boundary of AdS, and therefore one must relax the divergenceless condition. It would be interesting to investigate how different bulk fields couple to V and how they source specific thread configurations. One specific example that might be of interest is the case of two disjoint intervals, in the case that the entanglement wedge is disconnected [54]. For this system, the leading 1/N corrections to entanglement entropy gives the leading term in the mutual information I(A, B).
• Tensor networks. Finally, we remark that the bit thread picture of holographic entanglement entropy, as well as the proposal to compute entanglement of purification, were motivated in part by notions of cMERA and other holographic models tensor networks [55][56][57][58]. An interesting avenue for further research would be to consider the explicit thread configurations constructed in this paper and extrapolate them back to the realm of tensor networks. For example, it would be natural to associate the magnitude |V | to a density of tensors in a continuous network on AdS. The associated integral lines might in turn represent paths of optimal displacement along the RG scale. If so, do the geodesic and maximally packed flows constructed here have a specific interpretation?
We hope to come back to some of these points in the near future.