Chaos of QCD string from holography

TL;DR

This paper investigates chaos in large N_c, strongly coupled QCD strings using holography, measuring Lyapunov exponents and identifying conditions under which chaos emerges.

Contribution

It demonstrates that QCD strings exhibit chaos at strong coupling via holographic duality and quantifies this chaos through Lyapunov exponents.

Findings

Chaos appears when boundary force exceeds a certain threshold.

The chaos bound increases with quark-antiquark separation.

Chaos originates from the endpoints of the QCD string.

Abstract

It is challenging to quantify chaos of QCD, because non-perturbative QCD accompanies non-local observables. By using holography, we find that QCD strings at large $N_c$ and strong coupling limit exhibit chaos, and measure their Lyapunov exponent at zero temperature. A pair of a quark and an antiquark separated by $L_q$ in the large $N_c$ QCD is dual to a Nambu-Goto string hanging from the spatial boundary of the D4-soliton geometry. We numerically solve the motion of the string after putting a pulse force on its boundaries. The chaos is observed for the amplitude of the force larger than a certain lower bound. The bound increases as $L_q$ grows, and its dependence is well approximated by a hypothesis that the chaos originates in the endpoints of the QCD string.

Figures

Equations

$$d{s}^2$$

$$+{\left(\frac{\mathcal{R}}{U}\right)}^{3/2}\frac{d{U}^2}{f(U)}+{\mathcal{R}}^{3/2}{U}^{1/2}d{\Omega }_4^2,$$

$$f(U)=1-\frac{U_{\text{KK}}^3}{U^3}.$$

$$\delta \tau =\frac{4\pi }{3}\frac{{\mathcal{R}}^{3/2}}{U_{\text{KK}}^{1/2}}=\frac{2\pi }{M_{\text{KK}}},$$

$${\mathcal{R}}^3=\frac{1}{2}\frac{g_{\text{YM}}^2{N}_c{l}_s^2}{M_{\text{KK}}},U_{\text{KK}}=\frac{2}{9}{g}_{\text{YM}}^2{M}_{\text{KK}}{l}_s^2.$$

$$U(r)=U_{\text{KK}}(1+{\tan }^{2}r).$$

$$d{s}^2=\frac{4}{3}\lambda l_s^2\frac{1}{{\cos }^{3}r}\left[\frac{M_{\text{KK}}^2}{9}(-d{t}^2+d{\vec{x}}^2)+\frac{d{r}^2}{1+{\cos }^{2}r+{\cos }^{4}r}\right],$$

$$S_{\text{NG}}=-\frac{1}{2\pi {\alpha }^{\prime }}\int d\tau d\sigma \sqrt{-h},$$

$$S_{\text{NG}}$$

$$E=-\frac{2}{\mathcal{T}}{S}_{\text{NG}}[\tilde{X}].$$

$$E_0=-\frac{2}{\mathcal{T}}\left(\frac{\lambda }{6\pi }\mathcal{T}{\int }_{\pi /2}^{0}dr\frac{1}{{\cos }^{3}r}\frac{2}{\sqrt{1+{\cos }^{2}r+{\cos }^{4}r}}\right).$$

$$S_{\mathrm{NG}}=$$

$$+6{\int }_{r_0}^{\pi /2}\frac{dr}{{\cos }^{3}r\sqrt{1+{\cos }^{2}r+{\cos }^{4}r}}].$$

$$L_q{M}_{\text{KK}}=\frac{2\cos r_0}{\sin r_0\sqrt{1+{\cos }^2{r}_0+{\cos }^4{r}_0}}.$$

$$r_0=\frac{2}{\sqrt{3}}\frac{1}{L_q{M}_{\text{KK}}}+\frac{4}{9\sqrt{3}}\frac{1}{(L_q{M}_{\text{KK}}{)}^3}+\cdots .$$

$$r=r_0+\alpha (t)+\beta (t)\frac{x}{L_q}(-L_q/2\le x\le L_q/2)$$

$$S=\frac{-2\lambda M_{\text{KK}}}{9\pi }\int dt\left({\mathcal{L}}_{\text{radial}}+{\mathcal{L}}_{\text{bottom}}\right)$$

$${\mathcal{L}}_{\text{radial}}=$$

$$+{\int }_{r_0+\alpha +\beta /2}^{\pi /2}\frac{d\sigma }{{\cos }^3\sigma \sqrt{1+{\cos }^2\sigma +{\cos }^4\sigma }},$$

$${\mathcal{L}}_{\text{bottom}}=\frac{M_{\text{KK}}}{3}{\int }_{-L/2}^{L/2}d\sigma \frac{1}{{\cos }^{3}r}$$

$$\times {\left[1-\frac{9}{M_{\mathrm{KK}}^2}{\left(\dot{\alpha }+\dot{\beta }\frac{\sigma }{L}\right)}^2\frac{1}{1+{\cos }^{2}r+{\cos }^{4}r}\right]}^{1/2}$$

$$\times {\left[1+\frac{9}{L_q^2{M}_{\text{KK}}^2}{\beta }^2\frac{1}{1+{\cos }^{2}r+{\cos }^{4}r}\right]}^{1/2}$$

$$S=$$

$$V\equiv$$

$$+a_{111}{\alpha }^3+a_{122}\alpha {\beta }^2$$

$$+b_{111}\alpha {\dot{\alpha }}^2+b_{212}\beta \dot{\alpha }\dot{\beta }+b_{122}\alpha {\dot{\beta }}^2.$$

$$E_{\text{chaos}}\equiv \frac{2b_{11}^3}{a_{111}^2}$$

$$E_{\mathrm{chaos}}=\frac{3\lambda M_{\mathrm{KK}}}{98\pi }(L_q{M}_{\mathrm{KK}}{)}^3\left(1+\mathcal{O}\left(\frac{1}{(L_q{M}_{\mathrm{KK}}{)}^2}\right)\right),$$

$$\displaystyle S_{\text{NG}}=-\frac{\lambda}{6\pi}\int dudv\frac{1}{\cos^{3}R}\bigg{[}$$

$$\displaystyle\quad+\frac{4R_{,u}R_{,v}}{1+\cos^{2}R+\cos^{4}R}\bigg{]}.$$

$$T_{,uv}$$

$$R_{,uv}$$

$$+\frac{3}{8}\tan R(1+{\cos }^{2}R+{\cos }^{4}R)(-T_{,u}{T}_{,v}+{\vec{X}}_{,u}\cdot {\vec{X}}_{,v}),$$

$${\vec{X}}_{,uv}$$

$$C_u$$

$$C_v$$

$$X_1(\tau ,\sigma =\delta )=\frac{L_q}{2}+ϵ\alpha (\tau ;\Delta \tau ),$$

Full text

††preprint: OU-HET-1000

I Introduction: QCD chaos

How chaotic is QCD? — a question which is simple but unanswered, should drive the understanding of our universe based on quantum field theories. It is challenging to define the extent of chaos for QCD, because QCD is truly quantum while the popular measure of chaos, the Lyapunov exponent, is defined classically. Analyses based on weakly coupled picture Matinyan et al. (1981, 1986); Muller and Trayanov (1992); Biro et al. (1995); Gong (1994); Kunihiro et al. (2010); Muller and Schafer (2011); Iida et al. (2013); Tsukiji et al. (2016) and on out-of-time ordered correlators Maldacena et al. (2016) (which define a quantum chaos) suggest a QCD chaos at high temperature, but what about the usual picture of hadronic phase of QCD?

As lattice QCD, a popular strategy to study non-perturbative nature of QCD, still lacks a way to follow time dependence necessary to analyze any chaos, we need some other way. The holography, or the AdS/CFT correspondence Maldacena (1999), is suitable for the purpose. Taking a large $N_{c}$ limit and a strong coupling limit, QCD is approximated by a classical gravity dual, while keeping the quantum nature and the time dependence of QCD. In this paper, we analyze chaos of a quark antiquark pair by using the holography. We find a condition for the chaos to occur, and draw a phase diagram of the QCD chaos.

The study of chaos in the AdS/CFT was initiated in Pando Zayas and Terrero-Escalante (2010). While the chaos of chiral condensate in QCD was studied in Hashimoto et al. (2016); Akutagawa et al. (2018) via the holography, the physical excitation of QCD at low energy is non-local. Wilson loops, and pairs of quark-antiquark connected by an open Wilson loop, are the low energy physical degrees of freedom of QCD and their quantization provides the hadronic world. Spectra of hadrons exhibit quantum chaos Pascalutsa (2003), so, we need to locate the origin of the chaos of QCD, and measure the extent of the QCD chaos, based on the non-local Wilson loops.

Here we have to remind the readers of the fact that a Nambu-Goto (NG) closed string in three spatial dimensions, a phenomenological model of glueballs in the large-$N_{c}$ QCD, is integrable. Then, what is the origin of the QCD chaos? Naively, we can expect two possible origins: one is the boundary of the NG string, which is the quark, and the other is the thickness of the QCD string which has not been taken care of for the three-dimensional NG string. The question can be addressed in holography, because the QCD string corresponds to a NG string in the higher-dimensional spacetime in the gravity dual, and its static nature, such as the quark boundaries and the thickness, has been well-studied. We here provide a detailed analysis of a quark antiquark pair in motion, and locate the origin of the QCD chaos.

Through the AdS/CFT, The Wilson loop in QCD is identified with a NG string Maldacena (1998); Rey and Yee (2001) hanging down from the boundary of confining geometry Witten (1998a) which is considered to be a dual to a pure 4-dimensional Yang-Mills theory. The $q\bar{q}$ potential is the free energy of the string. Since the geometry has the bottom of the spacetime, the hanging string has the part sitting at the bottom, which provides the QCD string tension, and the parts connecting the bottom and the boundary, which correspond to the quarks in a gluon cloud. The motion of the QCD string is caused by a pulse force acting on the infinitely massive quarks, and we solve numerically the motion of the NG string in the geometry. The chaotic Lyapunov exponent is observed when the strength of the pulse force exceeds a certain bound. We study the dependence of the bound on the interquark distance, and find that the $q\bar{q}$ pair is less chaotic for larger interquark distances.

Our numerical result is explained well by a popular effective picture of the quarks connected by a long QCD string, assuming that the QCD string motion is integrable while the endpoint regions (the quarks with a gluon cloud) are chaotic. It suggests that the chaos of QCD string originates in its endpoints. The chaos of motion of closed NG string in various geometry has been studied Basu and Pando Zayas (2011); Basu et al. (2012); Pando Zayas and Reichmann (2013) (see also Basu and Ghosh (2014); Giataganas et al. (2014); Giataganas and Sfetsos (2014); Bai et al. (2016); Asano et al. (2015a); Ishii and Murata (2015, 2016); Basu et al. (2017)), which corresponds to the chaos due to the thickness of the QCD string. Our study about the quark antiquark pair shows a different origin of the QCD chaos.

The organization of this paper is as follows. First, in Sec. II, we review the static Wilson loop in the holographic QCD, and introduce our coordinate system in the bulk. In Sec. III, a rectangular NG string is introduced as a toy model, and its fluctuation analysis is presented to show the existence of the chaos and the chaos energy bound. Our main numerical study of the NG string in motion in holography is presented in Sec. IV. There we find Lyapunov exponent of the string motion, and draw a phase diagram of chaos, as a function of $\epsilon$ which is the magnitude of the pulse force and $L_{q}$, the interquark distance. The chaotic behavior is observed in the interquark force which is an observable of QCD. Sec. V is for a discussion to locate the origin of chaos. We introduce a simple effective picture of an open QCD string and discuss its chaos, to fit the numerical result of the phase diagram obtained in Sec. IV. We conclude that the chaos originates in the endpoints of the QCD string, the quarks. Sec. VI is for a summary and discussions. App. A provides details of numerical calculations. App. B calculates the formula for the interquark force.

II Confining geometry and $q\bar{q}$ potential

In this section, we first review the confining geometry Witten (1998a) and compute the static $q\bar{q}$ potential through the AdS/CFT, based on the dictionary Maldacena (1998); Rey and Yee (2001). The string configuration serves as an initial one upon which an external pulse force is put to produce a time-dependent motion of the QCD string, later in Sec. IV.

The D4-soliton background holographically corresponds to a five dimensional super Yang-Mills theory on a non-supersymmetric circle, giving a four-dimensional pure Yang-Mills theory at low energy Witten (1998a). The background is an example of confining geometries which has the bottom of the spacetime. Let us first obtain the static NG string configuration hanging down from the boundary of the spacetime, to calculate the expectation value of the Wilson loop in the Yang-Mills theory.

The D4-soliton background is of the following form Witten (1998a); Kruczenski et al. (2004):

[TABLE]

The coordinates $x^{\mu}$ and $\tau$ are the directions along the D4-branes, and the $\tau$ direction is compactified on $S^{1}$. The coordinate $U$ is a radial direction transverse to the D4-branes. To avoid a conical singularity at $U=U_{\text{KK}}$, the period of the $\tau$ direction must be

[TABLE]

so $1/M_{\text{KK}}$ is the radius of $S^{1}$. Parameters in the metric can be expressed by those of the dual gauge theory as

[TABLE]

The motion of the NG string studied in Sec. IV often goes through the tip of the geometry $U=U_{\rm KK}$, thus we need a coordinate system which does not have a coordinate singularity there. The new coordinate $r$ is introduced as

[TABLE]

Then, the metric becomes

[TABLE]

where $\lambda=g_{\text{YM}}^{2}N_{c}$ is the ’t Hooft coupling. Since we are interested in QCD we do not consider $\tau$ and $\Omega_{4}$ directions in the following, and so here we have omitted them. In this coordinate, the bottom of the D4-soliton $r=0$ is totally regular. The asymptotic boundary of D4-soliton is now $r=\pi/2$, which was $U\to\infty$ in the $U$ coordinate.

The $q\bar{q}$ potential is given by the energy of a static NG open string in the geometry of the gravity dual Maldacena (1998); Rey and Yee (2001). We consider a Wilson loop with the quark-antiquark separation $L_{q}$, and take an ansatz that the string is extended in the $x^{1}$-$r$ plane and the endpoints of the string are located at $x^{1}=\pm L_{q}/2$. The NG action in the geometry (6) is

[TABLE]

where $h=\det(h_{ab})$ and $h_{ab}$ is an induced metric on the worldsheet. We take the static gauge: $(\tau,\sigma)=(t,r)$ and then the static solution is provided as $x^{1}=X_{1}(r)$.111 We will use small letters for target space coordinates and capital letters for functions specifying the string position.

For numerical calculations, we choose a unit $M_{\text{KK}}=3/2$. The NG action now becomes

[TABLE]

where ${}^{\prime}=\partial_{r}$ and $\mathcal{T}=\int dt$. The integration region is $\pi/2\leq r\leq r_{\rm center}$ where $r=r_{\rm center}$ is the point of the bottom of the hanging string, which should solve the equation $X(r_{\rm center})=0$ due to the parity symmetry $X(r)=-X(-r)$ following from our boundary condition. Solving the equations of motion, we obtain static configurations of string for each quark separations $L_{q}$.

Fig. 1 shows the configurations of the static string in the $x^{1}$-$r$ plane. When the quark-antiquark separation becomes larger, the tip of hanging string sticks to the bottom of the geometry $r=0$. This actually implies that confining potential appears.

Let us evaluate the $q\bar{q}$ potential holographically. Considering the on-shell NG action $S_{\text{NG}}[\tilde{X}]$, where $\tilde{X}(r)$ is a static solution to the equation of motion, the $q\bar{q}$ potential is given by

[TABLE]

This has a divergence which stems from the infinitely long string hanging from the boundary, but it can be naturally understood as the infinite quark mass. Subtracting the contribution of that, the $q\bar{q}$ potential turns out to be $E-E_{0}$, where

[TABLE]

The quantity in the parenthesis is (8) with $X^{\prime}=0$, except that the integration region is $\pi/2\leq r\leq 0$.

Fig. 2 shows the relation between the quark-antiquark separation $L_{q}$ and the $q\bar{q}$ potential $E-E_{0}$. When $L_{q}$ is large, the potential becomes linear in $L_{q}$, which means it is a confining potential. This is totally consistent with the study developed in Brandhuber et al. (1998); Greensite and Olesen (1998); Kinar et al. (2000).

In the next section, we consider a toy model of the motion of the string and study a chaos bound, and in Sec. IV, we investigate numerically the full time-dependent dynamics of the string in the D4-soliton background and interquark force in the gauge theory.

III Toy model of string in motion

Before getting into the full numerical simulation of the NG string in motion, we here first study the motion of a toy model string to look intuitively how the chaos shows up in the motion. The toy model assumes the shape of the string and its fluctuation modes: the shape of the toy string is rectangular, and the fluctuation modes are only of two types, one is the motion keeping the rectangular shape, and the other is the motion giving a linear slope at the bottom of the rectangular string, see Fig. 3.

As is seen in comparison to the actual shape in Fig. 1, this toy model could capture some intrinsic feature of the motion of the NG string. In fact, the rectangular string hanging down from the boundary has been used in many literature to mimic the QCD string holographically, and it was used in Hashimoto et al. (2018) to show a universal chaos behavior near black hole horizons.

We here show that the toy model has no chaos when the total energy of the string is small, while for a larger energy the chaos appears. We estimate the lower energy bound of the chaos in the model, and the bound is used in the next section for intuitively understanding the simulation results. The existence of the chaos bound itself is easy to understand: For a very small fluctuation around the static string shape, the motion is that of a harmonic oscillator, so there is no chaos. On the other hand, if one puts a larger energy, any modes are excited and interacting with each other, generally inducing chaos. We are interested in the energy lower bound and its dependence on the distance between the quarks, $L_{q}$.

To obtain the fluctuation action of the toy string model, first we determine the static stable configuration. Denoting the location of the bottom of the rectangular string as $r=r_{0}$ with $-L_{q}/2\leq x\leq L_{q}/2$, the NG action is

[TABLE]

Extremizing of this action with respect to $r_{0}$, we obtain the relation between $r_{0}$ and $L$ for the static stable rectangular string,

[TABLE]

In particular, for a large interquark distance $L_{q}M_{\text{KK}}\gg 1$, this relation is rephrased as

[TABLE]

Let us proceed to obtain the fluctuation action. We include the lowest two modes, which are represented by the following linear shape of the bottom of the string,

[TABLE]

The fluctuation modes $\alpha(t)$ and $\beta(t)$ deform the bottom of the string. With this shape it is straightforward to obtain the NG action,

[TABLE]

where the radial and the bottom parts of the action of the string are given by

[TABLE]

and

[TABLE]

with $r=r_{0}+\alpha+\beta\sigma/L$ substituted for the last expression. We expand the total action (15) to the third order in the fluctuations $\alpha(t)$ and $\beta(t)$. The result is

[TABLE]

The coefficients $a$’s and $b$’s are functions of $r_{0}$, namely, of $L_{q}$. The “potential” term $V(t)$ in general includes time-derivative terms.222 Even if we make a field redefinition of the form $\alpha\to\alpha+\alpha^{2}+\alpha\beta$ and $\beta\to\beta+\alpha\beta+\beta^{2}$ (which respects the $x$-party transformation $(\alpha,\beta)\to(\alpha,-\beta)$, we cannot absorb the third order terms including time-derivatives (which are $\alpha\dot{\alpha}^{2}$, $\beta\dot{\alpha}\dot{\beta}$ and $\alpha\dot{\beta}^{2}$).

Generically, for the chaos to occur, the interaction terms (the cubic terms) in $V(t)$ need to contribute. For small fluctuation, only the quadratic terms in $V(t)$ (which are mass terms for $\alpha(t)$ and $\beta(t)$) provide the full dynamics and it is just a set of harmonic oscillators. When the fluctuation is larger, the cubic interaction term contributes, and the chaos emerges. To estimate the typical value for the energy lower bound of the chaos to emerge, we pick up the terms of $\alpha(t)^{2}$ and $\alpha(t)^{3}$ in $V(t)$ and obtain the energy at which the values of these two terms are equal to each other, under the condition $\dot{\alpha}(t)=0$. This energy is

[TABLE]

We plot this chaos bound in Fig. 4. We find that the chaos energy bound diverges for $L_{q}\to\infty$ or $L_{q}\to 0$, while it takes its lowest value around $L_{q}\sim 1/M_{\text{KK}}$.

This behavior of $E_{\rm chaos}$ is naturally understood, because in the limit $L_{q}\to\infty$ or $L_{q}\to 0$, the system is expected to reduce to an integrable model. For example, in the former limit $L_{q}\to\infty$, the string is straight and resides at the bottom of the D4-soliton geometry, and the coefficients of the interaction $V(t)$ is suppressed as $1/L_{q}$, therefore the chaos disappears. In fact, in this limit $L_{q}\gg 1/M_{\rm KK}$ the chaos energy lower bound is calculated as

[TABLE]

and it diverges as $\sim L_{q}^{3}$.

The lessons from the toy model are that the chaos should appear at an energy above some nonzero value, and that the energy lower bound for the chaos diverges as $L_{q}\to\infty$. We will find that these are exactly seen in the full numerical simulations presented in the next section.

IV Chaos of interquark force

In this section, we explore the full dynamical motion of string by numerical simulation to examine the chaotic motion. For this purpose, we employ the numerical techniques to study dynamical string developed in Ishii and Murata (2015, 2016). The detail of numerical calculations are summarized in appendix A. To induce the nonlinear dynamics of the string, we instantly move the position of the string endpoints at the boundary of the D4-soliton geometry. In the gauge theory, this corresponds to an instantly forced motion of the quarks: a small deformation of the Wilson loop along the time direction. This produces a nonlinear dynamics of the gluon flux tube induced by the motion of quark and antiquark pair, afterwards.

To perform the numerical calculation, we again employ the $r$ coordinate, so the metric is (6). Now, as a world-sheet coordinate system, we take double null coordinates $(u,v)$ and specify the string configuration by $t=T(u,v),\,r=R(u,v),\,\vec{x}=\vec{X}(u,v)$. The condition on the induced metric $h_{ab}$ for $u,v$ to be null coordinates is given by $h_{uu}=h_{vv}=0$. Then, we have $-h=h_{uv}^{2}-h_{uu}h_{vv}=h_{uv}^{2}$. Thus in the double null coordinate, the Lagrangian for the string is proportional to $h_{uv}$, and working in the unit $M_{\text{KK}}=3/2$, the NG action becomes

[TABLE]

From this action, we obtain the evolution equations of the string:

[TABLE]

The double null conditions give constraints

[TABLE]

They are conserved by time evolution: $\partial_{v}(\cos^{3}RC_{u})=\partial_{u}(\cos^{3}RC_{v})=0$. We impose them at the initial surface and time-like boundaries of the string worldsheet and solve time evlution based on Eqs.(23-25). (The numerical technique for solving evolution equations are summarized in appendix A.2.)

Using the residual coordinate transformations, $u\to F(u),\,v\to G(v)$, we put the boundaries of the worldsheet at $u-v=0$ and $u-v=\pi$. As an initial condition, we take a static string configuration obtained in Sec. II. (See appendix A.1 for the detail of the numerical construction of the initial data.) In the unit $M_{\text{KK}}=3/2$, static string configurations form a one-parameter family of initial conditions for $r_{\text{center}}$, where $r_{\text{center}}$ denotes the initial $r$-coordinate at the tip of the string. This $r_{\text{center}}$ is one-to-one correspondent to the interquark distance $L_{q}$. Here we use the static solution with the initial condition $r_{\text{center}}=0.2$ (corresponding to $L_{q}=2.884$) to demonstrate the simulation of the dynamics.

To induce the nonlinear dynamics of the string, we impose a time-dependent boundary condition on string endpoints. Introducing time and spatial coordinates on the worldsheet as $\tau=u+v$ and $\sigma=u-v$, we consider the following forced motion (“quench”) of the string endpoints along the $X_{1}$ direction:

[TABLE]

where $\alpha(\tau;\Delta\tau)$ is defined by

[TABLE]

In our numerical simulation, we introduce a small cutoff $\delta$ near string endpoints and set our numerical domain in $\delta\leq\sigma\leq\pi-\delta$. As long as we take a small enough cutoff $\delta$, it may not matter to our results as shown in appendix A.3.333 This cutoff might be understood as an integrable deformation of dual field theories, see Smirnov and Zamolodchikov (2017); Cavaglia et al. (2016); McGough et al. (2018); Kraus et al. (2018); Chakraborty (2018). There are two parameters $\epsilon$ and $\Delta\tau$, which are the amplitude and the time scale of the quench. One can check that $\alpha(\tau;\Delta\tau)$ is $C^{\infty}$ in all $\tau$ and has a compact support in the region $0\leq\tau\leq\Delta\tau$.

Boundary conditions for the other variables at time-like boundaries are $X_{2}=X_{3}=0$ and $R=R_{\textrm{ini}}$, where $R_{\textrm{ini}}\equiv R(\tau=0,\sigma=\delta)$ is the initial value of the $R$ at the boundary. Because of the trivial boundary conditions of $X_{2}$ and $X_{3}$, they are identically zero throughout time evolution. The boundary value of $T$ is determined by constraints $C_{u}=C_{v}=0$. (See appendix A.2 for the detail.) With these boundary conditions, the string motion is $\mathbb{Z}_{2}$-symmetric under $x^{1}\to-x^{1}$. Fig. 5 shows the string configuration in the $x^{1}$-$r$ plane for $\epsilon=0.4$ and $\Delta\tau=4$. In this figure, we took a time slice of bulk coordinate as $t=T(u,v)$. Fig. 5 contains snapshots of string at early and late times. At early times, string profiles seem smooth. On the other hand, at late times small spatial deformations are observed. We also monitored violation of constraints (26, 27) and found that they are sufficientlly small. (See appendix A.3).

To observe its chaos, we focus on the tip of the string, $r_{\text{center}}=R(\tau,\sigma=\pi/2)$, since the string motion is completely $\mathbb{Z}_{2}$-symmetric due to our boundary condition and the quench. Fig. 6 shows time dependence of $r_{\text{center}}(t)$. Chaos means a sensitivity to the change of the initial conditions. To explore the sensitivity of the string motion, we consider a linear perturbation: $T\to T+\delta T,R\to R+\delta R$ and $\vec{X}\to\vec{X}+\delta\vec{X}$. We numerically solve the linear evolution equations for $(\delta T,\delta R,\delta\vec{X})$ on the time dependent background $(T(u,v),R(u,v),X_{1}(u,v))$. Initial conditions are [math] for all variables and boundary conditions are $\delta X_{1}(\tau,\sigma=\delta,\pi-\delta)=\pm\alpha(\tau;\Delta\tau=4)$, and $(\delta R,\delta X_{2},\delta X_{3})|_{\sigma=0,\pi}=0$. The boundary conditions for $\delta T$ at $\sigma=\delta,\pi-\delta$ are again determined by linearized constraints $\delta h_{uu}=\delta h_{vv}=0$.

The results of time evolution of the linear perturbations are as follows. Fig. 7 shows a time evolution of $\delta r_{\text{center}}(t)$, the tip of the string (located at $X_{1}=0$ due to the $\mathbb{Z}_{2}$ symmetry), for $\epsilon=0.4$. The horizontal axis is the bulk time coordinate $t=T(\tau,\sigma=\pi/2)$. We can find an exponential growth of the initial perturbation, which implies chaos. Fitting the amplitude, we obtain the positive Lyapunov exponent as $\lambda_{L}\simeq 0.09$ .

In the numerical simulations, we also observed that for small enough $\epsilon$ chaos does not occur, which means that there may be a chaos threshold bound of $\epsilon$ for each initial configuration given by $L_{q}$ (corresponding to the initial condition for $r_{\text{center}}$). We investigate the bound by running the simulation for different values of the parameters: quark separation $L_{q}$ and $\epsilon$. Our final “phase diagram” of the chaos of the QCD string is shown in Fig. 8. For each fixed quark-antiquark separation $L_{q}$, we numerically solve the full dynamical motion of string with different amplitude of quench $\epsilon$ and study whether the motion is chaotic. Below the solid line in Fig. 8 the motion is regular, while above the line chaos appears.

The phase diagram (Fig. 8) shows the following two important behavior: First, there exists a lower bound for the magnitude of the boundary pulse force $\epsilon$, for the chaos to occur. Second, the bound is a function of the interquark distance $L_{q}$, and it grows as $L_{q}$ grows. This in particular means that long strings are less chaotic. The shape of the bound described in Fig. 8 appears to be consistent with what we obtained in the rectangular string toy model in the previous section, Fig. 4. We shall investigate more on this behavior of the chaos bound in the next section, by using a physical model, to locate the origin of the chaos.

Finally, let us provide a prediction about a QCD quantity. We can also observe positive Lyapunov exponents for an observable of the gauge theory. When the string endpoint does not move, the AdS/CFT tells us that the force acting on the quark and the antiquark in the gauge theory is given by

[TABLE]

where $\lambda$ is the ’t Hooft coupling appearing as an overall coefficient in (22). The derivation of (31) is given in App. B. Fig. 9 shows the time evolution of $\langle\delta\vec{F}(t)\rangle$ and it implies the sensitivity of the interquark force to an initial perturbation. $\langle\delta\vec{F}(t)\rangle$ grows exponentially and its Lyapunov exponent is consistent with that of $\delta r_{\text{center}}(t)$. We find chaos of the interquark force via the AdS/CFT: the force in large $N_{c}$ pure Yang-Mills theory is generically sensitive to initial perturbations.

V Possible origin of the chaos

We found in the numerical simulation that the NG string in the confining geometry shows chaos, when the energy of the string exceeds some lower bound which is a function of the interquark distance $L_{q}$. In this section we argue why this behavior appears, based on a simple argument.

First of all, in the rectangular string model in Sec. III, the chaos shows up as a result of the different oscillation modes at the bottom of the string. On the other hand, it is known that straight string is integrable. This leads us to suspect that the origin of the chaos should be at the boundaries of the string. Let us consider a string model shown in Fig. 10: Two quarks are connected by a long straight string in the three-dimensional space. On the string any wave can propagate, and the motion is integrable. The wave will hit the boundary which is a quark. The boundary is not a point, but a region of the QCD scale. Since the string propagation part is integrable, any chaos, if exists, should originate in the boundary regions. We naively assume that when the magnitude of a wave hitting the boundary region exceeds some threshold value $\epsilon_{0}$ the chaos emerges. The wave amplitude will decay while it propagates, and so, the system with a larger interquark distance $L_{q}$ is expected to be less chaotic.

To quantify this physical model, we solve a motion of the wave propagating on the straight string. If the NG string sits at the bottom of the geometry, the fluctuation of the string obeys the wave equation

[TABLE]

The mass can be obtained by the analysis of the fluctuation of the straight string. A typical solution with a momentum $k_{0}$ larger than the mass scale, $k_{0}\gg M_{\rm KK}$, is

[TABLE]

where $f(k)$ is centered at $k=k_{0}$. Expanding this for small $M_{\rm KK}/k_{0}$, we obtain

[TABLE]

which means that the amplitude of the fluctuation decays along the propagation on the string. The timescale for the fluctuation to reach the other side of the string is estimated as $t\sim L-L_{0}$ where $L_{0}/2$ is the size of the boundary region which is expected to be the QCD scale. The typical momentum $k_{0}$ is estimated as $k_{0}\sim\pi/(2\Delta\tau)$ for the initial kick in our numerical simulation. Using these, the lower bound for the chaos is given by

[TABLE]

This expression shows that a larger $L_{q}$ makes the chaos diminished.

By this analytic expression (35), we can fit our numerical lower bound of the chaos, Fig. 8. Our numerical simulation uses $M_{\text{KK}}=3/2$, $\Delta\tau=4$. We find that choosing $\epsilon_{0}=0.064$ and $L_{0}=1.2$ fits the numerically obtained bound qualitatively, see Fig. 11. The obtained value, $L_{0}/2\sim 0.6$, roughly coincides with $1/M_{\rm KK}$ which is the QCD scale of the model.

From this argument, we find that a physical picture consistent with the results of the numerical simulation is a quark model in which an integrable string connect two boundaries whose size is of the QCD scale, and the boundary region produces chaos if the input wave exceed a certain threshold amplitude. The chaos originates in the boundaries of the QCD string, the constituent quarks.

VI Summary and Discussion

In this paper, we studied chaos and time evolution of interquark force by using the AdS/CFT correspondence. We performed a full nonlinear numerical simulation of the dynamics of a NG string in the confining geometry in the gravity side. The AdS/CFT translates the chaos of the NG string to the chaos of the interquark force. We found that the interquark force in large-$N_{c}$ four-dimensional pure Yang-Mills theory is generically sensitive to initial perturbations, and it is actually chaotic.

Our numerical calculation of the string in the D4-soliton background enabled us to analyze the full dynamical motion in details, and the Lyapunov exponent was obtained. Using the AdS/CFT dictionary, we further obtained the Lyapunov exponent of the interquark force. Normally, time-dependence of gauge-invariant non-perturbative observables of QCD is quite difficult to compute, thus, our results provide a theoretical prediction: the dynamics of the non-perturbative Yang-Mills gauge theory may be generically chaotic.

Our numerical simulations have two adjustable parameters: the interquark distance $L_{q}$ and the strength $\epsilon$ of the impulse force on the quarks to make them start moving. By area-bombing the parameter space, we obtained a phase diagram of the chaos, Fig. 8. It exhibits a unique picture: there exists a lower bound of $\epsilon$ for the chaos to occur, and the bound grows as $L_{q}$. This feature can be understood if the chaos originates in the constituent quark sectors (which are the boundaries of the QCD string), as provided in Sec. V with a simple model.

We provided a prediction of the Lyapunov exponent for the interquark forces. We hope we can confirm the exponent by some other direct calculations of QCD. Recently, the gradient flow techniques have been applied to lattice QCD simulations and the energy-momentum tensor on the lattice was defined through a flow equation Suzuki (2013). By using these techniques, the three dimensional distribution of energy-momentum stress tensor in $SU(3)$ gauge theory is non-perturbatively computed Yanagihara et al. (2018). However, the lattice QCD analyses are still only for static observables, and it is difficult to follow the time dependence. Nevertheless, it would be beneficial to compare the structure of the lattice QCD string with the holographic QCD string and find some difference, to locate possible origin of chaos qualitatively.

Our study focused on light modes of the large $N_{c}$ QCD, which are mesons and glueballs, while heavy nonlocal excitations exist: baryons and nuclear resonances. It would be important to quantify chaos of large $N_{c}$ baryons and nuclei and compare them with that of mesons and the QCD strings to find any difference in origin. Again, holography can help analyzing the chaos of the single or multiple baryon(s). They are known to be dual to D-branes called baryon vertices Witten (1998b) in the gravity side, so the motion of the baryons are well-approximated by a dimensionally reduced Yang-Mills theories Hashimoto (2009). Based on the classical 1-dimensional Yang-Mills analyses Matinyan et al. (1981, 1986) and on their D0-brane interpretation Gur-Ari et al. (2016); Asano et al. (2015b), or more detailed ADHM-like matrix model formulation Hashimoto et al. (2010) and its quantum states Hashimoto et al. (2019), it is possible to quantify the chaos of baryons. Since it is known that nuclear resonances follow quantum chaos Haq et al. (1982), finding out random matrix-like behavior from the classical holographic baryons would be interesting.

The chaos in the gravity side has been studied in the context of black hole horizons and the infinite redshift. The universal chaos bound discovered in Maldacena et al. (2016) is $\lambda\leq 2\pi T$ for large $N_{c}$ system with a finite temperature $T$, and it is proven that all observables in the large $N_{c}$ limit should obey this chaos bound for the quantum Lyapunov exponent defined by the out-of-time ordered correlators. Our case is at zero temperature, so, if we naively apply the chaos bound to the zero-temperature large $N_{c}$ QCD, any chaos is not allowed. This appears to contradict with our finding that the interquark force has a nonzero Lyapunov exponent and thus is chaotic — apparently there should be a loophole. The point is that the bound in Maldacena et al. (2016) was for local operators, while our observables are non-local, so the bound does not apply naively. Since non-Abelian gauge theories are always accompanied by non-local observables, it would be interesting to study how the quantum Lyapunov exponent of those non-local observables in generic gauge theories is theoretically observed, and how they play a role in determining the spectral/dynamical aspects of generic gauge theories.

Acknowledgements.

We would like to thank Tadakatsu Sakai, Motoi Tachibana, and Ryosuke Yanagihara for variable discussions. The work of K.H. was supported in part by JSPS KAKENHI Grants No. JP15H03658, No. JP15K13483, and No. JP17H06462.

Appendix A Numerical details

A.1 Initial data

As the initial data, we use the static string configuration. Here, we explain how to express the static solution in the double null coordinate $(u,v)$. Introducing $\tau=u+v$ and $\sigma=u-v$, we assume that the static solution is written as $T=\tau$, $\vec{X}=(X(\sigma),0,0)$ and $R=R(\sigma)$. In this assumption, Eq.(23) is automatically satisfied. Integrating Eq.(25) by $\sigma$, we obtain

[TABLE]

where $r_{\rm center}$ is the integration constant and ${}^{\prime}\equiv d/d\sigma$. Substituting above expression into constraints (26) and (27), we have

[TABLE]

At $R=r_{\rm center}$, we have $R^{\prime}(\sigma)=0$. Thus, $R=r_{\rm center}$ corresponds to the position of the tip of the hanging string. Note that this equation is regular at $R=\pi/2$ and well-behaved near the AdS boundary. On the other hand, near the tip of the hanging string, $R^{\prime}\sim\sqrt{r_{\rm center}-R}$. This is not a suitable form for the numerical integration around $R=r_{\rm center}$. From Eq.(24), we can obtain the other equation for $R(\sigma)$ as

[TABLE]

This can also be derived by differenciating Eq.(37) by $\sigma$. Contrary to Eq.(37), above equation is singular at $R=\pi/2$ but regular at $R=r_{\rm center}$. Therefore, in our numerical construction of the initial data, we integrate Eq.(38) from $R=r_{\rm center}$ to $R=(\pi/2+r_{\rm center})/2$. We then switch the equation to Eq.(37) and continue the integration from $R=(\pi/2+r_{\rm center})/2$ to $R=\pi/2$. Once we have the numerical solution of $R(\sigma)$, we also obtain $X(\sigma)$ integrating Eq.(36). As the result, we have the right half of the static string in Fig.2. We reparametrize the worldsheet coordinates as $\tau\to c\tau$ and $\sigma\to-c\sigma+c^{\prime}$ ($c$ and $c^{\prime}$ are constants) so that $R|_{\sigma=0}=\pi/2$ and $R|_{\sigma=\pi/2}=r_{\rm center}$ are satisfied. The left half of the static string can be easily generated by the $\mathbb{Z}_{2}$-symmetry: $R(\sigma)=R(\pi-\sigma)$ and $X(\sigma)=-X(\pi-\sigma)$. Then, time-like boundaries of the worldsheet are located at $\sigma=0,\pi$.

A.2 Time evolution

The original form of evolution equations (23-25) is numerically unstable. To stabilize time evolution, we eliminate $T_{,u}$ and $T_{,v}$ from Eqs.(23-25) using the constraint equations (26) and (27). Resultant equations are written in the form of

[TABLE]

where $\Phi=(T,R,\vec{X})$, $\hat{\Phi}=(R,\vec{X})$ and $\bm{F}$ is a non-linear function of its arguments.

We take uniform grid along $u$ and $v$ as in Fig.12. The grid points are explicitly written as $v=jh$ and $u=(i+j)h+\delta$ ($i=0,1,2,\cdots,N$, $j=0,1,2,\cdots$.) where $h=(\pi-2\delta)/N$ is the mesh size and $N$ is the number of grid points along the $u$-direction. Our numerical domain is in $\delta\leq u-v\leq\pi-\delta$ and $v\geq 0$. We introduced a small cutoff $\delta$ near the time-like boundaries of the worldsheet. If we set exactly $\delta=0$, the numerical simulation immediately breaks down and we cannot even see regular time evolutions.

Let us focus on points N, E, W, S and C in Fig.12. We can evaluate $\Phi$ and its derivatives at the point C with second-order accuracy in $h$ as $\Phi_{,uv}|_{\textrm{C}}=(\Phi_{\textrm{N}}-\Phi_{\textrm{E}}-\Phi_{\textrm{W}}+\Phi_{\textrm{S}})/h^{2}$, $\Phi_{,u}|_{\textrm{C}}=(\Phi_{\textrm{N}}-\Phi_{\textrm{E}}+\Phi_{\textrm{W}}-\Phi_{\textrm{S}})/(2h)$, $\Phi_{,v}|_{\textrm{C}}=(\Phi_{\textrm{N}}+\Phi_{\textrm{E}}-\Phi_{\textrm{W}}-\Phi_{\textrm{S}})/(2h)$ and $\Phi|_{\textrm{C}}=(\Phi_{\textrm{E}}+\Phi_{\textrm{W}})/2$, where $\Phi_{\textrm{N,E,W,S}}$ denote numerical values of $\Phi$ at points N,E,W,S. Substituting them into Eq.(39), we obtain the discretized version of the evolution equation. The equation determines $\Phi_{\textrm{N}}$ from $\Phi_{\textrm{E,W,S}}$. We use the Newton-Raphson method for solving the equation.

Once we have the initial data at $v=0$ and boundary data at $u-v=\delta,\pi-\delta$, we can determine the solution in our numerical domain by solving the discretized equation. As the initial data, we use the static string obtained in Sec.A.1. (So, the constraint (26) is satisfied at $v=0$.) At boundaries $u-v=\delta,\pi-\delta$, we do not change $R$ from its initial value: $R(\tau,\sigma=\delta)=R(\tau=0,\sigma=\delta)$ and $R(\tau,\sigma=\pi-\delta)=R(\tau=0,\sigma=\pi-\delta)$. We impose the Dirichlet conditions for $X_{1}$ as in Eqs.(28) and (29). To determine the boundary value of $T$, we consider points P, Q, and R in Fig.12. We can evaluate $\Phi$ and its $v$-derivatives at the point R as $\Phi_{,v}|_{\textrm{R}}=(\Phi_{\textrm{P}}-\Phi_{\textrm{Q}})/h$ and $\Phi|_{\textrm{R}}=(\Phi_{\textrm{P}}+\Phi_{\textrm{Q}})/2$. Substituting them into the constraint equation (27), we have the equation for $T_{\textrm{P}}$. By the similar way, using the other constraint (26), we obtain the left boundary value of $T$.

Substituting $\Phi\to\Phi+\delta\Phi$ into Eq.(39) and taking first order in $\delta\Phi$, we obtain the linear partial differential equation for $\delta\Phi=(\delta T,\delta R,\delta X)$. We also solve the evolution of the linear perturbation numerically. Its numerical procedure is completely parallel to that for the background.

A.3 Error analysis

As the measure of the numerical error, we monitior the violation of the constraints (26) and (27). We introduce the normalized constraint as

[TABLE]

where $\mathcal{N}_{u}$ and $\mathcal{N}_{v}$ are “scales” of constraints:

[TABLE]

We also add $1$ to the denominater of Eq.(40) for the case of $\mathcal{N}_{u}\simeq\mathcal{N}_{v}\simeq 0$. We further introduce the one dimensional function $C_{\textrm{max}}(v)$, which masures of the constraint violation on the fixed $v$-slice as

[TABLE]

Fig.13 shows $C_{\textrm{max}}(v)$ for $N=4000,8000,16000$. (The numerical integration by $N=4000$ broke down at $v\simeq 30$.) We considered the same setup as Fig.6. The cutoff near time-like boundaries is fixed as $\delta=0.01$. The constraint violation keeps small value ($C\textrm{max}\lesssim 10^{-3}$ for $N\gtrsim 8000$). We can also see $C_{\textrm{max}}\propto 1/N^{2}$. This is consistent with the fact that our numerical scheme has second order acuracy.

In Fig.13, we show the time dependence of the tip of the hanging string $r_{\textrm{center}}(t)$ for several values of the cutoff: $\delta=0.005,0.01,0.02$. The number of grid points are fixed as $N=8000$. We again considered the same setup as Fig.6. The dependence on $\delta$ is small and typical chaotic bahaviour of the string does not depend on the value of $\delta$. Based on the error analysis here, we show results in the main text of this paper for $N=8000$ and $\delta=0.01$.

Appendix B Interquark force from holography

B.1 Derivation

Here we derive the formula (31) giving a relation between the force acting on quarks in the gauge theory and the NG string in the gravity side. We follow the argument given in App. D of Ishii and Murata (2015).

We write the on-shell NG action as

[TABLE]

where $\tilde{X}$ is a solution of the equation of motion with the boundary condition $X(t,r\to\pi/2)=\vec{x}_{q}(t),\vec{x}_{\bar{q}}(t)$. The force acting on quarks in the gauge theory is holographically given by Ishii and Murata (2015)

[TABLE]

where $m$ is mass of the quark, $\vec{v}=\dot{\vec{x}}_{q}$ and $\gamma=(1-\vec{v}^{2})^{-1/2}$.

Now, let us evaluate $\delta S/\delta\vec{x}_{q}$ in the gravity side. Since the background metric is now given by (6), in the static gauge the NG action becomes

[TABLE]

where $\dot{}=\partial_{t}$ and ${}^{\prime}=\partial_{r}$. In what follows we always take our unit $M_{\text{KK}}=3/2$. Solving the equation of motion for $\vec{X}$ near the D4-boundary; $r=\pi/2$, we obtain an asymptotic expansion form of the solution as

[TABLE]

where $\vec{a}=\ddot{\vec{x}}_{q}$ and we have defined $\epsilon=\pi/2-r$.

To obtain the force, let us consider the variation of the action (46),

[TABLE]

where $\mathcal{L}$ is the integrand of (46) and we introduce a cutoff at $r=\pi/2-\epsilon$. Substituting the asymptotic solution (47) into (48), we obtain

[TABLE]

The last term $A$ involves complicated terms, but when we consider a probe approximation $\dot{\vec{x}}_{q}\to 0$, $A$ actually vanishes, so we do not care about $A$. From this expression we find that the quark mass $m$ corresponds to $\lambda/6\pi\epsilon^{2}$, which is divergent when $\epsilon\to 0$. Setting $m=\lambda/6\pi\epsilon^{2}$ and considering probe approximation $\dot{\vec{x}}_{q}\to 0$, we get the force acting on the quark

[TABLE]

where we have replaced $\vec{f}_{4}$ with $\partial_{r}^{4}\vec{X}/4!$.

B.2 Sensitivity of the interquark force

To numerically compute the sensitivity of the interquark force to initial perturbations, we can employ two procedures to do it. One is a direct calculation: In Eq. (50), change the worldsheet coordinate to double null $u-v$ coordinate and consider linear perturbations $\Phi\to\Phi+\delta\Phi$. Evaluating the linearized differential equation for $\delta\Phi=(\delta T,\delta r,\delta X)$ at $r=\pi/2$, we obtain the sensitivity of the interquark force $\langle\delta\vec{F}\rangle$. However, this is a little tough since the right hand side of Eq. (50) has four derivatives of $r$. So, we employ the other one which we explain in the following.

In Eq. (48), The integrand is explicitly written as

[TABLE]

Substituting a solution $\tilde{X}$ to the right hand side and evaluating it at $r=\pi/2$, this quantity is corresponding to infinitely heavy quarkmass and the interquark force. Changing to the double null coordinate,

[TABLE]

and plugging them into (51), we find that this can be more easily evaluated numerically since the right hand side includes only single derivatives. Considering linear perturbations $\Phi\to\Phi+\delta\Phi$ and evaluating it at $u=0$ or $u=\pi$, we obtain the sensitivity to perturbations of the interquark force:

[TABLE]