Driven black holes: from Kolmogorov scaling to turbulent wakes

General relativity governs the nonlinear dynamics of spacetime, including black holes and their event horizons. We demonstrate that forced black hole horizons exhibit statistically steady turbulent spacetime dynamics consistent with Kolmogorov's theory of 1941. As a proof of principle we focus on black holes in asymptotically anti-de Sitter spacetimes in a large number of dimensions, where greater analytic control is gained. We also demonstrate that tidal deformations of the horizon induce turbulent dynamics. When set in motion relative to the horizon a deformation develops a turbulent spacetime wake, indicating that turbulent spacetime dynamics may play a role in binary mergers and other strong-field phenomena.


INTRODUCTION
The recent observations of black hole mergers [1] add to the increasing evidence that black holes exist in nature. Whilst at early times black hole mergers are described by perturbative post-Newtonian physics [2], and at late times are described by Kerr plus a handful of the longest lived quasinormal modes [3], at intermediate times there is the exciting possibility that we are witnessing the full nonlinearity of general relativity [4].
A remarkably universal consequence of non-linear dynamics is turbulence, seen across a wide variety of systems that exhibit fluid-like behavior at the largest scales. The dynamics of the chaotic cascade of vortices across scales washes out all memory of how those vortices were created in the first place, resulting in universal characteristics. Turbulence is arguably universality par excellence, the consequences of which are witnessed in all corners of nature from galaxy formation to atmospheric dynamics to a cup of tea.
In this paper we examine the possibility that the dynamics of black holes exhibit similar universal features to fluids undergoing turbulent cascades. As is well known, the physics of black hole horizons and the dynamics of fluids are closely related, depending on the context. It was appreciated early on that black-hole horizons can be thought of as fluid membranes [5][6][7]. More recently, a one-to-one map between near-equilibrium black holes solutions in asymptotically anti-de Sitter (AdS) spacetimes and the solutions to conformal hydrodynamics with particular transport coefficients has been established [8][9][10]. Recent connections between near horizon dynamics and incompressible Navier Stokes equations were studied in [11], and moreover it has recently been established that the Stokes equations govern transport properties of in-homogeneous horizons [12]. We may therefore hope that the remarkable universality observed in turbulent cascades is also seen in the dynamics of black holes. This is our present goal.
The most widely celebrated results on the universality of turbulent cascades are captured by Kolmogorov's theory of 1941 [13,14] (K41). 1 Under similarity hypotheses for homogeneous isotropic turbulence, the statistical distributions of the velocity field in the inertial range depend only on the rate of transfer of kinetic energy within the cascade, ε. Dimensional analysis then reveals that the two-point functions of the velocity field v in momentum space -here written in terms of the kinetic energy spectrum, E(k) -takes a simple scaling form, while for higher n-point functions of the velocity field in position space, arranged into longitudinal structure functions, Here the angle brackets denote statistical averages. These results apply in any number of dimensions, so long as the underlying details of the dynamics meet the nontrivial test of the similarity hypotheses. We shall demonstrate that certain driven black hole spacetimes exhibit these universal features, namely we shall numerically demonstrate that (1) holds over 1-2 decades in momentum space, and (2) holds up to n = 10 over a decade.
These properties are demonstrated for horizons with dynamics restricted to two spatial dimensions (i.e. spacetimes with nontrivial dynamics in 3+1 dimensions). In two spatial dimensions the underlying dynamics which realise the scaling hypotheses and lead to (1) take the form of an inverse cascade, with power moving from a driving scale at some high-k to low-k over time. In two dimensions there is also a direct cascade which moves from the driving-k to the UV [17], and whilst we do see this occur, its properties will not be the focus of our work.
We focus on the simplest possible setting in which we can explore such turbulent horizons, as a proof of principle of the above properties. To this end we restrict our attention to horizons with planar (or toroidal) topology, rather than spherical. In asymptotically flat spacetimes planar horizons present an unstable starting point because of the Gregory-Laflamme instability [18], and are therefore unnatural objects to consider in our present goal of studying the dynamical response to a forcing term. We focus instead on asymptotically AdS spacetimes where planar horizons are intrinsically stable. Such black holes are also of interest as models of strongly interacting many body systems through the AdS/CFT correspondence [19]. Due to the universality of the mechanism of turbulence we may hope that our proof of principle examples shed light onto the universal dynamics of black holes in general, including those in asymptotically flat spacetimes.
In obtaining our results we do not utilise a hydrodynamic expansion, i.e. we do not carry out perturbation theory in gradients. Without this technical crutch, the connection between black holes and fluid-like dynamics is not readily apparent from the Einstein equations. However, the connection is once more seen if one treats the inverse spacetime dimension 1/D as a perturbative parameter according to the seminal constructions of [20][21][22][23].
At large D there is a separation of scales between the black hole size, r 0 and the region occupied by a nontrivial gravitational potential, r 0 /D. A separation of scales signals an effective theory, which can be constructed by analytically solving the integrals for radial evolution. What remains is a set of constraint equations in D − 1 dimensions. It is these constraints that resemble fluid-dynamics equations. Crucially, however, even though they are perturbative in 1/D, they are exact in gradients. Therefore, far from a mere technical simplification, the 1/D expansion allows us to directly connect black holes with the turbulent behaviour of a class of fluid-like equations that are exact in gradients. We thus note that these equations may be valuable in their own right as a natural candidate for the study of turbulent behaviour, as compared to those obtained in an arbitrarily truncated gradient expansion (such as the Navier-Stokes equations). We return to this point in the discussion.
The driving we consider is obtained by adding forcing terms directly to the black hole equations, F i (t, x), in a way that injects vorticity consistent with the requirements of statistical homogeneity and isotropy. Due to the nature of our particular setup, we are also afforded the opportunity to drive the horizon fluid by turning on a deformation to the gravitational potential at the boundary of AdS, γ tt (t, x), using a generalisation of the large-D equations derived in [24]. Previous work has explored decaying turbulent dynamics of black holes that result after starting from unstable initial conditions [25][26][27][28], where [28] also utilised the large D expansion as we do here. To distinguish the work we present here, we do not require starting with unstable initial data, and the driving allows us to achieve a quasistationary 2 turbulent regime that can be compared with the predictions of K41. Motivated by the holographic connection there has also been a focus on turbulence in 2 + 1 conformal hydrodynamics [29][30][31][32] where analyses of the energy spectrum and comparisons to (1) are made. Quasi-normal mode resonances of a rapidly spinning Kerr black hole have been argued to result in a phenomenon resembling an inverse turbulent cascade in 2+1 dimensions [33].

LARGE D BLACK HOLE DYNAMICS IN ADS
We shall now record the most salient points regarding the effective theory that describes the dynamics of black branes in asymptotically AdS D spacetimes at large D, derived in [20][21][22][23] (see also [34]). In taking the large D limit we focus on the near horizon region of the black brane, so that the resulting theory effectively describes how the near horizon deformations evolve in time. Furthermore, a mathematical simplification arises that allows us to solve the constraints in the radial direction, so that the basic variables can be readily related to the energy and momentum density of a fluid. This effective theory was later extended in [24] by introducing a general class of boundary conditions which induce changes in the gravitational potential on the boundary of the nearhorizon region. These map to sources for the dual stressenergy tensor in the dual picture. Among the class of deformations derived in [24], here we will only consider the one corresponding to adding a source for the energy density.
Our action is simply Einstein-Hilbert with negative cosmological constant in D dimensions where Λ = −(D − 1)(D − 2)/2. We choose a coordinate system adapted to the black brane; we split the spacetime coordinate X α into X α = {t, r, x i , y a }, where r is the coordinate transverse to the brane, which also plays the role of the holographic coordinate, t, x i and y a are the coordinates along the black brane (i = 1..p, a = 1..ñ), so that D =ñ + p + 2. The reason for the split between x i and y a is that we restrict only to dynamics in a subset p of the boundary spatial directions, in other words, we dimensionally reduce on anñ-torus and keep only the zero modes. As in [20], we take the large D limit by takingñ → ∞ keeping p fixed. This is facilitated by choosing a metric ansatz of the form As shown in [22][23][24], it is consistent to solve the Einstein equations with the following ansatz as a perturbative expansion in where R = rñ. Here, a, p i are the energy and momentum density of the black brane, while γ tt is a deformation of the AdS-boundary metric. As such, we can think of it as providing an external gravitational potential, or providing an adjustment to the gravitational environment in which the black hole lives. 3 The ansatz (4) can be thought of as a radial foliation of spacetime. As a consequence, the Einstein equations split into evolution equations in the radial direction and constraints on r = const surfaces. The radial equations can be solved order by order in 1/ñ, so we are left with the constraints. At leading order in 1/ñ, these are where in the presence of the boundary metric deformation we have F i = a 2 ∂ i γ tt and where indices are raised and lowered with the flat metric δ ij . Equations (9), (10) correspond to conservation equations associated to time and spatial translation invariance, and behave in accordance with expectations for a viscous fluid, complete with sound and shear-diffusion modes as detailed in [23]. These equations will serve as the basis of our analysis. All that remains is to provide the relations between these gravitational variables (a, p i ), in which we carry out numerical evolution, and a set of fluid variables: the fluid energy density is a, whilst the velocity is v i = (p i −∂ i a)/a, see [24] for details. Note that this large D limit results in a set of non-relativistic equations.

KOLMOGOROV SCALING
In this section we present the key results of our work, that when appropriately driven by a homogeneous and isotropic forcing function, the large-D equations of general relativity in AdS (9)-(10) exhibit turbulent behaviour consistent with the predictions of K41 theory. We verify the agreement with K41 by carrying out 256 independent random realisations and using them to compute statistical properties of the velocity field, v i .
To achieve the conditions required for K41, namely homogeneity, isotropy and driving, we replace F i on the right hand side of (10) by an explicit forcing function, instead of utilising γ tt , which we set to zero. While γ tt is such a forcing function, it appears as the derivative of a scalar and cannot be directly used to supply a source of vorticity. 4 We shall consider flows in the presence of nontrivial γ tt later. We adopt periodic boundary conditions for a torus of size L × L. Details of the choice of F i , the numerical methods used and their implementation are given in the appendix.
To acquaint the reader, we first illustrate the pattern of vorticity, ω = ij ∂ i v j , obtained during a single realisation at late times in figure 1 (left). The time-dependence of this picture resembles that of a liquid of vortices moving in larger coherent structures. Larger structures can be seen in this picture in accordance with power beginning to accumulate on the largest scales available, i.e. the torus size. This part of the spectrum, subject to finitesize effects, is expected to lie outside the inertial range.  The inertial range is seen to grow over time, from the driving scale downwards, i.e. an inverse cascade. The inertial range approaches two decades-worth of K41 scaling before meeting the finite size of the box. At this time power begins to accumulate at low k and the scaling is destroyed there. Our simulations thus demonstrate quasistationary turbulence. Interestingly, since we are working on a torus, the structure begins to resemble that of a square vortex-antivortex lattice at late times, but continually fed from high k. This phenomenon also goes by the term 'energy condensation' [35].
Next we demonstrate that the longitudinal structure functions (2) are also seen in our simulations, from n = 2 up to n = 10 in figure 2 (right). For this position-space calculation we pick a single direction in which we compute velocity differences,x, average over each row labelled by y, and then over 256 realisations. The consistency with K41 is indicative of the absence of intermittency corrections; we are driving the system homogeneously and not providing any window of opportunity for intermittent laminar flow to develop.
So far we have provided evidence for the presence of K41 scaling in the form of power spectra and structure functions. Finally, we provide further support for why K41 predictions are seen here. A crucial part of the K41 analysis is that there is a dominant scale in the problem, ε, the rate of energy transfer, which is taken to be a constant. This governs how quickly kinetic energy is passed between vortices of different scales. Given that our inertial range is growing over time, exhibiting an indirect cascade from the injected k-scale downwards, we must be providing a constant injection of kinetic energy via the F i forcing. This is indeed the case, as is clearly shown in figure 3. We see a constant growth of v 2 with time (integrated over the torus), whilst the total enstrophy (Ω = d 2 xω 2 ) remains approximately constant. 5

TURBULENT WAKES
So far in this work we have examined universal turbulent dynamics of black hole horizons by explicitly adding a set of custom forcing terms F i to the Einstein equations which directly source vorticity. Whilst such terms are the most convenient for achieving homogeneous isotropic turbulence and comparing to K41, it is desirable to have a fully gravitational realisation of turbulent dynamics. In fact our setup is capable to also address this question directly. To this end, in this section we consider the case where the forcing is instead provided by a gravitational deformation of the AdS conformal boundary metric, γ tt (t, x). The associated equations for such deformations simply result in a forcing term F i = a∇ i γ tt /2 as in (10). We emphasize, however, that this means that 5 For a discussion of the approximate constancy of Ω for unforced equations see [28]. Note that here we are forcing the equations and the approximate constancy of Ω emerges from the detailed dynamics. As further support we also show the approximate constancy of the total enstrophy of the system in red.
the forcing term is entirely physical, and results from an explicit gravitational source. The fact that we can implement such a scheme is another illustration of the power of our setup for the study of gravitational turbulence.
For this simulation we set initial conditions where the fluid is in uniform motion, and then quench γ tt from zero to a symmetric Gaussian profile at a fixed location x on the torus. The resultant dynamics, shown in figure 4, clearly demonstrates that this deformation develops a turbulent gravitational wake on the horizon. Note that this example is not statistically homogeneous nor isotropic, and moreover it is decaying over time, and so we do not expect it to fall into the class of turbulent flows meeting the basic requirements to be described by K41. 6 We note that since drag force varies qualitatively between turbulent and laminar wakes we expect that this tidally-induced turbulence can impact the dynamics of strong-field gravitational processes. For example, the mechanism we have proposed may be of relevance to the dynamics of quark-gluon plasmas via AdS/CFT. To add further relevance to this scenario, the γ tt calculation we have performed may be taken in the same spirit as studies of near-extremal Kerr black holes and their near-horizon regions using CFT [36]. There, when a massive body falls into the near horizon region it appears as a source term in the CFT, as pointed out in [37]. In our case, the γ tt may be viewed as the gravitational deformation due to such a massive body, and the appearance of a turbulent wake in this context may affect the plunge dynamics and associated waveforms for BH-BH or BH-NS mergers. Of course the extrapolation of our results to asymptotically flat spacetimes is not straightforward, particularly with regards to identifying an appropriate hierarchy of scales, and a direct computation would be required to confirm its astrophysical relevance.

DISCUSSION
In this paper we have worked at strictly infinite spacetime dimension, D, though we considered dynamics constrained to 2 + 1 of them. There are however reasons to expect that results we obtained continue to hold in lower dimensions. First, the dimensional analysis behind K41 scaling is dimension independent; indeed, k −5/3 is predicted and seen in a range of 2 + 1 and 3 + 1 dimensional scenarios, as well as 2 + 1 dynamics of an infinitedimensional system as studied here. This universality occurs despite the clear differences in the dynamics that underpin the cascades; in 3 + 1 the cascade is direct, whilst in 2 + 1 there is also an inverse cascade. Thus whilst the detailed dynamics may change substantially as we lower the number of dimensions, we anticipate that K41 remains robust. Second, the map between gravity and hydrodynamics holds in all spacetime dimensions D.
We also highlighted the large-D equations as a potentially useful model in the study of turbulence in general. These equations are simultaneously dissipative and exact in gradients, by virtue of the parametric control afforded by the 1/D expansion. This should be contrasted with the usual treatment of hydrodynamics in a dynamical setting, where one typically truncates at a finite order and treats the resulting system of equations as exact -a procedure which fundamentally changes the theory. As a consequence of this change one discovers physically undesirable qualities such as instabilities and acausal behaviour both for relativistic theories [38], and non-relativistic theories [39]. Furthermore, one must of course also verify post-hoc that the solution remained a good approximation within the framework of a perturbative gradient expansion. None of these concerns apply to our setup.
Finally we emphasise that whilst the large-D equations of motion appear fluid-like, these are the radial constraint pieces of a full solution to the Einstein equations. The radial evolution equations were solved analytically as part of the construction, and so our solutions correspond to full black hole solutions to the Einstein equations. Thus the K41 behaviour we illustrate here from the fluid perspective is naturally encoded -through the a, p i variables -as geometric data in these spacetimes. In a more general setting the appropriate geometric data corresponding to the fluid observables may be difficult to identify, in which case it may be helpful to consider approaches to visualising horizon vorticity, for example [40].

Numerical Method
The equations we are solving in this case are given by (9) and (10) with the right hand side of (10) supplemented by the additional term Notice that this enters the equations in such a way to produce vorticity without directly sourcing a. We consider the system on a torus with The force F is constructed in such a way as to drive the system isotropically at a fixed energy scale, m, and consists of random combinations of all the Fourier amplitude coefficients whose corresponding wave vectors lie close to a circle of fixed radius in wave vector space. Specifically, given a discretisation scheme for the spatial directions, F is given by where k i is a set of M vectors in Fourier space that was sampled over an annular region | k (i) | = m ± δm in an isotropic way. The coefficient A is a fixed overall amplitude controlling the strength of the forcing. At times that are integer multiples of ∆t the mode coefficients c (i) are drawn from a normal distribution with zero mean and variance ∆t, normalised such that while the angles φ (i) are drawn from a uniform distribution. These random variables are assigned at times that are integer multiples of an interval ∆t. For the times in between, i.e. from t 1 = n∆t for some n ∈ Z + and t 2 = t 1 + ∆t the forcing function is interpolated, for t 1 ≤ t < t 2 . This choice ensures the squarenormalisation is maintained during the interpolation, with uncorrelated cross-terms cancelling after averaging. We typically used ∆t = 10δt, where δt is the time step used for the numerical evolution. A couple of comments are now in order. First of all, the benefit of forcing the system at a particular energy scale, m, is that there will be a sharp peak on the power spectrum at that scale, while the remaining of the spectrum will not be contaminated, allowing easier identification of scaling behaviours at larger and smaller energy scales. Furthermore, interpolating over the random amplitudes in this way allows us to use deterministic time-stepping techniques instead of stochastic ones.
Let us now move on to discuss the details of the numerical method used for solving these equations. We utilise a uniformly spaced discretisation of x, y with N x , N y grid points respectively, taking N x = N y = N . We adopt fourth-order finite difference approximations of the derivative operators. The variables a, p x and p y are then evolved forward in time using fourth-order Runge-Kutta time-stepping (RK4), subject to periodic boundary conditions. In addition we add a sixth order Kreiss-Oliger dissipation term [41]. This is done by replacing the time derivative operator with where h = L/N is the grid spacing. Note that this term approaches to zero faster than the error in the spatial finite differences as h → 0. At the highest resolutions considered, N = 1024, we have also performed numerical evolutions with η = 0. The difference between η = 0 and η = 0 is only visible near the UV in the power spectrum, as expected. The properties of the solutions away from the UV are not affected by η.
The last piece of information needed for the time evolution is to specify the initial conditions. For simplicity, we consider homogeneous initial conditions for the evolution functions, namely a = 2 , p i = 0.
In the particular evolutions discussed we considered N = 1024, M = 100, η = 0.4, δt = 6.4 × 10 −4 L, A = 0.005, and m = 100. In order to observe turbulence, the associated Reynolds numbers, Re ∼ L, should be large, and in our simulations we use L = 2×10 5 . Random numbers are generated using the Mersenne twister algorithm [42], and normally distributed variables are obtained using a Box-Muller transform [43]. We work with double floating-point precision.

Overview of GPU implementation
Given the relative simplicity of the explicit update steps we have found the use of GPUs to be of significant utility. The methods used are standard and do not warrant a detailed exposition, however a high-level overview of the step-by-step procedure used may be of value. This is given below.
Seed the pseudo-random number generator uniquely for each run. Allocate memory on the device to store, for all variables a, p x , p y : their spatial derivatives, first time derivatives, and their values at four RK4 intermediate steps. Allocate memory to store 100 pairs of k x , k y values specifying the location of driving points in momentum space, as well as two buffers for the associated amplitudes and phases (two buffers are required because we interpolate between two sets of random variables over time). Initialise the allocated memory with initial data for the run, the driving point locations, initial random amplitudes and phases (all generated on the CPU). Then, for each complete time step: 1. Check if the random amplitudes and phases need updating at this time step. If so, cycle between the two buffers, compute new random values on the CPU and copy them to the device.
2. For each RK4 intermediate step: • Compute spatial derivatives in the x-direction (of the intermediate RK4 variable appropriate for this step). This is performed by splitting the N × N grid into N thread blocks, one for each y-value, and using N threads in each thread block to perform the computations at each x. In this way we can set the shared memory for each thread block to be the entire y-column, extended to include an extra 8 ghost points for periodicity. See, for example, [44].
• Transpose and repeat for all y-derivatives (mixed second derivatives do not appear in the equations).
• Compute time derivatives using the previously computed values through the equations of motion.
• Add F i terms to the time derivatives as required. These are constructed from interpolating between F i values evaluated using initial and final amplitude/phase buffers.
• Fill the next RK4 intermediate buffer using the time derivatives. For the final RK4 intermediate step, the next RK4 buffer written is the first RK4 buffer.
3. Rarely, copy the field values (i.e. the first RK4 buffer) from the device and write to disk. In practise we found it invaluable to use the visualisation software VisIt [45] and so we also write an appropriate 'brick of values' descriptor file to accompany each binary data file.
In this way the bulk of the evolution is performed on the device itself. The exception to this is also the slowest part of this procedure, namely the device synchronisation bottlenecks for step 1.