Dynamical meson melting in holography

We discuss mesons in thermalizing gluon backgrounds in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \mathcal{N} $\end{document} = 2 super-symmetric QCD using the gravity dual. We numerically compute the dynamics of a probe D7-brane in the Vaidya-AdS geometry that corresponds to a D3-brane background thermalizing from zero to finite temperatures by energy injection. In static backgrounds, it has been known that there are two kinds of brane embeddings where the brane intersects the black hole or not. They correspond to the phases with melted or stable mesons. In our dynamical setup, we obtain three cases depending on final temperatures and injection time scales. The brane stays outside of the black hole horizon when the final temperature is low, while it intersects the horizon and settles down to the static equilibrium state when the final temperature is high. Between these two cases, we find the overeager case where the brane dynamically intersects the horizon although the final temperature is not high enough for a static brane to intersect the horizon. The interpretation of this phenomenon in the dual field theory is meson melting due to non-thermal effects caused by rapid energy injection. In addition, we comment on the late time evolution of the brane and a possibility of its reconnection.

1 Introduction and summary

Introduction
The gauge/gravity duality [1][2][3] enables us to tackle strongly-coupled gauge theories by using classical gravity. The deconfined phase in the finite-temperature gauge theories is described by the gravity dual in the presence of black holes [4]. This viewpoint has been applied to studying characteristics of the quark-gluon plasma, observed at the RHIC and the LHC experiments. It has been predicted that in the strongly coupled plasma, the shear viscosity over entropy is universally very small [5][6][7][8][9][10]. The gravity computations are related with fluid dynamics [11][12][13]. An advantage of using the gravity dual is that time-dependent systems can be easily handled. Thermalization and its out-of-equilibrium dynamics are of interest in real-world processes in experiments. In the gravity dual, thermalization is described by black hole formation [14]. 1 The timescale for the thermalization is typically fast, and the plasma evolution can be understood from hydrodynamic computations [18][19][20][21][22][23]. On the other hand, very-far-from-equilibrium effects can be significant in the very beginning of the fast thermalization. We would like to consider such a situation in QCD using the gravity dual.
In this paper, we will study time-dependent dynamics of mesons in fast thermalization in the gauge/gravity duality. In particular, we would like to focus on the initial-stage effects in the thermalization. For that purpose, we will do numerical computations to solve the partial differential equations describing time evolution in the gravity dual. As the setup, we will consider the N = 2 SQCD realized by D-branes.
The N = 2 SQCD provides a simple ground for catching qualitative features of QCDlike theories. This theory is realized by N c D3-branes and N f D7-branes, where the D7-branes supply fundamental quarks. In the strong coupling in the large-N c limit, the D3-branes are replaced with the AdS geometry, and we probe the curved background with D7-branes in the limit of N f ≪ N c [24]. In this paper, we consider the case when N f = 1. 2 In finite temperatures, the gravity background is the AdS-Schwarzschild black hole with flat horizon [4]. There is no confinement/deconfinement phase transition in the color sector of the N = 2 SQCD, and the gluons are deconfined in finite temperatures. In the flavor sector, however, if we introduce a nonzero quark mass to explicitly break the global U(1) symmetry of the theory, we can find two phases with or without stable mesons, where the parameter is the ratio of the quark mass and the temperature.
The phases are characterized by the embeddings of the D7-branes into the black hole geometry [25]. There are two kinds of embeddings when the system is static. Let us consider units where the quark mass is unity; the parameter is then the temperature determined by the black hole size. In low temperatures, the D7-branes do not touch the black hole horizon, and the fluctuations on the D7-branes give stable mesons [26]. This case is called the Minkowski embedding. In high temperatures, the D7-branes end on the black hole, and the quasi-normal modes on the D7-branes correspond to melting mesons in the thermal JHEP04(2014)099 background [27,28]. This case is called the black hole embedding. The phase transition has been studied in detail in [29][30][31][32][33][34]. 3 Our interest is what happens when we consider dynamical embeddings of the D7-brane in dynamical backgrounds. To address this question, we will compute time evolutions of the D3/D7 system by numerically solving the equations of motion, which are given by partial differential equations. 4 We consider dynamical embeddings of the probe D7-brane in a D3brane background corresponding to a system thermalizing from zero to finite temperatures.

Setup
Let us introduce our setup in a more concrete manner. In this paper, we take the Vaidya-AdS 5 × S 5 spacetime as the gravity background. In the gravity dual, thermalization is realized by the gravitational collapse and the subsequent black hole formation. The Vaidya-AdS metric is one of the simplest geometry of this kind, obtained for instance as a result of planar-shell collapse [14] or the collapse of null dust. 5 This collapse in the bulk gravity will be identified with a spatially-homogeneous injection of a fixed amount of energy per unit volume in the boundary field theory. We focus on the simplest case of the Vaidya-AdS metric to clarify basic aspects of the dynamical behaviors of the holographic mesons, though in general we could consider more complicated thermalization patterns such as those realized by local heating caused by point particle infalling into the bulk black hole [42] or an anisotropic wake induced by a time-dependent deformation of the boundary metric [43]. The thermalization time scale is typically fast for the Vaidya-AdS metric [44][45][46], and hence it would be suitable to focus on the effect of rapid thermalization.
The time evolution we consider is sketched as follows. In the Vaidya-AdS 5 × S 5 spacetime, it is convenient to use the ingoing Eddington-Finkelstein time V as the bulk time coordinate, and V coincides with the time of the field theory at the AdS boundary. Initially for V < 0, the spacetime is locally pure AdS. We then inject energy from the AdS boundary at V = 0 to form the black hole that settles down to the Schwarzschild-AdS 5 black hole in the final state. The initial static embedding of the D7-brane before the energy injection is given by the analytic solution in the pure AdS spacetime. After the injection, the brane starts to change its embedding shape. We will numerically compute the evolution of the branes. Figure 1 gives schematic illustration when the final mass of the black hole is sufficiently large and the brane intersects with the black hole horizon.
Here, we would like to emphasize subtleties in treating dynamical backgrounds of this kind. When we discuss non-stationary phenomena in the gravity dual description, especially when the event horizon exists, we should pay careful attention to the causality and non-uniqueness of the time foliation of the bulk spacetime, unlikely to static or stationary cases. In our study, such subtleties associated with the bulk causality become manifest when we try to answer questions such as the following: • How should we determine the phase of the brane embeddings at a given time? JHEP04(2014)099 Figure 1. The Penrose diagram of the Vaidya-AdS 5 spacetime. The initial data for the bulk brane configuration is given by the static embedding in the pure AdS. The null dust is injected to generate the bulk black hole, and H + is the event horizon of the black hole. The world-volume of the brane is shown by the light gray region. In this figure, the final state of the D7-brane embedding is the black hole embedding.
• How can we know the moment at which the phase transition takes place in the boundary theory?
In static cases, the former question becomes almost trivial: we can distinguish the BH and Minkowski embeddings just by seeing whether the brane intersects with the event horizon or not. In general dynamical cases, however, this naive consideration is no longer suitable. Firstly, there is no definitely preferred time foliation, and this fact implies that the brane configuration at a given bulk time depends on bulk time foliation and is not unique. In addition, the part of the brane world-volume captured in the event horizon cannot be observed by any observers outside the black hole, since that part lies in the causal future of those observers. The information of that part does not affect the dynamics of the other part of the brane outside the black hole, as well as time-dependent boundary observables determined by the brane dynamics. This observation tells us that some naive notions based on the static cases are no longer rigorous in the dynamical cases. We address this issue in more detail in section 4. In the gauge/gravity duality, it is the dual boundary operators that contains information of the boundary field theory. In our study, we will compute the time evolution of the operators by solving the equation of the motion of the brane. This approach guarantees that the results are not affected by bulk time foliations.

Summary
Below, we summarize properties of the dynamical D7-brane embeddings clarified by our numerical analysis. Our model possesses two parameters, which are the final black hole mass and the mass injection speed.

JHEP04(2014)099
We find there are three cases of the dynamical D7-brane embeddings depending on these parameters. Two of them are analogous to static embeddings. If the final black hole mass is sufficiently small, the brane never intersects with the event horizon, and the embedding is the Minkowski embedding with non-decaying oscillations. In the dual field theory, the quark condensate c(V ) oscillates around the equilibrium value. This is the sub-critical case in which the phase transition does not occur since the final temperature is sufficiently low. If the final mass is sufficiently large, the brane intersects with the horizon, and relaxes to the black hole embedding. This is the super-critical case in which the final temperature is sufficiently high for the phase transition to occur and the mesons are melting.
We find that another kind of dynamical embedding appears when the final black hole mass is slightly smaller than the phase transition value. There is no static BH embeddings in this regime. Despite this fact, the brane dynamically touches the horizon if the injection is sufficiently fast: the rapid injection induces a large brane motion, and it drives the brane into the horizon. Roughly speaking, a part of the brane falls into the black hole because of inertia. We call it the overeager case since the brane is made "overeager" to plunge into the black hole under the influence of the dynamical environment.
After intersecting, the brane configuration on constant V surfaces moves in order to escape from the bulk black hole since the final equilibrium configuration in the future cannot be the black hole embedding. The brane configuration, however, tends to be singular within a finite time as the intersection locus approaches the pole where the S 3 wrapped by the D7brane shrinks to zero size. It is expected that if stringy or quantum effects are taken into account there, the brane will reconnect near that locus, and then go back to the Minkowski embedding. We show schematic illustration of this process in figure 2(c), together with the sub-critical and super-critical cases in figures 2(a) and 2(b), respectively. We argue that this phenomenon may be the gravity dual of quarks recombining into mesons in the boundary field theory.
The remaining of this paper is organized as follows. In section 2, we will review the static embeddings to fix notations. In section 3, we describe the setup for our dynamical computation. Results of numerical computations are shown in section 4. Section 5 is devoted to discussions on the boundary theory interpretations and future directions. Technical details for numerical computations are provided in appendices.

Static embeddings
We start from reviewing the static embeddings of the D7-brane to the Schwarzschild-AdS 5 × S 5 spacetime in order to fix notations. The metric of the Schwarzschild-AdS 5 × S 5 background can be given as and L is the AdS radius. The AdS boundary is at z = 0, while the event horizon is at z = L 2 /r h . We use the same worldvolume coordinates as the target space coordinates themselves such that σ 0 , σ 1 , · · · , JHEP04(2014)099  σ 7 = (t, z, x 3 , Ω 3 ). The embedding of the D7-brane is specified by its position (φ, ψ) in the transverse space as a function of z: φ = Φ(z) and ψ = 0. Note that we can set ψ = 0 without loss of generality thanks to the U(1)-symmetry generated by ∂ ψ . The induced metric on the D7-brane is given as The embedding of the D7-brane is determined by the Dirac-Born-Infeld (DBI) action. In the absence of the world-volume field strength, the action is where h = det h ab . The equation of motion is a second-order ordinary differential equation for Φ(z), It is also convenient to introduce new bulk coordinates (w, ρ) = L 2 z −1 sin φ, L 2 z −1 cos φ . The new embedding function can then be a function of ρ, w = W (ρ). In practical numerical computations, we solve the equation for W (ρ) obtained by rewriting eq. (2.4) in terms of W (ρ).

JHEP04(2014)099
The asymptotic behavior of the field at the AdS boundary gives information of the dual field theory operators. Near the AdS boundary ρ = ∞, the asymptotic solution behaves as 6 where the two integral constants m and c are related to the quark mass M q and condensate O m as [25,[29][30][31][32]35] where ℓ s is the string length and g s is the string coupling constant. O m is the quark bilinear associated with its supersymmetric counterparts (see refs. [34,47] for further details).
Hereafter, we ignore the proportional constants and refer to m and c as the quark mass and the quark condensate for notational brevity.
In the bulk, we also need to impose boundary conditions at where the D7-brane terminates. The boundary conditions depend on the topology of the embedding, that is, whether the D7-brane intersects the black hole horizon or not. When the brane intersects with the black hole, we impose regularity at the horizon. The asymptotic solution is then obtained as where ρ H (0 < ρ H ≤ r h ) is the value of ρ at the horizon found in the (ρ, w)-plane. On the other hand, when the brane does not intersect with the black hole, the brane terminates at the point where S 3 wrapped by the brane shrinks to zero size at a pole of S 5 : Φ = π/2, corresponding to ρ = 0. The regular asymptotic solution at ρ = 0 is given by where W P ≡ W (0) > r h is the value of W at the pole. Solving the equations of motion by using either the boundary condition (2.7) or (2.8) when the brane intersects with the black hole or not, respectively, we obtain a sequence of embeddings. In figure 3, we show the solutions for several parameter values. For m/r h ≃ 0.92, we can see that both the black hole and Minkowski embeddings are allowed, and they are characterized by different values of the quark condensate c for a given m. We plot c/r 3 h as a function of r h /m in figure 4. c/r 3 h becomes multivalued for 1.0802 < r h /m < 1.0913. This implies that there is a first-order phase transition between the Minkowski and black hole embeddings, and it can be confirmed by comparing the free energy [32,34]. 6 In the notation of ref. [34], W is expanded as W = 2 −1/2 r hm + 2 −3/2 r 3 hc /ρ 2 + · · · . In this notation, the quark mass and condensate are written as Mq = √ λTm/2 and Om = − √ λN f NcT 3c /8, where T ≡ r h / πL 2 is the Hawking temperature and λ is the 't Hooft coupling.   In this section, we consider the very-far-from-equilibrium dynamics of the probe D7-brane in a dynamical spacetime focusing on the thermalization process. Thermalization in the boundary theory is realized in the gravity dual by the black hole formation due to gravitational collapse in the bulk AdS space. In this paper, we use the Vaidya-AdS spacetime as the dynamical background with the black hole formation. Hereafter, we set units where the AdS radius is unity, L = 1. The Vaidya-AdS 5 metric is given by

JHEP04(2014)099
where the mass function M (V ) is a free function representing the Bondi mass (density) of this spacetime. This metric is an exact solution of the following five-dimensional Einstein equation with null dust, where k µ denotes the ingoing null vector defined by k µ dx µ = −dV . The null dust injected from the AdS boundary infalls into the bulk, and then the event horizon is formed. When the mass function M (V ) stops depending on V , namely, when the energy injection ceases, the spacetime is locally isometric to the Schwarzschild-AdS 5 with a planar horizon. In our calculations, we chose a C 1 function for M (V ) as where M f is the final mass of the spacetime, and ∆V determines the duration of the energy injection. Note that the final radius of the event horizon, r h , is given by M f = r 4 h . Uplifting the solution to 10-dimensions, we obtain Vaidya-AdS 5 × S 5 spacetime: We consider the probe branes on this 10-dimensional background geometry. Figure 1 gives schematic illustration of our setting, as explained in section 1.2.

Evolution equations of the D7-brane
Having explained the background spacetime, we proceed to consider evolution equations of the D7-brane to solve dynamical embeddings numerically. Dynamics of a probe Dbrane is described by the DBI action. It has the diffeomorphism invariance with respect to transformations of world-volume coordinates. We should choose some useful coordinates, namely fix the gauge, for solving the evolution equations numerically. In order to derive equations of motion taking advantage of coordinate transformations on the brane, here we would like to begin with the Polyakov type action instead of the DBI action. Of course, once we know the appropriate ansatz, we could simply start from the DBI action to derive the equations of motion. The result is the same as that obtained from the Polyakov type action since the computation is on-shell. The Polyakov type action for the D7-brane in the absence of the world-volume field strength is given by where X µ is the brane collective coordinate and g µν is the metric in the target space. Note that λ is an arbitrary non-zero constant. Variating the action with respect to X µ and γ ab ,

JHEP04(2014)099
we obtain the equations of motion of the D7-brane as 7 where D a is the covariant derivative with respect to γ ab , and Γ µ ρσ is the Christoffel symbol with respect to g µν . The second equation implies that γ ab can be regarded as the induced metric on the brane up to a factor 6λ −1 . One finds that the wave equations (3.7) for X µ (σ) on the world-volume with the induced metric γ ab are evolution equations, and the equations (3.8) which determine components of the induced metric are constraint equations. In the following, we will study the D7-brane embeddings on the Vaidya-AdS background spacetimes solving eqs. (3.7) and (3.8).
There are eight world-volume coordinates σ a (a = 0, 1, · · · , 7) on the brane. For six of them, we use the target space coordinates themselves as σ 2 , · · · , σ 7 = ( x 3 , Ω 3 ) because we assume that the brane has the same symmetry as the background bulk spacetime in those directions. For the other two coordinates, we introduce the double null coordinates σ 0 , σ 1 = (u, v), and then the brane collective coordinates are parametrized as where we set ψ = 0 without loss of generality because of the U(1)-symmetry generated by ∂ ψ . The (u, v)-coordinates are chosen so that the ∂ u and ∂ v are null generators. Then, the brane induced metric is written as In the following calculation, it is convenient to use a new variable Ψ Note that there are still residual coordinate freedoms on the world-volume, This residual gauge freedom will be fixed by boundary conditions and an initial condition. Substituting eqs. (3.9) and (3.10) into eq. (3.7), we can obtain evolution equations for V (u, v), Z(u, v) and Ψ(u, v). 7 Using eq. (3.8) and eliminating the auxiliary field γ ab from the action (3.6), we obtain Dirac-Nambu-

JHEP04(2014)099
The equations we will solve are summarized as follows. The evolution equations for V , Z and Ψ are obtained as From the requirements that the (u, u) and (v, v)-components of eq. (3.8) vanish, we obtain the constraint equations, Note that the evolution equations guarantee the conservation of the constraints: Consequently, once we set an initial data satisfying the constraint equations, we may solve only the evolution equations as the constraints will be kept satisfied. The details of the numerical method to solve these equations are summarized in appendix A. Eliminating u and v, we can regard Ψ as a function of Z and V . Near the AdS boundary Z = 0, Ψ(V, Z) is expanded as The constant m and the function c(V ) are related to the quark mass and condensate as in eqs. (2.6). Solving the evolution equations, we can determine the time dependence of the quark condensate c(V ).

Boundary conditions and initial data
In general, two timelike boundaries appear on the world-volume of the brane: one is the AdS boundary Z = 0, and the other is the pole Φ = π/2 at which the radius of S 3 wrapped by the D7-brane shrinks to zero. When a numerical domain to be solved contains these boundaries, we need to impose boundary conditions there. Since, however, the evolution

JHEP04(2014)099
equations become singular at each boundary, we should fix the location of the boundary in the world-volume (u, v)-coordinates for numerical convenience.
As mentioned previously, we have the residual gauge freedom, which is a transformation from the original null coordinates to other ones (3.12). This is nothing but a conformal transformation between two-dimensional flat spacetimes. Using it, one can fix the location of the AdS boundary on the world-volume coordinates as u = v, i.e. Z| u=v = 0, or that of the pole as u = v + π/2, i.e. Φ| u=v+π/2 = π/2, or both. They become coordinate conditions on the world-volume domain.
Note that fixing the residual gauge is equivalent to choosing a certain conformal factor, which behaves as a harmonic function obeying a (1 + 1)-dimensional free massless field equation. 8 Hence, if we give a set of coordinate conditions at the boundaries and the initial surface, they will provide necessary boundary and initial conditions for the conformal factor. Those conditions will fix the residual gauge completely, and the null-coordinate system will be automatically determined on the computational domain.

Boundary conditions at the AdS boundary
Now, we derive boundary conditions at the AdS boundary. Since we will focus on observables there, the AdS boundary should be always contained in our computational domain and located at u = v throughout our calculations. The boundary conditions for Z and W are trivial: Z| u=v = 0 and Ψ| u=v = m. We choose the quark mass as time-independent, namely m is constant. Remaining boundary conditions on the AdS boundary can be derived by requiring regularities of the evolution equations at Z = 0. We expand the variables near the boundary as Substituting them into the evolution equations (3.13)-(3.15), we obtain the asymptotic solution as where · ≡ d/dv. They giveV Hence, integrating Z ,u | u=v along the AdS boundary, we can determine the time evolution of V | u=v . (The second equation is used to determine Z ,u at the boundary. For details, see appendix A.) Note that while the regularities yield V ,u | u=v = 0 or V ,v | u=v = 0, we have taken V ,u | u=v = 0 because the time direction in our setup leads to V ,v | u=v > 0. We can confirm that V ,u | u=v = 0 and V ,v + 2Z ,v | u=v = 0 (it is equivalent to the first equation of (3.21)) are consistent with the constraint equations at the AdS boundary. 8 We denote a conformal transformation asγ ab = e ω γ ab . In two dimensions, the scalar curvature is where ∇ is the covariant derivative with respect to γ ab . Now, since bothγ ab and γ ab are flat, namely R = 0 andR = 0, the conformal factor must be determined by a harmonic function ω satisfying ∇ 2 ω = 0.

Boundary conditions at the pole
Next, we consider the boundary conditions at the pole, at which the radius of S 3 wrapped by the D7-brane shrinks to zero. The brane may intersect with the event horizon in our calculations, and in that case we do not need to impose the boundary conditions at the pole since the AdS boundary is not in the causal future of the pole on the brane. Strictly speaking, when our computational domain does not contain this boundary, we do not need any boundary conditions there. If we need the boundary conditions, we will locate the pole Φ = π/2 at u = v +π/2 by using the residual gauge. Hence, one of the boundary conditions is ZΨ| u=v+π/2 = π/2. The others are given by regularities of the evolution equations at Φ = π/2 as at u = v + π/2. They are merely Neumann boundary conditions at the pole.

Initial data
Finally, we explain the initial data for our calculations. Before the energy injection V < 0, the background spacetime is pure AdS, and the brane is static there. We have the exact solution of the static embedding given by Z = m −1 sin Φ, i.e., Ψ = Z −1 arcsin(mZ). Solving the constraint equations (3.16) and (3.17), we obtain where φ is a free function corresponding to the residual coordinate freedom, and V ini is an integration constant. At the initial surface v = 0, the solutions become (3.24) where we have set φ(0) = 0, and then V (0, 0) = V ini < 0 is the initial time at the AdS boundary.
The residual gauge φ(u) corresponds to the coordinate freedom for u on the initial surface, and then we can choose an appropriate function for φ(u) arbitrarily. When the brane does not intersect with the event horizon, we can simply choose φ(u) = u. In this case, the pole Φ = π/2 can be fixed at If the brane intersects with the event horizon at the initial surface, or if the brane is close to the event horizon, the above choice φ(u) = u is not good. The reason is as follows. Each u-constant line represents an out-going null geodesic on the brane worldvolume. Roughly speaking, δφ ≃ φ(u + δu) − φ(u) is a proper distance between two null geodesics described by u + δu and u at the initial surface. If these out-going null geodesics pass through the vicinity of the event horizon towards the AdS boundary, the interval of each arrival time at the AdS boundary, δV ≃ V (u + δu, u + δu) − V (u, u), will become very large compared to the proper distance δφ at the initial surface. 9 This behavior essentially JHEP04(2014)099 originates from the gravitational red-shift effect near the horizon. As a result, if we choose φ(u) = u, V blows up rapidly at late time and then the numerical calculations easily break down.
To avoid this problem, we fine-tune φ(u) so that V and v are synchronized at the AdS boundary, i.e., V ,v | u=v = 1 (or V | u=v = v + V ini ). This implies that dφ/du becomes exponentially small in the vicinity of the event horizon. In appendix A, we explain how to fine-tune φ(u) in our numerical calculations.
Note that the condition V ,v | u=v = 1 is another boundary condition at the AdS boundary imposed in addition to the one that fixes the boundary position at u = v. This means that we do not fix the location of the pole Φ = π/2 on the world-volume coordinates. However, this is not a problem because the pole is not contained in our computational domain and then we do not need to impose the boundary conditions there.

Results
In this section, we show results of our numerical calculations. We use units in which m = 1. 10 The parameters are the final mass of the black hole (in other words, the final temperature) and the duration of the energy injection, which are respectively described by r h and ∆V .

Sub-critical case
In this subsection, we show numerical results for the small final mass case in which the D7brane does not touch the horizon and remains in the Minkowski embedding after the black hole formation. This is the sub-critical case in which the final temperature is sufficiently low and the phase transition does not occur. Since there is no dissipation when the brane is in the Minkowski embedding, excitations on the brane caused by dynamics of the bulk geometry will remain without decay. 11 Figure 5 shows the time evolution of the quark condensate c and the Bondi mass M with respect to the boundary time V . Note that all of these quantities are well-defined at the AdS boundary and are directly related with observables in the boundary theory. The Bondi mass M corresponds to the energy density of the CFT (gluon) fluid. We can recognize seemingly periodic oscillations after the energy injection. These correspond to excitations of stable mesons, namely normal modes of the brane fluctuations. Indeed, the power spectrum of c(V ) shown in figure 6 obviously indicates that the excitations have discrete spectrum, similarly to the case of linear perturbations on the static Minkowski embedding. Note that we have solved non-linear evolution equations of the brane rather 10 In our system, there is a scaling symmetry: ∆V → k∆V , m → m/k, and c → c/k 3 , where k is a constant. Using the scaling symmetry, we can set m = 1 without loss of generality. 11 In this paper, since we consider the probe D7-brane assuming N f /Nc ≪ 1, effects of its backreaction are not present. If one takes into account sub-leading effects in terms of N f /Nc, it is expected that the excess of the energy gained by the injection will be radiated through closed string modes such as gravitons over a long time period. . It is obvious that the excitations have discrete spectrum, namely normal modes. The first peak in the spectrum is located at ω = 2.8, and it corresponds to the lowest normal mode of the D7-brane oscillations.
than linear perturbations around an equilibrium brane configuration, and the discrete spectrum is a nonlinear version of the normal mode spectrum for linear fluctuations. Figure 7(a) shows time evolutions of c(V ) for various ∆V . It turns out that the shorter the injection duration ∆V is, the harder the condensate c(V ) oscillates. This implies that the larger non-adiabaticity is, the larger the excitations on the brane are. To study the fluctuations in more detail, we show the total power of fluctuations of the condensate, figure 7(b). We calculated c(V ) for 0.1 ≤ ∆V ≤ 25 setting the final horizon radius to r h = 0.5, and obtained P tot (∆V ) from c(V ) using V 0 = 30 and V 1 = 50. We confirmed that the numerical results are almost independent of V 1 as long as we keep it larger than 50. The numerical results imply that P tot (∆V ) shows a damped oscillation whose amplitude behaves as ∆V −4 , and appears to converge into the finite maximum value in the ∆V → 0 limit.
Before moving to the next subsection, we make some comments on the relationship of these results to a harmonic oscillator model. We find that features of P tot (∆V ) can be well  reproduced by a harmonic oscillator with an external force described by as the external force f (V ). The total oscillation power in the final state of this model is given by P tot (∆V ) = C(1+cos(ω∆V )) (π 2 −ω 2 ∆V 2 ) 2 , where C is a normalization constant independent of ∆V . 12 Fitting this analytic formula to the numerical data, we obtain ω = 2.8. This value coincides with the lowest normal mode frequency of the brane oscillations (see figure 6).
Precise behavior of P tot (∆V ) depends on the function form of M (V ). As an example, in figure 8 we show P tot (∆V ) for M (V ) given by . the external force applied to the D7-brane. It may be interesting to pursue these ideas in more detail. Also, it may be interesting to study further on behaviors in the ∆V → 0 limit to give an estimate on the largest amplitude realized by a rapid energy injection. Our numerical results and the analytic results based on the harmonic oscillator model in this context suggest that P tot (∆V ) converges to finite constants in the ∆V → 0 limit. It may be interesting to generalize results of ref. [48], which focuses on abrupt quenches in the holographic setup, to the D3/D7 system and compare them with our results.

Super-critical case
In this subsection, we show numerical results for the large final mass case in which the D7brane is drawn into the event horizon and its configuration changes from the Minkowski embedding to the BH embedding. This is the super-critical case in which the final temperature is sufficiently high for the phase transition to occur.
In figure 9, we show a typical behavior of the quark condensate c and the black hole mass M in this case, setting the energy injection duration to ∆V = 1. The quark condensate c gradually evolves as the black hole mass increases due to the injection, and eventually settles down to the final value by V ∼ 4. This final value coincides with c of the static BH embedding for the same bulk black hole mass. It turns out that a quasi-normal mode is excited in the ring-down phase. This fact implies that the phase transition to the BH phase has taken place and the brane comes to equilibrium.
In figure 10, to focus on the brane world-volume in the bulk, we show snapshots of the brane configuration on time slices characterized by V = constant. Note that we define the bulk coordinates on the time slices as (w, ρ) = z −1 sin φ, z −1 cos φ . The radius of the event horizon is evolving larger, and we can see that the brane eventually intersects with the horizon. At late time the configuration of the brane will coincide with the static BH embedding.
When we consider dynamical phenomena in the bulk, especially when an event horizon exists, we should pay careful attention to the causality. In the static case, two meson phases of the boundary field theory can be distinguished by seeing whether the brane configuration at any moment is intersecting with the horizon or not. However, such a naive definition can not be used in dynamical cases by the following reason. A temporal configuration of the brane, which is the intersection of the brane world-volume and a certain time slice of the bulk, depends on the given time slice. Since there is no definitely preferred time slice in general spacetime, we cannot uniquely identify the phase of the brane embedding according to the temporal configuration at that time. To see it, let us consider the following example for time slices determined by V = const. In the current case, the event horizon does exist even before the injection starts and when the bulk spacetime is still locally pure AdS. The brane can enter the horizon within this early time domain. It means that the brane can intersect with the event horizon even at V < 0 (before the injection starts!). However, the brane configuration in that time domain is given by the exact solution in the pure AdS, and it will be obviously regarded as the Minkowski embedding if static.

JHEP04(2014)099
In the spirit of the gauge/gravity duality, to determine the phase of the brane embedding, we should look at the quark condensate c rather than the temporal configuration of the brane on a certain time slice. If we consider the causality in the bulk further, we notice that observers outside the black hole (including those on the AdS boundary) cannot know whether the brane intersects the horizon or not, since the intersection locus lies on the event horizon and the exterior observers can never obtain information on that locus. This implies that an attempt to determine the phase based on the temporal brane configuration requires information inaccessible from the boundary. On the other hand, the behavior of the quark condensate c is defined only by the information outside the horizon, since this quantity is governed by the brane equations of motion that respect the bulk causality. Indeed, using the quark condensate c, we can examine the phase of the brane embedding without any problem, as shown in figure 9. From the viewpoint of the boundary, the phase transition takes place when the brane is drawn into the vicinity of the black hole (namely, when the gravitational red-shift at the brane becomes strong) rather than the inside of that.

Overeager case
As a new significant feature in the dynamical evolution, we find that the brane can fall into the event horizon even in a parameter region where any equilibrium BH embedding does not exist. We would like to refer to this situation as the overeager case, as explained in section 1.3. Figure 11 shows the time evolution of quark condensate c and the black hole mass M in such an overeager case.
The time evolution of c(V ) markedly differs from the previous cases: there are neither periodic oscillations nor relaxation to the equilibrium value. Also, we show snapshots of the brane configuration at early time (V ≤ 0.6) and late time (V ≥ 0.8) in figures 12 and 13, respectively. We take r h = 1.06 as the final radius of the black hole, for which only the Minkowski embedding exists as the static brane configuration (see figure 4). However, the brane configurations in figures 12 and 13 indicate that the brane intersects with the event horizon. This fact is reflected in the non-periodic behavior of c. We also notice that the brane keeps moving even after the bulk spacetime has settled down to the final state.
We can understand this bulk phenomenon in the following way. It is not forbidden that a part of the brane dynamically intersects with the event horizon even if any static BH embedding does not exist. 13 Once some part of the brane oversteps the horizon, that part can never come out from the horizon by the definition of the event horizon. Since no equilibrium configuration as the BH embedding exists, the brane has no choice but to remain dynamical with intersecting the horizon.
It is worth noting that it depends on the energy injection duration ∆V whether this overeager case is realized or not. As we mentioned in section 4.1, the excitations on the brane become larger as the energy injection becomes faster. Therefore, it is expected that 13 In a special case, this statement can be easily proved as follows. When ∆V is infinitesimal, namely for the fastest injection, the radius of the event horizon becomes r = r h just at V = 0 and it will be the final radius. On the other hand, if the quark mass is m < r h , the tip of the static brane is located at r = m (< r h ) before the injection V < 0. This implies that the portion of the brane near the tip resides within the event horizon at V = 0. Since there is no equilibrium BH embedding for 1 < r h /m 1.08, the brane must be overeager in such cases. if the energy injection is sufficiently slow, the brane can not be overeager. In fact, we can see this expectation is true. Figure 14 shows time evolution of c in a case of slow injection. This reveals that the quark condensate c is periodically oscillating around the equilibrium value, and its behavior is different from the one in the overeager case shown previously. This result implies that the brane is the Minkowski embedding with periodic oscillations in the sub-critical cases. Thus, non-adiabaticity of the energy injection plays an important role in the overeager effect.
Are the mesons melting in the overeager case? As shown in figure 11, the quark condensate c(V ) shows a damped oscillation around the smooth time-evolving component, while the normal-mode oscillations present in the sub-critical case is absent in this case. Also, if we focus on the snapshots of the brane configuration on constant V surfaces, we find that the brane touches the bulk black hole for a time scale relatively longer than that of the thermalization (see figure 13). These facts suggest that the meson is melting in the  overeager case, though the temperature of the bulk black hole is lower than the critical one. We argue the dual field theory implications of this overeager case in section 5.1.

Final fate for the overeager case
What will be the final brane configuration of the time evolution in the overeager case? As mentioned in the previous subsection, there is no corresponding static BH embedding. As explained in appendix A, we can solve the evolution equations only outside of the event horizon and the red-shift prevents us from computing in the region very close to the horizon. For these reasons, we do not have whole numerical solutions of the brane configuration at late time. In this subsection, to infer a final fate of the overeager case, we extrapolate numerical solutions of the brane configuration using the Padé approximation, and estimate late-time evolutions of the brane. Let us define e ≡ (r min − r h )/r h where r min is the minimum value of r ≡ 1/z in our numerical data at each bulk time. This parameter e measures the distance of the extrapolation. We use numerical data satisfying e < 0.03 in this subsection. Figure 15 We set parameters as r h = 1.06 and ∆V = 0.5. Our numerical data are plotted by black dots and their inter/extrapolation are depicted by solid curves. We see that, at V = V crit ≃ 5.15, the brane position on the horizon reaches the pole, at which the radius of the S 3 wrapped by the brane becomes zero. In figure 15(b), we plot the radius of the S 3 on the horizon, R S 3 ≡ cos Φ| r=r h , as a function of V . In this figure, numerical data with e < 10 −2 and 10 −2 < e < 3 × 10 −2 are shown by black and white circles, respectively. From the figure, we can see that the radius of S 3 becomes small as time increases, and it seems to reach zero within finite time. Since the extrapolation of the brane shape employed here may cease to be a good approximation near the horizon and the pole when e becomes large, the precise motion of the brane at late time could be different from that described here in that case. However, these results with a help of the extrapolation imply that the brane would reach the pole in the neighborhood outside the horizon around V ∼ V crit and the brane embedding would be singular. 14 Classical description of the brane would break down at such a singular point and a singularity outside the horizon could be seen from the boundary.
If we take into account stringy effects, there is a possibility that the brane reconnects near the pole around V ∼ V crit . When the radius of the S 3 wrapped by the brane becomes small, we may expect that brane will reconnect. Let us estimate the condition for the 14 Originally, the tip of the brane is located at the pole and the embedding is regular there. If the embedding is regular at V = V crit , the tip of the brane must be just at the pole on the horizon. However, the tip of the brane has fallen into the black hole and a regular tip cannot be at the pole on the horizon. This implies that other parts of the brane have concentrated on the pole and the embedding would be singular.

JHEP04(2014)099
brane reconnection. The radius of the S 3 is given by R S 3 = L cos Φ| reconnect ≃ Lǫ where we defined ǫ = π/2 − Φ| reconnect . Note that we have restored the AdS curvature scale L, which is the radius of S 5 too. The condition for the reconnection is given by R S 3 ℓ s , and it can be rewritten as ǫ λ −1/4 , where λ is the 't Hooft coupling defined by λ = L 4 /ℓ 4 s . This inequality indicates that, if we take into account the stringy (finite λ) effect, emergence of the singularity may be avoided by the brane reconnection. In a similar manner, the reconnection may occur by taking bulk quantum effects into account via membrane nucleation. Once the reconnection occurs, the part of the brane connected to the boundary would be separated from the black hole and oscillate around the equilibrium configuration, which is the static Minkowski embedding determined by the bulk black hole mass.

Boundary theory interpretations
We discuss boundary theory interpretations and implications of our numerical results, especially focusing on the overeager case. In this case, the brane touches the horizon temporarily and tends to be singular in a finite time duration. If we take into account stringy or quantum effects in account, the brane may reconnect near the singular point, and go back to the Minkowski embedding with oscillations induced by the reconnection. We consider this process in the boundary theory point of view below.
In the overeager case, the meson is melting even though the final temperature of the gluon medium (represented by the thermal D3 background) is lower than the meson melting temperature. One possible explanation for this phenomenon would be as follows. A rapid energy injection drives the ambient gluon medium into a non-thermal state. In such a state, the distribution function of the gluons will be deformed from the thermal one governed by its temperature. If the deformed distribution function possesses sufficient amount of high-energy components, gluons would break bound states of quarks and melt the mesons. It would be interesting if we could examine such a scenario from the bulk point of view. For example, the spectrum of the bulk Hawking radiation corresponds to the distribution function of the gluon medium mentioned above, and can be studied using techniques of refs. [49,50]. Such a study would provide useful information for this issue.
At the late time of the overeager case, the gluon medium approaches thermal equilibrium and high-energy gluons decrease. Since the final temperature is lower than the meson melting temperature, some of the dissociated quarks and antiquarks will recombine into mesons giving off the latent heat into the mesons and the gluon medium. The recombination of mesons may correspond to the reconnection of the D7-brane in the gravity side.

Open questions and future directions
We have clarified some basic aspects of dynamics in the D3/D7 system with the aid of numerics in this paper, but we still have many open questions to be addressed. Below, we list them along with ideas for future directions of the study.

JHEP04(2014)099
• Applications to other systems In this paper, we have developed a numerical scheme that realizes a robust time evolution for a long duration to study the probe brane on the dynamical background spacetime. This technique and its generalizations may have various applications to study non-linear dynamics within the context of the gauge/gravity duality.
One possibility is to apply to more general brane systems. An example is the black hole-string system like those studied in ref. [51], which corresponds to the quarks moving in the thermal plasma medium. Our scheme will enable to consider the effect of the dynamical background and relaxation phenomena therein. It will be interesting also to consider brane systems that realize more realistic field theories, in particular those realizing confinement, such as the D4/D6 system [35] or the Sakai-Sugimoto model [52]. We did not include the world-volume field strength and external fields in this paper, and their effects deserve further investigation (see e.g. ref. [36] for the effect of external magnetic field).
Yet another future direction is to consider higher-cohomogeneity cases. For example, ref. [42] considered a similar melting process realized by a rapidly-moving point particle modeled by the Aichelberg-Sexl solution in the bulk. Starting from a system whose temperature is slightly below the melting temperature and using the nearhorizon approximation, they shown that the mesons melt in a spatially localized region whose size is larger than the scale naturally expected from the injected energy E and the temperature difference ∆T from the melting temperature by a factor of O((E/N 2 c ∆T ) 3/2 ). This behavior is qualitatively similar to that of the overeager case, in the sense that the it becomes easier to melt mesons in the non-equilibrium regime. They focused on the regime where ∆T /T ≪ 1, and then it will be interesting to see how this behavior changes in general cases. Also, it may deserve to study if the phase corresponding to the overeager case, in which the mesons melt at temperature lower than the critical one, persists in this local heating case.

• Brane motion and reconnection
We observed that, in the overeager case, the intersection locus of the brane and bulk horizon moves toward the pole, where a brane reconnection is expected to occur. It would be interesting to clarify the brane dynamics after the reconnection assuming an explicit reconnection process or taking into account effects beyond the DBI action. Dynamical processes involving a brane and a black hole have been studied in, e.g., refs. [53][54][55][56] in the context of braneworld models, and these previous results may give us implications to the physics after the brane reconnection or the quark recombination. Also, it would be important to understand influence of the effects beyond the DBI action, such as the curvature corrections arising from non-zero brane thickness [57][58][59][60][61].
One possible way to analyze the brane dynamics in a simplified manner is to focus on the near-horizon region. Ref. [62] studied static brane configurations in the nearhorizon region. They constructed a family of static solutions in that region and JHEP04(2014)099 clarified its basic properties. Also ref. [42] focused on the near-horizon region in the D3/D7 system, and clarified the brane motion just after the local heating as explained above. It may be interesting to study the brane dynamics in this setup, especially focusing on the regime where the reconnection occurs.

• Brane induced geometry and fluctuations
We may read out some extra pieces of information from the induced geometry on the D7-brane. Fluctuations of the brane is governed by the induced geometry, since the intrinsic curvature gives rise an effective potential for those fluctuations. By studying the relationship between the intrinsic geometry and the (quasi-)normal mode frequency of the brane fluctuations, we may achieve a better understanding on the time dependence of c(V ) from the bulk point of view.
In this context, it is interesting to focus on the overeager case, where the brane turns from the BH embedding into the Minkowski embedding at some moment. A naive expectation is that the decay rate of the quasi-normal modes becomes smaller as the transition time approaches, since it becomes zero after the transition. The intrinsic curvature tends to be singular in this process, and it may cause such a time dependence through the effective potential for the fluctuations. It will be interesting to test this conjecture by comparing the time dependence of the brane intrinsic curvature and the brane fluctuations. Also, this transition can be viewed as a black hole formation and subsequent evaporation in the brane induced geometry. It may deserve a further consideration if we can relate such an idea with the time dependence of the brane fluctuations and c(V ). For this purpose, it would be fruitful to study more about the condensate fluctuation and its relaxation time scale for the overeager phase and also for the BH phase, as refs. [63,64] in the context of holographic superconductors.
In figure 6, we observed that the discrete spectrum of c(V ) in the Minkowski embedding shows a exponential decay at higher frequency region. We may be able to read out a characteristic scale from this exponential decay. It might be interesting to elaborate this idea further, though this kind of exponential decay in the spectrum could be simply due to basic properties of the decomposition of the initial perturbations into the normal modes at high frequency band (see e.g. ref. [65] for behavior of Fourier coefficients in the high-frequency limit).

JHEP04(2014)099
of the National Strategic Reference Framework (NSRF) under "Funding of proposals that have received a positive evaluation in the 3rd and 4th Call of ERC Grant Schemes". SK is partially supported by the JSPS Strategic Young Researcher Overseas Visits Program for Accelerating Brain Circulation "Deeping and Evolution of Mathematics and Physics, Building of International Network Hub based on OCAMI". The work of NT is supported in part by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, and JSPS Grant-in-Aid for Scientific Research 25·755.

A Numerical methods
In this appendix, we explain our numerical methods to solve the time evolution of the D7brane. We need two kinds of numerical schemes depending on whether the brane intersects with the event horizon or not. Before and after the brane intersects with the event horizon, we use numerical method in subsection A.1 and A.2, respectively.

A.1 Before the brane intersects with the event horizon
In this section, we explain our numerical method before the brane intersects with the event horizon. 15

A.1.1 Evolution scheme
The time evolution of the D7-brane is described by non-linear wave equations (3.13)-(3.15). We will use the finite-difference method for solving them numerically. It turns out that the time evolution based on these equations is numerically unstable. To stabilize it, we eliminate V ,u and V ,v in the right hand side of the equations using the constraints (3.16) and (3.17) as where each sign in the front of the square root is determined by the boundary conditions at Z = 0, namely V ,u = 0 and Z ,u > 0, or V ,v > 0 and Z ,v < 0. As a result, the evolution equations can be schematically rewritten as where X = (V, Z, W ) andX = (Z, W ). We can realize a stable numerical time evolution if we use these equations as the evolution equations. 16 As shown in figure 16, we discretize u and v coordinates with the grid spacing h. Now, we focus on the points N, E, W, S 15 If the brane intersects with the event horizon at the initial surface, which occurs when the initial background temperature is sufficiently large, we use the method in section A.2 from the beginning. 16 We found the method to stabilize the numerical calculation by trial and error, and the reason of the stabilization is unclear at this moment. It will be useful to clarify the stabilization mechanism, especially when we try to apply this numerical method to more generalized systems. One way to analyze it would be to show the well-posedness of the evolution equations, employing arguments similar to those in ref. [66][67][68][69].  and C in the figure. We can evaluate the function X and its derivatives at point C with second-order accuracy as where the index of X denotes each point N, E, W or S at which X is evaluated. Substituting the above expressions into eq. (A.2), we obtain nonlinear equations for X N . Solving the nonlinear equations by the Newton-Raphson method or the predictor-corrector method, we can determine unknown value of X N from X E , X W and X S which have been already given.

A.1.2 Boundary conditions at the AdS boundary
One of the boundaries in our numerical domain is the AdS boundary, Z = 0. As we explained in section 3.3, we fix the location of the AdS boundary at u = v by using a residual coordinate freedom throughout our calculations. Suppose that the point C is located on u = v. The boundary conditions for Z and Ψ leads to Z N = 0 and Ψ N = 1 (where the quark mass m is unity), shortly. For V , the boundary conditions (3.21) are discretized as

JHEP04(2014)099
where the point E is a dummy grid point introduced in order to impose the boundary conditions. From Z N = Z S = 0 on the boundary, we obtain This equation determines the time evolution of V at the AdS boundary.

A.1.3 Boundary conditions at the pole
The other boundary is the pole, Φ = π/2, at which the radius of S 3 wrapped by the D7brane shrinks to zero. We set the location of the pole at u = v + π/2. For Z and V , we have the Neumann boundary conditions as at u = v + π/2. (If we introduce Cartesian coordinates σ = u − v and τ = u + v, they are rewritten as Z ,σ | u=v+π/2 = V ,σ | u=v+π/2 = 0.) Now, we focus on the points P, Q and R in figure 16. Since V and Z satisfy the Neumann boundary conditions, we can approximate them by quadratic functions as V = V P + a(u − v − π/2) 2 and Z = Z P + b(u − v − π/2) 2 near the pole. We can determine V P , Z P , a and b from values of V and Z at the points Q and R. As a result, we obtain These expressions are essentially the backward differences of V and Z at the point P of second-order accuracy in the spatial direction. For Ψ, we have Ψ P = π(2Z P ) −1 as the boundary condition. Note that we need to know the values of Z, V and Ψ at the points Q and R prior to imposing these boundary conditions the point P. It is possible if we solve the evolution equations initially along v-direction and then proceed to the next step along u-direction.

A.1.4 Initial data and time evolution
We impose initial data as eq. (3.24) on the initial surface v = 0. Since we have located the AdS boundary and the pole at u = v and u = v + π/2, respectively, the function φ(u), which is the residual coordinate freedom, should be φ : [0, π/2] → [0, π/2]. For simplicity, we will choose φ(u) = u here. Now, we can solve the time evolution in the current numerical scheme until the brane intersects with the event horizon. However, if the brane has intersected with the event horizon, V will blow up rapidly near the horizon along v-direction and our numerical calculation will be broken immediately after the intersection. Therefore, we should monitor values of Z and V at the pole for each v-slice. If the pole is inside the horizon, Z| u=v+π/2 > Z EH (V | u=v+π/2 ), we pause the numerical integration at a intermediate surface v = v int and switch to the other scheme explained in the following subsection. Then, we will restart to integrate taking the numerical solution on the intermediate surface as the initial data for the new scheme. We denote the solutions on the intermediate surface as V int (u) = V (u, v int ), Z int (u) = Z(u, v int ) and Ψ int (u) = Ψ(u, v int ). On the intermediate

JHEP04(2014)099
surface, the u-coordinate is defined in v int ≤ u ≤ v int + π/2. The lower and upper bounds of the inequality correspond to the AdS boundary and the pole, respectively. Since the solutions on the intermediate surface are obtained just on grid points, we interpolate them by polynomials and regard them as functions in u ∈ [v int , v int + π/2] to construct the initial data for the new scheme.

A.2 Time evolution after the brane intersects with the event horizon
Here, we explain the numerical method after the brane intersects with the event horizon. The evolution scheme and the boundary conditions at the AdS boundary are the same as shown in section A.1.1 and section A.1.2. At the pole, we do not need to impose any boundary conditions because it is already inside the event horizon and is no longer contained in the current numerical domain.
As the initial data, we use the functions V int (u), Z int (u) and Ψ int (u), obtained in the last subsection. If the brane has intersected with the event horizon at an initial surface or if the brane will intersect shortly after the energy injection, we use eq. (3.24) as the initial data from the beginning. As we mentioned before, V at the AdS boundary blows up if the brane is close to the event horizon. To avoid it, we will adopt a coordinate system to synchronize, at the AdS boundary, the coordinate time on the world-volume with the asymptotic time V using the residual coordinate freedom at the initial surface. We consider the coordinate transformation on the initial surface as u = φ(ū). Then, the v-coordinate is also transformed as v = φ(v) to locate the AdS boundary atū =v. We choose this free function φ so that V andv are synchronized up to a constant at the AdS boundary, i.e., V | u=v =v + V 0 .
Hereafter, we explain how to determine the φ which synchronizes V | u=v andv. As in figure 17, we take the domain ofū-coordinate asū ≥ 0 and discretize it asū i = ih (i = 0, 1, 2, . . .). We denote φ(ū i ) = φ i and set φ 0 = v fin . We also denote boundary values of V at each grid point as V 0 , V 1 , V 2 , . . .. We determine the numerical solution in the order of P 1 → P 2 → P 3 → · · · . (Roughly speaking, we regard u as the "time" coordinate.) Firstly, we solve the evolution equations at P 1 using a trial value φ trial 1 . 17 Then, at the point P 1 , we obtain an asymptotic time V trial 1 at the AdS boundary. Since we would like to synchronize V andv at the AdS boundary (namely, V 1 − V 0 = h), the appropriate choice of φ 1 should be For this choice of φ 1 , we recalculate the evolutions and obtain Next, we solve the evolution equations at P 2 and P 3 with a trial value φ trial 2 and obtain an asymptotic time V trial 2 at the AdS boundary. Since we would like to impose V 2 − V 1 = h, we should choose φ 2 as 17 In our calculation, we used φ trial  Then, we recalculate them and obtain V 2 = V 1 + h + O h 2 . By repeating this procedure, we can synchronize V andv at the AdS boundary. Note that, in this choice of free function φ, our numerical domain does not touch the event horizon since, along aū-constant surface, the function V approaches a finite value at the AdS boundary and the surface should be outside the horizon. (Intuitively speaking, observers at the AdS boundary cannot see the event horizon within finite asymptotic time.) In principle, we can solve the time evolution forever (V | u=v → ∞) by the current scheme in contrast to the previous scheme in section A.1. However, especially in the overeager case, numerical error accumulate in actual calculations and it becomes an obstacle for an numerical evolution over a long time duration. In a typical set up, our numerical calculation is reliable in V | u=v − V 0 10 at least, which is sufficient for our purpose.

B Error analysis
We test our numerical code by (i) checking if a static initial data remains static in the time evolution by our code, and (ii) showing second-order convergence of the constraints. Our code passes these tests as shown below.

B.1.1 Vacuum bulk
When the bulk spacetime is the pure AdS spacetime without a black hole, the static initial data is exactly given by eq. (3.24). We check if our time evolution code keeps it to be static and the quark-condensate c(V ) to be zero.
We show |c(v)| in this case in figure 18(a) for grid number N = 100, 400, 1600. For each N , |c(v)| is approximately given by an constant component with periodic pulse noise, whose detailed structure is shown in figure 18(b). These sharp peaks are generated by time evolution near the pole, where the evolution equations become singular.  Figure 18. |c(v)| for the exact static initial data (3.24) in the pure AdS bulk spacetime.c is the average value of |c(v)| obtained by numerically fitting |c(v)| to a constant over 0 < v < 10. We measure the peak height from the average by max 0<v<10 ( |c(v)| −c ).  We show the N dependence of the constant component and the peak height around it in figure 18(c), where we defined the average value of |c(v)|,c, by numerically fitting to a constant over 0 < v < 10, and the peak height by max 0<v<10 |c(v)| −c . Fitting the numerical data, we find that both the constant component and the peak height around it converge to zero as we increase N by power-law with powers p = −1.94 and p = −1.43, respectively.

B.1.2 Bulk with a black hole
We also test if static initial data remains static in the time evolution given by our code even when a bulk black hole exists. For this purpose, we construct static initial data by integrating eq. (2.4) numerically, and studied time dependence of c(v) provided by that.
We show c(v) for the static Minkowski embedding by eq. (2.4) for the black hole mass M = 0.500 (the horizon radius r h = 0.841) in figure 19(a). Similarly to the previous case, c(v) is approximately composed of constant component and periodic pulse noise around it (see figure 19(b)). We calculate the average value of the quark condensatec by fitting c(v) to a constant over 0 < v < 10 for each N . This average value depends on N , and its asymptotic value, c ∞ , is estimated by fittingc(N ) = c ∞ + aN p numerically. We

B.2 Constraint violation
In this section, we check that the constraint violation is maintained to be sufficiently small in the dynamical calculations shown in section 4. We show that the constraint violation is sufficiently small throughout the time domain discussed therein, and also that it shows second-order convergence when the resolution is raised.

B.2.2 Dynamical black hole embedding
Below, we show the constraint violation for the dynamical black hole embeddings discussed in sections 4.2 and 4.3.    Figure 22(a) showsC 1 at V = 0.5, 5.5, 6.0, 7.0, where Φ on the horizon approaches π/2 at V ≃ 5. 5 We can see thatC 1 grows to O(10 −4 ) at Φ ≃ 0.97 × π/2 for N = 100. Figure 22 figure 23(c). Fitting shows that they converge to zero by power law with p = −1.98 and −2.02, respectively. It implies that the constraint violation up to V ≤ 8 shows second-order convergence to zero for higher resolution.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.