Non-linear perturbation of black branes at large $D$

The Einstein equations describing the black-brane dynamics both in Minkowski and AdS background were recently recast in the form of coupled diffusion equations in the large-$D$(imension) limit. Using such results in the literature, we formulate a higher-order perturbation theory of black branes in time domain and present the general form of solutions for arbitrary initial conditions. For illustrative purposes, the solutions up to the first or second order are explicitly written down for several kind of initial conditions, such as a Gaussian wave packet, shock wave, and rather general superposed sinusoidal waves. These could be the first examples describing the non-trivial evolution of black-brane horizons in time domain. In particular, we learn some interesting aspects of black-brane dynamics such as the Gregory-Laflamme (GL) instability and non-equilibrium steady state (NESS). The formalism presented here would be applicable to the analysis of various black branes and their holographically dual field theories.


Introduction
It is important to understand the dynamics of higher-dimensional black objects, since it tells us much about the nature of higher-dimensional gravitational theories and their holographically dual quantum field theories. The strong non-linearity of gravity, however, usually prevents us from understanding the dynamical properties of black objects beyond the linear-perturbation regime without highly sophisticated skills of numerical computation. The Gregory-Laflamme (GL) instability [1], which is a universal instability of higher-dimensional black objects, is a good example to see the above situation. Though the GL instability in the nonlinear regime is quite interesting, its analysis needs sophisticated skills of numerical relativity [2]. While there exists a semi-analytic higher-order perturbation method [3], it seems applicable only to static problems.
Recently, Emparan, Suzuki, and Tanabe showed that the Einstein equations describing the horizon dynamics of black branes in both Minkowski and Anti-de Sitter (AdS) background are recast in the form of coupled non-linear diffusion-type equations when the number of spatial dimensions is large [4]. This result provides us with a unique approach to the non-linear dynamics of black objects in higher dimensions. The authors indeed showed that the unstable black strings converge to nonuniform black strings (NUBSs), which had been predicted to happen above a critical dimension [5], by solving the diffusion equations numerically with a few lines of Mathematica code. It is added that the blackfold approach [6] is also thought to serve as a powerful approach to analyze the evolution of GL instability.
Once the simple diffusion equations were obtained [4], it is natural to ask if the non-linear properties of black-brane dynamics can be understood analytically. In this paper, we develop a systematic non-linear perturbation theory of asymptotically flat and AdS black branes, allowing the perturbations to be dynamical. Using the Fourier and Laplace transformation to solve the partial differential equations (PDEs), the perturbation equations are solved order by order for given arbitrary initial conditions up to the integration associated with the inverse transformation.
While the formulation is so general that it would be applicable to various problems, we pick up several examples as the initial conditions, which are a Gaussian wave packet, a step-function like shock configuration, and quite general discretely superposed sinusoidal waves. For these examples, the integration associated with the inverse transformation is completed up to the first or second order, and the properties of solutions are examined. Through these examples, one will see the validity of formalism itself and some unknown, or yet-to-be-confirmed, non-linear properties of black-brane dynamics. In particular, in the case of asymptotically flat black branes, an interesting non-liner property of GL instability resulting from the mode-mode coupling is unveiled at the second order. In the case of shock propagation on asymptotically AdS black branes, the analytic description of non-equilibrium steady state (NESS), which was recently discussed in the Riemann problem of relativistic fluid mechanics and field theories [7], is presented. This paper is organized as follows. In Sec. 2, the asymptotically flat black branes are investigated. In Sec. 2.1, we present the perturbation equations for asymptotically flat black branes and their general form of solutions. In Sec. 2.2, we apply the general result to the Gaussian wave packet. In Sec. 2.3, we consider the discretely superposed sinusoidal waves. In Sec. 3, we consider the nonlinear perturbation of asymptotically AdS black branes. Here, the formulation and applications are presented in parallel with Sec. 2, but a new example of initial condition, the step-function like shock, is investigated in Sec. 3.3. Section 4 is devoted to conclusion. Throughout this paper, we follow the notations in Ref. [4].

Asymptotically flat black branes 2.1 Perturbation equations and general form of solutions
In the large-D(imension) approach, the horizon dynamics of vacuum black branes without a cosmological constant are described by two functions, m(t, z) and p(t, z), where t is time and z is the spatial coordinate along which the horizon extends [4]. m and p represent the mass and momentum distributions along the horizon, respectively. m → +0 corresponds to the pinching off of the horizon. The equations of motion for these quantities take form of coupled non-linear diffusion equations, A uniform black-brane solution corresponds to m(t, z) ≡ 1 and p(t, z) ≡ 0. Since we are interested in the dynamical deformation of such a uniform solution, we introduce one-parameter families of m(t, z) and p(t, z), and expand them around the uniform black-brane solution, where ǫ is a constant parameterizing the families. Substituting these expansions into Eqs. (2.1) and (2.2), we obtain the equations of motion at O(ǫ ℓ ) (ℓ ∈ N), where the dot and prime denote the derivatives with respect to t and z, respectively. The righthand side of Eq. (2.7), ψ ℓ (t, z), which we call a source term, is a polynomial of the lower-order perturbations and their first spatial derivatives, For example, the source terms for ℓ = 2 and ℓ = 3 are given by In the rest of this section, we are looking for the general form of solutions to the perturbation equations (2.6) and (2.7), combining the Fourier and Laplace transformations (see, e.g., [8]). A similar technique is found to be used in Ref. [9,10] to analyze the higher-order perturbation of surface-diffusion equation, which is a single non-linear PDE.
Before starting to solve Eqs. (2.6) and (2.7), let us introduce the notations associated with the Fourier and Laplace transformations. For a given function, say f (t, z), we shall denote its Fourier transformation with respect to z byf (t, k), and its Laplace transformation with respect to t by the corresponding capital letter F (s, z). Namely, Then, a capital letter with a bar denotes a Fourier-Laplace transformation as (2.14) In addition, we define two kind of convolutions, With the notations introduced above, the Fourier-Laplace transformed version of Eqs. (2.6) and (2.7) are written as coupled algebraic equations in a matrix form .
By simple algebra, the inverse matrix A −1 is found to be decomposed into two parts, 20) . Here, we have used Equations (2.23) and (2.24) are exactly what we wanted, namely, the general form of solutions to the perturbation equations (2.6) and (2.7) for given arbitrary initial conditions, m ℓ (0, z) and p ℓ (0, z) (ℓ ∈ N). It depends on the problem which expression, (2.23) or (2.24), is easer to compute.
For all examples considered in this paper, Eq. (2.23) seems easer to compute. Using Eq. (2.22), one can easily obtain which is useful when one uses expression (2.24). In principle, one can obtain the arbitrary-order solutions, m ℓ (t, z) and p ℓ (t, z), by computing the right-hand side of Eq. (2.23) or (2.24) order by order. However, as the source term ψ ℓ (t, z) becomes complicated as ℓ increases, it is fare to say that to obtain the solutions analytically until arbitrary order is impossible in general. In addition, when the initial condition is a complicated function, even the first-order solution can be impossible to obtain analytically. Namely, the inverse Fourier transformation in Eq. (2.23) cannot be performed analytically in such a case. In the rest of this section, we shall consider two examples of initial conditions, for which the first few-order solutions are analytically obtainable.

Gaussian wave packet
As the first example, we adopt the Gaussian wave packet as the initial perturbation given to the asymptotically flat black brane. For this perturbation, the inverse Fourier transformation at the first order in Eqs. (2.23) and (2.24) can be computed analytically. Since the unperturbed black brane is unstable, the perturbation of course grows unboundedly and nothing unexpected happens in this sense. However, one can see how to use the general results obtained in the previous section and the validity of the method. In particular, comparing the perturbative solution with a full numerical solution, the perturbative solution turns out to be effective even for finite-amplitude dynamics. In other words, the convergence of ǫ-expansion (2.4) and (2.5) is rapid enough at least for this example.
Let us assume that the situation where the black brane is given O(ǫ)-perturbation taking form of a Gaussian wave packet, where β, b (> 0), and z 0 are real constants. Obviously, b and z 0 parameterize how much the wave packet spatially extends and the central position of the wave packet, respectively. β is just a normalization constant to give ∞ −∞ m 1 (0, z)dz = β. The initial perturbation of p 1 is given by the spatial derivative of m 1 for simplicity, though their initial conditions can be independent in nature.
The Fourier transformation of the above initial conditions arē (2.29)
One can compute such quantities as Substituting these results into Eq. (2.23), we obtain the first-order solutions We believe that this is the first example describing the non-trivial evolution of GL instability analytically in time domain, which is realized by virtue of the large-D method and perturbation theory developed in this paper.
Like the diffusion phenomenon of a Gaussian wave packet according to an ordinary diffusion equation, solutions (2.32) have temporal decay factor, which behaves as 1 √ b 2 +2t , and spatial de- The sinusoidal parts represent spatial oscillation with time-dependent wavelength 2π(b 2 +2t) t , which interestingly asymptotes to a universal value 4π as t → +∞. What crucially different from the ordinary diffusion is that the solutions temporally grow exponentially due to factor exp[ t 2 2(b 2 +2t) ]. It is stressed that this exponential growing happens eventually irrespective of b, which characterizes the extension of the initial wave packet. Substituting z = z 0 into Eq. (2.32), we can see the time dependence of the peak height, to diverge eventually. Thus, while the GL instability is generally said to be a long-wavelength instability, the initial perturbation taking form of a Gaussian wave packet necessarily grows exponentially however the 'scale' of perturbation b is small. The reason is that the Fourier spectrum of any Gaussian wave packet necessarily contains the GL mode k ∈ (−1, 1) \ {0} (see Sec. 2.3.1). This is quite reasonable but might be a somewhat interesting point.
Three-dimensional plots of 1 + m 1 (t, z) and p 1 (t, z) are presented in Figs. 2(a) and 2(b), respectively. One can observe the growth and oscillation described above. In addition, snapshots of 1 + m 1 (t, z) and p 1 (t, z) at selected moments, compared with numerical solutions, are presented in Figs. 2(c) and 2(d), respectively. Compared with the full numerical solutions, which are obtained by directly solving original equations (2.1) and (2.2) with the same initial conditions, i.e., m(0, z) = 1 + m 1 (0, z) and p(0, z) = p 1 (0, z), one can observe that the first-order solutions almost completely capture the qualitative features of the full solution during the time domain considered. Note that the deviation from the full solution, however, becomes large as the time proceeds, which results in the divergence of m and p.

Notes on second-order perturbation
The comparison between the first-order solution and full solution above tells us that O(ǫ 2 ) perturbations are negligible during the amplitudes of m and p are O(1) in the current example despite the ordinary expectation that the perturbation becomes invalid for such a large amplitude. We will see in Sec. 3 that O(ǫ) approximation is more accurate for the Gaussian perturbation to the asymptotically AdS brane than the present case.
For the aim to see the non-linear effects at the second order, it is natural to assume that the initial perturbation at the second order vanishes, m 2 (0, z) = p 2 (0, z) = 0. The reason is that m 2 (t, z) and p 2 (t, z) are composed of two independent parts as seen in Eq. (2.23): one is the contribution from initial perturbation m 2 (0, z) and p 2 (0, z), and the other is that from the source term ψ 2 (t, z). The former clearly has the same time dependence as the first term from Eq. (2.23). Namely, if we prepare the Gaussian wave packet as the initial condition of second-order perturbation, the secondorder perturbation evolves in the exactly same way as the O(ǫ) perturbation described above. Only the latter, the contribution from the source term, can have a different time dependence from the first-order perturbation. This will be seen explicitly in Sec. 2.3.
From the reason described above, we should assume that the initial perturbations vanish at O(ǫ 2 ), m 2 (0, z) = p 2 (0, z) = 0. Then, we see from Eq. (2.23) that what to compute at O(ǫ 2 ) is only the inverse Fourier transformation of the convolution between e sσ(k)t and the spectrum of source termψ 2 (t, k). Unfortunately, however, such a convolution in the present example involves the Gauss error function, not written in terms of elementary functions. Thus, it seems difficult to obtain the second-order solutions analytically, and therefore we stop the analysis on this example here.
and non-perturbative solution m(t, z) (resp. p(t, z)) obtained by numerically solving Eq. (2.1) and (2.2). The blue-dashed curve represents initial configuration m(0, z) (resp. p(0, z)). The green, red, and black solid curves represent the first-order solutions at t = 3.3, 6.7, and 10, respectively. The (green, red, and black) dashed curves represent the full numerical solutions at the corresponding time.

Superposed sinusoidal waves
As the second example, we consider the situation where the black brane is initially given an O(ǫ)perturbation being a superposition of an arbitrary number of sinusoidal waves. This example is simple but interesting enough to see what happens in the non-linear regimes. We set the following initial conditions where a n , k n , and ϕ n (n = 1, 2, · · · , N ) are real constants. It is noted that the right-hand side of Eq. (2.35) is not written as the general form of Fourier series expansion of a periodic function. However, choosing appropriate wave number k n and phase ϕ n , and taking summation over n from 0 to infinity, rather than from 1 to N , Eq. (2.35) can cover the Fourier series expansion of arbitrary piecewise continuous periodic function. The assumption that the second and higher order perturbations vanish initially are adopted from the same reason as the previous example in Sec. 2.2.
The Fourier transformations of the above initial configurations arē In the rest of this section, we shall compute the right-hand side of Eq. (2.23) order by order for these initial conditions.

First-order solutions
Since we have no source term at O(ǫ), ψ 1 ≡ 0, we see from Eq. (2.23) that what to compute is only the inverse Fourier transformation of the initial spectra,m 1 (0, k) andq 1 (0, k), multiplied by e sσ(k)t . These are easily computed to give where we have used s σ (−k) = s −σ (k). Substituting these results into Eq. (2.23), we obtain the first-order solutions, Since we are considering linear equations of motion, there is no mode-mode coupling appearing in non-linear regime, and therefore Eqs. (2.43) and (2.44) have simple interpretation. The factor of cos(k n z + ϕ n ) in Eq. (2.43) represents the time-dependent amplitude of the initially given mode cos(k n z + ϕ n ). Each mode evolves independently according to its growth or damping rate determined by e s + (kn)t and e s − (kn)t . From the concrete form of s ± (k) in Eq. (2.22), one can see that if k n ∈ (−1, 0) (resp. k n ∈ (0, 1)), such a mode grows exponentially due to e s − (kn)t (resp. e s + (kn)t ), which represents the GL instability.

Second-order solutions
Since we assume that the initial perturbations vanish at O(ǫ 2 ), m 2 (0, z) = p 2 (0, z) = 0, we see from Eq. (2.23) that what to compute is the inverse Fourier transformation of the convolution between e sσ(k)t and the Fourier spectrum of source termψ 2 (t, k). Using Eqs. (2.10) and (2.44), such a quantity is computed and written down in a simple form as by defining a function of time, a n a n ′ k n ′ [C a n a n ′ k n ′ [C Note that m(t, z) = 1 + m 1 (t, z) + m 2 (t, z) and p(t, z) = p 1 (t, z) + p 2 (t, z) with Eqs. Since the initial perturbations are assumed to vanish at O(ǫ 2 ), m 2 (0, z) = p 2 (0, z) = 0, the above O(ǫ 2 ) solutions contain only the contribution from the source term ψ 2 = −2p 1 p ′ 1 . If one prepares for non-vanishing initial conditions at the second order, its contribution is simply added to the above solution, but such a contribution will exhibit no interesting behavior since it has the time dependence similar to that of O(ǫ) solution.
In general, the multiple summation of any quantity with two indices n,n ′ T nn ′ can be decomposed as n,n ′ T nn ′ = n T nn + n<n ′ (T nn ′ + T n ′ n ). Here, n<n ′ represents the summation over all n and n ′ satisfying 1 ≤ n < n ′ ≤ N . Using this decomposition, one can decompose the multiple summation in Eqs. (2.47) and (2.48) as (2.50) The first term of Eqs. (2.49) and (2.50) represents the self-interference of each mode k n (n = 1, 2, · · · , N ). On the other hand, the second and third terms represent the interference between k n and k n ′ (n < n ′ ). The non-linear source term involves the mode-mode coupling, which is absent at the linear order. This coupling excites the terms of cos[(k n ± k n ′ )z] and sin[(k n ± k n ′ )z] in Eqs. respectively. Thus, the second-order solutions exhibit a variety of dispersion given by the exponents of quantities (2.51) and (2.52).

Notes on Gregory-Laflamme instability
Let us consider the meaning to investigate the higher-order perturbations from the stability point of view. The asymptotically flat black brane we consider here is essentially unstable. Namely, as seen in Sec. 2.3.1, if the initial perturbation contains any mode of which wave number k n ∈ (−1, 1) \ {0}, such a mode grows unboundedly. However, we consider the black brane in the large-D limit, namely, above the critical dimension [5]. Thus, the GL instability initially grows but it gradually damps in non-linear regimes, and eventually the horizon converges to non-uniform configuration [4]. It is pointed out that the second-order perturbation cannot stabilize the first-order instability because the first-order perturbation is treated as the fixed background, which appears as the source term, when we solve the second order. Nevertheless, one might expect if there appears the sign of the convergence to the non-uniform horizon at the second order. Although such a sign can appear at the second order, we cannot catch the sign in the result (2.49) and (2.50) unfortunately.
On the other hand, a black brane that is linearly stable can become unstable at the second order. If the initial perturbation does not contain any unstable mode, the initial perturbation will damp exponentially at linear level, as seen in the results of Sec. 2.3.1. However, the second-order perturbation involves various time dependence as seen in Eqs. (2.51) and (2.52). In order to see directly this situation, let us focus on a simple case as follow.
The above phenomenon is an interesting aspect of the GL instability, which was revealed for the first time by the present non-linear perturbation theory in time domain. It is intuitively understandable. The superposition of the two modes forms the beat. For simplicity, assume a 1 = a 2 ( = 0) in Eq. (2.35), then the superposed wave is written as This exhibits the fast spatial oscillation with the large wave number k 1 +k 2 2 which is enveloped by the slow oscillation with the small wave number k 1 −k 2 2 , which is called the beat phenomenon especially when the difference of the wave numbers is rather small k 1 − k 2 ≪ k 1 + k 2 . This slow oscillation is nothing but the origin of the GL instability at the second order. In Fig. 3, we present the three-dimensional plots of m(t, z) = 1 + m 1 (t, z) + m 2 (t, z) and p 1 (t, z) + p 2 (t, z) and their snapshots at selected moments. One can observe that the beat formed by the superposition of two modes at t = 0. As soon as the dynamics starts, such an initial wave rapidly damps as predicted by the first-order perturbation. As the time proceeds, however, the waves of which scale is the same order as that of the beat begin to grow and eventually diverge.

Asymptotically AdS black branes
In this section, we consider the non-linear perturbation of the asymptotically AdS black branes in the large-D limit. The governing equations of motion are almost the same as those in the asymptotically flat case except for a sign of one term. Thus, the formulation of the perturbation theory proceeds completely in parallel with the asymptotically flat case.
What different from the asymptotically flat black brane in the previous section is that the AdS black branes are stable: they are not suffered from the GL instability at least linearly. In addition, the gravitational phenomena in AdS background are interpreted as corresponding phenomena in the dual field theories via the AdS/CFT dictionary, and therefore have much more applications than the asymptotically flat case. In fact, we will apply the result of general perturbation theory to the problem of shock-wave propagation, which has been discussed in the context of AdS/CFT [7], in Secs. 3.3 and 3.4 in addition to the Gaussian wave packet and general superposed sinusoidal waves.

Perturbation equations and general form of solutions
For the asymptotically AdS neutral black branes in the large-D limit of general relativity, the mass and momentum distribution, m(t, z) and p(t, z), obey Eq. (2.1) and the following equation [4,7] with the domain (2.3), where the source term ψ ℓ is the same as those in the asymptotically Minkowski case, Eqs. (2.8)-(2.11).
Performing the Laplace and Fourier transformation in Eqs. (2.6) and (3.2), we obtain a couple of algebraic equations, which is written as Eq. (2.17) with A replaced by the following matrix, The inverse matrix D −1 is decomposed into two parts, After this decomposition, one can perform the inverse Laplace transformation L −1 in Eq. (2.19) with A replaced by D. Then, the general form of solution is given by .
While it depends on the chosen initial condition which expression between (3.7) and (3.8) is easier to compute, Eq. (3.7) is solely used in the rest of this paper. Using Eq. (3.6), the inverse Fourier transformation of e sσ(k)t is easily computed as This will be useful when one uses expression (3.8).

Gaussian wave packet
As in Sec. 2.2, we investigate the Gaussian wave packet as the initial perturbation given to the asymptotically AdS black brane. Compared with the full-order numerical solution, the linear-order approximation turns out to be rather accurate approximation in this case.

First-order solutions
For the initial perturbation given by the Gaussian wave packet, Eqs. (2.26) and (2.27), one can compute the following quantities Substituting these quantities into Eq. (3.7), we obtain the first-order solution, where cosh Ξ := b 2 + 3t Like the diffusion phenomenon of a Gaussian wave packet according to an ordinary diffusion equation, solutions (3.12) have temporal decay factor, which behaves as 1 √ b 2 +2t , and spatial decay factor exp[− (z−z 0 ) 2 2(b 2 +2t) ]. What crucially different from the ordinary diffusion is that the solution temporally damps rapidly due to factor exp[− t 2 2(b 2 +2t) ]. Three-dimensional plots of 1 + m 1 (t, z) and p 1 (t, z) are presented in Figs. 4(a) and 4(b), respectively. One can observe the fast damping described above. In addition, snapshots of 1+m 1 (t, z) and p 1 (t, z) at selected moments, compared with numerical solutions, are presented in Figs. 4(c) and 4(d), respectively. Compared with the full numerical solutions, which are obtained by directly solving original equations (2.1) and (3.1) with the same initial conditions, i.e., m(0, z) = 1 + m 1 (0, z) and p(0, z) = p 1 (0, z), one can observe that the first-order solution completely captures the full solution throughout the time domain considered. In other words, the higher-order perturbations are negligible, meaning that the ǫ-expansion Eqs. (2.4) and (2.5) converges rapidly for this example.

Shock wave
Here, let us consider the step-function like shock as the initial perturbation to the asymptotically AdS black brane. The propagation of this kind of shock is known as the Riemann problem in fluid mechanics. This classic problem attracts attentions recently in relativistic hydrodynamics since it makes us understand the non-equilibrium physics of quantum field theories. See the introduction of Ref. [7] for a brief but nice review for the recent development.
Assume that the black brane is given O(ǫ) perturbation as follows, where α is a real constant and sgn denotes the sign function, .

(3.15)
The assumption that p 1 initially vanishes is adopted to reproduce a situation considered in Ref. [7] (see the left panel of Fig. 5 in [7]). The Fourier transformation of the above initial conditions arē x 0 e −x 2 dx is the Gauss error function. Here, we have used the following formula (3.20) Using the fact that the Gauss error function is an odd function, one can immediately show that m 1 (t, −z) = −m 1 (t, z) and p 1 (t, −z) = p 1 (t, z) hold. Namely, m 1 (t, z) and p 1 (t, z) are spatially odd and even functions, respectively. In Fig. 5, three-dimensional plots of 1 + m 1 (t, z) and p 1 (t, z), where the parameter is chosen as α = −1/2, and their snapshots at selected moments are presented, compared with the nonperturbative numerical solutions obtained by directly solving Eqs. with α = −1/2, respectively. The comparison between first-order solution (c) 1 + m 1 (t, z) (resp. (d) p 1 (t, z)) and full-order numerical solution m(t, z) (resp. p(t, z)), obtained by solving Eqs. (2.1) and (3.1). The blue-dashed curve represents initial configuration m(0, z) (resp. p(0, z)). The green, red, and black solid curves represent the first-order solutions at t = 7.3, 14, and 22, respectively. The (green, red, and black) dashed curve represents the full numerical solution at the corresponding time.
condition similar to Eq. (3.14). Since the sign function is difficult to treat in numerical computations, we replace sgn(z) by tanh(cz) with a large positive number c (say, c = 50), based on the fact that lim c→+∞ tanh(cz) = sgn(z). In Fig. 5, one can see that as soon as the dynamics begins, the discontinuity separates into two shock fronts moving to left and right, where the former and latter are called the rarefaction wave and shock wave, respectively [7]. The point is that the rarefaction and shock is interpolated by an expanding plateau region with a non-zero constant flux, called the non-equilibrium steady state (NESS).
While it is interesting that we obtained the semi-analytic results (3.19) and (3.20) describing the shock propagation and NESS, they are not satisfactory from some viewpoints. In Fig. 5, one can observe that the full solution becomes asymmetric with respect to z → −z, though the linear solutions (3.19) and (3.20) continue to be symmetric. Thus, for example, the value of m at the NESS, which is always unity at O(ǫ), deviates between the first order and full solutions. These suggest that the higher-order perturbations are necessary to fill the above gaps. It seems impossible, however, to obtain the second-order solutions analytically since the first-order solution involves the error function. Thus, we will return to the Riemann problem in the next section as the example of the superposed sinusoidal waves.

Superposed sinusoidal waves
We consider the initial condition which is the superposition of sinusoidal waves like Eqs. (2.35)-(2.37). For later purpose, however, let us assume the O(ǫ) initial momentum vanishes which is assumed instead of Eq. (2.36). The first-and second-order solutions satisfying initial conditions (2.35)-(2.37) are presented in Appendix A.

First-order solutions
For the initial perturbation which is the superposition of sinusoidal waves (2.35) and (3.21), we obtain where we have used s σ (−k) = s −σ (k). Substituting these results into Eq. (3.7) and using the concrete form of the dispersion relation s σ (k), we obtain the first-order solutions, a n e −k 2 n t cos(k n t) cos(k n z + ϕ n ), (3.24) a n e −k 2 n t sin(k n t) sin(k n z + ϕ n ). (3.25) These results (3.24) and (3.25) tell us that the initial perturbation necessarily exhibits the damped oscillation for arbitrary non-zero wave number k n ∈ R \ {0} (n = 1, 2, · · · , N ), showing the black brane to be linearly stable.

Second-order solutions
Since we assume that the initial values vanish at O(ǫ 2 ), m 2 (0, z) = p 2 (0, z) = 0, we see from Eq. (3.7) that what to compute is the inverse Fourier transformation of the convolution between e sσ(k)t and the Fourier spectrum of source termψ 2 (t, k). Using Eqs. (2.10) and (3.25), such a quantity is computed and written down in a simple form as by defining a function of time, Substituting the above result (3.26) into Eq. (3.7), we have the second-order solutions as a n a n ′ k n ′ [F a n a n ′ k n ′ [F Note that m(t, z) = 1 + m 1 (t, z) + m 2 (t, z) and p(t, z) = p 1 (t, z) + p 2 (t, z) with Eqs. (3.24), (3.25), (3.28), and (3.29), represent O(ǫ 2 ) approximate time evolution of the initial perturbation, which takes the form of superposed sinusoidal waves (2.35), (3.21), and (2.37). These approximate solutions are rather general in the sense that initial condition (2.35) is general.
Since the initial perturbations are assumed to vanish at O(ǫ 2 ), m 2 (0, z) = p 2 (0, z) = 0, the above solutions contain only the contribution from the source term ψ 2 = −2p 1 p ′ 1 . If one prepares for non-vanishing initial conditions at the second order, its contribution is simply added to the above solution, but such a contribution will exhibit no interesting behavior since it has the time dependence similar to that in O(ǫ) solution.

Shock wave
We re-investigate the Riemann problem in Sec. 3.3, by choosing appropriate parameters (a n , k n , ϕ n ) in initial condition (2.35). In order to do so, the sign function, sgn(z), in initial condition (3.14) is replaced by the following function, where L is a positive constant and the periodic extension to entire R with period 2L is assumed.
Since this function has the following Fourier series expansion, The comparison between second-order solution (c) 1 + m 1 (t, z) + m 2 (t, z) (resp. (d) p 1 (t, z) + p 2 (t, z)) and full-order numerical solution m(t, z) (resp. p(t, z)), obtained by solving Eqs. (2.1) and (3.1). The blue-dashed curve represents initial configuration m(0, z) (resp. p(0, z)). The green, red, and black solid curves represent the second-order solutions at t = 83.3, 167, and 250, respectively. The (green, red, and black) dashed curve represents the full numerical solution at the corresponding time.
the parameters in initial condition (2.35) should be a n = 4 (2n − 1)π , k n = (2n − 1)π L , θ n = − π 2 (3.32) for n = 1, 2, · · · , N with the limit N → ∞. Taking large L and focusing on the spatial region around the center z = 0, there must be no difference in the dynamics during a finite interval of time between ones using sgn(z) and sgn L (z), due to the (non-relativistic) causality encoded in the equations of motion. In Figs. 6(a) and 6(b), we present the three-dimensional plots of full-order numerical solutions of m(t, z) and p(t, z), respectively, starting from the initial condition m(0, z) = 1 + m 1 (0, z) and p(0, z) = 0, where m 1 (0, z) is given by Eqs. Unlike the O(ǫ) approximation in Sec. 3.3, the O(ǫ 2 ) approximation presented here captures the spatially asymmetric features under the reflection z → −z of full solution. Furthermore, the error of values of m and p at the plateau of NESS between the O(ǫ 2 ) solution and full solutions are within a few percent (1.2% and 5.0% for m and p, respectively, for the above choice of parameters), which were inadequately large for the O(ǫ) approximation in Sec. 3.3. This is by virtue of the second-order solutions, which appropriately take into account the 'back-reaction' of the first-order solution through the source term.
Due to the Gibbs phenomenon, there appears the oscillation near the jump in the initial shock. Such oscillations stem from the artificial cutoff of the Fourier series expansion. As the time proceeds, the 'horns' as the remnant of Gibbs oscillations can be observed to propagate near the front of rarefaction and shock waves, and to be amplified. This artifact due to the cutoff will not be a problem when one treats the solution exactly.

Conclusion
In the large-D(dimension) limit of general relativity, the Einstein equations describing the horizon dynamics of asymptotically flat (resp. AdS) black branes are written in the form of coupled diffusion equations (2.1) and (2.2) (resp. (2.1) and (3.1)) [4,7]. While these equations are much simpler than the original Einstein equations, there is a non-linear term, making it difficult to solve exactly. Therefore, we have formulated the perturbation theory in this paper making it possible to obtain analytic results on the black-brane dynamics.
The metric functions m(t, z) and p(t, z), which represent the mass and momentum distributions in the direction of horizon z, respectively, were expanded around a uniform black-brane solution with formal small parameter ǫ as Eqs. (2.4) and (2.5). Then, the perturbative equations of motion were obtained, and solved order by order using the Laplace and Fourier transformation with respect to t and z, respectively. As the result, the general form of solutions m ℓ (t, z) and p ℓ (t, z) in the flat (resp. AdS) background were obtained (2.23) (resp. (3.7)), in which the inverse Fourier transformation of initial spectram ℓ (0, k) andp ℓ (0, k), and the spectrum of source termψ ℓ (t, k), being a polynomial of lower-order perturbations, were left to be computed.
As the example of initial conditions of perturbation, the Gaussian wave packet was considered for both the asymptotically flat and AdS black branes in Secs. 2.2 and 3.2, respectively, and the first-order solutions were written down explicitly. The resulting dynamics from this initial condition was not surprising itself. Namely, the wave pack grows and damps rapidly for the asymptotically flat and AdS black branes, respectively, which is expected from their known stability. A remarkable point revealed by this example is that the first-order solution captures the features of full-order solution rather accurately even for a finite amplitude. Thus, the convergence of expansion by ǫ (i.e., amplitude) is rather rapid and can be used for finite-amplitude perturbations for a certain class of problems.
Only for the asymptotically AdS black brane, the step-function like initial condition was considered in order to investigate the shock propagation in Sec. 3.3, and the first-order solution was explicitly written down. While the solution captures the emergence and propagation of NESS (non-equilibrium steady state) qualitatively, the first-order solution is not enough to reproduce the properties of NESS such as the values of metric functions and asymmetry of the full solution.
The discretely superposed sinusoidal waves, which can be the Fourier series expansion of an arbitrary piecewise continuous periodic function, were considered for both the asymptotically flat and AdS black branes in Secs. 2.3 and 3.4, respectively, and the first-and second-order solutions were written down explicitly. For the black brane in the flat background, the non-trivial feature of GL (Gregory-Laflamme) instability was revealed. Namely, the mode-mode coupling at the second order can make the perturbation grows even if the initial perturbation damps at the first order. For the black brane in the AdS background, the shock propagation considered in Sec. 3.3 was re-investigated. Thanks to the second-order contribution, the values and asymmetric features of the full solution were reproduced, illustrating the usefulness of the formalism in this paper.
There are many things to do by applying and generalizing the formulation presented in this paper. For instance, (i) it would be interesting to investigate further why the sign of convergence to non-uniform black string (NUBS) of asymptotically flat black string does not appear even in the second-order perturbation. (ii) One is able to investigate the various type of shock-wave propagation such as those discussed in Ref. [7] if one slightly changes the ansatz of expansion (2.4) and (2.5), which is straightforward. (iii) Including 1/D corrections and charges of background black branes [11,12,13] would increase the problems to be worked by our formalism.