Going beyond ER=EPR in the SYK model

We discuss generalizations of the TFD to a density matrix on the doubled Hilbert space. We suggest that a semiclassical wormhole corresponds to a certain class of such density matrices, and specify how they are constructed. Different semi-classical profiles correspond to different non-overlapping density matrices. We show that this language allows for a finer criteria for when the wormhole is semiclassical, which goes beyond entanglement. Our main tool is the SYK model. We focus on the simplest class of such density matrices, in a scaling limit where the ER bridge is captured by chords going from one space to another, encoding correlations in the microscopic Hamiltonian. The length of the wormhole simply encodes the extent these correlations are eroded when flowing from one side to the other.


Introduction
A key development in the holographic study of black holes was Maldacena's identification of the thermofield double state with the eternal black hole [1]. The thermofield double state (TFD) is an entangled state in two copies H L ⊗ H R of a quantum system, This state can be prepared by a Euclidean path integral over an interval of length β/2 in Euclidean time. In a holographic system, an uncharged black hole provides a bulk saddle-point for these Euclidean boundary conditions; in situations where this is the dominant saddle-point, on continuation to the Lorentzian picture the eternal maximally extended black hole then provides the dual description of the state (1.1). This connection is a key motivating example for the idea that entanglement is related to connection of the bulk geometry [2]. There are two-point functions involving operators in the two copies of the system O L O R which have non-zero values in this state; in the connected bulk geometry, influences can propagate from one boundary to the other, while in the boundary quantum theory, the two-point functions can be non-vanishing in the absence of any interaction due to the entanglement of the state [3]. We would like to know more generally what states in H L ⊗ H R can have such a wormhole description. Entanglement is clearly necessary (in its absence, correlations O L O R in non-interacting systems would vanish), but it has been argued that it is not sufficient. In [4,5], it was argued that in a generic state with the same degree of entanglement as the TFD, the correlations of simple operators between the two boundaries is exponentially small, so generic entangled states do not have a geometric wormhole description. Some studies have considered small departures from the TFD, finding that these lead to longer wormholes. This was first explored in [6], who considered a local operator insertion on the TFD evolved in Lorentzian time. A related investigation in [7] considered operator insertions in the Euclidean path integral defining the TFD state. We will review these studies in the next section.
Our aim in this paper is to further explore when states in two copies of a holographic system have a bulk wormhole description. We will argue that as we move away from the TFD, it is natural to consider mixed states on H L ⊗H R . 1 For small departures from TFD, this is motivated by noting that the geometry obtained in [6] or [7] doesn't depend on the details of the operators considered, so the bulk geometry can be related to a density matrix obtained by averaging over the operators. Density matrices are usually written as a sum over states in H L ⊗ H R ; we refer to this as the state frame of the density matrix. The density matrices corresponding to different bulk geometries should have disjoint support in the space of quantum states. Different states drawn from the same density matrix have the same gravitational profile, but differ in microscopic data.
We will argue that in considering two-sided correlations like O L O R , it is also useful to rewrite the same density matrix as a sum of products of pairs of operators, one acting on H L and the other on H R , We will refer to this way of writing as the operator frame of the density matrix. We can then view the density matrix as an object which takes an operator on H R and gives us a corresponding operator on H L ; that is we think of the density matrix in terms of operator conversion from one side to the other. If this operator conversion is close to diagonal, then two-sided correlations like O L O R will be of the same order as single sided correlators (like O R O R ). Thus, we argue that the density matrices of interest are those which 1) average in an appropriate way over microscopic details, 2) have disjoint support in the space of quantum states, and 3) have an approximately diagonal operator conversion for simple operators for which we expect to have gravitational duals. Going to the density matrix is essential in order to be able to formulate the last criteria, which goes beyond just entanglement. These general ideas are developed in section 2.
We will explore these questions explicitly in the context of the SYK model [10][11][12], and in a related model where the Majoranas are replaced by Pauli matrices (see for example [13]). The SYK model is a quantum mechanical system with N Majorana fermions, with an all-to-all interaction of p fermions at a time, where the strength of each such p-interaction is an independent random Gaussian. In the Pauli matrices model there are N spin 1/2 Hilbert spaces with a general random k-local interaction made out of Pauli matrices. In the large N limit, the SYK model has a nearly-conformal symmetry in the IR, and fluctuations around it are described by a Schwarzian effective action, which is the same dynamics that describes JT gravity on NAdS 2 [14,15]. This provides a useful context for studying holographic relations where we can do some calculations directly in the quantum mechanical system. One of our main aims is to study semiclassical generalizations of the TFD from the boundary perspective, in contrast with earlier works which have focused primarily on calculation in the bulk holographic picture. We will mainly do calculations in a double scaled limit of SYK [13,16], where we take p and N to infinity with p 2 /N fixed, using the chord diagram approach to calculations in this limit discussed in [17,18], but the main conclusions carry over to the usual SYK large N limit. The models and this calculational approach are reviewed in section 3.
Part of going from a state on the doubled Hilbert space to desntity matrices on it, has to do with averaging over microscopic information (keeping macroscopic data of the wormhole, such as its length, fixed). An important clue for how to do so is the following. After averaging over couplings, the SYK model has an O(N ) symmetry under rotation of the basis of fermions, ψ i = O j i ψ j . When we consider two copies of the SYK model, because we have a single average over the couplings, the averaged system is invariant under a diagonal O(N ) symmetry which rotates both the left and right fermions simultaneously. There is a similar symmetry in the random spin model. We will argue that this rotation is an example of the kind of microscopic details we should average over in constructing duals of bulk wormholes -two configurations that differ by an O(N ) rotation will look the same gravitationally. Thus we will consider O(N )-invariant density matrices ρ, or more precisely, density matrices where the only terms which break the O(N ) symmetry are explicit insertions of the Hamiltonian. In section 4 we set up the specific class of density matrices we consider, and in section 5 we use the chord diagram approach to calculate two-sided correlation functions in this class of density matrices. We show that in a suitable limit the results obtained from our microscopic calculation reproduce the holographic calculation of correlations in a shockwave perturbation of TFD in [6].
In the following section 2 we summarize the main points of the paper in more detail.
2 Summary of main points -from TFD to density matrices on the doubled Hilbert space We wish to consider generalisations of the thermofield double state from a microscopic perspective, and explore which ones could be dual to a semiclassical wormhole geometry in AdS spaces [19]. Our main aim is to study this explicitly in the context of the SYK model, where calculations in the microscopic quantum theory are (somewhat) tractable. But before setting up the details of our study in the SYK context, in this section we will set out the kind of generalisation we want to consider in a more general context and make contact with related work. We will consider states and density matrices in a bipartite quantum system with Hilbert space H L ⊗ H R . We consider theories with decoupled dynamics, so the Hamiltonian is H = H L + H R ; the connection between the two copies of the system comes only from entanglement in the quantum state. 2 Our main points are that 1) for the purposes of understanding the gravitational dynamics of the wormhole it is useful to go to a density matrix on H L ⊗ H R , and 2) there is a simple characterization of whether a density matrix corresponds to a semiclassical wormhole, which is not strictly the entanglement between the Right and Left spaces. From now on when we say density matrix we will refer to a density matrix on H L ⊗ H R .

A short review of deformations of the TFD
In [6], small modifications of the TFD state were studied, introducing a perturbation by acting on one side with an operator W (t) = e −iHtw Oe iHtw , where O is some local operator. This "timefold" operator corresponds to evolving the TFD state at t = 0 back into the past for a time t w , inserting O and then evolving forward in time by t w back to t = 0 to produce a modified state. In the limit of large t w , the dual of this construction is a geometry with a lightlike shock propagating from the boundary at early times into the black hole. The back-reaction of this shock lengthens the wormhole. This investigation was extended to consider multiple shocks in [24].
A related investigation in [7] considered what they termed the "partially entangled thermal state" (PETS) whereŌ n,m are the matrix elements of some local operatorŌ between energy eigenstates and Θ is an anti-unitary operator, for example CPT. The construction of this state is essentially a Euclidean variant of the previous construction: it's given by a Euclidean path integral with evolution over a period β R /2 in Euclidean time, followed by insertion of the operatorŌ, and further evolution over a period β L /2 in Euclidean time, as depicted in figure 1. This state was introduced as an interesting generalization of both the TFD state and the "thermal pure states" considered in [25]. In [7], these states were studied from the spacetime dual perspective using JT gravity: it was argued that the dual of such states is a bulk geometry with a particle emitted into the bulk by the insertion of the operatorŌ. They considered the regime where this particle is heavy (the dimension of the operatorŌ scales with the central charge in the holographic large central charge limit) so the emission of this particle back-reacts on the trajectory of the "boundary particle" in JT gravity, modifying the geometry. The modified geometry looks like a black hole from the perspective of each asymptotic region, but with a separation between the two horizons, producing a lengthened wormhole, see figure 1. In [7], the entanglement entropy and entanglement wedges of this bulk geometry were investigated using the replica trick, showing that the entanglement entropy is given by the smaller of the two horizon entropies. On the left, the construction of the PETS by a Euclidean path integral, and the dual bulk geometry. The t = 0 slice of the Euclidean geometry provides initial data for the Lorentzian geometry on the right.
A similar lengthening of the wormhole was also obtained in [5] when the state in the doubled theory was slightly rotated away from thermofield double, and in [26] when entangling a non-gravitating system and a gravitating (with JT gravity) one, and including the backreaction on the geometry.

A convenient purification
Usually, when going from a single to the two-sided discussion of a black hole, one purifies the thermal density matrix with another copy of the same Hilbert space. For us it will be more convenient to use a slightly different purification. Starting from the thermal density matrix on a single space ρ(β) = e −βE i |i i|, we purify it using an additional copy of the Hilbert space and a Hamiltonian with the same spectrum as the original Hamiltonian. This is true for another copy of H but it is equally true if we take H † for the second copy, and we will take our purifying space to be the latter.
That is, we just purify the density matrix by a state on H † L ⊗ H R which we write as This choice of purification arises naturally from a Euclidean path integral perspective. This is true in any dimension, but in the context of 2D gravity, this is particularly simple. Two sided black holes can be created by cutting a Euclidean space, whose boundary in the past is an open 1D surface (as in figure 1 on the left). This 1D boundary is can be chosen to have an orientation (any of the two will do). It is then natural to associate different orientations to the edges of this 1D surface, which translates into doubling H on one side by an H † on the other. One example of the usefulness of this choice is that we can modify a wormhole by the insertion of operators on this surface, and in this purification prescription the concatenation of operators action is straightforward and does not require any anti-unitaries. We can then generate a host of other states, such as states discussed in Fig. 15 of [7]: There is a one-to-one correspondence between this description using H † L ⊗ H R and the description where we purify using H L ⊗ H R (adopted for example in (2.1)) by mapping from H † to H as where Θ is an anti-unitary isometry, for example CP T . The need for this anti-unitary is because the time-evolution on H L and H † L are reversed in relation to each other. We will work in the H † L ⊗ H R purification scheme.

Going to density matrices on the doubled Hilbert space
The examples in section 2.1 share the feature that they focus on pure states on the doubled Hilbert space, which we take to be H † L ⊗ H R . The goal of this paper is to show that certain questions clarify significantly when going further to density matrices on H † L ⊗ H R . The idea is that deforming away from the thermofield double involves introducing an excitation in the bulk spacetime, but the excited wormhole geometry does not depend on many of the details of how we do it; that is, the same gravitational profile can be generated in various ways, and we can hence identify this geometry with some density matrix on H † L ⊗ H R . A concrete example would be taking the state (2.1) of [7] and averaging over the operator. The bulk geometry in [7] depends only on the mass of the particle, that is on the operator dimension, and not on the details of which particular operator we insert. It therefore seems natural to consider a duality not between the bulk spacetime and an individual pure state of the form (2.1), but between the spacetime and a density matrix ρ = i c i |Ō i Ō i |, where we average over all operators with conformal dimensions in some window, with some appropriate weights.
The reverse direction is more interesting: suppose that we have a way of associating a density matrix to each geometric profile of a wormhole. Given a density matrix ρ in some Hilbert space, here H † R ⊗ H L , it defines a probability on the states in that Hilbert space 3 . We can then consider generic states by this measure. I.e., there would be states whose probability will decrease like 1/e S , where S in the entropy of the density matrix, and there would be states where the probability will be much smaller. The former are the generic states in the density matrix. The dual of any of these generic states will look like the same gravitational wormhole.
Working with density matrices in the doubled Hilbert space will simplify the discussion, as it removes the dependence on some of the microscopic details. More interestingly it highlights a criterion for the existence of a semiclassical wormhole which is complementary to entanglement, which we discuss next.

The state frame and the operator frame
We will argue next that we can give a criteria (beyond entanglement) for when the wormhole is semi-classical in terms of this density matrix. The criteria reflects the more detailed structure of the algebra of operators in the theory.
There are two complementary ways of viewing the density matrix, corresponding to different perspectives on the emergence of a wormhole in the bulk. A density matrix We can write the density matrix as a sum over pure states in H † L ⊗ H R . We will call this the state frame As each state entagles left and right d.o.f, it emphasizes the entanglement structure between left and right. Alternatively we can write the density matrix as a sum of a product of operators acting on H † L and on H R . We will call this the operator frame, and it emphasizes the correlators O L O R between operators on the two boundaries. The relation between the sum of states and the correlator is well understood for the TFD, but in the generic case, we argue that the existence of a semiclassical wormhole is better addressed from the latter perspective. Just requiring entanglement (the first point of view) does not imply large correlation functions.
In more detail, the operator frame is obtained instead by expanding (2.5) as a sum of products of terms, one in H R ⊗ H † R and one in H † L ⊗ H L -i.e. a sum of products of operators, one acting on H R and one on H † L : where O a R are a list of operators acting on H R , andÕ b L the similar list of operators acting on H † L . Given a rich enough set of operators we can expand any ρ like this 4 with some coefficients C ab . In this form, the density matrix implies an operator conversion between the two spaces, to which we turn next.

Operator conversion/pairing
A common probe of the connectedness of the bulk spacetime geometry is to consider the two-point function of a simple operator O L O R ; if the dimension of the operator O is much larger than one (but not scaling with the central charge, so we can neglect its back-reaction), the boundary correlator can be calculated in the bulk using the geodesic approximation, and the correlator will measure the (regularised) length of a geodesic connecting the two boundaries. Thus, such two-point functions being order one 5 , and being independent (at leading order) of the details of the operator considered, is a useful diagnostic of the existence of a bulk dual.
Consider now (2.7). This object gives a pairing of an operator from the right and from the left. If we compute a two point function then Concrete examples of this are discussed in section 4. A semiclassical, weakly coupled, wormhole corresponds to a density matrix which has C ab ∝ δ ab , as will be explained later.
Another way of phrasing it is as a conversion of the operators in H R into operators on H L . A density matrix on the doubled Hilbert space literally plucks an operator from space R and inserts it into space L. To see this take (2.7) and consider it as a map of operators The interpretation is that if we start with an operator in R and move it through the wormhole, we will obtain the image O L on the left hand space. For example if we start with the particle that corresponds to O R and move it through the wormhole it will correspond to the superposition of particles encoded in the image O L . Clearly here again we see that a weakly coupled, semiclassical wormhole, is characterized by C ab being close to diagonal (for any O R that corresponds to a low energy field). The size of the diagonal C ab tells us what is the length of the wormhole.

Invariance
We argued before that the key advantage of going to a density matrix is that it clumps together many pure states, so we can avoid considering features of the pure state which are not relevant to the gravitational profile of the wormhole. In cases like the SYK model this can be made more precise. The theory has an O(N ) action which rotates the fermions, and it is broken by the random couplings in the Hamiltonian. After we average over couplings, this symmetry is restored. If we consider two copies of the theory, since we consider a single average over couplings, it is only the diagonal O(N ) symmetry which acts on both copies which is restored. Since gravity captures averaged quantities, it is sensitive to objects that are invariant under O(N ), so it is natural to consider density matrices which are invariant under O(N ) -for example, we can average (2.1) over all operators related by the O(N ) rotations. Note that this is only possible once we consider density matrices; the Hilbert space H † L ⊗ H R has just one singlet state which is invariant under the diagonal O(N ) symmetry (we will denote it by |s ). But there are non-trivial O(N ) invariant density matrices, allowing us to focus on the gravitational data.
Technically, we will not restrict strictly to O(N ) invariant density matrices; we will also consider density matrices where the O(N ) symmetry is broken only by explicit insertions of factors of the Hamiltonian. For example, the infinite temperature limit of the thermofield double state is the O(N ) invariant singlet state |s , but in |T F D T F D| = e −βH/2 |s s|e −βH/2 , the O(N ) symmetry is broken by the insertions of the Hamiltonian. We are interested in objects which generalise the thermofield double state so we allow breaking of O(N ) by the Hamiltonian in our more general density matrices as well.
This can probably be generalized to higher dimensions, with the O(N ) symmetry restrictions replaced by algebraic structures in the theory such as the OPE, but we leave this for future work.

The models, symmetry and chord diagrams
In this section we review the two models that we use -the SYK Majorana fermion [10][11][12] and the Pauli matrices model [13,17,18]. We will be interested in the case of the two sided black hole and hence the space is doubled. We will then review the chord diagram techniques for these models, in the double scaling limit. Our main conceptual points do not rely on the double scaling limit, but it will be useful computationally.

Two decoupled SYK models and their symmetry
The Sachdev-Ye-Kitaev model is a quantum mechanical system in 0+1 dimensions, constructed of N Majorana fermions with all-to-all interactions, and random (disordered) couplings. Denote the Majorana fermions by χ i , i = 1, · · · , N . These satisfy the algebra For some integer p ∈ 2N, we define the ordered multi-index I ≡ (i 1 , · · · , i p ) (with 1 ≤ i 1 < i 2 < · · · < i p ≤ N ), and the Majorana string χ I = χ i 1 · · · χ ip . The Hamiltonian of the system is given by where J is the coupling constant, but by choosing the appropriate scale we set this to one. The sum is over all possible multi-indices of length p. The J's are taken to be Gaussian random variables with zero mean and normalized variance. We denote the expectation value over these by · J , such that we have This is a slightly different normalization than, say, in [12] but going between the conventions is straightforward.
Since we would like to study the dual of the wormhole geometry, which has two boundaries, we are going to deal with two decoupled copies of this model, correlated by the fact that we have the same random coupling on both sides. We label the two boundaries as L and R (left and right), and mark it as a subscript. The Hamiltonians are the same on the two spaces, i.e., and we take the Majorana fermions on the two sides to anticommute with each other, i.e., {χ L , χ R } = 0. We will be intersted in states in the doubled theory. There is an independent O(N ) symmetry rotating each set of Majorana fermions, which is broken by the Hamiltonian. The average over the random couplings restores the diagonal O(N ) symmetry which rotates both χ i L and χ i R at the same time (it is only the diagonal symmetry which is restored as the couplings in the two theories are correlated).

Random spin model
The random spin model in 0+1 dimensions is a close relative of the SYK model, where spin 1/2 degrees of freedom assume the role of the Majorana fermions. The relevant operators that appear in the Hamiltonian are then Pauli matrices acting on the various spins, i.e., σ a i i , i = 1, · · · , N , and a i = 1, 2, 3. Let e = (i 1 , · · · , i p ) be a vector of length p of distinct integers defining a subset of the N sites (again we take 1 ≤ i 1 < · · · < i p ≤ N ), and let a = (a 1 , · · · , a p ) be a second vector of length p, with a i = 1, 2, 3. Denoting the pair (e, a) by I, we define σ I = σ (e,a) = σ a 1 i 1 · · · σ ap ip . The Hamiltonian is now where again the sum runs over all possible choices for I, and J sets an energy scale which we normalize to 1. As in section 3.1, the J I 's are taken to be independent random Gaussian variables with zero mean and unit variance, such that we have (3.3).
In the IR at finite p, the model was argued on quite robust grounds to not support a non-fermi liquid phase [27]. In fact, a related model, at fixed p, is reviewed in [28], and exhibits a plethora of complicated IR behaviors. However, we will be interested in the double scaling limit p ∝ √ N where it shows the same universal dynamics as the double scaled SYK model which has nearly conformal physics in the IR.
We take two decoupled random models (3.5) but with the same random coupling, as in the Majorana SYK case. The Hamiltonian we are interested in is therefore 6 (3.6) For later reference, we can clump together degrees of freedom associated with spin i from the left and from the right Hilbert space into a single four dimensional Hilbert space, which we denote as H i .
The symmetry of the system after averaging over the random coupling is SU (2) N S N , where each SU (2) acts on H i (which is made out of an singlet and a triplet under this action), and S N is the permutation of the N species of Pauli matrices on the left and the right boundaries at the same time.

Probe operators
Once we map the set of states that correspond to some wormhole, we would like to start probing it. A set of convenient probes are two sided correlation functions, but we still need to choose of which operators.
The operators that we will choose as probes are random operators. A random operator M of lengthp is defined to be where the couplingsJ I are again random Gaussians with zero mean and normalized variance, where · J denotes the ensemble average. Notice that the coefficientsJ I are independent of the coefficients J I . We could consider instead the operators to be "monomials" σ I or Ψ J for a specific I or J -nothing would have changed in that case, since we can view the random operator above as computing the properties of a generic monomial.
There is, however, a deeper reason for choosing such operators. Consider the case that we have some higher dimensional theory that flows to a (near) extremal background of the form AdS 2 × M , and suppose that this IR can be thought of as an SYK model in some variables (which may not be related to the original UV d.o.f). The set of observables correspond to fields defined in the UV of the entire background. The local energy-momentum is one of them, but it is typically part of a whole tower of similar ("single trace") operators. If the Hamiltonian in the full theory flows to some random Hamiltonian on the near extremal BH degrees of freedom, then we can expect the same to be true of all other UV operators. So the natural probes are random operators in a statistical class similar to the Hamiltonian, except that that we allow them to have different quantum numbers, which in this case is just the conformal dimension/length of the operatorp (more precisely, the conformal dimension of the operator in the IR is p/p, where p is the length of the Hamiltonian).

The double scaled limit and chord diagrams
In this work, we will be interested in computing quantities of the form Z = Tr(e −βH ) J , where the subscript J means averaged over the ensemble of random couplings. We will also be interested in the 2-sided correlation functions in the density matrices discussed above.
This computation will prove to be especially simple in the limit in which we take N → ∞, and scale p ∼ √ N . We refer to this limit as the double scaled limit. Namely, we keep fixed the parameter We can define analog quantities to λ, q for probe operators other than the Hamiltonian byλ pp N Random spin model, (3.10) wherep is the length of the new operator, and for both models we defineq,˜ bỹ q = e −λ = q˜ ,˜ =p p . (3.11) The main advantage of working in this limit is that it enables us to use chord diagrams in order to compute various quantities. These computations have already been performed in [17,18], and the following summary of the methods used will serve as the basis for the computation in the rest of the work.
As an example, we'd like to compute the partition function 7 . First we Taylor Tr(H k ) J , and then we evaluate each of the moments Since the coefficients J I are random Gaussians, upon taking the ensemble average we Wick contract them. Together with the cyclic structure of the trace, this allows us to represent the moment m k as a sum of chord diagrams, as in figure 2. The chord diagram then determines the arrangement of repeating Majorana fermions or Pauli matrices inside the trace. A little combinatorics is then needed to show that in the double scaled limit, the value of every such chord diagram is given by the number of chord intersections, so where CD(k) represents chord diagrams with k nodes, and #intersections is the number of pairwise intersection of chords in the diagram. This sum can be evaluated using a transfer matrix method in which one goes linearly along the circle, adding a new chord or closing a chord (with appropriate weight) at each node. A full explanation of the chord techniques is given in appendix D. The advantage of this limit is that, in our application, it "combinatorizes" the geometric relation between the two universes. In the one side case, generally, the chords are the objects that carry correlations of the different terms in the Hamiltonian from one time to another, i.e, from one boundary point to another in the dual GR picture. When we have two spaces we can have chords stretching between pieces of the Hamiltonian on the two different sides. This is a combinatorial manifestation of a bridge in spacetime that forms between the two systems.

Microscopics of wormhole density matrices
In this section we discuss the density matrices that are relevant to the Pauli matrices model (sections 4.1 to 4.3) and to the Majorana case (section 4.4). We discuss their form in the state frame and in the operator frame, and their entanglement, which is easily read form one frame, vs. their length, which is easily read from the other frame.

Density matrices in the random spin model
Let us set up the states and density matrices we want to consider in the random spin model. We will begin by setting up our conventions with a single spin, and then go on to consider N spins and the limit of large N .
4.1.1 Hilbert space structure of the two qubit system Basis of pure states: As a toy model, let us study a double copy of a system with a single spin. We think of the two copies as living on the left and the right boundaries. Remembering that our random spin model had the symmetry SU (2) N S N after ensemble averaging, we will construct SU (2) invariant density matrices.
We will start with pure states. The doubled Hilbert space of a single spin has four states, which split into a singlet and a triplet under the diagonal SU (2). Pure states on H † L ⊗ H R can be represented as operators from H L to H R (or operators on H since H L,R are isomorphic). In the present case, the identification of states with operators is where m takes values in −1, 0, or 1, corresponding to the three Pauli matrices τ a . An element g of the diagonal SU (2) acts on the operators as O → gOg † , so 1 is a singlet and τ a are a triplet under the action of SU (2).

Invariant density matrix
We are interested in the density matrix which is invariant under the diagonal SU (2). The most general invariant density matrix is and positivity of the density matrix implies 0 ≤Â ≤ 1. This object lives in be Pauli matrices and the identity, defined in the following spaces (4.3) Using this notation, we can write the correspondence 1 2 so the density matrix becomes This is the form of the density matrix in the state frame.
The operator frame We can translate this to the crossed channel using the Fierz identity This gives us the density matrix in the crossed channel notation This is the form of the density matrix in the operator frame.
In terms of the variable A, the density matrix has the singlet with probability 1+3A 4 , while each triplet has a probability of 1−A 4 .

Invariant density matrix of the large N random spin model (single shot)
Let us now move to our main interest, which is the N spin system. We take N copies of the spin-1/2 system, i.e., H = (C 2 ) N . The symmetry of the two-sided system here is, as we have already explained, G = SU (2) N S N , where each SU (2) is the diagonal symmetry on the two copies of H i , and S N is the permutation of fermion species. The simplest G-invariant density matrix is attained by taking the product of N copies of the single spin density matrix above, with the same value A for each spin: Our investigations will focus on this case, which we will refer to as single shot density matrices. In appendix B, we consider a more general class of G-invariant density matrices, where we take different values of A for different spins, and argue that (4.8) is preferred as the largest entropy density matrix for fixed wormhole length. Thus, this represents the broadest coarse-graining over the microstates, and is a natural candidate for a dual of a smooth semi-classical geometry, which does not distinguish between microstates. In section 4.3, we discuss the generalization to multi shot density matrices, which we leave for future work.
To make contact with the thermofield double, we need to include thermal factors on the left and on the right to suppress contributions from high energy states, considering the density matrix where Z normalizes it to have trace 1. We will use G when discussing these specific density matrices, and we will reserve ρ for discussing density matrices in general.
|s i is the product of invariant singlet states, which is the purification of the infinite temperature thermal density matrix (i.e., the identity matrix on H R ). Thus adding these thermal factors converts it into the thermofield double state with temperature β = (β L + β R )/2, as also discussed in section 2.6. For A = 1 the density matrix depends only on β L + β R , but for a general A it depends on both separately.
Thus, (4.9) is an interesting direction to generalize the thermofield double. The SU (2) N S N symmetry is broken only by the explicit insertions of the Hamiltonian. As discussed in section 2.3, this density matrix can be thought of as taking the PETS state (2.1) and averaging over a class of operators. We will argue below that (4.9) are good candidates for the dual of semi-classical wormholes on the gravity side, where we keep the gravitational profile of the wormhole and forget about the microscopics of the particles that created it.

Properties of the single shot density matrices
We now discuss several properties of the single shot density matrix: we comment on a matching property for A = 1, calculate the entropy of these density matrices, briefly discuss the length of the corresponding bulk wormhole -arguing that for A < 1 the wormhole gets longer, as in previous studies of generalisations of the thermofield double -and argue that the overlap of these density matrices for different values of A vanishes in the large N limit.
TFD and matching property: This means that if we act with H R on G then we can convert it into an action of H L . This can be viewed as the technical reason why the thermal state has zero length in the interior of the wormhole. In fact, this is more general. The length of the wormhole corresponding to G(1) is zero for any operator, which means that something like (4.10) should hold for any operator and not just H. We will refer to this the matching property. It is easy to verify that for any Pauli matrix, This means that any operator on H R made out of a string of Pauli matrices can be matched with (or fully converted into) an operator on the left Hilbert space. This is done with relative strength 1, which translates into a zero length wormhole. When we discuss the thermofield double with finite temperature, we conjugate G(1) with e −β L,R H L,R . The core of G(1) still matches operators in left and right but now there is an additional contribution of the thermal suppression factors which we interpret as a contribution in Left and Right separately outside the horizon.
This matching property is distinct from the general operator conversion we discussed in section 2.5, and is special for G(1) and the thermofield double. When we consider G(A = 1), we will not have such a matching between operators on the left and the right. We can always use the operator conversion property introduced in section 2.5 to push an operator on the right through the density matrix to an operator on the left, but we can't in general find a (simple) operator acting on the density matrix on the left whose action matches that of an operator on the right as in (4.10).
Indeed, already for a general pure state in H † L ⊗ H R such a matching may not be possible; matching an operator on the left to an operator on the right requires sufficient entanglement, and will fail for example for product states. But when we consider density matrices, the matching problem is more severe: a relation like (4.10) is a four-index tensor equation. If we try to solve it for the left operator O L given a right operator, we have far fewer variables than equations, and there is no general reason to expect a solution to exist. As we consider A = 1, we have a superposition of pure states in the density matrix; if the solution of matching for different states in the superposition is different, no solution will exist for the density matrix itself.
Entropy From (4.5), the single species density matrix G(A) contains the singlet with probability 1 −Â and each basis of the doubletÂ/3,Â = 3 4 (1 − A). This means that the entropy of the large N density matrix is If we take A to be fixed in the large N limit, the entropy scales like N . This is too large because summing over e cN states is expected to correspond to additional horizons embedded in the background. In the explicit calculations, we will study the models in a double scaling limit of large N, p, with fixed p 2 /N . We will find (see next section) that the density matrices of interest in the double scaling limit have which scales like √ N log(N ). The same argument holds qualitatively for the finite p model (see section 5.2).
The origin of this entropy is the following: Recall that in order to obtain a wormhole we are inserting some source which backreacts on the geometry. The entropy is just the entropy associated with this source. We will see below that a also determines the length of the wormhole, so in this case there is a straightforward relation between the entropy of the states that can generate a specific geometry and its length. This might be related to how one identifies complexity in the bulk [29,30], but for this one would have to make a more precise relation between complexity and the entropy in (4.13), as well as study it for more general density matrices.
When we turn on temperature, G(A) gets dressed to become G(A, β L , β R ). Since the dressing is by conjugating G(A) the basic structure of the gluing remains the same. Actually, a formal number of singlets and triplets can be defined in the dressed operator, by computing the expectation value of an appropriate projection operator but we will not pursues this further in the current paper. We expect however that the entropy will pick up some β L , β R dependence but that the N scaling will remain the same.
Length of the wormhole We will argue that the length of the throat is proportional to − log(A). A more precise discussion of the length of the wormhole will be given in section 5, but for now let us give the gist of the argument.
On the gravity side, if we have an operator which corresponds to a particle of mass m, then a two sided correlator (at time zero on both sides) will receive a contribution of the form where L is the length of the wormhole, and M L (M R ) is the insertion of the operator on the left (right). For simplicity we will take m 1 such that mR AdS ∼ ∆, the dimension of the operator.
We realize an operator with conformal dimension ∆ as a random operator of length p = ∆ · p, as discussed in 3.3. We would like to identify a contribution of the form (4.14) to the two sided correlator when evaluated with G(A), A = 1. I.e., in the QM side we compute and focus on terms which look like the RHS of (4.14). We explain the computation a bit more in section 5, but the structure is clear. To get a non-zero answer, we need to have the same monomial of Pauli matrices both on M L and M R -say that it is characterized by an index setĨ (of lengthp). For each index i 1 , ...ip in the setĨ we need to use the Aσ L σ R term in that index to get a non-zero answer. This means that we get a contribution of the form (4.16) We equate this with e −mL to obtain that L/R AdS ∝ −p · log(A). This is the case for any p. Thus, we see that taking A < 1 decreases the value of the two-sided correlator, corresponding to a longer wormhole in the gravitational dual. In the double scaling limit we consider in the explicit calculations in the next section, A = 1 − a/p for fixed a, and hence L/R AdS ∼ a.
Overlaps of G(A)'s: Providing the gravitational data for any semi-classical gravitational background will determine a density matrix. But since different semiclassical wormholes correspond to different quantum states, then it should be that the overlap in the support of two such measures is small. One can verify that this is the case. Consider G(A 1 ) and G(A 2 ) for A 1 = A 2 . Since A determines the length then these two density matrices corresponds to wormholes with different length -i.e, they correspond to different semiclassical geometries that are macroscopically different. Recall that for each species G(A i ) = α i P s + β i P t , where P s,t are projection operators on the singlet/triplet (α i + 3β i = 1), then a reasonable measure of the overlap of the density matrices 8 is given by in the large N limit for (α 1 , β 1 ) = (α 2 , β 2 ). Actually, in the double scaling limit we have to be a little more careful since α i ∼ 1 − a i /p, but the overlap still goes to zero.

Transition matrices and multi shot spaces
Transition Matrices: Our starting point was to consider an object in ρ ∈ (H † L ⊗H R )⊗ (H L ⊗ H † R ) either as a pairing of operators, one on H R and one on H † L (the operator frame) or as a density matrix on H † L ⊗H R (the state frame). The semiclassical nature of a wormhole manifests itself simply in the first approach, and the entanglement structure in the second approach.
In the state frame ρ = α c α |ψ α ψ α | where |ψ α are some states in H † L ⊗ H R . The sum over α averages all the ways that we build a wormhole with a fixed given gravitational profile. It is clear, however, that we can just as well build a mixed state transition matrix (or Pseudo-entropy [31]) where |φ α are a distinct set of states than |ψ α . The only requirement that we impose is that (4.18) has no explicit violations of the symmetry other than by insertions of the Hamiltonian.
The simplest example is the following. Consider (4.9) or (4.21). For β 1 = β 2 it describes, in Minkowski space, a symmetric L ↔ R configuration with a massive particle sitting in the middle. This is a density matrix of the form α c α |ψ α ψ α |. Suppose we want to change its initial entry point and its final point -this would be of the form (4.18). This would then be given by the expression where t f , t p parameterize the entry/deprature point of the massive particle.
Multishot wormholes: The single shot density matrix (4.9) and its generalization (4.19) can be graphically depicted as in the first line of figure 3. The labels on the external lines encode which Hilbert space is associated with them. The ⊗ on the external legs indicate an action of some e −β i H , where β i can be complex.
If we want to have two, more or less parallel, heavy particles generating a wormhole, we can use the expression in the second line of that figure (the arrangement of the external legs is the same). We can think about it as if we introduced another Hilbert space of states that live in the interior of the wormhole -let's call it I. We can then think about the density matrix as a single shot linking R to I and then another shot linking I to L. The two shots are invariant (under SU (2) N S N or under SO(N )) but they are linked by insertions of a Hamiltonian evolution, which is the only source of symmetry breaking.
The latter is the simplest configuration. Clearly the most general transition matrix -and we conjecture, the most general semicalssical gravitational background that corresponds to a transition matrix -is given by a complicated network of shocks connected between themselves in different ways -for example the third line of that figure. In this case it describes 5 shots, some of which cross others. A general gas of particles in a wormhole will be described by a dense net of shots connected to each other in different ways, which encode the initial and final states of the wormhole. We are going to leave the study of these networks for future work.

The Majorana SYK model
So far we have been working on the random spin model introduced in section 3.2. Let us now consider the Majorana model of section 3.1 and discuss it in a parallel way. In this section, we define a class of density matrices for the SYK model, which turns out to have the similar property as the ones we've discussed before.
As in the random spin model, we are interested in density matrices invariant under the symmetry of the two-sided system. For the Majorana model, this symmetry is a diagonal O(N ) which rotates the left and right spinors simultaneously. Finding invariant density matrices is again equivalent to splitting H † L ⊗ H R into irreducible representations of this group. This is a cleaner representation theory problem than in the spin model case, and we have worked out the general invariant density matrices in appendix C. We find there that a convenient basis for the space of invariant density matrices is provided by the (4.20) In the large N limit, density matrices with different values of A have disjoint support, and we take the density matrices G M (A) for different A as our candidates for the dual of gravitational wormholes, analogous to the G(A) in the random spin model defined in (4.8).
We can add a temperature by defining the density matrices we represent the Majorana fermions in terms of Pauli matrices. The convenient representation for us is Here, σ j,L lives on the Hilbert space with the (2j − 1)-th and (2j)-th left Majorana fermions, and likewise for the right Majorana fermions. By using this representation, we have This is the density matrix in the operator frame.
In the state frame (or in order to understand how ρ(A) can be prepared using the Euclidean path integral) we repeat the Fierz identity exercise as in the Pauli matrices and write the density matrix as where τ a are defined the same way as in Sec. 4.1.1. The entropy of this density matrix can also be computed, which is Special values of A and TFD states It is immediate from the eigenvalue decomposition of ρ(A), the density matrix is pure at A = ±1. This is in contrast to the random spin model, where A = 1 corresponds to a pure state, while the other end of the spectrum, A = −1/3, was not a pure state, due to three states in the triplet. Let us start from A = 1. This is, in the Euclidean preparation, having the unit operator inserted in the thermal half cycle. This therefore corresponds to a TFD state at infinite temperature. We can also think about it in the following way, that we have the relation and since G M (A = 1) = |Ψ + Ψ + | is pure, this gives us the defining relation for the TFD state [20], The case A = −1 is similar, but with a little twist. This corresponds to inserting the fermion parity operator between the two half spaces when preparing the state (still as in figure 1), and this has the effect of changing the definition of parity between the left and the right boundary. This operator which is yet another choice of the TFD state given in [20]. The difference between A = ±1 is simply the difference in definitions of the TFD state, if we twist the choice of the CP T operator using the fermion number, in the definition of the TFD, |TFD = n |n ⊗ CP T |n .
PETS interpretation Equation (4.28) also gives us access to what operators we need to insert if we want to construct this space a-la PETS states [7]. This means that the density matrix can be understood as inserting the operator 1 with probability (1+A) 2 4 , τ 1 or τ 2 with probability 1−A 2 4 , and the fermion parity operator with probability . For the density matrix to be unitary, we need −1 ≤ A ≤ 1. This is favourable, since the number of these operators inserted in the path integral follows the multinomial distribution, similar to the Pauli matrix case. For example, the probability of finding operator ψ for K times follows the probability of and the distribution has a peak at large N .

Explicit correlators in two-sided systems
In this section we will evaluate the partition function and two sided correlation functions in our proposed (unnormalised) density matrices in the random spin model, i.e., with G(A; β L , β R ) defined in (4.9). We would like, for example, to see how the length of the wormhole depends on A. We will set up the machinery mainly in the double scaling limit (and will also comment on the scaling of A in the regular large N limit). We will focus on the random spin model -the result of the Majorana SYK model is exactly parallel.

Integrating over the random couplings
We will begin with the partition function By expanding everything as a power series in H L,R , we will generally have to compute the moments In appendix A we detail how G, H L and H R are contracted in order to implement the trace. The upshot is that we can rearrange the insertions of H L to bring it to the ordinary form of a matrix product where a lower index is contracted with an upper index of a matrix to the right. We will denote this by Tr with no subscript. If we have some insertions of an additional operator then we need to keep track of their location. For example, Next, we proceed as usual by using the fact that J I JĨ J = δ I,Ĩ . We then have where k = k L + k R . In the last equality we changed to an ordinary trace over H, by flipping the order of operators in space L. We have written the sum as a sum over all the possible pairings, and over the composition of the index sets. We denote the set of all chord diagrams with k nodes as CD(k), such that the sum over k pairs of I's can be rewritten as pairings of I = CD(k) . We can now apply the chord technology as in appendix D: the σ matrices on the right "bubble" chords in H R . When we reach G, we are typically left with some chords open, i.e., an odd number of Pauli matrices in some of the species. We now need to take G into account and carry out the trace on H R . Whenever there is an odd number of Pauli matrices in a species, we need to use the A σ R σ L (for that species index) which will convert it into an incoming single σ L . So the role of G(A) is to convert 9 , in the species needed, an index from H R to the same index in H † L . After G(A) converts chords from space R to space L we then continue the bubbling process in space L. Of course various weights, associated with A, are incurred during this conversion -we will discuss them in a moment. This is illustrated in figure 4.  (G(A)σ R,I 1 σ R,I 2 σ R,I 3 σ R,I 2 σ L,I 3 σ L,I 1 σ L,I 4 σ L,I 4 ) = Tr (σ L,I 4 σ L,I 4 σ L,I 1 σ L,I 3 G(A)σ R,I 1 σ R,I 2 σ R,I 3 σ R,I 2 ), contributing to m 4,4 . Note that the ordering of operators in space L is with respect to the their ordering in the ordinary trace. This is the 3rd line of (5.5).
In this diagrammatic language, we can now present the different weights associated with chord intersections.
Weights of chord diagrams in the double scaling limit As discussed in appendix D we can neglect the overlap of three index sets and treat each overlap of two index sets as independent. We can now start computing the contribution from two chords, and then multiply each contribution to evaluate a specific chord diagram.
Let I 1 , I 2 , ... denote the index sets associated to the different chords. We decompose the Hilbert space into site-wise components which we denote by the index i. Traces are always ordinary traces. • If i appears both in two chords -I 1 and I 2 Similarly, we list all the possible contribution from Hilbert space H i when i appears both in I 1 and I 2 : In order to get the total contribution of the two chords, we can multiply contributions from each of the species's Hilbert spaces. An index set contains p indices, and so a chord running inside space 1 or 2 gets a factor of 1, and a chord running from space 1 to 2 gets a factor of A p , where p is the number of species that don't appear in any other chord, as one can see from (5.7) and (5.8).
In the double scaling limit p differs from p by some additive finite amount (whereas they both scale like √ N ). We now also set A = 1 − a/p. Hence the final result in the double-scaling limit is This can conveniently be summarized in figures as (5.13) We depicted chords weighing 1 in blue, while the ones weighing A p in green. The red line is there to remind us that we are computing the chord diagram in the presence of G(A). Now that we assigned a value to each chord of two types above, we discuss the contribution coming from the intersections between two chords: as in appendix D, we can change variables in (5.5) from I 1 , · · · , I k/2 to the index overlaps m ij , with i, j = 1, · · · , k/2, and change the measure accordingly. As discussed there, the measure is the Poisson distribution with parameter p 2 /N , and the factor N p −k/2 ensures this distribution is properly normalized. This means that in order to find the contribution of a pair of chords, we take the relevant factor from (5.9) for each overlapping index, and the total number of indices follows the Poisson distribution. The results are summarized in figure 5. Figure 5. Factors of chord configurations in the two space construction. We compute the relevant factors by multiplying the factors from (5.9) by the number of times they appear, which follows a Poisson distribution, and setting A = 1 − a p , as we're working in the double scaled limit.
For example, the third diagram is given by where we've used the fact that in the double scaled limit p 2 /N is finite, as well as A = 1 − a p .
In conclusion the moment m k L ,k R is given by where CD(k L , k R ) are chord diagrams with 2 regions involving k L , k R nodes respectively, and by crossings we mean chords that cross from left to right.
Transfer matrix method Next we would like to use the transfer matrix method in order to compute the moment (5.14). To do so we use the same auxiliary Hilbert space and transfer matrix defined in appendix D. The similarity between (D.9) and (5.14) means that we need to do only small modifications to the moment computed by the transfer matrix (D.12). Indeed, the new ingredient is B #crossings , which can be accounted for simply by an operator insertion between the two regions. This means that we can compute the moment m k 1 ,k 2 by The insertion BN is to account for the number of green chords, and T is the transfer matrix defined in (D.11).
The moment (5.15) was previously computed in [17], where it arises in the computation of the two-point function of a random operator in a single copy of the system: if M B is a random operator of length p , the two-point function involves which was shown in [17] to be given by (5.15) with p = − p log B λ . This is in line with our expectation that the partition function of the two sided space is given by an insertion and extraction of a massive operator when taking the Euclidean boundary to be a circle (as in figure 1).
This correspondence can be generalized to include higher point functions. For example one just takes the 2nd line in equation (5.4) and converts it into a chord prescription and we take O L = O R , an Hermitian random operator of lengthp, which is replaced by a chord whose index set is of this length. This can also be written as a four-point function in a single SYK copy We will borrow the formulas for this expression from [18]. In a similar manner, the n-point function evaluated in a doubled system over the density matrix G(A) can be mapped to n + 2 point functions evaluated in a single copy.

Size of correlators and length of wormholes
We would like to get some sense of the relation between the length of the wormhole and the parameter A appearing in G(A, β L , β R ) in (4.9). The computation was basically done in section 4.2, but we have just now derived all the factors in a more systematic way.
Consider a two-sided two-point function of an operator O made from a sum over strings ofp Pauli matrices with random coefficients as defined in (3.7). Given that the operator is made out ofp symbols, its IR conformal dimension isp/p. We will place this operator in the Left and Right Hilbert spaces at times t R = t L = 0. First take β L = β R = 0. Then as in the discussion in the previous section, each Pauli matrix that gets converted from Left to Right by G(A) gives a weighting by A, so the conversion of an operator withp Pauli matrices gives Adding temperature to the system (by going to ρ(A, β L , β R )) will slightly modify the calculation, but we still need to have that the chord of the probe operator goes from Right to Left, hence we will still have a contribution like (5.20). There will be a modified prefactor that depends on the β's, but thep scaling of the RHS of (5.20) and the qualitative dependence on B will remain similar.
On the gravity side, the two-point function behaves as e −mL , where L is the length of the wormhole, and for heavy operators the mass m is related to the conformal dimension by m = ∆/R AdS =p/pR AdS . Hence we have When A = 1, then the wormhole is of length 0 as we expect, and the wormhole becomes longer as we move away from the thermofield double. This is consistent with the observation above for the entropy of ρ(A) that as A moves away from 1 the entropy of the density matrix increases -as the length of the wormhole increases, then more states can be accommodated in the wormhole (with the same gravitational profile).
The formulas given above are for the double scaling limit. In the finite p case, one might worry that the situation is different because B = A p is finite for any value of A and not just A close to one (p is a fixed integer and need not be particularly large). We would like to argue however that a semiclassical wormhole still exists only for A close to 1 (as some power of 1/N ). The reason is that for finite p, one obtains semiclassical gravity by going to low temperatures, which means a high power of H if we use the moment method. In fact one has to take the number of H's to scale like a positive power of N . To have a semiclassical wormhole we similarly require a large number of chords going between the right and left spaces. I.e, the number of chords that go from left to right scale like n LR ∼ N α , for some α > 0. In order to have a finite B n LR we need again to take A = 1 − a/N α .

Partition function and two point functions
By using the transfer matrix and the expression for the moments given above, we can compute the partition function and the 2-point functions of the two-sided random spin system. Because of the discussion above, we can just cite the result given in [17] to get the final result.
Partition function For the partition function, the result is (a k ; q) n , (e ±iθ ; q) ≡ (e +iθ ; q)(e −iθ ; q) Double-sided two-point function For the double-sided two-point functions we can cite the results for the crossed four-point function Remember thatq = e −2 pp N , as defined in (3.11). The basic hypergeometric series 8 W 7 is defined in (E.5). We can also check that it reduces to the partition function when we take M → 0.

Low energy limit
The formulas above are quite general, and can be generalized further to any correlator across the wormhole, and in fact, using a discussion similar to 4.3, can be generalized to particles entering from the past singularity or leaving through the future singularity. Here we will have a more modest goal of just testing our machinery using the two sided correlator in the shock wave [6]. The shock wave is a particularly convenient background as the two-sided correlation function was evaluated for any length/shockwave strength, and it plays a key role in the study of quantum chaos in black holes.
We will first review some formulas that have to do with the low energy limit and then we will specialize them to the shock wave in the next subsection.

Useful formulae
The low-energy limit we take is the same as in [18], where one takes q → 1 and the lowenergy limit at the same time in a consistent way, and is known to reproduce the result of the Schwarzian quantum mechanics, or equivalently, of JT gravity. We summarize below all the necessary formulae without derivation, for which the reader is referred to [18] and references therein.
First take q = e −λ and take λ → 0. At the same time, take This E(θ) is the energy eigenvalue of the SYK model in the double-scaling limit, as discussed in the previous sections. In the following, for the sake of simplicity, we shift E(θ = π) to 0. Also, in the λ → 0 limit, we have where B = q B . This B can be understood as the operator dimension of the fictitious operator M B introduced in the previous subsection. In particular, in the double-scaling limit A = 1 − a p , we have Here a = O(p 0 ), but can still depend on λ non-trivially. We have also defined M by usingq ≡ q M , and it is the dimension of the random operator, whose two-point function we compute later. We will eventually use the saddle-point approximation to evaluate the forthcoming integrals -the consistency of the saddle-point approximation is that we are in the low temperature limit, which is In principle, this should be checked every time one does the saddle-point approximation, but we refer the reader to [17] and proceed.

The low-energy limit and the shock wave
Double-sided two point function We first compute the low energy limit of the unnormalised double-sided two point function. We begin with expression (5.25), for the two sided correlator with R and L insertions att L,R + t d . We will eventually take t d to infinity, but we shift the dependence on t d into a time evolution of M B . This shifts the entry and exit point of the defect in the past and future singularities. In other words, we have the unnormalised density matrix The expression that we get is therefore (with a small cyclic rearrangement of the terms) We now take the low energy limit using the formulas presented above. where This is equivalent to [32]. As in the single-sided case, we can a posteriori justify that ω 1,2,B are much smaller than M (for further justification see [17]). The integral over u can be done using the contour integral -at large M , it can been shown that the only relevant poles are at u = 0 and at u = i(y 2 + y 4 − y 1 − y 3 ) [32]. More precisely, the other poles are suppressed for any t d , and between these two poles, one of them is dominant depending on whether t d goes to ∞ or −∞. This corresponds to the two possible orientations of the shockwave.
We first analyze the u = 0 pole.
where as always ν B,1,2 ≡ ω B,1,2 2 √ M . The ν B integral can be performed by using B , ν B ν 1 , ν 2 at the saddle point and then using the contour integration, We get The u = i(y 2 + y 4 − y 1 − y 3 ) pole can likewise be evaluated. It is much smaller then u = 0 when we take t d → ∞ (it is the dominant pole when t d → −∞).
We can now do the M integral, which is localised at the saddle-point 39) and the final expression becomes Finally, we need to normalize by the partition function. Since we get the partition function by taking M = 0, and hence taking the saddle where E 1 = E 2 = 0 in the integral in the above expression, we simply get It therefore follows that the normalised two-point function in the low-energy limit becomes which should be used in the next subsection, where we compare with the gravity computation. A comment is in order. Our result indicates that the normalised two-point function does not depend on the difference between β 1 and β 2 . This is natural, because we demanded that M is much smaller than B . In the dual gravity picture, changing β 1 − β 2 corresponds to changing the location of the end-of-the-worldline created by a particle of mass M [32]. In the approximation we are using, this is too light to affect the two-point function that we are computing. The effect will be visible at subleading order in the saddle-point approximation, since the cancellation of 2 cos π between the partition function and the unnormalised two-point function is due to the approximation M , ν B ν M in the ν B integral.

Comparison to JT gravity
In the previous section we started from a particle going from the past region of the black hole to the future region. This is slightly different than the construction in [6] where one throws in the particle from the spatial (field theory) boundary. Still, in the limit t d → ∞ the particle grazes the horizon and we expect it to converge to the shock wave geometry. The two-point function that we would like to compare is determined by the renormalised geodesic distance d of two boundary points in the shockwave geometry, where m probe is the mass of a probe particle. Since JT gravity is an s-wave reduction of the three-dimensional Einstein gravity [33], we can use the result from [6], which is where the operators are inserted at timest L,R on the two sides, and T BH = R BH /2πl 2 AdS . α shock quantifies the strength of the shock and is given by α shock = (E/M BH )e 2πT BHtw wheret w is the time of shock and we take it to infinity as we take E/M BH → 0 keeping α shock fixed 10 . By Fourier transforming in terms oft L,R , this can be rewritten as This expression is of the same form as (5.44), which supports the claim that the density matrix (5.32) is dual to the shockwave geomtry in JT gravity. The parameters between the two are matched byω 2πT e 2πT t d (5.48) 10 The computation was done in the theory of three-dimensional gravity, but since they considered the shockwave to be S-wave, we can do the dimensional reduction and then the result is the same as in JT gravity.

Summary and Outlook
We have studied generalizations of the thermofield double state, looking for suitable density matrices on a doubled Hilbert space that could correspond to smooth semiclassical wormholes in a dual gravitational description. We argued that in generalising the thermofield double state, it is natural to consider density matrices on the doubled Hilbert space, averaging over microscopic details that do not affect the gravitational profile in the bulk dual. Studying density matrices also makes it possible for us to translate from a state frame description to the operator frame description, which allowed us to formulate a criterion for the existence of a wormhole in the dual description in terms of the operator conversion properties in the operator frame description.
We have considered these questions explicitly in the context of the SYK model and a related random spin model. These models have a diagonal symmetry after averaging over the random couplings -O(N ) for the SYK model and SU (2) N S N for the random spin model -and we argued that we should consider density matrices invariant under these symmetries up to explicit insertions of the Hamiltonian. We focused our considerations on a one-parameter family of invariant density matrices in each model, (4.9) for the random spin model and (4.21) for the SYK model.
The minimal length of the wormhole is simply encoded in the parameter of the density matrix (and can also be simply related to its entropy). The ER bridge is encoded by correlations of the random terms of the Hamiltonians on the two sides, and the length is encoded by the extent to which correlations are suppressed.
For the random spin model, we calculated the two-sided correlation function for probe operators in this density matrix directly in the spin system, and found a nice correspondence with previous gravitational calculations of shock wave deformations of the eternal black hole in a low energy limit. The result is actually the same for the Majorana once we reduce the latter to the same chord prescription, as in [18]. This provides a proof of principle calculation demonstrating that our proposed density matrices have the expected properties to correspond to semiclassical wormholes.
There are number of potential directions for further development: • It would be interesting to further study the properties of our density matrices and understand their dual gravitational description in more detail.
• The density matrices we considered are "single shot": they correspond to a single localised perturbation in the wormhole in the bulk. It would be interesting to extend our study to the "multi shot" transition matrices introduced in section 4.3 and relate to previous discussion of wormholes with multiple shocks in the bulk.
• From the point of view of the invariance property, there are many possible generalizations of the simple density matrices we considered. The most general invariant density matrices probably do not have a simple gravitational description, but there may be other classes which do; it would be valuable to explore this space further.
• These density matrices correspond to wormholes with two boundaries. It would be interesting to also understand the general description of multi-boundary wormholes [34][35][36], in particular if we could shed some further light on the nature of the entanglement structure in the dual of these geometries. In the SYK or spin models, it is straightforward to set up density matrices on a Hilbert space with k factors, but it is less clear in this case what are the natural density matrices to consider. Our operator conversion and invariance discussions do not have a straightforward extension.
An operator from H † to itself is given by However, both O and U are really defined in the same class of objects and can be identified. In short, an operator on H can be used to define an operator on H † when we act with the other index on the vector of coefficients 11 .
2. Returning to our density matrices, G now carries the indices 3. When we compute a two sided correlator of some Hermitian operator, then we are given such an operator on H and we know what to insert in the right Hilbert space. As we discussed before, an operator on H R can be used just as well to define an operator on the left side Hilbert space H So suppose we are given some Hermitian operator on H R . In an AdS/CFT context, in the bulk it would correspond to a real field. Now take the field and move it to the left side. The operator that we get is the operator just defined. For example consider the Hamiltonian on H R , and it has a series of eigenvectors and eigenvalues Then the Hamiltonian as it acts on H † L is where w i is the dual basis of ψ i w i a ψ a j = δ i j (A.8) 11 If we have a preferred basis and use it to define anti-unitary transformation which identifies vectors in H and H † , then the definition above is the transpose in that basis -but we do not need to use such a map in our conventions. So the prescription above indeed keeps the eigenvalues of the Hamiltonian (or any other Hermitian operator).

Time evolution for H † looks a bit unusual in these conventions
Note that the identification (A.5) is such that it identifies O a b (t) = U a b (−t) as expected from the overall invariance of the background when A = 1.

B Motivating restriction to G(A)
In our discussion of the Pauli spin model, we focused on the simple invariant density matrix introduced in (4.8), which is simply the product of invariant density matrices for each of the spins with the same value of A. In this appendix, we consider generalising this to allow different values of A for different spins; this gives an SU (2) N invariant density matrix We can obtain SU (2) N S N invariance by summing over all permutations of the indices with equal weight. So we consider This gives us an N parameter family of SU (2) N S N invariant density matrices. Writing (B.2) more explicitly in terms of states, it is Taking the product, there are N K 3 K states with K of the modes in triplet states. All these states have the same probability p K , We now argue that setting the A i equal maximises the entropy of the density matrix, for fixed strength of the correlation between the two boundaries, providing some additional motivation for this choice.
For simplicity, we will give an explicit calculation where we consider just two possible values forÂ i , that is there are two valuesÂ andB which appear in N A , N B of the species (N A + N B = N ). We define α A,B = N A,B /N . It should be straightforward, although probably messy calculationally, to extend the argument to more general cases.
We want to consider varyingÂ −B holding fixed the "length of the wormhole" in the dual description. We take the size of double-sided two-point correlations as a proxy for this length, as in the bulk, these are related to the length of geodesics between the two boundaries, as seen in (5.45). Let M i be a random boundary operator of length p, according to the definition in (3.7). Consider now the double sided 2-pt function M 1 (t 1 )M 2 (t 2 ) . On the SYK side we can use appendix D to see that in the double scaled limit Thus, our condition is that we varyÂ −B holding fixed We consider a density matrix (B.7) There are N K 3 K states with (N − K) states in the singlet and K states in the triplet. After summing over permutations, all these states enter the density matrix ρ with probability p K =p K F (z), wherẽ (B.8) The entropy for the density matrix is We consider varyingÂ,B, holding W = A N A B N B fixed. Let's consider the variation around some common valueÂ =B =Ĉ. At linear order,Â =Ĉ +α −1 A γ,B =Ĉ −α −1 B γ, and the first derivative of the entropy is We have and dF dz (z = 1) = KN A N F (z = 1). This gives, at γ = 0, ∂p K ∂B = 0 at γ = 0, and hence ∂S ∂γ = 0; γ = 0 is an extremum of the entropy.
To see whether it's a minimum or a maximum, we need to consider the second derivative. First, work out the dependence on γ to second order: ifÂ =Ĉ +α −1 Now we correct our previous formula for the first derivative to After some calculation, we find This gives a negative second derivative, indicating a maximum. This gives a motivation for preferring the simple density matrix (4.8). We should note however that the most general SU (2) N S N invariant density matrix is more complicated than (B.2). In the basis of singlet and triplet states, the general density invariant density matrix is block diagonal, as in (B.3), but the states with K triplets do not form an irreducible representation of the symmetry, so the most general density matrix is not proportional to the identity on the space of states with K triplets. The irreducible representations of S N on the space of states with K triplets are labelled by Young tableaux with up to two rows, and the second row has up to L boxes, where L = min(K, N − K). So there are L + 1 irreps in the space of states with K triplets, and the general SU (2) N S N invariant density matrix is where ρ K,I are the projectors onto the irreps. For N even, this has 1 4 (N +2) 2 parameters p K,I , subject to the overall constraint K I p K,I = 1. 12 We will not explore this more general family of invariant density matrices here; we see no reason to expect these more refined objects have a nice geometrical dual.

C O(N ) invariant density matrices in the Majorana model
We here discuss the symmetry-invariant density matrices for the Majorana model. The effective symmetry for the Majorana SYK model is O (N ). This is the symmetry of rotating the index of Majorana fermions in each space at the same time. One may also require that the density matrix be symmetric under the exchange of spaces R and L.
Representation structure Let us now briefly remind the reader the representation structure of O(N ) (we will use the notations of [37]). We can write the total space as a tensor product of the left and right spaces. Restricting for concreteness to even N , then each Hilbert space is the sum of the two spin representations of O(N ) with different chiralities. We will refer to them as S and S . So in total the doubled Hilbert There is one word of caution, however. Since χ L and χ R anticommutes rather than commutes, χ R also acts and flips the chirality of the left Hilbert space. Instead, (−) F L χ R commutes with χ L and hence acts trivially on the left Hilbert space.
Next consider an SO(N )-invariant density matrix on this space. First recall that In order to build an O(N ) invariant density matrix, we need to sum up all the states in the same representation with same weight. Namely, the general density matrix can be written as a linear sum of projectors on each representation, where ρ i denote the projectors on N + 1 different representations of O(N ) included in the tensor product, while a i are arbitrary positive coefficients, such that Tr[ρ] = 1.

Basis of O(N ) invariant density matrices
Although ρ i comprise a set which spans the possible density matrices, there is more convenient option, which is similar to the Aσσ that we used before. Consider the operator We will see that we can achieve any density matrix (C.4) by taking linear combinations of Q n . First of all, Q 2 is related to a quadratic Casimir 14 C 2 of SO(N ). Recalling that the SO(N ) generators Σ ij are given by 14 The Pauli case G(A)'s can also be written using the quadratic casimir of SU (2).
For the antisymmetric representation of rank 0 ≤ r ≤ N , C 2 is C 2 = N r − r 2 , which means One should also note that Q 2 is the same for rank r and rank N − r representations, and that Q 2 = 0 for the rank N/2 rep, irrespective of self-dual or anti-self-dual representations. Now, since Q itself is an O(N ) invariant operator, this can also be written as a sum of projection operators onto individual representations. Also, since Q has a single left fermion and a single right fermion we understand that it changes the chirality of the left and the right spinor at the same time. We see that Q exchanges the two representations of the same dimensionality which were not distinguished by Q 2 , This means that there is a basis in which on a given representation of SO(N ), and upon diagonalization, Q distinguishes the two representations of the same dimensionality with the sign of their eigenvalues. We now conclude that We used the Hermiticity of Q here. This means that Q n spans the complete basis for the O(N ) invariant density matrix, for 0 ≤ n ≤ N , thus establishing our claim above.

Convenient basis of O(N ) invariant density matrices We have seen that the
Since Q N +1 ⊆ Sp{Q 0 , · · · , Q N }, such a function is always of the form where Tr[ρ] = 1. One can also take a more convenient basis where N k=0 b k = 1 to ensure Tr[ρ] = 1. We can also rewrite this by using where A k = tanh(α k ), so that The choice of A k is arbitrary, but it is usually convenient to take it so that A k = −1+ 2k N . In the limit of N → ∞, the basis of O(N ) invariant density matrices becomes ρ(A), and we can take A continuously from −1 to 1.

D Crash course on chord diagrams
In this appendix we will review the chord diagram techniques for the models introduced in 3, in the double scaling limit.
Consider the moments of the Hamiltonian (this is a re-iteration of (3.12)) Let us now see how the double scaled limit allows us to compute it. We will proceed with the random spin model as it will be useful for later discussion, but the discussion is very similar in the Majorana SYK model, and performed in full in [18]. We plug the Hamiltonian (3.5) into (3.12), in order to get Due to the Gaussian distribution (3.3), the expectation value over the coefficients is given by a sum over Wick contractions. This in turn means that the moment m k is given by all possible traces involving k/2 operator strings σ I , each of which appears twice in m k , as Tr(σ I 1 σ I 2 · · · σ I 1 · · · ). (D.3) Each term in the sum over Wick contractions can be represented using a chord diagram. Let each of the k operators σ I define a node on a circle. Each node is labelled by an index j = 1, · · · , k. We then connect the nodes in pairs, to designate which pairs have identical sets of sites I i . See figure 2 for an example of a chord diagram. Now focus on a specific term in the double sum, i.e, a specific chord diagram (Wick contraction) with some specific choice of indices. We compute the trace by commuting the trace on each of the species. The obstruction to doing so is that some of the sites i appear in more than a single index set I; If two chords do not intersect they contribute Tr(σ I 1 σ I 1 σ I 2 σ I 2 ), (D.4) whereas if they intersect they give a factor proportional to Tr(σ I 1 σ I 2 σ I 1 σ I 2 ).

(D.5)
If there is a non-trivial overlap I 1 ∩ I 2 = ∅ then these factors will be different.
In the double scaled limit described above, [13] showed that there are two major simplifications in computing this sum over chord diagrams: 1. The number of overlapping indices between any two index sets m ij ≡ |I i ∩ I j | is Poisson distributed with parameter p 2 /N .

2.
With probability 1, the intersection of any three index sets vanishes, namely |I i ∩ I j ∩ I k | = 0, for i = j = k. This statement is summarized in lemma (9) there, and subsequent discussion.
These two statements are just a consquence of the double scaling p ∝ √ N , and they hold in both the models described above 15 .
These simplifications imply the following. Consider now the sum over all index sets in some specific chord diagram, namely for some specific pairing, and start carrying out the traces in each of the species's Hilbert spaces.
If a species appears only in a single chord, then the trace contributes 1 3 · 3 a Tr(σ a σ a ) = 1 tom k . Next, due to property (2) above, we can change the variables in the sum over I 1 , · · · , I k/2 to the size of the overlap between pairs of index sets m ij = |I i ∩ I j |, with i, j = 1, · · · , k/2, along with the measure which is the probability of having a given overlap. Due to property (1) above, we know that this is exactly the Poisson distribution with parameter p 2 /N . The N p −k/2 factor precisely turns counting of appearances of a certain type in the sum into probabilities of such events. Each overlap for a given intersection gives a factor of where CD(k) represents chord diagrams with k nodes.
Transfer matrix method Next we will use some linear algebra in order to compute the weighted sum over all chord diagrams (D.9). Consider cutting the circle open at some point and going sequentially along the line. We define the Hilbert space H aux , which is spanned by {|n } ∞ n=0 , along with the diagonal inner product n|n = δ nn . We can think of |n as a state representing n open chords, and a vector in the Hilbert space will be denoted by n≥0 v n |n .
Define T : H aux → H aux the Transfer matrix on H aux . We think of T as acting on a state |n by opening a new chord or closing an existing one, see figure 7. We can reproduce the sum (D.9) if we decide that: 1. T always opens a new chord below all existing chords. This means that chords cannot intersect when they open.
2. Whenever a chord closes and intersects another chord, it does so with a factor of q.

(E.10)
If the last condition is not satisfied on the entire imaginary line, the contour should be indented according to a usual Mellin-Barnes prescription (i.e. separating poles going to right from going to the left). When phrased in terms of gamma and q-gamma functions, the above integral representation is immediately seen to reduce in q → 1 − limit to a Mellin-Barnes representation of the corresponding (undeformed) well-poised hypergeometric 7 F 6 (1), i.e. a Wilson function [39] (up to appropriate Pochhammer factors).