Non-linear black hole dynamics and Carrollian fluids

The dynamics of black hole horizons has recently been linked to that of Carrollian fluids. This results in a dictionary between geometrical quantities and those of a fluid with unusual properties due its underlying Carrollian symmetries. In this work we explore this relation in dynamical settings with the interest of shedding light on either side by relevant observations. In particular: we discuss how the null surface where the Carrollian fluid evolves is affected by its behavior; that the fluid's equilibration properties are tied to teleological considerations; the connection of higher derivative contributions as both source of energy and dissipation for the fluid and the non-linear behavior of black holes. This latter point, connects with discussions of non-linear modes in the relaxation to equilibrium of perturbed black holes.


Introduction
The understanding of gravity around black holes has long been the subject of scrutiny.These objects have been a tremendous source of inspiration and central to questions ranging from astrophysics to field theories.At the same time, they have long presented puzzles that have often benefited from insights from other branches of physics.In particular, the membrane paradigm [1][2][3] connecting the physics of black holes with that of relativistic fluids at a timelike surface in the vicinity of the event horizon shed light on particular aspects of the region's behavior.However, freedom in the choice of such surface could obscure somewhat particular connections.In the context of black holes in asymptotically AdS spacetimes, this issue was largely addressed by choosing the AdS boundary -a unique surface in the spacetime-and a clean fluid/gravity correspondence was established [4][5][6][7].This enabled, in particular, motivating the possibility of turbulent behavior in General Relativity with a negative cosmological constant [7][8][9].In the opposite direction, insights on possible completion of relativistic hydrodynamics and the analysis of transport coefficients, have been discussed from the gravitational side (e.g.[10][11][12]) In the context of asymptotically flat spacetimes special surfaces to attempt a similar connection are future null infinity I + or the black hole event horizon.Since both are null surfaces, special care must be taken to deal with their singular nature.In the case of the event horizon, such a program has recently been carried out, leading to the discovery that gravity around the black hole is dual to Carrollian fluids [13].These fluids appear as the limit of relativistic hydrodynamics as the speed of light goes to zero.Even though in that limit a single particle is not allowed to have any dynamics [14], interacting particles retain some dynamics in the limit [15,16].The physical nature of this Carrollian hydrodynamics remains puzzling; in particular the notion of entropy and a microscopic description of such fluids are still work in progress (e.g.[17,18]).In the current work, we further explore this connection by focusing on the dynamical aspects.In particular, we study the black hole region's behavior both geometrically and in terms of a (Carrollian) fluid description.To this end, we focus on the relaxation to equilibrium of a black hole that is perturbed by a pulse of gravitational waves -physically motivated by scenarios describing the after-merger stage of binary black hole collisions.We find that the non-linear dynamics results in an excitation of shorter wavelength modes which occurs prior to the perturbation -related to the teleological character of the horizon.Depending on the strength of the perturbation, appreciable non-linear energy exchange is induced and short wavelength modes can become even more relevant than those modes excited by the perturbation.To understand this behavior, we analyze the structure of the equations and note the Carrollian fluid dynamics has a non-trivial/non-linear relation with the "teleological" volume (the horizon's surface) it occupies.To unearth this effect, we connect to the bulk spacetime thus revealing this intricate connection.In particular, this allows for motivating a definition of gravitational Reynolds number to help distinguish between laminar and turbulent regimes.Further, since our studies on this topic are intrinsically linked with the behavior of horizons, and in particular their non-linear behavior, we draw from also a complementary way to explore dynamical horizons.Particularly, the transition from non-linear to linear regimes.This work is organized as follows: section 2 revisits the geometry of null surfaces and how can they naturally be endowed with a Carrollian structure.Section 3 establishes the gravitational dynamics of this geometry, and its mapping to Carrollian fluids by properly identifying geometric and fluid degrees of freedom.In section 4 we describe our nonlinear perturbative set-up and the equations that will determine the horizon's dynamics.We introduce the quantities that we will monitor during the evolution in section 4.2 and discuss our numerical implementation in section 5. Finally we analyze the results of our simulations in section 6.We end with final observations in section 7. We also include in the appendices more details on Carrollian geometry as well as further potential connections with hydrodynamical behavior that would be explored elsewhere.
In what follows, we restrict to 3+1 spacetime dimensions, indices with upper case letters are intrinsic, 2D, indices A, B = 2, 3 and lower case ones denote spacetime indices a, b = 0 . . .3.

Geometry of null surfaces
Let H be a null surface, and assume it is located at some adapted coordinate ρ = 0.It is always possible to find a chart containing the surface such that the metric is written in the form [19,20]: where v is the advanced null time, ρ parameterizes the "distance" to the surface, and x A are the angular coordinates on the sphere, A, B = 2, 3. We expand the metric in powers of ρ, where κ, U A , Ω AB and λ AB are functions of (v, x A ).Note this gauge is a restricted version of the one chosen in [21]; their geometrical characterization also applies for our choice, upon making the necessary identifications.We next define the geometric variables that will be used to characterize the evolution of the null surface.
The null surface has a null normal vector l a that is also tangent to it.We introduce another null normal vector to the surface, n a , normalized such that n a l a = −1.These two vectors are given by 3) It is useful to characterize the geometry of the null surface in terms of quantities whose definition do not depend on the coordinates, but that acquire a simple expression in terms of our coordinate choice.Notice that the induced metric on the surface H is given by h ab = g ab + 2l (a n b) , which yields Therefore, for each null normal m = {l a , n a }, which we debote with a upperscript in between brackets, we can construct its deviation tensor.This is defined as the projection onto the null surface of its covariant derivative that we decompose into its trace (the expansion) and its traceless symmetric part (the shear). 1 These are given now by where the dot (˙) denotes the coordinate derivative with respect to v, and Ω = det Ω AB .Finally, the last variable that we will use is the Hajiceck one-form, which is given by We can provide a geometric meaning to all the quantities involved in the metric expansion (2.2).The parameter κ is the non-affinity of our parameterization, which coincides with the surface gravity for stationary black holes since l a ∇ a l b = κl b .The shift term U A is related to the Hajiceck one-form.The angular metric Ω AB describes the intrinsic metric of the null surface, whereas the term λ AB describes the angular metric right outside of the surface itself, and therefore takes into account ambient effects.The geometry of a null surface can be endowed with a Carrollian structure.This observation comes from the realization that the diffeomorphisms ξ that preserve the gauge (2.1), once restricted to the null surface, are just which are the infinitesimal version of the Carrollian diffeomorphisms [23].
A Carroll structure is given by a surface H, a degenerate metric on that surface h AB and a preferred generator, A = h A b l b such that l A h AB = 0.Although this structure can be defined intrinsically for any given null surface (e.g.[24][25][26], it is perhaps easier to make the identification by comparing the near-horizon geometry (2.1) with the Randers-Papapetrou parameterization [27], given by (2.9) Then, the null surface of interest H is endowed with a Carroll structure with degenerate geometry Ω AB and the following identifications: with α the lapse function and b A the shift or temporal connection.More importantly, the distance to the surface ρ plays the role of the speed of light c 2 .This motivates the following intuition: A timelike surface located at some distance ρ to the horizon has a relativistic fluid dual as per the usual membrane paradigm [1][2][3].Therefore, as one takes the limit in which the surface becomes null ρ → 0, then one must take the ultrarrelativistic or Carrollian limit c 2 → 0 on the dual fluid.

Dynamics of null surfaces
The membrane paradigm [1] identifies the gravitational dynamics on the stretched horizon to the evolution of a relativistic fluid in 2 + 1 dimensions.The stretched horizon is defined to be a timelike surface, sitting at some distance outside of the event horizon of a given black hole spacetime.If one insists on finding a duality with the dynamics of a relativistic fluid at the level of the event horizon, the quantities involved in this process diverge.These can be regularized by introducing a family of fiducial inertial observers [2].However, this procedure depends explicitly on the foliation of stretched horizons chosen, and therefore it is not universal.One can avoid the divergences by identifying the degrees of freedom at the level of the event horizon with those of a Carrollian fluid, as it was shown in [13].In this section we frame this construction in a different way, not requiring the use of a foliation of stretched horizons.The well-defined gravitational action of a spacetime M with a null boundary H is given by [28]: where we have fixed units such that 8πG = 1.The variations with respect to the boundary degrees of freedom were shown in [29] to be equivalent to the conservation of the null stress-energy tensor given by where , where l a is the null generator of the surface, and W = W a a is its trace.Notice how this stress energy tensor is similar in structure and in construction to the Brown-York stress energy tensor [30,31], although that is constructed for a timelike surface.A discussion about the relation between this stress energy tensor and the null limit of a timelike surface is given in [29].Also, from the full spacetime point of view, this stress energy tensor is obtained through a suitable projection with h ab .This projector is not unique but is invariant with respect of scalings, see [32].Regardless, once a pair of null vectors {l, n}is defined a consistent decomposition is made, and employed to study a black hole's behavior and explore its dynamics from several angles as we do in this work.
Gravitational dynamics on the null boundary are given by the conservation of this null stress energy tensor: where R ab is the Ricci tensor and we are assuming a vacuum spacetime.These equations give directly the Raychaudhuri equation (the b = 0 component) and the Damour-Navier-Stokes equation (the b = 2, 3 components), also referred to sometimes as the Hajiceck equation [3]. 2 These equations are given, in terms of the geometric variables, by where ∇ A here is the covariant derivative associated to the metric Ω AB .Notice that these equations are just a restricted version of the constraint equations derived in [21], which also allow for a geometric description [33].By defining everything from the variational principle adapted for a null surface we have avoided introducing a foliation of stretched horizons.

Mapping to the Carrollian fluid
A Carrollian fluid can be described as the c → 0 limit of a relativistic fluid [26,34].Upon taking the limit, the fluid is described by a similar set of hydrodynamical variables: the energy density e, the pressure p, a pair of heat currents J A and π A , and a pair of dissipation tensors Σ AB and Π AB .It could seem that by introducing another pair of heat current and dissipation tensors we are increasing the number of degrees of freedom of the system.This is not the case, since these two additional terms enter in a different order in the expansion in powers of c of the fluid, and must satisfy two constraints: Here θ (l) and N (l) AB are the expansion and shear of the Carroll structure where the fluid is defined. 3The dot denotes the time derivative ẏ = ∂ v y.We have also introduced the Carrollian acceleration [34] and the Carrollian derivative ∇A (see appendix A).Then, the Carrollian fluid is described by the following two evolution equations In the above AB is the Carrollian vorticity (see next section and appendix A).Now one realizes that these equations are just the gravitational equations (3.4) in disguise.
Identifying terms is straightforward, and yields the following fluid/gravity dictionary: e = θ (l) , p = −κ, In obtaining this dictionary, notice that the second pair of heat flux and dissipation tensor Σ AB = J A = 0 vanish.The vanishing of J A is due to the identification with the Raychaudhuri equation, and then the constraints (3.5) necessarily imply the vanishing of Σ AB .A very similar result is obtained by mapping the gravitational stress energy tensor (3.2) to the stress energy tensor of a Carrollian fluid [26,32].The main difference lies on the ambiguity in assigning contributions to pressure and the bulk viscosity ζ term of the dissipation.Their dissipation tensor is required to be traceless, which implicitly assumes vanishing bulk viscosity.This is the case, e.g. in the fluid/gravity correspondence for asymptotically AdS spacetimes, where the fluid retains a conformal symmetry.However in the case of the horizon there are no symmetry arguments which would result in a vanishing bulk viscosity.Therefore, one would expect, in a similar fashion as in the identifications done in the membrane paradigm works [2], a non-vanishing bulk viscosity, thus we adopt this point of view.
Our result also yields a very similar dictionary to the one first proposed in [13], but with some terms that depend on the gradients of the surface gravity absent.The difference was discussed in [29], where it was explained how those terms appear as additional terms due to the regularization needed to take the null limit, when starting from a timelike surface.Since we do not consider a timelike stretched horizon at any time, these terms do not appear.
The resulting correspondence maps the dynamics on a black hole horizon to those of a Carrollian fluid with very particular properties.Importantly, we highlight the pressure is negative, and the energy density vanishes in equilibrium.At this point we find it relevant to note that there is an ambiguity in the identification: we could have chosen to re-define the dictionary (3.8) with a negative energy density, positive pressure and positive bulk viscosity.Notice that the sign of the shear viscosity η would remain unchanged.However, even though a fluid with negative pressure poses a challenging question about its interpretation it has a precedent in the membrane paradigm [1].On the other hand a negative energy density seems rather unphysical.A thermodynamical formalism for Carrollian fluids could help clarify this, as well as the question on the entropy current for Carrollian fluids.
Before ending this section, we highlight an important subtlety.While it is clear that a horizon can be mapped to a Carrollian fluid, in the dynamical regime, there are ingredients of such mapping that are yet to be clarified and which are inescapably linked to the very volume where such a fluid "lives".That is, its backreaction on it, as well as the influence of the embedding spacetime are important to close the system and studying its behavior.Furthermore, for the specific case of a horizon, the Carrollian fluid's dynamics is also strongly determined by the teleological character of this surface.In particular, the requirement that in the asymptotic future θ (l) → 0 which, in turn, implies θ (l) > 0 at early times in the dynamical case.In what follows, we perform a perturbative analysis of the spacetime in the neighborhood of the horizon to illustrate this relation and motivate further connections.

Perturbative set-up
So far, the previous analysis has been purely kinematical and the established dictionary, while interesting in its own right can potentially be misinterpreted if translating from fluid considerations lightly.Namely, in the case of fluids dissipative effects are encoded by the shear and expansion appearing in the tensor Π AB ; however, for Carrollian fluids dual to black hole horizon these need not be the case as the horizon is affected by energy accretion.Such process certainly affects the resulting fluid dynamics in a non-linear fashion, which is not straightforwardly appreciated in the equations (3.7).To elucidate this behavior, we here consider a horizon perturbed by incoming gravitational waves, and treat the system perturbatively to sufficiently high order to keep track of its behavior.We note that in full non-linearity one should study this problem by tackling the full set of equations in a framework making contact with the horizon's behavior.A transparent way to this end was presented in [19,35].Alternatively, insights can be gained through the dynamical horizon approach [36,37] or via a suitable second order perturbative approach [38,39].However, we here want to keep a close connection with the Carrollian fluid structure and find it convenient to proceed as follows.
We will look at a simplified situation: a spherically symmetric black hole perturbed by a gravitational wave. 4The geometry will be described by ab , where 1 is some perturbative parameter and g Sch ab the Schwarzschild metric, given by where q AB is the usual unit sphere metric.For consistency, the perturbative expansion is carried to O(2) as the 'energy' e ≡ θ (l) is a quantity of that order outside equilibrium.The perturbation is also expanded in the gauge (2.1), so we define where all terms (k, h A , c AB , s AB ; K, H A , C AB , S AB ) are functions of (v, x A ) and we assume that the angular part is traceless to each order, q AB c AB = q AB C AB = q AB s AB = q AB S AB = 0.The source of perturbations is encoded by the trace-free (first and second order in tensors {s AB , S AB }.They encode perturbations incoming from the bulk.We will define these perturbations with a suitable, physically motivated, profile.The rest of the perturbation degrees of freedom are fixed by solving the Einstein equations in vacuum R ab = 0. We observe that the R vv and R vA components give just the Raychaudhuri and Damour equations (3.4).Moreover, as it was observed in [19,21], the R ρb components constrain the higher order terms in the ρ expansion of the metric.Finally, we split the R AB equations into its trace and trace-free symmetric parts: where ∇ A is the covariant derivative associated to the metric Ω AB , and R AB and R are the Ricci and scalar curvatures associated to that connection.The above geometric variables are the ones defined in section 2, associated to the choice of null basis (2. 3) The T-F superscript denotes the trace-free symmetric part.We now expand all relevant terms up to second order in .First, for convenience, we write the Ricci scalar as The expansions and shears can then be expressed as, with X = c AB (c AB + s AB ).We can then re-express Einstein equations and solve for each relevant order, obtaining: From the trace of R AB : (4.5) and from the trace-free part of R AB : ) These equations make the connection of quantities intrinsic to the horizon and the behavior of the embedding space evident.In particular, this implies the Carrollian fluid dynamics can be conveniently interpreted from this point of view.For instance, to leading order, ), the shear associated to the Hajicek field and in the second line we replaced to leading order with the shear associated to n a .Clearly, Π AB contains contributions of dissipation as well as source of energy.As well, the equation of state to leading order results: where D A is the covariant derivative associated to the metric on the sphere q AB .Notice that for any tensorial quantity T we can write ∇ A T = D A T + Γ T + . . .for some terms Γ depending on derivatives of c AB .The above equation indicates that the equation of state depends on both the metric of spheres at v = const and divergence of the Hajicek field.
The evolution equations, to order two in , is given by Eq. (4.6) and the following: where we have introduced the variable ϑ = θ (l) , since it is the only expansion that appears explicitly in our system as a dynamical variable.We choose to express them in terms of spinweighted quantities for convenience when implementing them numerically.In order to do this, we introduce a complex dyad q A on the sphere such that q AB = q (A qB) , and normalized with q A qA = 2, q A q A = qA qA = 0.Then, the Hajiceck vector field is described by the complex function H = q A H A , and symmetric rank-2 tensors T AB , {c AB , C AB , s AB , S AB }, by the complex function 2T = q A q B T AB .The covariant derivatives on the sphere can also be projected into this dyad to construct spin raising and lowering operators, dubbed ð and ð, as discussed in [41,42].
To sum up, we started with a set of first order {k, h A , c AB } and second order {K, H A , C AB } dynamical degrees of freedom, and some source terms that describe perturbations incoming from the bulk {s AB , S AB }.We have used a subset of Einstein equations (R AB = 0) to eliminate algebraically the corrections to the surface gravity k and K in terms of other dynamical variables as well as derive evolution equations for {c AB , C AB }.To leading order in and in spin-weighted form, the resulting equations are: (4.10) The second order terms satisfy the evolution equations: (4.11)For convenience, we set 2m = 1; further, as we will perturb the system in the longwavelength regime, we drop terms containing third derivatives assuming a quasispherical regime of sorts [43].
In order to complete the description of the evolution of the system, we need to provide boundary conditions and a prescription for the incoming gravitational wave, which is encoded in the variable s.For the latter, we give a compact support pulse of gravitational waves, with frequency and decay time associated to quasi-normal modes in Schwarzschild [44].For the former, we note that unless asymptotic boundary conditions in the future are given such that ϑ → 0 the null surface will either expand to null infinity or converge towards the singularity.Therefore, we demand that ϑ, c, and C are zero at asymptotic time v → ∞.We complement this with choosing a vanishing Hajiceck field at v → ∞, since we assume that asymptotically the metric is given by the Schwarzschild black hole in "non-rotating" coordinates.This amounts to fixing h = H = 0 in the boundary.

Non-linear considerations
Before moving on to the numerical implementation of the relevant equations, we note that resumming the equations for { ḣa , ḢA } one can see the non-linear structure of the Hajicek equation is, schematically, of the form: This equation reveals features hidden at first sight in equations (3.7),(4.9).Specifically, the non-linear structure of the underlying equation for Hajicek field (or the momentum flux) and a relevant consequence.Considering derivatives scale as ∇ ∼ 1/L, where L is the relevant scale, as It is tempting to propose a Reynolds number defined as which bears a close resemblance to the standard expression in fluids as well as the one suggested in [45].A new ingredient here is the presence of c AB which accounts for its role in defining the volume on which the dynamics develops.

Monitoring quantities
We will focus on a few quantities to explore the horizon's (=fluid's) behavior.First, the black hole's area.This can be obtained from the total expansion by direct integration (recall that from the Carrollian fluid dual perspective the expansion is the fluid's energy density).We also monitor the expansion and its angular structure.Out of the evolution variables {θ, h A , H A , c AB , C AB } we will monitor the re-summed combinations H A = h A + H A and C AB = c AB + C AB .Notice that while in the unperturbed regime, H A is a gauge choice, in the non-trivial regime its dynamics can help discern nonlinear behavior, as argued in the previous subsection.In particular, a memory imprint will result.In the case of relativistic or Navier-Stokes hydrodynamics, the vorticity associated to the fluid's velocity is quite informative.In a Carrollian fluid however, the dual to such velocity is the null vector l a which, by construction does not have vorticity and by virtue of it being null, does not carry a similar physical content.However the Carrollian acceleration π A or the Hajicek field H A provide a sense of energy flux.Thus, and recalling in the Eckart frame such flux governs, in particular, entropy growth; we also monitor the behavior of the vorticity associated to H A , defined as This vorticity coincides with the Carrollian vorticity defined in [46].Moreover it also coincides with the Carrollian vorticity defined for a Carrollian fluid in [34] up to a term that is sub-leading in perturbative order.Further, equation (4.7) highlights the role of the Hajicek field's shear as a viscous contribution.Geometrically, the vorticity is given by the imaginary part of the Newman-Penrose Ψ 2 .The duality motivates an interesting physical interpretation to a purely geometrical quantity and has also been highlighted in the context of asymptotically AdS spacetimes and the fluid-gravity correspondence on the spacetime timelike boundary [47,48].

Numerical Implementation
We evolve equations (4.10), (4.11) using the method of lines by implementing a 4th order Runge-Kutta time integrator (in Python), with time-spacing h.We use a (spectral) spinweight decomposition to discretize the angular directions.The system is initialized at v = v 0 sufficiently late in time where we impose the final conditions given by θ (l) = H = 0.
In terms of spin-weighted spherical harmonics, we express the main variables as: ,m (x A ), ,m (x A ). (5.1) The spatial derivatives, which we already have written in terms of spin raising and lowering operators, take a very simple form since [49]: ,m .
(5.2) Then, the system of equations becomes a system of ODEs for each angular mode of each variable.The coefficients in that system are given by integrals of three spin-weighted spherical harmonics, which can be efficiently computed from the 3j symbols: where s 1 + s 2 + s 3 = 0. We use the implementation of the 3j symbols in the sympy package [50] to efficiently compute these integrals.The angular integrals on the sphere required to construct scalars such as the area are calculated using the quadpy package, and we use the numba package [51] to speed-up the simulations by writing optimized machine code.The complete implementation makes use of the numpy framework [52].
Finally we consider a source term describing an in-coming perturbation pulse followed by a decay motivated by quasinormal modes of a Schwarzschild black hole [44].This is given, in detail, by where we introduce an auxiliary function 5.5) in order to make the source term C 1 .Here v S is the impinging time of the source, which for simplicity we will always set to zero v S = 0, σ the width of the Gaussian, ω and τ the real and imaginary part of the quasinormal mode frequencies, respectively, and s ,m the excitation coefficients.The explicit source profile is shown in figure 1.Notice the source chosen interacts with the black hole on a timescale commensurate with its size, as a result one is not in the adiabatic regime and modes absent in the source can be populated by non-linear interactions (see e.g.[53]).
The system of equations present two practical numerical challenges.The first relates to their character which contains both first and second spatial derivatives which constrain ∆t rather strongly.The second one relates to the "aliasing problem", which can contaminate relevant low modes of the solution with high modes aliased down to lower modes.To address this, we adopt the standard "two-thirds" rule.That is, we regard the highest 1/3 max modes as "superflous" in a sense, zero-ing them at the end of each time step.
We observe that the Eqs.(4.11) are only non-linear in the first order in perturbations variables {h, c}.For these variables, only the (2, ±2) and the (3, ±3) modes are excited, since those are the modes directly excited by the source.As a consequence, only modes coming from quadratic terms in these variables will be excited in our problem: the higher modes will appear with = 6.We will fix then max = 6, and include an additional 3 ghost modes to avoid aliasing errors.We have checked explicitly that increasing the number of "physical" modes does not change the solution.
We also check the convergence of the time-integrator in figure 2. We compute the Cauchy convergence order through r = log 2 ( ϑ 2h − ϑ h 2 / ϑ h − ϑ h/2 2 ) and the result obtained, consistent with 4th-order convergence, is shown in figure . 2 Most of the simulations presented here use an initial step-size of h = 0.02M , and excitation coefficients given by s 2,±2 = , s 3,±3 = 0.5 × .The source is convoluted with a Gaussian profile to smoothen it out, attaining its maximum amplitude at v S = 0 with standard deviation σ = 5M .

Numerical Results
The first quantity to note is the evolution of the area, which to second perturbative order is given by Now we recall that the dynamics of ϑ is governed by the evolution of the first order terms {h, c}.The solution to equations (4.10) have exponential growing modes in both directions of evolution.Since we here demand θ → 0 in the asymptotic future, such modes are only relevant to the past.Their growth, in particular, will break the approximation adopted here, so we define a "stoppage" time v 0 as the time where the area of the black hole is 80% of its asymptotic (future) value.figure 3 illustrates this behavior.
To better appreciate the resulting dynamics, we show the mode-by-mode contribution to the total Hajiceck field H = h+ H and the metric perturbation C = c+ C in figure 4, as well as the evolution of the expansion ϑ (dual to the energy density of the fluid) in figure 5. We observe modes with different ( , m) than the ones directly excited by the source are populated.This behavior is a consequence of the non-linear character of Eq. (4.11).Note that close to our stopping the evolution, non-linear modes become comparable to main (linear) perturbing ones.
We measure the decay rate τ ,m for each mode both for the Hajiceck and the angular metric perturbation degrees of freedom.One would expect a clear genealogy of the modes, where the presence of a (2, 2) mode in the linear perturbation h to H induces quadratic contributions, e.g. a (4, 4) mode.This genealogy, in turn, manifests in the decay rates.For instance, τ (4,4) = 2τ (2,2) (see also [39,[54][55][56]).Indeed, we show in figure 6 that this is the case; the decay rates for extracted modes are consistent with the intuitive expectation with the exception of the (4, 0) mode.To confirm the obtained values are largely insensitive to the time-window used for the rates extraction, we vary it v ∈ [v S , v ], where v S is the time at which the source attains its maximum v the final time included in the linear regression.Figure 7 illustrates the decay times measured are consistent in such window with the exception of the m = 0 mode.This is related to the fact that the observed τ H (4,0) = 2τ (2,2) .Notice also the the non-linearly excited are essentially in lockstep with the linear ones.
For further scrutiny, we extract the amplitude of different modes vs. time and their  dependency on the strength of the perturbation ( ). Figure 8 illustrates some relevant modes of the Hajiceck field H extracted at 3 different times.Interestingly, at the order we work, there is no qualitatively distinct behavior in the scaling of the non-linear modes in the "decaying" stage of the perturbation considered as the strength of the perturbation is changed.Also, we note that higher values of m have a somewhat larger amplitude.At even earlier times however, stronger perturbations imply a more rapid failure of the perturbative expansion as already evident in figure 3. We can also connect this behavior where we have taken as characteristic wavelength L = 1/ for each mode, and recalled that H = h + H and C = c + C are the re-summed versions of the perturbations to the Hajiceck vector field and the angular metric.This Reynolds number also shows an exponential behaviour, and becomes larger than Re > 1 prior to the stop the evolution as shown in figure 9.This also indicates non-linear effects become important during the  evolution.
The impact on non-linearities also reflects in changes to vorticity (ω AB ) structure.Figure 10 illustrates how at late times the (spin weighted) vorticity field weakens, and its angular structure is dominated by the (2, 2) mode.However at earlier times it becomes stronger and the presence of higher modes induces a more complicated angular structure.This behavior is akin to the one observed in perturbed Schwarzschild AdS black holes within studies of fluid-gravity correspondence [48].Throughout the evolution total integral of the vorticity S 2 ω AB = 0 vanishes, as noted in appendix B.
Finally, for a more graphical way to visualize the impact of these non-linear features, we show an explicit embedding diagram of the horizon in figure 11.In order to do so we find an isometric embedding of the metric q AB + C AB in 3-dimensional space, assuming that 1.At late times the horizon becomes gradually more spherical, as imposed by the teleological conditions.At early times, however, it shows features with shorter wavelengths (high modes).As mentioned, our perturbative scheme is not able to fully capture the earlier dynamics of these modes.However, it is clear that for the kind of perturbations here considered (perturbations with a characteristic timescale commensurate with the black hole crossing time), these short-wavelength features naturally arise.A full non-linear analysis adapted to the null surface of the horizon, akin to the one proposed in [35] could be ideal to further explore these observations.As a final comment, we note the non-linear modes generated depend non-trivially on the perturbation considered, as the number of "channels" for non-linear interactions depend on this.For instance, when the perturbation is just a (l, m) = (2, 0) mode, the most relevant mode non-linearly generated is the (4, 0) though it is appreciably smaller than, e.g. the (4,4) mode generated with an equal amplitude combination of a (2,2) and (2,-2) perturbation at similar stages of the evolution.In particular the value of Re is < 1 even up to the point when we stop our simulation in the former, while it becomes > 1 in the latter.

Final Comments
In this work we have taken first steps to explore dynamical aspects of black holes and Carrollian fluids with a two-fold goal: On one hand, to explore Carrollian dynamics in settings relevant to event horizons.In this context, we identified several salient features: (i) the fluid evolves in a time-dependent domain that is affected by the fluid's dynamics, and higher order (in derivatives) contain both sources of dissipation and energy, the latter connected to the behavior of the embedding space.(ii) thermalization is linked to teleological boundary conditions.Both a negative bulk viscosity, and a "wrong" sign of dissipative terms in the evolution of the Hajicek flow appear as dooming the dynamics.However, this is just a manifestation of the "instability" associated to null surfaces, which typically will either expand to infinity or contract towards the singularity.Certainly, the above considerations imply that the insights drawn for Carrollian fluids in the gravitational context arguably apply to a rather special, sub-class of fluids with such a symmetry. 5Last, and as noted in the context of fluid-gravity correspondence, the full Einstein equations can potentially provide guidance for a completion of Carrollian hydrodynamics, in the regime of an "open system" with energy exchange with the embedding spacetime and backreacting on the volume in which it evolves.
On the other hand, to explore concepts of hydrodynamics and how they can motivate new insights on the non-linear behavior of horizons.In particular, the connection between rate of energy change, vorticity and the flow of energy towards IR/UV regimes; and definition of Reynold's number.In the former, we find a relation can be established though it is unclear what message can be derived from it as it is certainly different from the one found in, e.g.Navier-Stokes equations and standard results can not be assumed [58].For the latter, we can use similar qualitative arguments as in the fluid case to motivate a Reynolds number for perturbed black holes.However, potential connections with phenomena understood in hydrodynamics must be taken with care.Negative bulk viscosity indicate instabilities on the fluid side while in the context of black hole horizon such behavior is ruled out through teleological boundary conditions.As well, the rate of change of the Hajicek field's energy is governed by a significantly distinct expression than what is found in hydrodynamics.Of course this need not be surprising; here the "fluid" backreacts and defines the volume in which it exists and Carrollian symmetries are different from those in standard and relativistic hydrodynamics.Our first steps here have been focused on exploring first connections for future work in this area.
From a broader point of view, we were also motivated by exploring non-linear gravitational phenomena, and its impact/imprint on dynamical black hole horizons.Our definition of Reynolds number seem to capture when the non-linear regime becomes relevant.In addition, our solving for the black hole dynamics backwards in time, provides a rather natural way to assess such a regime.In particular, our perturbation included only leading modes and integration to the past highlighted the generation of non-linear modes and how they become commensurate with the perturbing ones.This behavior is particularly relevant in discussions on the onset of linear behavior in black hole decay to equilibrium, and the role of different mechanisms and different seeds inducing non-linear mode generation, e.g.[39,45,53,56].This is the Carrollian derivative appearing in the fluid equations (3.5)-(3.7).An interesting property is that it acquires an effective torsion when acting on scalar fields: We refer to AB as the Carrollian vorticity.If we understand the vector field b A as defining a Carrollian frame, then the vorticity associated to that frame would be naturally defined as where ϕ A is the acceleration, given by (A.9) B On vorticity and enstrophy

B.1 Energy of Hajicek field
In analogy to relevant discussion in fluids described by the Navier-Stokes equations, let us define the following "energy" associated to the vector field H A .To leading order, Its rate of change is given (to leading order) by (B.2) In the case of Navier-Stokes equations, the rate of change of the energy of the fluid flow can be related to the enstrophy, defined in terms of the flow's vorticity w AB as In the context of horizon's dynamics, two proposed definitions for the vorticity have been presented.One coincides with the imaginary part of the Newman-Penrose scalar Ψ 2 [46].The other one, proposed in [60] includes a further contribution so as to make it Carrollcovariant.However, both expressions coincide at order 2 .For completeness, these are given by w F S = ∂ Plugging this into the evolution of the energy density ρ yields: This bound is not akin to the one found in the Navier-Stokes case.For instance, in the incompressible case, one arrives at ρ = −2νZ (with ν the kinematic viscosity).Could this be interpreted as the fluid having a negative kinematic viscosity?An analog expression can be derived for an energy associated to the Carrollian field π A .The implications of these relations are yet unclear, their consequences will be explored elsewhere.

B.2 Vorticity conservation
Notice that the Hajiceck equation (3.4) implies that Thus, the integral of the vorticity remains constant; since we adopt H A = 0 asymptotically, it implies this integral vanishes.We have checked this condition is satisfied in our numerical solution.
Last, w AB is not a Carrollian 2-tensor, but a Carrollian vorticity can be defined, The conservation law (B.7) is still satisfied once every derivative is promoted to its Carrollian counterpart, and this quantity has the additional advantage of being the ultrarrelativistic limit of the vorticity of a relativistic fluid [34].

Figure 1 .
Figure 1.Source modes with excitation coefficients s 2,±2 = 2s 3,±3 = 10 −10 .The source is asymmetric in time, since the initial growth is governed by a Gaussian pulse whereas the late time behavior by quasinormal modes.

Figure 2 .
Figure 2. Convergence order measured with h = 0.04M, 0.02M, 0.01M and max = 6 with an additional 3 ghost modes.The result is consistent with the expected 4th-order rate.

Figure 3 .
Figure 3. Evolution of the area of the black hole (normalized with respect to the area at the time where the evolution starts v = 70M ) for different values of .The perturbative scale is controlled by the amplitude of the excitation coefficients of the source s ,m .

Figure 4 .
Figure 4. Modes evolution i the dynamical variables H = h+ H (left) and C = c+ C (right).The dashed lines correspond to the modes that are directly excited by the source ((2, ±2) and (3, ±3)) whereas the solid lines correspond to modes excited to second-order.The perturbative scale is = 10 −12 in this case.

Figure 5 .
Figure 5. Modes evolution of θ.The perturbative scale is = 10 −12 in this case.The (0, 0) mode drives the black hole's area change.

Figure 7 .
Figure 7. Dependence of the decay time τ H ( ,m) as a function of the final time of the fit v , for different modes, as indicated in the legend.

Figure 9 .
Figure 9. Evolution of the Reynolds number defined in Eq. (6.2) for different perturbative scales.In all cases we see that the Reynolds number becomes larger than 1 before we stop the evolution.

Figure 10 .
Figure 10.Evolution of the angular structure of the vorticity 2ω = q A qB ω AB as a function of the angular coordinates (θ, φ) at different times.The colormap (not normalized at different times) represents the value of the vorticity at each point.

Figure 11 .
Figure 11.Isometric embedding diagram of the induced metric in the horizon at different times.The colormap represents the radius of the embedding sphere (normalized).Short wavelength perturbations become relevant at early times.