A Random Unitary Circuit Model for Black Hole Evaporation

Inspired by the Hayden-Preskill protocol for black hole evaporation, we consider the dynamics of a quantum many-body qudit system coupled to an external environment, where the time evolution is driven by the continuous limit of certain $2$-local random unitary circuits. We study both cases where the unitaries are chosen with and without a conserved $U(1)$ charge and focus on two aspects of the dynamics. First, we study analytically and numerically the growth of the entanglement entropy of the system, showing that two different time scales appear: one is intrinsic to the internal dynamics (the scrambling time), while the other depends on the system-environment coupling. In the presence of a $U(1)$ conserved charge, we show that the entanglement follows a Page-like behavior in time: it begins to decrease in the middle stage of the"evaporation", and decreases monotonically afterwards. Second, we study the time needed to retrieve information initially injected in the system from measurements on the environment qudits. Based on explicit numerical computations, we characterize such time both when the retriever has control over the initial configuration or not, showing that different scales appear in the two cases.


I. INTRODUCTION
In the past decade, quantum information ideas have become increasingly relevant in high energy physics, especially in connection to the black hole information paradox [1][2][3][4][5]. In this context, a particularly fruitful line of research was initiated by the seminal work by Hayden and Preskill [6], where the authors studied how quantum information is released from a black hole, under the assumption that it is not destroyed during the evaporation process. Their study suggested that information could be released in a time which is much shorter than the black hole lifetime, and related to the time needed for localized information to spread, or scramble, over all the degrees of freedom.
These considerations provided an obvious motivation for a systematic study of information scrambling and the related concept of many-body quantum chaos, also due to the subsequent conjecture by Sekino and Susskind that black holes are the fastest scramblers in nature [7,8]. In turn, this led to the development of several measures of information spreading and chaos, including out-of-time-ordered correlation (OTOC) functions [9][10][11][12][13][14] (historically introduced in the context of disordered superconductors [15]), and the tripartite mutual information defined in Ref. [16].
In this work, motivated by the recent technical advances in the study of RUCs, and inspired by the Hayden-Preskill evaporation protocol, we consider the dynamics of a quantum many-body qudit system coupled to an external environment, where the time evolution is driven by the continuous limit of certain 2-local RUCs. These consist of qudits nonlocally coupled, but with only two of them interacting at a time. This setting allows us to study quantitatively the contribution of the environment and internal dynamics on the scrambling of information. Furthermore, we consider a modified tensor network model with U (1) charge conservation, which evaporates to a unique vacuum state, instead of reaching the maximally entangled state. This provides a more realistic toy model of evaporating black hole in flat space, for which the entropy after the Page time eventually decreases to zero [50,51]. The U (1) charge conservation is an analog of the energy conservation. Pictorial representation of the model introduced in Sec. II. We consider a system S and an environment E consisting of N and M qudits respectively. The evolution is driven by the continuous limit of a random quantum circuit which implements a fast-scrambling dynamics for S with a tunable coupling between S and E. At each infinitesimal time step ∆t a random unitary operator Ui,j is applied to randomly chosen qudits in S with probability p1 = N λ1∆t, while a swap W l,m between a qudit in the system and one in the environment (randomly chosen) is applied with probability p2 = λ2∆t.
In the rest of this paper, we focus on two aspects of the dynamics. First, we study analytically and numerically the growth of the second Rényi entropy of the system, highlighting the implications of conservation laws and the emergence of two different time scales: one is intrinsic to the internal dynamics (the scrambling time), while the other depends on the system-environment coupling. Second, following Hayden and Preskill [6], we study the time needed to retrieve information initially injected in the system from measurements on the environment qudits, and how this depends on the knowledge of the initial configuration of the system.
In the past years, several works have appeared discussing ideas and techniques related to those of the present paper. First, we note that our setting differs from those studied in Refs. [52][53][54][55][56][57][58][59][60] in the context of measurement-induced phase transitions. Indeed, in our model no projective measurement is taken, and we consider instead an environment which is eventually traced over in our calculations. A similar setting was studied in Ref. [61], but there the authors considered random global Hamiltonians, with no notion of local interactions. Next, quantum mechanical evaporation protocols displaying some analogy with our setting were investigated in Refs. [62,63] for an SYK model [11,64] coupled to an external environment (see also [65]). However, the dynamics studied in these works is not Brownian, and is analyzed by means of the Keldysh formalism.
We also mention that very recently the effects of decoherence on information scrambling has been analyzed in Ref. [83] within a quantum teleportation protocol related to the setting of this paper, see also Ref. [84] for an experimental implementation. Furthermore, we note that the Hayden-Preskill protocol with a U (1) conserved charge has been studied before in Ref. [85], where global random unitary transformations (instead of k-local circuits) were considered. Finally, two papers closely related to the present article appeared very recently. First, a random quantum circuit model for black hole evaporation was studied in Ref. [86], but there the authors focused on a different setup and quantities . Second, analogously to our work, the emergence of a Page curve in a unitary toy model for a black hole has also been shown in Ref. [87], based on recently-developed concepts of many-body quantum chaos. However, in this work we focus on a specific microscopic model which is different from the one studied in Ref. [87], and employ different techniques in our calculations.
The rest of this manuscript is organized as follows. In Sec. II we introduce our model, while in Sec. III we analyze the growth of the entanglement both in the case of Haar-scrambled local unitary evolution (Sec. III A) and in the presence of a U (1) conserved charge (Sec. III B). The retrieval of quantum information initially injected in the system is studied in Sec. IV, while we report our conclusions in Sec. V. Finally, the most technical aspects of our work are consigned to a few appendices.

II. THE MODEL
We start by introducing the model studied in the rest of this work, which is pictorially depicted in Fig. 1. We consider two sets of N and M d-level systems (qudits), denoted respectively by S (the system) and E (the environment). The Hilbert spaces associated with S and E are then We anticipate that in our calculations we will always take the limit M → ∞, corresponding to the physical situation where the number of degrees of freedom in the environment is much larger than in the system.
Motivated by the Hayden-Preskill evaporation protocol [6], we would like to construct a quantum circuit which implements a fast-scrambling dynamics for S and with a tunable coupling between S and E. Let us begin by considering a discrete process, and divide the time interval [0, t] into n steps t j = (j/n)t, so that t j − t j−1 = ∆t = t/n. At each time step, the system evolves according to the following rules: 1. with probability p 1 , two qudits in S, placed at random positions i and j, interact. We model this process by the action on h S of a unitary operator U i,j , chosen out of a suitable random ensemble; 2. with probability p 2 ≤ 1 − p 1 , one qudit in S and one qubit in E at random positions are swapped. This models the simplest possible interaction between S and E.
Note that at each time step the system is not evolved with probability 1 − p 1 − p 2 . The random choice of interacting qudits should be considered as "fixed once chosen": as we will see later, this means that when multiple replicas of the system are considered, the circuit is always identical in each copy. The above rule defines a quantum circuit with discrete time steps. It is convenient to take a continuous limit of the former, which allows us to simplify some aspects of the computations. In order to do so, we choose the probability p 1 and p 2 to scale with the time interval ∆t as where λ 1 , and λ 2 are two positive real numbers. Note that while both p 1 and p 2 are proportional to ∆t, they have a different dependence on N . As we will comment on again later, this ensures that the internal time scales are much shorter that those related to the interaction with the environment, as it is assumed within the Hayden-Preskill protocol [6]. With the above choices, expectation values of observables computed at time t display a well defined limit for ∆t → 0 (namely n → ∞), yielding a continuous dynamics for S ∪ E. Importantly, we will be interested in the limit of an infinitely large environment, which will then play the role of a "qudit" reservoir. In the discrete dynamics, it is enough to choose the number M of environment qudits to be M N t/∆t, so that M → ∞ in the continuous limit. In the rest of this work, we will focus on the computation of averaged physical quantities: at each time step this amounts to averaging over all the possible choices of pairs of qudits and of gates U i,j , with the proper probability distribution. For a given fixed time t, this is equivalent to averaging over all the realizations of allowed quantum circuits. A crucial point is that each individual realization corresponds to a unitary evolution. In particular, if the initial state of S ∪ E is pure, it will remain so for any realization, and its von Neumann entanglement entropy will remain zero for all times (and so will its average over realizations).
Finally, regarding the ensemble of two-qudit gates U i,j , we will consider two distinct physical situations. In the first one, the internal dynamics is "maximally chaotic", namely each gate U i,j is drawn out of a Haar distribution. In the second situation, we assume a locally conserved U (1) charge, namely we choose each gate U i,j to preserve the U (1) sectors in the product h S , as it was done in Refs. [39,42] for the case of spatially local RUCs.

III. THE ENTANGLEMENT GROWTH
In this section we study the entanglement growth for a subsystem K ⊂ S, which is naturally quantified by means of the von Neumann entanglement entropy where ρ K (t) is the density matrix reduced to the subsystem K. Denoting by {|j } d−1 j=0 , a basis for the local Hilbert spaces h S ,h E , we will assume that both the system and the environment are initialized in product states, denoted by |Ψ S 0 and |Ψ E 0 respectively. In particular we set, for finite M , while we will consider different initial product states for S. Note that by construction there is no entanglement between S and E at time t = 0. Despite the importance of the von Neumann entanglement entropy, it is known that the latter is difficult to obtain in the setting of RUCs [35]. For this reason, in the following we focus on the related Rényi-2 entropy; more precisely, we will compute S (2) We note that S K (t) is not the averaged second Rényi entropy, as the disorder average is taken inside the logarithm. In fact, Eq. (5) is proportional to the logarithm of the averaged purity P K , which is defined as However, for large N one expects the effect of fluctuations in the disorder to be small, so that the behavior of S K (t) should be qualitatively the same as the averaged Rényi-2 entropy [38].
Let us now define K = S \ K and rewrite where X K is a swap operator exchanging the two copies of K, while in the last line, "tr" represents the trace over the entire Hilbert space. From this expression it is clear that S . In order to compute the latter, it is convenient to recall the Choi-Jamiolkowski mapping which allows us to interpret the operator ρ S (t) ⊗ ρ S (t) (defined on the tensor product of two "replicas" H S ⊗ H S ) as a state in H ⊗4 S . In particular, we define where we introduced the maximally entangled state with In the following, we label with 1 and 2 the Hilbert spaces of the two replicas associated with ρ S (t) in Eq. (8), and with1 and2 the other two. Accordingly, the Hilbert space corresponding to the four replicas is Finally, we also define Within this formalism one can recover the value of the purity using where Eq. (13) can be verified straightforwardly by expanding the scalar product. We note that when the initial state |Ψ S 0 is a product state and invariant under arbitrary permutations of qudits in H S , then the initial state |ρ S (0) ⊗ ρ S (0) , is invariant under permutation of qudits in H S . As it will be clear from the subsequent discussion, this is also true for the evolved state |ρ S (t) ⊗ ρ S (t) : accordingly, the value of the purity P K (t) only depends on the cardinality of K, k = |K|, and not on which sites belong to K and we may write P k (t) = P K (t).
The formalism above allows us to write an equation describing the evolution of the state E [|ρ S (t) ⊗ ρ S (t) ] under the continuous RUC introduced in Sec. II. In particular, in the limit M → ∞, we derive in Appendix A where L is a super operator (the Lindbladian) acting on H S , which reads with In order to proceed further, we need to specify the probability distribution for the two-qudit unitary gates U i,j , which in turn determines the average in Eq. (17). As we already anticipated, we focus on two different physical situations. First we consider the case where U i,j are Haar-distributed over the group U (d 2 ), which corresponds to a maximally chaotic evolution. Second, we consider random gates U i,j with a block structure determined by the presence of a U (1) charge, as done for local RUCs in Refs. [39,42]. The two cases are treated separately in the next subsections.

A. Random Brownian circuit without conservation law
As we have anticipated, we start by choosing the unitary gates U i,j to be Haar distributed over U (d 2 ). In this case, the average in Eq. (17) can be computed easily, and we have (see for instance Refs. [37,38]) where Furthermore, throughout this section we initialize the system in the product state With the above choices, one can now plug the explicit expression (18) into (16) and solve, at least numerically, Eq. (15). Unfortunately, the exact numerical solution to Eq. (15) is difficult to obtain for large values of N , as the dimension of H S grows exponentially with the system size. Luckily, in the present case the problem can be considerably simplified due to permutation symmetry between different qubits, and one does not need to solve Eq. (15) directly. Instead, based on Eq. (13), it is possible to derive the following system of differential equations for n = 0, . . . , N and with the convention P −1 (t) = P N +1 (t) ≡ 0. Here P n (t) is the purity for a subsystem with n qudits, while the initial conditions [corresponding to the state (20)] are P n (0) = 1 , n = 0, . . . , N .
We note that Eq. (21) represents a rare example where an explicit result for the dynamics of the Rényi entropy can be obtained for open systems [88]. Since its derivation is rather technical, we reported it in Appendix B.     It is important to comment on this result. First, we note that setting λ 2 = 0 in Eq. (21), we recover the same set of equations (up to prefactors) that was derived in Ref. [28] for a Brownian Hamiltonian evolution. Thus, the internal dynamics driven by the RUC defined in Sec. II is qualitatively equivalent to a continuous Brownian Hamiltonian evolution. This observation allows us to apply directly some of the results of Ref. [28] to our model. In particular, it was shown in Ref. [28] that the system (21) leads (for λ 2 = 0) to the emergence of a time scale which is logarithmic in N . More precisely, let us call t * (k) the amount of time needed before the purity of a subsystem of size k becomes less than (1 + δ)2 −k , where δ is a small positive real number, and 2 −k is the purity of a maximally mixed state. Then, for 0 < κ < 1 fixed, it was shown that t * (κN ) ∼ ln(N )t * (1) . In our case, due to the choice made in Eq. (1), we have that t * (1) has a constant limit for N → ∞, so that t * (κN ) ∼ ln(N ) for large N . In Ref. [28] this was defined as the scrambling time of the system. Note that, following later developments, the scrambling time is now usually defined as the time needed for OTOCs to decay to zero. However, the latter was shown to be also logarithmic in the system size N for the Brownian Hamiltonian evolution of Ref. [28], see Ref. [33], so that, up to prefactors, they can be identified in our model.
The features of the entanglement dynamics for λ 2 = 0 discussed above are illustrated in Fig. 2, from which the emergence of a time scale logarithmic in N is manifest.
Next, note that for λ 1 = 0 Eq. (21) predicts the purity of any subsystem to remain constant, namely P n (t) ≡ 1 for all values of n. This is due to the fact that, in each realization of the quantum circuit, S remains a pure state, since the evolution only amounts to an exchange of qudits |1 and |0 between S and E. On the other hand, when both λ 1 , λ 2 = 0, the entanglement growth is non-trivial. In the following, we present our results based on the numerical solution of Eq. (21).
In Fig. 3(a) we report the numerical values of S n (t) as a function of the subsystem size n, for different times t and λ 2 = 0. We can immediately appreciate that the effect of the environment is to increase the entanglement of S, even though the environment itself consists of a product state. This is due to following mechanism: if j is a qudit in S, the internal dynamics will generate entanglement between j and S \ j. When j is swapped with a qudit in E, the latter becomes entanglement between S \ j and E. As a consequence, S does not remain in a pure state, and its entanglement grows in time. We also see that the Rényi entropy of K and S \ K is not equal anymore, since the larger of the two can accommodate more entanglement with E.
It is particularly interesting to follow the time evolution of a subsystem K larger than half of the system size, as displayed in Fig. 3(b). We see that there are two relevant time scales that characterize its qualitative behavior: for short times, the Rényi entropy S K (t) starts to increase with a constant slope up to a time t p , at which saturation occurs (the indices s and p stand for "scrambling" and "Page" respectively: the use of these names will be justified in the next section). We can interpret the increase of S   t < t s as mainly due to the internal scrambling dynamics. Based on this picture, we expect t s ∼ ln(N ), while, due to the normalization choice in Eqs. (1) and (2), t p t s for large N .
To verify this, we have computed numerically the time derivative of S K (t), from which the emergence of different regimes is manifest, cf. Fig. 4. We see that for t < t s the derivative is large and increases with N , while for t s < t < t p it approaches a constant s λ2 as N → ∞. It is not straightforward to compute s λ2 directly from Eq. (21): indeed, while at short times the r.h.s of Eq. (21) is dominated by the term proportional to λ 1 , for t ∼ t s the absolute value of the latter becomes comparable to the term proportional to λ 2 and both contribute in a non-negligible way to s λ2 . Nevertheless, we can make the conjecture .
In order to motivate Eq. (23), we consider the case K = S, so that only the term proportional to λ 2 in Eq. (21) is non-vanishing. In the limit λ 2 → 0, one can make the assumption that that after a time t > t s ∼ ln(N ) the system is almost maximally scrambled. Then, for large n one would get P n (t) d −N +n , so that and so d dt S (2) Remarkably, we found that Eq. (23) is in perfect agreement with the numerical solution to Eq. (21) for arbitrary values of λ 1 and λ 2 , and also for general K ⊂ S (with |K| > N/2), suggesting that it should be possible to derive it rigorously from Eq. (21).
We can estimate t s precisely, by defining it as the amount of time needed in order for dS (2) K (t)/dt to become smaller than s λ2 + ε, where ε is a positive small number. We see clearly from Fig. 4(b) that t s ∼ ln(N ), as we also verified with a quantitative fit. On the other hand, one can define analogously t p to be the amount of time needed in order for dS (2) K (t)/dt to be smaller than a small positive constant, and as it is clear from Fig. 4(a), one has t p ∼ N , so that indeed t p t s . In summary, the above analysis shows that in the presence of both internal dynamics and system-environment interaction, two distinct time scales emerge: one can be associated with the internal scrambling time t s , with t s ∼ ln(N ), while the other, t p , depends on the interaction with E, and for the RUC constructed in Sec. II we have t p ∼ N .  In the previous subsection we have seen that the second Rényi entropy for a subsystem K ⊂ S grows always monotonically with time, even if E is initialized in the product state |0 ⊗M . On the other hand, in a unitary black hole evaporation process, one expects that the entanglement follows a "Page-like" behavior in time [51]: namely it initially grows but starts to decrease in the middle stage of the evaporation, and eventually vanishes when the black hole evaporates completely 2 . This difference between a black hole and random tensor networks originates from the absence of energy conservation in the latter. In the long-time limit, the black hole returns to a vacuum state since its energy leaves with the radiation, while the random tensor network model approaches a random state with large entanglement entropy between the system and the bath.
It is difficult to introduce energy conservation in tensor network models, but it is possible to introduce a U (1) charge conservation, which plays a similar role. When the bath is infinitely large and initialized in the zero (i.e. lowest ) charge "vacuum" state, the black hole charge will gradually decrease and approach zero in the final state. As long as the zero-charge state is unique, the black hole entropy will eventually vanish in the long-time limit (for a very different approach that achieves similar phenomena, see Ref. [67]).
We implement a dynamics with a U (1) conserved charge by imposing that the two-qudit unitary gates U i,j have some special structure, as done in Refs. [39,42] for the case of spatially local circuits. For the rest of this section, we will focus on the case of qubits, namely d = 2. Then, following [39,42] we consider gates of the form where the first and last blocks are 1 × 1 and the second block is a 2 × 2 Haar-random unitary matrix. Since the interaction with the environment is driven by swap gates, Eq. (26) defines a dynamics conserving the charge where the charge operator is on each site.
Averaging over unitary gates of the form (26) introduces additional computational difficulties with respect to the case of Haar-distributed operators. In particular, by exploiting the results derived in Ref. [39], Eq. (18) has to be replaced by Here |I ± Qj Q k are states living in the tensor-product of two local sites of the four replica space. In terms of single-site states, they can be written as where, for the case of qubits, the Greek indices take value in {0, 1}, while d 0 = d 2 = 1 and d 1 = 2.
The form of Eq. (29) makes the computations considerably more involved. In particular, one can not derive a set of N + 1 differential equations for the purity, and a different strategy is needed to obtain S K (t) efficiently. Luckily, one can exploit an observation of Ref. [39]. Namely, the states |I ± Qj Q k can be written in terms of the following 6 states [39] |0 ≡ |0000 , Once again, we stress that the states |0 , |1 , |A , |B , |C and |D live in a single local space of the four replicas. This means that the evolution dictated by the averaged gates U j,k effectively takes place in a Hilbert space H eff is a matrix acting on the space C 6 ⊗ C 6 . The above consideration becomes particularly powerful when combined with the underlying permutational symmetry of the operator L and of the initial state |ρ S (0) ⊗ ρ S (0) . Indeed, this allows us to exploit a logic which is similar to the one developed in Ref. [34], and obtain an efficient scheme to compute the evolution of the system in a numerically exact fashion.
We start by introducing the following class of permutationally invariant states on the space H eff S |n 0 , n 1 , n A , n B , n C , n D = 1 √ N !n 0 !n 1 !n A !n B !n C !n D ! π∈S N π |0 ⊗ · · · ⊗ |0 n0 ⊗ · · · ⊗ |D ⊗ · · · ⊗ |D n D π −1 . (35) Importantly, we can rewrite these states by introducing a set of bosonic creation operators as [34] |n 0 , n 1 , n A , n B , n C , n D = 1 where [a j , a k ] = a † j , a † k = 0 and a j , a † k = δ j,k , while |Ω is a vacuum state. One of the advantages of the bosonic representation is that the operator L, the initial state |ρ S (0) ⊗ ρ S (0) and the vector |W K defined in Eq. (14) admit a simple expression in terms of the a-operators. Since we won't make use of them in the following, we report them in Appendix C, to which we refer the interested reader.
Since both the initial state and the Lindbladian L are invariant under arbitrary permutation of qubits, the states (36) form a basis of the Hilbert space in which the dynamics takes place. Crucially, the corresponding dimension is and thus grows only polynomially (rather than exponentially) with N . In practice D perm is still very large for the values of N considered in the previous subsection. Nevertheless, we were able to perform numerically exact calculations       up to N = 80. This was done by implementing the matrix corresponding to L in the vector basis (36), and then computing the evolved state |ρ S (t) ⊗ ρ S (t) solving the system of differential equations encoded in Eq. (15). Note that in this way we did not need to diagonalize exactly the matrix associated with L, which would be unfeasible for N = 80. In rest of this section, we report the numerical results obtained by following the above procedure.
We first consider the case where the system is initialized in the product state (20), and study the time evolution of the Rényi entropy S  Fig. 5. We immediately see that the qualitative behavior is different from the Haar-scrambled case, since S  κN (t) reaches its maximum grows linearly with N/λ 2 . We call t p the Page-time [51]. After t p , we see from  κN (t) decreases and approaches zero exponentially fast, with an exponent that does not appear to depend on κ.
Besides its non-monotonic behavior, S K (t) displays another qualitative difference. Indeed, the initial state of S, defined in Eq. (20), is a fixed point for the internal dynamics. Hence, at short times, one can not distinguish clearly the contribution of the internal scrambling, since the initial growth of S (2) K (t) is only due to system-environment coupling (there is no evolution if λ 2 = 0). For this reason, we consider in Fig. 6 the Rényi entropy S Note that this state is not invariant under permutation of qubits. Accordingly, we consider a protocol where not only we sample over different realizations of the RUC, but we also take an average over all the initial product states obtained by permuting the qubits in (38): namely, we average over all product states where half of the qubits are initialized to |0 , and the rest are set to |1 . It is straightforward to see that in the four replica space H S the state |ρ S (0) ⊗ ρ S (0) (obtained after averaging over the initial configurations) is indeed permutationally invariant, and we can employ the approach explained above. As expected, we see from Fig. 6(a) a separation of time scales for the initial state (38). In order to make this more transparent, and following the previous section, we report in Fig. 6(b) the time derivative of the Rényi-2 entropy S (2) κN (t). Although the results are now plagued by larger finite-N effects, we can see the same qualitative behavior displayed in Fig. 4 for the Haar-scrambled dynamics. In particular, after a time t s ∼ ln N the derivatives appear to approach a plateau, and it remains approximately constant for t s < t < t p , where t p ∼ N is the Page time.
Finally, in order to push further the analogy between our model and a unitary black hole evaporation process, it is interesting to study the time evolution of the system charge The computation of Q S (t) can be carried out using the very same techniques outlined above for the second Rényi entropy. In fact, we note that the calculations are simpler, since they only involve two replicas, instead of four. In particular, it turns out that the average charge can be obtained as the solution to a system of (N + 1) 2 coupled linear differential equations, which can be easily treated numerically. Since no additional complication arises, we omit the details of the computation here, and only report our final numerical results. These are displayed in Fig. 7, where we also show S N (t) as a function of Q S (t). We note that the dimension of the Hilbert space associated with a given integer value Q of the charge is N Q . Of course, the evolved system state will have nonzero projection onto different sectors of the charge at a given time. Nevertheless, we can define an effective "black hole" Hilbert space dimension associated with the averaged charge as We see that the behavior of D BH (t) depends on the initial state chosen for S. If the system is initialized as in Eq. (20), then the effective Hilbert space dimension will first increase, and then decrease after Q S (t) reaches the value N/2.
On the other hand, if the system is initialized in the state (38), the effective Hilbert space decreases monotonically, as one would expect in a more realistic unitary black hole evaporation process. Before leaving this section, we comment on the choice (4) for the initial state of the environment. As we have mentioned, this is motivated by an analogy with the black hole evaporation process, where (4) plays the role of the "global" vacuum state. However, it is natural to wonder what would happen if E is initialized, instead, in a random product state. In this case, a non-vanishing charge would be pumped into the system at each time interval, so one would expect that the qualitative behavior of the Rényi entropy would be similar as in the evolution without the U (1) conserved charge. We have verified by explicit numerical calculations that this is indeed the case. In particular, if E is initialized in a random product state we observe no monotonic decrease of the Rényi entropy after the Page time.

IV. RETRIEVAL OF QUANTUM INFORMATION
In this section we finally discuss how the RUC introduced in Sec. II provides a microscopic model for the informationretrieval protocol studied by Hayden and Preskill [6], and allows us to investigate quantitatively several aspects of the latter. We start by briefly reviewing the setting of Ref. [6], and then proceed to present our results.
We recall that the information stored in a black hole is emitted in the form of Hawking radiation [2], so that one can ask what is the minimum amount of time that is needed before such information can be recollected from   , we initialize the system S (the "black hole") in a product state. A qudit A is initially injected into the black hole, and a third party (C) holds a reference system, namely an ancilla which is maximally entangled with the former at time t = 0. The system is in contact with E (representing the exterior of the black hole, including Hawking radiation) and evolved by the RUC introduced in Sec. II. In the second setting [subfigure (b)], the retriever (B) initially holds a copy of the black hole, namely S is initialized in a maximally entangled state with a set of ancillary qudits. measurements performed outside of the black hole. In order to make contact with our model, we interpret S as the black hole, while E consists of all its exterior degrees of freedom (hence including Hawking radiation). Following Ref. [6], we then imagine that Alice injects a qudit A into the system at time t = 0, and that a third party C (Charlie), holds a reference qudit which is maximally entangled with the former. The system is in contact with E and evolved by the RUC introduced in Sec. II. Finally, we imagine that Bob wants to recover information on the injected qudit by only performing measurements outside of S. Depending on the initial configuration of the system, the ability to faithfully do so after a given time t is captured quantitatively by the mutual information between different sets of qudits, as we now explain.
First, let us consider the setting pictorially depicted in Fig. 8(a): in this case, Bob has no control over the initial configuration of the system S, which is initialized in a given product state at time t = 0. The capability of recovering information on the injected qudit by measurements on E is quantified by the mutual information which tells us how much information can be extracted from the reference qudit C by accessing those in E. In particular, if I (a),[C,E] (t) is close to its maximal value, then Bob can faithfully recover the information initially injected into S. Note that in Eq. (41) we used an index (a) to distinguish the two settings in Fig. 8. As usual, due to computational limitations, in the following we will not compute the quantity in Eq. (41), but rather its Rényi-2 version, namely where S   Fig. 8(b), and increasing values of N . For both plots, the evolution is driven by the maximally chaotic RUC of Sec. II (without conserved charges), where we set λ1 = 1, λ2 = 2 and chose d = 2.
In the second setting, displayed in Fig. 8(b), we imagine instead that the black hole formed long ago, and that Bob has been collecting its emitted Hawking radiation ever since. Accordingly, by the time the qudit A is injected, the black hole S is in a maximally entangled state with the previously emitted radiation system, which is under Bob's control [region B in Fig. 8(b)]. In this case, Bob can also access these qudits, together with those in the environment E, and his capability to recover the initially injected information is quantified by I (b),[C,E∪B] (t). Accordingly, analogously to the previous case, in the following we will compute its Rényi-2 counterpart It turns out that the formalism introduced in the previous section is adequate to compute numerically the mutual information in Eqs. (42) and (43), for both RUCs without and with a conserved U (1) charge. To see this, we can exploit the fact that the Rényi entropy of a subsystem K is equal to that of its complementary one (with respect to the whole space), and rewrite Each individual entropy in the r.h.s. of the above equations can be computed by exploiting the approach in Sec. III A (for the maximally chaotic RUC) and in Sec. III B (in the presence of a conserved U (1) charge). In particular, in each case we can map the problem onto the computation of the time evolution in a four replica space H S , where the dynamics is driven by the Lindbladian operator (16). The only difference with respect to the steps presented in the previous section is in the initial state and purity vector W K |, which have to be modified for each individual term in the r.h.s of Eqs. (44) and (45). Since these calculations do not present additional difficulties, we report them in Appendix D, and in the rest of this section we present our final results. We begin by discussing Fig. 9, where we report data for the maximally chaotic RUC (no conserved charge). Subfigures (a) and (b) correspond to the two different settings discussed above, and display respectively I In both cases, the mutual information has a monotonic behavior, although with qualitative differences. In the first case, it reaches its maximum value in a time which is clearly proportional to the system size N . Interestingly, we see that after a time scale of the order of the scrambling time t s ∼ ln N , the mutual information reaches a small non-zero value, which, however, is seen to decrease with the system size N . We can interpret this as follows: after the scrambling time, Bob is able to only reconstruct a small amount of the initially injected information, and needs to wait for a time proportional to the black hole size in order to retrieve all of it.
Conversely, we see from Fig. 9(b) that the information retrieval is much faster in the case Bob holds a copy of the black hole [cf. Fig. 8(b)]. In particular, from the plot we clearly see that the mutual information reaches its maximum value after a time which is logarithmic with the system size N , namely the scrambling time. We have repeated the same calculations for a RUC with a conserved U (1) charge, and reported our results in Fig. 10. We see that the functions I (2) (a),[C,E] (t) and I (2) (b),[C,E∪B] (t) display the same qualitative features. It is interesting to note that, in the setting corresponding to subfigure (a) of Fig. 8, the value of the mutual information after the scrambling time is larger than that in the maximally chaotic case, although still vanishing for N → ∞. This is intuitive: the presence of conservation laws constrains the Hilbert space that can be explored by the system, hence generally increasing the knowledge on its state.
In energy conserved systems, the Lyapunov exponent measured by OTOC growth generically depends on temperature. Usually a slower scrambling occurs at lower temperature T , and an upper bound of 2πT has been proven for a particular regularized version of the OTOC [13]. The analog of temperature dependence in our model is the charge dependence of the information retrieval time. We expect that when the charge is closer to 0 or 1, the Hilbert space size is smaller, leading to a similar effect as reducing temperature in energy conserved system. For this purpose, we study the mutual information growth for states with different charge. To this end, we first consider the protocol depicted in Fig. 8(a), but we initialize the system S in the product state As we have clarified after Eq. (38), we actually consider averages over all the initial states obtained by permuting different qubits in Eq. (46), namely over all the product states with n qubits initialized to |0 and (N − 1 − n) 12. Pictorial representation of the third settings considered in Sec. IV. We initialize the system S (the "black hole") in a product state, except for a qubit A, randomly chosen, which is maximally entangled with an ancillary one, denoted by C. The system E ∪ S \ A is evolved with the RUC with a U (1) conserved charge for a time t1. After that, E ∪ S is evolved with the same RUC, for a time t2. to |1 (where the last qubit N corresponds to A, and is entangled to the ancilla C). This allows us to exploit the permutational symmetry in the four-replica space, and proceed following the very same steps outlined above to obtain numerically exact data for the mutual information I We report our results in Fig. 11, for different values of n, and a fixed system size N = 60. In subfigures (a) and (b) we report data for decreasing values of the initial charge Q = N −1−n, respectively larger or smaller than N/2. As we have already pointed out, in the former case the effective Hilbert space dimension (40) has a non-monotonic behavior, whereas in the latter case it is monotonically decreasing, as one would expect in a more realistic unitary evaporation protocol. This is reflected in the fact that, at short times, the two plots display a different qualitative behavior as n increases: in subfigure (a), I (a),[C,E] (t) decreases as n increases, while the opposite happens in subfigure (b). In subfigure (c) we report instead the logarithm of the difference between the maximum value 2 of the mutual information and I (2) (a),[C,E] (t) at late times. The plot shows the emergence of an exponential decay, which starts first for smaller initial charge (a larger initial charge takes longer to evaporate).
Next, we consider the protocol reported in Fig. 8(b). In this case, the initial state is given by a maximally entangled state between S and B. This has a non-vanishing projection over all the charge sectors, so we can not vary arbitrarily its charge, as for Fig. 8(a). For this reason, we consider a different setting, which maintains some of its features, but allows us to tune the initial charge of S. This is depicted in Fig. 12. The idea is to initialize the system in a product state, and let the RUC generate an entangled state between S and E. After the Page time, S is approximately a maximally mixed state in a certain charge sector (with charge decreasing in time), as we discussed earlier. At time  t = t 1 , we introduce a new qubit which is maximally entangled with an ancillary one, denoted by C. After that, the dynamics of E ∪ S is dictated by the same RUC for a time t 2 . We are interested in the retrieval of this qubit in E.
Thus we study the mutual information where the index (c) here is used to distinguish this protocol from those in Fig. 8. As usual, we average over the choice of the qudit A: this allows us, once again, to rely on the permutational symmetry in the four-replica space, and exploit the exact same techniques developed so far to efficiently simulate the dynamics (cf. Appendix D 3). We report our numerical results in Figs. 13 and 14, which we now discuss. First, Fig. 13(a) displays the mutual information I  Fig. 10(a)], so that t 1 < t p for the data reported in Fig. 13(a) . In this case, we see that I (c),[C,E] (t 2 ) saturates faster as time increases, which is what we expect. Indeed, for t 1 smaller than the Page time t p , the RUC increases the entanglement between S and E, so a retriever accessing E at time t = t 1 has more control over the configuration of the "black hole" when the extra qubit is injected.
At t 1 ∼ t p , S and E will be maximally entangled within a given charge sector. Thus, the retriever should be able to faithfully recollect information on the injected qubit after the scrambling time t s . However, since the charge is conserved, the portion of the Hilbert space that can be explored during the dynamics is smaller that 2 N . For this reason, we expect t s ∝ log S, where S = ln D BH (t p ) and D BH (t p ) is the effective dimension defined in Eq. (40). Unfortunately, we can not reach large enough system sizes to test this statement quantitatively.
Next, we report in Fig. 13(b) the mutual information I (c),[C,E] (t 2 ) for t 1 > t p and fixed initial charge Q. The plot shows that as t 1 increases the mutual information saturates more slowly, which is due to the fact that the entanglement between S and E decreases for t 1 > t p . In this respect, it is particularly simple to understand the limit t 1 → ∞: in this case the configuration of S at time t = t 1 will be extremely close to the vacuum, and there will be essentially no scrambling of information in S, leading to an extremely slow saturation of I (2) (c),[C,E] (t 2 ) . From Fig. 13(b) we can also extract the dependence of the scrambling time for information injected at time t 1 on the system Rényi-2 entropy at time t 1 , namely S N (t 1 ) . Here the scrambling time t s (t 1 ) is defined as the value of t 2 at which the mutual information reaches the value 2 − δ, where δ is some small positive number. This is reported in Fig. 14(a), where we chose δ = 0.2. From the plot it is clear that t s (t 1 ) is a monotonically decreasing function of S (2) N (t 1 ) for t 1 > t p , as we already discussed above. Finally, Fig. 14(b) shows I (2) (c),[C,E] (t 2 ) for different values of the initial charge Q, for fixed t 1 ∼ t p (Q = 39) (the Page time depends on the initial charge). In this case, we see that I (2) (c),[C,E] (t 2 ) is decreasing with Q, which is what we expect: if the initial charge is small, then the corresponding Page time is short. So, for Q < 39 and a given time t 1 > t p (Q = 39), the entanglement between S and E will be small, leaving the retriever with little control over the configuration of S when the extra qubit is injected.

V. CONCLUSIONS
In this work, we have considered the dynamics of a quantum many-body qudit system coupled to an external environment, where the time evolution is driven by the continuous limit of certain 2-local random unitary circuits. We have shown that the growth of the second Rényi entropy displays two different time scales that are related to the internal information scrambling and the interaction with the environment. Furthermore, we have characterized the qualitative differences that emerge choosing the unitaries to be Haar-distributed or with a conserved U (1) charge. In the latter case, we have shown that the entanglement displays a Page-like behavior in time, where it begins to decrease in the middle stage of the "evaporation". Finally, we have shown that our model provides a microscopic realization of the Hayden-Preskill protocol for information retrieval, studying quantitatively the time evolution of the mutual information between different subsystems. The conserved U (1) charge provides a tunable effective Hilbert space size, and allow us to study the charge dependence of scrambling dynamics.
The RUC considered in this work can be enriched in a number of ways. For instance, we have always considered the limit where the environment has an infinite number of qudits, that are non-interacting with one another. One could wonder whether the qualitative features described in this work are modified by considering an environment with a finite number of qudits, possibly with a non-trivial internal dynamics.
Next, it would be extremely interesting to consider the growth of local operators [9][10][11][12][13][14] in our setting. While the effect of decoherence on the latter has been already considered in the literature [83], our model provides an ideal playground where numerical and analytic results can be derived for large values of N , and the implications of conservation laws explored in detail. We plan to go back to these questions in future investigations.
Finally, when compared to holographic duality, our model gives us a toy model for the boundary dynamics. It would be interesting to use a tensor network approach to describe bulk degrees of freedom, and study the entanglement wedge structure. We wish to write down an evolution equation for the state |ρ S (t) ⊗ ρ S (t) . To this end, we start with the discrete version of the quantum circuit introduced in Sec. II. Choosing a time t j fixed, we focus on an individual realization of the circuit. This defines a global unitary transformation on S ∪ E which we denote by U (t j ). Then, we have (A1) The operator U (t j + ∆t) is obtained from U (t j ) by applying a suitable unitary operator. In particular, according to the evolution described in Sec. II, we have three possibilities: • with probability 1 − p 1 − p 2 no unitary is applied at time t j , so that U (t j + ∆t) = U (t j ); • with probability p 1 a unitary between j and k is applied, so that U (t j + ∆t) = U j,k U (t) • with probability p 2 a swap exchanges one qudit in S and one qudit in E.
We can now take the average over all possible realizations. We note that the average can be taken independently at each time step, so that, due to the above considerations, the r. h. s. of Eq. (A1) splits into the sum of three contributions The first, corresponding to no unitary applied, is trivial Next, C 2 can be easily determined, since the action of U j,k , for j, k ∈ S, commutes with tracing over E. We obtain The term C 3 is more complicated, because it couples the system S and the environment E. However, it can be computed explicitly in the limit M → ∞. Indeed, let us denote by j and k the qudits in S and E respectively that are swapped at time t j . Assuming M N t j /∆t, we have a negligible probability that qudit k in the environment has interacted before with S. Hence, we can assume k to be in its initial configuration |0 k , and hence having no entanglement with the rest of the qudits in E. Under this assumption (which becomes exact in the limit M → ∞), it is straightforward to compute where |I + j was introduced in (10). Putting all together and scaling the probabilities p 1 and p 2 with ∆t and N as defined in (1) and (2) results in with the final Lindbladian (16) as L. In the limit ∆t → 0 we recover the differential equation (15).
Appendix B: Derivation of the system of differential equations for the purity in the Haar-scrambled case In the maximally chaotic case, we do not need to evaluate directly (15) to obtain |ρ S (t) ⊗ ρ S (t) . Instead, we can derive the system (21) of N + 1 coupled differential equations for the purities P n = W n |ρ S ⊗ ρ S (see (13)) for subsystems of size n.
To this end, we insert the Lindbladian (16) into the equation (13) defining the purity, Next, the action of W n | from (14) onto the Lindbladian L from (16) with U i,j from (18) can be computed. Using the identities and keeping in mind that for W n | = W K | only the size |K| = n of the region matters, this results in by considering separately the three sets of terms in the sum 1≤j<k≤N where the j'th and k'th site of W n | consist of I ± | j , I ± | k with the signs +, + or −, − or opposite, respectively. The differential equation (21) for the purities P n (t) then easily follows.
In this section we discuss in more detail the formalism introduced in Sec. III B, and derive a set of formulas that are needed for numerical implementations. We start by showing how to write operators in terms of the bosonic a-operators. First, we notice that one simply has as can be explicitly checked by comparing the action of the two sides on any state. From this, is follows for x, y, z, t = 0, 1, A, B, C, D. One can now prove a general formula, which can be directly applied for implementing the effective Hamiltonians appearing in the main text. Let us consider j<k x,y,z,t Γ x,y |x j ⊗ |y k Λ z,t z| j ⊗ t| k =: ( * ) , where Γ x,y = Γ y,x and Λ z,t = Λ t,z are symmetric matrices. We can rewrite In the second term, the parenthesis that multiplies Γ x,y Λ z,t is antisymmetric under simultaneous exchange x ↔ y, z ↔ t. Since Γ x,y Λ z,t is instead symmetric, the sum is zero. Accordingly, we have where we used Eq. (C2). Finally, we show how to write symmetrized states in terms of bosonic a-operators. For this we consider a general state described by coefficients c i,z , i ∈ {1, . . . , N }, z ∈ {0, 1, A, B, C, D}, which we symmetrize: where c I,z = i∈I c i,z . From the above general formulas, it is now straightforward to rewrite the Lindbladian (16), with the choice (29), in terms of bosonic operators, together with the states relevant for our computations. In particular, we derived j<k U j,k = 1 2 α=0,1,A,B,C,D a † α a † α a α a α + α=A,B,C,D a † 0 a † α a 0 a α + a † 1 a † α a 1 a α + 1 3 a † 0 a † 1 a A a B + a † A a † B a 0 a 1 + a † 0 a † 1 a C a D + a † C a † D a 0 a 1 and N j=1 |0, 0, 0, 0 j I + | j = N j=1 |0 j ( 0| j + 1| j + A| j + B| j ) = a † 0 (a 0 + a 1 + a A + a B ) .
Furthermore, it follows from Eq. (C6) that Eq. (14) can be rewritten in terms of bosonic modes as where k = |K|. Note that Eq. (C9) actually corresponds to symmetrizing over all possible sets K of k elements. This is correct, since we are interested in the overlap (13), and the state |ρ S (t) ⊗ ρ S (t) is invariant under arbitrary permutations.
Finally, let us consider the initial state (20). It is immediate to see that this corresponds to the state In this section we provide all the necessary details to obtain the results on information retrieval presented in Sec. IV. For ease of presentation and numerical efficiency, we restrict to qubits.

Scenario (a)
Let us begin with scenario (a), in which the black hole is in an initial product state except for one qubit A [cf. Fig. 8(a)]. Rather than the initial product state (20) for the system S, the initial state is now an entangled state with the two-replica Jamiolkowski representation After time evolution, we may extract the purities necessary for the mutual information (44) similarly to (13), but the vector |W is now defined on systems S and C. In particular, for the various Rnyi entropies needed, we have In order to perform the calculation of the purities P X (t) = W X |e −Lt |ρ S∪C (0) ⊗ ρ S∪C (0) , (D6) we symmetrize over system S, including the location choice of qubit A entangled to C. Due to the projection onto I ± | C , the sum over s may be restricted to s ∈ {0, 1, A, B, C, D}.
In the maximally chaotic case, we can derive and use the differential equation (21) as in section III A, with initial conditions where K ⊂ S, k = |K| as usual.
For the case with conservation laws, the symmetrization allows us to express all the states in the four replica space in the bosonic formalism as in section III B, within which we can numerically compute the purities. For this, one needs to write down explicitly an expression for |ρ S∪C (0) ⊗ ρ S∪C (0) . To this end, let us generalize the case considered in Eq. (D1), by considering instead the case corresponding to the initial state (46), where we sum over all the possible permutations of qubits. Then, following the technical derivations in the previous section, it is possible to derive |ρ S∪C (0) ⊗ ρ S∪C (0) = 1 4 s∈{0,1,A,B,C,D} Note that here we have N −1, and not N , appearing in the second exponent, because one qubit is maximally entangled with C, so only N − 1 qubits in S are in a product state.

Scenario (b)
Now let us move to scenario (b), in which the black hole is maximally entangled to a retriever B, except for one qubit A, that is maximally entangled to C [cf. Fig. 8(b)]. Here, the initial state is an entangled state which reads The |W vectors for the purities involved in the mutual information (43) are as in (D5) with an additional |I + B,j for each j ∈ B, since B is never within a region we compute the purity of. Therefore we may directly evaluate sj ∈{0,1} N |s j S,j I + |s j B,j = |0 + |1 + |A + |B and use the simplified initial state after restriction of s and symmetrization of S as above. For the evolution with charge conservation, this bosonic formalism is again the basis for our numerical calculations. Note, finally, that in the case of Haar-scrambled evolution, we can again use the differential equation (21), where the initial conditions are now

Scenario (c)
In order to implement the two-step protocol depicted in Fig. 12, it is crucial to remember that the interaction with the bath is Markovian. First, we simply evolve for time t 1 the initial pure (and symmetrized) state (46) of N − 1 qubits in a given charge sector. Then, we add a qubit s maximally entangled to the ancilla. Symmetrizing its position, this amounts to the following change of basis vectors of the system state: |n 0 , n 1 , . . . , n D S → 1 4 s∈{0,1,A,B,C,D} n s + 1 N |n 0 + δ s0 , n 1 + δ s1 , . . . , n D + δ sD S ⊗ |s C .
The rest of the protocol is then analogous to scenario (a) for time t 2 and the initial mixed state above.