On shape dependence of holographic mutual information in AdS4

We study the holographic mutual information in AdS4 of disjoint spatial domains in the boundary which are delimited by smooth closed curves. A numerical method which approximates a local minimum of the area functional through triangulated surfaces is employed. After some checks of the method against existing analytic results for the holographic entanglement entropy, we compute the holographic mutual information of equal domains delimited by ellipses, superellipses or the boundaries of two dimensional spherocylinders, finding also the corresponding transition curves along which the holographic mutual information vanishes.


Introduction
Entanglement entropy has been extensively studied during the last decade and its important role in quantum gravity, quantum field theory and condensed matter physics is widely recognized. Given a quantum system in its ground state |ψ and assuming that its Hilbert space can be decomposed as H = H A ⊗ H B , one can introduce the reduced density matrix ρ A ≡ Tr B ρ by tracing over H B the density matrix ρ = |ψ ψ| of the whole system. Here we focus on a bipartition of the Hilbert space associated with a separation of a spatial slice into two complementary regions. The entanglement entropy is the Von Neumann entropy associated with ρ A , namely S A ≡ −Tr A (ρ A log ρ A ), and it measures the entanglement between A and B. In the same way, one can introduce ρ B ≡ Tr A ρ and S B . Since ρ is a pure state, we have that S B = S A . Understanding the dependence of S A on the geometry of the region A is an important task.
Let us consider a conformal field theory in D + 1 dimensions at zero temperature in its ground state. The entanglement entropy S A between a D dimensional spatial region A and its complement B can be written as an expansion in the ultraviolet cutoff ε, where the leading divergence is S A ∝ Area(∂A)/ε D−1 + . . . [1,2]. This behaviour is known as the area law for the entanglement entropy and ∂A is sometimes called entangling surface. When D = 1 and the domain A is an interval, ∂A is made by its two endpoints and the area law is violated because the leading divergence is logarithmic. In particular, S A = (c/3) log( /ε) + const, where c is the central charge of the model [3,4].
By virtue of the holographic correspondence [5][6][7] (see [8] for a review), the entanglement entropy S A of a conformal field theory in a D + 1 dimensional Minkowski spacetime can be also calculated from its dual gravitational model defined in a D + 2 dimensional asymptotically anti-De Sitter (AdS) spacetime whose boundary is the spacetime of the original conformal field theory. In the regime where it is enough to consider only classical gravity, the holographic prescription to compute the entanglement entropy is [9,10] (1.1) where G N is the D + 2 dimensional Newton constant and A A is the area of the codimension two minimal area spacelike surfaceγ A at some fixed time slice such that ∂γ A = ∂A. Sinceγ A reaches the boundary of the asymptotically AdS D+2 spacetime, its area A A is divergent and therefore it must be regularized through the introduction of a cutoff ε in the holographic direction, which corresponds to the ultraviolet cutoff of the dual conformal field theory. The leading divergence of (1.1) as ε → 0 provides the area law of the entanglement entropy. The covariant generalization of (1.1) has been proposed in [11] and it has been extensively employed to study holographic models of thermalization. Recent reviews on entanglement entropy in quantum field theory and holography are [12][13][14].
The minimal area surfaces anchored on a given curve defined on the boundary of AdS D+2 occur also in the holographic dual of the expectation values of the Wilson loops [15,16]. Nevertheless, while the bulk surfaces for the Wilson loops are always two dimensional, for the holographic entanglement entropy they have codimension two. Thus, when D = 2 the minimal surfaces to compute for the holographic entanglement entropy (1.1) are the same ones occurring in the gravitational counterpart of the correlators of spacelike Wilson loops.
As for the dependence of A A on the geometry of ∂A, analytic results have been found for the infinite strip and for the sphere when D is generic [9,10]. Spherical domains play a particular role because their reduced density matrix can be related to a thermal one [17]. When D = 2, the O(1) term in the expansion of S A as ε → 0 for circular domains provides the quantity F , which decreases along any renormalization group flow [18][19][20]. Some interesting results have been found about A A for an entangling surface ∂A with a generic shape [21][22][23][24][25][26][27][28], but a complete understanding is still lacking.
When A = A 1 ∪ A 2 is made by two disjoint spatial regions, an important quantity to study is the mutual information It is worth remarking that S A1∪A2 provides the entanglement between A 1 ∪ A 2 and the remaining part of the spatial slice. In particular, it does not quantify the entanglement between A 1 and A 2 , which is measured by other quantities, such as the logarithmic negativity [29][30][31][32]. In the combination (1.2), the area law divergent terms cancel and the subadditivity of the entanglement entropy guarantees that I A1,A2 0. For two dimensional conformal field theories, the mutual information depends on the full operator content of the model [33][34][35][36]. When D 2, the computation of (1.2) is more difficult because non local operators must be introduced along ∂A [37][38][39].
The holographic mutual information is (1.2) with S A given by (1.1). The crucial term to evaluate is S A1∪A2 , which depends on the geometric features of the entangling surface ∂A = ∂A 1 ∪ ∂A 2 , including also the distance between A 1 and A 2 and their relative orientation, being ∂A made by two disjoint components. It is well known that, keeping the geometry of A 1 and A 2 fixed while their distance increases, the holographic mutual information has a kind of phase transition with discontinuous first derivative, such that I A1,A2 = 0 when the two regions are distant enough. This is due to the competition between two minima corresponding to a connected configuration and to a disconnected one. While the former is minimal at small distances, the latter is favoured for large distances, where the holographic mutual information therefore vanishes [40][41][42]. This phenomenon has been also studied much earlier in the context of the gravitational counterpart of the expectation values of circular spacelike Wilson loops [43][44][45][46]. The transition of the holographic mutual information is a peculiar prediction of (1.1) and it does not occur if the quantum corrections are taken into account [47]. A similar transition due to the competition of two local minima of the area functional occurs also for the holographic entanglement entropy of a single region at finite temperature [48][49][50].
In this paper we focus on D = 2 and we study the shape dependence of the holographic entanglement entropy and of the holographic mutual information (1.1) in AdS 4 , which is dual to the zero temperature vacuum state of the three dimensional conformal field theory on the boundary. This reduces to finding the minimal area surfaceγ A spanning a given boundary curve ∂A (the entangling curve) defined in some spatial slice of the boundary of AdS 4 . The entangling curve ∂A could be made by many disconnected components. When ∂A consists of one or two circles, the problem is analytically tractable [9,10,15,[51][52][53][54][55]. However, for an entangling curve having a generic shape (and possibly many components), finding analytic solutions becomes a formidable task. In order to make some progress, we tackle the problem numerically with the help of Surface Evolver [56,57], a widely used open source software for the modelling of liquid surfaces shaped by various forces and constraints. A section at constant time of AdS 4 gives the Euclidean hyperbolic space H 3 . Once the curve embedded in H 3 is chosen, this software constructs a triangular mesh which approximates the surface spanning such curve which is a local minimum of the area functional, computing also the corresponding finite area. The number of vertices V , edges E and faces E of the mesh are related via the Euler formula, namely V − E + F = χ, being χ = 2 − 2g − b the Euler characteristic of the surface, where g is its genus and b the number of its boundaries. In this paper we deal with surfaces of genus g = 0 with one or more boundaries.
The paper is organized as follows. In §2 we state the problem, introduce the basic notation and review some properties of the minimal surfaces occurring in our computations. In §3 we address the case of surfaces spanning simply connected curves. First we review two analytically tractable examples, the circle and the infinite strip; then we address the case of some elongated curves (i.e. ellipse, superellipse and the boundary of the two dimensional spherocylinder) and polygons. Star shaped and non convex domains are also briefly discussed. In §4 we consider ∂A made by two disjoint curves. The minimal surface spanning such disconnected curve can be either connected or disconnected, depending on the geometrical features of the boundary, including the distance between them and their relative orientation. The cases of surfaces spanning two disjoint circles, ellipses, superellipses and the boundaries of two dimensional spherocylinders are quantitatively investigated for a particular relative orientation. Further discussions and technical details are reported in the appendices.

Minimal surfaces in AdS 4
Finding the minimal area surface spanning a curve is a classic problem in geometry and physics. In R 3 this is known as Plateau's problem. A physical realization of the problem is obtained by dipping a stiff wire frame of some given shape in soapy water and then removing it: as the energy of the film is proportional to the area of the water/air interface, the lowest energy configuration consists of a surface of minimal area. In this mundane setting, the requirement of minimal area results into a well known equation where H = k i i /2 is the mean-curvature given by the trace of the extrinsic curvature tensor k ij = e i,j ·N , with N the surface normal vector, e i a generic tangent vector, such that the surface metric tensor is h ij = e i · e j , and ( · ) ,i = ∂ i ( · ).
The metric of AdS 4 in Poincaré coordinates reads where the AdS radius has been set to one for simplicity. The spatial slice t = const provides the Euclidean hyperbolic space H 3 and the region A is defined in the z = 0 plane. According to the prescription of [9,10], to compute the holographic entanglement entropy, first we have to restrict ourselves to a t = const slice and then we have to find, among all the surfaces γ A spanning the curve ∂A, the one minimizing the area functional where U A is a coordinate patch associated with the coordinates (u 1 , u 2 ) and h = det(h ij ). We denote byγ A the area minimizing surface, so that A[γ A ] ≡ A A provides the holographic entanglement entropy through the Ryu-Takayanagi formula (1.1). Since all the surfaces γ A reach the boundary of AdS 4 , their area is divergent and therefore one needs to introduce a cutoff in the holographic direction to regularize it, namely z ε > 0, where ε is an infinitesimal parameter. The holographic dictionary tells us that this cutoff corresponds to the ultraviolet cutoff in the dual three dimensional conformal field theory. Considering z ε > 0, the area A[γ A ] and therefore A A as well become ε dependent quantities which diverge when ε → 0. Important insights can be found by writing A A as an expansion for ε → 0. When ∂A is a smooth curve, this expansion reads where P A = length(∂A) is the perimeter of the entangling curve and o(1) indicates vanishing terms when ε → 0. When the entangling curve curve ∂A contains a finite number of vertices, also a logarithmic divergence occurs, namely The functions F A , B A and W A are defined through (2.4) and (2.5). They depend on the geometry of ∂A in a very non trivial way. We remark that the section ofγ A at z = ε provides a curve which does not coincide with ∂A because of the non trivial profile ofγ A in the bulk. As the area element in AdS 4 is factorized in the form dA = du 1 du 2 √ h/z 2 , a surface in AdS 4 is equivalent to a surface in R 3 endowed with a potential energy density of the form 1/z 2 . By using the standard machinery of surface geometry (see §A), one can find an analog of (2.1) in the form whereẑ is a unit vector in the z direction. The relation (2.6) implies that, in order for the mean curvature to be finite, the surface must be orthogonal to the (x, y) plane at z = 0: i.e.ẑ · N = 0 at z = 0. As a consequence of the latter property, the boundary is also a geodesic ofγ A (see §A).

Simply connected regions
In this section we consider cases in which the region A is a simply connected domain. We first review the simple examples of the disk and of the infinite strip, which can be solved analytically [9,10]. In §3.1 we numerically analyze the case in which A is an elongated region delimited by either an ellipse, a superellipse or the boundary of a two dimensional spherocylinder, while in §3.2 we address the case in which ∂A is a regular polygon. In §3.3, star shaped and non convex domains are briefly discussed. If A is a disk of radius R, the minimal area surfaceγ A is a hemisphere, as it can be easily proved from a direct substitution in (2.6). Taking N = r/|r|, with r = (x, y, z) and |r| = R, one findsẑ · N = z/R, hence H = −1/R, which is the mean curvature of a sphere whose normal is outward directed. The area of the part of the hemisphere such that ε z R is Comparing this expression with (2.4), one finds that F A = 2π in this case. It is worth remarking , as peculiar feature of the disk, that in (3.1) o(1) terms do not occur. A special case of (2.6) is obtained when the surface is fully described by a function z = z(x, y) representing the height of the surface above the (x, y) plane at z = 0. In this case and (2.6) becomes the following second order non linear partial differential equation for z (see §A for some details on this derivation) with the boundary condition that z = 0 when (x, y) ∈ ∂A. The partial differential equation (3.3) is very difficult to solve analytically for a generic curve ∂A; but for some domains A it reduces to an ordinary differential equation. Apart from the simple hemispherical case previously discussed, this happens also for an infinite strip A = {(x, y) ∈ R 2 , |y| R 2 }, whose width is 2R 2 . The corresponding minimal surface is invariant along the x axis and therefore it is fully characterized by the profile z = z(y) for |y| R 2 . Taking Equivalently, the infinite strip case can be studied by considering the one dimensional problem obtained substituting z = z(y) directly in (3.2) [8][9][10]. Since the resulting effective Lagrangian does not depend on y explicitly, one easily finds that z 2 1 + z 2 ,y is independent of y. Taking the derivative with respect to y of this conservation law, (3.4) is recovered, as expected. The constant value can be found by considering y = 0, where z(0) ≡ z * and z ,y (0) = 0. Notice that z * is the maximal height attained by the curve along the z direction. Integrating the conservation law, one gets where Γ is the Euler gamma and 2 F 1 is the hypergeometric function. Thus, the minimal surfaceγ A consists of a tunnel of infinite length along the x direction, finite width R 2 along the y direction and whose shape in the (y, z) plane is described by (3.5). Considering a finite piece of this surface which extends for R 1 R 2 in the x direction, whose projection on the (x, y) plane is delimited by the dashed lines in the bottom panel of Fig. 1, its area is given by [9,10,15,16] where ε z z * . Comparing (2.4) with P A = 4R 1 and (3.6), one concludes that F A = s ∞ R 1 /R 2 . In order to compare (3.6) with our numerical results, we find it useful to construct an auxiliary surface by closing this long tunnel segment with two planar "caps" placed at x = ±R 1 , whose profile is described in the (y, z) plane by (3.5), with a cutoff at z = ε. These regions are identical by construction and their area (see §D.1) is given by A cap = 2R 2 /ε − π/2 + o (1). Thus, the total area of the auxiliary surface reads where the coefficient of the leading divergence is the perimeter of the rectangle in the boundary (dashed curve in Fig. 1). It is worth remarking that this surface is not the minimal area surface anchored on the dashed rectangle in Fig. 1. Indeed, in this case an additional logarithmic divergence occurs (see §3.2). Since in the following we will compute numerically A A for various domains keeping ε fixed, let us introduce In Fig. 2 the values of F A for the surfaces discussed above are represented together with other ones coming from different curves that will be introduced in §3.1: the black dot corresponds to the disk (see (3.1)), the dotted horizontal line is obtained from (3.6) for the infinite strip, while the dashed line is found from the area (3.7) of the auxiliary surface.
x y z x y R 1 R 2 Figure 1: Top panel: Minimal surfaces constructed by using Surface Evolver where the entangling curve ∂A is a circle with radius R = 1 (red), an ellipse (orange), a superellipse (3.9) with n = 8 (purple) and the boundary of a spherocylinder (green) with R 1 = 3R 2 . The cutoff is ε = 0.03 and only the y 0 part of the minimal surfaces has been depicted to highlight the curves provided by the section y = 0. Bottom panel: In the (x, y) plane, we show the superellipses with R 1 = 3R 2 with n = 2 (orange), n = 4 (blue), n = 6 (magenta) and n = 8 (purple), the circle with radius R 1 (red curve) and the rectangle circumscribing the superellipses (dashed lines). The green curve is the boundary of the two dimensional spherocylinder with R 2 = 3R 1 .

Superellipse and two dimensional spherocylinder
The first examples of entangling curves ∂A we consider for which analytic expressions of the corresponding minimal surfaces are not known are the superellipse and the boundary of the two dimensional spherocylinder, whose geometries depend on two parameters. The two dimensional spherocylinder nicely interpolates between the circle and the infinite strip.
In Cartesian coordinates, a superellipse centered in the origin with axes parallel to the coordinate axes  Figure 2: Numerical data for F A , defined in (3.8), corresponding to domains A which are two dimensional spherocylinders or delimited by superellipses. Here ε = 0.03. In the main plot R 2 = 1, while in the inset, which shows a zoom of the initial part of the main plot in logarithmic scale on both the axes, we have also reported data with R 2 = 2. The horizontal dotted black line corresponds to the infinite strip (3.6) and the dashed one to the auxiliary surface where the sections at x = ±R 1 have been added (see (3.7)). The red and blue dotted horizontal lines come from the asymptotic result (C.10) evaluated for n = 2 and n = 3 respectively.
is described by the equation where R 1 , R 2 and n are real and positive parameters. The curve (3.9) is also known as Lamé curve and here we consider only integers n 2 for simplicity. The special case n = 2 in (3.9) is the ellipse with semi-major and semi-minor axes given by R 1 and R 2 respectively. As the positive integer n increases, the superellipse approximates the rectangle with sides 2R 1 and 2R 2 . When R 1 = R 2 , the curves (3.9) for various n are known as squircles because they have intermediate properties between the ones of a circle (n = 2) and the ones of a square (n → ∞). In the bottom panel of Fig. 1, we show some superellipses with R 1 = 3R 2 , the circle with radius R 1 included in all the superellipses and the rectangle circumscribing them. In order to study the interpolation between the circle and the infinite strip, a useful domain to consider is the two dimensional spherocylinder. The spherocylinder (also called capsule) is a three dimensional volume consisting of a cylinder with hemispherical ends. Here we are interested in its two dimensional version, which is a rectangle with semicircular caps. In particular, the two dimensional spherocylinder circumscribed by the rectangle with sides 2R 1 and 2R 2 is defined as the set S ≡ D ∪ C + ∪ C − , where the rectangle D and the disks C ± are The perimeter of this domain is P A = 2πR 2 + 4(R 1 − R 2 ) and an explicit example of ∂S with R 2 = 3R 1 is given by the green curve in the bottom panel of Fig. 1. When R 1 = R 2 , the curve ∂S becomes a circle, while for R 1 R 2 it provides a kind of regularization of the infinite strip. Indeed, when R 1 → ∞ at fixed R 2 the two dimensional spherocylinder S becomes the infinite strip with width 2R 2 . Let us remark that the curvature of ∂S is discontinuous while the curvature of the superellipse (3.9) is continuous. Moreover, the choice to regularize the infinite strip through the circles C ± in (3.10) is arbitrary; other domains can be chosen (e.g. regions bounded by superellipses) without introducing vertices in the entangling curve. A straightforward numerical analysis allows to observe that a superellipses with n > 2 intersects once the curve ∂S in the first quadrant outside the Cartesian axes.
In Fig. 2 we show the numerical data for F A , defined in (3.8), when A is given by the domains discussed above: disk, infinite strip, two dimensional spherocylinder and two dimensional regions delimited by superellipses. In particular, referring to the bottom panel of Fig. 1, we fixed R 2 and increased R 1 . For the two dimensional spherocylinder, this provides an interpolation between the circle and the infinite strip. Surface Evolver has been employed to compute the area A A and for the cutoff in the holographic direction we choose ε = 0.03. Below this value, the convergence of the local minimization algorithm employed by Surface Evolver becomes problematic, as well as for too large domains A, as discussed in §B.
When R 1 = R 2 , we observe that F A for the squircles with different n > 2 increases with n. For large R 1 /R 2 , the limits of F A /(R 1 /R 2 ) for the domains we address are finite and positive. The values of these limits associated with the superellipses are ordered in the opposite way in n with respect to the starting point at R 1 = R 2 and therefore they cross each other as R 1 /R 2 increases. We remark that the curve corresponding to the two dimensional spherocylinder stays below the ones associated with the superellipses for the whole range of R 1 /R 2 that we considered. In Fig. 2 the horizontal black dotted line corresponds to the infinite strip (see (3.6)) while the dashed curve is obtained from the auxiliary surface described above (see (3.7)). The latter one is our best analytic approximation of the data corresponding to the two dimensional spherocylinder.
Focussing on the regime of large R 1 /R 2 , from Fig. 2 we observe that the asymptotic value of F A /(R 1 /R 2 ) for the two dimensional spherocylinder is very close to the one of the auxiliary surface obtained from (3.7) and therefore it is our best approximation of the result corresponding to the infinite strip. This is reasonable because the two dimensional spherocylinder is a way to regularize the infinite strip without introducing vertices in the entangling curve, as already remarked above. As for the minimal surfaces spanning a superellipse with a given n 2, in §C an asymptotic lower bound is obtained (see (C.10)), generalizing the construction of [28]. In Fig. 2 this bound is shown explicitly for n = 2 and n = 3 (red and blue dotted horizontal lines respectively). Since this value is strictly larger than the corresponding one associated with the infinite strip (see (3.6)), we can conclude that F A /(R 1 /R 2 ) for the superellipse at fixed n does not converge to the value s ∞ associated with the infinite strip.

Polygons
In this section we consider the minimal area surfaces associated with simply connected regions A whose boundary is a convex polygon with N sides. These are prototypical examples of minimal surfaces spanning entangling curves with geometric singularities. For quantum field theory results about the entanglement entropy of domains delimited by such curves, see e.g. [58][59][60].
The main feature to observe about the area A A of the minimal surface is the occurrence of a logarithmic divergence, besides the leading one associated with the area law, in its expansion as ε → 0. We find it convenient to introduce Since (2.5) holds in this case, we have that When ∂A is a convex polygon with N sides, denoting by α i < π its internal angle at the i-th vertex, for  Figure 4: Left: Section of the minimal surfaces anchored to an equilateral triangle (red, magenta and purple points), a square (blue points) or an octagon (green points) inscribed in a circle, as indicated in the inset by the black line. The continuos lines are z = ρ/f 0 (α N ), where f 0 (α) is found from (3.15) with N = 3 (red), N = 4 (blue) or N = 8 (green). The dashed black curve is the hemisphere corresponding to the circle circumscribing the polygons at z = 0 (dashed in the inset), while the dashed grey horizontal line corresponds to the cutoff ε = 0.03. Right: A zoom of the left panel around the origin, placed in the common vertex of the polygons. For the triangle, three different values of ε ∈ {0.03, 0.02, 0.01} has been considered to highlight how the agreement with the analytic result improves as ε → 0. the coefficient of the logarithmic term in (2.5) we can write (3.12) The function b(α) has been first found in [61], where the holographic duals of the correlators of Wilson loops with cusps have been studied, by considering the minimal surface near a cusp whose opening angle is α. Notice that (3.12) does not depend on the lengths of the edges but only on the convex angles of the polygon. Further interesting results have been obtained in the context of the holographic entanglement entropy [52,62]. Introducing the polar coordinates (ρ, φ) in the z = 0 plane, one considers the domain {|φ| α/2 , ρ < L}, Figure 5: The quantity B A in (3.11) with A A evaluated with Surface Evolver when the entangling curve ∂A is either an isosceles triangle whose basis has length (top panel) or a rhombus whose side length is (bottom panel). Here ε = 0.03. The black continuous curves are obtained from (3.12) and (3.16).
where L 1. By employing scale invariance, one introduces the following ansatz [61] z = ρ f (φ) , (3.13) in terms of a positive function f (φ), which is even in the domain |φ| α/2 and f → +∞ for |φ| → α/2. Plugging (3.13) into the area functional, the problem becomes one dimensional, similarly to the case of the infinite strip slightly discussed in §3. Since the resulting integrand does not depend explicitly on φ, the corresponding conservation law tells us that (f 4 + f 2 )/ (f ) 2 + f 4 + f 2 is independent of φ. Thus, the profile for 0 φ < α/2 (the part of the surface with −α/2 < φ 0 is obtained by symmetry) is given by 14) being f 0 ≡ f (0). When f → ∞, we require that the l.h.s. of (3.14) becomes α/2 and, by inverting the resulting relation, one finds f 0 = f 0 (α). In this limit the integral in (3.14) can be evaluated analytically in terms of elliptic integrals Π and K (see §E for their definitions) as follows Notice that when f 0 → 0 we have α → π, which means absence of the corner, while α → 0 for f 0 → ∞.
An interesting family of curves to study is the one made by the convex regular polygons. They are equilateral, equiangular and all vertices lie on a circle. For instance, a rhombus does not belong to this family. Denoting by R the radius of the circumscribed circle and by N the number of sides, the length of each side is = 2R sin(π/N ) and all the internal angles are α N ≡ N −2 N π. When N → ∞ we have that α N → π and the polygon becomes a circle. Thus, the area of the minimal surface spanning these regular polygons is (2.5) with P A = N and It is interesting to compare the analytic results presented above with the corresponding numerical ones obtained with Surface Evolver. Some examples of minimal surfaces anchored on curves ∂A given by a polygon are given in Fig. 3, where the triangulations are explicitly shown. In Fig. 4 we take as ∂A an equilateral triangle, a square and an octagon which share a vertex and consider the section of the corresponding minimal surfaces through a vertical plane which bisects the angles associated with the common vertex, as shown in the inset of the left panel. Focussing on the part of the curves near the common vertex, we find that the numerical results are in good agreement with the analytic expression (3.14). It would be interesting to find analytic results for the profiles shown in the left panel of Fig. 4.
By employing Surface Evolver, we can also consider entangling curves given by polygons which are not regular, as done in Fig. 5, where we have reported the data for B A (defined in (3.11)) corresponding to the area of the minimal surfacesγ A when ∂A is either an isosceles triangle (top panel) or a rhombus with side (bottom panel) . These examples allow us to consider also cusps with small opening angles. The size of the isosceles triangles has been changed by varying the angles α adjacent to the basis. Thus, the limiting regimes are the segment (α = 0) and the semi infinite strip (α = π). As for the rhombus, denoting by α the angle indicated in the inset, its limiting regimes are the segment (α = 0) and the square (α = π). The cutoff in the holographic direction has been fixed to ε = 0.03 (see the discussion in §B). Increasing the size of the polygon improves the agreement with the curve given by (3.12) and (3.16), as expected, because ε/P A gets closer to zero. Moreover, the agreement between the numerical data and the analytic curve gets worse as α becomes very small.
In Fig. 6 we report the data for B A found with Surface Evolver for regular polygons with various number N of edges. The agreement with the curve given by (3.12) and (3.16) is quite good and it improves for larger domains.
It is worth emphasizing that, for entangling surfaces ∂A containing corners, the way we have employed to construct the minimal surfaces with Surface Evolver (i.e. by defining ∂A at z = ε) influences the term W A in the expansion (2.5) for the area, as already remarked in [61].
It could be helpful to compute the length P ε of the curve defined as the section at z = ε of the minimal surface anchored on the long segments of a large wedge with opening angle α, which has been introduced above. From (3.13) we find that, in terms of polar coordinates whose center is the projection of the vertex at z = ε, this curve is given by ρ = εf (φ). Being L 1, we find that P ε reads (3.17) where α ε α is defined by the relation L = εf (α ε /2) and in the last step a change of variable has been performed. It is easy to observe that α ε < α. Considering the integral in the intermediate step of (3.17), one notices that it diverges because of its upper limit of integration (see the text below (3.13)), while the lower limit of integration gives a finite result, providing a contribution O(ε) to P ε . The expression of ∂ f φ can be read from the integrand of (3.14), finding that where the finite term has been found numerically. As a cross check of the finite term, we observe that f 0 = 0 when α = π (see below (3.15)), as expected. Thus, we can conclude that P ε = 2L + O(ε), being P A = 2L the length of the boundary of the wedge at z = 0. Notice that, performing this computation for the minimal surface anchored on a circle of radius R, which is a hemisphere, one finds that P ε = 2πR + O(ε 2 ).
Let us remark that P ε is not related to the regularization we adopt in our numerical analysis, as it can be realized from the right panel of Fig. 4. Indeed, in order to analytically the profiles given by the numerical data in the right panel of Fig. 4 the ansatz (3.13) cannot be employed and a partial differential equation must be solved.

Star shaped and non convex regions
The crucial assumption throughout the above discussions is that the minimal surfaceγ A can be fully described by z = z(x, y), where (x, y) ∈ A. Nevertheless, there are many domains A for which this parameterization cannot be employed because there are pairs of different points belonging to the minimal surfacesγ A with the same projection (x, y) / ∈ A in the z = 0 plane. In these cases, being the analytic approach quite difficult in general, one can employ our numerical method to find the minimal surfaces and to compute their area. The numerical data obtained with Surface Evolver would be an important benchmark for analytic results that could be found in the future.
An interesting class of two dimensional regions to consider is given by the star shaped domains. A region A at z = 0 belongs to this set of domains if a point P 0 ∈ A exists such that the segment connecting any other point of the region to P 0 entirely belongs to A. As for the minimal surface anchored on a star shaped domain A, by introducing a spherical polar coordinates system (r, φ, θ) centered in P 0 (the angular ranges are φ ∈ [0, 2π) and θ ∈ [0, π/2]), one can parameterize the entire minimal surface. Thus, we have ρ = r sin θ and z = r cos θ, being (ρ, φ) the polar coordinates of the z = 0 plane. Some interesting analytic results about these domains have been already found. In particular, [22] considered minimal surfaces obtained as smooth perturbations around the hemisphere and in [23] the behaviour in the IR regime for gapped backgrounds [63] has been studied. Our numerical method allows a more complete analysis because, within our approximations, we can find (numerically) the area of the corresponding minimal surface without restrictions.
In Fig. 7 we show a star convex domain A delimited by the red curve at z = 0, which does not contain vertices, and the corresponding minimal surfaceγ A anchored on it. Notice that there are pairs of points belonging toγ A having the same projection (x, y) / ∈ A on the z = 0 plane. It is worth recalling that in our regularization the numerical construction of the minimal surface with Surface Evolver has been done by defining the entangling curve ∂A at z = ε.
In order to give a further check of our numerical method, we find it useful to compare our numerical results against the analytic ones obtained in [22], where the equation of motion coming from (2.3) written in polar coordinates (r, φ, θ) has been linearized to second order around the hemisphere solution with radius R, finding where the r 1 (θ, φ) and r 2 (θ, φ) are given by [22] being k ∈ N and µ ∈ R two parameters of the linearized solution. The minimal surface equation coming from (2.3) is satisfied by (3.18) at O(a 2 ). Notice that r 1 (θ = 0, φ) = r 2 (θ = 0, φ) = 0, which means that the maximum value reached by the linearized solution along the z direction is R, like for the hemisphere.
Neglecting the O(a 3 ) terms in (3.18), one has a surface spanning the curve r(π/2, φ) ≡ R 2 (φ) at z = 0, z Figure 9: Minimal surfaces constructed with Surface Evolver corresponding to non convex domains at z = 0 delimited by the red and blue curves, which are made by arcs of circle centered either in the origin or in the points identified by the black dots. The green and magenta curves are sections of the minimal surfaces anchored on the red and the blue curves respectively.
which reads In Fig. 8 we construct the minimal surfaces providing the holographic entanglement entropy of some examples of star shaped regions A delimited by (3.21) where R and µ are kept fixed while a takes different values, taking the φ = π/4 section of these surfaces (see also the green curve in Fig. 7). Compare the resulting curves (the solid ones in the main plot of Fig. 8) with the corresponding ones obtained from the second order linearized solution (3.18) (made by the empty circles), we observe that the agreement is very good for small values of a/R and it gets worse as a/R increases, as expected. Our numerical method is interesting because it does not rely on any particular parameterization of the surface and this allows us to study the most generic non convex domain. In Fig. 9 we show two examples of non convex domains A which are not star shaped: one is delimited by the red curve and the other one by the blue curve. We could see these domains as two two dimensional spherocylinders which have been bended in a particular way. Constructing the minimal surfacesγ A anchored on their boundaries and considering their sections given by the green and magenta curves, one can clearly observe that some pair of points belonging to the minimal surfaces have the same projection (x, y) / ∈ A on the z = 0 plane, as already remarked above. An analytic description of these surfaces is more difficult with respect to the minimal surfaces anchored on the boundary of star shaped domains because it would require more patches. z Figure 10: Minimal surface constructed with Surface Evolver for a domain A = A 1 ∪ A 2 delimited by two disjoint and equal ellipses at z = 0 (blue curves). Here ε = 0.03 and the minimal surface is anchored on ∂A defined at z = ε, according to our regularization prescription. The minimal surface has (V, F ) = (18936, 37616) (the number of edges E can be found from the Euler formula with vanishing genus and two boundaries). Only half surface is shown in order to highlight the curves given by the two sections suggested by the symmetry of the surface.

Two disjoint regions
In this section we discuss the main result of this paper, which is the numerical study of the holographic mutual information of disjoint equal domains delimited by some of the smooth curves introduced in §3.1. For two equal disjoint ellipses, an explicit example of the minimal surface whose area determines the corresponding holographic mutual information is shown in Fig. 10.
Let us consider two dimensional domains A = A 1 ∪A 2 made by two disjoint components A 1 and A 2 , where each component is a simply connected domain delimited by a smooth curve. The boundary is ∂A = ∂A 1 ∪∂A 2 and the shapes of ∂A 1 and ∂A 2 could be arbitrary, but we will focus on the geometries discussed in §3. Since the area law holds also for S A1∪A2 and P A = P A1 + P A2 , the leading divergence O(1/ε) cancels in the combination (1.2), which is therefore finite when ε → 0.
Considering the mutual information (1.2) with the entanglement entropy computed through the holographic formula (1.1), we find it convenient to introduce I A1,A2 as follows where G N is the four dimensional Newton constant. Since ∂A 1 and ∂A 2 are smooth curves, from (2.4) and (3.8) we have In the following we study I A1,A2 when ∂A is made either by two circles ( §4.1.2) or by two superellipses or by the boundaries of two two dimensional spherocylinders. Once A 1 , A 2 and their relative orientation have been fixed, we can only move their relative distance. A generic feature of the holographic mutual information is that it diverges when A 1 and A 2 become tangent, while it vanishes when the distance between A 1 and A 2 is large enough.

Circular boundaries
In this section we consider domains A whose boundary ∂A is made by two disjoint circles. The corresponding disks can be either overlapping (in this case A is an annulus) [51][52][53] or disjoint [54,55].

Annular regions
Let us consider the annular region A bounded by two concentric circles with radii R in < R out . The complementary domain B is made by two disjoint regions and, since we are in the vacuum, S A = S B . The minimal surfaces associated with this case have been already studied in [51,53] as the gravitational counterpart of the correlators of spatial Wilson loops and in [52] from the holographic entanglement entropy perspective.
In §D.2 we discuss the construction of the analytic solution in D dimensions for completeness, but here we are interested in the D = 2 case. Because of the axial symmetry, it is convenient to introduce polar coordinates (ρ, φ) at z = 0. Then, the profile of the minimal surface is completely specified by a curve in the plane (ρ, z).
A configuration providing a local minimum of the area functional is made by the disjoint hemispheres anchored on the circles with radii R in and R out . In the plane (ρ, z), they are described by two arcs centered in the origin with an opening angle of π/2 (see the dashed curve in Fig. 23). Another surface anchored on ∂A that could give a local minimum of the area functional is the connected one having the same topology of a half torus. This solution is fully specified by its profile curve in the plane (ρ, z), which connects the points (R in , 0) and (R out , 0). Thus, we have two qualitatively different surfaces which are local minima of the area functional and we have to establish which is the global minimum in order to compute the holographic entanglement entropy. Changing the annulus A, a transition occurs between these two types of surfaces, as we explain below. This is the first case that we encounter of a competition between two saddle points of the area functional.
The existence of the connected solution depends on the ratio η ≡ R in /R out < 1. As discussed in §D.2, a minimal value η * can be found such that for 0 < η < η * only the disconnected configuration of two hemispheres exists, while for η * < η < 1, besides the disconnected configuration, there are two connected configurations which are local minima of the area functional (see Fig. 23). In the latter case, one has to find which of these two connected surfaces has the lowest area and then compare it with the area of the two disconnected hemispheres. This comparison provides a critical value η c > η * such that when η ∈ (η c , 1) the minimal surface is given by the connected configuration, while for η ∈ (0, η c ) the minimal area configuration is the one made by the two disjoint hemispheres.
Let us give explicit formulas about these surfaces by specifying to D = 2 the results found in §D.2 (in order to simplify the notation adopted in §D.2, in the following we report some formulas from that appendix omitting the index D). The profile of the radial section of the connected minimal surface in the plane (ρ, z) is given by the following two branches where, by introducingz ≡ z/ρ, the functions f ±,K (z) are defined as follows (from (D.11)) (4.4) ∆A ε η c η * Figure 11: Left panel: Radial profiles of the connected surfaces anchored on the boundary of an annulus A which are local minima of the area functional. Comparison between the section of the surfaces constructed with Surface Evolver (black dots) and the analytic expressions reported in §4.1.1. While the external radius is kept fixed to R out = 1, for the internal one the values R in = 0.38 (red), 0.5 (green) and 0.7 (magenta) have been chosen. The cutoff is ε = 0.03 and, according to our regularization prescription, ∂A has been defined at z = ε in the numerical construction. Right panel: The sign of ∆A establishes the minimal area surface between the connected surface and the two disjoint hemispheres. The black curve is obtained from (4.12) by varying K > 0 and it is made by two branches joining at η = η * , where the lower one corresponds to the connected solution which is not the minimal one between the two connected ones. The data points have been found with Surface Evolver for various annular domains. Notice that in the left panel η < η c only for the red curve.
The integral occurring in f ±,K can be computed in terms of the incomplete elliptic integrals of the first and third kind (see §E), finding where we have introduced The matching condition of the two branches (4.3) provides a relation between η η * and the constant K, namely (from (D.13)) where K(m) and Π(n, m) are the complete elliptic integrals of the first and third kind respectively. The relation (4.7) tells us η = η(K) and κ ∈ [1/ √ 2, 1]. As discussed in §D.2, where also related figures are given, plotting this function one gets a curve whose global minimum tells us that η * = 0.367. From this curve it is straightforward to observe that, for any given η ∈ (η * , 1), there are two values of K fulfilling the matching condition (4.7). This means that, correspondingly, there are two connected surfaces anchored on the same pair of concentric circles on the boundary which are both local minima of the area functional. We have to compute their area in order to establish which one has to be compared with the configuration of disjoint hemispheres to find the global minimum.
Performing the following integral up to an additive constant (from (D.20) for D = 2) Plotting the O(1) term of this expression in terms of K, it is straightforward to realize that the minimal area surface between the two connected configurations corresponds to the smallest value of K.
As for the area of the configuration made by two disconnected hemispheres, from (D.23) one gets We find it convenient to introduce ∆A ≡ A dis −A con , which is finite when ε → 0. In particular, ∆A → 2π∆R as ε → 0, where ∆R is (D.27) evaluated at D = 2. From (4.10) and (4.11), we have Considering as the connected surface the one with minimal area, the sign of ∆A determines the minimal surface between the disconnected configuration and the connected one and therefore the global minimum of the area functional. The root η c of ∆A can be found numerically and one gets η c = 0.419 [45,51]. Thus, the connected configuration is minimal for η ∈ (η c , 1), while for η ∈ (0, η c ) the minimal area configuration is the one made by the disjoint hemispheres. By employing Surface Evolver, we can construct the surface anchored on the boundary of the annulus at z = 0 which is a local minimum, compute its area and compare it with the analytic results discussed above. This is another important benchmark of our numerical method.
In the left panel of Fig. 11 we consider the profile of the connected configuration in the plane (ρ, z). The black dots correspond to the radial section of the surface obtained with Surface Evolver, while the solid line is obtained from the analytic expressions discussed above. Let us recall that the triangulated surface is numerically constructed by requiring that it is anchored to the two concentric circles with radii R in < R out at z = ε and not at z = 0, as it should. Despite this regularization, the agreement between the analytic results and the numerical ones is very good for our choices of the parameters. It is worth remarking that, when η η * and therefore two connected solutions exist for a given η, Surface Evolver finds the minimal area one between them. Nevertheless, it is not able to establish whether it is the global minimum. Indeed, for example, the red curve in the left panel of Fig. 11 has η * < η < 1 and therefore the corresponding surface is minimal but it is not the global minimum. Instead, considering an annulus with η < η * , even if one begins with a rough triangulation of a connected surface, Surface Evolver converges towards the configuration made by the two disconnected hemispheres.
In the right panel of Fig. 11 we compare the values of ∆A obtained with Surface Evolver with the analytic curve from (4.12), finding a very good agreement. Numerical points having η * < η < η c are also found, for the reason just explained.

Two disjoint disks
In this section we consider domains A made by two disjoint disks by employing the analytic results for the annulus reviewed in §4.1.1 and some isometries of H 3 . This method has been used in [64] for the case of a circle, while the case of two disjoint circles has been recently studied in [54,55]. The analytic results found in this way provide another important benchmark for the numerical data obtained with Surface Evolver.
Let us consider the following reparameterizations of H 3 , which correspond to the special conformal transformations on the boundary [64] x , When z = 0 in (4.13), the maps (x, y) → (x,ỹ) are the special conformal transformations of the Euclidean conformal group in two dimensions. These transformations in the z = 0 plane send a circle C with center c = (c x , c y ) and radius R into another circle C with centerc = (c x ,c y ) and radius R which are given bỹ . (4.14) Notice that the centerc is not the image of the center c under (4.13) with z = 0. Moreover, when c is such that the denominator in (4.14) vanishes, the circle is mapped into a straight line [64]. z Figure 13: Two examples of minimal surfaces (constructed with Surface Evolver) corresponding to A made by two disjoint and equal disks (∂A is given by the red and blue circles). Only half of the surfaces is shown in order to highlight their section through a plane orthogonal to z = 0 and to the segment connecting the centers. This section provides a circle whose radius and center are given in (4.20). In this figure ε = 0.03, the red circles have radius R = 1 and the distance between their centers is d = 2.16, while for the blue ones R = 0.75 and d = 1.68.
Considering two concentric circles at z = 0 with radii R in < R out , their images are two different circles at z = 0 which do not intersect. In order to deal with simpler expressions for the mapping, let us place the center of the concentric circles in the origin, i.e. c = (0, 0). By introducing η ≡ R in /R out < 1 for the initial configuration of concentric circles centered in the origin and denoting by out | the radii of the circles after the mapping, the distance between the two centers reads where β 2 ≡ |b| 2 R 2 in . Thus, η and β fix the value of the ratioδ ≡ d/ R 1 . The final disks are either disjoint or fully overlapping, depending on the sign of the expression within the absolute value in the denominator of (4.15). In particular, when β 2 ∈ (η 2 , 1) the two disks are disjoint, while when β 2 ∈ (0, η 2 ) ∪ (1, +∞) they overlap. As for their ratioη ≡ R 1 / R 2 , we find (4.16) Notice thatη → 1/η > 1 for β 2 → ∞. Thus, given η and β, the equations (4.15) and (4.16) provideδ and η. By inverting them, one can write η and β in terms ofδ andη. The system is made by two quadratic equations and some care is required to distinguish the various regimes. When the disks after the mapping are disjoint, i.e. η 2 < β 2 < 1, an interesting special case to discuss is R 1 = R 2 , namely when the disjoint disks have the same radius R = R in /(1 − η) = R out /(η −1 − 1), being R in < R out the radii of the two concentric circles at z = 0 centered in the origin. Settingη = 1 in (4.16), one finds that it happens for β 2 = η, i.e. |b| 2 = 1/(R in R out ). The distance corresponding to this value of β can be found from (4.15) and it is given by where the root η(δ) < 1 has been selected andδ > 2 must be imposed in order to avoid the intersection of the two equal disks.
is chosen by fixing the initial and final configurations of circles at z = 0, the transformations (4.13) for the points in the bulk are fixed as well and they can be used to map the points belonging to the minimal surfaces spanning the initial configuration of circles. In particular, let us consider a circle given by (R cos φ, R sin φ, z ) for φ ∈ [0, 2π), lying in a plane at z = z parallel to the boundary. This circle is mapped through (4.13) into another circle C whose radius is given by 17) and whose centerĉ ≡ (ĉ x ,ĉ y ,ĉ z ) has coordinateŝ (4.18) Setting z = 0, R = R and R = R in (4.17) and (4.18), the expressions in (4.14) with c = (0, 0) are recovered. The circle C lies in a plane orthogonal to the following unit vector where 2z |b| R/R < 1, as can be easily observed from (4.17).
In the top left panel of Fig. 12 we consider as initial configuration the annulus at z = 0 for some given value of η and the corresponding connected minimal surface in the bulk anchored on its boundary, which has been discussed in §4.1.1. The transformation (4.13) with β = √ η maps this surface into the connected surface anchored on two equal and disjoint circles (bottom right panel in Fig. 12). It is interesting to follow the evolution of the former surface into the latter one as β ∈ [0, √ η] increases: in Fig. 12 we show two intermediate steps where the surfaces are qualitatively different and they correspond to different regimes of β separated by β = η. For 0 < β < η the disks at z = 0 are still overlapping but they are not concentric (top right panel of Fig. 12). Within this range of β, the radius of the largest disk, which is R out /|1 − β 2 /η 2 |, increases with β and it diverges when as β → η. When η < β √ η, instead, the disks at z = 0 are disjoint and the images of the initial surface through (4.13) are shown in the bottom panels of Fig. 12, where the surface on the left has η < β < √ η, while the one on the right corresponds to the final stage of disjoint equal disks (β = √ η). In Fig. 12 the mapping preserves the color code and we have highlighted the green circle because in the top left panel it corresponds to the circle at z = z m along which the two branches given by (4.3) match, as imposed by the condition (4.7). When β = √ η, this matching circle is mapped into the vertical one shown in the bottom right panel, whose radius R v and whose coordinate z v > R v of its center along the holographic direction are given respectively by where R is the radius of the two equal disjoint disks written above andz m is a function of η (see (4.4) and (4.7)). In Fig. 13 we show two examples of minimal surfaces constructed with Surface Evolver which provide the holographic mutual information of two equal disjoint disks. Considering the section of these surfaces through a vertical plane which is orthogonal to the boundary and to the line passing through the centers of the disks, we find a good agreement with (4.20). As for the finite part of the area, once η and β have been written in terms ofη andδ by inverting (4.15) and (4.16), the limit ε → 0 of either ∆A or I A1,A2 (depending on whether the final disks are either overlapping or disjoint respectively) is given by the r.h.s. of (4.12), where κ = κ(η) is obtained through the numerical inversion of (4.7), being η = η(δ,η) found above.
The special case of two equal disjoint disks corresponds toη = 1 andδ = (1 + η)/ √ η, and therefore the limit ε → 0 of I A1,A2 depends only on the parameterδ, as expected. The relationδ = (1 + η)/ √ η can be used to find the critical distance d c between the centers beyond which the holographic mutual information vanishes and also the distance d * > d c beyond which the connected surface does not exist anymore. They correspond to η c and η * respectively and, in particular, one getsδ c = 2.192 andδ * = 2.256. In order to check that the surfaces obtained through (4.13) are local minima of the area functional, one can compare the analytic results found as explained above against the corresponding surfaces constructed by Surface Evolver. In Fig. 14 we have performed this check for a section profile: the black dots come from the surface obtained as in the bottom right panel of Fig. 12 (notice that the black dots reach z = 0), while the red curve is the section of the corresponding surface constructed by Surface Evolver (see also the red curves in Fig. 10 for a similar construction with different A). In Fig. 15 we have performed another comparison between the analytic expressions and the numerical data of Surface Evolver by computing the holographic mutual information of a domain A made by two equal disjoint disks. The black triangles have been found by mapping the black curve for the annulus in the right panel of Fig. 11 (which is given by the r.h.s. of (4.12)) through η = η(δ) found above. The agreement with the corresponding data obtained with Surface Evolver (red curve) is very good. Notice that, as already observed for the annulus in §4.1.1, also in this case Surface Evolver finds a surface which is a local minimum of the area functional, even if it is not the global minimum. Let us conclude by emphasizing that, while this numerical method is very efficient in finding surfaces which are local minima for the area functional when they exist, it is not suitable for studying the existence of a surface with a given topology.

Other shapes
In §4.1.2 we have considered the holographic mutual information of two disjoint circular domains, for which analytic results are available. When A = A 1 ∪ A 2 is not made by two disjoint disks, analytic results for the corresponding holographic mutual information are not known and therefore a numerical approach could be very useful. Here we employ Surface Evolver to study I A1,A2 (defined in (4.1)) of disjoint regions delimited by some of the smooth curves introduced in §3.1.
The holographic mutual information of non circular domains depends on the geometries of their boundaries, on their distance and also on their relative orientation. Independently of the shapes of ∂A 1 and ∂A 2 , once the domains and their relative orientation have been fixed, the holographic mutual information vanishes Figure 16: Holographic mutual information of two equal and disjoint domains delimited by ellipses (top panels) or superellipses with n = 4 (bottom panels), which are defined by R 1 and R 2 (see the bottom panel of Fig. 1 and (3.9)), while d is the distance between their centers. The relative orientation is like in Fig. 10. Left panels: Density plots for I A1,A2 whose zero provides the corresponding transition curve (solid black line) in the plane (d/R 2 , R 1 /R 2 ). The straight vertical line indicates the transition when A is made by two equal and disjoint infinite strips whose width is 2R 2 and the distance between their central lines is d. Right panels: I A1,A2 in terms of d/R 2 for various fixed values of R 1 /R 2 indicated by the horizontal dashed lines in the corresponding left panel, with the same color code. The lower curves (orange) in the right panels correspond to the squircles (R 1 = R 2 ) with n = 2 (top) and n = 4 (bottom) and therefore they reproduce the red and orange curves in Fig. 15 respectively. The data reported here have been found with R 2 = 1 and some checks have been done also with R 2 = 2. Figure 17: Holographic mutual information of two equal and disjoint two dimensional spherocylinders oriented like the two ellipses in Fig. 10. The parameters R 1 and R 2 specify the domains (see the bottom panel of Fig. 1 and (3.10)) and d is the distance between their centers. The same notation and color coding of Fig. 16 has been adopted.
when the distance between A 1 and A 2 is large enough. The critical distance d c beyond which I A1,A2 = 0 depends on the configuration of the domains. This transition occurs because, for a generic distance d between the centers of A 1 and A 2 , the global minimal area surface comes from a competition between a connected surface anchored on ∂A and a configuration made by two disconnected surfaces spanning ∂A 1 and ∂A 2 , which are both local minima. Beyond the critical distance between the centers, the disconnected configuration becomes the global minimum and therefore I A1,A2 vanishes.
In Fig. 10 we show an example of a connected surface constructed with Surface Evolver where ∂A is made by two equal and disjoint ellipses at z = 0. Let us recall that in our numerical analysis we have regularized the area by defining ∂A at z = ε, as discussed in §B. In the figure, we have highlighted two sections of the surface suggested by the symmetry of this configuration of domains, which are given by the red curves and by the green one.
We have constructed minimal area connected surfaces also for configurations of equal disjoint domains with other shapes and in Fig. 14 we have reported the corresponding curves obtained from the section giving the red curves in Fig. 10. The red curves in Fig. 14 are associated with circular domains and they can be recovered analytically (black dots), as explained in §4.1.2. Instead, for the remaining curves analytic expressions are not available and therefore they provide a useful benchmark for analytic results that could be found in the future.
Besides the profiles for various sections, Surface Evolver computes also the area of the surfaces that it constructs. Considering a configuration of disjoint domains with given shapes and relative orientation, we can compute I A1,A2 while the distance d between their centers changes. In Fig. 15 we show the results of this analysis when ∂A 1 and ∂A 2 are squircles (i.e. (3.9) with R 1 = R 2 ≡ R). As for their relative orientation, drawing the squares that circumscribe ∂A 1 and ∂A 2 , their edges are parallel. Since I A1,A2 0, the critical distance d c corresponds to the zero of the various curves and I A1,A2 vanishes for d d c . Thus, I A1,A2 is continuos with a discontinuous first derivative at d = d c . The points found numerically which have I A1,A2 < 0 correspond to connected surfaces that Surface Evolver constructs but they are not the global minimum for the area functional because the disconnected configuration is favoured for that distance. z z Figure 18: Minimal surfaces obtained with Surface Evolver for a domain A = A 1 ∪ A 2 made by the interior of two disjoint and equal squares. All the squares have the same size but the relative orientation of A 1 and A 2 is different in the two panels.
Once the relative orientation has been chosen, a configuration of two equal and disjoint squircles is completely determined by two parameters: the distance d between the centers and the size R of the squircles. Instead, when A 1 and A 2 are two equal two dimensional spherocylinders or equal domains delimited by two disjoint superellipses and the relative orientation has been chosen, we have three parameters to play with: the distance d between the centers and the parameters R 1 and R 2 which specify the two equal domains (see the bottom panel of Fig. 1). In Fig. 16 we show I A1,A2 for two disjoint domains delimited by ellipses and superellipses with n = 4, whose relative orientation is like in Fig. 10. In the left panels, the black thick curve is the transition curve along which the holographic mutual information vanishes, while the continuos straight line identifies the transition value corresponding to two disjoint infinite strips [42]. Comparing the transition curve in the top left panel with the one in the bottom left panel, it is evident that the one associated with the superellipses having n = 4 is closer to the value corresponding to the infinite strips than the one associated with the ellipses. In Fig. 17 we study I A1,A2 for a domain A made by two equal and disjoint two dimensional spherocylinders. In this case the transition curve is closer to the line corresponding to the transition for two infinite strips with respect to the transition curves of Fig. 16. Nevertheless, from our data we cannot conclude that the transition curve for the two dimensional spherocylinders approaches the value corresponding to the infinite strips as R 1 /R 2 → ∞. It would be interesting to have further data and some analytic argument to understand whether some bounds prevent the transition curves to approach the value associated with the infinite strips for R 1 /R 2 → ∞. Let us remark that the lowest curves (orange) in the right panels of Figs. 16 and 17 correspond to disjoint squircles with n = 2 (i.e. circles) or n = 4 and therefore they reproduce the red and the orange curves of Fig. 15. Configurations of domains having smaller values of d than the ones shown in the plots provide unstable numerical results.
By employing Surface Evolver, we could also study the holographic mutual information of disjoint domains whose boundaries contain corners. In particular, one could take both A 1 and A 2 bounded by polygons, but also A 1 bounded by a smooth curve and A 2 by a polygon. In Fig. 18 we show the minimal area surfaces corresponding to ∂A made by two equal and disjoint squares having different relative orientation. As discussed in §3.2, when ∂A has vertices a further logarithmic divergence occurs after the area law term in the ε → 0 expansion (see (2.5)). If the coefficient of the logarithmic divergence in (2.5) is additive, i.e. B A1∪A2 = B A1 + B A2 for two disjoint regions, then the holographic mutual information is finite. An expression like (3.12) with the sum extended over the vertices of both the components of ∂A is additive, leading to a finite I A1,A2 . Also for these cases we could find plots similar to Figs. 16 and 17 but the curves would not be suitable for a comparison with an analytic formula because of the regularization procedure that we have adopted. Indeed, in our numerical computations ∂A is defined at z = ε and this regularization affects the O(1) term in (2.5) [61], as already mentioned in the closing part of §3.2. Figure 19: Minimal surface corresponding to three disjoint and equal red circles in the plane z = 0 (the z axis points downward). This surface has 13147 vertices and 26624 faces, while the number of edges is given by Euler formula with vanishing genus and 3 boundaries. This kind of surfaces occurs in the computation of the holographic tripartite information for the union of three disjoint disks.

Conclusions
In this paper we have studied the area of the minimal surfaces in AdS 4 occurring in the computation of the holographic entanglement entropy and of the holographic mutual information, focussing on their dependence on the shape of the entangling curve ∂A in the boundary of AdS 4 .
Our approach is numerical and the main tool we have employed is the program Surface Evolver, which allows to construct triangulated surfaces approximating a surface anchored on a given curve ∂A which is a local minimum of the area functional. We have computed the holographic entanglement entropy and the holographic mutual information for entangling curves given by (or made by the union of) ellipses, superellipses or the boundaries of two dimensional spherocylinders, for which analytic expressions are not known. We have also obtained the transition curves for the holographic mutual information of disjoint domains delimited by some of these smooth curves (see Figs. 15,16 and 17), providing a solid numerical benchmark for analytic expressions that could be found in future studies. We focused on these simple examples, but the method can be employed to address more complicated domains.
Besides the fact that the surfaces constructed by Surface Evolver are triangulated, a source of approximation in our numerical analysis is the way employed to define the curve spanning the minimal surface. Indeed, once the cutoff ε > 0 in the holographic direction has been introduced to regularize the area of the surfaces, the numerical data have been found by defining ∂A at z = ε. It would be interesting to understand better this regularization with respect to some other ones and also to decrease ε in a stable and automatically controlled way in order to get numerical data which provide better approximations of the analytic results.
There are many possibilities to extend our work. The most important ones concern black hole geometries and higher dimensional generalizations. An interesting extension involves domains A made by three or more regions (see [65] for some results in two dimensional conformal field theories and [66][67][68] for a holographic viewpoint). In Fig. 19 we show a minimal surface anchored to an entangling curve made by three disjoint circles. The area of this surface provides the holographic entanglement entropy between the union of the three disjoint disks and the rest of the plane, which is the most difficult term to evaluate in the computation of the holographic tripartite information [66]. Another important application of the numerical method employed here involves time dependent backgrounds modelling the holographic thermalization [11,[69][70][71][72][73][74].
Surface Evolver is a useful tool to get numerical results for the holographic entanglement entropy, which can be used to test analytic formulas that could be found in the future. Figure 20: Example of a typical evolution obtained by Surface Evolver in the case of a circular boundary. The initial configuration consists of an octagonal prism composed of 40 triangles (left). The shape is then optimized and refined as described in §B, finding the final configuration given by the rightmost surface, which consists of 10240 triangles and yields F A = 1.99843π whereas F A = 2π is the exact value from the analytic result (3.1). In this example the radius of the circle is R = 1 and ε = 0.03.
with t the tangent vector of ∂γ A , κ n = t ,s (with s the arc lenght) and κ n and κ g the normal and geodesic curvature respectively. Since ∂γ A lies on the z = 0 plane andẑ · N = 0 at z = 0, then N = ±n where the choice of the sign is conventional. By virtue of (A.5) this implies that κ g = 0. Thus ∂γ A is a geodesic over γ A .
An interesting consequence of the previous statement is that the total Gaussian curvature of the surface is constant, regardless the shape of the boundary in the z = 0 plane. The Gauss-Bonnet theorem tells us thatˆγ where K G is the Gaussian curvature and χ is the Euler characteristic. Since κ g = 0 in our case, we havêγ Let us recall that the Euler characteristic is χ = 2 − 2g − b, where g is the genus of the surface and b is the number of its boundaries.

B Numerical Method
The numerical results presented in §3 and §4 have been obtained with Surface Evolver [56,57]. This is a multipurpose shape optimization program created by Brakke [56] in the context of minimal surfaces and capillarity and then expanded to address generic problems on energy minimizing surfaces. A surface is implemented as a simplicial complex, i.e. a union of triangles. Given an initial configuration of the surface, the program evolves the surface toward a local energy minimum by a gradient descent method. The energy used in our calculations is the H 3 area function given in (2.3). The initial configuration is preferably very simple and contains only the least number of triangles necessary to achieve a given surface topology (Fig. 20). A typical evolution consists in a sequence of optimization and mesh-adjustment steps. During an optimization step, the coordinates of the vertices are updated by a local minimization algorithm (conjugate gradient in our case), resulting in a configuration of lower energy. The topology of the mesh (i.e. the number of vertices, faces and edges) is not altered during minimization. A mesh-adjustment step, on the other hand, consists of a set of operations whose purpose is to render the discretized surface smooth and uniform. These operations can be broadly divided in two class: meshrefinements and mesh-repairs. In a mesh-refinement operation a finer grid is overlaid on the coarse one. This is obtained, for instance, by splitting a triangle in four smaller triangle obtained by joining the mid points of the original edges. In a mesh-repair operation, the triangles that are too distorted compared to the average are eliminated. This operation can change the topology of the mesh and possibly also the topology of the surface which can then breakup into two or more connected parts. This happens, for instance, in the ! " "# ! " "" ! " "$ ! " "#% Figure 21: The quantity F A (see (3.8)) computed with Surface Evolver for ellipses having R 1 = 2R 2 (see the bottom panel of Fig. 1), for various R 2 and ε. When ε/R 2 is too small, our numerical data are not stable. The fitted value on the vertical axis is 3.728. case of the surfaces described in §4. As explained, the minimal surface spanning a disconnected boundary curve can be either connected or disconnected depending on the shape of the boundary. Evolving an initially connected surface in the regime of geometric parameters where the only stable solution is disconnected causes the surface to form narrow necks and eventually pinches off once the triangles around the necks become too squeezed.
Due to the divergence of the area element dA = √ h/z 2 du 1 du 2 at z = 0, the boundary curves used in the numerical work have been defined on the plane z = ε. In order to maximize the accuracy of the numerical solution, it is preferable to choose value of ε that is much smaller than any other length scale in the problem and yet large enough to allow the convergence of the optimization steps. With this goal in mind, we have adopted an empirical selection criterion based on the following procedure. Let ∂γ A be an ellipse and let R 1 and R 2 = R 1 /2 be the semi-major and semi-minor axes. Using Surface Evolver we have calculated the finite part of the area F A for various choices of ε and R 1 . In the limit of ε → 0 the ratio F A /R 1 is expected to approach a finite value, but from the data shown in Fig. 21 we see that for ε/R 2 < 0.02, the accuracy of the numerical calculation starts to drop. Based on this numerical evidence we have set in most of our numerical calculations ε/R = 0.03, where R is the typical length scale of the boundary. It is worth remarking that in our numerical computations it is easier (namely the evolution is more stable) to deal with smaller values of ε/R by increasing R than by decreasing ε. Smaller values of ε/R obtained by decreasing ε keeping R fixed can be achieved by setting up ad hoc evolutions, tailored for a specific type of boundary shape. This has been done only for the triangles in Fig. 4, while in the remaining figures we have increased R keeping ε = 0.03 fixed. Nevertheless, for ε fixed, numerical instabilities are encountered when R is too large as well. The values of ε/R adopted in our numerical calculations have been chosen to guarantee both stable evolutions and a satisfactory precision to compare the data with the analytic results, when they are available.
Other alternative methods are available to construct minimal surfaces. A popular one by Chopp [76] consists of evolving the surface level sets under the surface mean curvature flow. A variant of this method has been employed in [50] to study minimal surfaces in the Schwarzschild-AdS D+2 background.

C Superellipse: a lower bound for F A
In this appendix we provide a lower bound for the quantity F A (see (2.4)) associated with the entangling curves ∂A given by the superellipses (3.9), that we have discussed in §3.1.
If A is a simply connected domain without corners in its boundary, let us consider a surface γ * A anchored on ∂A, but different fromγ A , and such that A[γ * Beingγ A the minimal area surface anchored on ∂A, it is immediate to realize that F * A < F A . Here we consider the superellipses (3.9), whose perimeter is given by where the integration variablex = x/R 1 as been employed. Let us adapt to this case the choice of the trial surface suggested in [28] for the ellipse, namely we consider γ * A such that any section along the x direction provides the profile of the infinite strip whose width is given by y(x) obtained from (3.9), i.e.
Given the symmetries of the superellipse, we are allowed to restrict ourselves to x > 0 and y > 0. From (D.2) for D = 2, we construct the trial surface γ * A by requiring that we have that any section at x = const is given by where the integration variable Z ≡ z/z * has been employed and z * (x) has been introduced by taking z * in (3.5) with s ∞ defined in (3.6) and replacing R 2 with y(x) defined in (C.2). From (C.3), it is straightforward to show that y(0,x) = y(x) and this guarantees that the trial surface is anchored on the superellipse (C.2). The occurrence of the cutoff ε in the holographic direction influences the integration domain along the x direction. In particular, by employing (C.2) and (C.3), the requirement z * (x) ε becomesx x ε , wherẽ Plugging (C.3) inside the area functional, being y written in terms of x and z, we get Computing (C.5) analytically is too hard, but one can check that the area law is satisfied. When ε → 0, from (C.4) we have thatx ε = 1 + O(ε n ). In this limit, the most divergent term of M ε (x) comes from the limit of integration ε/z * (x) and it can be found by considering an integration on the interval [ε/z * (x), a], where Z is infinitesimal if a 1. The remaining integral provides O(1) terms. For Z → 0 we have that C(0) = 1 and therefore the leading term in (C.5) is given by where P A given in (C.1) can be recognized after (C.3) and (C.2) have been employed. We are not able to find F * A analytically but it can be obtained numerically as F * given by (C.5), getting a lower bound for F A associated with the superellipse.
It is interesting to consider F * A in the limit of a very elongated superellipses, namely when R 1 /R 2 → ∞. This means that (C.5) must be studied in the double expansion ε → 0 and R 2 /R 1 → 0. Assuming that the order of this two limits does not matter, let us set R 2 /R 1 = 0 in the expressions of M ε (x) in (C.6) and expand it for small ε, finding where z * (x) is given in (C.3). By plugging (C.8) into (C.5) and expanding the resulting expression for ε → 0, we have that Notice that, from (C.1), one can observe that P A = 4R 1 1 + o(1) when R 1 /R 2 → ∞. We conclude that the leading term of F * A as R 1 /R 2 → ∞ reads F * A = πs ∞ n sin(π/n) When n = 2, the result of [28] is recovered, as expected. Moreover, the expression (C.10) in the special cases of n = 2 and n = 3 has been checked in Fig. 2 against the data obtained with Surface Evolver (see respectively the red and the blue dotted horizontal lines), finding a good agreement. Notice that the expression in the r.h.s. of (C.10) is strictly larger than the value of F A corresponding to the infinite strip (see (3.6)), which is approached as n → ∞.

D.1 Sections of the infinite strip
In this section we discuss the computation of the area of the domain identified by an orthogonal section of the minimal surfaces associated with the infinite strip. The metric of AdS D+2 in the Poincaré coordinates reads ds 2 = − dt 2 + dz 2 + dx 2 1 + · · · + dx 2 D z 2 . (D.1) Considering an infinite D-dimensional strip on the spatial slice t = const extended along the x 2 , . . . , x D directions whose width is given by 2R 2 , i.e. |x 1 | R 2 , the minimal area surface associated with this domain is characterized by the profile z = z(x 1 ). Because of the symmetry of the problem, z(x 1 ) is even and therefore we can restrict to 0 x 1 R 2 . The profile is obtained by solving the following differential equation [9,10] where z * is the maximum value of z, which is reached at x 1 = 0. A way to get an orthogonal section of the infinite strip is defined by x 2 = · · · = x D = const. Then, one considers the two dimensional region enclosed by the profile z(x 1 ) and the cutoff z = ε in the plane (x 1 , z). The domain along the x 1 axis is |x 1 | R 2 − a, where a is defined by z(R 2 − a) = ε. Its area readŝ where (D.2) has been employed. Another section of the infinite strip to study is defined by x i = const for some 2 i D and |x j | R 1 for j = i. In this case we are interested in the volume of the D dimensional region enclosed by the profile z(x 1 ) and z = ε, whose projection on the z = 0 hyperplane is included within the section of the infinite strip we are dealing with. It is given bŷ For any D, a minimal value η * > 1 occurs, which is shown in the inset. Given a value η ∈ (η * , 1), two values of K correspond to it, providing two different radial profiles (see an example for D = 2 in Fig. 23).
which correspond to two different parts of the profile. As for the integration constant K, it must be strictly positive becausez = 0 corresponds to the boundary z = 0, which is included in the range of z. The domain forz is 0 z z m , wherez m is the first positive zero of the polynomial under the square root in (D.10). For D = 2 we are lead to solve a biquadratic equation, which givesz 2 m = K + K(K + 4) /2. Notice that z m → 0 when K → 0.
The differential equation (D.10) can be solved through the separation of the variables. In particular, from the r.h.s. of (D.10), we find it convenient to introduce Imposing that these two branches match at the point P m , whose (ρ, z) coordinates are (ρ m , z m ≡ z(ρ m )), where z m has been found above, we get the following relation Sincez m depends on K, from (D.13) we get a relation between η and K, which is represented in Fig. 22 for 2 D 7. The first feature to point out about (D.13) is the existence of a minimal value for η that will be  Figure 23: Radial profiles in the (ρ, z) plane for the connected surfaces anchored on the boundary of the same annulus A having R in < R out . They correspond to local minima of the area functional and they are characterized by the two different values of K associated with the same η. These connected surfaces are obtained through (D.12) and (D.13), where different colours are used for the various branches. The dashed curves represent the two concentric hemispheres anchored on ∂A as well. The continuous grey curves are the paths in the (ρ, z) plane of the points P 0 , P m and P * as K ∈ (0, ∞). Here D = 2, R in = 0.43, R out = 1 and the values of K are K = 0.81 (global minimum) and K = 2.05 (local minimum). Comparing the area of the two connected surfaces, we find that the one having minimal area has P * closer to the boundary. denoted by η * > 0. For instance, we find η * = 0.367, η * = 0.542 and η * = 0.643 for D = 2, D = 3 and D = 4 respectively (see the inset in Fig. 22 for other D's). Then, for any η * < η < 1, there are two values of K giving the same η, while for 0 < η < η * connected solutions do not exist. The two different K's associated with the same η * < η < 1 provide two different radial profiles and therefore two connected surfaces having the same ∂A. In order to find the global minimum of the area functional, we have to evaluate their area. Through a numerical analysis, one observes that z m is an increasing function of K.
Given R in and R out , besides the two connected surfaces having the same η but different K, we have also another surface γ A which is a local minimum for the area functional (D.6) such that ∂γ A = ∂A: it is made by two disjoint concentric hemispheres in the bulk with radii R in and R out which are anchored on the boundaries of the concentric spheres in the boundary (see the dashed curves in Fig. 23). The area of a hemisphere of radius R in the bulk anchored on the boundary of a sphere with the same radius at z = 0 can be found by integrating (D.6) for 0 ρ R − a, where z(ε) ≡ a, finding

E Elliptic integrals
When D = 2, the integrals encountered in §3.2 and in §D.2 can be computed analytically in terms of elliptic integrals. Here we report their definitions for completeness, following [77] (notice that Mathematica adopts the same notation).
The incomplete elliptic integrals of the first, second and third kind are defined respectively as follows which are the complete elliptic integrals of the first, second and third kind respectively.