Deep Quantum Geometry of Matrices

TL;DR

This paper uses machine learning to find accurate wavefunctions for matrix quantum mechanics, revealing how emergent geometries like fuzzy spheres behave quantum mechanically and their entanglement properties.

Contribution

It introduces a variational quantum Monte Carlo approach with deep generative flows to study gauge invariant states in matrix quantum mechanics, including supersymmetric extensions.

Findings

Recovered semiclassical fuzzy sphere states

Probed collapse of geometries in quantum regimes

Observed boundary-law entanglement entropy in large N limit

Abstract

We employ machine learning techniques to provide accurate variational wavefunctions for matrix quantum mechanics, with multiple bosonic and fermionic matrices. Variational quantum Monte Carlo is implemented with deep generative flows to search for gauge invariant low energy states. The ground state, and also long-lived metastable states, of an $\mathrm{SU}(N)$ matrix quantum mechanics with three bosonic matrices, as well as its supersymmetric `mini-BMN' extension, are studied as a function of coupling and $N$. Known semiclassical fuzzy sphere states are recovered, and the collapse of these geometries in more strongly quantum regimes is probed using the variational wavefunction. We then describe a factorization of the quantum mechanical Hilbert space that corresponds to a spatial partition of the emergent geometry. Under this partition, the fuzzy sphere states show a boundary-law entanglement entropy in the large $N$ limit.

Equations

$$H_B=\mathrm{tr}\left(\frac{1}{2}{\Pi }^i{\Pi }^i-\frac{1}{4}[X^i,X^j][X^i,X^j]+\frac{1}{2}{\nu }^2{X}^i{X}^i+i\nu {ϵ}^{ijk}{X}^i{X}^j{X}^k\right).$$

$$[X^i,X^j]=i\nu {ϵ}^{ijk}{X}^k.$$

$$H=H_B+\text{tr}\left({\lambda }^{†}{\sigma }^k[X^k,\lambda ]+\frac{3}{2}\nu {\lambda }^{†}\lambda \right)-\frac{3}{2}\nu (N^2-1).$$

$$H=$$

$$+i{f}_{ABC}{\lambda }_A^{\alpha †}{X}_B^k{\sigma }_{\alpha }^{k\beta }{\lambda }_{C\beta }+\frac{3}{2}\nu {\lambda }_A^{\alpha †}{\lambda }_{A\alpha }-\frac{3}{2}\nu (N^2-1),$$

$$Q_{\alpha }=\left(-i\frac{\partial }{\partial X_A^i}+i\nu X_A^i-\frac{i}{2}{f}_{ABC}{ϵ}_{ijk}{X}_B^j{X}_C^k\right){\sigma }_{\alpha }^{i\beta }{\lambda }_{A\beta },{\bar{Q}}^{\alpha }=(Q_{\alpha }{)}^{†},$$

$$\{Q_{\alpha },{\bar{Q}}^{\alpha }\}=4H.$$

$$G_A=-i{f}_{ABC}\left(X_B^i\frac{\partial }{\partial X_C^i}+{\lambda }_B^{\alpha †}{\lambda }_{C\alpha }\right).$$

$$|M\rangle\equiv\sum_{r=1}^{D}\prod_{a=1}^{R}\Big{(}\sum_{\alpha=1}^{2}\sum_{A=1}^{N^{2}-1}M_{A\alpha}^{ra}\lambda_{A}^{\alpha\dagger}\Big{)}|0\rangle,$$

$$(U\psi )(X)=f(U^{-1}XU)|(UM{U}^{-1})(U^{-1}XU)⟩,$$

$$f(X)=\tilde{f}(\tilde{X}),|M(X)⟩=|U\tilde{M}(\tilde{X})U^{-1}⟩\text{ where }X=U\tilde{X}{U}^{-1}.$$

$$d^{3(N^2-1)}X=\Delta (\tilde{X})d^{2(N^2-1)}\tilde{X},$$

$$\Delta (\tilde{X})\propto \prod _{i\ne j=1}^N\left|{\tilde{X}}_{ii}^1-{\tilde{X}}_{jj}^1\right|\prod _{i=1}^{N-1}\left|{\tilde{X}}_{i(i+1)}^2\right|.$$

$$E_{\theta }=⟨{\psi }_{\theta }|H|{\psi }_{\theta }⟩=\int dX|{\psi }_{\theta }(X){|}^2{H}_X[{\psi }_{\theta }]={\mathbb{E}}_{X\sim |{\psi }_{\theta }{|}^2}[H_X[{\psi }_{\theta }]].$$

$${\theta }_{t+1}={\theta }_t-\alpha {\nabla }_{{\theta }_t}{E}_{{\theta }_t},$$

$${\nabla }_{\theta }{E}_{\theta }={\mathbb{E}}_{X\sim p_{\theta }}[{\nabla }_{\theta }{H}_X[{\psi }_{\theta }]]+{\mathbb{E}}_{X\sim p_{\theta }}[{\nabla }_{\theta }\left(\ln p_{\theta }(X)\right)\left(H_X[{\psi }_{\theta }]-E_{\theta }\right)].$$

$$F=A_{\theta }^m\circ \tanh \circ A_{\theta }^{m-1}\circ \tanh \circ \cdots \circ \tanh \circ A_{\theta }^1.$$

$$q(y)=p(x){\left|\det \frac{Dy}{Dx}\right|}^{-1}=p(F^{-1}(y))|\det DF{|}^{-1}.$$

$$q(x)=\prod _{i=1}^M{p}_i(x_i;F_i(x_1,\dots ,x_{i-1})).$$

$$R=\sum _{A\alpha }{\lambda }_A^{\alpha †}{\lambda }_{A\alpha },[R,H]=0,$$

$$p_i(x_i)=\sum _{r=1}^s{k}_r^i\frac{{\beta }_r^i}{2{\alpha }_r^i\Gamma (1/{\beta }_r^i)}{e}^{-(|x_i-{\mu }_r^i|/{\alpha }_r^i{)}^{{\beta }_r^i}},\sum _{r=1}^s{k}_r^i=1.$$

$$p_i(x_i>0)=\sum _{r=1}^s{k}_r^i\frac{({\beta }_r^i{)}^{{\alpha }_r^i}}{\Gamma ({\alpha }_r^i)}(x_i{)}^{{\alpha }_r^i-1}{e}^{-{\beta }_r^i{x}_i},\sum _{r=1}^s{k}_r^i=1.$$

$$r=\sqrt{\frac{1}{N}\mathrm{tr}(X_1^2+X_2^2+X_3^2)},$$

$$[X^i,X^j]=i\nu {ϵ}^{ijk}{X}^k.$$

$$X^i=\nu J^i,$$

$$r^2=\frac{1}{N}\sum _{i=1}^3\mathrm{tr}(X^i{)}^2=\frac{{\nu }^2(N^2-1)}{4}.$$

$$\sum _{i=1}^3[J^i,[J^i,{\hat{Y}}_{jm}]]=j(j+1){\hat{Y}}_{jm},[J^3,{\hat{Y}}_{jm}]=m{\hat{Y}}_{jm}.$$

$$X^i=\nu J^i+\sum _{jm}{y}_{jm}^i{\hat{Y}}_{jm},$$

$${\omega }^2=0$$

$${\omega }^2=j^2$$

$${\omega }^2=(j+1{)}^2$$

$${\Delta E|}_{\text{bos}}=\frac{|\nu |}{2}\sum |\omega |=\frac{4N^3+5N-9}{6}|\nu |.$$

Full text

1. 1 Introduction

2. 2 The mini-BMN model

3. 2.1 Representation of the fermion wavefunction

4. 2.2 Gauge invariance and gauge fixing

5. 3 Architecture design for matrix quantum mechanics

6. 3.1 Parametrizing and sampling the gauge invariant wavefunction

7. 3.2 Benchmarking the architecture

8. 4 The emergence of geometry

9. 4.1 Numerical results, bosonic sector

10. 4.2 Semiclassical analysis of the fuzzy sphere

11. 4.3 Numerical results, supersymmetric sector

12. 5 Entanglement on the fuzzy sphere

13. 5.1 Free field with an angular momentum cutoff

14. 5.2 Fuzzy sphere in the mini-BMN model

15. 6 Discussion

16. A Geometry of the gauge

17. B Evaluation of observables

18. C Semiclassical analysis of the fuzzy sphere

19. D Training and tuning

20. E Entanglement of free fields on a sphere

1 Introduction

A quantitative, first principles understanding of the emergence of spacetime from non-geometric microscopic degrees of freedom remains among the key challenges in quantum gravity. Holographic duality has provided a firm foundation for attacking this problem; we now know that supersymmetric large $N$ matrix theories can lead to emergent geometry [1, 2]. What remains is the technical challenge of solving these strongly quantum mechanical systems and extracting the emergent spacetime dynamics from their quantum states. Recent years have seen significant progress in numerical studies of large $N$ matrix quantum mechanics at nonzero temperature. Using Monte Carlo simulations, quantitatively correct features of emergent black hole geometries have been obtained, e.g. [3, 4, 5]. To grapple with questions such as the emergence of local spacetime physics, and its associated short distance entanglement [6, 7], new and inherently quantum mechanical tools are needed.

Variational wavefunctions can capture essential aspects of low energy physics. However, the design of accurate many-body wavefunction ansatze has typically required significant physical insight. For example, the power of tensor network states, such as Matrix Product States, hinges upon an understanding of entanglement in local systems [8, 9]. We are faced, in contrast, with models where there is an emergent locality that is not manifest in the microscopic interactions. This locality cannot be used a priori; it must be uncovered. Facing a similar challenge of extracting the most relevant variables in high-dimensional data, deep learning has demonstrated remarkable success [10, 11, 12], in tasks ranging from image classification [13] to game playing [14]. These successes, and others, have motivated tackling many-body physics problems with the machine learning toolbox [15]. For example, there has been much interest and progress in applications of Restricted Boltzmann Machines to characterize states of spin systems [16, 17, 18, 19].

In this work we solve for low-energy states of quantum mechanical Hamiltonians with both bosons and fermions, using generative flows (normalizing flows [20, 21, 22] and masked autoregressive flows [23, 24, 25] in particular) and variational quantum Monte Carlo. Compared with spin systems, the problem we are trying to solve contains continuous degrees of freedom and gauge symmetry, and there is no explicit spatial locality. Recent works have applied generative models to physics problems [26, 27, 28] and have aimed to understand holographic geometry, broadly conceived, with machine learning [29, 30, 31]. We will use generative flows to characterize emergent geometry in large $N$ multimatrix quantum mechanics. As we have noted above, such models form the microscopic basis of established holographic dualities.

We will focus on quantum mechanical models with three bosonic large $N$ matrices. These are among the simplest models with the core structure that is common to holographic theories. The bosonic part of the Hamiltonian takes the form

[TABLE]

Here the $X^{i}$ are $N$ by $N$ traceless Hermitian matrices, with $i=1,2,3$. The $\Pi^{i}$ are conjugate momenta and $\nu$ is a mass deformation parameter. The potential energy in (1.1) is a total square: $V(X)=\frac{1}{4}\operatorname{tr}\left[\left(\nu\epsilon^{ijk}X^{k}+i[X^{i},X^{j}]\right)^{2}\right]$. The supersymmetric extension of this model [32], discussed below, can be thought of as a simplified version of the BMN matrix quantum mechanics [33]. We refer to the supersymmetric model as ‘mini-BMN’, following [34]. For the low energy physics we will be exploring, the large $N$ planar diagram expansion in this model is controlled by the dimensionless coupling $\lambda\equiv N/\nu^{3}$. Here $\lambda$ can be understood as the usual dimensionful ’t Hooft coupling of a large $N$ quantum mechanics at an energy scale set by the mass term (cf. [35]).

The mass deformation in the Hamiltonian (1.1) inhibits the spatial spread of wavefunctions — which will be helpful for numerics — and leads to minima of the potential at

[TABLE]

In particular, one can have $X^{i}=\nu J^{i}$ with the $J^{i}$ being, for example, the $N$ dimensional irreducible representation of the $\mathfrak{su}(2)$ algebra. This set of matrices defines a ‘fuzzy sphere’ [36]. There are two important features of this solution. Firstly, in the large $N$ limit the noncommutative algebra generated by the $X^{i}$ approaches the commutative algebra of functions on a smooth two dimensional sphere [37, 38]. Secondly, the large $\nu$ limit is a semiclassical limit in which the classical fuzzy sphere solution accurately describes the quantum state. In this semiclassical limit, the low energy excitations above the fuzzy sphere state are obtained from classical harmonic perturbations of the matrices about the fuzzy sphere [39]. See also [40] for an analogous study of the large-mass BMN theory. At large $N$ and $\nu$, these excitations describe fields propagating on an emergent spatial geometry.

By using variational Monte Carlo with generative flows we will obtain a fully quantum mechanical description of this emergent space. This, in itself, is excessive given that the physics of the fuzzy sphere is accessible to semiclassical computations. Our variational wavefunctions will quantitatively reproduce the semiclassical results in the large $\nu$ limit, thereby providing a solid starting point for extending the variational method across the entire $N$ and $\nu$ phase diagram. Exploring the parameter space, we find that the fuzzy sphere collapses upon moving into the small $\nu$, quantum regime. We will consider two different ‘sectors’ of the model, with different fermion number $R$. The first will be purely bosonic states, with $R=0$. The second will have a $R=N^{2}-N$. In this latter sector, the fuzzy sphere state is supersymmetric at large positive $\nu$, so we refer to this as the ‘supersymmetric sector’. In the bosonic sector of the model the fuzzy sphere is a metastable state, and collapses in a first order large $N$ transition at $\nu\sim\nu_{\text{c}}\approx 4$. See Figs. 2 and 3 below. In the supersymmetric sector of the model, where the fuzzy sphere is stable, the collapse is found to be more gradual. See Figs. 6 and 7. In Fig. 8 we start to explore the small $\nu$ limit of the supersymmetric sector.

Beyond the energetics of the fuzzy sphere state, we will define a factorization of the microscopic quantum mechanical Hilbert space that leads to a boundary-law entanglement entropy at large $\nu$. See (5.14) below. This factorization at once captures the emergent local dynamics of fields on the fuzzy sphere and also reveals a microscopic cutoff to this dynamics at a scale set by $N$. The nature of the emergent fields and their cutoff can be usefully discussed in string theory realizations of the model. In string-theoretic constructions, fuzzy spheres arise from the polarization of D branes in background fields [41, 42, 43, 44]. A matrix quantum mechanics theory such as (1.1) describes $N$ ‘D0 branes’ — see [32] and the discussion section below for a more precise characterization of the string theory embedding of mini-BMN theory — and the maximal fuzzy sphere corresponds to a configuration in which the D0 branes polarize into a single spherical D2 brane. There is no gravity associated to this emergent space, the emergent fields describe the low energy worldvolume dynamics of the D2 brane. In this case, the emergent fields are a Maxwell field and a single scalar field corresponding to transverse fluctuations of the brane. In the final section of the paper we will discuss how richer, gravitating states may arise in the opposite small $\nu$ limit of the model.

2 The mini-BMN model

The mini-BMN Hamiltonian is [32]

[TABLE]

The bosonic part $H_{B}$ is given in (1.1). The $\sigma^{k}$ are Pauli matrices. The $\lambda$ are matrices of two-component SO(3) spinors. It can be useful to write the matrices in terms of the $\mathfrak{su}(N)$ generators $T_{A}$, with $A=1,2,\ldots,N^{2}-1$, which obey $[T_{A},T_{B}]=if_{ABC}T_{C}$ and are Hermitian and orthonormal (with respect to the Killing form). That is, $X^{i}=X^{i}_{A}T^{A}$ and $\lambda^{\alpha}=\lambda^{\alpha}_{A}T^{A}$.111The $ijk$ and $ABC$ indices are freely raised and lowered. Lower $\alpha\beta$ indices are for spinors transforming in the $\mathbf{2}$ representation of SO(3), while upper indices are for $\bar{\mathbf{2}}$. We will not raise or lower spinor indices. The full Hamiltonian can then be written

[TABLE]

where $\lambda_{A}^{\alpha\dagger}\equiv(\lambda_{A\alpha})^{\dagger}$ and $\{\lambda_{A}^{\alpha\dagger},\lambda_{B\beta}\}=\delta_{AB}\delta^{\alpha}_{\beta}$ are complex fermion creation and annihilation operators. This Hamiltonian is seen to have four supercharges

[TABLE]

that obey

[TABLE]

States that are invariant under all supercharges therefore have vanishing energy.

Matrix quantum mechanics theories arising from microscopic string theory constructions are typically gauged. This means that physical states must be invariant under the $\mathrm{SU}(N)$ symmetry. In particular, physical state are annihilated by the generators

[TABLE]

2.1 Representation of the fermion wavefunction

The mini-BMN wavefunction can be represented as a function from bosonic matrix coordinates to fermionic states $\psi(X)=f(X)|M(X)\rangle$. Here $X$ denotes the three bosonic traceless Hermitian matrices. The function $f(X)\geq 0$ is the norm of the wavefunction at $X$ while $|M(X)\rangle$ is a normalized state of matrix fermions. A fermionic state with definite fermion number $R$ is parametrized by a complex tensor $M^{ra}_{A\alpha}$ such that

[TABLE]

where $|0\rangle$ is the state with all fermionic modes unoccupied.

The definition (2.6) is parsed as follows: for any fixed $r$ and $a$, $\eta^{ra\dagger}=\sum_{\alpha A}M^{ra}_{A\alpha}\lambda_{A}^{\alpha\dagger}$ is the creation operator for the matrix fermionic modes, where $A$ runs over some orthonormal basis of the $\mathfrak{su}(N)$ Lie algebra and $\alpha=1,2$ for two fermionic matrices. Then $\prod_{a}\eta^{ra\dagger}|0\rangle$ is a state of multiple free fermions created by $\eta^{\dagger}$. The final summation over $r$ in (2.6) is a decomposition of a general fermionic state into a sum of free fermion states. Such a representation is seen to be completely general (but not unique) if we have the number of free fermion states $D$ sufficiently large.

For purely bosonic models, $|M(X)\rangle$ is simply the phase of the wavefunction.

2.2 Gauge invariance and gauge fixing

The generators (2.5) correspond to the following action of an element $U\in G=\mathrm{SU}(N)$ on the wavefunction:

[TABLE]

that is, the group acts by matrix conjugation. The wavefunction is required to be invariant under the group action, i.e. $U\psi=\psi$ for any $U\in G$.

Gauge invariance allows us to evaluate the wavefunction using a representative for each orbit of the gauge group. Let $\widetilde{X}$ be the representative in the gauge orbit of $X$. Gauge invariance of the wavefunction implies that there must exist functions $\widetilde{f}$ and $\widetilde{M}$ such that

[TABLE]

The functions $\widetilde{f}$ and $\widetilde{M}$ take gauge representatives as inputs, or may be thought as gauge invariant functions. The wavefunction we use will be in the form (2.8). The functions $\widetilde{f}$ and $\widetilde{M}$ will be parametrized by neural networks, as we describe in the following section 3.

We proceed to describe the gauge fixing we use to select the representative for each orbit, as well as the measure factor associated with this choice. The $\mathrm{SU}(N)$ gauge representative $\widetilde{X}$ will be such that

$X^{i}=U\widetilde{X}^{i}U^{-1}$ for $i=1,2,3$ and some unitary matrix $U$. 2. 2.

$\widetilde{X}^{1}$ is diagonal and $\widetilde{X}^{1}_{11}\leq\widetilde{X}^{1}_{22}\leq\ldots\leq\widetilde{X}^{1}_{NN}$. 3. 3.

$\widetilde{X}^{2}_{i(i+1)}$ is purely imaginary with the imaginary part positive for $i=1,2,\ldots,N-1$.

The third condition is needed to fix the $\mathrm{U}(1)^{N-1}$ residual gauge freedom after diagonalizing $X^{1}$. The representative $\widetilde{X}$ is well-defined except on a subspace of measure zero where the matrices are degenerate. Then $\widetilde{X}$ can be represented as a vector in $\mathbb{R}^{2(N^{2}-1)}$ with a positivity constraint on some components. The change of variables from $X$ to $\widetilde{X}$ leads to a measure factor given by the volume of the gauge orbit:

[TABLE]

with

[TABLE]

Keeping track of this measure (apart from an overall prefactor) will be important for proper sampling in the Monte Carlo algorithm. The derivation of (2.10) is shown in Appendix A.

3 Architecture design for matrix quantum mechanics

In this work we propose a variational Monte Carlo method with importance sampling to approximate the ground state of matrix quantum mechanics theories, leading to an upper bound on the ground state energy. The importance sampling is implemented with generative flows. The basic workflow is sketched as follows:

Start with a wavefunction $\psi_{\theta}$ with variational parameters $\theta$. In our case $\theta$ will characterize neural networks. 2. 2.

Write the expectation value of the Hamiltonian to be minimized as

[TABLE]

In the mini-BMN case $X$ denotes three traceless Hermitian matrices (indices omitted) and $H_{X}[\psi_{\theta}]$ is the energy density at $X$. Notationally $\mathbb{E}_{X\sim p(X)}$ is the expectation value, with the random variable $X$ drawn from the probability distribution $p(X)$. 3. 3.

Generate random samples according to the wavefunction probabilities $X\sim p_{\theta}(X)=|\psi_{\theta}(X)|^{2}$, and evaluate their energy densities $H_{X}[\psi_{\theta}]$. The variational energy (3.1) can then be estimated as the average of energy densities of the samples. 4. 4.

Update the parameters $\theta$ (via stochastic gradient descent) to minimize $E_{\theta}$:

[TABLE]

where $t=1,2,\ldots$ denotes the steps of training and the parameter $\alpha>0$ sets the learning rate. The gradient of energy is estimated from Monte Carlo samples:

[TABLE]

The method is applicable even if the probabilities are available only up to an unknown normalization factor. 5. 5.

Repeat steps 3 and 4 until $E_{\theta}$ converges. Observables of physical interest are evaluated with respect to the optimal parameters after training.

In the following we discuss details of parametrizing and sampling from gauge invariant wavefunctions with fermions. Technicalities concerning the evaluation of $H_{X}[\psi_{\theta}]$ are spelled out in Appendix B. More details concerning the training are given in Appendix D. Benchmarks are presented at the end of this section.

3.1 Parametrizing and sampling the gauge invariant wavefunction

We first describe how gauge invariance is incorporated into the variational Monte Carlo algorithm. As just discussed, an important step is to sample according to $X\sim|\psi(X)|^{2}$. From (2.8), for a gauge invariant wavefunction $|\psi(X)|^{2}=|\widetilde{f}(\widetilde{X})|^{2}$. However, in sampling $\widetilde{X}$ we must keep track of the measure factor $\Delta(\widetilde{X})$ in (2.10). This is done as follows:

Sample $\widetilde{X}$ according to $p(\widetilde{X})=\Delta(\widetilde{X})|\widetilde{f}(\widetilde{X})|^{2}$. 2. 2.

Generate Haar random elements $U\in\mathrm{SU}(N)$. 3. 3.

Output samples $X=U\widetilde{X}U^{-1}$.

The correctness of this procedure is shown in Appendix A.

Conversely at the evaluation stage, $\psi(X)$ can be computed in the following steps for gauge invariant wavefunctions (2.8):

Gauge fix $X=U\widetilde{X}U^{-1}$ as discussed in the last section. 2. 2.

Compute $\widetilde{M}(\widetilde{X})$ and $\widetilde{f}(\widetilde{X})$. Details of the structure of $\widetilde{M}$ and $\widetilde{f}$ will be discussed below. 3. 3.

Return $\psi(X)=\widetilde{f}(\widetilde{X})|U\widetilde{M}(\widetilde{X})U^{-1}\rangle$ according to (2.8).

We now describe the implementation of $\widetilde{M}$ and $\widetilde{f}$ as neural networks. The basic building block, a multilayer fully-connected (also called dense) neural network, is an elemental architecture capable of parametrizing complicated functions efficiently [12]. The neural network defines a function $F:x\mapsto y$ mapping an input vector $x$ to an output vector $y$ via a sequence of affine and nonlinear transformations:

[TABLE]

Here $A^{1}_{\theta}(x)=M^{1}_{\theta}x+b^{1}_{\theta}$ is an affine transformation, where the weights $M^{1}_{\theta}$ and the biases $b^{1}_{\theta}$ are trainable parameters. The hyperbolic tangent nonlinearity then acts elementwise on $A^{1}_{\theta}(x)$.222We experimented with different activation functions; the final result is not sensitive to this choice. Similar mappings are applied $m$ times, allowing $M^{i}_{\theta}$ and $b^{i}_{\theta}$ to be different for different layers $i$, to produce the output vector $y$. The mapping $F:x\mapsto y$ is nonlinear and capable of approximating any square integrable function if the number of layers and the dimensions of the affine transformations are sufficiently large [45].

The function $\widetilde{M}(\widetilde{X})$ is implemented as such a multilayer fully-connected neural network, mapping from vectorized $\widetilde{X}$ to $\widetilde{M}$ in (2.6), i.e., $\mathbb{R}^{2(N^{2}-1)}\to\mathbb{R}^{DR\,2(N^{2}-1)}$. The implementation of $\widetilde{f}(\widetilde{X})$ is more interesting, as both evaluating $\widetilde{f}(\widetilde{X})$ and sampling from the distribution $p(\widetilde{X})=\Delta(\widetilde{X})|\widetilde{f}(\widetilde{X})|^{2}$ are necessary for the Monte Carlo algorithm. Generative flows are powerful tools to efficiently parameterize and sample from complicated probability distributions. The function $\widetilde{f}(\widetilde{X})=\sqrt{p(\widetilde{X})/\Delta(\widetilde{X})}$, so we can focus on sampling and evaluating $p(\widetilde{X})$, which will be implemented by generative flows.

Two generative flow architectures are implemented for comparison: a normalizing flow and a masked autoregressive flow. The normalizing flow starts with a product of simple univariate probability distributions $p(x)=p_{1}(x_{1})\ldots p_{M}(x_{M})$, where the $p_{i}$ can be different. Values of $x$ sampled from this distribution are passed through an invertible multilayer dense network as in (3.4). The probability distribution of the output $y$ is then

[TABLE]

The masked autoregressive flow generates samples progressively. It requires an ordering of the components of the input, say $x_{1},x_{2},\ldots,x_{M}$. Each component is drawn from a parametrized distribution $p_{i}(x_{i};F_{i}(x_{1},\ldots,x_{i-1}))$, where the parameter depends only on previous components. Thus $x_{1}$ is sampled independently and for other components, the dependence $F_{i}$ is given by (3.4). The overall probability is the product

[TABLE]

When $p_{i}(x_{i})$ are chosen as normal distributions, both flows are able to represent any multivariate normal distribution exactly. Features of the wavefunction (such as polynomial or exponential tails) can be probed by experimenting with different base distributions $p_{i}(x_{i})$. Choices of the base distributions and performances of the two flows are assessed in the following benchmark subsection and also in Appendix D. We will use both types of flow in the numerical results of section 4.

3.2 Benchmarking the architecture

In [34] the Schrödinger equation for the $N=2$ mini-BMN model was solved numerically. Comparison with the results in that paper will allow us to benchmark our architecture, before moving to larger values of $N$. In [34] the Schrödinger equation is solved in sectors with a fixed fermion number

[TABLE]

and total SO(3) angular momentum $j=0,1/2$. We do not constrain $j$, but do fix the number of fermions in the variational wavefunction.

The variational energies obtained from our machine learning architecture with $R=0$ and $R=2$ are shown as a function of $\nu$ in Fig. 1. We take negative $\nu$ to compare with the results given in [34], which uses an opposite sign convention.333There is a particle-hole symmetry of the Hamiltonian (2) via $\nu\to-\nu$, $\lambda\to\lambda^{\dagger}$, $\lambda^{\dagger}\to\lambda$ and $X\to-X$. The masked autoregressive flow yields better (lower) variational energies. These energies are seen to be close to the $j=0$ results obtained in [34]. The variational results seem to be asymptotically accurate as $|\nu|\to\infty$, while remaining a reasonably good approximation at small $\nu$. Small $\nu$ is an intrinsically more difficult regime, as the potential develops flat directions (visualized in [34]) and hence the wavefunction is more complicated, possibly with long tails. In the ‘supersymmetric’ $R=2$ sector, where quantum mechanical effects at small $\nu$ are expected to be strongest, further significant improvement at the smallest values of $\nu$ is seen with deeper autoregressive networks and more flexible base distributions, as we describe shortly. Analogous improvements in these regimes will also be seen at larger $N$ in Sec. 4.3 and Appendix D.

In Fig. 1 the base distributions $p_{i}(x_{i})$, introduced in the previous subsection, are chosen to be a mixture of $s$ generalized normal distributions:

[TABLE]

Here the $k^{i}_{r}$ are positive weights for each generalized normal distribution in the mixture. In (3.8) the $k^{i}_{r}$, $\alpha^{i}_{r}$, $\beta^{i}_{r}$ and $\mu^{i}_{r}$ are learnable (i.e. variational) parameters. For autoregressive flows these parameters further depend on $x_{j}$, with $1\leq j<i$, according to (3.4).

Due to the gauge fixing conditions 2 and 3 in section 2.2, some components $x_{i}$ are constrained to be positive. In the normalization flow this is implemented by an additional map $x_{i}\mapsto\exp(x_{i})$. For the autoregressive flows we have a more refined control over the base distributions; in this case, for components $x_{i}$ that must be positive, we draw from Gamma distributions instead:

[TABLE]

Where again the $k^{i}_{r}$, $\alpha^{i}_{r}$ and $\beta^{i}_{r}$ depend on $x_{j}$, with $1\leq j<i$, according to (3.4).

In Fig. 1 we have shown mixtures with $s=1,3,5$ distributions. The number of layers in (3.4) has been increased with $s$ to search for potential improvements in the space of variational wavefunctions. As noted, the only improvement within the autoregressive flows in going beyond one layer and one generalized normal distribution is seen at the smallest values of $\nu$ with $R=2$. On the other hand, the gap between the variational energies of the two types of flows in Fig. 1 suggests that the wavefunction is complicated in this regime, so that the more sophisticated MAF architecture shows an advantage. The recursive nature of the MAF flows means that they are already ‘deep’ with only a single layer. The complexity of the small $\nu$ wavefunction should be contrasted with the fuzzy sphere phase at large positive $\nu$ discussed in the following section 4 and shown in e.g. Figs. 2 and 3 below. The wavefunction in this semiclassical regime is almost Gaussian, and indeed the NF(1, 1) and MAF(1, 1) flows give similar energies when initialized near fuzzy sphere configurations. The NF architecture in fact gives slightly lower energies in this regime, so we have used normalizing flows in Figs. 2 and 3 for the fuzzy sphere.

The numerics above and below are performed with $D=4$ in (2.6), so that the fermionic wavefunction $|M(X)\rangle$ is a sum of four free fermion states for each value of the bosonic coordinates $X$. In Appendix D we see that increasing $D$ above one lowers the variational energy at small $\nu$, indicating that the fermionic states are not Hartree-Fock in this regime.

4 The emergence of geometry

4.1 Numerical results, bosonic sector

The architecture described above gives a variational wavefunction for low energy states of the mini-BMN model. With the wavefunction in hand, we can evaluate observables. We will start with the purely bosonic sector of the model (i.e. $R=0$). Then we will add fermions. An important difference between the bosonic and supersymmetric cases will be that the semiclassical fuzzy sphere state is metastable in the bosonic theory but stable in the supersymmetric theory.

Figure 2 shows the expectation value of the radius

[TABLE]

for runs initialized close to a fuzzy sphere configuration (solid) and close to zero (open). For large $\nu$ a fuzzy sphere state with large radius is found, in addition to a ‘collapsed’ state without significant spatial extent. Below $\nu_{\text{c}}\approx 4$, the fuzzy sphere state ceases to exist.

The nature of the transition at $\nu_{\text{c}}$ can be understood from the variational energy of the states, plotted in Figure 3. The bosonic semiclassical fuzzy sphere state is seen to be metastable at large $\nu$, as the collapsed state has lower energy. For $\nu<\nu_{\text{c}}$ the fuzzy sphere is no longer even metastable. We will gain a semiclassical understanding of this transition in section 4.2 shortly.

Figures 2 and 3 show that the radius and energy of the fuzzy sphere state are accurately described by semiclassical formulae (derived in the following section) for all $\nu>\nu_{\text{c}}$. In particular this means that $E/N^{3}$ and $r/N$ are rapidly converging towards their large $N$ values. Figure 4 further shows that the probability distribution for the radius $r$ becomes strongly peaked about its semiclassical expectation value at large $\nu$.

Analogous behavior to that shown in Figures 2 and 3 has previously been seen in classical Monte Carlo simulations of a thermal analogue of our quantum transition [46, 47, 48]. These papers study the thermal partition function of models similar to (1.1) in the classical limit, i.e. without the $\Pi^{2}$ kinetic energy term. The fuzzy geometry emerges in a first order phase transition as a low temperature phase in these models. We will see that in our quantum mechanical context the geometric phase is associated with the presence of a specific boundary-law entanglement.

4.2 Semiclassical analysis of the fuzzy sphere

The results above describe the emergence of a (metastable) geometric fuzzy sphere state at $\nu>\nu_{\text{c}}$. In this section we recall that in the $\nu\to\infty$ limit the fluctuations of the geometry are classical fields. For finite $\nu>\nu_{\text{c}}$ the background geometry is well-defined at large $N$, but fluctuations will be described by an interacting (noncommutative) quantum field theory.

In the large $\nu$ limit, the wavefunction can be described semiclassically [39, 40]. We will now briefly review this limit, with details given in the Appendix C. These results provide a further useful check on the numerics, and will guide our discussion of entanglement in the following section 5.

The minima of the classical potential occur at:

[TABLE]

These are supersymmetric solutions of the classical theory, annihilated by the supercharges (2.3) in the classical limit, and therefore have vanishing energy. The solutions of equations (4.2) are

[TABLE]

where the $J^{i}$ are representations of the $\mathfrak{su}(2)$ algebra, $[J^{i},J^{j}]=i\epsilon^{ijk}J^{k}$. We will be interested here in maximal, $N$-dimensional irreducible representations. (Reducible representations can also be studied, corresponding to multiple polarized D branes.)

The $\mathfrak{su}(2)$ Casimir operator suggests a notion of ‘radius’ given by

[TABLE]

Indeed, the algebra generated by the $X^{i}$ matrices tends towards the algebra of functions on a sphere as $N\to\infty$ [37, 38]. At finite $N$, a basis for this space of matrices is provided by the matrix spherical harmonics $\hat{Y}_{jm}$. These obey

[TABLE]

We construct the $\hat{Y}_{jm}$ explicitly in Appendix C. The $j$ index is restricted to $0\leq j\leq j_{\text{max}}=N-1$. The space of matrices therefore defines a regularized or ‘fuzzy’ sphere [36].

Matrix spherical harmonics are useful for parametrizing fluctuations about the classical state (4.3). Writing

[TABLE]

the classical equations of motion can be perturbed about the fuzzy sphere background to give linear equations for the parameters $y^{i}_{jm}$. The solutions of these equations define the classical normal modes. We find the normal modes in Appendix C, proceeding as in [39, 40]. The normal mode frequencies are found to be $\nu\omega$ with

[TABLE]

Recall that $1\leq j\leq j_{\text{max}}=N-1$. The three different sets of frequencies in (4.7) correspond to the group theoretic $\mathfrak{su}(2)$ decomposition $j\otimes 1=(j-1)\oplus j\oplus(j+1)$. Here $j$ is the ‘orbital’ angular momentum and the $1$ is due to the vector nature of the $X^{i}$. We will give a field theoretic interpretation of these modes shortly. The modes give the following semiclassical contribution to the energy of the fuzzy sphere state

[TABLE]

This energy is shown in Figure 3. The scaling as $N^{3}$ arises because there are $N^{2}$ oscillators, with maximal frequency of order $N$. This semiclassical contribution will be cancelled out in the supersymmetric sector studied in section 4.3 below.

The normal modes (4.7) can be understood by mapping the matrix quantum mechanics Hamiltonian onto a noncommutative gauge theory. The analogous mapping for the classical model has been discussed in [49]. We carry out this map in Appendix C. The original Hamiltonian (1.1) becomes the following noncommutative $\mathrm{U}(1)$ gauge theory on a unit spatial $S^{2}$ (setting the sphere radius to one in the field theory description will connect easily to the quantized modes in (4.7)):

[TABLE]

The noncommutative star product $\star$ is defined in the Appendix and

[TABLE]

where the derivatives generate rotations on the sphere $L^{i}=-i\epsilon_{ijk}x^{j}\partial_{k}$ and $[f,g]_{\star}\equiv f\star g-g\star f$. In (4.9) and (4.10) the vector potential $a^{i}$ can be decomposed into two components tangential to the sphere, that become the two dimensional gauge field, and a component transverse to the sphere, that becomes a scalar field. This decomposition is described in Appendix C. The normal modes (4.7) are coupled fluctuations of the gauge field and the transverse scalar field. The zero modes in (4.7) are pure gauge modes, given in (4.11) below. In (4.10) the effective coupling controlling quantum field theoretic interactions is seen to be $1/(N\nu)^{3/2}$. The extra $1/N$ arises because the commutator $[a^{i},a^{j}]_{\star}$ vanishes as $N\to\infty$, see Appendix C. Corrections to the Gaussian fuzzy sphere state are therefore controlled by a different coupling than that of the ‘t Hooft expansion (recall $\lambda=N/\nu^{3}$).

The $\mathrm{SU}(N)$ gauge symmetry generators (2.5) are realized in an interesting way in the non-commutative field theory description. We see in Appendix C that upon mapping to non-commutative fields, the gauge transformations become

[TABLE]

Here $\bm{n}$ is the normal vector and $y(\theta,\phi)$ a local field on the sphere. The first term in (4.11) is the usual $\mathrm{U}(1)$ transformation. The second term describes a coordinate transformation with infinitesimal displacement $\bm{n}\times\nabla y$. Indeed, it is known that non-commutative gauge theories mix internal and spacetime symmetries, which in this case are area-preserving diffeomorphisms of the sphere [50, 51]. The emergent $\mathrm{U}(1)$ non-commutative gauge theory thereby realizes the large $N$ limit of the microscopic $\mathrm{SU}(N)$ gauge symmetry, as area-preserving diffeomorphisms [37, 38].

The fluctuation modes about the fuzzy sphere background allow a one-loop quantum effective potential for the radius to be computed in Appendix C. The potential at $N\to\infty$ is shown in Fig. 5. At large $\nu$ the effective potential shows a metastable minimum at $r\sim N\nu/2$. For $\nu<\nu^{\text{1-loop}}_{\text{c},N=\infty}$ this minimum ceases to exist. The large $N$, one-loop analysis therefore qualitatively reproduces the behavior seen in Figs. 2 and 3. The quantitative disagreement is mainly due to finite $N$ corrections. The transition is only sharp as $N\to\infty$.

4.3 Numerical results, supersymmetric sector

We now consider states with fermion number $R=N^{2}-N$. The fuzzy sphere background is now supersymmetric at large positive $\nu$ [32]. The contribution of the fermions to the ground state energy is seen in Appendix C to cancel the bosonic contribution (4.8) at one loop:

[TABLE]

In Figure 6 the variational upper bound on the energy of the fuzzy sphere state remains close to zero for all values of $\nu$. Figure 7 shows the radius as a function of $\nu$. Probing the smallest values of $\nu$ requires a more powerful wavefunction ansatz than those of Figs. 6 and 7. We will consider that regime shortly.

In contrast to the states with zero fermion number in Figure 3, here the fuzzy sphere is seen to be the stable ground state at large $\nu$. However, the fuzzy sphere appears to merge with the collapsed state below a value of $\nu$ that decreases with $N$. This is physically plausible: while the classical fuzzy sphere radius $r^{2}\sim\nu^{2}N^{2}$ decreases at small $\nu$, quantum fluctuations of the collapsed state are expected to grow in space as $\nu\to 0$. This is because the flat directions in the classical potential of the $\nu=0$ theory, given by commuting matrices, are not lifted in the presence of supersymmetry [52]. Eventually, the fuzzy sphere should be subsumed into these quantum fluctuations. This smoother large $N$ evolution towards small $\nu$ (relative to the bosonic sector) is mirrored in the thermal behavior of classical supersymmetric models [53, 54].

Indeed, exploring the small $\nu$ region with more precision we observe a physically expected feature. In Fig. 8 we see that as $\nu$ decreases towards zero, the radius not only ceases to follow the semiclassical decreasing behavior, but turns around and starts to increase. The variance in the distribution of the radius is also seen to increase towards small $\nu$, revealing the quantum mechanical nature of this regime. These behaviors (non-monotonicity of radius and increasing variance) are expected — and proven for $N=2$ — because the flat directions of the classical potential at $\nu=0$ mean that the extent of the wavefunction is set by purely quantum mechanical effects in this limit.

The small $\nu$ regime here is furthermore an opportunity to test the versatility of our variational ansatz away from semiclassical regimes. In Appendix D we see that for small $\nu$ MAFs achieve much lower energies than NFs. Increasing the number of distributions in the mixture and the number $D$ of free fermions states in (2.6) further lowers the energy. These facts mirror the behavior we found in our $N=2$ benchmarking in Sec. 3.2 at small $\nu$, increasing our confidence in the ability of the network to capture this regime for large $N$ also. The error in a variational ansatz is, as always, not controlled and therefore further exploration of this regime is warranted before very strong conclusions can be drawn. We plan to revisit this regime in future work, to search for the possible presence of emergent ‘throat’ geometries as we discuss in Sec. 6 below.

5 Entanglement on the fuzzy sphere

In this section we will see that the large $\nu$ fuzzy sphere state discussed above contains boundary-law entanglement. To compute the entanglement, one must first define a factorization of the Hilbert space. For our emergent space at finite $N$ and $\nu$ the geometry is both fuzzy and fluctuating, and hence lacks a canonical spatial partition. The fuzziness of the sphere is captured by a toy model of a free field on a sphere with an angular momentum cutoff. Recall from the previous section 4 that the noncommutative nature of the fuzzy sphere amounts to an angular momentum cutoff ${j_{\max}}=N-1$. We will start, then, by defining a partition of the space of functions with such a cutoff.

5.1 Free field with an angular momentum cutoff

Consider a free massive complex scalar field $\varphi(\theta,\phi)$ on a unit two-sphere with the following Hamiltonian:

[TABLE]

Here $\pi$ is the field conjugate to $\varphi$. We impose a cutoff $j\leq{j_{\max}}$ on the angular momentum, rending the quantum mechanical problem well-defined. The fields can therefore be decomposed into a sum of spherical harmonic modes:

[TABLE]

The ‘wavefunctional’ of the quantum field $\varphi(\theta,\phi)$ is then a mapping from coefficients $a_{jm}$ to complex amplitudes. The ground state wavefunctional of the Hamiltonian (5.1) is

[TABLE]

To calculate entanglement for quantum states a factorization of the Hilbert space $\mathcal{H}=\mathcal{H}_{1}\otimes\mathcal{H}_{2}$ is prescribed. To motivate the construction of such a factorization in the fuzzy sphere case, we now review a general framework of defining entanglement in (factorizable) quantum field theories. In quantum mechanics, a quantum state is a function from the configuration space $Q$ to complex numbers, and the Hilbert space of all quantum states is commonly the square integrable functions $\mathcal{H}=L^{2}(Q)$. In quantum field theories, the space $Q$ is furthermore a linear space of functions on some geometric manifold $M$, and thus an orthogonal decomposition $Q=Q_{1}\oplus Q_{2}$ induces a factorization of $\mathcal{H}=L^{2}(Q_{1})\otimes L^{2}(Q_{2})$, which can be exploited to define entanglement.

To define entanglement it then suffices to find an orthogonal decomposition of the space of fields on the fuzzy sphere. Without an angular momentum cutoff, i.e. with ${j_{\max}}\to\infty$, there is a natural choice for any region $A$ on the sphere, which sets $Q_{1}$ to be all functions supported on $A$, and $Q_{2}$ all functions supported on $\bar{A}$, the complement of $A$. Any function $f$ on $M$ can be uniquely written as a sum of $f_{1}\in Q_{1}$ and $f_{2}\in Q_{2}$, where $f_{1}=f\chi_{A}$ and $f_{2}=f(1-\chi_{A})$. Here $\chi_{A}$ is the function on the sphere that is 1 on $A$ and 0 otherwise. Note that the map of multiplication by $\chi_{A}$, $f\mapsto f\chi_{A}$, acts as the projection $Q_{1}\oplus Q_{2}\to Q_{1}$. Conversely, given any orthogonal projection operator $P:Q\to Q$, we can decompose $Q=\mathrm{im}P\oplus\ker P$.

When the cutoff ${j_{\max}}$ is finite, multiplication by $\chi_{A}$ will generally take the function out of the subspace of functions with $j\leq{j_{\max}}$. However, we can still do our best to approximate the projector $P_{A}^{\infty}$ of multiplication by $\chi_{A}$, as defined in the previous paragraph, with a projector $P^{j_{\max}}_{A}$ that lives in the subspace with $j\leq{j_{\max}}$. Formally let $Q^{j_{\max}}$ be the space of functions on the sphere spanned by $Y_{jm}(\theta,\phi)$ with $j\leq{j_{\max}}$. Define the orthogonal projector $P_{A}^{j_{\max}}:Q^{j_{\max}}\to Q^{j_{\max}}$ to minimize the distance $\|P_{A}^{j_{\max}}-P_{A}^{\infty}\|$. The projector $P^{j_{\max}}_{A}$ annihilates all functions in the orthogonal complement of $Q^{j_{\max}}$, when viewed as an operator acting on $Q^{\infty}$. It is convenient to choose $\|\cdot\|$ to be the Frobenius norm, and in Appendix E an explicit formula for $P^{j_{\max}}_{A}$ is obtained.

The projector $P^{j_{\max}}_{A}$ then defines a factorization of the Hilbert space $L^{2}(Q^{j_{\max}})=L^{2}(\mathrm{im}\,P^{j_{\max}}_{A})\otimes L^{2}(\ker P^{j_{\max}}_{A})$ for any region $A$, and entanglement can be evaluated in the usual way. In particular, the second Rényi entropy of a pure state $|\psi\rangle$ on a region $A$ is

[TABLE]

where $x_{A}=Px$ and $x_{\bar{A}}=(I-P)x$ are integrated over $\operatorname{im}P$ and $\ker P$, for $P=P_{A}^{j_{\max}}$, and $x_{A}$ and $x_{\bar{A}}$ can be more compactly combined into a field $x$ with $j\leq{j_{\max}}$. Note that the various $x$’s in (5.4) denote functions on the sphere.

The projector $P^{j_{\max}}_{A}$ is found to have two important geometric features:

The trace of the projector, which counts the number of modes in a region, is proportional to the size of the region. Specifically, at large ${j_{\max}}$, $\operatorname{tr}P^{j_{\max}}_{A}\propto j_{\text{max}}^{2}\left|A\right|$ as is seen numerically in Fig. 9 and understood analytically in Appendix E.

The second Rényi entropy defined by the projector follows a boundary law. At large ${j_{\max}}$, with the mass fixed to $\mu=1$, the entropy $S_{2}\approx 0.03\,{j_{\max}}\left|\partial A\right|$ as is seen numerically in Fig. 10 and understood analytically in Appendix E.

This boundary entanglement law in Fig. 10 is of course precisely the expected entanglement in the ground state of a local quantum field [6, 7]. As the cutoff ${j_{\max}}$ is removed, the entanglement grows unboundedly.

The partition we have just defined can now be adapted to the fluctuations about the large $\nu$ fuzzy sphere state in the matrix quantum mechanics model. We do this in the following subsection. Intuitively, we would like to replace the $j(j+1)+\mu^{2}$ spectrum of the free field in the wavefunction (5.3) with the matrix mechanics modes (4.7). Recall that the matrix modes are cut off at angular momentum $j_{\text{max}}=N-1$.

5.2 Fuzzy sphere in the mini-BMN model

Now we address two additional subtleties that arise when adapting the free field ideas above to the mini-BMN fuzzy sphere. Firstly, the mini-BMN theory is an $\mathrm{SU}(N)$ gauge theory. It is known that entanglement in gauge theories may depend upon the choice of gauge-invariant algebras associated to spatial regions [55]. Different prescriptions correspond to different boundary or gauge conditions [56]. However for a fuzzy geometry, the boundaries of regions and gauge edge modes are not sharply defined. To introduce the fewest additional degrees of freedom, we choose to factorize the physical Hilbert space, instead of an extended one [57, 58], to evaluate entanglement in the mini-BMN model. This is similar to the ‘balanced center’ procedure in [55], where edge modes are absent.444It should, nonetheless, be possible to identify meaningful $\mathrm{SU}(N)$ ‘edge modes’ that would reproduce the edge mode contribution of the emergent Maxwell field. This is an especially interesting question in the light of the fact that the microscopic $\mathrm{SU}(N)$ gauge symmetry also acts as an area-preserving diffeomorphism on the emergent fields in (4.11). This is left for future work.

Secondly, the emergent fields include fluctuations of the geometry itself. The factorization that we have discussed in the previous subsection is tailored to a region on the sphere, and does not need to approximate a spatial region in other geometries. The partition is even less meaningful in non-geometric regions of the Hilbert space. The variational wavefunction we have constructed can be used to compute entanglement for any given factorization of the Hilbert space, but it is unclear that preferred factorizations exist away from geometric limits. In this work we will focus on the entanglement in the $\nu\to\infty$ limit where the fields are infinitesimal, and hence do not backreact on the spherical geometry. In this limit the factorization is precisely — up to issues of gauge invariance — that of the free-field case discussed in the previous subsection.

The matrices corresponding to the infinitesimal fields on the fuzzy sphere are, cf. (4.6),

[TABLE]

which should be thought of as living in the tangent space at $X^{i}=\nu J^{i}$. At large $\nu$ the wavefunction is strongly supported on the classical configuration and hence in this limit the infinitesimal description is accurate. Gauge transformations then act as

[TABLE]

where $\epsilon$ is infinitesimal and $Y$ is an arbitrary Hermitian matrix. The $\epsilon[Y,A^{i}]$ term is omitted in (5.6) as it is of higher order. Gauge invariance of the state is manifested as

[TABLE]

Physical states are wavefunctions on gauge orbits $[A^{i}]$, the set of infinitesimal matrices differing from $A^{i}$ by a gauge transformation (5.6). Similarly to the discussion of free fields above, a partition of the space of gauge orbits is specified by a projector $P$. We will now explain how this projector is constructed. Given a projector $P^{\prime}$ acting on infinitesimal matrices $A^{i}$, a projector acting on gauge orbits can be defined as

[TABLE]

However, for $P$ to be well-defined, $P^{\prime}$ must preserve gauge directions:

[TABLE]

for any $A^{i}$, $Y$ and some $Y^{\prime}$ dependent on $Y$. Let $V$ be the subspace of gauge directions:

[TABLE]

then (5.9) is equivalent to the requirement that $P^{\prime}(V)\subset V$. The strategy for finding the projector $P$ is to solve for the projector $P^{\prime}$ that minimizes $\|P^{\prime}-\chi_{A}\|$ subject to the constraint that (5.9) is satisfied. Then $P$ is defined via $P^{\prime}$ as in (5.8).

The problem of minimizing $\|P^{\prime}-\chi_{A}\|$ for orthogonal projectors $P^{\prime}$ such that $P^{\prime}(V)\subset V$ is exactly solvable as follows. The condition that $P^{\prime}(V)\subset V$ is equivalent to imposing that $P^{\prime}=P_{V}\oplus P_{V_{\perp}}$, where $P_{V}$ is some projector in the subspace $V$ and $P_{V_{\perp}}$ in its orthogonal complement $V_{\perp}$. And $\|P^{\prime}-\chi_{A}\|$ is minimized if and only if $\|P_{V}-\left.\chi_{A}\right|_{V}\|$ and $\|P_{V_{\perp}}-\left.\chi_{A}\right|_{V_{\perp}}\|$ are both minimized. Via the correspondence between matrix spherical harmonics $\hat{Y}_{jm}$ and spherical harmonic functions $Y_{jm}(\theta,\phi)$ in Appendix C, both of these minimizations become the same problem as in the free field case, with a detailed solution in Appendix E.

The second Rényi entropy, in terms of gauge orbits, is evaluated similarly to (5.4):

[TABLE]

where $\Delta$ are measure factors for gauge orbits and $\psi_{\text{inv}}([A])=\psi(\nu J+A)$. Recall that $\psi$ is gauge invariant according to (5.7). The formula (5.2) as displayed does not involve any gauge choice. However, there are some gauges where evaluating (5.2) is particularly convenient. The gauge we choose for this purpose, which is different from that in section 2.2, is that $A\in V_{\perp}$, i.e., the fields are perpendicular to gauge directions. In this gauge measure factors are trivial and the projector is simply $P_{V_{\perp}}$ that minimizes $\|P_{V_{\perp}}-\left.\chi_{A}\right|_{V_{\perp}}\|$:

[TABLE]

where $\psi_{\perp}(A)$ is defined as $\psi(\nu J+A)$ for $A\in V_{\perp}$.555We can find a gauge transformation $U\in\mathrm{SU}(N)$ mapping any matrices $X^{i}$ into this perpendicular gauge as follows. We are looking for $\widetilde{X}^{i}=UX^{i}U^{-1}$, such that $\widetilde{X}^{i}-\nu J^{i}\in V_{\perp}$. This means that $\sum_{i}\operatorname{tr}\left([Y,J^{i}]^{\dagger}(\widetilde{X}^{i}-\nu J^{i})\right)=0$ for any Hermitian matrix $Y$. Equivalently, $\sum_{i}\operatorname{tr}\left(J^{i}[Y,\widetilde{X}^{i}]\right)=0$ for any $Y$. This is achieved by numerically finding the $U$ that maximizes the overlap $\sum_{i}\operatorname{tr}\left(J^{i}UX^{i}U^{-1}\right)$.

The bosonic fuzzy sphere wavefunction can be written in the $\nu\to\infty$ limit as follows. As in (4.6), the perturbations can be decomposed as $A^{i}=\sum_{a}\delta x_{a}\sum_{jm}y^{i}_{jma}\hat{Y}_{jm}\,,$ where the $y^{i}_{jma}$ diagonalize the potential energy at quadratic order in $A$ so that $V=\frac{\nu^{2}}{2}\sum_{a}\omega_{a}^{2}(\delta x_{a})^{2}+\cdots$ (see Appendix C). The wavefunction is then, analogously to (5.3),

[TABLE]

The frequencies are given by (4.7), excluding the pure gauge zero modes. Using this wavefunction, the Rényi entropy (5.2) can be computed exactly and is shown as a solid line in Fig. 11. As $N\to\infty$ these curves approach a boundary law

[TABLE]

Here $\left|\partial A\right|=2\pi\sin\theta_{A}$ is again the circumference of the spherical cap $A$ (in units where the sphere has radius one, consistent with the field theoretic description in (4.9)). The result (5.14) is the same as that of the toy model in Fig. 10, with $j_{\text{max}}$ now set by the microscopic matrix dynamics to be $N-1$.666A (simpler) instance of entanglement revealing the inherent graininess of a spacetime built from matrices is two dimensional string theory [59, 60]. This regulated boundary-law entanglement underpins the emergent locality on the fuzzy sphere at large $N$ and $\nu$. Recall from the discussion around (4.9) that there are only two emergent fields on the sphere: a Maxwell field and a scalar field. The perpendicular gauge choice we have made translates into the Coulomb gauge for the emergent Maxwell field, cf. the discussion around (4.11) above. The factor of $N$ in (5.14) is due to the microscopic cutoff at a scale $L_{\text{fuzz}}\sim L_{\text{sph}}/N$.

Previous works on the entanglement of a free field on a fuzzy sphere involved similar wavefunctions but a different factorization of the Hilbert space, which was inspired instead by coherent states [61, 62, 63, 64]. Those results did not always produce boundary-law entanglement. Here we see that the UV/IR mixing in noncommutative field theories does not preclude a partition of the large $N$ and large $\nu$ Hilbert space with a boundary-law entanglement.

We can also evaluate the entropy (5.2) using the large $\nu$ variational wavefunctions, without assuming the asymptotic form (5.13). The results are shown as dots in Fig. 11. However, we stress that only the $\nu\to\infty$ limit has a clear physical meaning, where fluctuations are infinitesimal. The variational results are close to the exact values in Fig. 11, showing that the neural network ansatz captures the entanglement structure of these matrix wavefunctions.

The results in this section are for the bosonic fuzzy sphere. The projection we have introduced in order to partition the space of matrices can be extended in a similar, but more involved, way to factorize the fermionic Hilbert space.

6 Discussion

We have seen that neural network variational wavefunctions capture in detail the physics of a semiclassical spherical geometry that emerges in the mini-BMN model (2.1) at large $\nu$. Away from the semiclassical limit, the spherical geometry either abruptly or gradually collapses towards a new state. In Fig. 8 we saw that in the ‘supersymmetric’ sector this new state was characterized by an increase in both the expectation value and quantum mechanical variance of the radius as $\nu\to 0$. To understand the physics of this process, and to start thinking about the nature of the collapsed state as $\nu\to 0$, it is helpful to consider the string theoretic embedding of the model.

The mini-BMN model can be realized in string theory as the description of $N$ D-particles in an AdS4 spacetime. Let us review some aspects of this realization [32]. The parameter

[TABLE]

Here $L_{\text{AdS}}$ is the AdS radius, $L_{\text{s}}$ is the string length and $g_{\text{s}}$ is the string coupling. The proportionality in (6.1) depends on the volume, in units of the string length, of internal cycles wrapped by the branes in the compactification down to AdS4. In particular, the mass of a single D-particle goes like $1/g_{\text{s}}$ times the wrapped internal volume. The strength of the gravitational backreaction of $N$ coincident D-particles is then controlled by $G_{N}\cdot N/g_{\text{s}}$. Here $G_{N}\sim g_{\text{s}}^{2}$ is the four dimensional Newton constant, where we have suppressed a factor of the volume of the compactification manifold. Therefore, if we keep the AdS radius fixed in string units, gravitational backreation becomes important when $g_{\text{s}}N\sim N/\nu^{3}\gtrsim 1$. Up to factors of the volume of compactification cycles, this is equivalent to the statement that the dimensionless ’t Hooft coupling $\lambda=N/\nu^{3}$, introduced below (1.1), becomes large.

For $N/\nu^{3}\lesssim 1$, then, the D-particles can be treated as light probes on the background AdS spacetime. The fuzzy sphere configuration describes a polarization of the D-particles into spherical ‘dual giant gravitons’. From the string theory perspective, this polarization is driven by the 4-form flux $\Omega\sim 1/L_{\text{AdS}}$ supporting the background AdS4 spacetime. Together with the discussion in the previous paragraph on the strength of the gravitational interaction, we can write the heuristic relation $N/\nu^{3}\sim\text{gravity}/\text{flux}$. At large $\nu$ the flux wins out and semiclassical fuzzy spheres can exist, but at small $\nu$ gravitational forces cause the spheres to collapse. The entanglement and emergent locality that we have described in this paper is that of the polarized spheres, whose excitations are described by the usual gauge fields and transverse scalar fields of string theoretic D-branes.

For $N/\nu^{3}\gg 1$ it is possible that the strongly interacting, collapsed D-particles will develop a geometric ‘throat’, in the spirit of the canonical holographic correspondence [1]. It is not well-understood when such a throat would be captured by the mini-BMN matrix quantum mechanics. The variational wavefunctions that we have developed here provide a new window into this problem. In particular, we hope to investigate the small $\nu$ collapsed state in more detail in the future, with the objective of revealing any entanglement associated to emergent local dynamics in the throat spacetime. If the emergent dynamics includes gravity, there are two potentially interesting complications. Firstly, the entanglement of bulk fields may be entwined with entanglement due to the ‘stringy’ degrees of freedom that seem to be manifested in the Bekenstein-Hawking entropy of black holes as well as in the Ryu-Takayanagi formula [65, 66, 67, 68]. Secondly, and perhaps relatedly, it may become crucial to understand the ‘edge mode’ contribution to the entanglement, that we have avoided in our discussion here [69, 70].

More generally, the methods we have developed will be applicable to a wide range of quantum problems of interest in the holographic correspondence. The benefit of the variational neural network approach is direct access to properties of the zero temperature quantum mechanical state. Optimizing the numerical methods and variational ansatz further, and with more computational power, it should not be difficult to work with larger values of $N$. In addition to understanding the emergence of spacetime from first principles, it should also be possible to study, for example, the microstates and dynamics of quantum black holes.

Acknowledgements

It is a pleasure to thank Frederik Denef and Xiaoliang Qi for helpful discussions, Aitor Lewkowycz, Raghu Mahajan and Edward Mazenc for comments on the draft, and Tarek Anous for sharing his code with us. We also thank Zhaoheng Guo and Yang Song for collaboration on a related project. SAH is partially funded by DOE award de-sc0018134. XH is supported by a Stanford Graduate Fellowship. Computational work was performed on the Sherlock cluster at Stanford University, with the TensorFlow code for the project available online.

Appendix A Geometry of the gauge

Gauge invariant sampling

In the procedure of sampling bosonic matrices $X$ according to the wavefunction probability distribution $|\psi(X)|^{2}=|f(X)|^{2}$, it is asserted in the main text that $X\sim|f(X)|^{2}$ if we let $X=U\widetilde{X}U^{-1}$ where $U$ is a Haar random element in $\mathrm{SU}(N)$ and the representative of the gauge orbit $\widetilde{X}\sim\Delta(\widetilde{X})|\widetilde{f}(\widetilde{X})|^{2}$. A proof of this assertion, along with a more precise definition of the gauge orbit measure $\Delta$, is presented here.

To simplify notation, denote $\widetilde{X}\sim\widetilde{p}(\widetilde{X})$. If the random variable $X=U\widetilde{X}U^{-1}$, it follows the probability distribution

[TABLE]

where the integral over $\mathrm{SU}(N)$ is with respect to the normalized Haar measure, and $\delta$ is the Dirac delta distribution. For almost any $X_{0}$, there is a unique gauge representative $\widetilde{X}_{0}$, with a discrete set of $U_{i}\in\mathrm{SU}(N)$ ($i=1,2,\ldots,N$), such that $U_{i}\widetilde{X}_{0}U_{i}^{-1}=X_{0}$. These unitaries differ by an overall phase (powers of $\exp(i2\pi/N)$). Hence

[TABLE]

where $J$ is the Jacobian determinant of the map $(\widetilde{X},U)\mapsto U\widetilde{X}U^{-1}$. As will be seen in the next subsection, $J(\widetilde{X},U)=J(\widetilde{X})$ does not depend on the unitary $U$. So if we assign

[TABLE]

and note $\widetilde{p}(\widetilde{X})=\Delta(\widetilde{X})|\widetilde{f}(\widetilde{X})|^{2}$,

[TABLE]

for a gauge invariant wavefunction (2.8). This is the desired result.

Derivation of the gauge orbit measure

From (A.3), the gauge orbit measure $\Delta$ is given by the Jacobian determinant $J$ of the map $X:(\widetilde{X},U)\mapsto U\widetilde{X}U^{-1}$. Recall that for a general mapping $F$ between smooth manifolds of equal dimension $S\to T$, the Jacobian determinant can be written in terms of the pullback of the volume form

[TABLE]

where $\omega_{S}$ and $\omega_{T}$ are volume forms on $S$ and $T$. That is, $J$ is the ratio of the volume element after and before the mapping. If $x_{i}$ and $y_{i}$ are two orthonormal coordinate systems at $x\in S$ and $y=F(x)\in T$, in terms of the wedge product,

[TABLE]

Therefore equation (A.5) can be expressed more explicitly as

[TABLE]

We would like to show firstly that $J(\widetilde{X},U)$ does not depend on $U$. Note that the map $X:(\widetilde{X},U)\mapsto U\widetilde{X}U^{-1}$ is equivariant with respect to the following actions of $G=\mathrm{SU}(N)$: for any $U^{\prime}\in G$, in the base space $U^{\prime}\cdot(\widetilde{X},U)=(\widetilde{X},U^{\prime}U)$, and in the target space $U^{\prime}\cdot X=U^{\prime}XU^{\prime-1}$. And the two actions preserve the volume forms, because the Haar measure is left invariant and the metric $\operatorname{tr}dX^{\dagger}dX$ is invariant under matrix conjugation. Hence the Jacobian $J(\widetilde{X},U)=J(\widetilde{X})$ is independent of $U$.

We will obtain the Jacobian by explicitly computing the pullback of the volume form at $X$. As the Jacobian does not depend on $U$, it is convenient to evaluate it at $U=I$. To further simplify the computation, we shall complexify the cotangent spaces, which does not change the Jacobian determinant. The $\mathfrak{su}(N)$ real Lie algebra is complexified to $\mathfrak{sl}(N)$, and the following basis $\{D_{i},E_{ij}\}$ of $\mathfrak{sl}(N)$ is employed. The basis is orthonormal with respect to the matrix inner product $\operatorname{tr}X^{\dagger}Y$:

For $1\leq i\leq N-1$, $D_{i}$ is a diagonal matrix with $(D_{i})_{jj}=1/\sqrt{i(i+1)}$ for $1\leq j\leq i$, $(D_{i})_{jj}=-(j-1)/\sqrt{i(i+1)}$ for $j=i+1$ and $(D_{i})_{jj}=0$ for $j>i+1$. 2. 2.

For $1\leq i,j\leq N$ and $i\neq j$, $E_{ij}$ is the matrix that has only one nonzero entry $(E_{ij})_{ij}=1$.

A general element in the complexified cotangent space of $\widetilde{X}$ is (with the gauge choice defined in the main text)

[TABLE]

where the superscript $i=1,2,3$ denotes three bosonic matrices. The equations (A.8) thus define a basis $\{d\widetilde{c}^{1}_{i},d\widetilde{c}^{2}_{i},d\widetilde{e}^{2}_{i(i+1)},d\widetilde{e}^{2}_{ij},d\widetilde{c}^{3}_{i},d\widetilde{e}^{3}_{ij}\}$ of the complexified cotangent space of $\widetilde{X}$.

The complexified cotangent space of $\mathrm{SU}(N)$ at $U=I$ is isomorphic to the Lie algebra $\mathfrak{sl}(N)$, so that (introducing basis forms $dc_{i},de_{ij}$):

[TABLE]

The differential of the map $X:(\widetilde{X},U)\mapsto U\widetilde{X}U^{-1}$ at $U=I$ is

[TABLE]

and the cotangent space of $X$ is complexified to three copies of $\mathfrak{sl}(N)$, so that (introducing basis forms $dc^{k}_{i},de^{k}_{ij}$):

[TABLE]

Substituting (A.8) and (A.9) into (A.10), recalling that $\widetilde{X}^{1}$ is diagonal, and equating the expressions for $dX^{1}$ we have

[TABLE]

Equating the expressions for $dX^{2}$ gives, with terms that drop out of the final result omitted:

[TABLE]

where the expression for $de^{2}_{ij}$ holds for $|i-j|\neq 1$ and the prefactor $i$ in the expressions for $de^{2}_{i(i+1)}$ and $de^{2}_{(i+1)i}$ is the imaginary unit. Subscripts are omitted if that term with any subscript is unimportant, e.g., $de$ means linear combinations of $de_{ij}$ for $1\leq i\neq j\leq N$. Similarly

[TABLE]

The Jacobian determinant $J$ is evaluated as, schematically,

[TABLE]

where $de^{1}_{ij}$ denotes $\bigwedge_{ij}de^{1}_{ij}$ for $1\leq i\neq j\leq N$, for example. Substitution of (A.12), (A.13) and (A.14) into the left-hand side of (A.15) yields a sum of wedge products of differentials. The wedge product is nonzero only if each factor on the right-hand side of (A.15) appears exactly once. Now observe that $de_{ij}$ already appears in $de^{1}_{ij}$ in (A.12), hence all $de_{ij}$ terms in other factors can be safely ignored.

With the $de_{ij}$ ignored, $de^{2}_{i(i+1)}\wedge de^{2}_{(i+1)i}$ is proportional to $dc_{i}\wedge d\widetilde{e}^{2}_{i(i+1)}$ for $i=1$, because for any differential $da$, $da\wedge da=0$. Then remaining factors of $dc_{1}$ and $d\widetilde{e}^{2}_{12}$ can be ignored. Next, for $i=2$, $de^{2}_{i(i+1)}\wedge de^{2}_{(i+1)i}$ must be proportional to $dc_{i}\wedge d\widetilde{e}^{2}_{i(i+1)}$ as well, up to terms that can be ignored. In the end we have (note that $\widetilde{X}^{2}_{i(i+1)}=-\widetilde{X}^{2}_{(i+1)i}$ is purely imaginary)

[TABLE]

Now terms with $dc_{i}$ can be ignored as well as they appear in (A.16). With the $dc_{i}$ and $de_{ij}$ ignored, $dc^{1}_{i}$, $dc^{2}_{i}$, $de^{2}_{ij}$ for $|i-j|\neq 1$, $dc^{3}_{i}$ and $de^{3}_{ij}$ on the left-hand side of (A.15) can be replaced by $d\widetilde{c}^{1}_{i}$, $d\widetilde{c}^{2}_{i}$, $d\widetilde{e}^{2}_{ij}$, $d\widetilde{c}^{3}_{i}$ and $d\widetilde{e}^{3}_{ij}$, respectively, in the light of (A.12), (A.13) and (A.14). The Jacobian is then a product of the factors in (A.12) and (A.16). Thus overall the gauge orbit measure is

[TABLE]

Appendix B Evaluation of observables

The physical observables that we are interested in fall into roughly three categories: (i) bosonic potentials; (ii) fermionic bilinears; (iii) casimirs of Lie group actions. Efficient numerical recipes for evaluating these observables via Monte Carlo simulation are discussed in this Appendix. Monte Carlo requires that the integrals are written as the average over samples $\mathbb{E}_{X\sim|f|^{2}}[\cdot]$.

Bosonic potentials are real functions of bosonic matrix coordinates $V(X)$, and they are straightforward to evaluate:

[TABLE]

Fermionic bilinears and casimirs are more elaborate to compute. The final results are (B.20) and (B.25) with detailed derivations presented below.

Fermionic bilinears

Expectation values of fermionic bilinears $B(\lambda^{\dagger},\lambda,X)$ are

[TABLE]

The problem is thus essentially to evaluate fermionic bilinears in the fermionic state $|M(X)\rangle$, which can furthermore be reduced to calcuating

[TABLE]

where $|M^{r}\rangle$ is the free fermion state

[TABLE]

The question is more generally formulated as follows: let $M$ be a complex matrix of size $R\times P$ and denote its corresponding free fermion state as

[TABLE]

then what are the matrix elements $\langle M^{\prime}|B(\lambda^{\dagger},\lambda,X)|M\rangle$? The starting point is the Slater determinant:

[TABLE]

and note that

[TABLE]

where the first term on the right-hand side can be evaluated from (B.6). The second term in (B.7) can be read as the overlap between free fermion states $\lambda_{q}^{\dagger}|M^{\prime}\rangle$ and $\lambda_{p}^{\dagger}|M\rangle$ and thus (B.6) is again applicable:

[TABLE]

A dummy variable $s$ is introduced for later convenience. Using (B.14) in (B.7)

[TABLE]

Differentiate both sides with respect to $s^{2}$ to obtain a more compact expression:

[TABLE]

where $\operatorname{adj}A=(\det A)A^{-1}$ is the adjucate of $A$. For an arbitrary bilinear $W$,

[TABLE]

Back to the original problem of calculating (B.3). Equation (B.17) is applicable if we regard the index $p$ in (B.5) as running over both the indices $\alpha$ and $A$ in (B.4). Define the overlap matrix

[TABLE]

then

[TABLE]

where the fermionic operators in the bilinear are replaced by complex matrices so that the expression is a complex number. Finally summing over $r$ and $s$,

[TABLE]

Casimirs

The observables discussed above do not involve derivatives. Derivatives show up in kinetic terms, for example, and can be understood in a geometric way. For an action of a Lie group $G$ on the wavefunction $\psi$, a casimir term can be defined as

[TABLE]

where the summation is over an orthonormal basis of the Lie algebra and

[TABLE]

As an example, consider the group of translations of bosonic coordinates $X\to X+\delta X$ that acts on the wavefunction as

[TABLE]

and thus in this case

[TABLE]

which is the usual kinetic term. If $G=\mathrm{SU}(N)$ with the adjoint action on matrices, the observable (B.21) is the casimir of the gauge group, and if $G=\mathrm{SO}(3)$ in the mini-BMN model, the observable measures the angular momentum quantum number of the state.

The summation and the integral in (B.21) are estimated from Monte Carlo samples as:

[TABLE]

where $f=|\psi|$, $|T_{A}|^{2}=\dim G$ means that the expectation value averages over all Lie algebra elements $T_{A}$ with norm $\sqrt{\dim G}$.

Appendix C Semiclassical analysis of the fuzzy sphere

Correspondence between matrices and fields on the emergent sphere

A mapping from any $N$-by-$N$ complex matrix $A$ to a function $f_{A}(\theta,\phi)$ is constructed as follows. The construction is motivated by the following principles: (i) the map $A\mapsto f_{A}(\theta,\phi)$ should be linear; (ii) the map should preserve the inner products:

[TABLE]

Here $\int d\Omega$ is the integral over a $4\pi$ solid angle; (iii) the map should preserve the $\mathfrak{su}(2)$ action:

[TABLE]

As in the main text, the $J^{i}$ are generators of the $N$ dimensional irreducible representation of $\mathfrak{su}(2)$ and the $L^{i}$ are generators for rotations of functions on a sphere:

[TABLE]

and $(x^{1},x^{2},x^{3})=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$.

Requirements (i) and (ii) can be accomplished by mapping an orthonormal basis of matrices to an orthonormal basis of functions on the sphere. In the light of (iii), we choose spherical harmonics $Y_{jm}(\theta,\phi)$ ($j\geq 0$, $|m|\leq j$) as the basis of functions:

[TABLE]

and they are orthonormal with respect to the inner product in (C.1):

[TABLE]

To construct matrix counterparts of spherical harmonics $\hat{Y}_{jm}$, we note that

[TABLE]

where $L^{\pm}=L^{1}\pm iL^{2}$, so (iii) requires (denote $J^{\pm}=J^{1}\pm iJ^{2}$)

[TABLE]

which fixes all the matrices $\hat{Y}_{jm}$ given $\hat{Y}_{j(-j)}$. The $\mathfrak{su}(2)$ representation further requires that $L^{-}Y_{j(-j)}=0$ and $L^{+}Y_{jj}=0$, which translates to the matrix side as $[J^{-},Y_{j(-j)}]=0$ and $[J^{+},Y_{jj}]=0$. Thus for some normalizing factor $C$,

[TABLE]

The matrix $J^{-}$ is nilpotent with order $N$: $(J^{-})^{N}=0$. Therefore the matrices in (C.8) are restricted to $j\leq N-1$. For $j\leq N-1$, the numerical factor $C$ is chosen such that

[TABLE]

The sign of $C$ is not fixed by the three requirements, and we pick $C>0$ in correspondence with spherical harmonics $Y_{j(-j)}\propto(x^{1}-ix^{2})^{j}$.

It is straightforward to verify that

[TABLE]

given the $\mathfrak{su}(2)$ algebra and eqs. (C.7) and (C.8). Hence the matrices $\hat{Y}_{jm}$ form an eigenbasis of adjoint actions of $J^{3}$ and the casimir $(J^{i})^{2}$, and are therefore orthogonal. They are normalized as well because of (C.9). The map $A\mapsto f_{A}(\theta,\phi)$ is then defined on the basis as $\hat{Y}_{jm}\mapsto Y_{jm}(\theta,\phi)$, fulfilling the requirements (i) to (iii).

Under the correspondence $\hat{Y}_{jm}\mapsto Y_{jm}(\theta,\phi)$, $N$-by-$N$ matrices describe fields on a sphere with angular momentum cutoff ${j_{\max}}=N-1$. Furthermore (C.1) connects matrix observables and averages of fields on the emergent sphere. For instance, the classical fuzzy sphere solution sets $X^{i}=\nu J^{i}$, and we would like to interpret $f_{X^{i}}(\theta,\phi)$ as coordinates $x^{i}$ of the point on the sphere at angle $(\theta,\phi)$. Thus according to (C.1), the radius of the emergent sphere (for irreducible representation $J^{i}$) is

[TABLE]

Noncommutative gauge theory on the fuzzy sphere

In the last subsection we have discussed the correspondence between matrix degrees of freedom and fields on the fuzzy sphere. Given that correspondence the matrix Hamiltonian (2) can be cast into a quantum field theory on the sphere. The caveat is that the fields on the sphere are not commutative, due to the noncommutative nature of matrix multiplication.

To be more precise, we define the ‘star product’ of the fields as induced from their corresponding matrix multiplications:

[TABLE]

where $\hat{f}$ and $\hat{g}$ are the matrix counterparts of functions $f(\theta,\phi)$ and $g(\theta,\phi)$ via the correspondence between matrix spherical harmonics and spherical harmonics on the sphere: $\hat{Y}_{jm}\leftrightarrow Y_{jm}(\theta,\phi)$. The prefactor is a result of the normalization (C.9).

The star product is associative but noncommutative. In particular, the commutator of scalar functions may not vanish. For example,

[TABLE]

where $[\cdot,\cdot]_{\star}$ is the commutator with the star product for multiplication. The structure constants $f$ in (C) are known to vanish as $1/N$ as $N\to\infty$ (see, e.g., the Appendix of [49]). The usual commutative product is recovered at $N=\infty$.

To repackage matrix degrees of freedom into emergent fields, expand the bosonic matrices around their classical values:

[TABLE]

where the $A^{i}$ are Hermitian matrices parametrizing fluctuations around the fuzzy sphere. Our re-writing of the Hamiltonian will be exact in $A$. The corresponding emergent fields $\widetilde{a}^{i}(\theta,\phi)$ are as follows:

[TABLE]

The conjugate momenta to the $A^{i}$ are

[TABLE]

obeying the canonical commutation relations $[A^{i}_{ab},(\Pi^{j}_{A})_{cd}]=i\delta^{ij}\delta_{ad}\delta_{bc}$. We will also want to introduce the momenta

[TABLE]

which obey

[TABLE]

The $\widetilde{\pi}^{i}$ therefore become the usual conjugate momenta when ${j_{\max}}=\infty$, where the summation in (C.18) becomes $4\pi\delta(\cos\theta-\cos\theta^{\prime})\delta(\phi-\phi^{\prime})$. Hermiticity of the matrices $A^{i}$ and $\Pi^{i}_{A}$ is manifested as reality of the fields $\widetilde{a}^{i}$ and $\widetilde{\pi}^{i}$.

Substituting (C.15), (C.16) and (C.17) into the matrix Hamiltonian, the kinetic terms are

[TABLE]

The bosonic potential in (1.1) can be written as a square:

[TABLE]

and substituting (C.14) into (C.20):

[TABLE]

The corresponding field is (recall (C.2) and (C.12))

[TABLE]

and the potential can now be written

[TABLE]

The fermionic potential in (2.1) is, in terms of $A^{i}$,

[TABLE]

Let $\widetilde{\psi}(\theta,\phi)$ be the fermionc field corresponding to $\lambda$, then (C.24) is recast into

[TABLE]

where $D^{k}\widetilde{\psi}\equiv iL^{k}\widetilde{\psi}+i\nu^{-1}[\widetilde{a}^{k},\widetilde{\psi}]_{\star}$.

Collect all three parts (C.19), (C.23) and (C.25), and rescale the fields

[TABLE]

The Hamiltonian for the emergent fields, which is equivalent to (2) for matrices, is then

[TABLE]

where

[TABLE]

The $\mathrm{SU}(N)$ gauge symmetry of the matrices leads to the noncommutative $\mathrm{U}(1)$ gauge symmetry of (C). Under an infinitesimal $\mathrm{SU}(N)$ gauge transformation parametrized by a Hermitian matrix $Y$, $\delta X^{i}=i[Y,X^{i}]$, $\delta\lambda^{\alpha}=i[Y,\lambda^{\alpha}]$, and thus by (C.14),

[TABLE]

Let $\widetilde{y}(\theta,\phi)$ be the field corresponding to the matrix $Y$, then the gauge transformation of the noncommutative fields is ($\bm{n}$ is the radial vector and fields should be considered as defined on the unit sphere)

[TABLE]

Recall the rescaling (C.26) and let $\widetilde{y}=y\sqrt{4\pi/N\nu^{3}}$,

[TABLE]

The first term in $\delta a^{i}$ is the usual $\mathrm{U}(1)$ transformation. The second term, which can be obtained from the algebra in (C), describes a coordinate transformation with infinitesimal displacement $\bm{n}\times\nabla y$ [38]. Indeed, it is known that non-commutative gauge theories mix internal and spacetime symmetries, which in this case are area-preserving diffeomorphisms of the sphere [50, 51]. The coordinate transformation in (C.31) is area-preserving because $\nabla\cdot(\bm{n}\times\nabla y)=0$.

In the commutative limit $\nu\to\infty$, the gauge field is decoupled from the fermions and the theory contains a U(1) gauge field on the sphere, with a real massive scalar and a massive Dirac fermion. To see more explicitly the field content of (C.27) in this limit, note that $\bm{L}=-i\bm{n}\times\nabla$ and $f^{ij}=\epsilon^{ijk}\left((\bm{n}\times\nabla)\times\bm{a}+\bm{a}\right)^{k}$ when $\nu\to\infty$ ($\bm{a}$ is the three-dimensional vector notation for $a^{i}$). We then obtain

[TABLE]

The scalar field $\varphi$ is the radial component of the gauge field, and we denote the U(1) gauge field on the sphere as $\bm{b}$:

[TABLE]

The U(1) curvature $f$ of the gauge field $\bm{b}$ defined on the sphere is

[TABLE]

and we have (after some vector calculus manipulations)

[TABLE]

Substituting (C.35) into (C.32), the commutative gauge theory can be rewritten as

[TABLE]

where $\pi^{a}$ and $\pi$ are the conjugate variables of $\bm{b}$ and $\varphi$, respectively, and $\bm{\sigma}$ is the vector of Pauli matrices. The fields in (C.36) should be thought as living on the unit sphere.

Fluctuation spectrum around the classical fuzzy sphere

The classical energy at the fuzzy sphere vanishes due to supersymmetry. In the following we analyze the spectrum of bosonic quadratic fluctuations near the fuzzy sphere configuration, and the spectrum of fermions, as the next order in a semiclassical expansion. The semiclassical correction to energy at this level is shown to be zero as well.

The bosonic potential in (1.1) can be written as a square:

[TABLE]

and quadratic fluctuations around a classical solution are given by

[TABLE]

where $\delta X^{i}=\sum_{a}\delta x_{a}Y^{i}_{a}$ and $Y^{i}_{a}$ are the normalized eigen-matrices:

[TABLE]

Here we specialized to the background solution $X^{j}=\nu J^{j}$.

To solve the eigenvalue equation in (C.39), expand $Y^{i}$ (subscript $a$ omitted) into a sum of matrix spherical harmonics $Y^{i}=\sum_{jm}y^{i}_{jm}\hat{Y}_{jm}$, and note

[TABLE]

For convenience introduce the $\pm$ basis: $y^{\pm}=y^{1}\pm iy^{2}$ and the indices must be raised with $g^{+-}=g^{-+}=2$ and $g^{33}=1$ (other entries are zero). In this basis $\epsilon_{+-3}=i/2$. Then (C.39) can be cast into equations for the coefficients $y^{3}_{jm}$ and $y^{\pm}_{jm}$:

[TABLE]

[TABLE]

Equations (C.41) and (C.42) consist of three linear equations with three variables $y^{3}_{jm}$, $y^{+}_{j(m+1)}$ and $y^{-}_{j(m-1)}$. For there to be nonzero solutions, the determinant must be zero:

[TABLE]

Hence for $0<j<N$, $|m|<j$, the eigenvalues are $\omega=0,-j,j+1$. The edge cases $|m|=j,j+1$ should be treated separately due to the additional constraint $y^{\pm}_{jm}=0$ if $|m|>j$. The eigenvalue equation at $m=\pm j$ is instead $\omega(\omega-j-1)=0$, and for $m=\pm(j+1)$ it is $\omega-j-1=0$.

The multiplicity of the eigenvalue $\omega=0$ is $N^{2}-1$, which accounts for the $\mathrm{SU}(N)$ gauge degrees of freedom. The other eigenvalues are $\omega=-j$ for $1\leq j\leq N-1$ with multiplicity $2j-1$ and $\omega=j+1$ for $1\leq j\leq N-1$ with multiplicity $2j+3$. The ground state energy of the bosonic oscillators (C.38) is therefore

[TABLE]

The spectrum of the fermionic bilinear is found similarly:

[TABLE]

Expand $\lambda^{\alpha}=\sum_{jm}y^{\alpha}_{jm}\hat{Y}_{jm}$ (note now $\alpha=\pm$ labels $\sigma^{3}=\pm 1$ basis). The equations are

[TABLE]

[TABLE]

The eigenvalue equations (C.46) and (C.47) have nontrivial solutions when

[TABLE]

so that for $0<j<N$ and $-j\leq m<j$ there are eigenvalues $\omega=j+3/2$ and $\omega=-j+1/2$. For $m=j$ or $m=-j-1$ the eigenvalue equation is instead $\omega-j-3/2=0$, as $y^{-}_{j(j+1)}=y^{+}_{j(-j-1)}=0$ is imposed.

So the eigenvalues for $0<j<N$ are $\omega=j+3/2$ with multiplicity $2j+2$ and $\omega=-j+1/2$ with multiplicity $2j$. For $\nu>0$ the $\omega=-j+1/2$ modes are occupied with a total number of fermions:

[TABLE]

And the fermionic energy for $\nu>0$ at this order is

[TABLE]

For $\nu<0$ the $\omega=j+3/2$ modes are occupied instead and the number of fermions is

[TABLE]

We see that supersymmetry requires different number of occupied fermions in the case of $\nu>0$ and $\nu<0$. The fermionic energy for $\nu<0$ is

[TABLE]

In either case (C.50) or (C.52) the energy is $-(4N^{3}+5N-9)|\nu|/6$, which exactly cancels the bosonic contribution (C.44). Hence the semiclassical correction to the fuzzy sphere energy is zero at this order, for the specific number of fermions (C.49) or (C.51).

One-loop effective potential and the estimate of $\nu_{\text{c}}$

In the main text we observe a first-order phase transition near $\nu_{\text{c}}\approx 4$ when the bosonic fuzzy sphere phase becomes unstable. Here we give an estimate of $\nu_{\text{c}}$ from the bosonic one-loop effective potential for the radius, at $N=\infty$.

We start with the bosonic potential (C.37) with matrix sources $S_{i}$:

[TABLE]

where the sources $S_{i}(\phi)$ are such that the local energy minimum is at $X^{i}=\phi J^{i}$. The parameter $\phi>0$ is proportional to the radius:

[TABLE]

The classical contribution to the energy (C.53) at $X^{i}=\phi J^{i}$ is

[TABLE]

Quadratic fluctuations of (C.53) around the local minimum give:

[TABLE]

The norm of the spin matrices $J^{i}$ scales as $N$, and hence to leading order in $N$:

[TABLE]

Diagonalizing this leading order piece as we did in the last subsection, the nonzero mode frequencies are now $\omega=-(j+1)\phi$ for $0<j<N$ with multiplicity $2j-1$ and $\omega=j\phi$ for $0<j<N$ with multiplicity $2j+3$. So, the one-loop quantum correction to the ground state energy is

[TABLE]

The one-loop effective potential $\Gamma(\phi)=E_{0}(S_{i}(\phi))+\frac{1}{2}\sum_{a}|\omega_{a}|-\operatorname{tr}S_{i}(\phi)\phi J^{i}$ is then

[TABLE]

where omitted terms are higher order in $N^{-1}$. The critical value of $\nu$ is estimated as when the second order derivative of $\Gamma(\phi)$ at the fuzzy sphere solution vanishes:

[TABLE]

It is clear in (C.59) that, at large $N$, the leading quantum correction to the classical solution is suppressed by $\nu^{-3}$. This shows that the large $\nu$ limit rapidly becomes classical. The critical $\nu_{\text{c}}$ estimated above is at $N=\infty$, where the transistion is sharp.

Appendix D Training and tuning

Training of the model is divided into three epochs, each of which consists of 5000 iterations. The learning rate is set to be $10^{-3}$ for iterations from 1 to 5000, $2\times 10^{-4}$ from 5001 to 10000 and $4\times 10^{-5}$ from 10001 to 15000. In each iteration the energy is evaluated from a batch of $10^{3}$ random samples, and while the Monte Carlo energy fluctuates among iterations, its average value converges. Some typical training histories are shown in Fig. 12.

The final energy of the trained variational wavefunction is evaluated from 5 million samples, with Monte Carlo uncertainties shown as error bars in Figs. 13, 14, 15 and 16. In these figures we compare performance of various architectures and observe that

• •

MAF obtains lower energies for small $\nu$ and NF has lower energies at larger $\nu$.

• •

The result does not significantly depend on the initialization for small $\nu$.

• •

In the supersymmetric sector the variational energy is close to zero (compared to a typical energy scale, say the bosonic energies).

• •

Consistent improvement is observed in MAFs if we increase the number of distributions in the mixture or $D$ as in the fermionic wavefunction. However, increasing the number of layers in neural networks does not improve the results.

Appendix E Entanglement of free fields on a sphere

Solution for the projector

We wish to solve the following optimization problem: find an orthogonal projection operator $P$ such that $\|P-Q\|$ is minimal given another Hermitian operator $Q$. We will now do this in the case that $\|\cdot\|$ is the Frobenius norm. In this case, diagonalize $Q=UQ^{\prime}U^{\dagger}$ such that $Q^{\prime}$ is diagonal with diagonal elements nonincreasing. Then $\|P-Q\|$ is minimized if and only if $\|P^{\prime}-Q^{\prime}\|$ is minimized and $P=UP^{\prime}U^{\dagger}$.

Firstly we search for $P^{\prime}$ that minimizes $\|P^{\prime}-Q^{\prime}\|$ in the subspace of projectors with fixed rank $r$. It is equivalent to maximizing $\operatorname{tr}(P^{\prime}Q^{\prime})$ by definition of the Frobenius norm. Let $F(V)=\operatorname{tr}(VP^{\prime}V^{\dagger}Q^{\prime})$ for unitary $V$. If $P^{\prime}$ maximizes $\operatorname{tr}(P^{\prime}Q^{\prime})$, $dF=0$ at $V=I$ for any $dV$ in the Lie algebra of the unitary group:

[TABLE]

If $Q^{\prime}$ is diagonal with distinct eigenvalues, (E.1) implies that $P^{\prime}$ should be diagonal as well. Then the $P^{\prime}$ that maximizes $\operatorname{tr}(P^{\prime}Q^{\prime})$ should be such that $(P^{\prime})_{ii}=1$ for $1\leq i\leq r$ and 0 otherwise, and the minimal value of $\|P-Q\|$ is

[TABLE]

The projector $P$ that achieves the minimum is unique when $Q^{\prime}$ has distinct eigenvalues; if $Q^{\prime}$ is degenerate, there may also be nondiagonal $P^{\prime}$ matrices that attain the minimal $\|P-Q\|$.

The second step is to minimize (E.2) with respect to the rank $r$. If $Q^{\prime}_{ii}\neq 1/2$, the rank should be the number of eigenvalues of $Q$ that are above $1/2$. The minimum is then

[TABLE]

When one half is among the eigenvalues, there are multiple $P$’s that minimize $\|P-Q\|$.

To summarize, let $Q=UQ^{\prime}U^{\dagger}$ such that $U$ is unitary and $Q^{\prime}$ is diagonal. Then the following $P$ minimizes $\|P-Q\|_{F}$ among orthogonal projectors:

[TABLE]

And this is the unique minimum if none of the eigenvalues of $Q$ is $1/2$.

Evaluation of the second Rényi entropy

As discussed in the main text, in the case where the configuration space $Q$ has a linear structure, an orthogonal decomposition $Q=Q_{1}\oplus Q_{2}$ induces a factorization of the Hilbert space $L^{2}(Q)=L^{2}(Q_{1})\otimes L^{2}(Q_{2})$. For any pure state $|\psi\rangle\in L^{2}(Q)$, the entanglement entropy is computed as $S(\rho_{1})$, where $\rho_{1}$ is the reduced density matrix of the subsystem $L^{2}(Q_{1})$. For numerical simplicity, we now focus on the Rényi entropy (of order $\alpha\geq 0$):

[TABLE]

The von Neumann entropy is recovered as the limiting case $\alpha\to 1$. And in the following consider $\alpha=2$ for concreteness; similar methods and arguments apply to the Rényi entropies of integer orders $\alpha\geq 2$.

The decomposition $Q=Q_{1}\oplus Q_{2}$ can be implicitly specified by an orthogonal projection operator $P:Q\to Q$, such that $Q_{1}=\operatorname{im}P$ and $Q_{2}=\ker P$. For a pure state $|\psi\rangle\in L^{2}(Q)$, the reduced density matrix $\rho_{1}$ is

[TABLE]

where $x,x^{\prime}\in Q_{1}=\operatorname{im}P$ and the integral is over the subspace $Q_{2}=\ker P$. Consequently the second Rényi entropy is

[TABLE]

To further simplify the integral, let $z=x+y\in Q$ and $z^{\prime}=x^{\prime}+y^{\prime}\in Q$, so that

[TABLE]

Thus the integral in (E.7) can be done over the full space $Q$ instead:

[TABLE]

Numerically the integral in (E.9) can be estimated by Monte Carlo:

[TABLE]

where in the square bracket, the overall normalization of the wavefunction is unimportant.

The integral in (E.9) is analytically tractable for Gaussian states:

[TABLE]

where $V$ is some positive definite matrix and $Z$ is the normalization factor. Up to numerical factors, for any positive definite matrix $A$,

[TABLE]

Substituting (E.11) into (E.9) and performing the integral using (E.12), for Gaussian pure states, one obtains

[TABLE]

where

[TABLE]

The factor of $\det S$ comes from the normalization $Z$ in (E.11). It is simpler to write

[TABLE]

where

[TABLE]

In the next subsection, geometric features of entanglement for free fields are understood analytically from the formulae (E.23) and (E.24).

Derivation of the geometric features of entanglement

Consider a free field on a sphere as in (5.1) with angular momentum cutoff $j\leq{j_{\max}}$. The ground state is a Gaussian state (E.11) with $V$ diagonal in the basis of spherical harmonic modes with eigenvalues $\sqrt{j(j+1)+\mu^{2}}$ and multiplicities $2j+1$. The projector $P$ is the one that minimizes $\|P-\chi_{A}\|$, with the region $A$ being a spherical cap with polar angle $\theta_{A}$. We would like to confirm the following numerical findings with analytic computations: as ${j_{\max}}\to\infty$, (i) $S_{2}\propto{j_{\max}}\sin\theta_{A}\propto{j_{\max}}|\partial A|$ and (ii) $\operatorname{tr}P\propto j_{\text{max}}^{2}\int_{0}^{\theta_{A}}\sin\theta d\theta\propto j_{\text{max}}^{2}|A|$.

To start, observe that from (E.23) naively we would expect $S_{2}\sim({j_{\max}})^{2}$ because of the trace, and thus if $S_{2}\sim{j_{\max}}$ it must be the case that the matrix $K$ is small. Hence it is reasonable to make the approximation

[TABLE]

In terms of matrix elements of the projector, (recall that $P^{\dagger}=P$ and $P^{2}=P$)

[TABLE]

where we have noticed that the projector preserves the $J^{z}$ quantum number because of the symmetry of region $A$. Also the eigenvalues of $V$ are approximated as $j$. Subleading terms will not modify the scaling as ${j_{\max}}\to\infty$, where $j$ is typically large.

For $j,j^{\prime}\ll{j_{\max}}$, the projector $P_{jm,j^{\prime}m}$ should converge to its value at infinite ${j_{\max}}$, which is the matrix element of multiplication by $\chi_{A}$:

[TABLE]

where $\chi_{A}$ restricts the $\theta$ integral to $[0,\theta_{A}]$. Up to numerical factors,

[TABLE]

where $P_{j}^{m}(x)$ are associated Legendre polynomials.

The asymptotic form of associated Legendre polynomials $P_{j}^{-m}(x)$ in the limit $j,m\to\infty$ with $\alpha=m/(j+1/2)$ fixed ($0<\alpha<1$) is given by the WKB formulae eqs. (3.28) and (3.30) in [71]: for $\beta=\sqrt{1-\alpha^{2}}$ and $\beta<x\leq 1$,

[TABLE]

while for $0\leq x<\beta$,

[TABLE]

where

[TABLE]

Let $x=\cos\theta$. At large $j$ the oscillating region of the integral in (E.28), where (E.30) holds, is $0<\alpha<\sin\theta$. Outside of this region, the Legendre polynomial is approximately (E.29), and hence exponentially small. We need therefore only consider the region where both Legendre polynomials are oscillating. In order to get the parametric dependence of observables right, we can furthermore restrict attention to $m\ll j,j^{\prime}$. In this limit $\beta\to 1$, $\alpha\to 0$ and hence

[TABLE]

So in this limit the integrand in (E.28) can be approximated as

[TABLE]

The terms $\cdots$ necessarily oscillate strongly at large $j,j^{\prime}$ and will not contribute to leading order. In the remaining term in (E.33), in contrast, the oscillations are slower when $j\sim j^{\prime}$. Performing the integral we obtain

[TABLE]

The lower limit of integration (at $m=[\min(j,j^{\prime})+1/2]\sin\theta$) can be ignored so long as $m\ll\min(j,j^{\prime})\sin\theta_{A}$. This is stronger than the previous assumption $m\ll j,j^{\prime}$. We can now use (E.34) to evaluate observables, using the fact that $P_{j(-m),j^{\prime}(-m)}=P_{jm,j^{\prime}m}$.

The Rényi entropy (E.26) is now (with $j_{\text{m}}=\min(j,j^{\prime})$)

[TABLE]

In the second line we used $j_{\text{m}}\sin\theta_{A}$ as a cutoff on the sum over $m$, to get an estimate of the scaling with $\sin\theta_{A}$. This is the boundary law entanglement that was observed numerically in the main text.

To get the rank of the projector one must treat the sum over $m$ a little more carefully. In particular, we refrain from taking $\alpha\to 0$, $\beta\to 1$. Keeping $\alpha=m/(j+1/2)$,

[TABLE]

Here $\cdots$ again denote terms that oscillate strongly in the large $j$ limit and are therefore subleading. The integrand in the second line is directly the non-oscillating part of (E.30) squared. At large $j_{\text{max}}$ we therefore have, approximating the sums as integrals and letting $\alpha=\sin\gamma$,

[TABLE]

The integrals are most easily done by exchanging the order of integration to $\int_{0}^{\theta_{A}}d\theta\int_{0}^{\theta}d\gamma$. This result shows that the rank of the projector goes like the area of the region on the sphere, as seen numerically in the main text. The prefactor in the final result (E.41) is easily restored by noting that when $\theta_{A}=\pi$, corresponding to the whole sphere, $\operatorname{tr}P\sim j_{\text{max}}^{2}$ at large $j_{\text{max}}$.