1 Introduction↩︎

Entanglement entropy plays a vital role in quantum information theory and quantum many-body systems, as it quantifies the degree of quantum entanglement in bipartite pure states. The comprehensive exploration is essential for understanding the underlying structure and the dynamics of quantum systems.

For a region \(\boldsymbol{A}\), entangled with its complement, the entanglement entropy is given by the von Neumann entropy, \(S_{\boldsymbol{A}} = -\text{Tr}(\rho_{\boldsymbol{A}} \log \rho_{\boldsymbol{A}})\), where \(\rho_{\boldsymbol{A}}\) denotes the reduced density matrix obtained by tracing out the degrees of freedom external to \(\boldsymbol{A}\). Within the framework of the AdS/CFT correspondence [1]–[3], classical spacetime geometry in the bulk is dual to a quantum state defined on the conformal boundary. The Ryu-Takayanagi (RT) formula [4], [5], and its covariant generalization [6], provide a geometric realization of entanglement entropy in this context: the entanglement entropy of a region \(\boldsymbol{A}\) on the boundary is proportional to the area of a minimal (or extremal) surface in the bulk that is homologous to \(\boldsymbol{A}\) and anchored on the boundary \(\partial\boldsymbol{A}\). Moreover, once the quantum fields in the bulk are considered, the RT surface should be generalized to a QES \(\gamma_{\boldsymbol{A}}\). The holographic entanglement entropy (HEE) of \(\boldsymbol{A}\) is calculated by extremizing the generalized gravitational entropy functional [7]–[9].

Recent breakthroughs in the black hole information paradox [10]–[15] have sparked growing interest in composite quantum systems comprising a \((d-1)\)-dimensional quantum system \(\partial\mathcal{B}\) coupled to a \(d\)-dimensional heat bath \(\partial\) [16]–[19]. This framework is commonly referred to as the boundary perspective. In its holographic dual, the \((d-1)\)-dimensional system \(\partial\mathcal{B}\) corresponds to a (quantum) gravity theory in a \(d\)-dimensional asymptotically Anti-de Sitter (AAdS) spacetime \(\mathcal{B}\) — a viewpoint often termed as the brane perspective. A key development in this context is applying the QES \(\gamma_{\partial\mathcal{B}}\) to describe the entanglement entropy of this \((d-1)\)-dimensional system \(\partial\mathcal{B}\) [10]–[15]. In this setting, quantum fields propagating in both the gravity and bath regions are typically described by the same conformal field theory (CFT), and their entanglement entropy are computed using the island formula [10], [20], which prescribes extremizing the generalized entropy over all candidate island regions and selecting the configuration that yields the minimal value. Direct computation of the entanglement entropy for quantum fields on curved spacetimes is often technically challenging [12]. However, when the quantum fields are described by a CFT with a large central charge, the framework of doubly holographic duality becomes applicable. This duality further maps the \(d\)-dimensional AAdS spacetime \(\mathcal{B}\) into a Planck brane embedded in a \((d+1)\)-dimensional classical bulk spacetime. The approach is commonly referred to as the bulk gravity perspective. In this scenario, the QES \(\gamma_{\partial\mathcal{B}}\) reduces to the standard RT surface \(\Gamma_{\partial\mathcal{B}}\), thereby enabling a purely geometric description of the entanglement entropy [21]–[25]. ¹

To investigate the entanglement structure between the \((d-1)\)-dimensional system \(\partial\mathcal{B}\) and its external environment, or equivalently the entanglement contribution from the quantum fields propagating in the curved spacetime background \(\mathcal{B}\), the RT surface corresponding to \((d-1)\)-dimensional regions \(\boldsymbol{A}\subset\partial\mathcal{B}\) must be constructed. However, for \(d\geq 3\), this task becomes significantly more challenging due to the complexity of identifying extremal surfaces near the Planck brane.

For a simple connected region \(\boldsymbol{A}\), the corresponding RT surface can be numerically constructed as a two-dimensional non-uniform extremal surface, by solving a system of partial differential equations in appropriately chosen coordinates [90]. This method offers a concrete characterization of the entanglement entropy associated with a single region \(\boldsymbol{A}\) on \(\partial\mathcal{B}\). In the semi-classical limit, the entanglement entropy of the boundary region \(\boldsymbol{A}\) in \(d=3\) is found to contain a leading contribution with linear divergence and a subleading term with logarithmic divergence. Remarkably, this expression corresponds to the area of the classical extremal surface \(\Gamma_{\boldsymbol{A}}\) in the bulk and admits two physically equivalent interpretations. From the brane perspective, the leading divergence in the entanglement entropy arises from the linear-law contribution associated with the entanglement between the brane CFT and the bath CFT, while the subleading divergence encodes the geometric contribution from the QES \(\gamma_{\boldsymbol{A}}\) localized on the brane. From the boundary perspective, the leading divergence reflects the volume-law entanglement between \(\boldsymbol{A}\) and the bath \(\partial\), whereas the subleading divergence arises from the entanglement between \(\boldsymbol{A}\) and the rest of the boundary \(\partial\mathcal{B}\).

For multiple disconnected regions, an analogous formula has not yet been established, despite its importance for studying quantum information measures such as \(n\)-partite information. In this work, we begin with the simplest case — BMI. Our goals are twofold: first, to analyze how BMI varies with the size of subregions; and second, to investigate the impact of quantum field entanglement on the brane on the holographic BMI. In this scenario, however, the aforementioned method becomes inadequate due to the difficulty in identifying a coordinate system capable of capturing disjoint configurations. To further investigate this mixed state entanglement of quantum fields on the brane, we adopt a shape optimization method using the Surface Evolver [91], [92], a numerical tool widely used in material science for minimizing surface energy via gradient descent techniques. This method constructs minimal surfaces by iteratively evolving an initial surface toward its area-minimizing configuration. It has been successfully applied in holographic settings [93]–[97], making it well-suited for analyzing entanglement structures.

The remainder of this paper is organized as follows. In Section 2, we introduce the doubly holographic setup. We begin by presenting three equivalent perspectives for the gravitational background with a single Planck brane. Subsequently, we describe three corresponding perspectives for entanglement entropy and BMI. Finally, we examine the universality of holographic entanglement entropy and BMI within the framework of double holography, particularly in the semiclassical gravity limit.

In Section 3, we perform a numerical analysis on the holographic entanglement entropy and BMI for various configurations. We also investigate the correction to BMI arising from quantum fields on the brane.

In Section 4, we briefly review the framework of RTNs and present a RTN-based interpretation of the correction term in the BMI.

Section 5 concludes the paper with a summary of our findings and a discussion of potential future directions.

2 The setup in double holography↩︎

Consider a general \((d+1)\)-dimensional asymptotic AdS spacetime truncated by a \(d\)-dimensional Planck brane \(\mathcal{B}\). This brane intersects with the conformal boundary \(\partial\) at infinity, forming a \((d-1)\)-dimensional junction \(\partial\mathcal{B}\) [10], [12], [21], [26]. This setup admits three equivalent descriptions within the framework of double holography.

2.1 The gravitational background↩︎

Within the doubly holographic framework [10], [21], the first description, referred to as the bulk gravity perspective, is governed by the \((d+1)\)-dimensional gravitational action: \[\begin{align} \label{eq:Action} I=&\frac{1}{16\pi G_N^{(d+1)}} \Bigg[ \int d^{d+1}x \sqrt{-g}\left(R+\frac{d(d-1)}{L^2}\right)+2\int_{\partial}d^{d}x\sqrt{-h_{\partial}}K_{\partial}\nonumber\\ &+ \int_{\mathcal{B}}d^{d}x\sqrt{-h}\left(K-\alpha\right)-\int_{\partial\mathcal{B}} d^{d-1}x \sqrt{-\Sigma}\, \theta_{0}\Bigg], \end{align}\tag{1}\] where \(h_{\partial}\), \(K_{\partial}\) denote the induced metric and the extrinsic curvature on the conformal boundary \(\partial\). \(h\), \(K\) correspond to the induced metric and the extrinsic curvature on the Planck brane \(\mathcal{B}\). \(\alpha\) characterizes the brane tension. The final term is the junction condition at the intersection \(\partial\mathcal{B}\), with \(\theta_{0}\) specifying the dihedral angle between the brane \(\mathcal{B}\) and the conformal boundary \(\partial\), and \(\Sigma\) representing the induced metric on the intersection \(\partial\mathcal{B}\). The tension on the brane can be further modified by introducing a Dvali-Gabadadze-Porrati term [98], [99], enabling controlled tuning of the Newton constant ratio between the brane and the bulk [21], [29], [90], [100].

This bulk gravitational description admits two equivalent dual interpretations:

Brane perspective

The \((d+1)\)-dimensional gravity is dual to semi-classical gravity on the brane \(\mathcal{B}\), coupled to CFTs living on both the brane \(\mathcal{B}\) and the heat bath \(\partial\) [10], [101].

Boundary perspective

The combined gravity-plus-bath theory is further dual to lower-dimensional-one quantum system \(\partial\mathcal{B}\) coupled to a heat bath \(\partial\) [10], [21], [26].

To realize dynamical gravity on the brane, we impose Neumann boundary conditions on the Planck brane \(\mathcal{B}\) (see also [28], [102], [103] for alternative viewpoints): \[\label{eq:BCSonBrane95General} K_{a b}-K h_{ab}+\alpha h_{ab}=0,\tag{2}\] where \(h_{ab}\) is the induced metric on the brane \(\mathcal{B}\).

For simplicity, we restrict our analysis to pure \(\text{AdS}_{4}\) spacetime with a single brane, corresponding to the vacuum state of the quantum fields living on the boundary. In Poincaré coordinates \((z, x, y, t)\), the metric is given by \[\begin{align} \label{eq:AdSmetric} ds^2=L^2\frac{-dt^2+dz^2+dx^2+d y^2}{z^2}, \end{align}\tag{3}\] where \(L\) is the AdS radius. The brane is translationally invariant along the \(y\)-direction – Fig. 1, with the position at \[x \tan \theta_{0} + z = 0.\] Substitute this constraint into 3 , the induced geometry on the brane is just an AdS\(_3\) spacetime, with the metric being given by \[ds^2=L^2\frac{-dt^2+dz^2+d y^2}{z^2},\] where \(L\to L \sin\theta_0\), \(t\to t/\sin\theta_0\) and \(y \to y/\sin\theta_0\). The Neumann boundary condition determining the brane tension (2 ) takes the following explicit form in AdS spacetime²: \[\label{eq:BCSonBrane95AdS} \alpha = \frac{2}{L}\cos \theta_{0},\tag{4}\]

With this geometry fully specified, we proceed in subsequent sections to construct the RT surfaces corresponding to boundary regions, which encode the BMI.

2.2 The entanglement entropy and BMI↩︎

From the boundary perspective, let us consider a \(1\)-dimensional spatial region \(\boldsymbol{A}\subset\partial\mathcal{B}\) within a \((2+1)\)-dimensional quantum system \(\partial\mathcal{B}\cup\boldsymbol{\partial}\) as shown in Fig. 2 (a). Its entanglement entropy admits two equivalent perspectives within the holographic framework:

1. Brane perspective – Fig. 2 (b): The entanglement entropy of \(\boldsymbol{A}\) is holographically given by the island formula \[\label{eq:EEfromBrane} \mathcal{S}[\boldsymbol{A}]= \min_{\gamma_{\boldsymbol{A}}}\left\{ \frac{\boldsymbol{Area}(\gamma_{\boldsymbol{A}})}{4G_{\boldsymbol{eff}}^{(3)}} + S_{\text{QFT}}(\mathcal{W}_{\boldsymbol{A}})\right\},\tag{5}\] where \(\gamma_{\boldsymbol{A}}\) denotes the quantum extremal surface on the AdS\(_3\) background. The first term on the r.h.s. \[\mathcal{S}_g[\boldsymbol{A}]:=\frac{\boldsymbol{Area}(\gamma_{\boldsymbol{A}})}{4G_{\boldsymbol{eff}}^{(3)}}\] captures the geometric contribution, while the second term \[\mathcal{S}_c[\boldsymbol{A}]:=S_{\text{QFT}}(\mathcal{W}_{\boldsymbol{A}})\] accounts for the entanglement of the quantum fields within the Q-EW \(\mathcal{W}_{\boldsymbol{A}}\). Here, \(G_{\boldsymbol{eff}}^{(3)}\) is the effective Newton constant on the brane \(\mathcal{B}\). The configuration of the quantum extremal surface \(\gamma_{\boldsymbol{A}}\) has been fully determined by the extremization procedure, and if multiple candidates \(\gamma_{\boldsymbol{A}}\) locally minimize the entropy, one should choose the one that achieves the global minimum.

2. Bulk gravity perspective – Fig. 2 (c): Both terms in (5 ) admit a second holographic dual description to a classical RT surface \(\Gamma_{\boldsymbol{A}}\) in the higher-dimensional one bulk, with the formula to be \[\label{eq:QES} \mathcal{S}[\boldsymbol{A}]=\min_{\Gamma_{\boldsymbol{A}}}\left\{\frac{\boldsymbol{Area}(\Gamma_{\boldsymbol{A}})}{4G_N^{(4)}}\right\}.\tag{6}\] The right-hand side term corresponds to the area of a classical extremal surface \(\Gamma_{\boldsymbol{A}}\) that extends into the bulk and the surface is bounded by the union of the region \(\boldsymbol{A}\) and the QES \(\gamma_{\boldsymbol{A}}\) on the brane. The entanglement entropy of \(\boldsymbol{A}\) is then obtained by selecting the surface \(\Gamma_{\boldsymbol{A}}\) associated with the globally minimal generalized entropy among all admissible QES \(\gamma_{\boldsymbol{A}}\) configurations.

Now consider a region \(\boldsymbol{A}=\boldsymbol{A_1}\cup\boldsymbol{A_2}\subset\partial\mathcal{B}\), composed of two disjoint, simply connected subregions \(\boldsymbol{A_1}\) and \(\boldsymbol{A_2}\), each residing on the boundary quantum system \(\partial\mathcal{B}\), respectively. The BMI between \(\boldsymbol{A_1}\) and \(\boldsymbol{A_2}\) is defined as \[\begin{align} \label{eq:mutualinformation} \mathcal{I}[\boldsymbol{A_1}:\boldsymbol{A_2}] &:= \mathcal{S}[\boldsymbol{A_1}] +\mathcal{S}[\boldsymbol{A_2}] -\mathcal{S}[\boldsymbol{A_1}\cup \boldsymbol{A_2}]. \end{align}\tag{7}\] This expression also admits two equivalent perspectives:

1. Brane perspective: The BMI naturally decomposes into geometric contribution and quantum field contribution separately as \[\begin{align} \mathcal{I}[\boldsymbol{A_1}:\boldsymbol{A_2}] =& \sum_{i=1}^{2}\Bigg\{\min_{\gamma_{\boldsymbol{A_i}}}\left[\frac{\boldsymbol{Area}(\gamma_{\boldsymbol{A_i}})}{4\;G_{\boldsymbol{eff}}^{(3)}}+S_{\text{QFT}}(\mathcal{W}_{\boldsymbol{A_i}})\right]\Bigg\} \nonumber\\ &-\min_{\gamma_{\boldsymbol{A_1}\cup\boldsymbol{A_2}}}\left[\frac{\boldsymbol{Area}(\gamma_{\boldsymbol{A_1}\cup\boldsymbol{A_2}})}{4\;G_{\boldsymbol{eff}}^{(3)}}+S_{\text{QFT}}(\mathcal{W}_{\boldsymbol{A_1}\cup\boldsymbol{A_2}})\right], \end{align}\] where the geometric contribution is given by \[\label{eq:AreaTerm} \mathcal{I}_g[\boldsymbol{A_1}:\boldsymbol{A_2}]:=\frac{\boldsymbol{Area}(\gamma_{\boldsymbol{A_1}})}{4\;G_{\boldsymbol{eff}}^{(3)}}+\frac{\boldsymbol{Area}(\gamma_{\boldsymbol{A_2}})}{4\;G_{\boldsymbol{eff}}^{(3)}}-\frac{\boldsymbol{Area}(\gamma_{\boldsymbol{A_1}\cup\boldsymbol{A_2}})}{4\;G_{\boldsymbol{eff}}^{(3)}},\tag{8}\] while the quantum field contribution is \[\label{eq:CorrectionTerm} \mathcal{I}_c[\boldsymbol{A_1}:\boldsymbol{A_2}]:=S_{\text{QFT}}(\mathcal{W}_{\boldsymbol{A_1}})+S_{\text{QFT}}(\mathcal{W}_{\boldsymbol{A_2}})-S_{\text{QFT}}(\mathcal{W}_{\boldsymbol{A_1}\cup\boldsymbol{A_2}}).\tag{9}\] Each QES \(\gamma(\mathcal{K})\), with \(\mathcal{K} \in \{\boldsymbol{A_1},\boldsymbol{A_2},\boldsymbol{A_1}\cup\boldsymbol{A_2}\}\) is determined by extremizing the corresponding entire entropy functional. The global minimum among all candidates is selected. Notably, the phase transition of BMI in this setting depends not only on geometric contributions, as in classical gravity, but also on the quantum field contributions.

2. Bulk gravity perspective: The combination of quantum extremal surfaces is further dual to the classical extremal surfaces in the higher-dimensional bulk as \[\mathcal{I}[\boldsymbol{A_1}:\boldsymbol{A_2}]= \sum_{i=1}^2\left\{\min_{\Gamma_{\boldsymbol{A_i}}}\left[\frac{\boldsymbol{Area}(\Gamma_{\boldsymbol{A_i}})}{4 G_N^{(4)}}\right]\right\} -\min_{\Gamma_{\boldsymbol{A_1}\cup\boldsymbol{A_2}}}\left[\frac{\boldsymbol{Area}(\Gamma_{\boldsymbol{A_1}\cup \boldsymbol{A_2}})}{4 G_N^{(4)}}\right],\] where \(\Gamma(\mathcal{K})\) denotes the RT surface of the subregion \(\mathcal{K} \in \{\boldsymbol{A_1},\boldsymbol{A_2},\boldsymbol{A_1}\cup\boldsymbol{A_2}\}\) in the bulk.

In the next subsection, we will further analyze the universal behavior of BMI in the semiclassical limit and summarize all associated formulas proved numerically in Sec. 3.

2.3 The universality in the semiclassical limit↩︎

In quantum information theory, BMI is a fundamental quantity for characterizing correlations in quantum systems, including quantum entanglement. Within the double holographic framework, BMI similarly serves to quantify entanglement between two black holes or between two SYK clusters [100]. In this section, we systematically study the BMI of codimension-three subregions. We begin by reviewing the semiclassical limit and the entropy formula for a single subregion. We then discuss the universal entanglement behavior and summarize the entropy and BMI formulas for multiple subregions in the doubly holographic setup, which will be obtained through numerical analysis in the following sections.

Before proceeding, it is crucial to note that for any boundary region \(\boldsymbol{A} \subset \partial\mathcal{B}\), the entanglement entropy undergoes a phase transition as a function of the brane parameter – in our case, the dihedral angle \(\theta_{0}\) (see [100] for more general cases). In particular, only in the subcritical regime (\(\theta_{0} < \theta_c\simeq0.64\)) does the entanglement entropy become non-vanishing, making this the only case meaningful for studying the entanglement properties of \(\boldsymbol{A}\) [16]. To assure this condition is met, our analysis mainly focuses on the semiclassical limit \(\theta_{0} \to 0\), which naturally lies within this subcritical regime. In this limit, the induced gravitational action on the brane takes the form [21]–[24]: \[I_{\text{eff}}=I_{b}+I_{\text{reg}},\] where \(I_{b}\) corresponds to the terms of the bulk action (1 ) on the brane: \[I_{b}= \frac{-1}{8\pi G_N^{(4)}} \int_{\mathcal{B}}d^{3}x\sqrt{-h}\alpha,\] and \(I_{\text{reg}}\) is given by [21], [104] \[I_{\text{reg}}=\frac{1}{16 \pi G_{N}^{(4)}} \int_{\mathcal{B}} d^d x \sqrt{-h}\left[\frac{4}{L}+L R_h \right]+\mathcal{O}[R_h]^2.\] Combining the above terms, the effective gravitational action on the brane can be expressed as \[\label{eq:HCgravity} I_{\text{eff}}=\frac{1}{16 \pi G_{\mathrm{eff}}^{(3)}} \int d^d x \sqrt{-h}\left[\frac{2}{\ell_{\mathrm{eff}}^2}+R_h\right]+\mathcal{O}[R_h]^2,\tag{10}\] with the effective Newton constant and AdS\(_3\) radius scale defined as \[\label{eq:effecitve95newton95constant} \frac{1}{G_{\mathrm{eff}}^{(3)}}=\frac{L}{G_{N}^{(4)}},\quad\text{and}\quad \frac{1}{l_\text{eff}^2}=\frac{2-\alpha L}{L^2}.\tag{11}\] In the semiclassical limit, the ratio of the central charges admits a bulk geometric interpretation through the relation as \[\begin{align} \label{eq:ratio95of95central95charges} \frac{c'}{c}& \simeq 6\sqrt{\frac{1}{2-2\cos \theta_{0}}},\nonumber\\ \text{with}\quad c' &:= \frac{3}{2}\frac{l_{\text{eff}}}{G_{\text{eff}}^{(3)}} \quad\text{and}\quad c :=\frac{L^{2}}{4G_{N}^{(4)}}, \end{align}\tag{12}\] where \(c'\) denotes the central charge of the boundary CFT on \(\boldsymbol{\partial} \mathcal{B}\), and \(c\) corresponds to the central charge of the bath CFT on \(\partial\).

For a one-dimensional single region \(\boldsymbol{A}=\boldsymbol{A_i}\subset \partial\mathcal{B}\), the corresponding entropy formula has been confirmed to take the form as [90] \[\begin{align} \label{eq:entanglemententropy} \mathcal{S}[\boldsymbol{A_i}]=&c \frac{l_i}{\epsilon}+\frac{c'}{3}\log \frac{l_i}{\delta}+\mathcal{F}_{\boldsymbol{A_i}}, \quad \text{with } i=1,2 \end{align}\tag{13}\] where \(l\) is the length of \(\boldsymbol{A_i}\), \(c\) denotes the central charge of the CFT\(_3\) on both the brane \(\mathcal{B}\) and bath \(\partial\), while \(c'\) represents the central charge of the CFT\(_2\) on the boundary \(\partial\mathcal{B}\). Moreover, \(\epsilon\) and \(\delta=\csc\theta_0 \epsilon\) are the UV-cutoffs of the CFT\(_3\) and CFT\(_2\), respectively [19], [21]. It is important to note that while the entire entropy formula in 13 is derived from computing the area of the classical extremal surface \(\Gamma_{\boldsymbol{A_i}}\), the geometric and quantum components can be isolated in the semiclassical limit as \[\mathcal{S}_g[\boldsymbol{A_i}] \to \frac{c'}{3}\log \frac{l_i}{\delta}, \quad\text{and}\quad\mathcal{S}_c[\boldsymbol{A_i}] \to c \;\frac{l_i}{\epsilon} + \mathcal{F}_{\boldsymbol{A_i}}.\] From the brane perspective, \(\mathcal{S}_g[\boldsymbol{A_i}]\) arises from the area of the QES \(\gamma_{\boldsymbol{A_i}}\), encoding the geometric contribution from the induced gravity on the brane, while \(\mathcal{S}_c[\boldsymbol{A_i}]\) captures additional entanglement from quantum fields within the Q-EW. The dominant contribution in \(\mathcal{S}_c[\boldsymbol{A_i}]\) arises from the first linear-law divergent term due to the entanglement between the brane CFT\(_3\) and the bath CFT\(_3\) near the induced AdS\(_3\) boundary, where the spacetime is infinitely stretched. Moreover, the term \(\mathcal{F}_{\boldsymbol{A_i}}\) is the finite, cutoff-independent contribution from the quantum fields deep into the brane within the Q-EW \(\mathcal{W}_{\boldsymbol{A}}\). From the boundary perspective, \(\mathcal{S}_g[\boldsymbol{A_i}]\) arises from the entanglement between \(\boldsymbol{A_i}\) and the remaining boundary \(\partial\mathcal{B}\), while \(\mathcal{S}_c[\boldsymbol{A_i}]\) mainly reflects the volume-law entanglement between \(\boldsymbol{A_i}\) and the bath \(\partial\).

Furthermore, when considering two disjoint intervals \(\boldsymbol{A_1}\) and \(\boldsymbol{A_2}\), the entanglement entropy of the union \(\boldsymbol{A} = \boldsymbol{A_1} \cup \boldsymbol{A_2}\) exhibits a phase transition depending on their separation as \[\begin{align} \label{eq:entropy2} \mathcal{S}[\boldsymbol{A_1} \cup \boldsymbol{A_2}] =\boldsymbol{min}\Bigg\{& c \frac{l_1 + l_2}{\epsilon}+\frac{c'}{3}\left(\log \frac{l_1}{\delta} + \log \frac{l_2}{\delta}\right) + \mathcal{F}_{\boldsymbol{A_1}} + \mathcal{F}_{\boldsymbol{A_2}}, \nonumber \\ & c\frac{l_1 + l_2}{\epsilon}+\frac{c'}{3}\left(\log \frac{l_1 + l_2 + a}{\delta} + \log \frac{a}{\delta}\right) + \mathcal{F}_{\boldsymbol{A_1} \cup \boldsymbol{A_2}} \Bigg\}, \end{align}\tag{14}\] where \(l_1\), \(l_2\), and \(a\) are the lengths of two subregions and their separation, respectively. The terms \(\mathcal{F}_\mathcal{K}\), with \(\mathcal{K} \in \{\boldsymbol{A_1}, \boldsymbol{A_2}, \boldsymbol{A_1} \cup \boldsymbol{A_2}\}\), also represent the finite, cutoff-independent entanglement from the quantum fields inside the wedge \(\mathcal{W}_{\boldsymbol{A_1}\cup\boldsymbol{A_2}}\). Recall that 14 is examined via the numerical analysis in the next section. For sufficiently separated subregions (large \(a\)), the QES of the union \(\boldsymbol{A_1} \cup \boldsymbol{A_2}\) decomposes into two disconnected QES surfaces corresponding to \(\boldsymbol{A_1}\) and \(\boldsymbol{A_2}\), respectively, and the entanglement entropy reduces to the sum of individual contributions as described by (13 ). In contrast, when the subregions are adjacent to each other (small \(a\)), the entropy exhibits a clear deviation from this behavior, signaling a nontrivial phase transition.

Since the UV-divergent part is always addictive as \[\operatorname{Area}(\boldsymbol{A_1}) + \operatorname{Area}(\boldsymbol{A_2}) = \operatorname{Area}(\boldsymbol{A_1}\cup\boldsymbol{A_2}),\] the BMI becomes finite and is governed entirely by the finite contributions. For two nearby subregions, the formula of BMI (7 ) reduces to \[\begin{align} \label{eq:semiMI} \mathcal{I}[\boldsymbol{A_1}:\boldsymbol{A_2}]=&\frac{c'}{3}\log \frac{l_1 l_2}{(l_1 + l_2 + a) a}+\mathcal{F}_{\boldsymbol{A_1} }+\mathcal{F}_{\boldsymbol{A_2}} - \mathcal{F}_{\boldsymbol{A_1} \cup \boldsymbol{A_2}}. \end{align}\tag{15}\] Within this expression, we further have \[\mathcal{I}_g[\boldsymbol{A_1}:\boldsymbol{A_2}]=\mathcal{S}_g[\boldsymbol{A_1}]+\mathcal{S}_g[\boldsymbol{A_2}]-\mathcal{S}_g[\boldsymbol{A_1}\cup\boldsymbol{A_2}]\to \frac{c'}{3}\log \frac{l_1 l_2}{(l_1 + l_2 + a) a},\] and \[\mathcal{I}_c[\boldsymbol{A_1}:\boldsymbol{A_2}]=\mathcal{S}_c[\boldsymbol{A_1}]+\mathcal{S}_c[\boldsymbol{A_2}]-\mathcal{S}_c[\boldsymbol{A_1}\cup\boldsymbol{A_2}]\to\mathcal{F}_{\boldsymbol{A_1} }+\mathcal{F}_{\boldsymbol{A_2}} - \mathcal{F}_{\boldsymbol{A_1} \cup \boldsymbol{A_2}}.\] From the boundary perspective, the first logarithmic term \(\mathcal{I}_g[\boldsymbol{A_1}:\boldsymbol{A_2}]\) denotes the BMI from the CFT\(_2\) on the boundary system \(\boldsymbol{\partial} \mathcal{B}\), while the second finite term \(\mathcal{I}_c[\boldsymbol{A_1}:\boldsymbol{A_2}]\) represents the correction from the CFT\(_3\). From the brane perspective, the first term \(\mathcal{I}_g[\boldsymbol{A_1}:\boldsymbol{A_2}]\) denotes the geometric contributions from the induced gravity on the brane, which equals the area of the QES, whereas the second correction term \(\mathcal{I}_c[\boldsymbol{A_1}:\boldsymbol{A_2}]\) therefore encodes the additional cutoff-independent contributions from the CFT\(_3\) on the brane. Interestingly, \(\mathcal{I}_c[\boldsymbol{A_1}:\boldsymbol{A_2}]\) is generally negative from the later numerical analysis. This is mainly because the Q-EW of the union region \(\boldsymbol{A_1}\cup\boldsymbol{A_2}\) is larger than the sum of those associated with \(\boldsymbol{A_1}\) and \(\boldsymbol{A_2}\) individually. We will elaborate on this behavior in the next section.

3 Numerical analysis↩︎

In this section, we begin by introducing the engineering software – “Surface Evolver”, which is employed to construct the RT surfaces relevant to our analysis. We then assess the numerical accuracy of the results obtained and verify the emergence of universal entanglement behavior. Finally, we explore the properties of BMI in the semiclassical limit.

3.1 Constructions of the RT surfaces↩︎

Traditionally, constructing the RT surface in the bulk involves solving two-dimensional partial differential equations, often facilitated by adopting suitable coordinate systems to regulate divergences near the boundary region \(\boldsymbol{A}\) and by imposing Neumann boundary conditions on the brane [90]. However, this coordinate-based approach faces inherent limitations when investigating BMI, as it becomes challenging to find appropriate coordinate charts for describing disjoint subregions. To address these challenges, we instead employ the shape optimization software “Surface Evolver” [91] to construct the minimal surface \(\Gamma_{\boldsymbol{A}}\). This method has been successfully applied in various asymptotically AdS spacetimes [93]–[97], offering a flexible and robust numerical framework for studying entanglement structures beyond symmetric or connected regions.

In Surface Evolver, any surface is represented as a collection of oriented triangular facets. Given a background metric, a boundary region \(\boldsymbol{A}\) is anchored at \(z = \epsilon\) to regulate divergences near the asymptotic boundary, and an initial trial surface is constructed – Fig.3 (a). The software then evolves this initial configuration toward a local minimum of the area functional via a gradient descent algorithm – see Fig. 3 (b) (see Appendix B of [93] for an overview of Surface Evolver). The output is a triangulated minimal surface \(\Gamma^\epsilon_{\boldsymbol{A}}\), anchored at \(z = \epsilon\), which approximates the true extremal surface \(\Gamma_{\boldsymbol{A}}\) (anchored at \(z = 0\)). The accuracy of this approximation improves with increasing triangulation resolution. For any given triangulated surface, the total area can be directly computed from the mesh data. For simplicity, in what follows we will not distinguish between the numerically evolved surface \(\Gamma^\epsilon_{\boldsymbol{A}}\) and the continuum surface \(\Gamma_{\boldsymbol{A}}\), and will refer to both as \(\Gamma_{\boldsymbol{A}}\).

Since this work primarily focuses on holographic BMI, we consider configurations where two disjoint subregions \(\boldsymbol{A_1}\) and \(\boldsymbol{A_2}\) (with \(\boldsymbol{A} = \boldsymbol{A_1} \cup \boldsymbol{A_2}\)) are both anchored on the boundary \(\boldsymbol{\partial} \mathcal{B}\). In general, the combined region \(\boldsymbol{A}\) exhibits two distinct entanglement phases: When the subregions are sufficiently close, the associated Q-EW \(\mathcal{W}\) is connected – Fig.3 (c). This phase signals strong correlations between \(\boldsymbol{A_1}\) and \(\boldsymbol{A_2}\). In contrast, when the subregions are far apart, the RT surface decomposes into a union of the individual RT surfaces corresponding to \(\boldsymbol{A_1}\) and \(\boldsymbol{A_2}\), each resembling the configuration shown in Fig.3 (b), indicating weak correlations. For visualization purposes, the RT surface shown in Fig. 3 (c) is rendered with a moderate resolution of approximately \(V \simeq 3000\) vertices. However, for precise area calculations, we use a higher-resolution discretization with \(V \simeq 13000\) vertices to ensure sufficient smoothness and numerical accuracy.

In the subsequent subsections, all entanglement-related quantities will be expressed in dimensionless form by normalizing with respect to the central charge \(c=\frac{L^2}{4 G_N^{(4)}}\). Specifically, we define \[\begin{align} \left\{\mathcal{I}[\boldsymbol{A_1}:\boldsymbol{A_2}],\mathcal{S}[\mathcal{\mathcal{K}}]\right\}= \frac{1}{c}\left\{ \mathcal{I}[\boldsymbol{A_1}:\boldsymbol{A_2}],\mathcal{S}[\mathcal{K}]\right\}, \end{align}\] where \(\mathcal{K} \in \{\boldsymbol{A_1},\boldsymbol{A_2},\boldsymbol{A_1}\cup\boldsymbol{A_2}\}\). Accordingly, the geometric and correction terms are also rescaled as \[\left\{\mathcal{I}_g,\mathcal{I}_c,\mathcal{S}_g,\mathcal{S}_c\right\}=\frac{1}{c}\left\{\mathcal{I}_g,\mathcal{I}_c,\mathcal{S}_g,\mathcal{S}_c\right\}.\] With these identifications, the effective central charge \(c'\) always appears in the dimensionless ratio \(c'/c\), which is determined by the dihedral angle via (12 ).

3.2 Entanglement entropy of disjoint subregions↩︎

In this subsection, we consider the entanglement entropy of a disconnected region \(\boldsymbol{A}=\boldsymbol{A_1}\cup\boldsymbol{A_2}\). For simplicity, we restrict our analysis to symmetric configurations with equal subregion lengths \(l_1 = l_2 = l\). Throughout, we set the UV cutoff to be \(\epsilon=0.01\), which is sufficiently small compared to all other relevant physical scales.

For two distant subregions, the Q-EW associated with the union region \(\boldsymbol{A_1}\cup\boldsymbol{A_2}\) becomes disconnected. In this regime, the total entanglement entropy is expected to reduce to the sum of the entropies of the individual components. Thus, we adopt the following fitting function as \[\label{eq:FitEEDis} \frac{1}{2}\mathcal{S}[\boldsymbol{A_1}\cup\boldsymbol{A_2}]=\mathcal{S}[\boldsymbol{A_1}]=\mathcal{S}[\boldsymbol{A_2}]=\frac{l}{\epsilon}+b_a \log \frac{l}{\delta}+b_{a1}.\tag{16}\] Here, \(b_a\) and \(b_{a1}\) are two fitting coefficients. We note that it has been shown that \(b_a\) asymptotically approaches the ratio \(c'/c\) in the semiclassical limit [90]. Therefore, the main aim here is to verify the convergence behavior of our numerical method and assess the robustness of the Surface Evolver-based approach. An illustrative example is presented in Fig. 4 (a), where the total entropy (blue curve) is decomposed into a geometric component (yellow) obeying the logarithmic law, and a correction component (green) satisfying the linear law. These numerical behaviors align well with the structure of the fitting function(16 ). Moreover, Fig. 4 (b) demonstrates that the fitting coefficient \(b_a\) converges to the theoretical ratio \(c'/c\) as \(\theta_0 \to 0\), thereby confirming the expected semiclassical behavior: \[\mathcal{S}_g\to b_a \log \frac{l}{\delta},\quad \text{and}\quad \mathcal{S}_c\to \frac{l}{\epsilon}+\mathcal{F}_{\boldsymbol{A}}, \quad \text{as}\quad \theta_0\to0.\]

For two adjacent subregions, the Q-EW \(\mathcal{W}_{\boldsymbol{A_1}\cup\boldsymbol{A_2}}\) of the combined region is connected. In this case, we propose the following fitting function for the numerical computation of the entanglement entropy as \[\label{eq:FitEECon} \mathcal{S}[\boldsymbol{A_1} \cup \boldsymbol{A_2}] = \frac{2l}{\epsilon} + b_s\left(\log \frac{2l + a}{\delta} + \log \frac{a}{\delta}\right) +b_{s1},\tag{17}\] where \(\{b_s,b_{s1}\}\) are fitting coefficients. The linear-divergent term arises from the entanglement of quantum fields near the induced AdS\(_3\) boundary. The logarithmic-divergent terms originate from the area of the QES – Fig. 3 (c), while other corrections are expressed by the finite term. Within this fitting structure, we identify the geometric and correction contributions as \[\mathcal{S}_g\to b_s\left(\log \frac{2l + a}{\delta} + \log \frac{a}{\delta}\right), \quad \text{and}\quad \mathcal{S}_c\to \frac{2l}{\epsilon}+\mathcal{F}_{\boldsymbol{A_1}\cup\boldsymbol{A_2}}, \quad \text{as}\quad \theta_0\to0.\] As illustrated in Fig. 4 (a), this decomposition clearly separates the divergent contributions. Similarly, in the semiclassical limit \(\theta_0\to0\), the fitting coefficient \(b_s\) converges to the central charge ratio \(c'/c\) – Fig. 5 (b). This observation leads to the conclusion that, regardless of whether the Q-EW is connected or disconnected, the coefficient of the logarithmic divergence consistently approaches \(c'/c\).

As a summary, our numerical analysis confirms the expected behavior of the entanglement entropy for disconnected regions in the semiclassical limit. In general, the total entropy is primarily contributed from two distinct sources: First, linear divergence arises from the quantum fields localized near the boundary of the induced gravitational region, reflecting the volume-law entanglement; Second, logarithmic divergence originates from the area of the QES, and corresponds to the entanglement entropy governed by induced gravity brane. However, there is an additional finite term that generally captures the contributions from quantum fields deep inside the Q-EW. These contributions are sensitive to the global structure of the QES and are not simply additive, depending on whether the Q-EW is connected or disconnected. In the following subsection, we will analyze how these finite-term contributions influence the BMI.

3.3 BMI of the disjoint subregions↩︎

In the previous analysis of the entanglement entropy (13 ) and (14 ), our attention was primarily focused on the divergent structures, while the finite contributions were largely neglected due to their sub-subleading nature and the difficulty in obtaining analytical expressions. However, in the case of BMI, all UV divergences are canceled by construction, leaving behind only finite terms. This feature makes BMI a divergent-free quantity, and simultaneously highlights the importance of reliable numerical methods, as these finite terms are generally inaccessible through purely analytical techniques.

Having already computed the relevant entropies, the BMI in the semiclassical limit can now be directly obtained via simple algebraic combinations, as described in (7 ). In this section, we proceed to numerically evaluate the BMI between two disjoint subregions \(\boldsymbol{A_1}\) and \(\boldsymbol{A_2}\), enabling us to isolate and analyze both the finite geometric contributions and the corrections of quantum fields.

We begin by analyzing the dependence of the BMI \(\mathcal{I}[\boldsymbol{A_1}:\boldsymbol{A_2}]\) on the dimensionless separation \(a/l\) between two subregions, as illustrated in Fig. 6 (a). When the subregions are close, the Q-EW of the union \(\boldsymbol{A_1} \cup \boldsymbol{A_2}\) remains connected, resulting in a non-vanishing (positive) BMI. As the separation increases, the BMI decreases approximately linearly – a behavior that sharply contrasts with predictions from classical gravity. Eventually, the Q-EW becomes disconnected, and the BMI drops to zero, signaling a phase transition in the underlying geometry.

Intriguingly, the transition from the connected phase to disconnected phase occurs at a larger separation \(a/l\) when the dihedral angle \(\theta_0\) decreases – Fig. 6 (b). On one hand, the region above the black curve corresponds to configurations where the Q-EW of \(\boldsymbol{A_1} \cup \boldsymbol{A_2}\) remains connected, indicating strong entanglement between these two subregions. On the other hand, below the curve, Q-EW becomes disconnected, signaling that the entanglement between \(\boldsymbol{A_1}\) and \(\boldsymbol{A_2}\) is weak. This trend suggests that in the semiclassical limit, the entanglement between two subregions is stronger and can persist over larger separations. This observation aligns with the interpretation that decreasing \(\theta_0\) corresponds to an increase in the degrees of freedom (d.o.f.) on the (1+1)-dimensional boundary quantum system \(\partial\mathcal{B}\) [29], [90], [100]. As the ratio \(c'/c\) increases, each subregion, such as \(\boldsymbol{A_1}\), possesses more internal d.o.f. to become entangled with \(\boldsymbol{A_2}\), thereby extending the range of significant correlations and shifting the transition point outward.

Furthermore, as shown in Fig. 7 (a), the geometric contribution \(\mathcal{I}_g[\boldsymbol{A_1}:\boldsymbol{A_2}]\) decreases monotonically with the increase of the separation. Remarkably, this geometric contribution consistently exceeds the total BMI \(\mathcal{I}[\boldsymbol{A_1}:\boldsymbol{A_2}]\), implying that the correction term is always non-positive: \[\begin{align} \mathcal{I}_c[\boldsymbol{A_1}:\boldsymbol{A_2}]=\mathcal{F}_{\boldsymbol{A_1}}+\mathcal{F}_{\boldsymbol{A_2}}-\mathcal{F}_{\boldsymbol{A_1}\cup\boldsymbol{A_2}}\leq0. \end{align}\]

At first glance, this phenomenon may appear counterintuitive, but it admits a natural interpretation from the brane perspective. Recall that each finite term \(\mathcal{F}_{\mathcal{K}}\) represents the cutoff-independent entanglement between the quantum fields contained within the Q-EW \(\mathcal{W}_{\mathcal{K}}\) and the environment. Therefore, these entanglements are naturally expected to be proportional to the area of the corresponding wedge \(\mathcal{W}_{\mathcal{K}}\) (\(\mathcal{K}=\boldsymbol{A_1}, \boldsymbol{A_2}, \boldsymbol{A_1}\cup\boldsymbol{A_2}\)). As illustrated in Fig. 7 (b), when the wedge \(\mathcal{W}_{\boldsymbol{A_1}\cup\boldsymbol{A_2}}\) is connected, its area always exceeds that of the union of two disconnected wedges \(\mathcal{W}_{\boldsymbol{A_1}}\cup\mathcal{W}_{\boldsymbol{A_2}}\). Consequently, the connected wedge \(\mathcal{W}_{\boldsymbol{A_1}\cup\boldsymbol{A_2}}\) contains a greater number of quantum fields. This difference directly leads to a negative value for the correction term in BMI as \[\begin{align} \label{eq:MI95QFT} \mathcal{I}_c[\boldsymbol{A_1}:\boldsymbol{A_2}]=S_{\text{QFT}}(\mathcal{W}_{\boldsymbol{A_1}})+S_{\text{QFT}}(\mathcal{W}_{\boldsymbol{A_2}})-S_{\text{QFT}}(\mathcal{W}_{\boldsymbol{A_1}\cup\boldsymbol{A_2}}). \end{align}\tag{18}\] Furthermore, as the separation between these two subregions increases while their sizes remain fixed (i.e., as \(a/l\) grows), the area of the connected wedge \(\mathcal{W}_{\boldsymbol{A_1} \cup \boldsymbol{A_2}}\) grows increasingly larger relative to the total area of the disconnected wedges \(\mathcal{W}_{\boldsymbol{A_1}}\cup\mathcal{W}_{\boldsymbol{A_2}}\). This growing amplifies the negativity of \(\mathcal{I}_c[\boldsymbol{A_1}:\boldsymbol{A_2}]\), causing the correction term to decrease monotonically with \(a/l\).

It is worth noting the fundamental distinction between the negative correction term \(\mathcal{I}_c\) in doubly holography and the non-negativity of BMI in a standard CFT, as shown in Fig. 8. In the scenario of a standard CFT – Fig. 8 (a), one always has the relation: \[\boldsymbol{Area}(\boldsymbol{A_1})+\boldsymbol{Area}(\boldsymbol{A_2})=\boldsymbol{Area}(\boldsymbol{A_1}\cup\boldsymbol{A_2}),\] which implies that in (7 ) the quantum fields contributing to \(\mathcal{S}[\boldsymbol{A_1}] +\mathcal{S}[\boldsymbol{A_2}]\) coincide with those entering \(\mathcal{S}[\boldsymbol{A_1}\cup \boldsymbol{A_2}]\). Since the former generally captures more entanglement, the BMI is manifestly non-negative. In contrast, within the framework of double holography – Fig. 8 (b), the Q-EWs satisfy the inequality as \[\boldsymbol{Area}(\mathcal{W}_\mathcal{A})+\boldsymbol{Area}(\mathcal{W}_\mathcal{B})\leq\boldsymbol{Area}(\mathcal{W}_{\mathcal{A}\cup\mathcal{B}}).\] This indicates that the quantum fields contributing to the entropy of the individual wedges \(\mathcal{W}\mathcal{A}\) and \(\mathcal{W}\mathcal{B}\) are fewer than those contributing to the combined wedge \(\mathcal{W}_{\mathcal{A} \cup \mathcal{B}}\). As a result, the correction term in the BMI, originating from cutoff-independent finite contributions, is generically non-positive. This discrepancy underlines a key qualitative difference between two holographic setups and highlights the nontrivial structure of entanglement in doubly holographic scenarios.

4 Interpretation from random tensor networks↩︎

From the brane perspective in double holography, the negative contribution to the BMI from bulk quantum fields arises from their volume-law entanglement entropy, which itself is a consequence of strong entanglement between the bulk fields and the radiation.

We reproduce this phenomenon using RTNs. As a toy model of holography, a RTN recovers the generalized-entropy formalism in the large bond-dimension limit and realizes a subsystem quantum error-correcting code with complementary recovery [105], [106].

In this section, we first briefly review the RTN construction and its entropy formula. We then show that a maximally mixed state on the bulk degrees of freedom, whose entropy scales with volume, yields a negative contribution to the BMI between boundary subsystems.

4.1 A review of holographic mapping from RTN↩︎

A RTN state is a random projected entangled pair state (PEPS) with boundary degrees of freedom. Consider a network with boundary, parameterized by \(\left(\mathcal{I},\partial,\mathcal{E}\right)\), where the network consists of a set of internal vertices \(\mathcal{I}\), a set of boundary vertices \(\partial\), and a set of edges \(\mathcal{E}\) connecting vertices. An internal vertex can be connected to more than one other vertex by edges. A boundary vertex is connected to a single internal vertex by one edge.

Although our construction is general, we mostly consider networks embedded in a two-dimensional manifold with negative curvature, as illustrated in Fig. 9, which is the typical case in the context of holography [105].

We associate a Hilbert space \(\mathcal{H}_x\) with each vertex \(x\in\mathcal{I}\cup\partial\) and factorize it into Hilbert spaces \(\mathcal{H}_{xy}\) associated with its neighboring vertices \(y\), together with a local bulk Hilbert space \(\mathcal{H}_{xb}\). Explicitly, \[\mathcal{H}_x=\left(\bigotimes_{y:\,\overline{xy}\in \mathcal{E}} \mathcal{H}_{xy}\right)\otimes \mathcal{H}_{xb},\] where \(y:\,\overline{xy}\in \mathcal{E}\) denotes a neighboring vertex \(y\) connected to \(x\) by the edge \(\overline{xy}\). The dimensions of these Hilbert spaces are denoted by \(d_x=\dim(\mathcal{H}_x)\), \(d_{xb}=\dim(\mathcal{H}_{xb})\), and \(d_{\overline{xy}}=\dim(\mathcal{H}_{xy})=\dim(\mathcal{H}_{yx})\), where we require \(\mathcal{H}_{xy}\) and \(\mathcal{H}_{yx}\) to have the same dimension.

From this factorization, we have \(d_x=d_{xb}\prod_{y:\overline{xy}\in\mathcal{E}}d_{xy}\). For any boundary vertex \(x\in\partial\), we require a trivial local bulk Hilbert space with \(d_{xb}=1\), so that \(\mathcal{H}_x=\mathcal{H}_{xy}\), where \(y\in\mathcal{I}\) is the unique internal vertex connected to \(x\). All local bulk Hilbert spaces together form the bulk Hilbert space \(\mathcal{H}_\mathcal{B}=\bigotimes_{x\in\mathcal{I}}\mathcal{H}_{xb}\), with total dimension \(d_\mathcal{B}=\prod_{x\in\mathcal{I}} d_{xb}\). Similarly, all boundary Hilbert spaces form the boundary Hilbert space \(\mathcal{H}_\partial=\bigotimes_{x\in\partial}\mathcal{H}_x\), with dimension \(d_\partial=\prod_{x\in\partial} d_x\).

We construct a bulk-to-boundary map based on the RTN. On each internal vertex \(x\in\mathcal{I}\), we associate a Haar-random state \(\left|U_x\right\rangle=U_x\left|0_x\right\rangle\in\mathcal{H}_x\), where each \(U_x\) is drawn from the circular unitary ensemble (CUE) on \(\mathcal{H}_x\), and \(\left|0_x\right\rangle\) is a fixed reference state. Due to the factorization of \(\mathcal{H}_x\), the state \(\left|U_x\right\rangle\) defines a random tensor with multiple indices associated with the edge Hilbert spaces \(\mathcal{H}_{xy}\), and one bulk index associated with \(\mathcal{H}_{xb}\).

On each edge \(\overline{xy}\in\mathcal{E}\), we place a maximally entangled state \[\left|\overline{xy}\right\rangle=d_{\overline{xy}}^{-1/2}\sum_{i=1}^{d_{\overline{xy}}} \left|i\right\rangle\otimes\left|i\right\rangle\in \mathcal{H}_{xy}\otimes\mathcal{H}_{yx},\] where \(\{\left|i\right\rangle\}\) are orthonormal bases of \(\mathcal{H}_{xy}\) and \(\mathcal{H}_{yx}\). The RTN is then defined as \[V = \left(\bigotimes_{x\in \mathcal{I}}\sqrt{d_x}\left\langle U_x\right|\right) \left(\bigotimes_{\overline{xy}\in \mathcal{E}}\left|\overline{xy}\right\rangle\right),\] where the factor \(\sqrt{d_x}\) is included to ensure the normalization condition 20 below. This defines a bulk-to-boundary map \(V:\mathcal{H}_\mathcal{B}\to\mathcal{H}_\partial\), such that a bulk state \(\rho_\mathcal{B}\) is mapped to a boundary state \[\begin{align} \label{eq:bdy95state} \rho_\partial=V \rho_\mathcal{B}V^\dagger. \end{align}\tag{19}\]

Denoting the ensemble average over the Haar measure as \[\overline{X}=\int_{\text{Haar}} X\prod_{x\in\mathcal{I}} dU_x,\] with normalization \(\overline{1}=1\), we have \[\overline{\left|U_x\right\rangle\left\langle U_x\right|}=\frac{1}{d_x}\mathbb{I}_x.\] It then follows that \[\begin{align} \label{eq:map95normalization} \overline{V^\dagger V}=\mathbb{I}_\mathcal{B}, \end{align}\tag{20}\] where \(\mathbb{I}_\mathcal{B}\) is the identity operator on \(\mathcal{H}_\mathcal{B}\).

To check whether \(V\) is approximately isometric, one may compute the purity [105] \[\frac{\overline{\mathrm{Tr}[(V^\dagger V)^2]}}{\overline{\mathrm{Tr}[V^\dagger V]^2}} \approx \frac{1}{d_\mathcal{B}} + \frac{1}{d_\partial},\quad d_{xb}\gg1,\;\forall x,\quad d_{\overline{xy}}\gg1,\;\forall\overline{xy}.\] When \(d_\partial\gg d_\mathcal{B}\), the first term dominates and \(V^\dagger V\approx \mathbb{I}_\mathcal{B}\), so that \(V\) becomes an isometric bulk-to-boundary map. Identifying the boundary Hilbert space \(\mathcal{H}_\partial\) with the boundary field theory, the bulk Hilbert space \(\mathcal{H}_\mathcal{B}\) describes weak fluctuations around a classical geometry, and the isometric map \(V\) realizes a one-directional holographic mapping from bulk fluctuations to boundary degrees of freedom [105]³.

The map \(V\) reproduces the Ryu-Takayanagi (RT) formula with bulk-state corrections in the large-dimension limit. For the boundary state 19 , the ensemble-averaged \(n\)-th Rényi entropy of a boundary subregion \(A\) is \[\mathcal{S}_n[A]=\frac{1}{1-n}\ln\frac{\overline{\mathrm{Tr}[\rho_A^n]}}{\overline{\mathrm{Tr}[\rho_A]^n}},\quad \rho_A=\mathrm{Tr}_{\bar A}\rho_\partial,\] where \(\bar A\) denotes the complement of \(A\). The replica partition functions can be mapped to a \(\mathrm{Sym}_n\) spin model and evaluated in the large-\(d_{xb}\) and large-\(d_{\overline{xy}}\) limits. For simplicity, we assume uniform dimensions \(d_{xb}=d_b\) and \(d_{\overline{xy}}=d_e\). The resulting Rényi entropy is [105] \[\begin{align} \label{eq:entropy95RTN} \mathcal{S}_n[A]=\min_{\gamma_A} \left(\left|\gamma_A\right|\ln d_e + S_n[W_A;\rho_\mathcal{B}]\right), \end{align}\tag{21}\] where \(\gamma_A\) is a cut homologous to \(A\), enclosing a bulk region \(W_A\) with \(\partial W_A=A\cup\gamma_A\). This reproduces the quantum extremal surface formula 5 and provides a RTN interpretation of bulk entropy contributions.

4.2 Negative bulk contribution to BMI↩︎

To mimic the double-holographic setup, we consider a bath system \(\mathcal{R}\) with Hilbert space \(\mathcal{H}_\mathcal{R}\) and a pure entangled state \(\left|\psi\right\rangle\in\mathcal{H}_\mathcal{B}\otimes\mathcal{H}_\mathcal{R}\), such that the reduced density matrix on the bulk is \(\rho_\mathcal{B}=\mathrm{Tr}_\mathcal{R}\left|\psi\right\rangle\left\langle\psi\right|\). When the bath \(\mathcal{R}\) is much larger than the bulk system \(\mathcal{B}\), the reduced state \(\rho_\mathcal{B}\) is expected to be highly mixed. The simplest case is \[\rho_\mathcal{B}=\frac{1}{d_\mathcal{B}}\mathbb{I}_\mathcal{B}.\] In this case, the reduced state on any bulk wedge \(W_A\) is also maximally mixed, and the Rényi entropy obeys a volume law, \[\begin{align} \label{eq:volume95law} S_n[W_A;\rho_\mathcal{B}]=\left|W_A\right|\ln d_b, \end{align}\tag{22}\] where \(\left|W_A\right|\) denotes the number of vertices in \(W_A\).

Following Subsec. 2.2, we consider two disconnected boundary regions \(A_1\) and \(A_2\) and compute their Rényi BMI in the large-dimension limit. Using 21 , we find \[\begin{align} \label{eq:MI95RTN} \mathcal{I}_n[A_1:A_2]&=(\left|\gamma_{A_1}\right|+\left|\gamma_{A_2}\right|-\left|\gamma_{A_1\cup A_2}\right|)\ln d_e + \mathcal{I}_{nc}[A_1:A_2], \\ \mathcal{I}_{nc}[A_1:A_2]&= S_n[W_{A_1};\rho_\mathcal{B}] +S_n[W_{A_2};\rho_\mathcal{B}]-S_n[W_{A_1\cup A_2};\rho_\mathcal{B}]. \end{align}\tag{23}\]

Assuming the hierarchy \(d_e\gg d_b\), the minimal cuts are determined purely by geometry, and the inequality \[\begin{align} \label{eq:gamma95ineq95RTN} \left|\gamma_{A_1}\right|+\left|\gamma_{A_2}\right|\geq \left|\gamma_{A_1\cup A_2}\right| \end{align}\tag{24}\] always holds. When \(A_1\) and \(A_2\) are sufficiently close, the minimal surface \(\gamma_{A_1\cup A_2}\) is connected, as shown in Fig. 9. In this case, using the volume law 22 , the bulk contribution to the BMI is non-positive, \[\mathcal{I}_{nc}[A_1:A_2]=(\left|W_{A_1}\right|+\left|W_{A_2}\right|-\left|W_{A_1\cup A_2}\right|)\ln d_b\leq 0,\] with equality in the disconnected phase.

By contrast, if the bulk entropy obeys an area law, \[S_n[W_A;\rho_\mathcal{B}]\sim \left|\gamma_A\right|\ln d_e,\] then 24 implies a non-negative bulk contribution, \[\mathcal{I}_{nc}[A_1:A_2]\sim (\left|\gamma_{A_1}\right|+\left|\gamma_{A_2}\right|-\left|\gamma_{A_1\cup A_2}\right|)\ln d_e\geq 0.\] We therefore conclude that the non-positive contribution of bulk quantum fields to BMI can be understood, within the RTN framework, as a consequence of the bulk being in a highly mixed, volume-law entangled state.

5 Conclusions and discussions↩︎

In this paper, we have investigated the BMI between two one-dimensional spatial subregions \(\boldsymbol{A_1}\) and \(\boldsymbol{A_2}\) on the (1+1)-dimensional boundary system \(\partial\mathcal{B}\), within the framework of double holography. In this setup, the boundary system \(\partial\mathcal{B}\) is coupled to a (2+1)-dimensional heat bath \(\partial\). To compute the BMI, we constructed the corresponding RT surface of a disconnected spatial region \(\boldsymbol{A}=\boldsymbol{A_1}\cup\boldsymbol{A_2}\), using a shape optimization program, “Surface Evolver”. This method begins with an initial oriented surface anchored on a spatial slice of AdS\(_4\) spacetime with a Planck brane \(\mathcal{B}\) and employs a gradient descent algorithm to evolve the surface into an extremal one with minimal area. Notably, Neumann boundary conditions are naturally implemented by constraining the RT surface’s boundary to lie on the brane.

Our analysis is restricted to a specific parameter regime \(\theta_0 \leq \theta_c\), where the Q-EW associated with any finite region \(\boldsymbol{A} \subset \partial\mathcal{B}\) remains finite. For \(\theta_0 > \theta_c\), by contrast, the size of the wedge collapses to zero, despite \(\boldsymbol{A}\) remaining finite [90].

We first validated the convergence and robustness of the numerical method by computing the entanglement entropy of a single spatial region on \(\partial\mathcal{B}\). The resulting entropy exhibits both a leading linear divergence and a subleading logarithmic divergence. The leading divergence originates from the brane CFT\(_3\) near the conformal boundary, while the logarithmic term arises from the geometric contribution of the area of the QES. The coefficient of this logarithmic term matches the central charge \(c'\) of the CFT\(_2\) on \(\partial\mathcal{B}\), consistent with the results in literature [90]. We then extended this analysis to the more technically involved case of a disconnected region \(\boldsymbol{A} = \boldsymbol{A_1} \cup \boldsymbol{A_2}\), finding a similar entropy structure in (14 ), again confirming the convergence of the coefficient to \(c'\) in the semiclassical limit.

The core of this study has focused on characterizing the entanglement structure between two codimension-three subregions by computing their BMI \(\mathcal{I}[\boldsymbol{A_1} : \boldsymbol{A_2}]\). In the semiclassical limit, the BMI can be decomposed into a geometric contribution \(\mathcal{I}_g[\boldsymbol{A_1} : \boldsymbol{A_2}]\) from the areas of the QESs, which is finite and logarithmic, and a correction term \(\mathcal{I}_c[\boldsymbol{A_1} : \boldsymbol{A_2}]\), stemming from quantum fields within the Q-EW. Notably, we have found that \(\mathcal{I}_g \geq \mathcal{I}\) and \(\mathcal{I}c \leq 0\). This negative correction arises because, although the UV components of the quantum fields within \(\mathcal{W}_{\boldsymbol{A_1}}\), \(\mathcal{W}_{\boldsymbol{A_2}}\), and \(\mathcal{W}_{\boldsymbol{A_1} \cup \boldsymbol{A_2}}\) are similar, the IR contributions differ significantly: the connected wedge \(\mathcal{W}_{\boldsymbol{A_1} \cup \boldsymbol{A_2}}\) contains obvious more quantum fields than the union of the individual wedges, leading to \(\mathcal{I}_c < 0\) as described by Eq.(9 ) and illustrated in Fig.8 (b).

Earlier work [90] conjectured the existence of an intermediate phase, in which the Q-EW on the brane is disconnected but the classical extremal surface in the bulk remains connected. However, our numerical analysis shows that such configurations are unstable: after a finite number of mesh refinements, they invariably converge to one of two stable phases discussed in this work. Therefore, these configurations cannot be regarded as genuine extremal surfaces. The underlying reason for the absence of this intermediate phase deserves further investigation.

Within the RTN framework, the negative contribution of bulk quantum fields to BMI acquires a transparent information-theoretic interpretation. In the large bond-dimension limit, the RTN entropy formula separates naturally into a geometric term determined by minimal cuts and a bulk entropy term associated with degrees of freedom inside Q-EWs. When the bulk is maximally mixed due to strong entanglement with a heat bath, its entropy obeys a volume law. As a consequence, the contribution from bulk fields to BMI becomes non-positive whenever the combined wedge contains strictly more degrees of freedom of fields than the individual wedges.

This analysis clarifies the physical origin of the negative correction term observed in double holography. From the brane perspective, fields on the brane are strongly entangled with the heat bath, leading to highly mixed states inside Q-EWs. The RTN model demonstrates that such volume-law entanglement generically reduces BMI by introducing a negative bulk entropy contribution. By contrast, if the bulk entropy followed an area law, the same RTN construction would yield a non-negative correction, emphasizing that the sign of the bulk term is controlled by the entanglement structure.

Looking forward, it would be of great interest to study other quantum information measures in the context of double holography. In general, any such measure may be decomposed into a geometric component from the induced gravity on the brane and a correction component from the quantum fields on the brane. The geometric contribution can be captured by the QES and reduces to classical behavior in the semiclassical limit. The correction component, involving quantum fields deep into the brane, is generally non-analytic and requires further numerical investigation.

Finally, an important extension would be to explore finite-temperature systems by introducing black holes into the bulk geometry. This involves solving the backreacted brane geometry, possibly via the DeTurck method [108], [109]. As is well-known, raising the temperature of the system will disrupt the long-range correlations. Consequently, both the QES and the quantum field contributions would be expected to undergo qualitative changes in their behavior, making this direction particularly compelling for future work.

Acknowledgments↩︎

We are grateful to Alexander Jahn and Yuan Sunfor the helpful discussions. Liu Yuxuan special thanks to Peiwen Cao for supporting his work. LYX is supported by the Natural Science Foundation of China under Grant No. 12405079, the Natural Science Foundation of Hunan Province, China (Grant No. 2025JJ60062), and Research start-up funds from the Central South University. YL is supported by the Natural Science Foundation of China under Grant No. 12275275. ZYX also acknowledges support from the berlin Quantum Initiative.

References↩︎

---

1. See related discussions in [19], [26]–[89].↩︎

2. Alternative boundary conditions have been discussed in [95].↩︎

3. A boundary-to-bulk isometry can be realized by enlarging the bulk Hilbert space [107].↩︎