Holographic Entanglement Entropy of Local Quenches in AdS$_4$/CFT$_3$: A Finite-Element Approach

TL;DR

This paper develops a finite-element numerical method to compute holographic entanglement entropy in higher-dimensional CFTs, revealing early-time agreement with the first law and complex late-time behavior for local quenches in AdS$_4$/CFT$_3$.

Contribution

It introduces a finite-element approach for calculating holographic entanglement entropy in higher dimensions and explores its dynamics after local quenches, highlighting deviations from the first law.

Findings

Early-time entanglement growth matches the first law for small conformal dimension.

Deviations from the first law occur at later times, showing sub-linear growth.

Large conformal dimension leads to qualitatively different time dependence than in 2D cases.

Abstract

Understanding quantum entanglement in interacting higher-dimensional conformal field theories is a challenging task, as direct analytical calculations are often impossible to perform. With holographic entanglement entropy, calculations of entanglement entropy turn into a problem of finding extremal surfaces in a curved spacetime, which we tackle with a numerical finite-element approach. In this paper, we compute the entanglement entropy between two half-spaces resulting from a local quench, triggered by a local operator insertion in a CFT$_3$. We find that the growth of entanglement entropy at early time agrees with the prediction from the first law, as long as the conformal dimension $\Delta$ of the local operator is small. Within the limited time region that we can probe numerically, we observe deviations from the first law and a transition to sub-linear growth at later time. In particular, the time dependence at large $\Delta$ shows qualitative differences to the simple logarithmic time dependence familiar from the CFT$_2$ case. We hope that our work will motivate further studies, both numerical and analytical, on entanglement entropy in higher dimensions.

Equations

$$|\Psi ⟩=\mathcal{N}{e}^{-\alpha H}O(x)|0⟩,$$

$$\text{d}{s}^2=\frac{R^2}{z^2}\left(\text{d}{z}^2-\text{d}{t}^2+\sum _{i=1}^{d-1}\text{d}{x}_i^2\right)\text{ .}$$

$$z(t)=\sqrt{t^2+{\alpha }^2}\text{ ,}{x}_1(t)=\cdots =x_{d-1}(t)=0\text{ .}$$

$$E=\frac{R}{z(t=0)}m=\frac{R}{\alpha }m\text{ .}$$

$$\text{d}{s}^2=-\left(R^2+r^2\right)\text{d}{\tau }^2+\frac{R^2\text{d}{r}^2}{R^2+r^2}+r^2\text{d}{\Omega }^2\text{ .}$$

$$T$$

$$W$$

$$Z$$

$$X_i$$

$$\text{d}{s}^2=-(R^2+r^2-\frac{M}{r^{d-2}})\text{d}{\tau }^2+\frac{R^2}{R^2+r^2-\frac{M}{r^{d-2}}}\text{d}{r}^2+r^2\text{d}{\Omega }_{d-1}^2\text{ ,}$$

$$M=\frac{8\Gamma (\frac{d}{2})G_N{R}^2}{(d-1){\pi }^{d/2-1}}m\text{ ,}$$

$$T_{tt}(\rho ,t)=\frac{M{\alpha }^3}{\pi G_{N}R}\frac{({\alpha }^2+t^2+{\rho }^2{)}^2+2t^2{\rho }^2}{{\left(({\alpha }^2+t^2-{\rho }^2{)}^2+4{\alpha }^2{\rho }^2\right)}^{5/2}}\text{ .}$$

$$S_A=\frac{|{\gamma }_A|}{4G_N}\text{ ,}$$

$$A=\left\{(t=t_0,z=z_0,x_1,\dots ,x_{d-1})|x_1\in {\mathbb{R}}_{>0}\wedge x_2,\dots ,x_{d-1}\in \mathbb{R}\right\}\text{ ,}$$

$$\tan {\tau }_i$$

$$\tan {\theta }_i$$

$$r_i$$

$$S_A=\frac{R}{4G_N}\left[\log \left(\frac{2}{s}{r}_1{r}_2\right)+\log \left(\cos \left(\frac{\sqrt{s}}{R}|{\tau }_2-{\tau }_1|\right)-\cos \left(\frac{\sqrt{s}}{R}|{\theta }_2-{\theta }_1|\right)\right)\right]\text{ ,}$$

$$\Delta S_A(t=0)=\frac{R}{2G_N}\log \left(\frac{\sin \left(\frac{\pi }{2}\sqrt{1-M/R^2}\right)}{\sqrt{1-M/R^2}}\right)\text{ .}$$

$$\Delta S_A(t)$$

$$\simeq \frac{R}{4G_N}\left[\log \left(\frac{\sin \left(\pi \sqrt{1-M/R^2}\right)}{\sqrt{1-M/R^2}}\right)+\log \frac{t}{\alpha }\right]\text{ .}$$

$$\text{d}{s}^2=-(R^2+r^2-\frac{M}{r})\text{d}{\tau }^2+\frac{R^2}{R^2+r^2-\frac{M}{r}}\text{d}{r}^2+r^2\text{d}{\theta }^2+r^2{\sin }^2\theta \text{d}{\varphi }^2\text{ .}$$

$$\tan \tau$$

$$r\cos \varphi \sin \theta$$

$$r\sin \varphi \sin \theta$$

$$r\cos \theta$$

$${\gamma }_A^{(0)}=\left\{(t=t_0,z,x,y=0)|x,z\in \mathbb{R}\wedge z>z_0\right\}\text{ .}$$

$$\Delta S_A=\frac{|{\gamma }_A|-|{\gamma }_A^{(0)}|}{4G_N}$$

$$|{\gamma }_A|$$

$$=\underset{|{\gamma }_A^{(0)}|}{\underbrace{\int {\text{d}}^2\xi \sqrt{\det G^{(0)}}}}+\underset{|{\gamma }_A^{(1)}|}{\underbrace{\frac{1}{2}\int {\text{d}}^2\xi \sqrt{\det G^{(0)}}\text{Tr}[G^{(1)}(G^{(0)}{)}^{-1}]}}+O(M^2)\text{ ,}$$

$$\Delta S_A\simeq \frac{|{\gamma }_A^{(1)}|}{4G_N}\text{ .}$$

$$|{\gamma }_A^{(1)}|=\frac{1}{2}{\int }_0^{\infty }\text{d}z\int {\text{d}}^{2}x\frac{R}{z}\mathrm{Tr}{G}^{(1)}\simeq \frac{4Mt}{R\alpha }\text{ .}$$

$$a{X}^2+b{Z}^2=r^2(a{\sin }^2\theta {\cos }^2\varphi +b{\cos }^2\theta )=\text{const ;}Y=r\sin \theta \sin \varphi =0\text{ ,}$$

$$\Delta S_A=\frac{R^2}{4G_N}\cdot \Delta \tilde{A}\left(\tilde{t},\tilde{M}\right)\text{ ,}$$

$$\tilde{M}=\frac{\Delta }{c_{3d}}\text{ .}$$

Full text

1. 1 Introduction

2. 2 Entanglement Entropy of Local Quenches

3. 2.1 Local quenches in AdS/CFT

4. 2.2 Half-space entanglement entropy

5. 2.3 Extension to AdS4/CFT3

6. 3 Numerical Studies in AdS4/CFT3

7. 3.1 Numerical surface extremization

8. 3.2 Computation of half-space entanglement entropy

9. 3.3 Numerical results

10. 4 Local Quenches and the First Law of Entanglement Entropy

11. 4.1 AdS3/CFT2

12. 4.2 AdS4/CFT3

13. 5 Summary and Discussion

14. A Finite Element Implementation

15. A.1 Area calculation

16. A.2 Triangulation and quadrilateralization

17. A.3 Numerical optimization

18. A.4 Normals

19. A.5 Time-slice constraints

20. A.6 Accuracy and error estimation

21. A.7 Performance

1 Introduction

The discovery of the AdS/CFT correspondence Ma ; Gubser:1998bc_Witten:1998qj has precipitated a number of new research directions in theoretical physics. One of the fields that has greatly benefitted from AdS/CFT involves the study of entanglement entropy. In particular, AdS/CFT provides a geometrical method to compute entanglement entropy Ryu:2006bv ; HRT (see CHM ; LM for its derivations). Although this approach had originally been developed from ideas related to the Bekenstein-Hawking formula for black hole entropy BKLS ; Sr , it has found applications within condensed matter physics and quantum information theory. Holographic entanglement entropy provides a useful method for certain strongly coupled quantum systems, so-called holographic CFTs, which are dual to classical gravity via AdS/CFT. The holographic approach has an advantage especially in dimensions larger than two111 Throughout this paper, we generally work in a relativistic setting, with dimensions refering to space-time if not otherwise noted. We also use natural units with $c=1$., as the analysis of entanglement entropy in interacting higher-dimensional CFTs is quite difficult.

The non-equilibrium dynamics of entanglement entropy are a field of extensive research (for a general review, see Eisert:2008ur ; Eisert:2014jea ). There has already been much progress on homogeneous excitations such as global quantum quenches, which can be analytically studied in 1+1-dimensional CFTs GQ . Holographic studies of entanglement entropy under global quenches AAL ; Ba ; HaMa ; Kundu:2016cgh have been successful even in higher dimensions. However, for local excitations in CFT, our knowledge of the behavior of entanglement entropy is highly limited, especially in higher dimensions. In this paper, we are therefore interested in holographic entanglement entropy for a specific class of locally excited states, referred to as local quenches222 Note that there is another class of local quenches where CFTs on two semi-infinite systems are instantaneously joined together Calabrese:2007 ; Calabrese:2016xau . We will not discuss this class of local quenches in this paper., in 2+1-dimensional CFTs using AdS${}_{4}/$CFT3. In the CFT description, we are considering excited states $|\Psi\rangle$ defined by acting with a local operator $O(x)$ on the CFT vacuum $|0\rangle$ in the manner

[TABLE]

where ${\cal N}$ is the normalization factor to unit norm. The parameter $\alpha>0$ provides a UV regularization as the literally point-like localized operator has infinite energy and is singular. We are interested in the time evolution of the entanglement entropy $S_{A}$ for the excited state $|\Psi(t)\rangle=e^{-iHt}|\Psi\rangle$ when we choose the subsystem $A$ to be the half-space. The excitation is located on the boundary between both half-spaces, thus producing additional entanglement between them. As time increases, a larger causal region is affected by the quench. Our main quantity of interest is the resulting growth of entanglement entropy $\Delta S_{A}(t)=S_{A}(|\Psi(t)\rangle)-S_{A}(|0\rangle)$ compared to the vacuum.

Previous analyses of $\Delta S_{A}$ for massless scalar fields have been performed in Nozaki:2014hna ; Nozaki:2014uaa ; Nozaki:2015mca ; Nozaki:2016mcy and it was found that the growth $\Delta S_{A}(t)$ approaches a finite positive constant at late time. This is interpreted as a system of entangled particles propagating at the speed of light (see the recent discussion Nozaki:2017hby ). The same behavior has been found for rational CFTs in two dimensions Nozaki:2014uaa ; He:2014mwa ; Chen:2015usa ; Caputa:2015tua ; Caputa:2016yzn ; Numasawa:2016kmo . Furthermore, a recent study of 1+1-dimensional orbifold CFTs found an exotic time evolution $\Delta S_{A}\propto\log\log t$ for irrational CFTs Caputa:2017tju . For other field theoretic progress on local quenches refer also to Shiba:2014uia ; Caputa:2014eta ; deBoer:2014sna ; Guo:2015uwa ; Caputa:2015waa ; Rangamani:2015agy ; David:2016pzn ; Sivaramakrishnan:2016qdv ; Numasawa:2016emc .

However, holographic results have so far been limited to the AdS${}_{3}/$CFT2 setup, where we can analytically compute $\Delta S_{A}$ Nozaki:2013wia ; Caputa:2014vaa . In this holographic description, the local excitation corresponds to a massive particle falling in AdS3, whose mass $m$ is related to the conformal dimension $\Delta$ of the local operator $O(x)$ in (1) via the standard relation $\Delta\simeq mR$, with $R$ being the AdS radius. The holographic results for 1+1-dimensional CFT show that $\Delta S_{A}\propto\log t$ under time evolution at late time Nozaki:2013wia ; Caputa:2014vaa . This time dependence has been precisely reproduced in Asplund:2014coa using a large $c$ CFT analysis. Such a behavior is assumed to stem from the chaotic nature of holographic CFTs, where the quasi-particle picture breaks down.

The main purpose of this paper is to conduct analogous holographic computations for AdS${}_{4}/$CFT3 (see also Rangamani:2015agy for different perspectives on this problem). Unfortunately, perturbative results on the AdS side are not useful for several interesting cases when we consider entanglement between two large regions. Thus, we would like to obtain an exact result for such limits and compare it with the lower-dimensional counterpart.

Our numerical approach relies on a finite-element optimization strategy which approximates extremal surfaces, required for the calculation of holographic entanglement entropy, by a discrete mesh. In contrast to previous studies using the finite-element method Fonda:2014cca , we require a method that does not restrict the extremal surface to a timeslice (when it becomes equivalent to a minimal surface), but can find more complicated space-like solutions. This is neccessary to tackle problems without timeslice constraints, such as the holographic local quench model, where no useful Killing symmetry of the time-dependent gravity background exists.

This paper consists of three parts. In section 2, we will review the holographic local quench model and its solutions in $\text{AdS}_{3}/\text{CFT}_{2}$. Section 3 will present the application of our numerical approach to the $\text{AdS}_{4}/\text{CFT}_{3}$ case and the computational results obtained with it. An interpretation of these results in terms of more general principles of entanglement entropy will be given in section 4. The details of our numerical method will be described in appendix A.

2 Entanglement Entropy of Local Quenches

2.1 Local quenches in AdS/CFT

A local excitation in a CFT can be described by a holographic dual consisting of a freely falling mass in AdS Poincaré coordinates Nozaki:2013wia , based on the construction Horowitz:1999gf . Pure (i.e. “empty”) AdSd+1 spacetime in Poincaré coordinates $(t,z,x_{1},\dots,x_{d-1})$ corresponds to

[TABLE]

In this notation, $R$ is the AdS radius. In such a geometry, a falling mass $m$ which is at rest at time $t=0$ at position $(z,x_{1},\dots,x_{d-1})=(\alpha,0,\dots,0)$ follows a trajectory

[TABLE]

Such a falling mass has a conserved energy $E$, which can be most easily evaluated as the rest energy at $t=0$:

[TABLE]

Note that $z(t)/R$ is equivalent to an energy scale, or inverse length scale, of pure AdS spacetime.

The induced metric following from the insertion of a mass into AdSd+1 spacetime is quite complicated. Fortunately, as found in Horowitz:1999gf , we can express it more conveniently by switching to global coordinates $(\tau,r,\theta_{1},\dots,\theta_{d-1})$ with a new time coordinate $\tau$, radius $r$ and $d-1$ angular coordinates $\theta_{i}$. For the pure AdS case, this corresponds to

[TABLE]

Here, $\text{d}\Omega$ is the standard angular differential on $S^{d-1}$, e.g. $\text{d}\Omega^{2}=\text{d}\theta_{1}^{2}+\sin^{2}(\theta_{1})\text{d}\theta_{2}^{2}$ for $d=3$. The new time coordinate $\tau$ is unitless. The general coordinate transformation between Poincaré and global coordinates is given by

[TABLE]

where we expressed the transformation in terms of the usual AdS embedding coordinates $(T,W,Z,X_{1},\dots,X_{d-1})$, in a $\mathbb{R}^{d,2}$ manifold. The expressions $\Omega_{1},\dots,\Omega_{d}$ are the Cartesian components of the unit vector in terms of the angular coordinates, e.g. $\Omega_{1}=\cos\theta_{1}\sin\theta_{2}$ for $d=2$. Note that the AdS boundary, located at $z=0$ in Poincaré coordinates, appears at $r=\infty$ in global coordinates.

The real parameter $\beta$, corresponding to an additional boost transformation, does not appear in the invariant $\text{d}s^{2}$ in the pure AdS case. This changes when we consider the additional mass $m$: If we choose $e^{\beta}={\alpha\over R}$, the trajectory (3) corresponds to a static point at $r=0$ in global coordinates. Thus, the induced metric in global coordinates is that of a static AdS black hole, independent of global time $\tau$. Explicitly, the full solution Witten:1998zw is given by

[TABLE]

where the parameter $M$ is related to the mass $m$ via

[TABLE]

where $G_{N}$ is the Newton constant (gravitational constant). We can use the coordinate transformation (6) to map the solution (7) back to Poincaré coordinates $(t,z,x_{1},\dots,x_{d-1})$. However, the resulting expression for the invariant $\text{d}s^{2}$ is rather involved and is therefore omitted here.

The falling mass in Poincaré AdSd+1 corresponds to a local excitation in the CFTd on the boundary. This is visualized in figure 1 for the case $d=3$: We show the component $g_{tt}$ of the induced metric in Poincaré coordinates as well as the corresponding boundary energy density, which is given by the $T_{tt}$ component of the stress energy tensor at $z\to 0$ (for details on the procedure, see Nozaki:2013wia ). Introducing the radial coordinate $\rho=\sqrt{x_{1}^{2}+x_{2}^{2}}$, the latter is given by

[TABLE]

As we can see, the falling mass in the AdS4 bulk creates complicated deviations from the pure AdS metric, which correspond to a radially symmetric excitation on the boundary starting at $t=0$. The peak of this shockwave has a width $\sim\alpha$ and approaches the speed of light at $t\gg\alpha$. This falling particle description can be applied to local quenches in dimensions $d$ larger than three, as well Nozaki:2013wia ; Nozaki:2014uaa .

The falling particle corresponds to a localized state in the AdS bulk, so we expect entanglement between any two regions on the peak of the shockwave. Thus, it is compelling to take a closer look at the entanglement properties of the local quench system.

2.2 Half-space entanglement entropy

The entanglement entropy $S_{A}$ of a subsystem $A$ of a CFT quantifies the entanglement between $A$ and its conjugate $B=A^{C}$ of the CFT space. By using the holographic entanglement entropy formula Ryu:2006bv ; HRT , it can be calculated as

[TABLE]

where $\gamma_{A}$ is an extremal hypersurface along the boundary of $A$ (i.e. $\partial\gamma_{A}=\partial A$) extending into the AdS bulk, and $|\gamma_{A}|$ is its area. For the AdSd+1/CFTd setup, $\gamma_{A}$ has $d-1$ dimensions.

For the case of a local quench, an interesting choice of $A$ is that of a half-space in Poincaré coordinates, given by

[TABLE]

where we consider the system at constant time $t_{0}$ and UV cutoff $z_{0}$. For $d=2$, the setup resembles the joining of two half-line subsystems at $x_{1}=0$, with $S_{A}$ characterizing the entanglement between both subsystems after they have been connected. For $d>2$, it corresponds to the entanglement between two half-spaces after a point-like excitation on the boundary between them.

In the $d=2$ case, $S_{A}$ can be analytically calculated using (10). Consider a compact region $A^{\prime}$ bounded by the points $(t_{0},z_{0},x_{1})$ and $(t_{0},z_{0},x_{2})$ in Poincaré coordinates $(t,z,x)$. Clearly, $A^{\prime}\to A$ as $(x_{1},x_{2})\to(0,\infty)$. In global coordinates, the boundary points correspond to:

[TABLE]

with $i=1,2$. Note that the constant $\alpha$ acts as a scaling factor in Poincaré coordinates, while $R$ corresponds to a radial rescaling in global coordinates. Also, $r_{1}$ and $r_{2}$ diverge as the points approach the AdS boundary at $z_{0}=0$.

The entanglement entropy for this subsystem $A$ is given by Nozaki:2013wia

[TABLE]

with $s=R^{2}-M$. For $s<0$, the expression can be analytically continued so that $S_{A}$ is real. We have omitted terms of order $O\left(r_{1}^{-2},r_{2}^{-2}\right)$ that vanish on the AdS boundary at $r_{1},r_{2}\to\infty$. In this limit, the first term is logarithmically divergent. Therefore, we rely on the finite quantity $\Delta S_{A}=S_{A}-\left.S_{A}\right|_{M=0}$ to describe the system, giving the entanglement entropy excited by the local quench.

The half-line subsystem we are interested in corresponds to setting $x_{1}=0$ and considering the limit $x_{2}\to\infty$. First, consider the $t=0$ case, where the entanglement entropy becomes

[TABLE]

As we are mainly interested in the region $t\gg\alpha$ where the two peaks of the CFT excitation are clearly separated, we can also expand 13 in ${\alpha\over t}$, yielding

[TABLE]

To produce the second line, we assumed $\sqrt{1-M/R^{2}}\cot\left(\pi\sqrt{1-M/R^{2}}\right)\ll{t\over\alpha}$. As $M=8G_{N}R^{2}m$ is a nonzero parameter, this simplification is valid at large time $t$. Note that the central charge is given by the standard formula $c={3R\over 2G_{N}}$ BrHe and the conformal dimension of the local operator is expressed as $\Delta={M\over 8G_{N}R}$. This leads to the relation $M/R^{2}=24\Delta/c$.

Thus, a local quench produces an asymptotically logarithmic time dependence in holographic 2-dimensional CFTs Nozaki:2013wia ; Caputa:2014vaa . Furthermore, the coefficient of logarithmic growth ${R\over 4G_{N}}={c\over 6}$ does not depend on $M$ (or equally the conformal dimension $\Delta$), that is, the size of the excitation. This logarithmic behavior was reproduced precisely in the 2-dimensional CFT analysis of Asplund:2014coa , using the large central charge method. It is curious to note that a similar logarithmic behavior has also been found in CFT calculations for a different class of local quenches created by joining two half-lines Calabrese:2007 , though the coefficient is different by a factor of two. Note also that the dependence on time $t$ in the form ${t\over\alpha}$ follows directly from the invariance of the system under $(t,z,x)\to(a\,t,a\,z,a\,x)$ and $\alpha\to a\,\alpha$ for any $a>0$.

It might be intriguing to compare the above results for holographic CFTs with those for integrable CFTs, where field theoretic computations can be done analytically. For 2-dimensional rational CFTs, including free CFTs, the growth of entanglement $\Delta S_{A}$ is finite and becomes a step function in the limit $\alpha\to 0$, which is simply explained by the behaviour of entangled particles propagating at the speed of light Nozaki:2014hna ; Nozaki:2014uaa ; He:2014mwa ; Chen:2015usa ; Caputa:2015tua ; Caputa:2016yzn ; Numasawa:2016kmo . Recently, another class of integrable 2-dimensional CFTs, namely orbifold CFTs, has been studied Caputa:2017tju . For irrational CFTs, an exotic time evolution $\Delta S_{A}\propto\log\log t$ has been found. Refer to Shiba:2014uia ; Caputa:2014eta ; deBoer:2014sna ; Guo:2015uwa ; Caputa:2015waa ; David:2016pzn ; Numasawa:2016emc for further results for local quenches in two dimensions.

In 3-dimensional CFTs, on the other hand, the only available results for local quenches created by the local operator insertions are those for the free field CFTs Nozaki:2014hna ; Nozaki:2014uaa ; Nozaki:2015mca ; Nozaki:2016mcy . In free field CFTs, we can again understand the evolution of entanglement entropy based on the picture of entangled particles propagating at the speed of light Nozaki:2014hna ; Nozaki:2014uaa ; Nozaki:2017hby . Therefore, it is desirable to explore many other examples of interacting CFTs in higher dimensions. Motivated by this, we focus on holographic CFTs in three dimensions (CFT3) in this paper. It is natural to ask whether our holographic 3-dimensional case exhibits the same logarithmic growth $\Delta S_{A}\propto\log t$ as the 2-dimensional one (2.2).

2.3 Extension to AdS4/CFT3

We will now consider the $d=3$ case of a holographic local quench. Again, we map the Poincaré coordinate description of a falling mass to global coordinates, where the induced metric (7) takes the form:

[TABLE]

According to (8), the mass parameter $M$ is now related to the real mass $m$ via $M=2G_{N}R^{2}m$. A point in global coordinates $(\tau,r,\phi,\theta)$ is translated to a point in Poincaré coordinates $(t,z,x,y)$ via

[TABLE]

The half-space subsystem (11) we are particularly interested in corresponds to a boundary $\partial A$ that is given by a line at $y=0$ in Poincaré coordinates. In the pure AdS case, i.e. for $M=0$, the extremal surface $\gamma_{A}$ is simply the $xz$ plane at constant time:

[TABLE]

Even for a finite $z_{0}$ the surface area $|\gamma_{A}^{(0)}|$ is infinite. For any nonzero $M$, however, the corresponding extremal surface $\gamma_{A}$ will only significantly differ from $\gamma_{A}^{(0)}$ in some local region around the mass. Thus, the growth of entanglement entropy

[TABLE]

is well-defined and finite.

Assuming only small deformations from $\gamma_{A}^{(0)}$ to $\gamma_{A}$, we can attempt to calculate $\Delta S_{A}$ perturbatively. The induced metric $G_{\mu\nu}$ resulting from projecting the full metric for $M>0$ on the plane $\gamma_{A}^{(0)}$ can be expanded in orders of $M$. The area of $\gamma_{A}$ can then be approximated by the expansion

[TABLE]

where $\xi$ denotes coordinates parameterizing $\gamma_{A}^{(0)}$. The first term of this expansion becomes irrelevant, as we can now write

[TABLE]

Evaluating this expression for the full $M>0$ metric on the $xz$ plane and considering the limit $t\ll\alpha$, we find that

[TABLE]

This perturbative approach suggests a linear growth of the entanglement entropy with time $t$ for small perturbations in the minimal surface. However, this result rests on the assumption that at small $M$, the growth of entanglement entropy is dominated less by changes in the shape of $\gamma_{A}$ than by changes in the background metric. Indeed this approximation corresponds to the first law relation of entanglement entropy, which can be applied only for small excitation, as we will explain in section 4.

Thus, we desire a non-perturbative approach to calculating the extremal surface area $|\gamma_{A}|$. However, an analytic computation holds considerable challenges. Even for the simple case of a disk-shaped subsystem $A$ around $(x,y)=(0,0)$, where (17) allows us to constrain the solution to constant $\tau$, computing $\gamma_{A}$ requires the evaluation of integrals that do not have an analytic expression. Solutions for the half-space subsystem are even more involved, motivatating the use of a numerical approach. Note that this situation differs from the AdS3 case, where we can compute the holographic entanglement entropy analytically. Similarly, analytical results can be obtained from field theoretic computations only in two dimensions. Therefore a non-perturbative AdS4 analysis using numerical tools allows us to make predictions which are impossible in any field theoretic analysis currently available.

3 Numerical Studies in AdS4/CFT3

3.1 Numerical surface extremization

Finding an optimal surface is a common numerical problem, often tackled using a finite element discretization. This discretization, usually a triangulation or quadrilateralization, approximates a continuous surface by a finite number of parameters that can be varied until a solution is found that optimizes a function of these parameters. By recursively refining the discretization, i.e. enlarging the parameter space of the optimization, a series of discretized approximations converging to the continuous solution is produced.

When searching for an extremal surface, the optimization function is the area of the surface itself. Such optimization methods are frequently used, and have even been applied to minimal surfaces in pure AdS4 spacetime Fonda:2014cca , which are related to ground state entanglement entropies. However, in our case of time-dependent backgrounds, we need to find extremal surfaces, i.e. space-like surfaces minimal with regard to space-like variations and maximal with regard to time-like ones. This introduces considerable complications compared to simple minimization problems. In particular, it imposes highly nonlinear constraints on the parameter space of discretized solutions, as each discretization element has to remain purely space-like.

The details of our implementation are shown in appendix A. In principle, it can be used to calculate extremal surfaces for a wide range of boundary conditions and 4-dimensional metrics.

3.2 Computation of half-space entanglement entropy

We calculate the growth of entanglement entropy $\Delta S_{A}=(|\gamma_{A}|-|\gamma_{A}^{(0)}|)/4G_{N}$ produced by a falling mass in the bulk AdS4 spacetime. As described in the previous section, the extremal surface $\gamma_{A}$ is computed for a given $M$, while $\gamma_{A}^{(0)}$ is the flat extremal surface for $M=0$. The numerical computation is performed in global coordinates $(\tau,r,\phi,\theta)$333 As explained in A.5, we actually use a modified global time coordinate $\tau^{\prime}$ which is constant on the boundary. with a radial coordinate $r$ and a static spherical horizon at $r=r_{\text{hor}}(M)$. $\gamma_{A}^{(0)}$ is equal to the $XZ$ plane, where $X=r\sin\theta\cos\phi$ and $Z=r\cos\theta$. Because the areas of both surfaces are divergent, a cutoff parameter is required. Notice that $\gamma_{A}$ turns into $\gamma_{A}^{(0)}$ at large $r$, far away from the horizon. Thus, we should choose our cutoff along $\gamma_{A}^{(0)}$, i.e. on the $XZ$ plane.

It is natural to introduce a geodesic distance $l$ from the horizon as a cutoff parameter. At $t=0$, the metric is completely isotropic, so the cutoff is simply a circular region at constant $r$ on the $XZ$ plane. The shortest proper distance betwen the cutoff region and the horizon is then given by $l$. Outside the cutoff, we assume $\gamma_{A}$ lies on the $XZ$ plane. The contribution to $|\gamma_{A}|-|\gamma_{A}^{(0)}|$ in this region can be computed using simple numerical integration. By extrapolating the convergence of $\Delta S_{A}$ with $l$, we can calculate the $l\to\infty$ limit. Our cutoff procedure is visualized in figure 2 for two cutoff distances $l_{1}$ and $l_{2}$.

For $t>0$, the metric becomes increasingly anisotropic. We therefore replace the circular cutoff by an elliptic one, given by

[TABLE]

with $a$ and $b$ chosen so that the distance between the horizon and the cutoff is equal to $l$ along both the $X$ and $Z$ axes (corresponding to $(\theta,\phi)=(\pi/2,0)$ and $(0,0)$). Due to the symmetries of our metric, geodesics are still straight lines along these axes. An example for the elliptic cutoff at two different $l$ is shown in figure 3, using both global and Poincaré coordinates.

In our definition of Poincaré and global coordinates, we introduced the AdS radius $R$ and the Poincaré distance $\alpha$ determining the starting point of the falling mass $m$. For our numerical purposes, we can set $R=\alpha=1$. As mentioned earlier, $\alpha$ acts as a scaling factor in Poincaré coordinates, thus setting $\alpha=1$ is equivalent to considering time dependence with respect to $\tilde{t}={t\over\alpha}$. For the global coordinates in which the computations are performed, $R=1$ corresponds to rescaling the metric to $\text{d}\tilde{s}^{2}={\text{d}s^{2}\over R^{2}}$ and using a unitless mass parameter $\tilde{M}={M\over R^{3}}={2G_{N}m\over R}$. Thus we can write the growth of entanglement entropy in AdS4/CFT3 in the form

[TABLE]

where $\Delta\tilde{A}(\tilde{t},\tilde{M})$ is a unitless function of unitless parameters, to be determined numerically. In the following sections, we usually omit the rescaled units and express the results directly in terms of $\Delta S_{A}$. Note that in terms of CFT quantities, we can write

[TABLE]

In this equation, $\Delta$ is the conformal dimension of local operator $O(x)$ (1) and $c_{3}\equiv{R^{2}\over 2G_{N}}$ is a conventional measure of the degrees of freedom of the 3-dimensional holographic CFT, which is a generalization of the central charge in 2-dimensional CFTs. Thus in the CFT language we can write (24) as $\Delta S_{A}={c_{3}\over 2}\cdot\Delta\tilde{A}\left(\tilde{t},{\Delta\over c_{3}}\right)$.

3.3 Numerical results

At time $t=0$, when the mass is at rest in Poincaré coordinates, we can ignore variations in $t$ and need only to find a spatially minimal surface on a timeslice. This is because time reversal symmetry of the trajectory of the mass leads to a symmetry in the minimal surface in $t$. In addition, setting $t=0$ leads to a metric that is isotropic in global coordinates.

A series of computed discretized minimal surfaces at $t=0$ for different values of $M$ is shown in figure 4. At large $M$, the minimal surface “wraps” closely around the coordinate horizon of our metric (7). As this metric is only an approximation to a horizon-less metric (like that of a star) corresponding to a proper pure state, this means that results are not completely physical at very large $M$.

The rotational invariance and time-slice constraint allow us to simplify the algorithm from appendix A considerably, as we only have to find the shape of a 1-dimensional profile curve in 2-dimensional space. Note that at $t\neq 0$, both simplifications break down and we have to use the full algorithm.

The quantitative results for $\Delta S_{A}$ at $t=0$ are shown in figure 5. For small $M$, the data points closely follow a linear function with slope ${\pi\over 4}$ (red line), which is due to the first law of entanglement entropy Bhattacharya:2012mi ; Blanco:2013joa ; Wong:2013gua as will be analyzed in the next section. At larger $M$, this function provides an upper bound to the entanglement entropy, a consequence of the positivity of relative entropy Blanco:2013joa ; Wong:2013gua . Including constants $R$ and $G_{N}$, and using the dimension $\Delta$ of the local operator, this bound is expressed as

[TABLE]

We will explore this bound in more detail in the next section.

By inspecting the numerical solutions in figure 4, we can also obtain an analytical form of $\Delta S_{A}(t=0)$ in the large $M$ limit. As the minimal surface wraps closer around the horizon with increasing $M$, we can approximate its shape by a half-sphere around the horizon. Thus, we expect an approximate behavior

[TABLE]

where we inserted the area $A_{\text{hor}}$ of the horizon at radius $r_{\text{hor}}\simeq M^{1/3}$. As we can see in figure 5, our approximation (orange line) is valid at large $M$, as expected.444 We ignored terms from the annulus region around the half-sphere, as well as the $M=0$ contribution to $\Delta S_{A}$. Together, these add a subleading contribution of $-{3\pi\over 8}M^{1/3}$. For comparison, consider the exact AdS3 result 14 (dashed line in figure 5). Expanding it in powers of $M$ yields a similar expression, where $A_{\text{hor}}$ is replaced by the horizon diameter $d_{\text{hor}}=2\pi M^{1/2}$.555 Note that in terms of physical mass $m$, the parameter $M$ varies with dimension $d$ according to 8. This suggests that in AdSd+1/CFTd, the initial half-space entanglement entropy of a local quench at large $M$ can be written as

[TABLE]

For the $t>0$ case, we need to compute extremal surfaces in the full 4-dimensional spacetime. While boundary time $t$ is constant, it may vary on the rest of the surface. In order to visualize the solutions, we color-code vertices in their local time coordinate. The change in shape from the $t=0$ solution is shown in figure 6 for $M=1$. The surfaces are no longer isotropic and the local surface time differs considerably from boundary time $t$. Around $\cos\theta<0$, which corresponds to the $z<\sqrt{\alpha^{2}+t^{2}}$ region in Poincaré coordinates between falling mass and AdS boundary, the surface dips into the local future. It returns to boundary time $t$ at its closest distance to the mass, and then dips into the past around $\cos\theta>0$.

With increasing time $t$ at the boundary, the surface itself becomes more time-dependent. For large $t$ and $M$, parts of the extremal surface cross the horizon of the Poincaré coordinate patch, forcing us to perform all computations in global coordinates.

The numerical results for the time evolution of $\Delta S_{A}$ are shown in figure 7 for a range of mass parameters $M$. The exact result for the AdS3/CFT2 case (dashed lines), given by (13), is shown for the same values of $M$. For easier comparison, it has been rescaled to match the AdS4/CFT3 data at $t=0$. We find clear deviations between them, especially at late time.

We consider two possible scenarios at large $t$: An asymptotic power law $\Delta S_{A}(t)\propto t^{p}$, and a logarithmic time dependence $\Delta S_{A}(t)\propto\log t$, the latter of which we observed in the AdS3 case.

To test for an asymptotic power law, we use an estimator $\hat{p}$ of the power coefficient $p$ of a supposed asymptotic time dependence $f(t)=a+b\,t^{p}$. Assume a triplet of data points $\left((e^{t_{0}-\Delta t},f_{1}),\;(e^{t_{0}},f_{2}),\;(e^{t_{0}+\Delta t},f_{3})\right)$ is given. If $f(t)$ matches these data points, the estimated power coefficent follows as

[TABLE]

By Gaussian error propagation, the associated absolute error is

[TABLE]

in terms of the absolute errors $\delta f_{k}$ of the function values. Estimators for $a$ and $b$ can be constructed analogously.

The results for $f(t)=\Delta S_{A}(t)$ are shown in figure 8. In the limit $t\to 0$, the power coefficient converges to $\hat{p}=2$ for all $M$. For $t\to\infty$, however, $\hat{p}$ strongly depends on $M$: At small $M$, an apparently unstable plateau close to $\hat{p}=1$ is visible, while at large $M$, there appears to be convergence to a small but nonzero power. At large $t$, a strict bound $\hat{p}<1$ is satisfied, i.e. the entanglement entropy increases sub-linearly. Note that the initial $\hat{p}=2$ growth and $\hat{p}=1$ plateau follow from the first law of entanglement entropy, to be explained in the next section. The late time sub-linear growth is a genuinely new behavior, which we present in this paper for the first time.

In order to test a possible logarithmic time dependence $\Delta S_{A}(t)=a+b\,\log{t}$ at large $t$, we simply check ${\text{d}\,\Delta S_{A}(t)\over\text{d}\log{t}}$ for convergence. This is plotted in figure 9. While the AdS3 result converges to a constant, the numerical data for the AdS4 case shows no sign of slowing down within the range of $t$ that we computed.

Unfortunately, our numerical approach cannot be extended to arbitrarily large $t$, as the increasing time dependence requires a larger finite element resolution, slowing down the algorithm considerably. Thus, we can cannot exclude the possibility that the apparent asymptotic power law breaks down at very large $t$, which might lead to a logarithmic growth in the end. However, our AdS${}_{4}/$CFT3 analysis shows clear differences to the AdS${}_{3}/$CFT2 case, which extend to at least ${t\over\alpha}\approx 50$.

In fact, there are a few theoretical hints suggestive of an asymptotically logarithmic time evolution $\Delta S_{A}(t)\propto\log t$. For instance, in an analysis of Renyi entanglement entropy for local quenches of holographic CFTs we find a logarithmic time evolution Caputa:2014vaa in any dimensions $d$ when $\Delta\ll c_{d}$, where $c_{d}=O(N^{2})$ describes the degrees of freedom of the CFT. However, this analysis breaks down in the von Neumann entropy limit $n=1$. Moreover, the tensor network description of a falling particle proposed in Nozaki:2013wia for 2-dimensional CFTs is straightworward to generalize to any dimensions, also leading to a logarithmic growth of entanglement entropy.

In contrast, the theoretical background of the initial time evolution leads to much less ambiguity, as our observations fully agree with the first law of entanglement entropy.

4 Local Quenches and the First Law of Entanglement Entropy

We shall now give an interpretation of some parts of our numerical results presented in the previous section. Consider the relative entropy $S(\rho_{1}|\rho_{2})$ between two different quantum states, defined as

[TABLE]

where $\rho_{i}$ is the density matrix corresponding to the $i$th state. The positivity of $S(\rho_{1}|\rho_{2})$ leads to a bound to the entanglement entropy change Blanco:2013joa ; Wong:2013gua :

[TABLE]

where $\Delta\langle H_{A}\rangle$ denotes the change in expectation value of the modular Hamiltonian $H_{A}=-\log\rho_{A}$. In particular, if we consider a $d$-dimensional CFT and choose the subsystem $A$ to be a round ball with radius $l$, we can explicitly write the change of $H_{A}$ between an excited states and the CFT vacuum as follows Blanco:2013joa ; Wong:2013gua :

[TABLE]

Here $x_{0}$ is the center of the round ball $A$ and $T_{tt}$ is the energy density.

If we consider an infinitesimally small excitation, then the leading linear order contribution saturates the inequality (32), leading to

[TABLE]

This is called the first law of entanglement entropy Bhattacharya:2012mi ; Blanco:2013joa ; Wong:2013gua , relating a change in energy (or energy density) to a change in entanglement entropy.

4.1 AdS3/CFT2

The energy density of local quenches Nozaki:2013wia can be obtained in the AdS3 setup as

[TABLE]

The subsystem $A$ is taken as an interval $x\in[0,2l]$ and the local excitation is situated at $x=0$. The first law (34) leads to

[TABLE]

We now take $l\to\infty$, extending $A$ to the entire half-line. We find

[TABLE]

Using the identity Gubser:1998bc_Witten:1998qj

[TABLE]

we can rewrite this result for AdS3/CFT2 ($d=2$) in terms of the mass parameter $M$:

[TABLE]

Thus, the first law predicts quadratic growth of $\Delta S_{A}$ at early time and a transition to linear growth at later time. However, in order to apply the first law we require $\Delta S_{A}\ll c$, where $c$ is the central charge Bhattacharya:2012mi . At small $t$, (37) turns this into the condition $\Delta\ll c$. At larger $t$, we get the additional constraint

[TABLE]

where we used the holographic relation ${\Delta\over c}={M\over 12R^{2}}$. If both constraints are met, our results from the first law relations can be trusted. On the other hand, when ${t\over\alpha}\gg{c\over\Delta}$ (late time zone), we have the logarithmic grow $S_{A}\simeq{c\over 6}\log t$, derived by both holographic and field theoretic methods Nozaki:2013wia ; Asplund:2014coa .

The exact results for $\Delta S_{A}$ in AdS3/CFT2, along with the AdS4/CFT3 computations, are shown in figure 8 in terms of a local power fit (29) with power coefficent $\hat{p}$. For a small mass parameter $M$, the first law requirement $\Delta S_{A}\ll c$ holds and we observe a quadratic growth ($\hat{p}=2$) at early time. In the intermediate region $1\ll{t\over\alpha}\ll{c\over\Delta}$ , the predicted transition to linear growth ($\hat{p}=1$) can be seen. With increasing $M$, the time dependence becomes logarithmic more quickly, with $\hat{p}$ converging to zero.

4.2 AdS4/CFT3

In the AdS4/CFT3 case ($d=3$) the energy density looks like

[TABLE]

with a radial coordinate $\rho=\sqrt{x^{2}+y^{2}}$. First, consider the configuration at small time $t\ll\alpha$, where the energy density reads

[TABLE]

We choose the subsystem $A$ to be a disk with a radius $l$ whose center is at the point $(x,y)=(l,0)$. At $t=0$, the change of the modular Hamiltonian is

[TABLE]

where in the final expression we used $l\to\infty$. This leads to the bound

[TABLE]

which reproduces the bound (26) that we observed in figure 5.

The calculation steps in (43) can extended to the $O(t^{2})$ term in (42), yielding

[TABLE]

as the bound to entanglement entropy growth at small $t$.

Next we assume $\alpha\ll t\ll l$. Then the integral of the first law (34) has a dominant contribution around $\rho\simeq t$. Thus it is approximated as follows:

[TABLE]

Hence we find

[TABLE]

The saturation of this inequality (i.e. the first law) reproduces the perturbative holographic result (22). This is to be expected, as both the first law and the perturbative approach are valid in the regime of small excitations, i.e. when (32) only contains subleading contributions with respect to changes in the density matrix describing the system. In our setup, this corresponds to small $M$, where the shape of the extremal surface is not greatly perturbed by the mass. Under these conditions, we expect the qualitative time dependence of $\Delta S_{A}$ at early and intermediate time $t$ to be equivalent to the AdS3/CFT2 case, i.e. with an initial quadratic growth at ${t\over\alpha}\ll 1$ and a linear behavior at intermediate time. Indeed, that is what we already saw for the numerical data in figure 8.

Even at $t=0$, the first law bound is not reliable if $\Delta$ is large. Thus, the behavior

[TABLE]

that we observed at large $M$ in figure 5 cannot be understood in terms of the first law.

The time-dependent bounds (45) and (47) (for $t\ll\alpha$ and $t\gg\alpha$, respectively) are shown in figure 10 alongside numerical data at small $M$. In their respective $t$ range, they provide an upper bound for all $M$ and saturate as $M\to 0$, as expected. As in the AdS3/CFT2 case, we expect the saturation to break down at very large $t$ as the entanglement entropy reaches its asymptotic, sub-linear growth.

5 Summary and Discussion

In this work, we studied the holographic entanglement entropy of an AdS4/CFT3 setup involving a local quench. To this end, we employed a numerical finite-element approach to calculate extremal surfaces in the AdS4 geometry dual to the CFT3. In particular, we considered the growth of entanglement entropy $\Delta S_{A}$ corresponding to a half-space subsystem $A$, on the boundary of which the local quench occurs. The strength of the quench is determined by the operator dimension $\Delta$ of the excitation, which is proportional to a mass $m$ that is falling freely in the bulk theory Nozaki:2013wia .

At time $t=0$ after the quench, the numerical data (see figure 5) is constrained by a bound

[TABLE]

which saturates at small $\Delta$, in agreement with the first law of entanglement entropy. At large $\Delta$, the extremal surface area is dominated by the region near the coordinate horizon around the mass, leading to the relation

[TABLE]

where $A_{\text{hor}}$ is the horizon area and $R$ the AdS radius. As a similar relationship can be derived for the AdS3/CFT2 case, a relation

[TABLE]

can be obtained for local quenches in holographic $d$-dimensional CFT in the large $\Delta$ limit.

At small $t$ and $\Delta$, our numerical results exhibit an initial quadratic bound that is again in agreement with the first law prediction

[TABLE]

At larger $t$, this turns into a bound linear in $t$:

[TABLE]

Again, both bounds are saturated at small $\Delta$ (see in figure 10).

Qualitative differences to the entanglement growth in 2-dimensional CFT become apparent at large $t$ and $\Delta$: Instead of quickly reaching an asymptotically logarithmic time dependence, we observe a growth similar to a power law $\Delta S_{A}\sim t^{p}$ (see figure 8 and 9) in the late time region of our analysis.

This observed time dependence may not be truly asymptotical, as our numerical methods cannot be extended to arbitrarily large time $t$. However, it is clear that significant deviations from the entanglement entropy in lower dimensions appear with increasing $t$. It is also interesting to note that the power $p$ (shown by an estimator $\hat{p}$ in figure 8) monotonically decreases toward zero as the time $t$ and the mass parameter $M$ (and thus $\Delta$) get larger. This suggests that the asymptotic limit is either $\Delta$-dependent, or slower than a power law, i.e. logarithmic in $t$. The latter case would be more consistent with arguments from Renyi entropies Caputa:2014vaa and tensor networks Nozaki:2013wia , but clearly, complete analytical studies of entanglement entropy evolution in higher-dimensional CFTs will be needed to explain our results. This would require further characterization of holographic CFTs in the language of field theory.

We would also like to mention that in our holographic approach Nozaki:2013wia , we regard the gravity dual of local quench as a massive heavy particle falling in the bulk AdS. This treatment precisely corresponds to the large $c$ approximation Asplund:2014coa . However, more precisely it should be described by a time-dependent classical solution made of a scalar field in the bulk AdS with gravitational backreactions as in Rangamani:2015agy . Such a detailed structure of massive excitation in the bulk may affect the late time evolution of holographic entanglement entropy, which is dual to the breakdown of the standard large $c$ approximation. It is a very interesting future problem to take this effect into account to see if the entanglement entropy in the late time limit can approach a finite constant or continue to grow forever.

In principle, the finite-element optimization strategy used in this paper can be applied to a large class of problems involving holographic entanglement entropy in AdS4/CFT3. While our current implementation (details in appendix A) only considers simply connected boundary regions, extensions to more complicated models are possible. As analytical calculations in CFT3 are notoriously complicated, our numerical technique can offer an alternative approach.

Acknowledgements

We would like to thank Vijay Balasubramanian, David Berenstein, Pawel Caputa, Jens Eisert, Mario Flory, Masairo Nozaki, Tokiro Numasawa, Erik Tonni and Herman Verlinde for useful discussions, and especially to Robert Myers for insightful comments. TT is supported by the Simons Foundation through the “It from Qubit” collaboration and JSPS Grant-in-Aid for Scientific Research (A) No. 16H02182. TT is also supported by World Premier International Research Center Initiative (WPI Initiative) from the Japan Ministry of Education, Culture, Sports, Science and Technology (MEXT). AJ is supported by the German Academic Scholarship Foundation (Studienstiftung des deutschen Volkes). His research exchange to Kyoto University was supported by a grant of the German Academic Exchange Service (DAAD). AJ is also grateful to Valentina Forini for supporting the thesis on which this work is partly based. We are very grateful to the long term workshop “Quantum Information in String Theory and Many-body Systems” held at YITP in Kyoto University and the IGST 2016 conference held at Humboldt University, where this work was partially conducted.

Appendix A Finite Element Implementation

This appendix covers the main features of our algorithm for finding the area of an extremal 2-surface with a given boundary in 3+1 spacetime dimensions. We implemented the algorithm in C++ without the use of any existing framework. Visualizations of the output surface data were produced with Mathematica.

We use a nomenclature customary in graphics programming: A discretization point is called a vertex, a line segment between two vertices is referred to as an edge, and a surface element is a face. Faces are typically triangles with three corner vertices or quadrilaterals, also called quads, with four. The full geometry is called a mesh.

In terms of data structures, a vertex is a collection of four floating-point numbers, one for each coordinate. An edge, a triangle, and a quadrilateral are collections of two, three, and four integers, respectively. Each integer serves as an index pointing to a vertex. Thus, only the vertices are dynamical objects, while edges and faces are derived objects used to calculate lengths and areas of the geometry.

The type of numerical approach used here is a finite element method, where a continuous problem is solved on a discretized mesh (for an introduction, see e.g. Zienkiewicz1977 ; Hughes2000 ). Shape optimization problems of this type are often encountered in engineering, as in the minimization of the strain on a component. In our case, the functional to be optimized is the area of the surface itself. However, extending this approach to surfaces in 4-dimensional curved spacetime leads to several complications, described throughout the rest of this appendix.

A.1 Area calculation

The primitive surface element for computing areas is the triangle. A triangle spanned by two differential vectors $\text{d}v^{\mu}$ and $\text{d}w^{\mu}$ on a Riemannian manifold with metric $g_{\mu\nu}$, has an area of

[TABLE]

where we used Lagrange’s identity for the wedge product $a\land b$ to get to the second step. If we parametrize the differential vectors as the derivative of a function $x^{\mu}$ along two coordinates $t$ and $u$, the result is equivalent to half of the familiar differential area element:

[TABLE]

where $G_{\mu\nu}$ is the induced metric on the surface parameterized by $x^{\mu}(u,v)$. In principle, we could integrate (A.1) over $u$ and $v$ to yield the exact area formula of a triangle in a spacetime with metric $g_{\mu\nu}$. However, this metric is usually too complicated to yield a compact analytical expression. Instead, we can assume that each face is small enough so that $g_{\mu\nu}$ varies little along its surface. Then we simply use (A.1) with $\text{d}v^{\mu}\to p_{1}^{\mu}-p_{3}^{\mu}$ and $\text{d}w^{\nu}\to p_{2}^{\nu}-p_{3}^{\nu}$ as the formula for the area of a triangle with corner vertices $p_{1}$, $p_{2}$ and $p_{3}$. As we subdivide the mesh into smaller and smaller elements, the error associated with this approximation gradually decreases.

Clearly, an area can only be defined for a space-like triangle, i.e. one for which (A.1) is real. However, the parameter space of all dynamic vertices allows for unphysical, time-like solutions. Approximating edges as differential vectors can also lead to configurations where (A.1) turns imaginary. While choosing a coordinate system that minimizes time dependence alleviates this problem (see A.5), it still forces us to use an optimization algorithm that avoids locally time-like solutions (outlined in A.3). First, however, we have to create a triangular discretization of the mesh we seek to extremize.

A.2 Triangulation and quadrilateralization

The boundary of the extremal surface is given as a closed chain of static vertices, chosen to approximate a continuous boundary function, between which we construct an initial mesh. After adding an additional vertex in the center of the boundary vertices and connecting every second boundary vertex with the center via an edge, we can define quads that produce a complete mesh.

The algorithm proceeds through $N_{i}$ iterations of first optimizing the mesh at a given discretization level and then refining it by subdividing the mesh into smaller quads. In order to calculate the area of a quad using (A.1), we need to divide it into a pair of triangles. However, depending on which quad diagonal is chosen as the separation between both triangles, there are two distinct ways of filling a quad. As the total surface has an orientation, these two solutions correspond to making the surface either locally convex or concave. As we subdivide the mesh into ever smaller quads and converge to a smooth surface, both solutions become equal. Therefore, we can choose the area of each quad as the average of the area of its two triangular fillings. From the perspective of the entire mesh, we average between an “outer” and “inner” triangulation that become more degenerate with each iteration step.

The subdivision process of an individual quad is visualized in figure 11. Each edge is divided into two, and an additional vertex is added in the center. Thus, each quad is divided into four smaller quads. The position of the new vertices is given by a weighted average of two (edge vertices) or four (center vertex) 4-vectors of the old vertices. The choice of weights depends on the problem: For calculating absolute areas, subdivision vertices should be chosen so that the new quads have roughly equal area. For the calculation of relative areas (i.e. between different metrics for similar boundary conditions), faster convergence can be reached by associating a greater weight to regions were deviations between both problems are largest.

Note than for the purposes of area calculation and subdivision, we treat edges as straight lines, a notion that depends on the choice of coordinates. Working with coordinate-independent geodesics is much more computationally demanding, if geodesic solutions are not explicitly known. As the edge length decreases exponentially with the number of iterations, the discrepancy quickly vanishes as long as the solution is sufficently smooth.

A full example of the quadrilateralization of a mesh with 12 initial boundary vertices is shown in figure 12. Note that boundary edges are subdivided so that the new vertices follow the continuous boundary function.

The use of a quadrilateralized mesh instead of a simpler triangulated one is preferable due to its stability with respect to subdivision. As quads are locally agnostic with respect to concavity, the optimization algorithm cannot get trapped in a local concave minimum. The problem is visualized in figure 13: When a coarse triangulated mesh is optimized around a sharply peaked solution, the result can have a jagged shape where concavity varies greatly between neighboring triangles. After subdividing the triangles to achieve better resolution, the triangles within the concave region would have to be “flipped” to make the surface locally convex again. This problem is altogether avoided by using quads instead of triangles.

A.3 Numerical optimization

We require a numerical method for finding a $D$-dimensional vector $\vec{x}=(x_{1},\dots,x_{D})$ extremizing a function $f(\vec{x})$, i.e. finding a point $\tilde{\vec{x}}$ with

[TABLE]

A textbook approach to this problem is Newton’s method666 The term “Newton’s method” is sometimes used exclusively for the original zero-point search method from which the extremum search method directly follows. (see e.g. BurdenFairesCh10_2 ; UeberhuberCh14_4_1 ; PolakCh1_4 ). This iterative method starts from an initial point $\vec{x}^{(0)}$, onto which a recursion step is applied. This formula is given by

[TABLE]

with the Hessian matrix $H$ defined by

[TABLE]

The formula (57) is computationally demanding for large $D$, as it requires an $O(D^{2})$ computation of $H$ as well as an $O(D^{3})$ matrix inversion. For our 4-dimensional problem with $N_{V}$ dynamic vertices describing the discretized geometry, we would set $D=4N_{V}$ for a naive implementation of Newton’s method. Clearly, such an approach is not computationally efficient, and a variety of approximate quasi-Newton methods are used. For example, in Fonda:2014cca the popular conjugate gradient method is used for similar problems on fixed time-slices.

Discretized surfaces in 4-dimensional spacetime live in a constrained region of the $4N_{V}$-dimensional parameter space we wish to probe numerically, as we have to avoid locally time-like solutions. On any quasi-Newton method that modifies several of the $4N_{V}$ degrees of freedom at once, this imposes a set of complicated dynamical constraints. To avoid this, we use Newton’s method in its exact form but apply it locally, making the constraints easier to implement.

The function to be extremized is the total area $A=\sum_{k}A_{\square_{k}}$ of all quads $\square_{k}$. While the total area depends on all $N_{V}$ dynamic vertices, varying the coordinates of one vertex $v^{\mu}_{i}$ alone only changes the area of quads that contain the vertex. In order to make the geometry locally extremal, we only need to apply Newton’s method to one vertex at a time. In other words, instead of a set of global conditions for extremal total area we use an equivalent set of local conditions

[TABLE]

where the sum runs over all quads $\square_{k}$ that contain the vertex $v_{i}$. The corresponding recursion step for Newton’s method turns from (57) into

[TABLE]

which is applied to each vertex 4-vector $v^{\mu}_{i}$ in the mesh. The Hessian is now only a $4\times 4$ matrix which can be inverted without much computational cost. Upon introducing normal vectors, the Hessian is further simplified to a $2\times 2$ matrix. If the recursion step leads to a locally time-like solution, the step size $\|v_{i}^{(n+1)}-v_{i}^{(n)}\|$ is reduced until the solution becomes space-like again.

The recursion step 60 is repeated until the local gradient norm $\|g_{i}\|$ is smaller than a threshold value.777 For an inhomogeneous discretization, it is preferable to apply the threshold to the gradient density $\|g_{i}\|/A_{v_{i}}$ instead, with $A_{v_{i}}=\sum_{\square_{k}\ni v_{i}}A_{\square_{k}}$ being the area of all quads containing the vertex $v_{i}$. By repeatedly optimizing all vertices, the gradient norm vanishes across the entire mesh, with $A$ converging to the extremal area.

This approach of successively optimizing each vertex to extremize its local environment is computationally efficient only if applied to a geometry that is already close to the exact solution. Otherwise, the number of local updates needed to fulfill condition (59) for all vertices $v_{i}$ will become strongly nonlinear in $N_{V}$. Recursively subdividing and optimizing avoids such problems: The first few iterations establish the rough shape of the solution by varying only a few dynamic vertices. Every subsequent iteration gradually refines the geometry, successively reducing the computational optimization cost per vertex.

Generally, Newton’s method could be applied to the $N_{V}$ vertices in any order. However, we find it is more efficient to sort the vertices according to their local gradient norm $\|g_{i}\|$ and always optimize the vertex with the largest gradient. After the vertex is optimized, the vertices in the local neighborhood have their entries in the gradient list updated, and the process is repeated. This gradient update only has a computational cost of $O(\log N_{V})$, but allows the algorithm to focus on the part of the mesh that is least extremal, speeding up convergence.

Our formulation of Newton’s method is not fully covariant: The gradient vector, i.e. the covariant derivative of the scalar area functional, is equivalent to an ordinary partial derivative, so the fixed points of the method are independent of the coordinate system chosen. The Hessian, however, is not. While fully covariant descriptions of Newton’s method exist Dedieu2003 , the computational cost of computing Christoffel symbols at every step does not make this approach favorable.

A.4 Normals

While there are now four dynamic parameters per vertex to optimize, we would like to constrain the optimization to directions orthogonal to the local surface. This is because moving the vertex along a tangent direction changes the resolution of the discretization, which can lead to clustering of vertices in one region of the mesh. To keep the discretization points at a distance, we restrict Newton’s method to directions along normal vectors.

An example for 3-dimensional vertex normals is shown in figure 14. Each dynamic (i.e. non-boundary) vertex lies along the corners of four quads, and each corner has an orthogonal direction given by the vector product of the two edges involved. Averaging those over all corners, an effective vertex normal is produced.

In curved 4-dimensional spacetime, we need to consider differential edges in the local neighborhood of each vertex. Furthermore, a 2-dimensional surface in four dimensions has a 2-dimensional space of 4-vectors orthogonal to each point on the surface, so normal vectors cannot be uniquely defined. We are therefore free to choose a pair of 4-dimensional (differential) normals $\text{d}n_{\mu}$ and $\text{d}m_{\mu}$ in the following manner:

[TABLE]

Here, $\text{d}v^{\mu}$ and $\text{d}w^{\mu}$ are two differential vectors spanning a corner of a quad, and $\epsilon_{\mu\nu\rho\sigma}$ is the 4-dimensional Levi-Civita symbol. The constant vector $\text{d}t^{\mu}$ is an arbitrary covariant time-like vector. Using differential forms, we can also write (61) more generally as

[TABLE]

where $\text{d}t$ is now a differential 1-form. In practice, we can simply set $\text{d}t=\delta_{0\mu}\text{d}x^{\mu}$.

In this construction, $\text{d}n_{\mu}$ is a space-like vector and $\text{d}m_{\mu}$ a time-like one. Thus, an extremal surface will be minimal along the former and maximal along the latter direction. For a surface at constant time, i.e. with $\text{d}v^{0}=\text{d}w^{0}=0$, the first normal $\text{d}n_{\mu}$ corresponds to the usual 3-dimensional normal vector with a vanishing time component while the second normal $\text{d}m_{\mu}$ only points in the time direction.

For computational purposes, normal vectors are simply vectors of unit length, as we are only interested in the direction of the normals. Due to this rescaling, we have omitted a factor $\sqrt{|\det g|}$ in (61) that would appear in a covariant definition of the Levi-Civita symbol.

A.5 Time-slice constraints

For strongly time-dependent solutions, it may not be possible during the first few iterations to find a low-resolution mesh that follows the general shape of the smooth solution. This can lead to unstable conditions after subdividing the quad mesh, such as time-like quads. This problem can be avoided using a series of preconditioning steps, during which we restrict the mesh to a time-slice compatible with the boundary conditions. When the mesh has become sufficiently smooth, the time-slice constraint is dropped.

A suitable choice of the metric $g_{\mu\nu}$ also increases the stability of the algorithm. Assume the initial coordinates are $(t,x,y,z)$. If the mesh boundary lies on some space-like surface given by $t=f(x,y,z)$, we can simply introduce a modified time coordinate $t\to t^{\prime}=t-f(x,y,z)$. In the new $(t^{\prime},x,y,z)$ coordinates, we can now construct an initial mesh at $t^{\prime}=0$. Under the new metric, tangent vectors along any point on this mesh will lie on the space-like surface, so we can perform a time-slice constrained optimization without encountering time-like solutions. After these preconditioning steps, we allow the algorithm to optimize along the $t^{\prime}$ coordinate, as well.

A.6 Accuracy and error estimation

The accuracy of the area of the extremized surface mainly depends on two parameters, the gradient convergence threshold and the number of iterations. The threshold value serves to terminate both the local steps of Newton’s method as well as the entire iteration, where in the latter case it is applied to the maximum gradient along the whole mesh. The absolute value of this threshold should be decreased with each iteration step, so that smoothness on increasingly smaller scales is achieved. For our purposes, we used a threshold gradient density of $10^{-2}$.

Given a sufficiently small gradient threshold, the area of the discretized surface converges exponentially with the number of iterations, as each iteration increases the area resolution by a factor of four. Denote as $\Delta A^{(i)}=|A^{(i)}-A^{(i-1)}|$ the difference between the extremal surface area at iteration $i$ and $i-1$. If $A^{(i)}$ converges exponentially with $i$, so does $\Delta A^{(i)}$. An example for such a convergence is shown in figure 15. Note that the convergence is more erratic during earlier iterations, as the extremal surface still changes shape considerably. After a few iterations, however, the convergence quickly becomes exponential. Taking the differences $\Delta A^{(N_{i})}$ and $\Delta A^{(N_{i}-1)}$ between the last three iterations, and assuming that $\Delta A^{(i)}$ continues to converge exponentially, we can estimate the absolute error $\delta A^{(N_{i})}=|A^{(\infty)}-A^{(N_{i})}|$ after $N_{i}$ iterations as

[TABLE]

where we sum up the projected steps of all further iterations. For example, the data in figure 15 leads to an estimated error $\delta A^{(10)}\approx 3\cdot 10^{-4}$. If $A^{(i)}$ does not converge monotonically (but $\Delta A^{(i)}$ still decreases exponentially), (63) gives an upper bound to the absolute error.

When computing the differences of surfaces in the limit of infinitely large boundaries, we also need to consider the results for different effective radii $l$. As long as $l$ is chosen to correspond to some proper length under the given metric, these results typically converge exponentially with $l$ as well, so we can compute errors similar to (63).

A.7 Performance

As the computational cost of one local update is constant, the performance of the algorithm scales linearly with the number of local updates required to reach the gradient convergence threshold. The iterative approach of optimizing the mesh at gradually higher resolutions implies that optimizing one vertex only affects the gradient of vertices within an effective local region $A_{\text{loc}}\propto\left({1\over 4}\right)^{k}$ of the entire mesh at the $k$th iteration. Because the number of vertices within $A_{\text{loc}}$ remains independent of $k$, the necessary number of local updates per vertex is constant as well. Thus, the performance of the algorithm should scale linearly in the number of vertices to be optimized, i.e. with $O(4^{N_{i}})$, $N_{i}$ being the number of total iterations. In reality, the actual performance drop is slightly higher, as the memory requirements also increase exponentially, reducing access times in any practical implementation. Note that processes such as normal construction and mesh subdivision have no noticeable impact on total runtime, as they are only executed once per iteration.

For producing the data presented in this paper, we used up to $N_{i}=10$ iterations, which corresponds to $\sim 10^{6}$ discretization points and as many quads. On a typical office CPU, one such calculation takes several days to complete and requires up to $\sim 1\mathrm{GB}$ of RAM.