Turbulent strings in AdS/CFT

We study nonlinear dynamics of the flux tube between an external quark-antiquark pair in $\mathcal{N}=4$ super Yang-Mills theory using the AdS/CFT duality. In the gravity side, the flux tube is realized by a fundamental string whose endpoints are attached to the AdS boundary. We perturb the endpoints in various ways and numerically compute the time evolution of the nonlinearly oscillating string. As a result, cusps can form on the string, accompanied by weak turbulence and power law behavior in the energy spectrum. When cusps traveling on the string reach the boundary, we observe the divergence of the force between the quark and antiquark. Minimal amplitude of the perturbation below which cusps do not form is also investigated. No cusp formation is found when the string moves in all four AdS space directions, and in this case an inverse energy cascade follows a direct cascade.


Introduction
The gauge/gravity duality [1][2][3] is successfully applied to investigating strongly coupled gauge theories. Through this duality, it is hoped that one can access their nontrivial aspects that are hard to be handled because of the strong coupling. In particular, among advantages in using gravity duals, it is worth notifying that dynamics in time-dependent systems can be powerfully computed from time evolution in classical gravity. Applications range over far-from-equilibrium dynamics governed by nonlinear equations, and using numerical techniques for solving them attracts much attention. For instance, physics of strongly coupled plasma of quarks and gluons at RHIC and LHC brings motivations to numerically study dual gravitational dynamics; a series of seminal works is in Refs. [4][5][6][7][8].
Far-from-equilibrium processes in the D3/D7 brane system dual to N = 2 supersymmetric QCD have been recently studied in Refs. [9][10][11][12][13], where partial differential equations for time evolution were solved numerically. As a phenomena characteristic in non-linear dynamics, it has been found that long time evolution of the D7-brane generates a singularity on the brane, and this is understood from the viewpoint of weak turbulence on the D7-brane: The energy in the spectrum is transferred from large to small scales [10][11][12]. Small scale fluctuations there correspond to excited heavy mesons in the dual gauge theory. The turbulent behavior of the D7-brane can be interpreted as production of many heavy mesons in the dual gauge theory, and the singularity formation is interpreted as deconfinement of such mesons.
To gain a deep insight into this kind of nonlinear dynamics in the gauge/gravity duality, in this paper, we will consider a string in AdS dual to the flux tube between a quark-antiquark pair in N = 4 super Yang-Mills theory, and fully solve its nonlinear time evolution with a help of numerical techniques. This setup corresponds to focusing on a Yang-Mills flux tube compared with the collective mesons described in the D3/D7 system. In addition, working in this setup is simpler than using the D3/D7 system and will provide a clear understanding of the turbulent phenomena and instabilities in probe branes in the gauge/gravity duality. To initiate the time evolution, we will perturb the endpoints of the string for an instant. Our setup is schematically illustrated in Fig. 1. The endpoints are forced to move momentarily and then brought back to the original locations. We will loosely use the terminology "quench" for expressing this process. This action introduces waves propagating on the string, and we will be interested in their long time behavior where nonlinearity in the time evolution of the string plays an important role. When we perform numerical computations, we will use the method developed in [9], which turned out to be efficient for solving the time evolution in probe brane systems. We are also motivated by the weakly turbulent instability found by Bizoń and Rostworowski [14]. Our setup may give a simple playground to study that kind of phenomenon.
The string hanging from the AdS boundary is one of the most typical probes AdS boundary in the gravity dual. This gives the gravity dual description of a Wilson loop corresponding to the potential between a quark and an antiquark [15,16]. Although the potential in a conformal theory is different from that in real QCD, linear confinement can be realized in nonconformal generalization [17]. In finite temperatures, a string extending to the AdS black hole corresponds to a deconfined quark [18,19], and this has been utilized for studying the behavior of moving quarks in Yang-Mills plasma [20][21][22][23][24][25][26]. Moving quark-antiquark pairs were also considered [27,28]. In farfrom-equilibrium systems, holographic Wilson loops have also been used as probes for thermalization [29,30]. Nevertheless, veiled by these applications to QGP, nonlinear (and non-dissipative) dynamics of the probe hanging string in AdS has not been shed light on so much. Some analytic solutions of non-linear waves on an extremal surface in AdS have been studied in Ref. [31]. Notice that that configuration corresponds to a straight string in the Poincare coordinates.
When it comes to nonlinear dynamics of a string, formation of cusps would be primarily thought of. In fact, it is well known in flat space in the context of closed cosmic strings that cusp formation is ubiquitous [32]. Cusp formation of fundamental strings ending on D-branes has been also found in Ref. [33]. We will turn our attention to whether there is such formation of cusps also in AdS. 1 The organization of the rest of this paper is as follows. We start from reviewing the static solution and linear perturbations in Section 2 and 3, where we introduce a parametrization convenient for our use. In Section 4, we explain the setup for our time dependent computations. We introduce four patterns of quenches that we consider, derive the evolution equations and the boundary conditions, prepare initial data, and explain measures for evaluating the time evolution. Sections 5, 6, and 7 are reserved for numerical results: In Section 5, we discuss two of the four quenches where the oscillations of the boundary flux tube are restricted to compression waves and therefore we call them longitudinal. We evaluate cusp formation, turbulent behavior in the energy spectrum, and the forces acting on the quark endpoints. In Section 6, we show results of a quench where the flux tube oscillates in one of its transverse directions. Finally, in Section 7, we examine the last quench where the motion of the flux tube is in all three spatial directions of the boundary 3+1 dimensions. Section 8 is devoted to summary and discussion. In appendices, we explain details for the numerical computations and derive some formulae used in the main text.

A review of the static solution
We briefly review the holographic calculation of the static quark-antiquark potential in the near horizon limit of extremal D3-branes [15,16]. The background metric is AdS 5 × S 5 , We consider a rectangular Wilson loop with the quark-antiquark separation L along the x 1 -direction where the quark and antiquark are located at x 1 = ±L/2. The dynamics of the string is described by Nambu-Goto action, where γ = det(γ ab ) and γ ab is the induced metric on the string. It is convenient to take a static gauge where the worldsheet coordinates (τ, σ) coincide with target space coordinates as (τ, σ) = (t, z). The static solution is then described by a single function x 1 = X 1 (z). 2 The Nambu-Goto action becomes where λ = 4πg s N c is the 't Hooft coupling. Solving the equation of motion of X 1 (z) gives the bulk string configuration. Let z = z 0 specify the bulk bottom of the string reached at x = 0, where the regular boundary condition ∂ x z = 0 is imposed. The embedding solution is given by where Γ 0 ≡ √ 2π 3/2 /Γ(1/4) 2 ≃ 0.599, and F and E are the incomplete elliptic integrals of the first and second kinds defined as Setting z = 0 in (2.4), we obtain L/2 = z 0 Γ 0 . In Fig. 2, we show the profile of the static string in the (z, x 1 )-plane. We will use this solution as the initial configuration for our time evolution. It is also known that the dependence of the potential energy on L is Coulomb. The energy is evaluated from the on-shell action, which in general diverges at the boundary, but this divergence can be regulated by comparing with the diverging energy of two strings straightly extending to the Poincare horizon. The regularized energy is then given by In this sense, the quark-antiquark potential in the AdS background does not correspond to the confining potential of real QCD. This, however, is considered as a handy playground for testing nonlinear evolution in the gauge/gravity duality.

Linear perturbation theory
Linearized fluctuations and stability of holographic quark-antiquark potentials have been studied in Refs. [35][36][37][38]. Here, we solve the linearized fluctuations of the static solution in coordinates convenient for our use. In (2.4), we parametrized the location of the string in terms of the z-coordinate. In this coordinate, however, the static solution X 1 (z) becomes multi-valued, and this may not be suitable for considering  linear perturbations. Instead, we introduce new polar-like coordinates (r, φ) in which the static embedding is expressed by single-valued functions, where we define the functions f and g as where sn(x; k) is a Jacobi elliptic function defined as the inverse function of F (x; k) given in Eq. (2.5): F (sn(x; k); k) = x. For k = i = √ −1, the Jacobi elliptic function has roots at x = β 0 n (n ∈ Z), where β 0 = π/(2Γ 0 ) ≃ 2.622. We find that there is a nice relation between f and g, With these functions f and g, the static embedding (2.4) is simply given by r = z 0 . The profiles of the functions f (φ) and g(φ) are shown in Fig. 3(a). We also depict the r = const. and φ = const. surfaces in the (z, x 1 )-plane in Fig. 3(b). Using t and φ as the worldsheet coordinates, we can describe the dynamics of the string in terms of three functions, We consider perturbations around the static solution, R = z 0 , X 2 = X 3 = 0, as  where χ i (i = 1, 2, 3) are dimensionless perturbation variables. We will refer to χ 1 as the longitudinal mode and χ i (i = 2, 3) as the transverse modes. Then, in the second order in χ, the Nambu-Goto action becomes where we define˙≡ ∂ t and ′ ≡ ∂ φ , and introduce h(φ) ≡ [(g/f ) ′ f ] 2 . To derive the above expression, we used the relation of f and g (3.4) and omitted the total derivative terms. The equations of motion for (χ 1 , χ 2 , χ 3 ) are (3.8) Operators H and H ′ are Hermitian under the inner products respectively. We denote the eigenvalues and eigenfunctions of H and H ′ as {ω 2 n , e n (φ)} and {ω ′ n 2 , e ′ n (φ)}, respectively. These are labeled by integers n ≥ 1 in ascending order of ω n and ω ′ n . The eigenfunctions are orthonormalized as (e n , e m ) = δ nm and (e ′ n , e ′ m ) ′ = δ nm . It is easy to check the linear stability of the static embedding against the longitudinal perturbation as In the same way, the stability against the transverse perturbations can be also checked, ω ′ n 2 ≥ 0. We numerically determine normal mode frequencies and eigenfunctions. These are plotted in Fig. 4. The spectra are asymptotically linear: ω n , ω ′ n ∝ n for n → ∞. In fact, from the WKB analysis, we can obtain z 0 ω n ≃ 2Γ 0 (n+1) for n → ∞ [35,36]. Our numerical results are consistent with the WKB approximation.

Non-linear dynamics of fundamental strings
The main focus of this paper is to study nonlinear dynamics of the string, where we make use of numerical techniques for solving the time evolution. In this section, the setup for this is prepared.

Setup
We consider AdS 5 × S 5 (2.1) as the background spacetime, and take the static solution (2.4) as the initial configuration. We then consider "quench" on the endpoints of the string: We move their positions momentarily and put them back to the original positions. A schematic picture of this setup is depicted in Fig. 1.
Let us denote the two endpoints of the string as x q (t) and xq(t), corresponding to the locations of the quark and antiquark, respectively. In this paper, we consider the following four kinds of quenches on x q (t) and xq(t): (i) Longitudinal one-sided quench: where α(t) is a compactly supported C ∞ function defined by The profile of this function is shown in Fig. 5. The flux tube vibrates in its longitudinal direction, and motions are not induced in the transverse directions by this quench. Thus, in this case, the motion of the string is restricted in (2+1)-dimensions spanned by (t, z, x 1 ).
(ii) Longitudinal Z 2 -symmetric quench:  (iii) Transverse linear quench: We shake one of the endpoints along x 2 -direction. String fluctuations in this direction are induced by the quench but not in x 3 -direction. Thus, the string oscillates in (3+1)-dimensions spanned by (t, z, x 1 , x 2 ).
In Fig. 6, we show schematic pictures of the quenches (i)-(iv). These patterns are chosen to represent typical string motions, particularly with different dimensionality. String dynamics is specified by two parameters ǫ and ∆t once we choose a quench type.
In this paper, we will take modest values for the amplitude of the quench (ǫ ∼ 0.01) since we are interested in nonlinear evolution starting from small deviation from the linear theory and focus on weak turbulence and cusp formation driven by the nonlinearity. For a large value of ǫ, we expect some effect similar to the overeager effect found in Ref. [9]: The string will be able to plunge into the Poincare horizon z = ∞ because of the strong perturbation. Such strong quenches will be studied in detail elsewhere [39].

Basic equations
To calculate the time evolution on the string worldsheet, we find it convenient to use double null coordinates. With worldsheet coordinates (u, v), the string position is parametrized as (4.6) Using these expressions with Eq. (2.1), we obtain the induced metric as (4.7) The reparametrization freedom of the worldsheet coordinates allows us to impose the double null condition on the induced metric as Notice that these conditions do not fix the coordinates completely: There are residual coordinate freedoms, These will be fixed by boundary conditions and initial data.
In the double null coordinates, the Nambu-Goto action (2.2) becomes where at the second equality we eliminate the square root in the action using the double null conditions (4.8). (Note that γ uv is negative.) From the action, we obtain the evolution equations of the string as (4.11) Using these, we find that the constraints (4.8) are preserved under the time evolution: Hence if we impose C 1 = C 2 = 0 on the initial surface and the boundaries, the constraints (4.8) are automatically satisfied in the whole computational domain. Our numerical method for solving the evolution equations is briefly summarized in Appendix A and its numerical error is estimated in Appendix B. For more details of the numerical method, see also Appendix A in [9].

Boundary conditions
On the worldsheet, there are two time-like boundaries that correspond to the two endpoints of the string attaching on the AdS boundary. We need boundary conditions there. Using the residual coordinate freedoms (4.9), we can fix the locations of the boundaries to u = v and u = v + β 0 . The boundary conditions for the spatial parts of the target space coordinates are given by We also need boundary conditions for T . To derive them, we solve Eq. (4.11) and (4.8) near the boundaries.
, we obtain the asymptotic solutions around σ = 0 as (4.14) where˙= d/dτ , v =ẋ q /ṫ 0 , a =v/ṫ 0 , j =ȧ/ṫ 0 and γ = 1/ √ 1 − v 2 . The asymptotic solutions near σ = β 0 can be given by replacing σ → β 0 −σ and x q → xq in the above expressions. From the second equation in Eq. (4.14), we have ∂ σ Z| σ=0 = ṫ 2 0 −ẋ 2 q . Using this, we obtain These equations determine the time evolution of T at the boundaries σ = 0, β 0 . Their numerical implementation is explained in Appendix A. For consistency, the speed of the quark endpoints during the quench should be slower than light, |v| < 1. Otherwise, the Lorentz factor γ becomes imaginary. Solving this condition, we obtain constraints for the quenches (i)-(iii), and for the transverse circular quench (iv), The parameter values examined in this paper satisfy these conditions.

Initial data
Before the quenches are applied, we assume that the string is static, namely, we use the static solution (2.4) as the initial configuration. For numerical computations, we need initial data written in the (u, v)-coordinates. Substituting Eq. (2.4) into Eqs. (4.8) and (4.11), we can express the static solution in terms of (u, v) as where φ 1 and φ 2 are arbitrary functions associated with the residual coordinate freedom (4.9), and the functions f and g are defined in Eqs. (3.2) and (3.3).
Locating the initial surface at v = 0 on the worldsheet, we parametrize the initial configuration as where we set the free functions φ 1 (u) = u and φ 2 (0) = 0 so that the boundaries are at u = 0, β 0 . If the string endpoints are not perturbed, our numerical calculations describe the static evolution of the exact solution (4.18).

Quantities for evaluation
In the non-linear dynamics of the string, we expect to observe formation of cusps. For detailed analyses, we will in particular use the following quantities for evaluation.

Cusp formation
Let us consider the profile of the string on a surface with t =constant. There is a cusp if the string does not change its target space position when the parameter on the string is varied. Using the u-coordinate as a parameter on the string in the t =constant surface, we obtain the conditions for the cusp formation as ∂ u X I | t=const = 0 where X I ≡ (Z, X i ) and i = 1, 2, 3. These conditions are rewritten as In our numerical computations, we monitor the roots of J I , which typically appear as curves on the (u, v)-plane. If these curves overlap at a point, we find cusp formation, and if the overlap continues in the time evolution, it is implied that the cusps continue to exist.
As an obvious corollary, if (4.20) is satisfied, there is a necessary condition that is also satisfied. This condition is conveniently utilized for a consistency check of the cusp formation detected by (4.20).

Energy spectrum in the non-linear theory
We will also study the energy spectrum of the non-linear fluctuations of the string because, from the spectrum, it is expected to find weak turbulence on the string. Once a dynamical solution (T (u, v), Z(u, v), X(u, v)) is calculated, we can convert it to the polar-like coordinates introduced in Eq. (3.1): , we can express the dynamical solution using target space coordinates (t, φ) as r = R(t, φ), As in Eq. (3.6), we define the non-linear version of the "perturbation" variableŝ χ 1 ,χ 2 ,χ 3 asχ We then decomposeχ 1 ,χ 2 andχ 3 with the eigenfunctions of the linear theory e n (φ) and e ′ n (φ), which were introduced below Eq. (3.8), aŝ Using the mode coefficients c n and c i n , we define the energy contribution from the n-th mode ε n (t) and the total energy in terms of the linear theory ε(t) as The quantify ε n is conserved in linear theory. Therefore, if we find time dependence in ε n , it is a fully non-linear effect. Note that the total energy ε defined in the linear theory is also time dependent in the non-linear theory although the time dependence is suppressed by the amplitude of the quench:ε/ε = O(ǫ). Since we consider only small ǫ in this paper (ǫ ∼ 0.01), we do not put emphasis on the time dependence of the total energy. In our actual numerical calculations, we take a cutoff at n = 50 for evaluating ε in Eq. (4.25). Its cutoff dependence is also not essential for our following arguments on the energy spectrum on each t-slice.

Forces acting on the heavy quarks
From the dynamical solution of the string, we can read off the time dependence of the forces acting on the quark and antiquark located at x = x q and xq. The force acting on the quark can be evaluated from the on-shell Nambu-Goto action as The same formula can be applied to the force acting on the antiquark F (t) by replacing x q → xq. The derivation of this expression is summarized in Appendix D. For our numerical analysis, it is convenient to rewrite this expression in terms of the (τ, σ)-coordinates. Using the asymptotic expansions (4.14), we obtain We need to extract the third order coefficient x 3 from our numerical data. For this purpose, it is convenient to define Using Eq. (4.14), we find that this function behaves as

Results for the longitudinal quenches
From this section to Section 7, we discuss results of numerical computations for our quenches (i)-(iv). In this section, we firstly treat the longitudinal one-sided quench (i) and then discuss the longitudinal Z 2 -symmetric quench (ii).

Cusp formation
We start from the longitudinal one-sided quench (4.1). In Fig. 7, we show snapshots of the time evolution of the string for ∆t/L = 2 and ǫ = 0.03. The left panel is just after the quench, t/L = 2, 2.2, 2.4. We observe that perturbations are induced on the string by the quench. Late time behavior is shown in the right panel for t/L = 7, 7.2, 7.4, where it is seen that cusps are formed on the string. The time for the cusp formation is evaluated as t/L ∼ 5 by using (4.20), while the snapshots are taken when it is easy to confirm the cusps visually. Although the solution plotted in the target space (z, x 1 )-plane looks singular on top of the cusps, the fields T (u, v), Z(u, v) and X 1 (u, v) are indeed smooth on the (u, v)-plane. Therefore, even after the formation of the cusps, the time evolution can continue without a breakdown of numerical computations. Physically, however, finite-N c effects can be important at cusp singularities, and the time evolution after the cusp formation may be considered as unphysical without such corrections. We will discuss possible finite-N c effects at cusps in section 8. For the present, the dynamics of cusps are discussed without taking into account these corrections. We find that the cusps appear as a pair. (The cusps at t/L = 7.0 are magnified in the inset of the right panel.) The propagating speeds of the two cusps are different, and they continue to propagate on the string. Following this example, we survey cusp formation by changing quench parameters using the method described in Section 4.5.1. Results are shown in Fig. 8 when ∆t/L = 2 is fixed and ǫ is varied. In the left panel, the times for the cusp formation for each ǫ are plotted, and the corresponding locations in the φ-coordinate are shown in the right panel. The time for the cusp formation becomes longer as the quench amplitude becomes smaller. It is also seen that the cusp formation times are mildly discretized, as well as the locations of the formation away from the boundary. This discretization indicates that the cusps might be difficult to form near the boundary.
The times and locations of the cusp formation tend to degenerate in the very early time before the first reflection of the initial perturbation wave at X = −L/2 around t/L ∼ 3.5. Features in this region seem to be influenced by the largeness of the quench and would be different from late time dynamics. Quenches with large amplitudes will be reported elsewhere [39].
The cusp formation time becomes longer as the amplitude gets smaller. A natural question then is if the cusps can be formed for any small amplitude ǫ. As ǫ decreases, however, the variations in the fields are tinier and tinier, and eventually it becomes numerically difficult to use the cusp condition (4.20) for finding the cusp formation. To obtain a reasonable inference, we fit results and extrapolate to smaller ǫ. A result is given in Fig. 9, where we fit data points 4 at ǫ ≥ 0.01 and t cusp /L ≥ 10 by a   polynomial 5 a ǫ + b ǫ 1/2 + c. From the extrapolation, we find that there is a critical value ǫ crit below which the cusp formation time would be infinity. In Fig. 9, we obtain ǫ crit ≃ 7.5 × 10 −3 . We repeat this procedure to estimate the critical ǫ for different ∆t/L and see how it changes. For each ∆t/L, we fit data by a ǫ + b ǫ 1/2 + c and extrapolate to the limit of infinite cusp formation time to read off the value of critical ǫ. Results are shown in Fig. 10. We find that the critical value scales as (∆t/L) 3 . A fit of our results is ǫ crit = 9.3 × 10 −4 (∆t/L) 3 . Note that this scaling may be altered if ∆t/L becomes very long and the wavelength of the induced wave is comparable with the consistent with the fit.  length of the hanging string. The critical amplitude in such a case, however, will be also big. For this reason, we do not focus on larger ∆t/L in this paper.
The cusp formation found here is similar to the weakly turbulent instability in the global AdS [14]: Small fluctuations in that AdS propagate between the boundary and the center and, eventually, the perturbations collapse into a black hole after several bounces. A difference between our cusp formation and that AdS instability is their critical amplitude of initial perturbations. The AdS instability occurs for arbitrary small initial perturbations while our cusps are formed only for ǫ > ǫ crit > 0. By a perturbative study in Ref. [14], it has been suggested that a commensurable spectrum is a necessary condition for the collapse by arbitrary small perturbations. As in Fig. 4, the linear spectrum of the string, ω n and ω ′ n , are commensurable only for n → ∞. Because of the non-commensurable spectrum of the string, we need finite perturbations for the cusp formation.

Energy spectrum in the non-linear theory
Given the formation of cusps, we look into the time dependence in the energy spectrum following the procedure described in Section 4.5.2. For samples, we focus on the following four cases: (a) ∆t/L = 2, ǫ = 0.005, (b) ∆t/L = 2, ǫ = 0.01, (c) ∆t/L = 2, ǫ = 0.03, and (d) ∆t/L = 4, ǫ = 0.07. With the parameters (a), cusps do not form on the string, while cusps are created for (b), (c) and (d). In Fig. 11, we show the time dependence of the energy spectra for these parameters. The dashed curves are the energy spectra computed in the linear theory; see Appendix C for the calculations. Although the spectra are defined only for integer n, for visibility of the plot these results are generalized to continuous n.
High frequencies in the energy spectra are suppressed when the cusps are not formed. In Fig. 11(a), the spectrum can be well approximated by the linear theory just after the quench (t/L = 2). Although it slightly deviates from the linear theory as the time increases, we do not find any remarkable change in the late time.
In contrast, the energy spectra show power law behaviors in the cases of cusp formation. In Fig. 11(b), although the spectrum can be well approximated by the linear theory just after the quench, we see the growth of the spectrum in high frequencies as time passes, apparently because of the nonlinearity in the time evolution equations. In Fig. 11(c) and 11(d), the spectra deviate from those in the linear theory even at t = ∆t, and this indicates that the nonlinearity evolves even in the quenching time, 0 < t < ∆t. For these three cases, we find a direct energy cascade: The energy is transferred to higher n-modes during the time evolution. Eventually, power law spectra are observed, and these behaviors persist until the time of cusp formation. (See magenta points.) 6 In particular, as seen in 11(b), the time for reaching the power law behavior can be rather earlier than that for the cusp formation, and once realized, the behavior lasts until the cusps are formed. 7 Hence, the cusp formation on the string can be regarded as a variation of weak turbulence [14]. 8 We fit the power law spectra by ε n ∝ n −a . Just before the cusp formation, we obtain a = 1.313 ± 0.171, a = 1.354 ± 0.128 and a = 1.432 ± 0.017 for (b), (c) and (d), respectively. The exponents are distributed around a ∼ 1.4.

Time-dependence of the forces acting on the heavy quarks
Following the procedure in Section 4.5.3, we compute the forces acting on the quark and antiquark as functions of time. Since the string motion is now restricted in the 6 After the cusp formation, the function R(t, φ) becomes multi-valued and the energy spectrum is ill-defined. Therefore, we only show spectra before the cusp formation. 7 The string is smooth before the cusp formation, and the energy spectrum ε n must fall off faster than any power law function as n → ∞. Hence, it is indicated that the power law spectrum will be no longer maintained at n ≫ 50, while many numerical efforts are necessary for computing the spectrum in such a region. 8 The AdS weak turbulence was found to be characterized by a Kolmogorov-Zakharov scaling [40]. (z, x 1 )-plane, only the x 1 -components of the forces are non-zero. In Fig. 12, we plot the forces F 1 (t) and F 1 (t) for ǫ = 0.005 and ∆t/L = 2. For these parameters, cusps do not form on the string. Figure 12(a) is for the early stage in the time evolution (0 ≤ t/L ≤ 20), and Fig. 12(b) for a long period (0 ≤ t/L ≤ 160). In Fig. 12(a), pulse-like oscillations are repeated at regular time intervals. The period corresponds to the timescale for the fluctuation induced by the quench (4.1) to go back and forth between the two boundaries. Because of the dispersive spectrum in the linear perturbation, the initially localized oscillations tend to spread as time increases as in Fig. 12(b), and there is no way for them to converge again. We can also see that F 1 (t) < 0 and F 1 (t) > 0 throughout the time evolution, and this implies that the force between the quark and antiquark is always attractive. For the parameters where the cusps form on the string, the time evolution of the forces is different from the previous example. In Fig. 13, we show F 1 (t) and F 1 (t) for ǫ = 1.0 × 10 −2 and ∆t/L = 2, where figures (a) and (b) are for 0 ≤ t/L ≤ 13 and for 0 ≤ t/L ≤ 30, respectively. In figure (b), we take the absolute values of the forces, and plot in the log scale in the vertical axis. We find that the pulse-like oscillations are getting sharp and amplified as the time increases. The forces can change the sign because of their large oscillations, and this implies that the force between the quark and antiquark can be repulsive temporarily. Eventually, when a cusp arrives at the boundary after the cusp formation, the force diverges.
A mathematical explanation why the forces diverge after the cusp formation is given as follows. Near the AdS boundary, the time dependent solution is well approximated by Eq. (4.14). The conditions for the cusp are then given by ∂ σ Z = 0 and ∂ σ X = 0 as discussed in Section 4.5.1. The latter condition is automatically satisfied near the boundary because of ∂ σ X ∼ σ, while the former givesṫ 0 = 0. In Eq. (4.27), the denominator hasṫ 0 , and it appears that the numerator does not cancel the zero in the denominator. Hence, it is natural that the force diverges when the cusp arrives at the boundary.

Z 2 -symmetric quench
When the two endpoints of the string are simultaneously quenched in the opposite directions with the same amplitude, there is a new contribution from the one-sided quench that the propagating waves collide at the Z 2 -symmetric point X = 0, and cusps are expected to form on the collision. This case is equivalent to impose the Neumann boundary condition at X = 0. Such a condition is typically imposed in probe D-brane embeddings. We would like to emphasize that understanding the difference between this case and the case without the Neumann condition would be  (4.20) to be satisfied for the first time, and the corresponding φ-coordinates are plotted in the right panel. We find that there are two kinds of cusp formation and wave collisions at the Z 2 -symmetric point: One is cusp formation by the wave collision, and the first cusp formation in this case is marked with red points in the plots, since in this case these cusps disappear once the colliding waves pass. The other is that cusps are formed on the propagating waves in the same way as in the one-sided quench, and the cusp formation for this case is marked with green triangles. The times for the cusp formation are clearly discretized, and the locations are concentrated to φ = β 0 /2. The small change in the formation time is because the propagating speed of the wave slightly varies for different ǫ.
Closely looking at the results in evaluating (4.20), we find that even numbers of cusps are created for the first case (red points). In the bulk coordinates, the cusp formation is seen at very close to φ = β 0 /2 and one might naively think that only one cusp was generated at a point. In the worldvolume results, however, when the waves collide at φ = β 0 /2, we find that a pair of cusps whose orientations are opposite are created, and then these cusps pair-annihilate shortly when the waves pass by. Hence, the cusps created by the collision do not propagate away from that point, and therefore these instantaneous cusps are not observed in the force at the boundary. In this case, other formation of cusps in the same way as that in the one-sided quench also happens afterward. In the cases of the green triangles in Fig. 14, the cusps are formed slightly before the collision point, and these cusps subsequently collide at the Z 2 -symmetric point. In fact, it is seen in Fig. 14(a) that these points go ahead of the cusp formation by the collisions. These cusps then continue traveling on the string, inferring that the waves are already magnified enough for forming cusps.
In the Z 2 -symmetric case, it is convenient to evaluate the time evolution of the worldsheet Ricci scalar at φ = β 0 /2 since many events of the cusp formation occur there. The Ricci scalar is given by For the static configuration, this monotonically changes from R = −2/ℓ 2 at φ = 0 to R = −4/ℓ 2 at φ = β 0 /2, and because of this coordinate dependence it might be desirable to compare the Ricci scalar at a fixed φ. The Ricci scalar diverges on top of a cusp, and the cusp formation condition (4.20) is consistent with the divergence of the Ricci scalar since γ uv in the denominator becomes zero when (4.20) is satisfied. Practically, as we use discretized computations, the Ricci scalar does not exactly become infinity, while at least it becomes huge. Results of the Ricci scalar representing the cusp formation at second, third, and forth collisions are shown in Fig. 15. In these results, the absolute value of the Ricci scalar suddenly becomes huge at the cusp formation, and come back to of order one when the waves pass. In our setup, the numerical evolution does not breakdown after the Ricci scalar diverges in contrast to the D3/D7 case [10][11][12], although finite-N c corrections may have to be taken into account.  6 Results for the transverse linear quench 6

.1 Cusp formation
In this section, we study the string dynamics induced by the transverse linear quench (4.4). With such a quench, the string moves in the (3+1)-dimensions spanned by (t, z, x 1 , x 2 ). Snapshots of string configurations under a quench with parameters ∆t/L = 2 and ǫ = 0.03 are shown in Fig. 16

Energy spectrum in the non-linear theory
We expect to see nonlinear origin for the cusp formation in the energy spectrum also in the case of the transverse linear quench. In Fig. 17, we show the time dependence of the energy spectrum when the parameters are ∆t/L = 2 and ǫ = 0.03. The time for the cusp formation is t/L = 14.45 for these parameters. Red points are just after the quench, and magenta points correspond to the time at cusp formation. We find the direct energy cascade as in Section 5.2, and eventually the spectrum obeys a power law until the time of cusp formation. Thus, also in the transverse linear quench, we find the turbulent behavior toward the cusp formation. Fitting the spectrum by ε ∝ n −a , we obtain a = −1.533 ± 0.133. The exponent has similar value as that for the longitudinal quench.

Time-dependence of forces acting on the heavy quarks
We turn to the time-dependence of the forces acting on the quark and the antiquark in the transverse linear quench. Since the motion of the string is in the (z, x 1 , x 2 )space, the x 1 -and x 2 -components of the forces can be non-zero. In Fig. 18(a), we show F and F as functions of time for ∆t/L = 2 and ǫ = 0.01, with which cusps do not appear on the string. Similar to the longitudinal quench, pulse-like oscillations are repeated at intervals. Although there are sharp peaks, they are always O(1) in units of λ −1/2 L 2 . We also find that F 1 (t) < 0 and F 1 (t) > 0. Thus, the force between quarks is always attractive. In Fig. 18(b), we show the absolute values of the forces for ∆t/L = 2 and ǫ = 0.03, where cusps are formed as seen in Section 6.1. The pulses are getting sharp and amplified as time increases, and eventually after the cusp formation, the forces diverge when the cusps arrive at the boundaries. We also monitoredṫ 0 , which is in the denominator in Eq. (4.27), at the boundaries, and found that it is consistent with zero at t/L ≃ 13.6 and 15.9 at x = x q and xq, respectively.

Results for the transverse circular quench
Finally, we study the string dynamics induced by the transverse circular quench (4.5).
The string moves in all (4+1)-dimensions spanned by (t, z, x 1 , x 2 , x 3 ). For the transverse circular quench, we did not find any cusp formation at least for modest parameters: around at ǫ ∼ 0.01 and ∆t/L ∼ 1. Nevertheless, we found an interesting behavior in the energy spectrum. In Fig. 19, we show the time dependence of the energy spectrum for parameters ∆t/L = 2 and ǫ = 0.02. In the early time evolution until t/L 14, there is a direct energy cascade: the energy is transferred from large to small scales, and eventually the spectrum obeys a power law at t/L ∼ 14. Fitting the numerical data at t/L = 14, we obtain ε n ∝ n −2.019±0.029 . For t/L 14, however, we find that this turns into an inverse energy cascade: The energy is transferred to the large scale. Thus, in contrast to the previous low dimensional cases, the power law once realized at an intermediate time t/L = 14 is not maintained in the late time.
In Fig. 20, we show snapshots of string configurations around the "turning point" of the energy cascade: t/L = 14, 14.2, 14.4. Since the string motion is in the (4 + 1)-dimensions, we project the string profile into (x 1 , x 2 , z)-and (x 1 , x 3 , z)-spaces. Although cusp-like points can be seen in the right figure, these are not real cusps: We find that although roots of J Z , J X 2 and J X 3 become close at t/L ≃ 14, J X 1 is not zero at the point. (See Eq. (4.20).) However, the perturbation variableχ 3 defined in Eq. (4.23) becomes cuspy, namely, its energy is transferred to the small scale. Hence the direct energy cascade appears in t/L 14. After t/L ≃ 14, the cuspy shape gets loose because of the dispersive spectrum in the linear perturbation.
In Fig. 21, we show the time dependence of the forces acting on the quark and  the antiquark for ǫ = 0.02 and ∆t/L = 2. We do not find the divergence of the forces. However, the forces are magnified until t/L ≃ 14 and can be repulsive at some time intervals ( F 1 > 0 and F 1 < 0) even though there is no cusp formation. The magnitude of the forces in the late time is not as big as that period, reflecting the looseness in the cuspy shape.

Summary and Discussion
We studied nonlinear dynamics of the flux tube between an external quark-antiquark pair in N = 4 SYM theory using the AdS/CFT duality. We numerically computed the time evolution of the string in AdS dual to the flux tube when we perturbed the positions of the string endpoints to induce string motions. We considered four kinds of quenches that were chosen to represent typical string motions: (i) longitudinal onesided quench, (ii) longitudinal Z 2 -symmetric quench, (iii) transverse linear quench, In the time evolution of the energy spectrum, we observed the weak turbulence, that is, the energy was transferred to the small scale, and the energy spectrum eventually obeyed a power law until the time of the cusp formation. The cusp formation occurred only when the amplitude of the quench was larger than a critical value, ǫ > ǫ crit , and the dependence of its magnitude on the quench duration ∆t was given by a simple form ǫ crit ∝ (∆t/L) 3 in small ∆t/L. When the cusps arrived at the AdS boundary, we observed the divergence of the force between the quark pair. For (iv), we found no cusp formation. Nevertheless, we observed a direct energy cascade and the power law spectrum for a while. However, in late time the direct cascade turned into an inverse energy cascade, where the energy was transferred to the large scale. There was no divergence of the force between the quark pair. How can we understand the weak turbulence of the string in view of gauge theory? Eigen normal modes e n of the fundamental string studied in section 3 can be regarded as the excited states |n of the flux tube in the gauge theory side. Hence, the fluctuating string solution, such as R(t, φ) = z 0 + n c n (t)e n (φ), corresponds to c n (t)|n (8.1) in the boundary theory, where |0 is the ground state. The weak turbulence implies that |c n (t)| with n ≫ 1 tends to increase as a function of time. Therefore, in the late time, the probability of observing highly excited states is high compared to the linear theory. Similar phenomenon has been found in the D3/D7 system dual to N = 2 supersymmetric QCD [10][11][12], where a direct energy cascade was found in the fluctuations on the D7-brane and regarded as production of many heavy mesons in the SQCD. In that paper, this phenomenon was referred as "turbulent meson condensation". Although the endpoints of the flux tube we considered are regarded as nondynamical and infinitely heavy quarks, the string turbulence found in this paper would be regarded as the microscopic picture of the turbulent meson condensation. We found cusp formation when the motion of the string is restricted in (2 + 1)-and (3 + 1)-dimensions. The divergence of the forces acting on the quarks is accompanied by the cusp singularities, where we expect that finite-N c effects will become important. These will contain quantum effects of the string, and such effects may resolve the cusp singularities and the divergence of the forces. Nevertheless, the cusp formation in the classical sense can give us observable effects: Finite-N c effects will also appear as gravitational backreactions. If these are taken into account, a strong gravitational wave will be emitted at the onset of the cusp formation. 9 For cosmic strings in flat spacetime, gravitational wave bursts from cusps have been studied in Ref. [42], and it has been found that their spectra obey a power law in the high-frequency regime. It would be nice to compute the gravitational waves from the string in AdS 5 and find the power law spectrum. The description in the dual field theory may be the power law spectrum in gluon jets from the flux tube.
When the motion of the string is in (4 + 1)-dimensions, we did not find cusp formation. Hence, in AdS 5 spacetime, the cusp formation on the string is not a general phenomenon but accidental one. This implies that the dual phenomenon to the cusp formation is not ubiquitous in the (3 + 1)-dimensional boundary field theory. However, if we consider the many-body system of quark-antiquark pairs, it is possible that the time evolution of some flux tubes happens to be restricted in lower dimensional spaces approximately, 10 and such flux tubes would be able to emit the gluon jets with the approximately power law spectrum, which is characteristic to the cusp formation. Besides, gravitational wave bursts in the presence of extra dimensions were discussed in [43,44]. Even though real cusps are not formed, there may be some gluon emission from cuspy shapes.
There are some future directions in our work. In this paper, we only considered modest values for the amplitude of the quench, ǫ ∼ 0.01. For a large value of ǫ, we expect that the string can even plunge into the Poincare horizon because of the strong perturbation. This will demonstrate a non-equilibrium process of breaking of the flux tube. It is also straightforward to take into account finite temperature effects in this process. It will be also interesting to consider the string motions in confined geometries [17]. In the theories dual to these backgrounds, the quarkantiquark potential is linear, and the presence of such potential may affect conditions for cusp formation. Studying nonlinear string dynamics in such geometries may give new insights into understanding the QCD flux tubes and non-equilibrium processes in realistic QCD.
Closed strings rotating in AdS and having cusps were constructed in Ref. [45]. Although our dynamical cusp formation on an open string is different from the existence of cusps in those steady solutions, it may be interesting to obtain useful information from such configurations. In [10], the universal exponent in the power law was deduced from a stationary solution called critical embedding in the D3/D7-brane system in the presence of a constant electric field, and results in time dependent computations obeyed that universal value. In our setup, we do not have a corresponding static cuspy configuration, but our power law exponents, distributed around 1.4, may be naturally understood from cuspy stationary strings in AdS.
Ultimately, it will be important to understand the mechanism relevant for the the turbulent behavior. For the integrable Wilson loops such as those in AdS 5 × S 5 , the turbulent behavior may be studied with the techniques of integrability.
In non-linear systems, there may be underlying chaos. In [46], closed strings moving in Schwarzschild-AdS background were studied from the viewpoint of chaos. It may be interesting if chaos is seen also in the motion of open strings and the turbulent behavior is understood, particularly in non-integrable situations. where we choose the positive signs for the square roots since we take ∂ u and ∂ v as future directed vectors. Eliminating T ,u and T ,v from Eq. (4.11), we obtain The evolution equations in these expressions are found numerically stable.
To numerically solve (A.2), we discretize the world volume (u, v)-coordinates with the grid spacing h as shown in Fig. 22. Let us denote the fields (T, Z, X) by Ψ. At a point C apart from the boundary, the fields and their derivatives are discretized with second-order accuracy as Discretization error is O(h 2 ). Substituting these into the evolution equations (A.2), we obtain nonlinear equations to determine Ψ N by using known data of Ψ E , Ψ W , and Ψ S . We use the Newton-Raphson method for solving the coupled nonlinear equations. The equations for the boundary time evolution (4.15) become coupled nonlinear equations of T N and X N . (Z N = 0 is trivially imposed.) The form at σ = 0 is where we used a relation for Z, Z E = −Z W , derived from the boundary condition Z ,uv = 0. At the other boundary σ = β 0 , Z W and x q in (A.4) are replaced with Z E and xq. These equations are also solved by using the Newton-Raphson method.

B Error analysis
In this section, we estimate errors in our numerical calculations. We definẽ These constraints should be zero for exact solutions. Hence, these can be nice indicators of our numerical errors. For visibility of the constraint violation, we introduce where we take the maximum value when we vary u on a fixed v surface. We also choose the bigger of the two constraints, |C 1 | and |C 2 |. Introducing an integer N such that the mesh size is given by h = β 0 /N, we plot C max (v) for several values of N in Fig. 23. We see that the constraint violation is small (C max ∼ 10 −3 even for N = 200) and behaves as C max ∝ 1/N 2 . This is consistent with the fact that our numerical method has the second order accuracy. In this paper, we mainly set N = 800. Then, the constraint violation is O(10 −4 ).

C Energy spectrum in the linear theory
In this appendix, we derive the energy spectrum induced in the linear theory of section 3 when a quench is added on the boundary. The equations of motion for the perturbation variables (χ 1 , χ 2 , χ 3 ) are given in Eq. (3.8). Here, we focus on the longitudinal mode χ 1 for simplicity and denote it as χ 1 = χ. Application to transverse modes is straightforward. For the quench, we consider the following boundary conditions for χ: where χ b (t) is the quench function assumed to have a compact support at 0 < t < ∆t. We also assume that the solution is trivial before the quench, χ(t ≤ 0, φ) = 0. We firstly consider a time independent solution to (3.8), A solution is given by is a constant for normalization. Near the boundaries the function S behaves as Let us introduce the quench χ b (t). Using the function S, we defineχ as The new variableχ satisfies trivial boundary conditions,χ(t, φ = 0) =χ(t, φ = β 0 ) = 0. The equation of motion forχ becomes where we used Eq. (C.2).
To solve the equation, we consider a Green's equation: where G is the Green's function. By using G, a special solution to Eq. (C.6) can be written in the formχ Operating t ′ +ǫ t ′ −ǫ dt on both sides of the Green's equation and taking the limit of ǫ → 0, we obtain the junction condition as For t < t ′ , we assume the Green's function is trivial: G = 0. For t > t ′ , the Green's function can be written as G = n [a n sin ω n (t − t ′ ) + b n cos ω n (t − t ′ )]e n (φ) , (C.10) From the continuity of G at t = t ′ we have b n = 0, and then from the junction condition (C.9), we find that a n satisfies n a n ω n e n (φ) = δ(φ − φ ′ ) . (C.11) Operating (e n , * ) to the above equation, we obtain a n = 1 ω n γ(φ ′ )e n (φ ′ ) . (C.12) Thus, the Green's function can be written as n ω −1 n sin ω n (t − t ′ ) γ(φ ′ )e n (φ ′ )e n (φ) (t > t ′ ) .

(C.13)
Since the Green's function is zero at t < t ′ , the special solution obtained from Eq. (C.8) is also zero before the quench, t < 0, and this is nothing but the solution we are looking for. After the quench t > T , the solution becomes (C.14) where S n ≡ (S, e n ). Note thatχ = χ after the quench. At the second equality, we integrated by parts twice. It is then straightforward to compute the energy spectrum. The mode coefficient c n = (χ, e n ) is computed by using (C.14) as c n (t) = ω n S n T 0 dt ′ χ b (t ′ ) sin ω n (t − t ′ ) , (C. 15) and then from Eq. (4.25) the energy spectrum for the longitudinal quench in the linear theory becomes ε n = √ λz 0 4π ċ 2 n + ω 2 n c 2 n = √ λz 0 ω 4 n S 2 n 4π |χ(ω n )| 2 .
(C. 16) where we defineχ(ω) = dtχ b (t)e −iωt . The energy spectrum does not depend on t as we expect. Taking into account the transverse modes, we obtain the energy spectrum as ε n = √ λz 0 4π ω 4 n S 2 n |χ(ω n )| 2 + and χ i b are the quench functions for χ i . The total energy is ε = ∞ n=1 ε n . The spectrum (C.17) is defined only for integer n. However, in Figs. 11, 17 and 19, we generalize ε n to a continuous number by interpolating ω n , S n , ω ′ n , and S ′ n and regarding them as function of continuous number n for visibility.

D Forces acting on the quark and the antiquark
In this appendix, we derive the formula for the forces acting on the quark endpoints (4.27). We denote the on-shell Nambu-Goto action as S[x q , xq], where x q and xq are the locations of the string endpoints regarded as the quark and the antiquark, respectively, at the AdS boundary: X(t, z → 0) = x q , xq. The on-shell action relates to the partition function of the boundary theory as In the field theory, the partition function is written as where φ represents the set of the fields in N = 4 super Yang-Mills theory and S SYM is its action. S q and Sq denote the actions for the quark and the antiquark, respectively. We regard x q (t) and xq(t) as external fields. The quark action is schematically written as where m is the quark mass and L int corresponds to the interaction term with SYM fields, and the velocity of the quark is introduced as v ≡ dx q /dt. We define the force acting on the quark as The variation of the on-shell Nambu-Goto action is then related to the field theory terms as δS δx q = −i δ δx q ln Z CFT = −m(γv) · + F , (D.5) where F = Z −1 DφF e iS SYM +iSq+iSq , γ = 1/ √ 1 − v 2 , and the dot in the upper right of the parentheses denotes · ≡ d/dt. In the first equality, we used the AdS/CFT duality (D.1). Now, we evaluate δS/δx q in the gravity side. For this purpose, it is convenient to use target space coordinates (t, z) as the world sheet coordinates. Then, the position of the string is specified by x = X(t, z), and the Nambu-Goto action becomes The equations of motion are given by where · ≡ ∂/∂t, ′ ≡ ∂/∂z, and ξ ≡ (1 −Ẋ 2 )(1 + X ′2 ) + (Ẋ · X ′ ) 2 . Solving these near the AdS boundary, we obtain where a = d 2 x q /dt 2 . Let us consider the variation of the on-shell action with respect to one of the endpoints of the string, x q , and the other endpoint is fixed, δxq = 0. By defining Lagrangian as S = dτ dσL, the variation of the action becomes δS[x q , xq] = dtdz δ(∂ a X) · ∂L ∂(∂ a X) + δX · ∂L ∂X = dtdz ∂ a δX · ∂L ∂(∂ a X) = − dt δX · ∂L ∂X ′ δS[x q , xq] = dt δx q · − √ λ 2πǫ (γv) · + 3 √ λ 2πγ (x 3 + γ 2 (v · x 3 )v) . (D.10) Hence, the upshot for δS/δx q is Comparing above expression with Eq. (D.5), we can see that the first term in (D.11) corresponds to a diverging quark mass m ∼ 1/ǫ. This is a natural consequence since we are considering an infinitely extended string. Setting m = √ λ/(2πǫ), we obtain the force acting on the quark from Eqs. (D.5) and (D.11) as where we replaced x 3 with ∂ 3 z X| z=0 /6. Note that there can be a finite difference between m and √ λ/(2πǫ), and an extra-term proportional to (γv) · may appear in Eq. (D.12). However, it can be eliminated by adding a local counter term proportional to √ 1 − v 2 in the quark action S p . The same formula can be applied to the force acting on the antiquark at xq.