Thermalization after holographic bilocal quench

We study thermalization in the holographic (1+1)-dimensional CFT after simultaneous generation of two high-energy excitations in the antipodal points on the circle. The holographic picture of such quantum quench is the creation of BTZ black hole from a collision of two massless particles. We perform holographic computation of entanglement entropy and mutual information in the boundary theory and analyze their evolution with time. We show that equilibration of the entanglement in the regions which contained one of the initial excitations is generally similar to that in other holographic quench models, but with some important distinctions. We observe that entanglement propagates along a sharp effective light cone from the points of initial excitations on the boundary. The characteristics of entanglement propagation in the global quench models such as entanglement velocity and the light cone velocity also have a meaning in the bilocal quench scenario. We also observe the loss of memory about the initial state during the equilibration process. We find that the memory loss reflects on the time behavior of the entanglement similarly to the global quench case, and it is related to the universal linear growth of entanglement, which comes from the interior of the forming black hole. We also analyze general two-point correlation functions in the framework of the geodesic approximation, focusing on the study of the late time behavior.

1 Introduction The description of thermalization and equilibration in closed quantum systems has been a long-standing theoretical problem. An interesting setting to study non-equilibrium physics in quantum systems is a quantum quench, when one prepares the ground state of a given Hamiltonian H 1 and then evolves it unitarily under the action of a different Hamiltonian H 2 = H 1 . The system is going to evolve to some pure state |ψ(t) after time t > 0. The system is said to thermalize if for any subsystem A the density matrix ρ A = TrĀ|ψ(t) ψ(t)| becomes equal to the thermal density matrix with a certain nonzero temperature. The entanglement entropy of the subsystem S(A) = −Tr A ρ A log ρ A and related quantities are very useful tools to probe subsystems which relax to thermal equilibrium as the full system evolves in time after the quantum quench. The AdS/CFT-correpsondence and holography [1,2] have given new tools for studying physics out of equilibrium in strongly coupled QFT. According to the holographic dictionary, thermal state in the boundary theory corresponds to a black hole in the bulk. The problem of studying the behavior of observables during thermalization after a quench in strongly-coupled quantum field theory can be treated in the leading order of the semiclassical approximation using non-stationary asymptotically AdS classical gravity solutions which describe formation of a black hole in the bulk. In holographic context, the entanglement entropy then can be calculated according to the Ryu-Takayanagi prescription [3] generalized to the general time-dependent case by Hubeny, Rangamani and Takayanagi [4], as area of the minimal surface anchored on the boundary region where the subsystem under the study is located. The most well studied holographic model of thermalization in CFT is the Vaidya dust shell collapse [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. This gravity dual models the global quench in CFT [20,21], when in the initial time moment a spatially uniform distribution of energy is injected into the system. Another known holographic model of global quantum quench is the end of world brane model [11]. The entanglement dynamics of the global quench in CFT were found to share many similarities between these two different models [11,17], hinting at the possibility of some universality of entanglement spreading at least in global quench situations.
Another type of quantum quenches is the local quench [21][22][23][24][25][26][27][28][29][30][31][32]. Such quenches have been studied holographically as perturbations of the zero temperature vacuum [22,24] as well as of thermal equilibrium state [26,27,[29][30][31][32][33]. In the present paper we study the thermalization in holographic (1 + 1)-dimensional compact CFT after a particular variation of a local quantum quench. This quench protocol, which we call the bilocal quench, is realized by simultaneously creating two high-energy excitations in the antipodal points on the cylinder. In the bulk the dynamics after this quench are described by the collision of two massless particles in the AdS 3 spacetime which leads to the formation of a static massive particle or a BTZ black hole [34]. We focus on the black hole formation case, which describes thermalization in the boundary theory after the quench.
The main feature of this model, which is not prominent in the Vaidya global quench model of thermalization, is the fact that the AdS 3 spacetime with two colliding particles which create a black hole is explicitly described as a topological quotient space of AdS 3 by a certain topological identification [34][35][36][37]. This simplifies dealing with boundary physical quantities which are expressed geometrically in the bulk as e.g. geodesic lengths by relating all geodesics in the quotient to certain auxiliary geodesics in the global AdS 3 spacetime. The bilocal quench setup is also interesting because the thermalization of a closed system after introduction of two high-energy local excitations is an attractive toy model for description of thermalization in systems such as quark-gluon plasma after a collision of two heavy ions [38].
Our main object of study in the present paper is behavior of the entanglement entropy of subsystems in the boundary CFT in the non-equilibrium regime after the quench described above. We perform the holographic computation of the entanglement entropy and mutual information in different subsystems after the bilocal quench, and we analyze the time dependence and spreading of the entanglement. We make a direct comparison to the global quench thermalization models, in particular the model based on the null dust collapse in the Vaidya-AdS spacetime. We find that because of lack of translational invariance in the initial state, the equilibration picture globally is substantially different from the picture given by the global quench. Specifically, subsystems which do not contain one of the initial excitations inside, exhibit thermal behavior of the holographic entanglement entropy right from the beginning of the time evolution. For subsystems, which do contain one of the excitations, however, the entanglement entropy demonstrates non-trivial non-equilibrium dynamics in many ways similar to the global quench situation, but with some substantial differences. Since the bulk spacetime is explicitly represented as a locally AdS 3 space with a topological identification, construction of HRT geodesics which calculate entanglement entropy becomes a purely geometrical problem. We discuss it in detail and make some observations about the loss of memory about the initial state upon equilibration of subsystems. We also study the leading behavior of two-point correlation functions in the framework of the geodesic approximation [35], including the long-time behavior.
The paper is organized as follows. In the section 2 we first set up our conventions and introduce the basic objects which are necessary for description of the bulk holographic dual to the thermalization after the bilocal quench. Then we describe the geometry of the AdS 3 spacetime with two colliding massless point particles which create a BTZ black hole and explain how it works as a bulk holographic dual. In the section 3 we study the boundary-to-boundary geodesics which are necessary for holographic computations in this bulk spacetime. We classify them, calculate the geodesic lengths and prove several statements about their behavior with respect to topological identifications generated by colliding particles. In the section 4 we use the results of the section 3 applied to the bulk spacetime described in the section 2.2 in BTZ coordinate patch to perform the holographic calculation of the entanglement entropy and mutual information and study the time dependence of entanglement in detail. In the section 5 we continue the holographic study of thermalization by analyzing the two-point correlation functions in the framework of the geodesic approximation. In the section 6 we recollect the results of the work and discuss their implications and future directions. We start the discussion by establishing the conventions and describing the basic objects which will help us then construct the holographic dual for bilocal quench in the boundary. On the gravity side, we deal with the pure AdS 3 spacetime, as well as with asymptotically locally AdS 3 solutions of 3D Einstein equations. Because 3D gravity is topological, solutions of gravitational equations with negative cosmological constant are global defects in AdS 3 . More precisely [39], they have the general form of AdS' 3 /Γ, where Γ -a discrete subgroup of the isometry group SL(2, R) 2 , and AdS' 3 is the subset of AdS 3 where Γ acts discretely. The objects in AdS 3 in which we are interested, namely point particles and black holes, are particular examples of such solutions. The AdS 3 spacetime has a simple geometry, which allows to use a unified framework to describe global defects in AdS 3 .
We begin with the description of pure AdS 3 space as a hypersurface in the 4dimensional flat spacetime R 2,2 . It is given by the quadratic equation (we set the AdS radius to 1): This quadric surface can be parametrized by coordinates, to which we will refer as global coordinates on AdS 3 : 2) x 1 = sinh χ cos φ , x 2 = sinh χ sin φ , The induced metric on the AdS 3 is given by Here χ ∈ [0, +∞) is the holographic coordinate, and other coordinates have ranges φ ∈ [0, 2π) and τ ∈ [−π, π]. The conformal boundary of the AdS 3 spacetime is located at χ → ∞. The spacetime can be visually represented as a cylinder together with its interior, where χ plays the role of a radial coordinate, τ is a coordinate along the vertical axis of the cylinder and φ is the angular coordinate. The global coordinates are most suitable for description of global defects in the bulk, since they keep the complete information about the topological identification associated with the given defect.
Another coordinate system which we will use is obtained by parametrization: x 0 = − r R cosh R ϕ , x 2 = r R sinh R ϕ , where t, ϕ ∈ R are the coordinates on the boundary, and r ∈ (R, +∞) is the radial coordinate. There is a coordinate singularity at r = R > 0. The metric in these coordinates has the form ds 2 = −(r 2 − R 2 ) dt 2 + dr 2 r 2 − R 2 + r 2 dϕ 2 , (2.5) We will refer to these coordinates as BTZ coordinates. In the AdS 3 spacetime, this patch covers only a part of the global AdS 3 . Note that the choice of the BTZ coordinate system is ambiguous. This ambiguity is the choice of the part of the global AdS 3 spacetime to cover with a BTZ patch. Namely, different choices of the patch can be implemented by changing the signs in front of the square root in the parametrization formulas.
As we will see, the BTZ coordinates are most natural for the holographic description of thermalization in CFT on a cylinder. However, description of topological defects in AdS 3 is more convenient in the global coordinates. Hence we will need the transformation formulas from the global coordinates to BTZ patch. They are obtained using the embedding coordinate parametrizations (2.2) and (2.4): To deal with classical solutions, which are quotients of AdS 3 , it is most convenient to use the algebraic representation of AdS 3 . The AdS 3 spacetime can be described as the SL(2, R) group manifold. We can treat points in AdS 3 as matrices: where the matrix basis is introduced In this notation the condition det X = 1 then gives the hypersurface equation (2.1). The physical quantities on the boundary which we are interested in are calculated from geodesics in the bulk. To study geodesics on quotients of AdS 3 , it is most convenient to work in terms of matrix notations. The geodesics in AdS 3 embedded into R (2,2) can be described as solutions of the Lagrangian [40,41]: where λ is a Lagrange multiplier. The geodesic length in AdS 3 can be expressed in terms of the scalar product in the embedding spacetime R (2,2) . Suppose that x and y are two points in the embedding space, and we denote their respective matrices defined according to (2.7) as X and Y . Then if points x and y belong to the AdS 3 hyperboloid, i. e. det X = det Y = 1, then it is true that The length of a spacelike geodesic between points x and y is expressed by formula and length of a timelike geodesic is given by the formula The isometry group of AdS 3 is the group SO(2, 2) SL(2, R) × SL(2, R), which acts on matrix X as follows: This group has an P SL(2, R) = SL(2, R)/Z 2 subgroup which corresponds to isometries which leave the origin of AdS 3 (which is represented by the unit matrix) fixed. It is realized by choosing u = g −1 = h as an element of the SL(2, R) group up to an overall sign. Then it can represent an isometry of AdS 3 which preserves the origin by acting on X via conjugation: Point-like objects in AdS 3 such as particles and black holes are obtained from empty AdS 3 by taking a topological quotient. The identification is defined by the isometry u acting on the SL(2, R) group manifold via conjugation (2.14). We will refer to the identification isometry u as the holonomy of the topological defect, in agreement with the discussion in [34]. Let us now proceed to concrete discussion of topological defects which we deal with in the present investigation.

Massless point particles in AdS 3
Our main ingredient for constructing the bulk spacetime is a couple of massless particles. A point particle in (2 + 1)-dimensional gravity produces a defect, which holnomy is determined by the momentum vector of the particle [42]. The most general form of a holonomy of a point particle with momentum p µ is given by B. Figure 1. A. Cartoon of propagation of a massless particle in the AdS 3 spacetime projected onto synchronous time slices of the AdS 3 cylinder in different moments of time. The particle moves from left (right) to right (left), and the spacetime is obtained by cutting out the wedge behind (in front of) the particle between the surfaces W ± . B. 3D plot of massless particle moving through the AdS 3 spacetime in global coordinates. The intersection of surfaces W ± is the worldline of the particle.
The condition det u = 1 then means that In the present work we focus on the case of massless particles, which means that we have to set p 2 = 0. Then from the equation above we have 1 u = 1. Thus, the holonomy of a massless particle is given by It produces an identification in AdS 3 , which glues together two surfaces W − and W + : and these surfaces intersect along the worldline of the particle, see Fig.1. The surfaces W ± intersect time slices of AdS 3 along the equal-time geodesics which we denote as w ± (see Fig.1A). Thus in the AdS 3 spacetime the massless particle cuts out the wedge between surfaces W − and W + . One can cut out the wedge either in front of the particle worldline, or behind the worldline. Or, equivalently, one can think that on the Fig.1 the particle moves either from left to right, or from right to left. Sometimes we will call the region space which is cut out by the identification, i. e. the complement of the fundamental domain to the global AdS 3 , as the dead zone, and we call the boundary of the fundamental domain as living space. The holonomy of a massless particle belongs to the parabolic conjugacy class, since |tr u| = 2 in this case. The fixed point of a parabolic holonomy is on the boundary of the H 2 [43]. That means that a massless particle can actually reach the boundary of the AdS 3 spacetime, and the turning points of its worldline at τ = π 2 + πn are located there. The motion of the particle is periodic with return points located at the boundary.

Maximally extended BTZ black hole
The bulk dual for the thermalization process in the CFT 2 is the creation of the BTZ black hole in the bulk. More specifically, in the present work we consider formation of the static BTZ black hole from point particle collisions. The black hole formed from matter in a dynamical process of some kind is dual to a pure state on the boundary, in contrast to the eternal (maximally extended) black hole, which is dual to the mixed thermal state on the boundary [44]. However, we will use the eternal black hole geometry as a reference point for description of the black hole formed from particle collisions.
The maximally extended BTZ black hole in global coordinates is described as an AdS 3 space quotient by a hyperbolic SL(2, R) element. We focus on the static case. The corresponding holonomy has the following form [34]: Here µ > 0 is a parameter related to the mass of the black hole. This holonomy generates an identification which identifies two surfaces V ± : In global coordinates these surfaces are defined by equations [34]: These surfaces intersect AdS 3 time slices along equal-time geodesics v ± , as shown in Fig.2. The maximally extended BTZ spacetime is defined as the region of the AdS 3 spacetime between the surfaces V ± . The part outside of this region is the dead zone which is cut out from the spacetime. The surface V + and V − do not intersect, except for τ = nπ, n ∈ Z, where they intersect along the horizontal diameter of the time slice disc. The spacetime is singular in these moments of time. The spacetime manifold as the region between the two surfaces has two boundaries. Holographically this means that the maximally extended BTZ black hole is dual to the thermofield double state on the boundary. The spacetime has an event horizon, which consists of two surfaces described by the equation cos φ tanh χ = cos τ . (2.22) Spacetime splits into four regions, and for τ ∈ [−π, 0] has the same global causal structure as the maximally extended AdS-Schwarzschild black hole, i.e. we have two external regions, to the left and to the right of both horizons, which are causally completely disconnected. At τ = −π we have the past singularity of the spacetime, and at τ = 0 we have the future singularity. Each of these two regions can be covered by a BTZ coordinate patch, with metric given by (2.5) with horizon located at R = µ π . The action of the holonomy (2.19) results in the identification ϕ ∼ ϕ + 2π. The length of the horizon equals 2µ = 2πR. The horizon radius is related to the mass of the black hole by the relation where G is three-dimensional Newton constant. The Hawking temperature of the black hole which equals the temperature in the dual theory is given by (2.24)

Black hole creation from particle collisions in AdS 3
Now let us discuss the picture of massless point particle collisions. We begin by setting up the stage in global AdS 3 as presented by Matschull in [34]. After that, we will make a transition to the BTZ coordinate patch in a similar way to [36], which gives the natural dual description of the thermalizing CFT on a cylinder. We will need both pictures for our analysis, the first one containing all the data we need for our holographic computations, and the second one for straightforward definition of holographic observables and temporal evolution in the boundary theory after the quench.

Black hole creation in global coordinates
An AdS 3 quotient spacetime contains a black hole if its total defect holonomy belongs to the hyperbolic conjugacy class. i.e. that it coincides with the BTZ holonomy (2.19) up to a coordinate transformation. The total holonomy of two colliding massless particles is a product of two holonomies of each particle. This product holonomy is not neccessarily hyperbolic, but it depends on the energies of the particles. The black hole creation threshold is thus can be expressed [34] as the condition: This translates into a lower bound of energy for colliding particles, which in itself produces a lower bound for the energy of excitations in the bilocal quench which would thermalize the boundary CFT. Since we are interested in thernalization after the quench, we consider only those particle collisions which create black holes. Thus the topological identification in the spacetime is constructed in such a way that we have two singularities with holonomies u 1 and u 2 corresponding to massless particles, and the total holonomy when circling around both particles must equal the BTZ black hole holonomy (2.19). The resulting spacetime can be obtained by making additional cutting and gluing in the maximally extended BTZ black hole spacetime.
More specifically that means that we have to choose two holonomies for particles u 1 and u 2 such that their product would equal u BTZ . The choice of holonomy of a massless particle is dictated by the choice of its momentum vector, according to (2.17). Suppose that two massless particles start from points φ = θ and φ = 0 from the boundary at τ = − π 2 . they move along the radial worldlines given by the equation At the moment of global time τ = 0, particles meet each other at the origin χ = 0, and the collision happens.
We can choose one holonomy, say u 1 , freely, and let the product constraint determine the other one. Since we can multiply in two different orders, there will be two possible choices for the holonomy of the second particle 2 : From these equations, one can find the parameters of the particles (see [34] for more details). We choose the momentum of the reference particle 1 such that it has the energy tan 1 and it moves along the radial direction, starting from the point φ = 0. The corresponding holonomy, from (2.17), reads (2.28) The second particle which starts from φ = θ will have the energy tan 2 . The equation (2.27) then dictates that the particle 2 moves with along the radial geodesic with angle φ = ±θ, where sin θ = tanh µ. It has the energy tan 2 = cosh µ coth µ 2 , and the first particle moves along the geodesic has the energy tan 1 = coth µ 2 . The last expression appears in many formulas in this work, so we introduce the notation: The resulting holonomy of the second particle reads Henceforth all parameters of the infalling particles are determined through the holonomy from the black hole mass parameter µ. Let us now recollect the kinematic data of the particles in the BTZ black hole creation process in global coordinates: • Particle 1: energy tan 1 = coth µ 2 , angle φ = 0 ; • Particle 2: energy tan 2 = cosh µ coth µ 2 , angle sin φ = ± tanh µ.
Having defined the holonomy of the particle 2 as a product of the other particle inverse holonomy with the black hole holonomy, we now can try to represent the geometry of the identification by this holonomy through the identifications corresponding to u 1 and u BTZ . One can show [34] that the second particle sits precisely on the intersection of a wedge face of the particle 1 with an identification geodesic of the BTZ black hole. The choice of the sign corresponds to the choice of the copy of the second particle with A.

B.
C. D. respect to the isometry of the first particle, plus corresponding to the copy located on the W + face, and the minus sign corresponding to the copy located on the W − wedge face. The holonomy of a defect can be thought of as an action of the identification which one encounters when moving along the closed contour with the defect inside. For example, when considering a time slice of AdS with a single particle, the identification cuts out a wedge with faces W ± (see Fig.1), which are identified by the action of the holonomy: The surfaces W ± are given by equation [34]: The forming BTZ black hole is represented by another holonomy, which identifies two surfaces V ± , which are described by the equation (2.21): Using these identifications, we can represent the action of u 2± defined as composition of the upper two holonomies by (2.27). That means that once we circle around the particle 2, we have to go through the identification V − → V + and through W + → W − (note that the u 2 enters in (2.27) as inverse) for any closed contour which lies inside a time slice and contains only the second particle. For the spacetime with two particles colliding into a black hole, we have to impose two more analogous requirements. So, if we circle along a contour containing both particles, we have to pass through the BTZ identification V − → V + . If we circle along the contour around the particle 1, we have to pass only through the identification W − → W + . Combining these requirements, one arrives to the conclusion that the geometry of identifications in the black hole rest frame looks as illustrated on Figs. 3 and 4.
Having described the collision picture in global coordinates, we now have to make some important remarks. First, the black hole which is formed in the collision is not an eternal one, it has only one external region with respect to the apparent horizon. The boundary of this spacetime has only one connected component, which holographically means that this black hole is dual to a pure state in the boundary, as expected in our quench scenario. This situation is very similar to black hole formation from the cloud of collapsing dust in the AdS-Vaidya metric. However, while in the latter case the pure state black hole is usually illustrated through a Penrose diagram, we have the precise picture on Fig.3 of the full spacetime in global coordinates similar to the Penrose diagram of a pure state black hole, but with more detail because we have no spherical symmetry. In particular, the future singularity in a Penrose diagram would correspond to the moment of the collision of particles τ = 0, when the spacetime in global coordinates shrinks into a singularity, see Fig.3D. However, we emphasize that there is much more information contained this picture than in Penrose diagram, because our spacetime is not just a cut out piece of AdS 3 , but a topological quotient. This simplifies the holographic calculation procedure, and yields some interesting details. For example, while the second external region of the BTZ black hole never becomes a part of the spacetime in the collision process and remains inside of the identification dead zone, it actually influences the behavior of holographic observables. This phenomenon will be pointed out precisely when we will discuss HRT geodesics which govern the behavior of holographic entanglement entropy. Another important point is that from the bulk point of view it is most intuitive to perform the diagnostic of black hole formation in the center of mass reference frame [34], where particles start from the opposite sides of the AdS 3 cylinder and move towards each other head-on. Unlike the black hole rest frame, the center of mass frame picture also covers the case of low energies, when a static massive particle is created instead of the black hole. However, the resulting holonomy in that picture (in the high energy case) is not equal to a holonomy of a static BTZ black hole, but it is related to u BTZ by a conjugation, which corresponds to the coordinate transformation from the black hole rest frame to the center of mass frame. However, while intuitively attractive, the center of mass picture is not a natural gravitational dual to the CFT 2 on a cylinder, because the living space is changing with time as wedges move, and it cannot be mapped straightforwardly to a cylinder by a simple coordinate transoformation, unlike the black hole rest frame picture. Nevertheless, the bulk spacetimes with defects which make the living space time-dependent were also studied in the context of holography, e.g. in case of a single moving particle [35,[45][46][47], colliding massless particles in center of mass frame [35,37], moving particles which orbit around the origin of AdS 3 [48].

Colliding particles in BTZ coordinates
We are finally ready to discuss the direct holographic dual to the CFT on a cylinder which equilibrates after the bilocal quench. We make transformation from the global coordianates to BTZ coordinates introduced in section 2.1. The collision of particles in BTZ coordinates in AdS 3 was discussed previously in holographic context in [35] and in context of near-horizon dynamics of black holes in [36].
The transformation formulas from global coordinates to BTZ coordinates are given by equations (2.6), and the metric is given by (2.5). We set the radius of the coordinate horizon R = µ π and we will express all quantities appearing from this point in terms of R, since it is proportional to the temperature. In this case one can show that the surface r = R coincides with the part of the surface of the horizon of maximally extended BTZ black hole given by equation (2.22) which bounds the patch covered by our parametrization in BTZ coordinates. Hence we see that the horizon of the black hole which is about to be formed will be located at r = R. The initial time slice in global coordinates τ = − π 2 , when particles start from the boundary, is mapped into the t = 0 time slice in the BTZ coordinates. Likewise, the t = ∞ time slice in the BTZ coordinates coincides with the horizon surface given by eq. (2.22). The embedding of finite-time BTZ time slices in the global AdS 3 cylinder is shown on Fig.5. Thus, the BTZ coordinates cover a patch which is outside of the horizon of the black hole, or in our case is outside of the apparent horizon of the colliding particles.
Now have to answer the question: how all the cutting and gluing on global AdS 3 performed in the previous section is reflected in the BTZ coordinate patch? In the global coordinates we have two sets of identification surfaces: BTZ identification V − → V + defined by equations (2.21) and the particle 1 wedge W − → W + . The BTZ black hole identification in BTZ coordinates leads to the periodic condition for the angular coordinate: This identification does not depend on the BTZ coordinate time t, and thus the living space in these coordinates is a cylinder, which is exactly what we want. We use transformations (2.6) to get the initial data for particles in the BTZ coordinates. Dividing Figure 6. Creation of the black hole by colliding particles in BTZ coordinates. Red surfaces are the faces of identification W ± introduced by colliding particles. As particles move towards one another, they asymptotically approach the black hole horizon, showed as the black cylinder (which is obscured behind the identification surfaces) in the center. Taking the boundary limit of this formula and substituting the above angle values for both particles in the global coordinates, one gets ϕ 1 = 0 for the particle 1 and ϕ 2 = ±π for the particle 2. We choose the fundamental domain in the BTZ coordinates as ϕ ∈ [−π, π]. In this case the identification in global coordinates φ = −θ ∼ φ = θ translates precisely into the identification ϕ = −π ∼ ϕ = π. Thus we arrive at a picture where two particles move towards each other head-on from the antipodal points of the cylinder, particle 1 moving along the ϕ = 0 and particle 2 moving along the ϕ = π ∼ −π worldline [36], see Fig.6. What is left is to describe the wedge cut out by the particle 1. To derive the equations for its faces, we transform the equation (2.32) in global coordinates to BTZ coordinates, which is the following: tanh χ sin(arctan(coth µ 2 ) ± φ) = − sin τ sin(arctan(coth µ 2 )) ; (2.36) We can expand the sine of sum into two terms and use the formula (2.35) for one term and an analogous formula for the second term. We arrive at the following expression for the W ± wedge faces in the BTZ coordinate patch (our choice of coordinate differs from that of Jevicki and Thaler [36] by a rescaling): These identification surfaces are anchored onto worldlines of particles given by equation and the horizon is located inside the dead zone between W − and W + . The 3D picture of particles moving towards each other in BTZ coordinates is shown in Fig.6, and the cartoon of time evolution is shown on Fig.7. Note that from the equation (2.39) it follows that particles cannot reach the horizon in finite time. This agrees with earlier observation that we will not see the emergence of the horizon in the BTZ coordinate picture in any finite time. Holographically, this means that the state in the dual theory will always remain pure. We conclude this section with a brief discussion of another property of thermalization which is prominent in our holographic description. A common feature of thermalization in a closed quantum system after unitary time evolution of a certain pure state is that at late times the system loses memory about any particular details of the initial state, keeping only information about the extensive characteristics of the initial state, such as total energy, total conserved charge, etc. In our case, these "details" are the initial locations and energies of the individual particles. From the shape of the identification wedge shown on Fig.6 we see that at late times the shape of the identification wedge gradually approaches the cylindrical shape, and cusps at the particle worldlines become smoother and smoother with time. One could say that as the time goes by, the bulk spacetime geometry gradually forgets about the parameters of the particles themselves. The thermal state is recovered in the limit of the infinite time, when the wedge completely falls onto the horizon, and complete rotational symmetry is restored. This kind of memory loss is not so evident in the global quench holographic duals, because by definition in the global quench scenario one deals with a translationally invariant initial state. That means that the details of the initial state are already smeared over the entire boundary time slice from the beginning (however the memory A.

B.
C. loss still absolutely can be observed in the time evolution of physical quantities such as HEE [8,9,19], which we will discuss later). We are certain that one could find the same property for situations where more than two particles create a black hole, or even more complex scenarios of thermalization with a non-homogeneous initial state.

Geodesics in AdS with colliding particles
To proceed with investigation of dynamics of entanglement and two-point correlation functions in the boundary dual of the AdS 3 spacetime with colliding particles, we need to study the geodesics in this spacetime with the endpoints located on the boundary. Generally speaking, we have a (locally) asymptotically AdS 3 spacetime with two conical singularities, moving along the lightlike worldlines. Geodesics can go from boundary to boundary directly, or they can wind around one defect, or around both of them.
To avoid possible confusions, we will refer to the geodesics of the first kind as direct geodesics, to the second kind as crossing geodesics (the meaning of the name will be clarified in the further discussion), and to the third kind as winding geodesics. The main task which we address in this section is to find all geodesics between two given boundary points in BTZ coordinates of colliding particles background, to calculate their lengths and to analyze what happens to geodesics when we evolve the system, that is move the boundary points along the time direction.
To calculate the lengths of geodesics, it is most convenient to use the SL(2, R) group formula (2.11). We are interested in geodesics between boundary points a and b. These points are parametrized by SL(2, R) matrices according to (2.7), where points in embedding space are parametrized by BTZ coordinates using (2.4): We set r a = r b = r 0 >> R as radial cut-off near the boundary. We'll also introduce the auxiliary notation which we will use throughout the rest of the paper: Using the formula (2.11), one can now find the length of a direct geodesic between spacelike-separated points a and b. In the limit r 0 R, it will have the form All holographic quantities which we consider in this paper are expressed through lengths of specific geodesics. However, in the length formula itself (3.4) there is no account for actual existence of the geodesic in the spacetime, since this formula is native to pure AdS 3 and not to a specific topological quotient which we consider in this paper. In order to make any holographic calculations correct in such spacetimes, one has to add to the length formulas the data about the interaction of geodesics with topological identifications. In the rest of this section, we are focusing on this issue in the case of AdS 3 spacetime with colliding particles described in the previous section. We will be using the parametrizations of geodesics in different coordinate systems, as well as isometry formulas from Appendix A and B.

Direct geodesics
It is known that between two given points at the boundary in the BTZ coordinates of the BTZ black hole spacetime one can construct one direct geodesic and an infinite number of geodesics which wind around the horizon (see e.g. [12,58,60]). Once we introduce the infalling particle topological identification described by the holonomy u 1 , some of those geodesics will cross the identification wedge in some manner. The subject of this subsection is to explain which points on the boundary can be connected by direct geodesics in the geometry described in sec. 2.2.2. In BTZ coordinates the identification wedge bisects the initial time slice, hence the geodesics between endpoints located on the same side to the collision line will behave differently compared to geodesics between the endpoints located to different sides of the collision line. Since the topological identification is realized by an isometry and the identification surfaces intersect the time slices along the pieces of boundary-to-boundary geodesics themselves, we can prove some statements about the behavior of the geodesics.
Proposition 1. Suppose that a and b are spacelike-separated points on the boundary such that either ϕ a , ϕ b ∈ (−π, 0), or ϕ a , ϕ b ∈ (0, π). Then there always exists a direct geodesic between these two points at any given moment of time and any time separation.
Proof. To prove this proposition, we observe that the geodesics w ± in BTZ coordinates shown on the Fig.7 are themselves parts of equal-time boundary-to-boundary geodesics. A direct geodesic between two given boundary endpoints can cease to exist if it somehow reaches w ± . The depth of the geodesic, i. e. minimum radial distance from the origin to the geodesic in the bulk, is given by the formula (B.11): where we set n = 0 since we are only interested in direct geodesics. We can directly compare the Γ + to the distance r ± from origin to w ± , which we can determine from the equations (2.38). Solving in terms of the radial variable, one gets (3.6) The r.h.s. is minimal at ϕ = ± π 2 , which is indeed clear from the symmetry of the wedge, see Figs.6-7. It gives the distance First, let us restrict ourselves to the case of equal-time geodesics, ∆t = 0. In this case r ± equals Γ + for ∆ϕ = π. Since by assumptions of the proposition we consider only points in upper or lower parts of the boundary, ∆ϕ < π, we have r ± < Γ + for all equal-time geodesics, as shown on Fig.8 in case of t = 0, when the wedge takes up the most space in the bulk. Now we only have to prove that r ± < Γ + for non-equal-time geodesics. From equations of geodesics (B.2-B.4) it is clear that a geodesic reaches deepest into the bulk in the moment t 0 = 1 2 (t a + t b ). Therefore it is this moment which makes sense as the edge case in (3.7) when a geodesic could possibly try to reach r ± . Further, Γ + depends on tanh R∆t 2 , whereas r ± depends on tanh Rt. It is true that t ≥ ∆t 2 , and the inequality is saturated when either t a or t b is zero. Both r ± and Γ + are decreasing functions of their respective temporal arguments, and their initial values coincide if ∆t = 2t 0 . The question that remains is how fast these functions decrease with time compared to each other. We have to compare two functions, which are proportional to the inverse of (3.5) and (3.7), respectively: where x = tanh 2 Rt ∈ [0, 1] is a variable and y = tanh 2 R π 2 < 1 is a constant. Obviously f (0) = g(0), but for x > 0 we have f (x) > g(x), since one can expand g(x) around x = 0 as follows: (both functions are monotonic). Thus we conclude that r ± < Γ + for any direct geodesic, and the proposition is proved . Now let us consider the situation when the direct geodesics not always exist, namely when the boundary segment between the endpoints crosses the collision line of particles. Suppose that a and b are spacelike-separated points on the boundary such that ϕ a , ∈ (−π, 0) and ϕ b ∈ (0, π), and ∆ϕ < π (∆ϕ > π). Then by definition the direct geodesic between points a and b exists if the geodesic does not intersect the identification wedge. When varying t a and/or t b , we observe that the edge case is when it intersects the worldline of the particle 1 (particle 2). A.
B. We have established the conditions of when the direct geodesics exist and when they do not. Further we will elaborate more on this in case of equal-time geodesics. For now we conclude this subsection by giving the expression of the regularized length of direct geodesic from (3.4). Taking the trace, one obtains the expression (B.15) with n = 0: This expression is valid only for ∆ϕ ≤ π. For ∆ϕ > π, the direct geodesic and the n = −1 winding geodesic change places, so in that case the length of the direct (minimal) geodesic is given by

Crossing geodesics
The crossing geodesic consists of two pieces of geodesics going from the endpoints to the identification surfaces. More specifically, suppose we have the endpoints located as follows: ϕ a ∈ [−π, 0), ϕ b ∈ [0, π), and t a , t b chosen in such a way that points are spacelike-separated. In this case the geodesic will consist of two pieces: a geodesic from a to point o ∈ W − and a geodesic from o * ∈ W + to point b, see Fig.9B 3 . The surfaces W ± are topologically identified, which is represented by the isometry, which acts according to the rule * : X → X * := u −1 1 Xu 1 ; (3.12) We will also need the inverse identification isometry, which is defined as follows: These isometries act on the identification surfaces W ± as follows: * : This enables us to use the geodesic image method [37,46,47] to find the length of the crossing geodesic. Since * is an isometry, we have On the other hand, we can define the inverse isometry # : Therefore, the length of the crossing geodesic can be found as or, equivalently, Thus, the length of the crossing geodesic between points a and b is equal to the length of the image geodesic from a to b # or from a * to b, and the crossing geodesic can be completely recovered from image geodesics and the identification surfaces. These image geodesics themselves are just regular geodesics in global AdS 3 spacetime. In our case of colliding massless particles in AdS 3 , we illustrate the behavior of image geodesics 4 in the black hole rest frame picture on the Fig.9B (BTZ identification is not shown on the picture). All information about the matter which produces the topological identification is encoded in the position of image points a * and b # . These points themselves can be generally located anywhere in AdS 3 . Because of this, one has to exercise caution when working in any coordinate patch which does not cover the AdS 3 spacetime globally, such as the BTZ coordinate patch. While the endpoints a and b belong to the BTZ coordinate patch, the image points, generally speaking, do not. Since the image points generically do not belong to the BTZ coordinate patch, we are in a tricky situation: while the pieces ao and o * b of a crossing geodesic do lie within the BTZ patch, image geodesics as a whole do not. However, image geodesics are the most convenient way to describe crossing geodesics, and we can make use of this machinery in global coordinates to prove some facts about crossing geodesics in BTZ coordinates.
We begin again with establishing the conditions of existence of a crossing geodesic is that image geodesics must intersect (different) identification surfaces. This ensures that the actual crossing geodesic will indeed run around the defect. Keeping this in mind, one can formulate some useful statements. Let us note that the above discussion of image geodesics representation for crossing geodesics is valid only for points located on different sides of the boundary. The following statement says that it is, in fact, the only set of situations when we encounter crossing geodesics.

Proposition 2.
There are no crossing geodesics between the endpoints located to the same side of the identification wedge.
Proof. This statement is not completely obvious in BTZ coordinates, but it is almost trivial in global coordinates. All the geodesics in this case which we deal with are boundary-to-boundary spacelike geodesics in global AdS 3 . A geodesic which starts from one side of the identification wedge (e.g. from the left on Fig.9B) and goes to W − comes out of the W + and has to go to the other side of the boundary, to the right in Fig.9B. Also a spacelike geodesic starting from the left cannot reach the W + surface first, before reaching the W − . An important point to note is that once a geodesic leaves the identification wedge W ± , it cannot enter it once more. The boundary-toboundary geodesics also cannot go through both V ± and W ± identifications, which is again evident from the black hole rest frame in the global coordinates. That means that winding geodesics cannot be also crossing, and vice versa. In the BTZ coordinates this is also clear from the fact that the part of winding geodesics which wraps around the horizon would have been lying inside of the cut out region close to the horizon, since the behavior given by formulas from the Appendix B dictates that winding geodesics always reach closer to the horizon than direct ones (more on that in the subsection 3.4), and the surfaces W ± can be considered as foliations of segments of direct geodesics in analogy to the considerations from the proof of the Proposition 1. These arguments imply the uniqueness of the crossing geodesic constructed from the image method as above, as well as the statement of the proposition. An important corollary is that the minimal spacelike geodesic connecting two points on the same side of the boundary is always the direct one. Now suppose that a and b are spacelike-separated points located on different sides of the boundary relative to the collision line. Then the crossing geodesic exists as long as image geodesics intersect the identification wedge 5 . The edge case in this situation is when o = o * belongs to the worldline of the particle. We have similar picture for direct geodesics in this case, which do not exist until they intersect the worldline. In the following subsection we begin to address this question in more detail in case of equal-time boundary endpoints.
We conclude this subsection by calculating the length of crossing geodesic in terms of BTZ coordinates of endpoints. As discussed above, it equals the length of a winding geodesic, either ab # or a * b. Suppose that endpoints a and b are parametrized as matrices according to (3.1). We write We again note that b # does not necessarily belong to the BTZ patch, however it does not matter since we do not have to calculate the actual coordinates of b # . To calculate the length, we take the trace: where u 1 is the holonomy of the particle given by (2.30). We introduce auxiliary notations Using the formula (2.30) for the holonomy and the matrix parametrization of the end- 5 Since the particle itself moves along the lightlike geodesic and image geodesics have to intersect W ± , image geodesics which constitute an existing crossing spacelike geodesic are always spacelike as well.
points (3.1), we come to the resulting expression: This expression has some symmetries: 1. The Z 2 -symmetry of the bulk background under reflection with respect to the collision line: 2. Replacing particle 1 with particle 2: 3. The symmetry between temporal and spatial coordinates of the center of the boundary segment on which the geodesic is anchored: t 0 ↔ ϕ 0 .
The first two symmetries are not surprising, but the last one is a somewhat unexpected unique feature of our bulk spacetime geometry. It is also worth noting that L cross is is a monotonically increasing function of t 0 and ϕ 0 .

ETEBA geodesics
The holographic entanglement entropy of a subsystem in the dual of a non-stationary bulk spacetime is calculated using the HRT prescription [4] as minimal surface anchored on the boundary region where the subsystem lives. In our case, this surface is a geodesic in BTZ coordinates anchored onto a segment [ϕ a , ϕ b ] on the boundary and with t a = t b = t 0 . In this subsection, we focus on such geodesics with equal-time boundary endpoints, which were labeled by Hubeny and Maxfield [12] as "equal-time-endpoint boundary anchored", or ETEBA geodesics. We will follow [12,13,19] and use this terminology. We are particularly interested in their behavior during time evolution. The HRT geodesic which computes the entanglement entropy is the minimal ETEBA geodesic which can connect a given pair of endpoints. For this reason, we also address the issue of existence of ETEBA geodesics of different types to know when they can and cannot participate in the HRT prescription.

Direct equal-time geodesics
These geodesics can be parametrized using equations (B.2,B.3,B.4). In the equal-time case, they have the form (we set λ 0 to zero): (3.28) From these equations it follows that these geodesics lie completely in the time slice t 0 , and their shape or length given by (3.10) does not depend on time. The question which we are interested in is at what times a direct geodesic exists between given ϕ a ∈ (−π, 0) and ϕ b ∈ (0, π). The direct geodesic does not exist when it crosses surfaces W ± . The worldlines of particles go into the bulk towards the horizon, and they are the closest points of the identification wedge to the boundary in any given time slice, see Fig.7. From the initial moment of time, the direct geodesic would have to intersect w ± and pass through the dead zone inside the wedge, therefore they do not exist for some time. However, as we evolve the system, the wedge gets smaller, and eventually it will get small enough to let the direct geodesic run clear of the identification fully inside the fundamental domain. From the shape of direct geodesics (which is constant under simultaneous evolution of both boundary points), shown e.g. in Fig.14B, it is clear that the worldlines of particles are the last points of the shrinking wedge which geodesics can touch before going completely out of the dead zone. The moment of emergence of a direct geodesic is illustrated in Fig.9A. The geodesic a 2 b 2 there touches the worldline. After this moment, the particles will move further into the bulk, the wedge will shrink more and the geodesic will lie completely in the fundamental domain of the identification, similarly to the smaller geodesic a 1 b 1 . Note that this argument holds both in case when we vary t 0 keeping ∆t = 0, and when we vary e.g. t a with fixed t b = 0. The first scenario is what relevant to the evolution of HEE, and the second scenario is relevant for time dependence of correlation functions. Now let us discuss the appearance of the direct geodesic in a bit more detail for the case of equal-time direct geodesics with ∆t = 0, since they are important for the calculation of holographic entanglement entropy. We can find explicitly the moment of time t cr = t a = t b = t 0 when the direct geodesic between points with given ϕ a and ϕ b , located according to the assumptions of the proposition, crosses the worldline of e.g. particle 1 with ϕ = 0. Suppose this happens at the timet in the point with radial coordinater. We use the parametrization of geodesics given by (3.25,3.26,3.27). The worldline equation is given by (2.39). Plugging it into (3.25) and solving in terms of λ, we find the value of the affine parameter at the intersection point as a function of t cr and ∆ϕ: Now we plug this into (3.26), requiring that ϕ(λ) = 0. This gives the equation, which we solve in terms of t 0 : cosh This is the moment of emergence of equal-time direct geodesic. Note that for symmetric segments, i. e. ϕ = 0, this expression reduces to The length of the direct equal-time geodesic is given by the expression (3.11) with t a = t b : This expression is time-independent.

Crossing geodesics with equal-time endpoints
The evolution of crossing ETEBA geodesics is the key to non-equilibrium dynamics of entanglement in our holographic bilocal quench setup. As explained in the previous discussion, the direct equal-time geodesics do not always exist. The following proposition establishes when crossing ETEBA geodesics exist and the chronology between vanishing of the crossing geodesic and emergence of the direct geodesic.
Proposition 3. For endpoints located as follows: ϕ a ∈ [−π, 0), ϕ b ∈ [0, π), and t a , t b chosen in such a way that points are spacelike-separated, it is true that: 1. For t a = t b = 0 the crossing geodesic always exists; 2. When increasing t a = t b = t, the crossing geodesic disappears only after the corresponding direct geodesic appears; 3. In the moment when the crossing geodesic disappears, L crossing (a, b) ≥ L direct (a, b). Proof. As mentioned earlier, we have to work in global coordinates to prove most of this proposition. To prove the point 1), we have to show that the image points are located in such a way that the image geodesics have to intersect the wedge at τ = − π 2 (which corresponds to t = 0 in BTZ coordinates. The global time slice τ = − π 2 is special because it is closed under the action of the isometries, that is (τ = − π 2 ) * ,# = − π 2 (see eqs. (A.9,A.12)). Therefore, we only have to prove that the angular coordinates of image points take desirable values. Specifically, if we have φ a ∈ (2π − θ, 2π) and where the angle θ is defined by sin θ = tanh πR . We use the formula (A.13) for the angle of the image point φ * a , setting in that formula χ → ∞ and τ = − π 2 : Fig.10. The coordinate of the image point is monotonically increasing function. We know that at τ = − π 2 there is a fixed point of the isometry located at φ = 0, which is where the particle sits, so φ # (0) = 0. On the other hand, by definition (3.13) we have θ ∈ W + ⇒ φ # (θ) = 2π − θ ∈ W − . By continuity and monotonicity, all image points with φ b ∈ [0, θ) will thus have angular . In other words, the actions of the #-isometry results in a rotation counter-clockwise with some angle less than 2π − 2θ. As a result, the φ # b will never end up in the interval (2π − θ, 2π). This is exactly what we needed to prove, and the argument for φ * a points is completely analogous. The point 2) concerns the time evolution of crossing ETEBA geodesics. From the transformation formulas from global to BTZ coordinates 2.6 we find that tanh Rt = coth χ cos τ cos φ ; (3.35) which means that on the boundary at χ → ∞ time evolution in BTZ coordinates corresponds to time evolution in global coordinates with angle-dependent rate. That means that once we fix the angles of endpoints in the initial moment, we can consider the evolution in global time to describe the evolution in BTZ time. Note, although, that in general case under the BTZ time evolution an equal-time geodesic on the initial time slice will be mapped to a geodesic with non-equal-time endpoints in global coordinates. The time coordinates of the endpoints will depend on their angular coordinates. However, in a special case of symmetric intervals ϕ a = −ϕ b the crossing ETEBA geodesic in BTZ coordinates will always remains an ETEBA geodesic in global coordinates. The statement 2) itself can be verified by plotting the geodesics in global coordinates and using the isometry formulas from Appendix A. The sample plots of geodesics are presented on the Fig.11. For symmetric endpoints, the direct geodesic appears precisely in the same moment when the crossing geodesic vanishes, as shown on Fig.11A. For general endpoints evolving in time, the direct geodesic appears before the crossing geodesic vanishes, as shown on Fig.11B. This argument can be strengthened by looking at isometry formulas and certain light cones. We know that the crossing geodesic vanishes after the image geodesics intersect the particle worldline in the same point o = o * . First, consider formulas (A.9,A.12) with χ → ∞: For fixed φ, it is clear that these are growing functions of t. Moreover, since the coefficient in front of the first term 1+2E 2 > 1, we observe that τ * ,# ≥ τ for τ ∈ [− π 2 , 0), where the inequality is saturated only in the initial moment. That means that generally Because of the continuity and monotonicity of geodesics in global time (see eq. (A.2)) that also means that  cone, to which the said geodesic belongs. The origin point of such light cone would be located on the boundary somewhere to the past of the geodesic. Since the * and # mappings are generated by parabolic Lorentz isometries, that means that both image geodesics a * b and ab # also belong to the same light cone. For a moment let us focus on the case of symmetric endpoints. In this case the particle lightlike worldline also belongs to this light cone in the moment of time when the given direct geodesic intersects it, which means that the image geodesic will also intersect it in the same moment, so we get precisely the picture shown on Fig.11A. For a case of general endpoints, the particle worldline intersects our imaginary light cone in the point of intersection of the direct geodesic and the worldline, when the direct geodesic appears, and then goes inside the light cone to the future, missing image geodesics. The time evolution with fixed angular coordinates of endpoints from that moment effectively means that we move the imaginary light cone upwards in time, while the worldline remains fixed. The intersection point between the imaginary light cone and the worldline will always move to the future, and inevitably will coincide with the point of intersection of two image A.
D. E. Figure 12. Evolution of ETEBA geodesics which compete in the HRT prescription during the particle collision process.
geodesics between themselves, eventually. This will be precisely the moment when the crossing geodesic vanishes. Thus the point 2) is proved. The point 3) can be seen from the geometric picture of geodesics, as illustrated in Fig.12D. The crossing geodesic at this point crosses the wedge in the position of the particle, and it consists of two AdS 3 geodesics joined together in that point in the bulk. Thus we have a curved triangle, all sides of which are spacelike geodesics in AdS 3 . From the Lorentzian AdS version of the triangle inequality for spacelike-separated points, it is therefore correct that.
which is precisely what we needed to prove. Note even though we work with geodesics which have diverging lengths we still can use the triangle inequality, since two corners of the triangle rest on the boundary, and we have identical divergences from both sides of the above inequality, understanding it in regularized sense.
The length of the ETEBA crossing geodesic is given by (3.25) with t a = t b = t: A.
B. C. Figure 13. Emergence of winding geodesics during the black hole creation. The blue geodesic has n = 1, and the green geodesic has n = −1. Here the disc depicts the bulk region r < 1, with the horizon at r = R = 0.3.
This is the primary formula we will use for study of non-equilibrium behavior of entanglement in our model.

Winding geodesics
The BTZ black hole solution admits infinitely many geodesics between two given spacelike-separated endpoints on the boundary. One of these geodesics is the direct geodesic and all other geodesics wind around the horizon. In BTZ coordinates all geodesics are parametrized by equations (B.2,B.3,B.4), direct and winding geodesics alike. For the direct geodesic |∆ϕ| = |ϕ b − ϕ a | < π, and for winding geodesics ∆ϕ > π.
Our holographic dual to the bilocal quench is the AdS 3 spacetime with colliding particles in BTZ coordinates. Locally, the geometry of this spacetime is identical to that of the BTZ black hole, so in principle the geodesic equations also admit winding geodesic solution. However, the topological identification can actually interrupt the existence of winding geodesics, just like in case of direct geodesics. So, in order to answer the question whether the winding geodesic exists, we have to check if it crosses surfaces W ± . A winding geodesic wraps around the horizon. It approaches to the horizon into radial distance given by eq. (B.11): where n = 0 ∈ Z is the winding number. From this formula it is clear that between two given endpoints all winding geodesics approach to the horizon closer than direct geodesics. Moreover, there is a strict hierarchy between the depths of windings: a winding with higher n always approaches closer to the horizon than any winding with lower n. Also, for smaller ∆ϕ winding geodesics approach the horizon closer.
The identification wedge shrinks with time according to equation (2.38) (see Fig.6). Thus we arrive to the following statement: Proposition 4. For given angular coordinates of endpoints ϕ a , ϕ b and the winding number n the corresponding winding geodesic exists only for t 0 = 1 2 (t a + t b ) >t n , wherẽ t 0 is the moment of time when the winding geodesic crosses synchronous slices w ± of the wedge faces at particle worldlines.
In other words, winding geodesics with given time separation between the endpoints will appear one by one as particle go deeper into the bulk towards the horizon. Geodesics with higher winding numbers will appear later than those with lower winding numbers. We illustrate the emergence of equal-time n = ±1 winding geodesics on Fig.13. On Fig.13A the wedge is still too large to accommodate both winding geodesics. Note that the direct geodesic at that time between two points on the lower half already exists from the beginning, according to proposition 1. On Fig.13B the identification wedge decreased in size enough to allow the blue winding geodesic to appear, but the green geodesic still intersects it. On Fig.13C even later moment of time is shown, when the green geodesic appears. In the initial moment t = 0 there are no winding geodesics, and in the limit t → ∞ the infinite amount of winding geodesic appears, which is the same situation as in the BTZ black hole spacetime. Note that these time scales of appearance of winding geodesics are generally much larger than the time scales of appearance of direct geodesics and vanishing of crossing geodesics. Holographically this means that winding geodesics are irrelevant for equilibration of entanglement, however they do contribute to late time behavior of correlation functions, as we will discuss later.
To conclude this subsection, we remind that length of a winding geodesic is given by the formula (B.15): The key observation is that lengths of winding geodesics between the given endpoints on the boundary are always larger than the lengths of corresponding direct and crossing geodesics.

Equilibration of entanglement
We now have at our disposal all tools needed to perform the holographic computation of the entanglement entropy in the boundary CFT after the bilocal quench and analyze its time dependence during the thermalization process. We start our analysis with calculation of the holographic entanglement entropy (HEE). We then use it to investigate the equilibration and spreading of entanglement in subsystems in the boundary theory which live on segments of the circle during the particle collision process. We also compute holographic mutual information and discuss different possibilities of its non-equilibrium behavior, depending on the location of subsystems.

Holographic entanglement entropy
Consider a subsystem in the boundary theory which lives on a segment L of the circle, which is bounded by points a and b. Then the entropy is calculated, according to the Ryu-Takayanagi proposal [3] generalized to the non-stationary case [4], as the minimal area of the codimension two surface in the bulk anchored on equal-time points a and b.
In an asymptotically AdS 3 spacetime, this surface is a geodesic, so we need to find the minimal geodesic between two given equal-time points on the boundary and calculate its length: Here G = 3 2c is the gravitational constant, and c is the central charge of the boundary CFT. In our case, the bulk background is set up explicitly as a quotient of the AdS 3 spacetime by a non-trivial identification. That means that the metric itself is stationary in our case, but the identification is non-stationary, which is what makes the entire spacetime non-stationary and requires the use of the HRT proposal, which generalizes usual Ryu-Takayanagi prescription to the non-stationary case. Our bulk spacetime is arranged in such a way, that depending on the position of the endpoints, the crossing geodesic either participates in the HRT prescription, or does not. These two situations describe qualitatively different behavior of entanglement.
In the BTZ black hole spacetime, which corresponds to the CFT at thermal equilibrium, the minimal geodesic is a direct geodesic. The length of such geodesic gives a result for HEE, and is obtained from (3.10) for small subsystems with size less than half of the circle, ∆ϕ ≤ π, by setting t a = t b . Introducing the UV cut-off in the boundary theory we have: This is the entanglement entropy of the thermal equilibrium state with temperature given by T = R 2π of a subsystem of size ∆ϕ = ϕ b − ϕ a < π. For large subsystems with ∆ϕ > π, one obtains the HEE from the expression (3.11) instead, which results in Since the geodesics which do not intersect W ± are identical to those in the BTZ black hole spacetime, direct equal-time geodesics always govern the equilibrium regime in our model. Meanwhile, the length of a crossing ETEBA geodesic (3.41) is time dependent because of the shape of the identification wedge, hence crossing geodesics must govern the non-equilibrium regime. In the further discussion we will focus on the small subsystems with ∆ϕ ≤ π. The results for small subsystems can be related to large subsystems if one keeps in mind the subtraction of 2π in the equilibrium HEE formula (4.4) and symmetries of the formula (3.41). Since the horizon never appears in any finite time, we do not have to worry about its contribution in the Ryu-Takayanagi prescription. Note also that the HEE of the complement of a subsystem is always given by the same geodesic as the HEE of the subsystem itself, because of the same reason. This is what we expect when considering the evolution of a pure state in the boundary CFT.

Equilibrium in the initial state
From the proposition 1 it follows that one can anchor a direct geodesic on segments of the boundary spatial circle which lie to the side of the collision line at any moment of time. Therefore the entanglement entropy of a subsystem located in either upper or lower (with respect to the collision line) semi-circle is maximal and is given by the expression (4.3), and is constant in time. Therefore we can come a conclusion that for subsystems which lie in between the initial excitations (orange subsystems on Fig.14) the single-interval HEE is exhibits constant equilibrium behavior. To study non-equilibrium dynamics in these subsytems, one might want to use more "fine-grained" observables, which would require information from the bulk beyond the minimal geodesic length. For example, one can consider Renyi entropies [23][24][25], entwinement [58], or contributions to holographic correlation functions from the non-minimal geodesics [37,46,47]. We will discuss the latter in the bilocal quench scenario in section 5.

Crossing geodesics and non-equilibrium regime
On the other hand, the subsystems with endpoints on different sides of the boundary with respect to the collision line are initially out of equilibrium. This is because in the initial moment the corresponding HRT geodesic is a crossing geodesic, and the evolution of crossing and direct HRT geodesics for such subsystems are described by the proposition 3. These subsystems contain one of the excitations on the boundary in the initial moment, like illustrated by blue subsystems in Fig.14. Let us restrict our attention without loss of generality to subsystems which contain the particle 1 moving along the ϕ = 0 worldline. Then we are interested in the case when ϕ a ∈ (−π, 0), and ϕ b ∈ (0, π). The evolution of the entanglement in such subsystems goes as follows, according to the proposition 3: • From the point 1) of the proposition, at t = 0, the HRT geodesic has to be a crossing ETEBA geodesic, as shown on , which we call the thermalization time of the subsystem. The points 2) and 3) of the proposition 3 ensure that the crossing geodesic still exists in that moment and the transition is continuous, as we will see below.
• At late times, for t > t , the behavior of HEE is governed by the direct equaltime geodesic, which corresponds to the equilibrium regime and is expressed by the formula (4.3). The point 3) of the proposition 3 ensures that the vanishing of the crossing geodesic does not influence the behavior of HEE (see Fig.12D-E).
Thus, the general formula for HEE of a crossing subsystem can be expressed as S(a, b|t) = min S non-eq (a, b|t) S eq (a, b) ; (4.5) where S non-eq is the contribution to HEE from the crossing ETEBA geodesic. It is obtained using the formula (3.41) with the boundary UV-cutoff introduced as in (4.2) for the length of the crossing ETEBA geodesic: where we remind that E = coth πR 2 , ∆ϕ = ϕ b − ϕ a and ϕ 0 = 1 2 (ϕ a + ϕ b ). At E = 0 one can recover the equilibrium result (4.3). The initial value at t = 0 is given by (4.7) As mentioned above, this quantity is smaller than the equilibrium value (4.3) for the same segment. From the formula (4.6) it is clear that the quantity S non-eq (a, b|t) grows monotonically with time, until the crossing geodesic vanishes. However, as we evolve the system with time, after the moment t = t cr the direct geodesic appears as well and starts competing in the HRT prescription, and takes over at t = t (a,b) * , realizing the saturation of HEE at equilibrium. Now let us discuss the time dependence of the HEE in more detail. We illustrate the typical time dependence of ∆S(t) = S non-eq (t) − S eq on Fig.15.
In the further discussion we assume that we deal with the formation of large black holes 6 with R > 1, or in other words when the temperature is higher than the AdS 3 Hawking-Page temperature, T > 1 2π . The formulae (4.6) and (4.3) are still valid for R < 1, but the bulk geometry with a small black hole gives a subdominant contribution to the CFT path integral and thus is not a proper holographic dual for a thermal state [60].
Early-time evolution. We can expand the expression (4.6) in time around t = 0. The leading terms of the expansion read: S non-eq (a, b|t) = S non-eq (a, b|0)+ c 6 (4.8) 6 Nevertheless, some illustrations are still presented with R < 1. This is done to make certain features more prominent on the availible scale.  where S non-eq (a, b|0) is given by the formula (4.7). This behavior is the same as the socalled pre-local equilibration growth [8,9], which appears in the Vaidya global quench scenario. On the Fig.15 it is shown that the HEE behavior is well approximated by the quadratic expansion until a time scale t = t 1 . We observe that We conjecture that this time scale has similar meaning as the "local equilibration" time scale in the Vaidya quench. However, the important difference of our case is that we do not make any assumptions about the size of the subsystem compared to the horizon radius, and due to this the time scale depends on ∆ϕ, unlike in AdS-Vaidya case in Poincare coordinates [9]. Intermediate regime and linear growth. For t > t 1 the behavior of HEE deviates from the quadratic growth. It starts approaching the linear regime, and reaches the linear asymptotic growth at a time scale t = t 2 . The time scale t 2 depends inversely on the horizon radius R, which is shown on the plot 16A, and thus the higher the temperature, the more prominent is the linear growth regime. The leading linear asymptotic behavior can be established by expanding the difference S non-eq − S eq in e −Rt . The result is the following 7 : The global quench models also exhibit the linear growth regime, which is evident from both CFT calculations [20,21] and holographic calculations [5, 6, 8-13, 16, 17, 19]. In the context of holographic thermalization global quench scenarios the linear growth regime is often referred to as entanglement tsunami [8][9][10]. It was found [8][9][10][11] that the linear behavior in global quench models is universal and can be expressed as Here s eq is the equilibrium density of HEE, A is the surface area of the boundary of the subsystem, and v E is the entanglement velocity, which depends only on the dimension of the spacetime. We can make contact with our case in a similar way to the discussion in [13], if we consider the case R ∆ϕ/2 1. In this case the equilibrium entanglement entropy (4.3) is given by pure area law (we omit the UV cutoff): from where we find that s eq = c 6 R. Taking into account that A = 2, since our subsystem is bounded only by two points, we find that the asymptotic linear behavior (4.11) should look like The universal result for global quenches in 2d CFT is v E = 1, and from comparing of (4.10) to (4.13) we see that this value for the entanglement velocity holds true for our bilocal quench scenario as well. Thus we obtain another argument for the fact that the notion of the entanglement velocity and its bounds is relevant not only for global quenches, but also for local quenches as well [29,30,32]. As it turns out, the linear growth regime in our case is directly related to the memory loss regime [8,9,19], since the function ∆S in the leading order only depends on the difference U = t = ∆ϕ 2 . Another interesting fact that the time scale t 2 is related to the crossing HRT geodesic going inside the horizon, which is an evidence for the fact that the linear growth of HEE is related to the HRT geodesics probing the interiors of the black hole. We discuss these observations in more detail later in subsection 4.1.4.  Thermalization. The transition to the equilibrium regime happens when the crossing geodesic and the direct geodesic have the same length. Hence thermalization time can be found from the condition (4.14) where we have emphasized the time dependence in the length of the crossing geodesic. Using the formula (3.41) for l. h. s. and the formula (3.32) for r. h. s., we find the expression for thermalization time: It is important that this moment in time is later than the time of emergence of the direct geodesic t cr given by (3.30), but it is also before the time when the crossing geodesic vanishes, since the point 3) of the proposition (3) directly states that the crossing geodesic vanishes at some time when its length has grown higher than the length of the direct geodesic. Thus we have the continuous transition from the non-equilibrium growth to saturation of HEE happening at t = t (a,b) * . The formula (4.15) dictates that larger subsystems thermalize slower, see plot 16B, which is expected. Also let us note that for symmetric intervals ϕ 0 = 0 one can obtain from (4.15) and (3.30) the result This is the same thermalization time as in the AdS-Vaidya quench model [5,6,8,13]. For a subsystem of the same given size ∆ϕ but non-zero ϕ 0 the thermalization time defined from (4.15) will be smaller than ∆ϕ 2 , as shown on the plot 15C. Now let us discuss the character of the transition to saturation. While the HEE itself is continuous, which is what expected from the thermalization models, particularly those based on the global quench scenario [5][6][7][8][9][10][11][12][13][15][16][17][18][19][20][21], the time derivative of the entanglement entropy is discontinuous at the transition point, which results in sharp transition to saturation. This fact is the key difference of non-equilibrium dynamics in our model from dynamics in holographic Vaidya [5-11, 13, 15-19] and end-of-world brane [11,17] quantum quench models in 2d CFT. However, the sharp transition to saturation, when the linear growth regime lasts all the way to thermalization, is remarkably similar to that of the quasiparticle picture of entanglement spreading in 2d CFT [21] and to equilibration of a strip subsystem after the global Vaidya quench in d ≥ 3 [8,17].

Emergent light cone
On the Fig.17A we plot ∆S as a function of time and ϕ = ∆ϕ/2 in the special case of centered subsystems, ϕ 0 = 0. The expression (4.6) in this case simplifies to S sym non-eq (ϕ|t) = From this picture and for the formula for the thermalization time t = ϕ (4.16) we can see that the entanglement spreads along an effective light cone. The sharp saturation, makes the light cone prominent, again hinting at similarity with the quasiparticle picture of entanglement spreading [21]. The effective light cone velocity is related to the butterfly velocity v B , which is the speed of propagation of quantum chaos in thermal state [16,[49][50][51]. In the setting of global quench in two-dimensional holographic CFT it is true that v LC = v B = v E = 1 [16]. From the plot 17A we observe that v LC = v E = 1 also holds in our case of the bilocal quench.
One can also verify that v LC = v B = 1 in our quench scneario using considerations analogous to those in section 5 of [16], which relate thermalization of perturbations with the concept of entanglement wedge reconstruction. The spreading of information on the boundary can be characterized by the rate of growth of a boundary region which contains the information about an initial excitation (for concreteness, we consider the particle 1). In analogy with [16], we assume that the full information about the infalling excitation is contained in the centered boundary segment for which the direct geodesic crosses the wordline of the particle in the given moment of time. The bulk subregion bounded by this segment and the direct geodesic which crosses the worldline of the particle is the entanglement wedge [52][53][54][55]. The corresponding boundary segment encodes all information about local physics in this bulk subregion [55]. In this case the rate of growth of this boundary segment is given by the formula for the thermalization time (4.16). From this formula it follows that v LC = 1.
Another point to note is that the dependence of S non-eq (0) as a function of ϕ 0 in (4.7) and the time dependence in (4.17) are identical. Moreover, if one considers the full entropy at t = 0 as a function of ϕ 0 , then the dependence will be the same as the time dependence of the entropy for ϕ 0 = 0. In the first case the saturation which happens for |ϕ 0 | > ∆ϕ 2 corresponds to the case when the segment completely lies in the upper or lower half of the boundary spatial circle and doesn't contain an initial excitation -as we discussed, all such subsystems are at equilibrium from the beginning. This symmetry between temporal and spatial coordinates of the segment appeared a symmetry 3) in the formula for the length of the crossing geodesic (3.25). In the case of HEE dynamics, this symmetry can be qualitatively explained using the quasiparticle model [21]. In 2d CFT after the quench the spreading entanglement can be approximately modeled as propagation of EPR pairss of particles, which are created in the moment t = 0 by the quench and propagate with the speed of light [11,14,21]. In our case these EPR pairs are created in two separate points: ϕ = 0 and ϕ = π. A subsystem with a given size ∆ϕ = 2ϕ will equilibrate once the quasiparticles will reach the its boundaries. If the subsystem is centered on the source of particles, i. e. ϕ = 0, then a given moment of time t the quasiparticles will reach the coordinate ±φ =t. On the other hand, this is the same situation as if particles were created at the point ±φ at t = 0. When t = ϕ, the particles reach the boundary of the segment, or, equivalently, the source of particles is at the boundary, thus equilibration happens.
This quasiparticle explanation hints that this symmetry is caused by the production of entanglement by a point source from a bilocal quench. Consequently, there is no analogue of such symmetry in any global quench model because there is always full translational invariance to work with, and EPR pairs are produced uniformly in every point of space. It is also unlikely that this symmetry would hold in a higher-dimensional generalization, because in general velocities which characterize the spread of entanglement are slower in higher dimensions [8,9,16,49]. However, we can expect that a CFT calculation of correlation functions after bilocal quench in 2d CFT using the techniques along the lines of [15,24] can help reveal its nature in relation to some deep symmetries in CFT such as modular invariance.

Entanglement tsunami and memory loss
As we mentioned in subsection 2.2.2, it is known that a closed system loses information about the details of the initial state upon thermalization. In the previos work [8-10, 17, 19] it was found that the memory loss phenomenon in the context of holographic global quenches is prominent in the late-time evolution of HEE. Namely, at late times for large enough subsystems it was found out that HEE in the leading order only depends on the difference of time and the subsystem size, and separate dependence on those variables are exponentially suppressed at late times. This is referred to as memory loss regime of HEE, and it happens shortly (in relative terms) before saturation in the global quench scenario. In this subsection we argue that the model of holographic bilocal quench also exhibits a property of memory loss for large subsystems, and this memory loss seems to be related with the black hole interior.
To observe it, we consider the expression for ∆S(t) = S non-eq (t) − S eq , where S non-eq is given by (4.6) and S eq is given by (4.3) in terms of coordinates of boundary points ϕ a , ϕ b instead of ∆ϕ and ϕ 0 . We find asymptotic expansion of ∆S, taking e Rϕa ∼ e −Rϕ b 1 and also considering the limit of late times, such that e −Rt → 0, but e R(t−ϕ b ) and e R(t+ϕa) are kept fixed. The asymptotic has the form So we see that indeed the late-time dynamics of entanglement of large systems are characterized by the functional dependence of time and endpoints of the subsystem as light cone-like variable combinations t−ϕ b , t+ϕ a . This asymptotic is dominant for any values of the horizon radius. For large horizon radius R > 1, one can expand further to obtain the following leading behavior: This expression is the linear growth regime discussed above (4.10), but with the additional requirement of large subsystem. Thus we conclude that for large subsystems during relaxation to equilibrium at high enough temperature the memory loss regime is basically the linear growth regime. On the Fig.17B we  lines become horizontal is where the distinguishable dependence on ϕ gets suppressed, which illustrates the memory loss regime. In other words, we observe that for large subsystems with boundaries far away from initial excitation at ϕ = 0 the entanglement propagates as a wave front with retarded coordinate t − ϕ, and for high temperatures this wave-like behavior can be identified as the entanglement tsunami linear growth regime. This picture is largely similar to the global quench story [8,9] in 2d.
We also can take into account that initially we have two excitations, and segments of similar size and with positions that mirror each other on the boundary circle will equilibrate synchronously. That means that on upper or lower half of the boundary circle we will have two waves of entanglement propagating in opposite directions. One can therefore make a speculation that the instant equilibrium of HEE in segments which belong completely to the upper or lower half of the circle can be caused by something similar to a standing wave of entanglement, which emerges as a result of clash of two waves of entanglement from initial excitations.

Linear growth and black hole interior
The peculiar feature of non-local observables such as HEE in holographic non-equilibrium processes is that they can probe the region inside the apparent horizon. More specifically, in the holographic global quench thermalization scenario it was understood [10,11] that HRT surfaces which calculate HEE in the boundary can probe the interior of the black hole, i. e. region of the spacetime inside the horizon. Moreover, it was found that the linear growth of HEE comes from growth of the piece of the HRT surface which lies inside the horizon. In the present subsection we discuss similar feature of the holographic bilocal quench model. We start by noting the fact that a crossing HRT geodesic which consists of two pieces of image geodesics, can reach inside the black hole. This happens because, as we discussed in section 3, the image geodesics can span the entire global AdS 3 cylinder, and not just the fundamental domain of the topological identification. Because of this and the shape of the identification wedge relative to the image geodesics (see Fig.9), the image geodesics can probe the interior of the black hole, which is located behind the horizon (see Fig.3). More important is the fact that the actual crossing geodesic, that is pieces ao and o * b of the image geodesics, can reach inside the horizon as well. In other words, a crossing geodesic probes the interior of the black hole when points o and o * lie inside the horizon. It appears that in our case the piece of the crossing geodesic lying inside the horizon is responsible for the linear growth of HEE. We were not able to prove this statement analytically like it was done in global quench scenario in [10,11], but we can establish this observation based on numerics. We also can obtain a bound for the subsystem size which separates subsystems which probe the black hole interior and thus exhibit the linear growth of HEE, from those which do not, in a special case.
To accomplish this, let us consider symmetric segments only, with ϕ b = −ϕ a = ϕ. For these segments, as discussed in the proof of the point 3) of proposition 3, the direct geodesic emerges in the same moment when the crossing geodesic vanishes, that is when o = o * . That means that we can use the condition o = o * as a condition of thermalization in global coordinates for symmetric segments. Because intervals are symmetric, it is true that L(a, o) = L(o * , b). But it is also true that L(o * , b) = L(o, b # ). That means that when o = o * , the point o bisects the image geodesics in half. In particular, As also discussed in the proof of the proposition 3, the functions τ b # (τ b ) and, conssequently, τ o (τ b ) are monotonically increasing faster than the value of τ b = τ a itself. That means that we can detect the enetering of the crossing geodesic inside the horizon by looking first at when the point o = o * will enter the horizon in global coordinates. The point o = o * is also located on the worldline of the particle 1, which enters the horizon at τ = − π 4 , which can be found from the worldline equation (2.26) and the horizon surface equation (2.22) [34]. Thus we come to a condition This is the moment of global time when a subsystem will thermalize exactly at the time when the deepest point of the HRT geodesic o = o * will just reach the horizon. The smaller subsystems which thermalize faster (because of the relation t = ϕ) will not probe the interior, and larger subsystems which thermalize longer will probe the interior. Keeping in mind that τ a , τ b # ∈ [− π 2 , π 2 ], we can rewrite this equation as where we used the formula (A.12) for τ # at the boundary χ b → ∞. Next, we use the equalities φ b = −φ a , τ b = τ a and transform the equation (4.21) into BTZ coordinates using the formulas (2.6: Then we substitute the expression (4.16) for thermalization time for symmetric subsystems t = ϕ, the definition of E = coth πR 2 , and then solve for ϕ. The result is (4.23) Figure 19. The grey region is the entanglement shadow, which appears when t > π 2 . The blue curve is the maximal RT geodesic. This is the lower bound for the size of the symmetric segment which will probe interior of a black hole of the given size at some point during its equilibration. We denote the corresponding thermalization time as t hor We plot time evolution of HEE for symmetric segments which do or do not probe the interior on the Fig.18. The blue curve corresponds to a segment that doesn't probe the interior since its size is smaller than ϕ hor . The green curve corresponds to the critical segment with ϕ = ϕ hor , which just touches the horizon in the moment of saturation. Red and brown curves all correspond to large segments that probe the interior and, as seen on the plot, exhibit the linear growth regime. Larger segments have longer linear growth. The time scale t 2 which signifies the start of the linear growth (see Fig.15) is the time scale when the crossing HRT geodesics penetrates the horizon. The formula (4.23) gives the lower bound for t 2 in case of symmetric segments for given R.
Thus we see that the probing of the black hole interior by crossing HRT geodesics is related to the linear growth of HEE. In the previous subsection we also discussed that for large temperatures the linear growth is related to the memory loss regime for a given subsystem. We see that in the bilocal quench model it is evident that loss of memory of the initial state at thermalization for a given subsystem is related to the problem of probing the black hole interior.

Scrambling time and entanglement shadow
Bigger subsystems thermalize longer, according to the formula (4.15) (see Fig.16B). On the other hand, among the segments of the same size ∆ϕ the segments which thermalize longest are the symmetric segments with ϕ 0 = 0. From these considerations it follows that there are two distinct subsystems among all subsystems which thermalize last: it is the symmetric segment with ϕ b = −ϕ a = π 2 and its complement. The thermalization time of these subsystems can be calculated from the formula (4.16) and equals to (4.24) At this time, all subsystems, which include half of total degrees of freedom in the boundary theory or less, are thermalized. This timescale is often referred to as scrambling time [56,57] in the context of relaxation of small perturbation of the thermal state. This is the time scale when the information about the initial state is scrambled so thoroughly that it cannot be restored from any small fraction of the total amount of degrees of freedom. Now let us ask the question: how the fact that after the scrambling time all small subsystems are thermalized is reflected in the bulk geometry? All small subsystems being at equilibrium means that their entnaglement entropy is governed by direct geodesics, which means that colliding particles are not probed by HRT geodesics anymore. In other words, beginning from t = t * , particles are located in a region of the spacetime which is not probed by the entanglement. This signifies that the entanglement shadow [58,59] appeared in the bulk in the moment t = t * , that is the region which is not probed by RT surfaces. This is a region between the horizon and the radius given by the depth of the direct geodesic with opening angle π or, equivalently, from the position of the particle at t = π 2 , see Fig.19. Using the latter definition and the worldline equation (2.39), we get Thus we get that the dimensionless radius of the entanglement shadow equals the energy of the particle 1 in global coordinates. This entanglement shadow is precisely the same as the entanglement shadow of the BTZ black hole spacetime [58,60]. Thus we observe that during the evolution of a pure state after the quench the information about the initial state gets scrambled when the ingoing matter falls into the entanglement shadow, which is associated with the forming black hole. While the true mixed thermal state is unreachable during the unitary time evolution, we do get the state, where all the initial data is completely scrambled over the entire system.  Holographically this means that the horizon does not form in finite amount of time, however the entanglement shadow appears after the scrambling time, which is finite for compact CFT.

Mutual information
In this section we consider the behavior of mutual information using HRT prescription. The mutual information for two disjoint regions is defined as: where S(A ∪ B) is a joint entropy, i.e. a minimal geodesic length between two possible types of geodesics. It means that the mutual information is zero if two regions are well separated. We compute the mutual information for different scales and locations of regions and analyze the time dependence. We are interested in equilibration of mutual information that can be described as reaching plateau behavior.
The mutual information in thermal equilibrium is either zero or a positive constant. We expect that the mutual information reaches equilibrium regime when all terms in the eq.(4.26) saturate. We present the most destinctive cases of mutual information behavior as function of time in the Fig.20. We can see that the equilibration comes at the different moments of time that we called t AB * that has been considered in several papers [13,18,23,27,28]. A value of this time scale depends on the location and size of the regions, but is bounded from above by the scrambling time of the entire system t * = π 2 (see the discussion above).
We have plotted the pictures of the mutual information in dependence of time up to scrambling time (4.24). The mutual information behavior in the system with two excitations is generally similar to the behavior in AdS-Vaidya geometry [13], however there are kinks that can be described differently which occur when the change of a regime happens in one of the terms in (4.26). In the Fig.20A. two kinks correspond to the thermalization of the small and large joined geodesics. In the Fig.20B. the first kink means that a small subsystem thermalizes as well as the second kink corresponds to the thermalization of the second subsystem. The mutual information in the last case, presented in the Fig.20C., has a peak that corresponds to the thermalization of the subsystem which contains the initial excitation.

Thermalizing correlators
Another way to probe thermalization is to study two-point correlation functions in the state in the boundary CFT which is dual to the bulk spacetime with a forming black hole. In this section we discuss non-equilibrium two-point correlation functions of a scalar operator O ∆ . We calculate the Lorentzian two-point correlation functions in the framework of the geodesic approximation [35], which holds for ∆ 1. We apply the geodesic approximation in the BTZ coordinates of AdS 3 spacetime with colliding massless particles. According to this prescription, to calculate the two-point correlator between points a and b on the boundary of a given asymptotically AdS spacetime, one has to sum over all geodesics between these points. The original formulation [35] is valid in the Lorentzian signature only when the points on the boundary are spacelikeseparated. In the case of timelike-separated points, the situation is more tricky, and a generalization is needed, e.g. see [7,47]. Because of this issue, we consider two cases separately. The pole structure of Lorentzian correlators is recovered by introducing the phase factors which result in a desired i -prescription. With that in mind, the schematic formula for the correlator between spacelike-separated points a and b in geodesic approximation reads: where index A stands for ret (retarded), F (Feynman) or W (Wightman) correlators. The lengths of the geodesics are appropriately renormalized by subtraction of diverging part. The limits of summation are defined by the concrete geometry of the defect in the AdS 3 spacetime. It can include a finite number of winding geodesics, like in case of a point particle [37,[45][46][47][48]61], or an infintite amount of winding geodesics, like in case of the BTZ black hole [62], see formula (C.1). In the general case, the set of values of the summation variable n depends on the choice of points a and b as well. However, one can always distinguish the dominating contribution coming from a minimal geodesic.

Leading contributions
In the section 3 we have established that if points a and b are spacelike-separated, one can always distinguish a minimal geodesic among the set of all geodesics between them -this is either direct geodesic, or the crossing geodesic, depending on the location of the endpoints. So for the most of this section we will only consider the leading contributions to correlators: Here Φ A is an appropriate factor for a given correlator dictated by the i -prescription in the Lorentzian signature. The explicit formulas for different Φ A are given by (C.3,C.4,C.5). Let us restrict ourselves to Wightman and Feynman correlators. Then for spacelikeseparated points Φ F,W = 1, and we omit these factors in this subsection. Analogously to HEE, we have two kinds of behavior of correlation functions depending on the location of the endpoints.
1. Suppose that endpoints are located to the same side of the collision line, that is either ϕ a , ϕ b ∈ [0, π] or ϕ a , ϕ b ∈ [−π, 0]. Then the minimal geodesic is the direct geodesic, as stated by proposition 1. The length of the direct geodesic is given by (3.10). We renormalize it by subtracting the 2 log r 0 R piece, and thus the correlator is thus given by This is the same result as in case of thermal equilibrium (C.1), and it is expected since in section 4 we discussed that with these endpoints HEE, which is also defined by the minimal geodesic length, is at its equilibrium value. 2. Suppose that endpoints are located on different sides of the boundary with respect to the collision line. More specifically, suppose that ϕ a ∈ [−π, 0] and ϕ b ∈ [0, π].
Then the minimal geodesic is either the direct geodesic, or the crossing geodesic. In the general case the lengths of these two geodesics are comparable, as evident from our discussion of HEE, so we take into account both these contributions on equal footing. However, as we discussed, generally crossing and direct geodesics for such choice of endpoints do not always exist, so we need to account for that as well. For this purpose we introduce auxiliary Θ-functions: The idea is that Θ = 1 if the corresponding geodesic exists, and Θ = 0 if it does not. The crossing geodesic exists when the image geodesics cross the identification wedge, and the direct geodesic exists when it does not cross the identification wedge and fully belongs to the fundamental domain in the bulk. Hence the definition. Now, using the formula for length of the crossing geodesic (3.25) and for length of the direct geodesic (3.10) and again subtracting the divergence, we can write the expression for the correlator as: This is a discontinuous function because of the fact that different geodesics not always exist. We will discuss the issue of these discontinuities in geodesic approximation in more detail in the next subsection. For now, let us consider the correlation function (5.6) in more detail in particular case of equal-time points t a = t b = t. In this case we can write Θ-functions as Heaviside step functions of time: The first equality reflects the fact that the direct geodesic exists after t = t cr given by formula (3.30). The second line means that the crossing geodesic exists until the time when the corresponding image geodesics cross W ± in the same point o = o * located on the worldline of the particle 1. The equal-time correlator then reads In the section 4.1.2 we have calculated the thermalization time t (a,b) * of a segment between a and b (4.15). From the perspective of correlation functions, this is the time when two terms in (5.9) exchange dominance: for t < t (a,b) * the first term is leading and we have non-equilibrium regime, and for t > t (a,b) * the second, equilibrium term dominates. Note that from proposition 3 it follows that t cr < t (a,b) * < t o=o * . That means that in general there are two small intervals in the end of non-equilibrium regime and in the beginning of the equilibrium regime when both terms contribute. This is a feature which arises in correlation functions, but the HEE is always defined by a single minimal geodesic and therefore is insensitive to such details.

Contribution of windings
In the section 3.4 we have discussed that in the bulk between the two given boundary points winding geodesics emerge at late times. These winding geodesics give subleading contributions to two-point functions in geodesic approximation, according to the general formula (5.1). The correlators dual to BTZ black hole spacetime include the contribution of infinite amount of winding geodesics (C.1), whereas in case of our bulk spacetime with particles the identification wedge prohibits the existence of any winding geodesics at least at early times. Therefore, it is interesting to look at the emergence of winding geodesics as subleading contributions to the correlation functions, as they tell us how thermalizing correlation functions approach the equilibrium correlation functions (C.1) at late times.
In the proposition (4) we have discussed that a winding geodesic between points a and b with the winding number n exists after the timet n , when it is able to wind around the horizon without intersecting W ± , see Fig.(13). We thus can define the Θ-functions which regulate the presence of winding contributions in the correlation functions: Θ winding (a, b; n) = 1 , if the n-th winding geodesic does not intersect W ± 0 , if the n-th winding geodesic intersects W ± . ; (5.9) where n = 0 is an integer. Note that Θ winding (a, b; |n|) = Θ winding (a, b; −|n|) (as can be seen from Fig.13). According to the proposition 4, we can simply write Θ winding (a, b; n) = θ(t 0 −t n ) ; (5.10) and also in the limit t 0 → ∞ it is true that Θ winding (a, b; n) = 1 for any n and any angular coordinates of endpoints a and b. Using these functions and the expression for the length of a winding geodesic (3.42), we conclude that the general form of the twopoint correlation functions with spacelike-separated points in geodesic approximation is the following, depending on the location of the endpoints.
• For ϕ a , ϕ b ∈ [0, π] or ϕ a , ϕ b ∈ [−π, 0], we have From these formulas (5.11-5.12) it is clear that in the limit t → ∞ one recovers the equilibrium two-point function (C.1). The geodesic approximation dictates that the subleading terms in the image sum appear one by one with time, as the particles move towards each other in the bulk. The presence of Θ-functions in the correlators obtained from geodesic approximation when there are multiple geodesics possible between two given endpoints are a common occurence in cases of locally AdS 3 spacetimes with singularities which admit multiple geodesics. A prime example for such occurence is the AdS 3 with a point particle [35,[45][46][47]61]. Our bulk spacetime is a special case of AdS 3 with two massless point particles. Correlators dual to AdS 3 with two particles were studied in global coordinates in [37,48], and they also exhibit such discontinuities. On the example of massive static particle in AdS 3 it was shown in [61] that geodesic approximation gives a continuous result which coincides with the correlators obtained from GKPW dictionary, if the spacetime is an AdS 3 orbifold. In the general case, however, the geodesic correlators are discontinuous, and these discontinuities are smoothened by corrections to the geodesic approximation, which is evident if one calculates correlation functions using e.g. full GKPW dictionary [61]. In the orbifold case, these corrections to the geodesic approximation are exactly zero.
In our case, the bulk spacetime itself is not an orbifold. However, as the time goes on, the spacetime geometry approaches that of the BTZ black hole spacetime, which is an oribfold. So one could say that in our case the bulk spacetime looks more and more like an AdS 3 orbifold. So if we assume that considerations from [61] can be extended for our case of the AdS 3 spacetime with two massless particles, then it can be expected that the corrections to the geodesic approximation will be diminishing at late times.

Timelike correlations
Now let us turn to the discussion of the time dependence of two-point correlation functions. We consider the case of endpoints a = (0, ϕ a ) and b = (t, ϕ b ) (with t > 0) and study the time dependence of correlation functions in the approximation of the leading term in geodesic approximation (5.2). First, let us note that when t is small enough such that t 2 − ∆ϕ 2 < 0, then the interval between the endpoints is spacelike and we can just use the geodesic approximation as discussed above. However, for timelike-separated points on the boundary t 2 − ∆ϕ 2 > 0 there is no real smooth geodesic between them. The geodesic approximation requires some continuation to the case of timelike-separated points.
This issue previously was tackled in the study of holographic global quench models [5][6][7] and other non-stationary backgrounds [45][46][47]. Several methods to continue the geodesic approximation to the timelike case were proposed in [7]: a specific analytic continuation, complexified geodesics and quasigeodesics. Here we use the method based on reflection mapping [46,47], which can be considered as a streamlined version of method of quasigeodesics adapted for topological quotients of AdS 3 in global coordinates.
The idea is the following [47]. Let us define a mapping R which acts on the boundary of AdS 3 in global coordinates: R : (τ, φ) → (τ + π, φ + π) ; (5.13) Figure 21. A reflection geodesic in relation to the wedge of the infalling particle.
This mapping has the key property that is useful for us. Consider a timelike interval ab on the boundary. If one acts with R on the point a, then the resulting interval a R b is spacelike. This is because in global coordinates φ ∼ φ + 2π. Now, since a R and b are spacelike-separated, we can consider the AdS 3 geodesic between points a R and b. We will call it the reflection geodesic. Its length can be found from the matrix formula (3.4). We parametrize the points a and b by matrices A and B, according to (2.7), by global coordinates (2.2), keeping in mind that they are located near the boundary at χ 0 → ∞: A = e χ 0 cos τ a + sin φ a sin τ a + cos φ a − sin τ a + cos φ a cos τ a − sin φ a ; (5.14) Using (5.13), we write the matrix parametrization for coordinates of reflected point a R : This formula is explicitly invariant and we will use it later for points a and b parametrized by BTZ coordinates as in ( In our case we deal with timelike-separated endpoints on the boundary, which in global coordinates are parametrized as a = (− π 2 , φ a ) and b = (τ, φ b ). To find the two-point correlation function between these points, we use the reflection prescription [47] and replace regular direct geodesic with the reflection geodesic for timelike-separated points: Note that we also take into account the Φ-factors of Lorentzian correlators, which are given in Appendix C by formulas (C.3,C.4,C.5). For our choice of endpoints, the reflection geodesic is a regular direct geodesic. This follows from the fact that by definition of reflection mapping τ a R = π. Meanwhile, in global coordinate the particle worldline is a light-like geodesic which goes from the point (0, 0) to the point (π, π) on the boundary, and the identification surfaces W ± are located beneath the worldline, see Fig.21 (the BTZ identification surfaces V ± are not shown for clarity). This means that all reflection geodesic will pass freely above the wedge, without intersecting it. Moreover, φ a R = φ a + π means that the point a R is located in the second external region with respect to the BTZ black hole identification, and thus the reflection geodesic is a direct geodesic that goes through the wormhole from the external region under consideration in our problem into the other exterior region which is located in the identification dead zone. Thus, the direct geodesic which goes through the wormhole is the only possible reflection geodesic. Now we can calculate its length using the formula (5.16) and parametrizing points a and b by BTZ coordinates according to (3.1). The result is (we subtract the diverging part): Note how it is different from the expression for the length of direct spacelike geodesic (3.10) by opposite sign of the expression under the logarithm, as expressed by the formula (5. 16. Now we substitute this expression into (5.17) to obtain the two-point correlation functions between timelike-separated points a = (0, ϕ a ) and b = (t, ϕ b ), t > 0: This expression establishes the fact that the leading temporal behavior of the two-point correlation functions after the quantum quench is the same as in thermal equilibrium (C.1) and does not carry any details about the infalling matter which forms the black hole. The same result was obtained in the case of Vaidya global quench [7]. Note, however, that if one takes the temporal coordinate of one endpoint before the quench, e.g. t a < 0, then the time dependence of the correlation function will be different [7,15].
In the context of our holographic bilocal quench model of colliding particles forming a black hole this calculation is out of the scope of this paper. However we can speculate that in such situation we would have to consider reflection geodesics which somehow pass through the identification W ± like spacelike crossing geodesic. Such reflection geodesic will carry information about the holonomy of the infalling particle and will lead to a more distinct time dependence of the correlation function.

Conclusions and discussion
In this paper we have investigated black hole formation from a collision of particles as a model of holographic thermalization in the boundary theory after the bilocal quantum quench in the framework of AdS 3 /CFT 2 -correspondence. This model holographically describes unitary time evolution of a non-homogeneous initial state with two excitations on the antipodal points of the CFT spacetime cylinder. We have probed the thermalization process by studying the non-equilibrium behavior of entanglement and two-point correlation functions. More specifically, we studied the dynamics of entanglement by using the time-dependent generalization of the RT prescription [3,4] to calculate the holographic entanglement entropy and mutual information, and investigating their time dependence, and we used the geodesic approximation and its extension to timelike case to study real-time two-point functions. The key results are the following.

1.
In the initial state of the boundary theory the subsystems which are located between the excitations exhibit constant thermal behavior of the entanglement entropy (4.3), since t = 0. Because we study the evolution of a pure state, their complements, which contain both initial excitations, are also thermalized, with HEE given by (4.4). We expect that this partial equilibrium in the initial state from the HEE perspective is one of the features of the bilocal quench in 2d specifically. One can think that this partial instant equilibration of the entanglement entropy is caused by the long-range entanglement induced just after the quench. The emergence of long-range entanglement is a characteristic feature of holographic local quenches, as discussed in [22] and also observed in [24]. It happens as a result of the fact that we use the time reversal-invariant classical bulk dynamics for holographic description of a non-reversible process, in which we instantly inject a large amount of energy. Meanwhile, the subsystems which at t = 0 contained one of the initial excitations, display non-trivial non-equilibrium dynamics of HEE (4.6) and two-point correlation functions (5.6). This picture suggests that the bilocal quench setup perhaps could work as a toy model for diffusion of two quantum fluids in a finite-volume vessel at finite temperature which are initially divided by two walls, which is analogous to a joining quench where two systems at the same finite temeprature are brought into contact in two points (for holographic consideration of joining quench of two thermal systems in a single point, see [32]). As explained in [24], the holographic description of the local joining quench is similar to the description of an operator local quench. In our case we observe the same on qualitative level.

2.
We have obtained the explicit formula (4.6) for the non-equilibrium behavior of the entanglement for subsystems which contained one of the initial excitations. The time dependence of HEE governed by the formula (4.6) is substantially similar to the global quench scenarios [5-9, 11-13, 15, 17-19, 21]. In our case the time evolution of HEE also display the early-time quadratic growth and the linear growth regime, see  This confirms the assumptions about universality of non-equilibrium growth of entanglement entropy with respect to the choice of initial out-of-equilibrium state in strongly-coupled systems proposed in previous work [8,9,16,17] for our special choice of initial state with two excitations. However, there is a strong difference in the character of transition to saturation: in our case the HEE undergoes a sharp transition to saturation with discontinuous first derivative, as the direct geodesic becomes dominating in the HRT prescription over the crossing geodesic. This behavior is similar to a first order phase transition 8 , which also happens in higher-dimensional global quench Vaidya models and in holographic models of formation of a charged black hole [8,17]. Contrary to this, in the global quench setup in two dimensions the saturation transition is smooth. The sharp transition to saturation also reveals that in our case we can observe a sharp emergent light cone (see Fig.17), which hints at similarity with quasi particle picture of entanglement spreading.

3.
Because of the similarities in non-equilibrium dynamics of HEE between the bilocal quench and global quench setups, we have observed that the universal characteristics of entanglement grwoth in the global quench models, namely the entanglement velocity v E and the emergent lightcone velocity v LC are relevant characteristics of the entanglement growth for the equilbration after the bilocal quench as well. We have confirmed that in our 2d model they equal to the speed of light, which is in agreement to their respective definitions in the case of 2d global quench. The fact that such velocity characteristics of entanglement propagation are meaningful in the case of certain local quenches was established in the work [29,30,32], where authors deal with local quenches which drive the system out of the initial thermal equilibrium. We have also discussed that the argument about the relation between v LC , v B and entanglement wedge subregion duality by Mezei and Stanford [16] also can be extended to our case. While it is significantly simplified by the fact that we work solely in (2 + 1)d bulk spacetime, and the RT surfaces are just direct geodesics, the interesting point is that in our case we deal not with small perturbations of the equilibrium state that fall into the black hole in the bulk, but with strong perturbations of vacuum which are responsible for the creation of the black hole and thermalization of the pure state. However, the bulk geometry of the particle creation itself in global coordinates (see Figs.3,4) looks similar to a setup when the particle 1 falls into a black hole. Such situtations have been considered in the work [26,27,29,30] as holographic duals to local quench in the thermal state. We can say that the bilocal quench holographic model is a model of thermalization of a pure state which shares similarities with both global quench models and local quench models of equilibration of perturbations of the thermal state.

4.
We have observed the relation between the scrambling time of the system t * = π 2 and the emergence of the entanglement shadow in the bulk. Namely, once the state evolves part the scrambling time, the identification wedge gets fully covered by the entanglement shadow, which cannot be probed by the entangling surfaces, which by that time are all direct geodesics and cannot reach into the bulk deeper than the radius of the entanglement shadow given by (4.25). In fact, after the scrambling time no subsystem can probe the infalling matter in the entanglement shadow, regardless of its particular features. That means that the same is true for e.g. an infalling Vaidya shell in AdS-Schwarzschild coordinates. We expect that the relation between the emergence of the entanglement shadow and scrambling time should be easily extendable to quantum quenches in higher-dimensional case.

5.
We have discussed that the linear growth of HEE at large distances from initial excitations is identified with the regime of memory loss and wave-like spreading of entanglement. We also have provided some evidence that, similarly to the holographic global quench case, the linear growth is governed by HRT geodesics that probe the interior of the forming black hole. Therefore, the behavior of HRT crossing geodesics hints that evidently there is a connection between the memory loss of details of the initial state during the thermalization and the interior of the black hole. We have also discussed that there are some hints that the thermalizing state loses detailed information about the initial configuration that can be seen purely from the geometry of the bulk spacetime. First, as we discussed in the end of section 2.2, the shape of the identification wedge, defined by holonomies of the colliding particles, approaches cylindrical at the limit t → ∞. The cusps at the particle worldlines, which break the rotational symmetry, are gradually smoothed as the state approaches thermal equilibrium.
The other point is related to which geodesics govern equilibrium and non-equilibrium regimes in our model. The non-equilibrium physics in our model is governed by crossing geodesics. The information about crossing geodesics is encoded in image geodesics. These geodesics are constructed using the identification isometry of the matter which forms the black hole, and can go anywhere in the global AdS 3 . Meanwhile, the equilibrium regime is described by direct geodesics which are completely restricted to the fundamental domain of the identification and do not carry any information about the infalling matter which describes the initial state. We can say that after the thermalization the system has lost all information about the infalling matter and about the orbit of its identification isometry which is located outside of the fundamental domain. All of these points lead us to a conclusion that we have a strong parallel between the memory loss during thermalization and information loss in the black hole formation, and that this parallel is prominent in the behavior of holographic entanglement entropy.
6. The holographic mutual information can exhibit a variety of different behaviors, depending on positions and sizes of the subsystems A and B. The most notable feature is the kinks in the time dependence of the mutual information, which happen when segments in the formula (4.26) thermalize. These kinks have the same origin as the sharp transition to saturation in the time dependence of HEE.

7.
Similarly to the HEE, the two-point correlation functions with spacelike-separated points in the minimal geodesic approximation exhibit non-trivial non-equilibrium behavior in the case when the collision line is located between the endpoints. However, the subleading contributions from winding geodesics appear with time for any choice of endpoints at late times, and they signify the approach to equilibrium of the system as a whole. The time dependence of the correlation functions after the quantum quench coincides with the thermal behavior, same as in case of the global quench. To verify this fact, we were able to use a particularly simple extension of the geodesic approximation to the timelike case based on the reflection mapping, which is suitable for dealing with AdS 3 topological quotients. It is remarkable that the reflection geodesics which provide the continuation for correlation functions have to go through the black hole interior into the second asymptotic region. This shows similarity between the reflection geodesics prescription and analytic continuation of Lorentzian amplitudes in the bulk background of the eternal BTZ black hole discussed in earlier work [67]. Contrary to the latter case, in our setup we deal with a pure state, and the second exterior region is completely unphysical, yet it still seems to play substantial role in description of late time behavior of holographic observables. Now let us discuss some possible future directions of the work.
(i ) First of all, we have not actually given the precise proof that the bulk spacetime with two colliding massless particles creating a black hole can indeed holographically describe the bilocal quench, realized by two operator excitations on the CFT side. To do that, one would need to calculate the 6-point correlation functions using CFT techniques and match the results to our holographic computation. We expect that one could use the method similar to that of [15], where such computation was performed for the boundary dual of the Vaidya global quench, and also to other earlier work, e.g. [24,66,68]. Specifically, one could consider a 6-point function of two light operators in the background generated by two "heavy" operators, and calculate it in the approximation of the vacuum Virasoro conformal block. The main difference from the calculation based on the monodromy method in [15] would be that in our case we have a finite amount (namely two) operators which produce excitations represented by massless particles in the bulk. Because of this, the limiting procedure which would allow for perturbative solution of monodromy equations at large central charge, will be different. The vacuum Virasoro conformal block of the correlation function can holographically be obtained from the length of the minimal geodesic [24,66,69], so we expect that such CFT computation would reproduce our holographic expressions for correlation functions (5.3) and (5.6) in the leading order of the semiclassical expansion. Also, throughout the paper we have been implicitly assuming that the energy of initial excitations must be high enough in order to form a BTZ black hole instead of just a static conical defect. These two cases are easily identified and separated, if one considers the collision of particles in the center of mass frame, as shown by the formula (2.25) (see [34]), but the black hole rest frame bulk geometry doesn't have a well-defined continuation to the case of the formation of a conical defect, so in our holographic computations this energy threshold is hidden. The CFT computation of correlators could clarify this issue as well.
(ii ) An interesting related question is study of the bilocal quench away from the holographic limit c → ∞. Specifically, it could be interesting to see how the entanglement scrambles in examples when the cental charge is finite. Such studies were conducted for local quench [23,24], as well as global quench [28] scenarios. The key observation in the latter case is that in some cases at finite central charge there are memory effects, which interrupt the equilibrium saturation of HEE and mutual information at late times. We could expect that in the bilocal quench scenario one would also observe similar memory effects for finite (small) c theories. Moreover, one could consider the N -local quench protocol generalization for more than 2 excitations in order to find out how the late time memory effects in HEE depend on the fraction of initial excitations located inside the given subsystem, or simply speaking how much memory do these memory effects carry about the initial state. However, of course the study of quenches produced by multiple localized excitations in a CFT at finite c seems to be a challenging problem.
(iii ) Another related question is the study of the 1/c corrections. We have observed that HEE in our model exhibits sharp transition to saturation. We could expect that the perturbative semiclassical corrections in 1/c would smooth out the transition. Even more interesting are the non-perturbative, e −c corrections. These are considered to be the corrections which should restore the information lost in black holes in holographic correlation functions [70]. In case of our setup that would mean that the e −c corrections to HEE would carry information about the topological identification even after the saturation. That fact could be used to investigate the question: just how much the entanglement entropy actually knows about the bulk geometry beyond the semiclassical approximation and out of the entanglement wedge?
(iv ) A possibly promising direction of further study is the higher-dimensional generalizations of the bilocal quench to the cases of AdS 4 /CFT 3 and AdS 5 /CFT 4 . The holographic dual for thermalization after bilocal quench in those cases should be the AdS spacetime with colliding shockwaves which create a black hole [38]. From the present work, as well as from previous work on local quenches in 3d [29,30] and bounds on the entanglement propagation [16,17] we can expect that the non-equilibrium dynamics of entanglement will be more rich in higher dimensions. Also, the bilocal quench in AdS 5 /CFT 4 could serve as a viable more realistic holographic model of thermalization in heavy ion collisions [38,65].

B Geodesics in the BTZ geometry
We start from the BTZ coordinate patch of the AdS 3 spacetime: where t ∈ [0, +∞), r ∈ (R, +∞), and for now ϕ ∈ R. The geodesics in this geometry are described by the following formulae [7]: The constants Γ ± . R are subject to the restriction 0 < Γ 2 − < R 2 < Γ 2 + ; (B.5) One can express the variables as a function of r ∈ [Γ + , ∞[: We are interested in geodesics with endpoints on the boundary. From the equations (B.3) and (B.4), the separation between two boundary endpoints (ϕ a , t a , r a = ∞) and (ϕ b , t b , r b = ∞) is given by In this paper we use BTZ coordinates in the spacetime with BTZ identification generated by the holonomy (2.19). In the BTZ coordinate patch it periodizes the angular coordinate: ϕ ∼ ϕ + 2πn, where n ∈ Z. We choose the fundamental domain such that the worldline of the first particle is in the middle of the domain, and the worldline of the second particle is at the periodically identified boundary, i. e. we set ϕ ∈ [−π, π).
Depending on the value of the difference of ∆ϕ = ϕ 2 − ϕ 1 , a geodesic in the BTZ patch of the pure AdS 3 space without topological defects is pulled back to either a direct (|∆ϕ| < π) geodesic, or a geodesic which winds around the horizon ((|∆ϕ| > π) on the topological quotient. The above formulas for the geodesics hold when pulled back to the BTZ coordinates, but with addition of the term 2πn to the ∆ϕ. With that in mind, let us now focus on geodesics in BTZ coordinates, taking into account all windings. The equations (B.9-B.10) can be used to express the integration constants Γ ± through the coordinates of the endpoints. From the equation (B.2) it is clear that the Γ + has the meaning of the closest distance along the radial coordinate to the horizon, into which the geodesic with given endpoints can approach, i. e. the depth of the geodesic in the bulk. The formula (B.11) shows that winding geodesics always approach to the horizon closer than direct geodesics. This fact is key to our understanding of dynamics of the pole structure of non-equilibrium correlators, as explained in the main text. The (regularized) length of a geodesic can be found to be ∆λ = λ + (r → r 0 ) − λ − (r → r 0 ) = 2 ln in terms of a regularized AdS boundary at r = r 0 . In the limit r 0 → ∞ this reduces to ∆λ = 2 ln 2r 0 (B.14) In terms of coordinates of the endpoints, this formula leads to the expression where r 0 is the near-boundary cut-off. From this formula, it is clear that the length of any winding geodesic is always larger than the length of any direct geodesic. Therefore, the winding geodesics in BTZ geometry will always give subleading contributions to holographic correlation functions for any placement of endpoints.

C Lorentzian correlators in thermal equilibrium
Consider real-time two-point correlation functions of a scalar operator O of scaling dimension ∆ in holographic (1 + 1)d CFT at finite temperature T > R 2π above the Hawking-Page threshold, so that R > 1. In this case the holographic dual bulk spacetime is the static BTZ black hole. The Wightman two-point correlator than reads [62,63]: (C.1) Let us assume for definiteness that t > 0. Then in analogy to the zero-temperature vacuum case of CFT on a cylinder [46,47,64], this expression can be rewritten as following sum over images: where the phase factor is introduced Φ ∆ W (R; t, ϕ; n) = exp (−iπ∆θ(cosh Rt + cosh[R(ϕ + 2πn)])) .

(C.3)
If initial points are spacelike-separated, then all phase factors in this expression equal to 1, and thus the Wightman correlator for spacelike-separated points can be interpreted as sum over all geodesics in the BTZ black hole geometry with metric (2.5). The n is the winding number. The leading contribution to the sum is given by the direct geodesic with n = 0, and winding geodesics give exponentially suppressed contributions to the correlator. If points are timelike-separated, then the contribution with n = 0 is interpreted as contribution from the reflection geodesic times the exponential factor e −iπ∆ . The same goes for contribution from image points. One can rewrite other real-time two-point functions in similar manner [7,46,47]. The difference is in the phase factors, which can be recovered using common QFT definitions of Green's functions through Wightman correlators. For example, for the retarded Green's function one obtains