Quasinormal modes of a semi-holographic black brane and thermalization

We study the quasinormal modes and non-linear dynamics of a simplified model of semi-holography, which consistently integrates mutually interacting perturbative and strongly coupled holographic degrees of freedom such that the full system has a total conserved energy. We show that the thermalization of the full system can be parametrically slow when the mutual coupling is weak. For typical homogeneous initial states, we find that initially energy is transferred from the black brane to the perturbative sector, later giving way to complete transfer of energy to the black brane at a slow and constant rate, while the entropy grows monotonically for all time. Larger mutual coupling between the two sectors leads to larger extraction of energy from the black brane by the boundary perturbative system, but also quicker irreversible transfer of energy back to the black brane. The quasinormal modes replicate features of a dissipative system with a softly broken symmetry including the so-called k-gap. Furthermore, when the mutual coupling is below a critical value, there exists a hybrid zero mode with finite momentum which becomes unstable at higher values of momentum, indicating a Gregory-Laflamme type instability. This could imply turbulent equipartitioning of energy between the boundary and the holographic degrees of freedom in the presence of inhomogeneities.


Introduction
Semi-holography introduces a way to model the complex dynamics of quantum field theories with asymptotic freedom. In this approach, one combines the perturbative description of weakly self-interacting ultra-violet degrees of freedom with a holographic description for those at lower energies and strongly self-interacting. This approach allows for flexibility in phenomenologically modelling open quantum systems, particularly in the case of a weakly self-interacting system coupled to a strongly self-interacting quantum critical bath [1][2][3].
The crucial ingredient in the semi-holographic construction is the "democratic" coupling between the two sectors which allows one to extract the low energy dynamics of the full system from the effective dynamics of the subsectors at any scale [4,5]. In this coupling scheme, the subsectors are subject to marginal/relevant deformations in their couplings and effective background metrics which are determined by the local operators of the other sector such that there is a local and conserved energy-momentum tensor of the full system in the physical background metric. The full dynamics can be obtained by solving the dynamics of both systems self-consistently in an iterative procedure [6][7][8]. Note that although one considers both subsectors at any scale, it is expected that the perturbative sector should dominate the ultraviolet behavior while the infrared behavior will be governed by the dynamics of the dynamical black hole horizon of the holographic sector. 1 Recently, the hydrodynamic attractor of such a hybrid system has been constructed by simplifying the description of both sectors to fluids [9]. It was found that the ratio of the energy densities of the strongly self-interacting to the weakly self-interacting sectors universally diverges as one approaches early proper time in Bjorken flow, confirming the expectation borne out of studies of perturbative QCD [10] that suggests such a bottom-up thermalization scenario. In explicit numerical simulations involving a black hole, it was found that there is irreversible transfer of energy from both the perturbative sector and the mutual interaction energy to the black hole whose apparent horizon grows monotonically. The rate of this irreversible transfer is very slow when the coupling between the subsectors is weak [8].
In this work, we present a simplified model of semi-holography which allows us to investigate the low energy dynamics and thermalization via the study of quasinormal modes (QNMs). In particular, we obtain robust understanding of why homogeneous thermalization of the full system can be parametrically slow, and how one can have inverse transfer of energy from the holographic to the weakly coupled sector at intermediate stages as observed in the simplified context of the hybrid hydrodynamic attractor. Another crucial question to which we gain new insights is whether one can have a mechanism in which there can be equipartitioning of energy between the holographic and the perturbative systems. We find instabilities that can lead to such a possibility in the presence of inhomogeneities.
The presence of a weakly broken symmetry plays a crucial role in our setup. We are also able to make connections with the broad literature of applicability of the quasihydrodynamic paradigm [11,12] in such systems [13] along with bounds on transport coefficients [14][15][16][17][18][19]. The novel feature of our model is that in addition to a diffusive Goldstone mode and a quasi-hydro mode, a third mode that is purely imaginary for small mutual coupling plays a crucial role in low energy dynamics creating potential instability towards turbulent/glassy dynamics.
The plan of the paper is as follows. In Section 2, we introduce our simplified semiholographic model explicitly detailing our motivation. In Section 3, we provide the technical details of the methodology of computation of QNMs. We then present our results for the homogeneous and inhomogeneous QNMs, and discuss their implications. In Section 4, we present the homogeneous non-linear dynamics of our model and investigate the parametrically slow thermalization. In Section 5, we discuss the questions raised by our results. In the appendices we provide supplementary information about the numerical implementation.

The simplified semi-holographic model
The construction of the semi-holographic framework is based on the following principles: 1. The non-perturbative part of the dynamics can be described by a strongly coupled large N holographic theory.
2. The interactions between the perturbative and holographic sectors can be described by the marginal and relevant deformations of the respective theories. The marginal and relevant couplings, and the effective background metric of each sector are promoted to ultralocal algebraic functions of the operators of the other sector in such a way that the full system has a local and conserved energy-momentum tensor in the physical background metric [4,5]. This coupling scheme is called democratic coupling [4].
3. The full dynamics should be solved self-consistently. This can be achieved via an iterative procedure in which the dynamics of each system is solved with couplings and the effective background metrics set to the values in the previous iteration until convergence is reached. The initial conditions should be held fixed throughout the iteration [6].
Numerous non-trivial examples explicitly demonstrate that the iterative procedure converges [7,8] with the present article serving as another such demonstration.

Review of the semiholographic glasma model
As an illustration, we consider a marginal scalar coupling between classical Yang-Mills theory and a holographic strongly coupled large N conformal gauge theory dual to classical -3 -Einstein's gravity with a negative cosmological constant. The dynamics of this system have been studied in [8] to gain insights into the possible non-perturbative dynamics of the color glass condensate, i.e. for understanding how the strongly interacting soft sector affects the initially overoccupied gluons at saturation scale which can be described by the classical Yang-Mills field equations.
The action for the full system in d spacetime dimensions is where χ(x) represents the deformation of the Yang-Mills coupling and h(x) is a source for a marginal operator in the holographic conformal gauge theory (CFT) with W CFT being the logarithm of its partition function. The inter-system coupling, β, has mass dimension −d. The equations of motion for the auxiliary fields χ(x) and h(x) lead to Substituting the above back in the action (2.1), we obtain (2.4) Finally, the holographic correspondence defines W CFT via where S grav is the renormalized on-shell action of the dual (d + 1)-dimensional gravitational theory which can be taken to be simply Einstein's gravity coupled to a massless dilaton field Φ dual to the CFT operator H(x). The non-normalizable mode of the dilaton, which specifies its boundary value (the boundary of the bulk spacetime is at r = 0), is identified with the source h(x) that couples to the dual operator H(x).
It is easy to see from the action (2.4) that the full system has a conserved energymomentum tensor which takes the form where t µν YM is the energy-momentum tensor of the Yang-Mills theory and T µν is that of the holographic sector. The full dynamics of the system was solved in [8]. Convergence -4 -was reached typically in 4 iterations as demonstrated by ∂ µ T µν = 0 being satisfied for all time to a very good accuracy. 2 We will present a simpler version of this set-up below and also the equations of motion explicitly. At this point, it could be mentioned that the most general democratic scalar couplings were found in [4].
The above model demonstrated that the energy in the Yang-Mills sector gets transferred completely to a growing black hole in the bulk holographic geometry for homogeneous initial conditions even if the bulk geometry is initially empty, i.e. a vacuum anti-de Sitter space with a vanishing dilaton. 3 At late times, both t µν Y M and the interaction term in the total energy-momentum tensor (2.7) decayed. The rate of transfer of energy to the black hole, as mentioned in the Introduction, was controlled by the mutual coupling β and was very slow for small β. It is to be noted that the model does not have any linear-coupling to the gauge field in the final thermal state. Since tr(F 2 ) vanishes at late times, the coupling to the bulk dilaton can only be quadratic in the gauge field. This does not allow us to relate the transfer of energy to the black hole to a quasi-normal mode of the full system easily. This motivates the simpler construction below.
The democratic coupling allows flexibility of constructing phenomenological models that capture the low energy dynamics of the full system based on effective descriptions of the subsystems. This has enabled understanding of the hydrodynamics of the composite system based on effective metric coupling of two fluids [5] and a preliminary study of the hybrid hydrodynamic attractor [9]. However, it is important to retain the dynamical black hole for capturing the infrared dynamics of the holographic sector in order to obtain the late time behavior and understand thermalization of the full system.
The effective metric coupling, unlike the scalar coupling discussed above, leads to hybridization of the thermal fluctuations of the black hole and the perturbative system. To understand hydrodynamization and thermalization in semi-holography, we should study the hybrid system of a gas of gluons described by kinetic theory coupled to the black hole by an effective metric coupling. The linearized hybrid modes were studied in the simpler version in [5] in which the black hole was substituted by a fluid. Here, we retain the black hole, but replace the kinetic theory by a massless scalar field, and the effective metric coupling by a linear scalar coupling. We find that the resulting simplified model retains many characteristics of the more complex models explored so far and gives several new insights. 2 The pile-up of numerical error breaks down the accuracy at very large time, however we could successfully extract the nature of the large time behavior. In particular, we were able to show the complete transfer of energy to the growing holographic black hole with both t µν Y M and the interaction term in the total energy T µν vanishing asymptotically at large time. 3 It was shown in [8] that we can start with such vacuum initial conditions by taking a suitable numerical limit in which the mass of an initial seed black hole is sent to zero.

Novel scalar semiholography
The simplified model introduces only a massless gauge-invariant scalar field at the boundary that couples to a black hole in the bulk via a (massless) bulk dilaton. Since such a field is gauge-invariant, we can couple it linearly to the holographic system unlike the gauge field which can couple only via tr(F 2 ), the energy-momentum tensor, etc. Therefore, instead of (2.4), we can consider the following action: with h(x) being the source of a marginal operator H(x) of the CFT. Note that here the inter-system coupling β has mass dimension −(d − 2)/2, different than the one introduced in (2.1). Furthermore, where S grav is the renormalized on-shell action of the dual (d + 1)-dimensional gravitational theory with a dilaton Φ whose boundary condition is given by (2.6). It is easy to see that the energy-momentum tensor of the full system is is the energy-momentum tensor of the boundary scalar field and T µν is that of the holographic CFT. It will be shown below that the explicit equations of motion of the full system directly implies the conservation of the full energy-momentum tensor with respect to the physical background metric η µν . Remarkably, the full energy-momentum tensor, unlike (2.7), is simply the sum of those of the two subsystems without an explicit interaction term. This feature will be helpful for us to deduce the dynamical consequences of the hybrid quasinormal modes. Linear semi-holographic couplings leading to such an energymomentum tensor of the full system have been explored in other contexts in [20,21]. 4 The equation of motion for the boundary scalar field from (2.8) is where we have used (2.2). The Ward identity (following from the diffeomorphism invariance) of the holographic theory implies that The linear couplings can be easily motivated also for fermions in the context of applications to condensed matter physics, since the boundary fermion is an electron which can be considered to be a gauge-neutral hadron made out of the partons of the holographic theory [1][2][3]22].
-6 -It is then easy to see from (2.10) and (2.11) that the full energy-momentum tensor given by (2.9) is indeed conserved, because Explicitly, the equations of motion for the metric and dilaton in the (d+1)-dimensional gravitational theory dual to the holographic strongly coupled large N CFT are Since all physical quantities of the gravitational theory should be measured in units of the AdS radius, 5 L, we set L = 1 for convenience. The generic solutions of the equations of motion have the following expansion in the Fefferman-Graham coordinates in which the radial coordinate is r, G rµ = 0, G rr = 1/r 2 and the boundary is at r = 0: The log terms above appear specifically for even d and capture the conformal anomaly of the dual CFT. To specify a unique solution, we need to specify the sources g µν (a.k.a. the boundary metric) and φ (0) aside from providing the initial conditions. For the semiholographic construction, these sources are determined by the gauge-invariant operators of the perturbative sector as discussed above. In the absence of effective metric couplings, so the boundary metric is the physical background metric. Furthermore, as implied by (2.8), The expectation values of the operators in the state of the CFT dual to the gravitational solution (determined now by the sources and the initial conditions) can be obtained from functional differentiation of the renormalized gravitational action [23][24][25]. The results are where X µν and ψ are local functionals of the sources of the theory, namely g (0) µν and φ (0) . The constraints of Einstein's equations imply two Ward identities, namely the conservation of T µν given by (2.11) and the trace condition, µν T µν = η µν T µν = 0. We will provide more explicit details in the ingoing Eddington-Finkelstein coordinates, which will be convenient for solving the dynamics numerically.
In what follows, we will explicitly analyze the hybrid quasinormal modes of this simplified semi-holographic system and study its non-linear dynamics. For numerical convenience, we will prefer to avoid the logarithmic terms in the radial expansion of the bulk fields. Therefore, we choose d = 3 for which we obtain the coupled system of a three-dimensional massless scalar field and a dynamical black hole with a dilaton field in AdS 4 .
From (2.10), (2.13) and (2.17) one can see that our simplified model has a global shift symmetry under which for a constant χ 0 . Note that this is a symmetry of the full non-linear theory. In the decoupling limit β = 0, the shifts of χ and Φ lead to two independent global symmetries. The coupling breaks these symmetries to the specific combination (2.20). Therefore, at finite but small β, the system has a quasi-hydro mode associated with a softly broken symmetry. It will be of fundamental interest to us as it will govern the relaxation dynamics of the system.
In what follows it is crucial that the diagonal shift-symmetry (2.20) is exact. Our model can therefore be interpreted also as a simple effective theory for a composite Goldstone boson interacting with a dissipative bath, where the underlying spontaneous symmetry breaking in the full (boundary plus bulk) system is not part of the model but takes place at a more fundamental level.
Generalized versions of our setup may help us to model the real time dynamics of the QCD axion and may be useful also for other phenomenological applications of holography. We will leave this for the future.

Quasinormal modes in semi-holography
QNMs are eigenfunctions of linearized perturbations which characterize the relaxation/growth of perturbations that govern the system away from a thermal equilibrium state. The thermal equilibrium state of the semi-holographic model discussed in the previous section is that of the bulk geometry being a static black brane with a vanishing or constant dilaton field Φ, while the boundary field χ also vanishes or is a constant. We will investigate whether this thermal background is stable against perturbations both linearly and non-linearly. In this subsection, we discuss how we can compute the hybrid quasinormal modes of this thermal equilibrium solution of the semi-holographic system, and save the nonlinear discussion for Sec. 4.1. For reasons mentioned before, we will consider a (2 + 1)-dimensional system (the -8 -bulk dual to the holographic sector is therefore (3 + 1)-dimensional). The description of our method will be modelled on the pedagogical account presented in [26] -we will highlight the crucial modifications brought in by the semi-holographic coupling. We also refer the reader to [27] for a comprehensive review of quasinormal modes of black branes.
The AdS 4 -Schwarzschild black brane dual to the thermal holographic sector takes the following form in the ingoing Eddington-Finkelstein coordinates where M is the Arnowitt-Deser-Misner (ADM) mass of the black brane and L is the AdS radius which we set to 1. The dual thermal equilibrium state has a temperature T = 3M where r = r h = M −1/3 is the radial position of the horizon.
Here we will focus on the hybrid fluctuations of the bulk dilaton and the boundary scalar field. At the linearized level, the metric perturbations will be exactly the same as those in the purely holographic case because the coupling to the boundary scalar is quadratic. The massless bulk dilaton field obeys the Klein-Gordon equation (2.13). With the background metric (3.1), a Fourier decomposition of the profile of the bulk dilaton according to Note that rotational symmetry of the black brane implies that the spectrum only depends on k = | k| 2 .
The bulk dilaton couples linearly to the boundary scalar χ, so it cannot be solved in isolation. It is more convenient to write the boundary equation in terms of the nonnormalizable mode which according to (2.17) should equal to −βχ(ω, k). The equation of motion for the boundary scalar field (2.10) in Fourier space can then be rewritten as where H(k, ω) can be obtained from the renormalized on-shell action (see Sec. 4.1 for more details) and is explicitly given by where φ (3) (ω, k) is the r 3 term of the near-boundary radial expansion of any solution of (3.4) that takes the form Physical solutions of the hybrid system corresponding to a causal response to perturbations must be ingoing at the horizon r = r h . A generic solution of (3.4) behaves as 8) and the ingoing boundary condition means setting c 2 = 0. The solution must also satisfy the semi-holographic boundary condition at r = 0 given by (3.5) and (3.6) which amounts to specifying a relation between φ (0) (k, ω) and φ (3) (k, ω), the two independent coefficients of the near-boundary radial expansion (3.7) of f (k, ω, r). For a given value of wave vector k, such solutions satisfying the boundary conditions at both r = 0 and r = r h can only exist for a discrete set of frequencies ω = ω i (k), which are the complex QNM frequencies.
We readily note that when β = 0, the boundary condition at r = 0 reduces simply to φ (0) = 0 as evident from (3.5) and (3.6). In this case, we get the usual conditions for the QNMs which require the linearized fluctuations to be normalizable. Since QNMs are intrinsic fluctuations of a system, we require them to exist source-free. This precisely implies φ (0) = 0 when the holographic system is decoupled from the boundary degrees of freedom. We thus reproduce the usual QNMs in the β → 0 limit, along with the ω = ±k modes of the decoupled boundary massless scalar field. Note that at finite mutual coupling β, the boundary conditions at r = 0 given by (3.5) and (3.6) still imply that the quasinormal mode fluctuations are intrinsic, i.e. they can exist without any external source. These equations impose the condition that the full hybrid system is not subjected to any external force.
The numerical method of determining the quasinormal modes by imposing both the ingoing boundary condition (3.8) at the horizon and (3.5) at the boundary has been discussed in details in Appendix A.

Homogeneous quasinormal modes
The homogeneous QNMs give us fundamental understanding of the relaxation dynamics of the system. In the decoupling limit, there exists two independent global symmetries, namely the constant shifts of the boundary and bulk scalar fields, which are broken to the specific combination (2.20) at finite value of β. Therefore, in the decoupling limit, there are two poles at the origin at zero momentum. At small β one of these is lifted but should be close to the origin. We call the latter the quasi-hydro mode ω Q for reasons to be discussed -10 -in the following subsection. Nonlinear simulations to be presented in Section 4 confirm that ω Q governs the homogeneous thermal relaxation of the full system. Explicitly, we find that at small values of β. The other (unlifted) pole which stays at the origin at k = 0 for any value of β will behave as a diffusion pole ω D ≈ −iDk 2 at small k and non-vanishing β. The diffusion constant D is negative above a critical value ofβ which is approximately 0.48 as discussed later.
Already the homogeneous quasinormal modes show a complex behavior. For the sake of notational convenience, we defineβ ≡ β √ T and use this dimensionless variable for the discussion. Note that the quasinormal modes will be of the general functional form In Fig. 1(a)-1(e), we have plotted the eight homogeneous (complex) quasinormal modes of the full system with lowest absolute values for various values ofβ. For small and nonvanishing values ofβ, there are three poles on the imaginary axis (aside from ω D which is at the origin for all values ofβ), namely 1. the quasi-hydro mode ω Q (plotted in red), which is parametrically close to the origin and is well approximated by (3.9) at smallβ, 2. a mode which we denote as ω G (plotted in green) that approaches the origin along the negative imaginary axis from −i∞ as the value ofβ is increased from zero, and 3. a mode which we denote as ω U (plotted in orange) that approaches the origin along the positive imaginary axis from +i∞ as the value ofβ is increased from zero.
In Fig. 1(b), corresponding toβ = 0.18, the red ω Q is slightly below the origin, whereas the green ω G and orange ω U are well separated from the origin on the negative and positive imaginary axes, respectively. Here, the poles in the decoupling limit [28], which have been plotted in Fig. 1(a), are shown again in gray color. The twin poles ω G and ω U clearly have no analogues in the decoupling limit. The remaining blue poles in Fig. 1(b) are on the lower half plane, and are only slightly displaced from their values in the decoupling limit shown in gray.
As evident from Fig. 1(b)-1(e), the unstable pole ω U stays on the positive imaginary axis for all non-vanishing values ofβ. It moves closer to the origin and attains a limiting value 1.03 × iπT asβ → ∞. This mode apparently implies an instability of the thermal state (corresponding to constant χ and Φ on a black brane geometry). In actuality, this only implies an instability over a short time scale as will be evident from our non-linear simulations presented in Section 4. Note, unlike the case of a closed system, a mode with Im ω > 0 may or may not imply instability in an open system. In our case, we -11 -  Fig. 1(a) shows the poles in the decoupling limit. Note that there are two poles at the origin corresponding to the independent constant shifts of χ and Φ. Fig. 1(b) shows that at small values ofβ one of the poles, namely ω Q (shown in red), is displaced slightly below the origin following (3.9), while two new poles ω G (shown in green) and ω U (shown in orange) appear from −i∞ and +i∞ respectively on the imaginary axis. Fig.1(c) shows that with increasingβ, ω G moves upward, while ω Q and ω U move downward on the imaginary axis. Eventually ω G and ω Q collide on the negative imaginary axis and transform into usual quasinormal modes as shown in Fig. 1(d). On the other hand ω U attains a limiting value on the positive imaginary axis asβ → ∞. As shown in Fig. 1(e), in the latter limit, all other QNMs (shown in blue) realign approximately on the same straight line on which the poles were located approximately in the decoupling limit but roughly at half-spacing.
-12 -do not have any instability in the homogeneous situation because of two reasons. The total conserved energy of the system shown in (4.15) is a sum of two non-negative terms, namely the boundary scalar kinetic energy and the black hole mass. Thus none of these can grow without bound in magnitude as they are bounded from both below and above. Furthermore, Birkhoff's theorem 6 guarantees that the homogeneous thermal state is the unique static solution where of course entropy cannot be produced. As the entropy given by the area of the apparent horizon grows monotonically (as explicitly verified in Sec. 4), the endpoint of evolution in the homogeneous case should be the static black brane.
So we can anticipate what is borne out by our non-linear solutions: for an arbitrary homogeneous perturbation about the thermal state, the unstable pole ω U governs the rapid transfer of some energy from the holographic sector to the boundary scalar field, which is followed by a slow, complete and irreversible transfer of energy back to the black hole over a timescale governed by Im ω Q at small β. The pole ω U thus does not signal an imminent transition to another phase (unlike e.g. the case of holographic superconductors [31]), only the propensity to process a perturbation in this particular way.
The pole ω G is associated with a Gregory-Laflamme type of instability at finite k as discussed later. 7 However, it remains on the lower half plane for all values ofβ at k = 0.
In Fig. 1(c), corresponding toβ = 0.4, we see that ω Q has moved down while ω G has moved up along the negative imaginary axis. At a sightly higher value ofβ, these poles collide on the negative imaginary axis, after which they move almost horizontally keeping the imaginary part almost unchanged as evident from Fig. 1 Thus ω Q and ω G are transformed to usual quasinormal mode poles for higher values ofβ.
Asβ is increased towards infinity, all poles at k = 0, except for ω D (which stays at the origin) and ω U (which goes to the limiting value on the positive imaginary axis), realign approximately on the same straight lines on the lower half plane along which the decoupled quasinormal poles were placed. Furthermore, in this limitβ → ∞, the poles are almost halfway in between the quasinormal mode poles of the decoupling limit. This is evident from Fig. 1(d) corresponding toβ = 10.
We provide a summary of the above discussion as a snapshot in Fig. 2. 6 The massless dilaton in our case has to be constant at the horizon for regularity, which implies it is constant everywhere for a static configuration. A constant dilaton has vanishing stress tensor and hence we obtain the AdS-Schwarzschild black brane solution. In the case of a holographic superconductor [29], there is a non-trivial potential for the scalar field and/or a non-trivial radial mass profile. The analysis of stable stationary configuration is more complicated but can be done via the method of Hollands and Wald [30]. 7 As we shall see, the fate of this additional pole when k is increased depends on the value of β. For sufficiently small β it always gives rise to a Gregory-Laflamme type instability, while at larger β the diffusion pole can take over this role. Then the label "G" simply stands for the color "green" in the plots. By contrast, ωU always refers to an unstable mode.

Quasinormal modes at finite momentum
The behavior of QNM frequencies with varying k at fixed β and T is even richer. Qualitatively different behaviors are observed for These are illustrated in Figs. 3-7. In the decoupling limit β = 0, we recover the propagating Goldstone modes ω = ±k of the boundary scalar and the usual complex quasi-normal modes of the bulk scalar at any value of k as mentioned before. As described below, even a small value of β changes the character of the Goldstone modes while two other non-trivial modes emerge as in the homogeneous case described previously.
The case ofβ 0.391: The representative case ofβ = 0.35 in this category is illustrated in Fig. 3. The pole ω D (shown in purple), which is at the origin at k = 0, becomes diffusive, Re ω π Im ω π (e) k = 2.54π Figure 3. Behavior of the first six QNMs in complex frequency plane with varying k atβ = 0.35.
As k evolves, the diffusion pole ω D shown in purple (and is at the origin at k = 0) collides with the quasi-hydro mode ω Q shown in red on the negative imaginary axis, and subsequently both of them transform to a pair of quasinormal modes with almost k-independent negative imaginary parts. The other mode ω G shown in green moves upwards on the imaginary axis, crosses the origin and collides with the downward moving ω U shown in orange above the origin. Subsequently both of them transform into a pair of unstable quasi-normal modes with almost k-independent positive imaginary parts.
i.e. it behaves as  at small k. We will discuss the dependence of the diffusion constant D onβ later. The quasi-hydro mode (shown in red) moves upwards along the negative imaginary axis with increasing k and eventually collides with ω D . We denote the value of k where this collision happens as k * , which is ≈ 0.555π for the chosen value ofβ. For k > k * , these two poles ω D and ω Q transform into a pair of QNMs which evolve almost horizontally (with constant imaginary parts). The momentum k * is critical because the dispersion relation of ω D giving the effective diffusive dynamics of the boundary massless field χ becomes non-analytic (with discontinuous first derivatives) at k = k * signalling the breakdown of an effective hydrodynamic description.
A representative characterisation of Re ω D (k) and Im ω D (k) forβ = 0.15 is presented in Fig. 4. The value of k * at which ω D and ω Q collide turns out to be 0.07π T , producing a pair of complex quasinormal poles which have almost k-independent imaginary part and with non-vanishing real parts Re ω D (k) ≈ ±|k − k * | δ as k → k * + (with > 0). We find that δ ≈ 1/2 and is independent of the choice ofβ in this range to a very good accuracy. As one would expect, Reω D (k)/k ≈ 1 as k → ∞. It is to be noted that we can get k * arbitrarily close to the origin by tuning β to smaller values. This feature that propagating modes with non-vanishing Re ω(k) exists for k > k * is called the k-gap [32]. In the context of phonon-like modes, this has been observed in [33,34].
Such phenomena of breakdown of hydrodynamics due to collision between a hydrodynamic and a quasi-hydrodynamic mode at a real and parametrically small value of momentum leading to a k-gap is a characteristic property of the collective modes of liquids [35]. In fact, our model interpolates between various qualitatively different behaviors about which we will have more to say in the following subsection.
The transformation of massless propagating Goldstone modes of the type ω = ±k -16 -(which is exactly how the modes of χ behave in the decoupling limit) to a pair of diffusion and a quasi-hydro modes has been observed before in other models of holography [36,37] where such modes are also obtained from an explicit and soft symmetry breaking [13,36,38] 8 and have been also discussed from an effective field theory point of view [40,41]. 9 For a recent review, see [42]. The independent shift symmetries are broken to the diagonal explicitly via the semi-holographic coupling in our model, while translation symmetry is explicitly broken in the models discussed in [42]. 10 However, in contrast to the latter models, ours can be thought of as an open quantum system with a finite total conserved energy if we consider the boundary scalar field as the system and the black hole with the dilaton field as the bath. Indeed our model has an additional feature, namely the presence of an additional modes ω G mode on the negative imaginary axis that also contributes to low energy dynamics. This is absent in usual setups. As evident from Fig. 3, the ω G mode moves upwards along the imaginary axis, and crosses the origin at a finite momentum k 0 . At k = k 0 , the system therefore has a Gregory-Laflamme type instability! Also when k is near k 0 , ω G cannot be excluded from the low energy description of the system.
As k is increased above k 0 , ω G collides with ω U (which moves downwards with increasing k) on the positive imaginary axis close to the origin, and both transform into a pair of unstable quasinormal modes with small imaginary parts as shown in Fig. 3. Subsequently, they move almost horizontally with almost k-independent positive imaginary parts.
The Gregory-Laflamme type instability [44] can have profound consequences for the dynamics of this system. Since the total conserved energy given by (4.15) is a sum of two non-negative terms, a repetition of our argument in the previous subsection, namely that the poles on the upper half plane could only lead to initial instabilities involving reverse transfer of energy from the black hole to the boundary scalar field, could have got through had there been no zero modes at finite k. The presence of a zero mode at k = k 0 may imply that the system may not be able to evolve to the static thermal configuration eventually and the final end point could be turbulent or glassy [45,46]. Unfortunately, non-linear simulation of this system in presence of inhomogeneities is difficult and we postpone this to a future work. In fact, unlike the usual Gregory-Laflamme instability of the black string, ours is intrinsically a (3 + 1)-dimensional gravitational problem.
The case ofβ ≈ 0.391: This is illustrated in Fig. 5. The collision of the diffusion pole ω D with the quasi-hydro pole proceeds as in the previous case. However, after these transform into a pair of stable quasinormal modes moving horizontally away from the imaginary axis with increasing k, they reverse back to the imaginary axis and once again 8 Remarkably, it has been shown in [39] that the k-gap is produced naturally also in holographic models with higher form fields. 9 In [36], actually the reverse transformation of a pair of complex poles to a pair of purely imaginary poles at higher values of k was reported. 10 A similar phenomenon of emergence of a diffusive pole in presence of time-translation symmetry breaking has been discussed in [43].
-17 -collide there. 11 Subsequently one of these poles moves downwards on the negative real axis, colliding with the pole ω G that moves upwards with increasing k, and thus producing two almost horizontally moving quasinormal modes. The other one moves upwards and produces the Gregory-Laflamme type instability as before. It is to be noted that for these values ofβ, there exists a region of value of k in which the three poles ω D , ω Q and ω G are on the negative imaginary axis and almost degenerate. 12 The case of 0.391 β 0.4425: The illustrative case ofβ = 0.4 is shown in Fig. 6. This is very distinct from the previous cases because the diffusive pole ω D never collides with the quasi-hydro pole. It initially behaves as a diffusion pole, but it starts moving upwards along the negative imaginary axis as the value of k is increased and eventually crosses the origin producing a Gregory-Laflamme type instability -there exists a finite value of k, namely k 0 , at which ω D (k 0 ) = 0 and Im ω D (k) > 0 for k > k 0 . The quasi-hydro pole ω Q moves downwards with increasing k, in contrast to the previous cases, and collides with the upward moving ω G pole, producing a pair of stable horizontally moving QNM poles.
The case ofβ 0.4425: The representative case ofβ = 0.5 is shown in Fig. 7. Firstly, even at k = 0, there exists no pole on the negative imaginary axis. We recall from the previous subsection that the ω Q and ω G poles are complex forβ 0.4425. Secondly, the diffusive pole ω D has a negative diffusion constant at small k forβ 0.48 and moves upwards on the positive imaginary axis to collide with the downward moving ω U pole. Consequently, forβ 0.48, there exists no finite value of k for which there is a quasinormal mode pole at the origin, and hence no Gregory-Laflamme type phenomena. We find that as we approachβ c ≈ 0.48 from below, k 0 scales like |β −β c | ρ with ρ ≈ 1/2. The negative diffusion constant could lead to clumping instabilities, which should be investigated via a numerical simulation of the inhomogeneous non-linear dynamics in future work. In the narrow range 0.4425 β 0.48, the diffusion constant is positive as in the previous case, but the value of k 0 at which ω D crosses the origin again moves towards zero asβ gets closer toβ c ≈ 0.48. It seems likely that the non-linear dynamics of the system is qualitatively different forβ ≥β c .

On the diffusion constant D and the Gregory-Laflamme momentum k 0
The mode ω D behaves as a diffusive mode at small k as discussed above. The plot of the dimensionless product of the diffusion constant times the temperature (DT ) as a function of the dimensionless mutual coupling (β) is presented in Fig. 8. We find that the diffusion constant D decreases monotonically withβ and changes sign atβ c ≈ 0.48 as discussed previously. Since the diffusive behavior exists for k k * (where k * is the momentum at 11 This second collision is closer to the phenomenon described in [36]. 12 It is possible that there exists a value ofβ around 0.391 on the negative real axis where the three poles ωD, ωQ and ωG coincide on the negative imaginary axis at a specific value k. -18 -    which ω D collides with ω Q ), and k * can be arbitrarily close to the origin for smallβ, it is very difficult to determine D at small values ofβ numerically.
It has been argued that the diffusion constant should satisfy an upper bound [17][18][19][47][48][49] in a wide class of many-body systems, i.e. 13 D v 2 * τ * , with v * = |ω * |/|k * |, and τ * = |ω * | −1 . (3.10) Above k * and ω * are the values of the momentum and frequency respectively at which the hydrodynamic description breaks down, i.e. they set the limits of the convergence of the gradient expansion. To be precise, k * is the value of momentum at which the hydrodynamic mode collides with a gapped mode or a branch point, and ω * is the value of the complex frequency at that point [12,[51][52][53]. Interestingly, k * can be complex and should be then determined by the analytic continuation of the hydrodynamic and nonhydrodynamic modes. It is expected that v * is essentially an effective state-dependent Lieb-Robinson velocity governing the ballistic growth of the operators at late time (see [54] and [17]). The inequality (3.10) is saturated typically in holographic theories and in other models such as SYK chains [16,55].
A lower bound on the diffusion constant has also been conjectured [14][15][16]19, 48] 14 to hold for many-body systems primarily inspired by the KSS bound [57] on η/s (which should be stated in terms of the product of the diffusion constant and the temperature more generally). In case of fermionic systems, v * has been identified with the Fermi velocity in fermionic systems [14] while the corresponding τ is the Planckian scattering time [14]. For holographic systems, the velocity is identified with the butterfly velocity v B and τ with the corresponding Lyapunov time τ L [15,16]. In the case that the dispersion relation of the diffusive mode is univalent over the entire complex z-plane with z ≡ k 2 except for a branch 13 See [50] for a discussion in the context of the Goldstone diffusivity. 14 See also [56].
-21 -point and at z = ∞, then according to [53], v * should indeed be the butterfly velocity. We will not have much to say about the lower bound in our model because we believe that we need an independent computation to establish the butterfly velocity and the Lyapunov time in our model (see below).
Forβ 0.391, the value of k * in our model is simply the (real) momentum at which the ω D collides with ω Q , and ω * = |ω D (k * )| = |Im ω D (k * )|. We find that indeed forβ 0.35, as shown on the right in Fig. 8. Thus the upper bound in (3.10) is half-saturated. In the regime 0.35 β 0.391, the stricter inequality (3.10) holds, i.e.
For instance, when β = 0.39/ √ T , we find that  Fig. 6 again), and 3. k * is the complex momentum at which ω D (k) collides with ω Q or another pole after analytic continuation.
Actually k * would be the smallest of these three possibilities.
In absence of an understanding of the analytic properties of ω D as a function of z ≡ k 2 , it is unclear how we can identify the butterfly velocity and Lyapunov exponent in our model; an independent computation following [58] could be necessary to settle this. At present, we cannot comment on the validity of the lower bound on the diffusion constant in our model. It is worth mentioning though that as shown in the section 4, the energy relaxation time is at smallβ. The inequality above follows from our previous discussion that for k < k * and smallβ, ω Q is purely imaginary and its imaginary part decreases with increasing k. This would imply that a lower bound similar to (3.10) but with the inequality reversed would be almost saturated if v B ≈ v * and τ eq ≈ τ L at smallβ.
It is not clear how a negative diffusion constant or even the vanishing of the diffusion constant could be reconciled with (3.10). 15 In the future, we would like to investigate this further.
The Gregory-Laflamme momentum k 0 at which ω G or ω D crosses the origin from the lower half plane also monotonically decreases withβ as shown in Fig. 9. As discussed in the previous subsection, k 0 goes to zero when the diffusion constant changes sign. Note that k 0 diverges in the limitβ → 0. Since in the latter limit ω G moves towards −i∞ at k = 0, it crosses the origin at higher and higher values of k 0 .
Another type of instabilities is known to appear in a weakly coupled plasma (and also glasma [59]) out of equilibrium: plasma instabilities involving gauge fields such as Weibel instabilities [60][61][62][63][64]. Those are clearly left out by our simplified model for the ultraviolet degrees of freedom. A full-fledged semi-holographic description of large-N c Yang-Mills plasmas should in principle also contain those. However, in the context of heavy-ion collisions it has been found that they tend to be too slow in their on-set to play a crucial role [63,65]. In particular, the bottom-up scenario of Ref. [10] may be qualitatively right even though it ignores plasma instabilities.
Let us also point out that plasma instabilities are qualitatively different from the instabilities we have found in the present simplified semi-holographic model. The former are present for a certain range 0 < k < k max with vanishing growth rate at k = 0, unlike the 15 The vanishing of the diffusion constant could be compatible with the ωG pole coming close to the origin at finite momentum for higher values ofβ. Also since open systems can have poles in the complex upper half frequency plane without implying instability of the thermal equilibrium state, one needs to reevaluate the bounds on transport coefficients arising from the analytic properties of the poles in such systems.
-23 -mode ω U . By contrast, the Gregory-Laflamme instabilities set in only above a nonzero value k 0 .
Before concluding this subsection, we would like to refer the reader to [66][67][68][69][70] for the computation of quasinormal modes with an examination of the diffusive behavior in purely holographic systems with a broken global symmetry. An analytic expression for the diffusion constant was obtained in [68] particularly in terms of thermodynamic data.

Emergence of conformality at infinite mutual coupling
In all our previous semi-holographic models, we have found emergence of conformality when the mutual coupling between the subsectors becomes infinite. At the level of quasinormal modes, the latter should imply that the quasinormal frequencies should behave as At fixed T and k, the above implies that all ω QNM should saturate to finite values in the limit β → ∞. We have indeed seen this feature in the case of k = 0 as shown in Fig. 2.
In Fig. 10, we have plotted the QNM poles at various values of β at T = 1 and k = 1.5. We readily notice that in the limit β → ∞ 1. all QNMs become independent ofβ saturating to finite values, and 2. they align themselves approximately on the same straight line as in the case of the decoupling limit, but placed approximately halfway between the latter poles (except for the two poles in the upper half plane).
We have verified that the above features hold at any value of k indicating the emergence of conformality at infinite mutual coupling.

Methodology of non-linear simulations
The full non-linear dynamics of the semi-holographic model can be readily computed based on the iterative procedure proposed in [6], and successfully demonstrated in the more complex case of classical Yang-Mills and dilaton plus black brane system in [8]. Here we implement a simpler version of [8] with more general initial conditions in the case of homogeneous non-linear dynamics.
-24 - The asymptotically AdS 4 metric representing the holographic sector can be assumed to be of the following form in the in-going Eddington-Finkelstein coordinates with the boundary being at r = 0. We -25 -can consistently set the anisotropy to zero, i.e. assume G xx = G yy . The bulk dilaton profile takes the form Φ(r, t) while the boundary scalar field χ(t) is only a function of time.
The gravitational equations of motion (2.13) reduce to the following nested set of partial differential equations: where d + is the directional derivative along the outgoing null radial geodesic, i.e.
Setting the boundary metric to be η µν and solving the equations of gravity (2.13) order by order in r near the boundary gives where denotes the time-derivative. Above, the residual gauge freedom r → r + f (t) of the Eddington-Finkelstein gauge is fixed by setting the coefficient of r in the expansion of A(r, t) to zero. The normalizable modes, namely a (3) (t) and φ (3) (t), remain undetermined in this procedure and need to be extracted from the full bulk solution that is determined by the initial conditions. The constraint of Einstein's equation implies that Performing holographic renormalization of the on-shell action and taking the functional derivative with respect to the sources, we obtain the expectation values of the energy momentum tensor and the scalar operator in the dual CFT state [8]. These turn out to be and respectively. Considering terms which are linear in φ (0) , the above reproduces the result for the linear fluctuations (3.6) when it is homogeneous. As a consistency check, we readily note that (4.11) reproduces the Ward identity of the CFT where we have used the key relation φ (0) = −βχ which sets the value of the non-normalizable mode in terms of the boundary scalar field. The equation of motion for the boundary field χ given by (2.10) reduces to (4.14) The conservation of the energy-momentum tensor of the total system (2.9) amounts to the total energy E tot , which is simply the sum of the boundary scalar field's kinetic energy and the ADM mass of the black hole, remaining constant. Indeed it is easy to verify using (4.13) and (4.14) that The non-linear dynamics of the full system is determined uniquely by the initial conditions for χ(t 0 ), χ (t 0 ), a (3) (t 0 ) and the initial profile of the bulk dilaton Φ(r, t 0 ) at the initial time t 0 . The iterative method of computing this numerically is discussed in Appendix B.

Results for the non-linear evolution of the homogeneous case
In this section we present our results for the non-linear simulations of generic but homogeneous initial conditions. As discussed before, each evolution is uniquely specified by the initial configuration of the bulk dilaton field Φ and the initial values of a (3) = −(1/2)E bh = −M (the ADM mass of the black hole), and χ and χ , i.e. the values of the boundary scalar field and its time-derivative. Note that due to the presence of the symmetry given by (2.20) we can set the initial value of χ (and therefore the boundary mode of the bulk dilaton field given by φ (0) = −βχ) to zero without loss of generality.
For the purpose of illustration, we consider two different types of initial conditions (ICs) at the initial time t = 0: IC1 : χ = 0, χ = 0, Φ(r) = r 5 e −r 2 , a (3) = −1;  These two initial conditions represent cases where the initial kinetic energy of the boundary scalar field is zero and non-zero, respectively. In both cases, the iterative method discussed -27 -  Figure 11. Boundary scalar field χ as a function of time for β = 0.1 and initial conditions given by (4.16). We see that χ slowly saturates to a constant value.
in Section 4.1 converges to a very good accuracy after four iterations, and the total energy (4.15) is conserved, i.e. is time-independent for a sufficiently long time which allows us to reliably draw our conclusions. A more detailed discussion about the numerical accuracy is provided in Appendix B.
For both initial conditions, as we have anticipated in Section 3.2, initially there is transfer of energy from the black hole to the boundary scalar field. This is followed by complete and irreversible transfer of energy to the black hole. When β is small, we also anticipated that the initial transfer of energy to the boundary should be rapid while the subsequent reverse transfer of energy back to the black hole should be slow. The final state is just the thermal state given by the black hole with a constant mass and with constant values of the boundary scalar and bulk dilaton fields. A plot of the boundary scalar field χ(t) with time is provided in Fig. 11 for β = 0.1 and initial conditions set by (4.16). The kinetic energy of the scalar field E kin (t) and that of the holographic sector E BH (t) = 2 − E kin (t) are provided in Fig. 12 with the same initial conditions note the total energy, E kin (t) + E BH (t), is 2). We find that the kinetic energy of the boundary scalar E kin fits perfectly to an exponentially decaying function, i.e. E kin (t) ≈ αe −γt with γ > 0 at late time.
From the QNM analysis, we can also anticipate γ quantitatively. At late time, we expect that where χ f is the final value of χ and γ Q is determined by the homogeneous (purely imaginary) quasi-hydro mode with the identification ω Q = −iγ Q , since ω Q is the QNM pole that is closest to the origin and is on the lower half plane. The kinetic energy of the scalar field -28 -   Higher values of β lead to larger energy extraction from the black brane by the boundary scalar field but then the energy is returned back to the black brane irreversibly and completely at a higher rate. The time t max at which the boundary scalar field attains its maximum energy E max is independent of β and is ≈ 2.9 as shown in the Table 1.
should then behave as E kin ≈ e −2γ Q t at late time and therefore γ = 2γ Q . Finally, since ω Q is given by (3.9) at small β, we obtain that γ ≈ 11.2πβ 2 T 2 (4. 18) with T determined by the final mass of the black hole (which is 1/2E BH ) via (3.2).
In Table 1, the values of α and γ have been computed for β ranging between 0.001 and 0.1 and initial conditions set by (4.16). First, we find that the value of γ satisfy (4.18) with remarkable accuracy. For the initial conditions set by (4.16), the final mass of the black hole should equal to its initial mass (which is unity) because the initial kinetic energy of χ is zero. Therefore, we obtain from (3.2) that the final temperature is 3/(4π) ≈ 0.239. From (4.18) with β = 0.1, we find γ ≈ 0.0201 which is precisely the value calculated in Table 1. The other values of γ are reproduced to the same degree of accuracy by (4.18). Fig. 13 confirms that γ and (4.18) scale like β 2 for fixed initial conditions. The time t max at which the kinetic energy of the boundary scalar attains its maximum value, E max , is independent of β. In Fig. 13 we see that E max also scales like β 2 for fixed initial conditions. We also note from Table 1 that The above scaling properties lead to a rather interesting result. For small β and initial conditions set by (4.16), we should have  Table 1. E kin , the kinetic energy of the boundary scalar, is fitted to function αe −γt from t = 10 to t = 18.5 for different values of β with the initial conditions set by (4.16). t max is the time when E kin attains its maximum value E max . The values of γ match with those predicted by (4.18) to a remarkably good accuracy. We note that t max is independent of β while α ≈ E max and both scale as β 2 . For fits see Fig. 13. We are using units in which the initial black hole mass is unity.
Since α and γ both scale as β 2 and t max is independent of β, we obtain that Thus the limit β → 0 is non-trivial. Furthermore, A is independent of β for small β to a very good approximation. This result implies that if the boundary system draws more energy from the black hole (bath), then it has to give it back to the black hole (bath) at a higher rate. It would be interesting to see if (4.21) could be a generic feature of out-ofequilibrium open quantum systems with the value of the limit possibly depending on the initial conditions.
The area of the apparent horizon acts as a proxy for the entropy of an out-of-equilibrium semi-holographic system as noted in [8]. We indeed find that the entropy grows monotonically although the black hole mass is non-monotonic as a function of time. In Fig. 14, we have plotted the radial position and area of the apparent horizon as a function of time for various values of β and initial conditions set by (4.16). We observe that the entropy saturates to its final thermal value more quickly for smaller values of the coupling β. Fig. 15 indicates that althoughĖ kin /E kin ≈ γ is setting the rate of energy exchange between the subsystems, the rate of growth of entropyṠ AH /S AH decays to zero at late times. This suggests that asymptotically there is an isentropic transfer of energy between the boundary scalar and the black hole. In fact, this is a proof of thermal equilibration because the latter implies that the rate of growth of entropy should vanish at late time.
We observe the same qualitative features for simulations of the system with initial  Table 1 for the corresponding β values in the table. Red dashed line is fitted function aβ 2 +b with a = 0.39445 and b = 3.791×10 −6 Figure 13. Scaling of α, γ and E max with semi-holographic coupling β for initial conditions set by (4.16). We are using units in which the initial black hole mass is unity.
conditions set by (4.17) in which the initial energy in the scalar field is non-vanishing (see Fig. 16). Of course γ, the rate of decay of the scalar kinetic energy at late time, is quantitatively the same as well. Furthermore, the maximum transfer of energy to the scalar sector, i.e. E max − E kin (t 0 ) with t 0 being the initial time, scales as β 2 and t max is almost independent of β. The rate of the monotonic growth of the entropy of the system goes to zero at late time confirming thermal equilibration. We have found that indeed that these features are present for generic initial conditions while the final mass of the black hole is simply determined by the fact that it is equal to the total initial energy of the system. Only the transient behavior at initial time depends on the initial conditions. Parametrically slow thermalization in non-conformal holography had been observed in [71,72]. However the mechanism discussed in these works was not related to a quasi-hydro mode but rather to the presence of two intersecting branches of black brane solutions.
-31 -  Although the apparent horizon first moves away from the boundary (at z = 0) and then towards it, the entropy grows monotonically.

Conclusions and outlook
Our simple semi-holographic model is fundamentally an open quantum system involving one preserved and one weakly broken global symmetry. For generic homogeneous initial conditions and weak inter-system coupling, an unstable mode implies rapid transfer of energy from the holographic sector to the massless field at the boundary, while the purely imaginary quasi-hydro mode governs a slow, irreversible and complete transfer of energy to the black brane at later stages. The entropy of the system, represented by the area of the apparent horizon of the black brane, grows monotonically although the black brane mass behaves non-monotonically. Higher values of the inter-system coupling leads to more extraction of energy from the black brane by the boundary scalar field, but then a quicker -32 -   Higher values of β lead to larger energy extraction from the black brane by the boundary scalar field but then the energy is returned back to the black brane irreversibly and completely at a higher rate. The time t max at which the boundary scalar field attains its maximum energy E max is independent of β. Thus the qualitative features are exactly like in the case of the initial condition (4.16) shown in Fig. 12.
irreversible and complete transfer of energy back to the black brane. Furthermore, the integral of the kinetic energy of the boundary scalar field with time for a vanishing initial value, remains finite even in the limit when the inter-system coupling vanishes. This feature deserves a more basic understanding with insights from non-equilibrium statistical mechanics.
The inhomogeneous dynamics is even richer. We find that for any value of the intersystem coupling, the mode at the origin becomes diffusive at finite momentum. At weak inter-system coupling, the quasi-hydro mode also moves up on the negative imaginary axis and collides with the diffusion pole producing a pair of complex poles as the momentum is increased. This results in the system having low energy propagating modes for momentum k > k c where k c is close to zero for small inter-system coupling. This feature is quite ubiquitous in dissipative systems with a softly broken symmetry and is called the k-gap [32,35]. Additionally, our system has a third pole which also moves from negative imaginary infinity along the negative imaginary axis with increasing momentum and produces an instability as it crosses the origin. This is similar to the Gregory-Laflamme instability [44]. At intermediate values of the coupling, these three poles are approximately degenerate for an intermediate range of momentum. At higher values of the inter-system coupling, the diffusion pole (instead of the third pole) reverses back and crosses the origin again as the momentum is increased producing the Gregory-Laflamme type instability, while the quasihydro mode collides with the third pole on the negative imaginary axis. At even higher values of the inter-system coupling, the diffusion constant of the diffusive mode becomes negative, and no pole crosses the real axis from the lower half plane at any value of the -33 -momentum. The Gregory-Laflamme momentum, at which one of the three modes has zero energy, exists only below a critical value of the inter-system coupling.
The model thus exhibits diverse behavior for different values of the inter-system coupling. Since the total conserved energy of the system is simply the sum of two non-negative terms, namely the kinetic energy of the massless scalar field at the boundary and the black brane energy, we have argued that in absence of zero modes with finite momenta, the unstable poles only imply short-term instability involving inverse transfer of energy from the holographic sector to the boundary scalar field. The dynamics is constrained by the facts that the energies cannot grow indefinitely and that the entropy represented by the total area of the apparent horizon should increase monotonically. The rate of growth of entropy becomes zero typically when the system reaches the thermal state represented by the static black brane geometry. However, there is no guarantee that the system has an entropy current generically although there exists a global monotonic entropy function. It is likely that the presence of zero modes at finite momenta can lead the system to turbulent or glassy final states even in the presence of an entropy current and the fate of the instability is also describable by a quasihydrodynamic theory. 16 In the former case, one may achieve equipartition of energy between the perturbative and holographic sectors at least for lower values of the inter-system coupling. We leave this issue for future investigation.
It would also be of interest to find an appropriate causal and consistent quasihydrodynamic effective theory 17 which can describe the interplay of the three modes that play a role in the low energy (macroscopic) dynamics of the system. Furthermore, we would like to have a better understanding of whether the diffusion constant saturates or satisfies conjectured bounds beyond those values of the inter-system coupling reported here. The latter would require us to study the applicability of quasihydrodynamics [13,40] in this system in more detail, and furthermore study the Lyapunov exponent and the butterfly velocity. This discussion needs to be reconciled with the change of sign of the diffusion constant with increased inter-system coupling. We leave this for the future.
We would also like to investigate whether the usual effective metric and scalar couplings of semi-holography [4,5] have similar features to the simple linear inter-system coupling described here. The numerical simulations done earlier demonstrate similar slow irreversible transfer of energy to the black brane [8] and could be indeed related to weak breaking of global symmetries via non-linearities. 16 If the infrared behavior of the system can be described by hydrodynamics, there is a generic expectation of the existence of an entropy current especially in holographic theories (see [73,74] for a review). A more general understanding of the existence of the entropy current has been led by the construction of an equilibrium partition function [75]. It is therefore of importance to understand if a quasihydrodynamic theory can describe the dynamics of the black brane horizon and thus the infrared behavior of the full system even in the presence of instabilities -see [46] for a hydrodynamic description of the fate of the Gregory-Laflamme instability. 17 In the context of holography, an early attempt has been made in [76,77] by adopting a consistent truncation of the full dynamics to that of the evolution of the energy-momentum tensor operator.

A Computation of the quasinormal modes
We can solve for the quasinormal frequencies of the semi-holographic system numerically by modification of the usual procedure described in [26]. In the usual case, one can reduce the problem of finding the quasinormal mode spectrum to a linear eigenvalue problem. In the semi-holographic case, it will be a cubic eigenvalue problem that smoothly reduces to the usual linear problem in the decoupling limit.
Following [26], let us define the dimensionless radial coordinate u = r/r h so that the boundary is at u = 0 and the horizon is at u = 1 (note from (3.2) that r h = 3/(4πT )). The domain of interest is u ∈ [0, 1]. As evident from the near boundary expansion (3.7), we can define a function g k,ω (u) via g(u) = 64 27 u 3 f (k, ω, u) − φ (0) (k, ω) − i ωφ (0) (k, ω) − k 2 r 2 2 φ (0) (k, ω) . (A.1) We also drop the k, ω subscripts in g k,ω (u) for notational simplicity. Crucially, g(u) is analytic at u = 1 owing to the ingoing boundary condition. The above definition of g(u) ensures a smooth decoupling limit in which β → 0. It is also convenient to define the dimensionless wave number q = k/(πT ) and the dimensionless frequency = ω/(πT ). We can readily obtain the differential equation for g(u) by substituting the above in (3.4).
We numerically approximate g(u) as a linear combination of the first M + 1 Chebyshev polynomials which are linearly mapped onto the domain the unknown expansion coefficients. We extract the coefficient q ij which multiplies the j th expansion coefficient in the equation for the i th grid point. We assemble a matrix Q = ||q ij ||, such that the complete set of equations takes the schematic form Q · coeffs + φ (0) P = 0 (A. 2) where P is a vector.
Since the equation for g(u) involves φ (0) (k, ω), we need to give an additional input, namely the boundary condition at r = 0 given by (3.5) and (3.6) which can be rewritten in the form where we have used φ (3) = (πT ) 3 g(u = 0) which follows from the defining equation (A.1) for g(u) and the asymptotic expansion (3.7), andβ ≡ β √ πT is the dimensionless mutual coupling. 18 Obviously, g(u = 0) is also a linear sum of the Chebyshev coefficients. We can then combine these coefficients and φ (0) (k, ω) into a column vector V , and also (A.2) and (A.3) together into a matrix equation of the form Q · V = 0. (A.4) The above system of linear equations for the elements of V will have solutions only for certain values of (the dimensionless frequency) for which det Q = 0. These values of will constitute the quasinormal spectrum of the full system. However, one can follow a strategy which is better than solving for det Q = 0 for determining the quasinormal mode spectrum.
We note that each element of Q is at most cubic in . 19 Therefore, we can write Q = Q 0 + Q 1 + 2 Q 2 + 3 Q 3 where Q 0 , Q 1 , Q 2 , Q 3 are independent of . Thus (A.4) is simply a cubic eigenvalue problem whose solutions give us the desired QNM for a given value of q, the dimensionless wave number.
We can readily solve a cubic eigenvalue problem by converting it into a generalized eigenvalue problem. To see this, we simply note that we can rewrite (A.4) in the form 18 Note that in d = 3 the boundary scalar field χ has mass dimension 1/2. Also βχ is dimensionless since φ (0) is the source of a marginal CFT operator. 19 The cubic term arises from the dependence of H on the third time-derivative φ (0) as explicit in (3.6).
In absence of semi-holographic coupling, Q depends only linearly on . where O is the null matrix and I is the identity matrix (with the same rank as Q). We can obtain the full QNM spectrum by using standard routines for solving the generalized eigenvalue problem of the type (A.5).

B Iterative procedure for computing the non-linear dynamics
The non-linear dynamics of the full system can be solved by the following iterative procedure with the input initial conditions for χ(t 0 ), χ (t 0 ), a (3) (t 0 ) and the initial profile of the bulk dilaton Φ(r, t 0 ). Note that it is not necessary to specify higher order time derivatives of χ at initial time although the equation for time-evolution of χ given by (4.14) is third order. Nevertheless the initial profile of the bulk dilaton Φ(r, t 0 ) needs to be consistent with the near-boundary radial expansion of χ, i.e. we need Φ(r, t 0 ) = −βχ(t 0 ) − βrχ (t 0 ) + f (r) with f (r) = O(r 2 ) at r = 0 so that it is consistent with (4.10) after setting φ (0) = −βχ.
The dynamics of the full system can be solved by the following iterative procedure. The initial conditions described above should be held fixed at each stage of the iteration.
1. We first solve the boundary scalar field equation (4.14) by setting its right hand side to zero. This implies that at the first step of the iteration χ(t) = χ(t 0 ) + χ (t 0 )t which gives the input to the holographic system by specifying the source for all time.
2. We proceed to solve the gravitational system of equations with the input source φ (0) as determined above, and the initial bulk profile of Φ and the initial value of a (3) . This can be achieved by utilizing the nested structure of the gravitational equations given by Eqs. (4.2)-(4.6) and employing the spectral method discussed in [78,79]. We use 30 Chebyshev grid points to capture the dependence on the holographic radial coordinate and a fourth order Adam-Bashforth time-stepping to evolve in time. We also choose a suitable radial cutoff. Solving the full non-linear system of equations, we can extract φ (3) (t) for from Φ(r, t).
3. Taking φ (3) (t) and χ(t) as inputs from the previous iteration we now compute the right hand side of the boundary scalar field equation (4.14). We solve this linear equation again with the source determined by the previous iteration and with the same initial conditions as before.
-37 -  Figure 17. The total energy is plotted for various iterations with initial conditions given by (4.16) with it 1 standing for first iteration, etc. We get convergence in 4 iterations as evident from the total energy being constant in the fourth iteration to a very good accuracy.