Hybrid Fluid Models from Mutual Effective Metric Couplings

Motivated by a semi-holographic approach to the dynamics of quark-gluon plasma which combines holographic and perturbative descriptions of a strongly coupled infrared and a more weakly coupled ultraviolet sector, we construct a hybrid two-fluid model where interactions between its two sectors are encoded by their effective metric backgrounds, which are determined mutually by their energy-momentum tensors. We derive the most general consistent ultralocal interactions such that the full system has a total conserved energy-momentum tensor in flat Minkowski space and study its consequences in and near thermal equilibrium by working out its phase structure and its hydrodynamic modes.


Introduction
The hydrodynamic analysis of heavy-ion collisions performed at RHIC and the LHC suggests that a droplet of strongly interacting matter is generated in the collisions. The value of the specific viscosity that best describes these data is very low [1,2], η/s 1, suggesting that the plasma is strongly coupled and does not have a description in terms of weakly interacting quasi-particles. This has encouraged much work in describing the plasma formed in terms of strongly coupled models, such as N = 4 super Yang-Mills theory as described by AdS/CFT duality [3] (and references therein).
While the low value of η/s implies that the system is strongly coupled, the collisions exhibit also hallmarks of weak coupling dynamics. In particular, it is seen that the hard components of high-p T jets go largely unmodified and resemble those created in p-p collisions [4]. This suggests that the medium formed in heavy-ion collisions cannot be strongly coupled at all scales and even if some of the modes are strongly coupled, others are weakly coupled. Even more strikingly, the interpretation of the observed long range rapidity correlations in p-A and high multiplicity p-p collisions through final state interactions, combined with no signature of jet quenching in these systems may be seen to suggest the presence of both strongly and weakly coupled modes.
The simultaneous presence of strongly and weakly coupled modes poses a theoretical challenge. In absence of any fully developed non-perturbative method to access real-time properties of QCD in the non-perturbative regime, we may attempt to model the nonperturbative modes using a theory that we can solve in the strong coupling limit, while discussing the perturbative sector in a weak coupling approximation. Corresponding attempts have been made in [5][6][7]. Such approaches in general pose the non-trivial question then, how the two sectors, described with different models describing different degrees of freedoms, be coupled.
A consistent coupling requires that the quantities that mediate the coupling should be well-defined in each theory and also gauge-invariant. In the context of jet-quenching, such couplings have been suggested for example in [5,8].
In a different attempt to formulate a generic coupling between the two subsectors for the study of collective dynamics and equilibration, a local coupling of all the marginal operators of the two subsectors was proposed [6,7], following previous examples of a semi-holographic framework where only part of the dynamics is described by gauge/gravity duality [9][10][11][12]. This includes in particular a coupling between the energy-momentum tensors of the two subsectors, which can be induced by deforming the boundary metric of a holographic sector.
Specifically, as a semi-holographic model of the early stages of heavy-ion collisions, the perturbative sector was assumed to be described by classical Yang-Mills equations as in the glasma effective theory [13] that describes the color-glass condensate initial conditions [14,15] of the deconfined gluonic matter liberated in the heavy-ion collisions, and the nonperturbative infrared sector by AdS/CFT, corresponding to strongly coupled N = 4 super-Yang-Mills theory. The toy model studied in [7] demonstrated that in this way a closed system with a conserved energy-momentum tensor in Minkowski space can be obtained. 1 In this work, we explore the implications of the "democratic" couplings proposed in [17] (and extensions thereof), where the effective metric of each subsystem depends on the energy-momentum fluctuations of the complements. Instead of far-from-equilibrium systems studied previously we concentrate on systems that are in equilibrium and the equilibration of near-equilibrium systems. Furthermore, we again restrict the couplings to that between the respective energy-momentum tensors of the subsystems.
We first address the question of what is the equilibrium state of two coupled conformal systems. As the system is assumed to be in thermal equilibrium, and the coupling depends only on the energy-momentum tensors of the subsystems, the microscopic features of the subsystems do not enter the discussion and therefore the results are generic for conformal subsystems and depend only on the properties of the coupling between the subsystems. We observe that requiring causality and ultraviolet completeness restricts the range the model parameters describing the coupling can take. In addition we find that the composite system -that breaks conformal symmetry due to dimensionful parameters of the coupling -exhibits a rich phase structure with a phase transition that takes the system from a sum of two separate conformal subsystems at low temperatures to a new emergent conformal system at high temperatures. As a function of the model parameters, this transition is either a cross-over or a first-order transition, and the two are separated by second-order critical endpoint with specific heat critical exponent α = 2/3.
Next we study in detail the collective dynamics of near-equilibrium systems. We first assume that each subsystem can be separately described as a conformal fluid in terms of first-order hydrodynamics. This assumption is generally valid if the length scale of the deviation of global equilibrium is sufficiently long for a well behaved gradient expansion and if no long-lived non-hydrodynamic modes are excited. Within this approximation we follow how linearized energy-momentum perturbations of the composite system approach global equilibrium and find a rich structure of two-fluid dynamics. In the shear sector we find that the overall viscosity interpolates between those of the subsystems and decreases with the coupling between the subsystems. In the sound sector we obtain two modes where only one is propagating with the thermodynamic speed of sound at large coupling. However both have attenuation vanishing with the square of momentum, implying that spatially homogeneous density perturbations of the individual subsystems are not attenuated and therefore more dynamics is required for the thermal equilibrium to be established between the two sectors. Indeed, this is in line with the findings in the semi-holographic toy model of Ref. [7], where also interactions beyond the ones between the energy-momentum tensors are needed for thermalization [? ].
Finally, we study to what extent non-hydrodynamic modes in one subsystem are attenuated because of coupling to the other dissipative subsystem.
The organisation of the paper is as follows. In Section 2, we describe the general setup, its motivation in the semi-holographic context, as well as the concrete mutual metric coupling and how a total energy-momentum tensor that is conserved with respect to the (Minkowski) background metric of the full system arises. In Section 3, we discuss the requirements of causality and UV-completeness and study the consequences of our couplings for the thermodynamics and phase structure of the full system. In Section 4, we study the hydrodynamic limit of the full system, and in Section 5 we further study the case when the weakly coupled system can be described by kinetic theory and the strongly coupled sector as a conformal fluid with appropriate transport coefficients.

Semi-holography and democratic coupling
We consider a dynamical system S in a fixed background metric g (B) µν (to be set to the Minkowski metric η µν eventually) which consists of two subsystems S 1 and S 2 .
In the previous approach to semi-holography [7], the full effective action S of system S was constructed as [7]: where S pert is the effective perturbative action for S 1 , and S 2 is represented by S hol , which is the holographic on-shell gravitational action in presence of sources. These sources, such as a non-trivial boundary metric g µν , are functionals of the gauge-invariant operators of the perturbative sector, and γ is a hard-soft coupling. The full conserved energy-momentum tensor calculated by varying the action with respect to g (B) µν cannot be written in terms of the effective operators of each sector and therefore the low energy dynamics of the full system cannot be readily derived from the coarse-grained descriptions of the individual subsystems. Furthermore, the way the two sectors are coupled is somewhat asymmetric. On the one hand, the coupling amounts to deforming the metric of S 2 by the energy-momentum tensor of S 1 . On the other hand, the energy-momentum tensor of S 2 enters via the equations of motion of S 1 . Nevertheless the main improvement of [6] made in [7] was that the full energy-momentum tensor is conserved, provided that the effective operators of S 2 satisfy a separate Ward identity and the equation of motion for the fields in S 1 are in effect.
Motivated by the semi-holographic approach in the democratic formulation [17], the two subsystems are assumed to have covariant dynamics with respect to individual effective metrics g µν andg µν , respectively. 2 Interactions between the two subsystems are introduced by promoting each effective metric to functions that are locally determined by the state of the complement system, The two subsystems are assumed to share the same topological space so that we can use the same coordinates for both of them (and thus the total system) 3 . Furthermore, the subsystems appear as closed systems with respect to their individual effective metrics, but they can exchange energy and momentum defined with respect to the actual physical background metric g (B) µν . Thus the effective metric tensors encode the interactions between the two subsystems.
The diffeomorphism invariance of the respective theories describing the two subsystems imply the Ward identities where ∇ and∇ refer to the covariant derivatives with respect to the different effective metrics with the Levi-Civita connections is the corresponding Levi-Civita connection, and the second equalities in (2.4) indicate that from the point of view of the actual physical background metric g µν the identities (2.3) actually imply that work is done on the respective subsystems by external forces. In what follows, we restrict the forms (2.2) of the effective metrics g andg (in a generally covariant manner) such that there exists a T µν for the full system that is locally conserved with respect to the physical background metric g µν , i.e., we can enforce the Ward identity for the total system: It turns out that the full energy-momentum tensor T µν is a functional only of the effective operators t µν andt µν of the two sectors. Hence one can readily construct effective descriptions of the full dynamics from the effective description of each sector. 4 The main advantage of our method in the context of phenomenology is that it works even when we cannot invoke action principles for the effective descriptions of one or both subsystems. The full dynamics is obtained by solving the subsystems in a mutually self-consistent way as has been illustrated in case of the vacuum state in a toy example [17].
In the present paper, we utilize this to construct the low energy phenomenology by considering appropriate effective description of each subsector. First we assume that both sectors are described by fluids. Then we describe the perturbative sector by an effective kinetic theory and the non-perturbative sector by a strongly coupled fluid. We will be able to find consistent solutions for the full thermal equilibrium and also study its linear perturbations. 3 Coordinate transformations would thus affect the background metric of the complete system and the effective metrics of the subsystems simultaneously. 4 Additionally such couplings can generate expectation values of high-dimensional irrelevant operators without the need of introducing a non-trivial irrelevant deformation of the respective theory [17]. This feature is needed for the cancellation of the Borel poles of the perturbative expansion.
As a general remark, the principle of democratic coupling can be extended to other couplings such as that between scalar operators O andÕ [17]. Let the theory describing the non-perturbative sector be also a (strongly coupled holographic) Yang-Mills theory with the couplingg YM whereas g YM is the coupling of the perturbative sector. These mutual deformations by scalar operators lead to the modified Ward identities (we turn off other couplings including the effective metric couplings for purpose of illustration) Then we may postulate a democratic coupling of the form: where g 0 YM andg 0 YM are constants. It is clear then that the above Ward identities along with (2.7) imply the existence of the conserved energy-momentum tensor T µ ν given by satisfying ∂ µ T µ ν = 0. In [17], the most general scalar couplings of the form have been explored and a toy construction has been done to illustrate how these "hard-soft" couplings (such as α) along with the parameters of the holographic classical gravity determining S hol can be derived as functions of the perturbative couplings in S pert via simple consistency rules. In the following subsection, we extend and correct the democratic effective metric type couplings set up in [17].

Consistent mutual effective metric couplings
We start the construction of the coupling rules between the two subsystems by demanding that the total system a conserved energy-momentum tensor T µν can be written for the total system S in the flat background metric (from now on we choose g (B) µν = η µν unless explicitly mentioned otherwise) while simultaneously satisfying the Ward identities of the two subsystems in their respective curved metrics where ∇ and∇ refer to the covariant derivatives with respect to the different metrics of the subsystems, with the corresponding Christoffel symbols (2.4).
For the rest of the paper, unless explicitly indicated otherwise, by t µ ν we will mean t µρ g ρν and by t µν we will mean g µρ t ρσ g σν , etc., with all lowering (and raising) of indices done by the effective metric (and its inverse) of the respective theory. The Ward identity of subsystem S 1 implies that i.e., and multiplied both sides of (2.12) with √ −g to obtain (2.13). Similarly, the Ward identity for subsystem S 2 implies that We require that both t µν andt µν are symmetric tensors. Using these Ward identities, it is straightforward to verify that the following local relations for the effective metrics where γ and γ are coupling constants (with mass dimension −4), allow us to construct a symmetric conserved tensor for the full system in flat space. From (2.13) and (2.15) it follows that Similarly it is easy to see that A symmetric and conserved total energy-momentum tensor T µν = η µρ T ν ρ = T µ ρ η ρν with ∂ µ T µν = 0 (also ∂ µ T µ ν = 0) can therefore be defined by We can easily generalize the above construction for a curved background metric g (B) µν instead of the Minkowski metric using the second identities in (2.4) which imply where we have used that √ −g/ −g (B) is a scalar under general coordinate transformations. With the help of these relations, one can readily see that the consistent coupling rules have the following general covariant forms , . (2.24) Then with the full conserved energy-momentum tensor is again given by (2.22), and it satisfies ∇ (B) µ T µ ν = 0 in the actual background where all degrees of freedom live. (Note T µν = T µ ρ g (B)ρν ). More general consistent couplings can be constructed if we permit higher powers of the energy-momentum tensors t µν andt µν together with new coupling constants carrying correspondingly higher inverse mass dimension. This is done in Appendix A (correcting and generalizing Ref. [17] in this respect); in the following we will restrict ourselves to the above two coupling terms with coupling constants γ and γ . 5 While the dimensionful coupling constants introduced here appear to be arbitrary at this stage, we shall see that certain restrictions appear when further physical requirements are imposed. In particular, the complete dynamics should be such that causality remains intact. This means that the effective lightcone speed in the subsystems should not exceed the actual speed of light defined by g (B) µν . At least in the following applications to equilibrium and near-equilibrium situations, we can confirm that with just the two terms of in the consistent coupling rules corresponding to γ and γ causality can be ensured -at arbitrary energy scales -by choosing γ > 0 and r ≡ −γ /γ > 1. Interestingly enough, a positive value of the tensorial coupling constant γ was also found to be required in the semi-holographic study of Ref. [7] in order that interactions lead to a positive interaction measure, E −3P = −T µ µ > 0, which is a feature of (lattice) Yang-Mills theories at finite temperature [18,19].

General equilibrium solution
We now assume that the full system S, living in a flat Minkowski background metric g (B) µν = η µν and composed of two sectors S 1 and S 2 that interact through mutually determining 5 Appendix A points out that in the most general set-up, where the total energy-momentum tensor satisfies thermodynamic consistency proven in Appendix B, each possible interaction term in the total energy-momentum tensor can be obtained via an appropriate coupling rule as a result of an interesting combinatoric identity. This is significant because a generic interaction term is not ruled out by any symmetry, and therefore it should indeed be reproduced by our way of introducing interactions via effective metrics. their effective metrics, has reached a homogeneous and isotropic equilibrium state with temperature T with total energy-momentum tensor Assuming furthermore that the subsystems S 1 and S 2 have also thermalized due to their internal dynamics taking place in their respective effective metrics, we expect that the latter will have a static, homogeneous and isotropic form for which we introduce the ansätze with constants a, b,ã,b to be determined self-consistently. 6 The energy-momentum tensors of the subsystems are then of the form i.e., with individual temperatures T 1 and T 2 . The simplest coupling rules (2.16) now read and these determine a, b,ã andb as functions of T 1 , T 2 and the coupling constants γ and 6 If one of the systems is to be described by gauge/gravity duality, the simple metric ansatz above is of course not pertaining to the bulk, but to the boundary of the gravity dual.
γ . The full energy-density and pressure following from (2.22) are In a thermal equilibrium for the full system as well as for its individual subsystems, the physical temperature T of S living in Minkowski space is given by the inverse of β 0 dτ , where τ is imaginary time and β its period. The temperature of the subsystem S 1 , which effectively lives in a metric with constant √ −g 00 = a, is then given by by the same token we have T −1 2 =ãT −1 . Hence, Thus T alone parametrizes the space of equilibrium solutions. Using the thermodynamic identities 1,2 + P 1,2 = T 1,2 s 1,2 , E + P = T S, (3.9) the result (3.6) implies showing that the total entropy density is the sum of the two entropy densities. Therefore, we identify the total entropy current as for s µ 1 = s 1 (T 1 )u µ , s µ 2 = s 2 (T 2 )ũ µ , and S µ = SU µ and U µ = (−1, 0, 0, 0). This indeed makes perfect sense in a general non-equilibrium situation. When each sector has an entropy current s µ 1,2 satisfying such that (3.14) In thermal equilibrium, we also need to have or, equivalently, dP/dT = S, for thermodynamic consistency. In Appendix B we prove this relation and the consistency of (3.8), (3.10) and (3.15) for the coupling discussed here as well as for the coupling rules that generalize (2.16). The mutual compatibility of the thermodynamic identities (3.9) and (3.15) of the full system with the global equilibrium condition (3.8) (along with the additivity of the total entropies that can be expected from the fact that each subsystem is closed in an effective point of view) provides a strong low-energy consistency check of our approach.

Causal structure of equilibrium solution
Since the causal structure of the dynamics taking place in the subsystems is dictated by the respective effective metrics only, causality in the full system, which is living in Minkowski space, is not guaranteed. For example, massless excitations from the point of view of subsystem S 1 with metric g µν = diag(−a 2 , b 2 , b 2 , b 2 ) propagate with velocity v = a/b with respect to the actual physical spacetime that the full system is occupying. (Recall that the two subsystems and the full system share the same topological space; the effective metrics of the subsystems just encode the effects of interactions between the two components of the full system.) At least for the above solution for the equilibrium configuration obtained in the case of the simplest coupling rules (2.16) we can ensure the absence of superluminal propagation by requiring that the tensorial coupling constant γ ≥ 0 together with P 1,2 ≥ 0 and 1,2 ≥ 0. To see this, take the sum of the first and second as well as of the third and fourth equation in (3.5). This leads to independent of γ , implying that the effective lightcones defined by the metrics g µν andg µν are contained within the lightcone defined by the background Minkowski metric.

Conformal subsystems
In the following we shall consider the case of two conformal subsystems. For example one may think of a gas of nearly free massless particles for S 1 coupled to a strongly interacting quantum liquid for S 2 . The energy-momentum tensors t µν andt µν are assumed to be traceless with respect to the effective metrics g µν andg µν , thus the equations of state of the two subsystems are then simply with constant prefactors n 1 and n 2 . The entropy is a simple expression in terms of the effective lightcone velocities v,ṽ associated with the effective metrics g µν andg µν , respectively, where Similarly we obtain for (3.5) is a dimensionless coupling constant that we shall use from now on in exchange for γ . Eliminating b andb yields the two equations Since causality implies 0 < v,ṽ < 1, we see that solutions exist for arbitrarily high T only when the denominator on the right-hand side of (3.22) is able to reach a zero and is positive. This is the case when r > 1, which thus turns out to be a necessary (as well as sufficient) condition for ultraviolet completeness for the simplest coupling rules (2.16); otherwise this model would exist only up to some finite value of T . As shown in Appendix C, the high-temperature behavior of the total system is governed by the fact that the metric factors a,ã, b,b asymptote to linear functions of the physical temperature T . Since the effective temperatures of the subsystems are given by T 1 = T /a and T 2 = T /ã, they stop growing together with T and and instead saturate at finite values proportional to γ −1/4 . For r = 2 figure 1 displays this behavior for equal and unequal subsystems, i.e., n 1 = n 2 and n 1 = n 2 , respectively.
Although the subsystems are conformal, when the two sectors interact, the full system in general is no longer conformally invariant. With the simplest coupling rules (2.16) and the resulting solution (3.6) we find Effective temperatures of the subsystems as a function of the physical temperature with r = 2 for equal and nonequal (n 2 = n 1 /10) subsystems (left and right panel, respectively). As the physical temperature increases, the effective temperature of the subsystems first increases in line with the former (the dotted line marks equality), but when T becomes larger than γ −1/4 , the effective temperatures asymptote to a limiting value. This limiting value is larger for the subsystem with fewer degrees of freedom.
Note that the term in square brackets in (3.24) is the square root of the denominator in (3.23); it is positive in the uncoupled case where v =ṽ = 1, and it cannot change sign for any finite value of γT 4 . Therefore, the conditions for causality γ > 0 and condition for ultraviolet completeness, r > 1, imply that the interaction measure E − 3P = −T µ µ is positive (as is the case with lattice QCD results), and that the full system approaches conformality at large temperature T .
While in general we have to resort to numerical evaluations, one can also derive perturbative expansions for all quantities (see Appendix C). For small couplings or for small temperature, 7 |γ|, |γ | T −4 , the resulting a,ã, v, andṽ are all close to unity, and thus E − 3P ≈ 24γn 1 n 2 T 8 , i.e., the full system approaches conformality at small temperature as expected.
This behavior can also be seen in the speed of sound (squared) of the full system, defined thermodynamically by which expanded up to third order in γT 4 reads (3.26) With conformal subsystems the dependence on r = −γ /γ appears only at third order; in quantities which only depend on v andṽ, as is the case for the entropy, also the third-order term is still independent of r.

Equal subsystems
For the special case n 1 = n 2 where v =ṽ, the numerical solution of (3.22) is displayed in Fig. 2 for various values of r > 1. 8 It turns out that for 1 < r < r c ≈ 1.1145 more than one solution exists. This corresponds to a phase transition that will be discussed in section 3.3.3. Concentrating first on the case r > r c , the behavior of the pressure and the interaction measure (divided by T 4 ) is shown in the left panel of Fig. 3 for a typical case (r = 2). Intriguingly, P/T 4 shows an increase somewhat reminiscent of the deconfinement crossover transition in QCD.
The speed of sound (squared) (3.25) shown in the right panel of Fig. 3 exhibits a pronounced dip, indicating a crossover as opposed to a phase transition as γ 1/4 T is increased from the conformal situation at γ 1/4 T = 0 to large values, where it asymptotes again to conformal value 1/3.
Since S/T 3 ∝ v −3 , the entropy increases from its interaction free value at T = 0, where v = 1, in parallel to the drop in v displayed in Fig. 2 In the case of two identical conformal subsystems, the relation between the effective lightcone velocity v and γT 4 is given by the roots of a polynomial equation of 9th degree (explicitly given in (D.2)), which in general can only be solved numerically. The asymptotic value of v is however determined by a simple quadratic equation which yields Evidently, the entire physical range 0 < v ∞ < 1 is covered as r varies between unity and infinity.
Left panel: Pressure (black line) and trace of the energy-momentum tensor (red) divided by T 4 , with the asymptotic value of the pressure indicated by the short dashed line; right panel: speed of sound squared (full black line) -both for n 1 = n 2 = 1 and r = 2. As γ 1/4 T increases from small to large values, a crossover between regimes with different values of P/T 4 takes place that is accompanied by a dip in the speed of sound which takes on a conformal value in both asymptotic regimes. At large γ 1/4 T and for sufficiently low values of r (including the case r = 2 at hand), the speed of sound in the full system turns out to be larger than the effective lightcone speed v of the subsystems (green dashed line: v 2 ).
Since the speed of sound c s approaches the conformal value 1/ √ 3 at large γ 1/4 T , for sufficiently small values of r (namely r < 7/3), c s can be larger than the effective lightcone speed v of the subsystems. This is no contradiction to causality, since besides dynamics within the subsystems, there is also collective dynamics between them. (In Section 4.2 this will be studied further.)

Unequal subsystems
For unequal systems one can show (using formulae (3.22) and (3.23)) that there exist solutions for v andṽ in the limit γT 4 → ∞ for any value of n 2 /n 1 and r > 1. They are given by the (sextic) equations which have a unique solution in the domain 0 < v ∞ ,ṽ ∞ < ∞ when r > 1. In the extreme limit that one of the systems completely dominates, say n 2 /n 1 → 0, the asymptotic effective lightcone velocity of the smaller system approaches zero,ṽ ∞ ∼ O((n 2 /n 1 ) 1/5 ), while the dominant system has the limit v ∞ → √ 1 − r −1 . In Fig. 4, the full numerical solution of the effective lightcone velocities is displayed for n 2 /n 1 = 1/10 and r = 2 as well as the resulting entropies of the two subsystems. While the smaller subsystem has a much larger relative growth of S/T 3 than the larger subsystem, the latter remains dominant. (Considering again the extreme limit n 2 /n 1 → 0, S 2 /S 1 changes from being of order n 2 /n 1 at low γT 4 to (n 2 /n 1 ) 2/5 at high γT 4 .) At the value r = 2 used in Fig. 4, the behavior of the speed of sound is similar to the case shown in Fig. 3. Again, there is a dip at the crossover between the regimes of  small and large γ 1/4 T , where c s asymptotes to the conformal value 1/ √ 3. In the case displayed in Fig. 4, now only one of the effective lightcone velocities, namelyṽ, falls below the (conformal) value of the speed of sound at large γ 1/4 T .

Phase transition
Perturbative expansions in the dimensionless parameter γT 4 turn out to work rather poorly and indeed have to break down for 1 < r < r c where multiple solutions appear at finite values of γT 4 , as shown in Fig. 2 for n 1 = n 2 . For 1 < r < r c ≈ 1.1145 in the case n 1 = n 2 and 1 < r < r c ≈ 1.25 for n 1 = n 2 , this corresponds to a first-order phase transition that turns into a second-order phase transition at r c .
In Fig. 5 pressure and entropy are plotted in the region around the first-order phase transition with n 1 = n 2 = 1 and r = 1.1. 9 The range in γ 1/4 T where the pressure has three solutions corresponds to the possibility of superheating or supercooling (depending on whether the phase transition is approached from higher or lower temperatures). This happens if one does not immediately switch to the thermodynamically preferred phase with higher pressure (lower free energy). The third solution which directly connects the endpoints of superheating and supercooling is always thermodynamically disfavored and cannot be accessed physically, because it comes with negative specific heat (corresponding to the part of the curve for the entropy with negative slope).
In Fig. 6 the effective temperature of the subsystems is shown for the same set of parameters. This shows a curious nonmonotonic behavior; at the phase transition the effective temperature jumps and approaches the asymptotic value from above as the physical temperature goes to infinity. In fact, although hardly perceptible in the left plot in Fig. 1,   the effective temperature also approaches the asymptotic value from above for r = 2 in the crossover region; only for r 2.048 (in the case of n 1 = n 2 ) the effective temperature eventually shows monotonic behavior.
At r = r c the phase transition becomes second-order with continuous pressure and entropy. In Appendix D the parameters of the second-order phase transition are obtained in closed form for n 1 = n 2 . In particular the critical exponent α that characterizes the behavior of specific heat, is obtained, with the result which is independent of n 2 /n 1 . It is thus different from any mean-field result, and it is also larger than the value in the Ising model (α ≈ 0.11) or in the polymer models (α ≈ 0.236), which are the largest values occurring in N vector models (for N = 1 and N = 0, respectively) [20]. The comparatively large value of α in (3.29) is curiously close to that obtained in the matrix model for deconfinement of Ref. [21], which yields α = 3/5. The qualitative features of the phase transition are the same for unequal conformal subsystems: For 1 < r < r c , the transition is first order, at r = r c the phase transition is second order, and for r > r c it is a crossover. Furthermore, the critical value r c shows a rather weak dependence on n 2 /n 1 , it lies in the narrow interval 1.119 . . . < r c < 1.25, and the critical exponent α at the second-order phase transition point r = r c is always 2/3 (for more details see Appendix D).

Massive subsystems
The simplest coupling rule (2.16) with r > 1 also makes sense for more general equations of state for the subsystems. In Fig. 7 we display the results for the speed of sound (squared) that is obtained by coupling two free Bose gases with various masses (for simplicity with r = 2, where only a crossover and no phase transition arises). When both subsystems have particles with mass, the speed of sound starts from zero at γT 4 = 0, and approaches the conformal value at large γT 4 . (When one or both components contain massless particles, c 2 s also starts from the value 1/3.) The way approximate conformality is approached at high T is again similar to the conformal case discussed above, although we cannot demonstrate this analytically as in Appendix C. The high-temperature behavior (at r > 1) is governed again by an asymptotically linear behavior of the metric coefficients a,ã, b,b ∼ T . Such a behavior is at least consistent with the (simplest) coupling rules (2.16): Once a,ã, b,b have grown sufficiently large, these equations are homogeneous of degree two in the metric coefficients, provided the effective temperatures T 1 , T 2 become constant, which is the case when a,ã ∼ T .
However, an important difference to the conformal case is that the trace-term ∆Kδ µ ν in the full energy-momentum tensor is no longer subdominant, but in fact needed to cancel the contributions to the trace of the full energy-momentum tensor at order T 4 . This is a consequence of the form (3.10) of the full entropy, S = s 1 (T 1 )b 3 + s 2 (T 2 )b 3 ∼ T 3 , together with thermodynamic consistency, S = dP/dT (which is proved in Appendix B for arbitrary equations of state of the subsystems). We expect that it is equally possible to couple more involved equations of state than gases of free massive particles with the simplest coupling rule and to obtain a UV-complete setup. However, we assume that the subsystems have self-interactions so that they are able to thermalize on their own.

Bi-hydrodynamics
In the following we investigate the linearized perturbations of the full hybrid system about thermal equilibrium for given values of the hydrodynamic transport coefficients within the two subsystems, i.e., parameterising their energy-momentum tensors to first order in the gradient expansion according to and similarly for the second subsystem with metricg µν .
Owing to the rotational symmetry of thermal equilibrium, the perturbations can be classified into three distinct sectors, which are called the shear, sound, and tensor channels. Each channel has distinct low energy characteristics. If we take the hydrodynamic limit in both sectors, only the shear and sound channels yield dynamic propagating modes with distinct forms of dispersion relations. The tensor channel in the bi-hydrodynamic limit consists only of a response local in space and time (without a pole) which is convenient for calculating the shear viscosity of the full hybrid system using the Kubo formula. For simplicity we will also analyze the case of conformal subsystems, and therefore we will set ζ 1 = ζ 2 = 0. Note ζ 1,2 do not affect the shear channel in any case.

Bi-hydro shear diffusion
In the shear sector, the velocity fields of both sectors point in the same direction but are orthogonal to the momentum (i.e., the direction of propagation) of a perturbation. Without loss of generality we may assume that the momentum k is in the z-direction and the velocity fields are in the x-direction. The (normalized) velocity fields in both sectors including the infinitesimal linearized perturbations then assume the form: The temperatures in both sectors remain unperturbed from their equilibrium values (in the shear channel). Furthermore, we can consistently assume that the effective metrics are: with the non-vanishing components of δg µν and δg µν being: Note that these metric perturbations preserve the norm of the velocity fields (4.2) at the linearized level. It follows then that the hydrodynamic energy-momentum tensors of the two sectors including the linearized infinitesimal perturbations will assume the forms: with the non-vanishing components of δt µν and δt µν being The hydrodynamic equations (2.13) and (2.15) of the two sectors in the two (selfconsistently perturbed) effective metrics are: These hydrodynamic equations automatically guarantee the conservation of the full energymomentum tensor at the linearized level provided the metric perturbations β,β, γ 13 and γ 13 are solved self-consistently in terms of the physical variables ν andν using the linearized version of the effective metric coupling equations. Once again we will assume that only the couplings γ and γ ≡ −rγ are non-vanishing.
detQ(ω(k), k) = 0, (4.13) and the corresponding eigenvectors involve a momentum dependent combination of ν and ν. It is to be noted that these modes are the intrinsic perturbations of the system which can exist without any extrinsic drive such as a perturbation to the fixed background metric g µν where the full system lives.
Of particular interest are the shear-diffusion modes whose dispersion relations assume the characteristic form: where the index I labels different solutions.
As discussed before, we can solve all equilibrium quantities as functions of T , γ and γ so that we can also express D I as functions of these variables. The perturbative expansions of the shear diffusion constants D I are given by: In the decoupling limit, γ 1/4 T → 0 (with fixed r), these two shear diffusion modes clearly reduce to individual shear diffusion modes of the subsystems 1 and 2; with nonzero coupling they instead involve both subsystems nontrivially. The propagating mode corresponding to the first diffusion constant D a involves velocity amplitudes with 10 10 Note that the combination of ν andν in the propagating mode is k−independent. This is so because each element in the matrix Q in (4.12) is O(k 2 ) at the leading order on-shell, i.e. when ω = −iD a,b k 2 + · · · . and therefore it is indeed localized mostly in the first subsystem when γT 4 is small. Similarly, the other propagating mode has and thus is localized mostly in the second subsystem for small γT 4 . For finite γT 4 , both these modes receive significant contributions from both subsystems (see Fig. 9). Furthermore, the dependence on γ of the perturbative expansions (4.15) start only at third order in the perturbative expansion -so this dependence is weak at small γT 4 . We also note that the perturbation expansion in γT 4 evidently breaks down when |κ 1 − κ 2 | γT 4 , irrespective of the values of n 1 and n 2 .
In the coincidence limit of κ 1 = κ 2 = κ, we instead obtain the following perturbative series where one of the diffusion modes turns out to be independent of γT 4 . The propagating mode corresponding to this γT 4 -independent diffusion constant has When n 1 = n 2 , i.e. when the two subsystems are identical, then the propagating mode is exactly given byν = ν (parallel and equal motion within the subsystems). In any case, this mode gets significant contributions from both subsystems even in the decoupling limit γ, γ → 0. The other propagating mode corresponding to the second diffusion constant D b in (4.18) is the following combination of ν andν where When n 1 = n 2 , this mode is exactly given by ν = −ν (anti-parallel and equal motion within the subsystems). This mode evidently gets significant contributions from both subsystems even in the decoupling limit γ 1/4 T → 0 (as long as |κ 1 −κ 2 | γ 1/4 T ). The nonperturbative dependence ofν/ν on γ 1/4 T is displayed in Fig. 9.
In order to define the shear viscosity of the full system, we can consider the tensor channel. Consider an extrinsic homogeneous perturbation such that the background metric in which the full system lives is perturbed by h µν (t) whose only non-vanishing component is h 13 (t). How does the full system respond? The coupling equations imply that the response involves homogeneous perturbations of the effective metrics γ 13 (t) andγ 13 (t). Furthermore, the hydrodynamic equations of the individual systems in the individual effective metrics (cf. (4.7)) imply that the velocity fields ν andν also vanish for homogeneous γ 13 and γ 13 up to second order in the derivative expansion along with the perturbations of the  Figure 9. The relation between the velocity amplitudes ν andν of the shear eigenmodes displayed in Fig. 8 in the form ξ ≡ 2 π arctan(ν/ν). A value of ξ = 0 or ξ = ±1 (with these two latter values to be identified) means that the mode is carried only by subsystem 1 or 2, respectively; ξ = 0.5 or ξ = −0.5 corresponds to exactly equal amplitudes with equal or opposite phase.
temperatures in each sector. Therefore, the linearized perturbations of the hydrodynamic energy-momentum tensors in the individual sectors assume the forms while the coupling rule (2.16) implies that the effective metric perturbations is determined by the extrinsic perturbation h 13 (t) according to the coupled linear equations Solving the above, we obtain (4.23) One can then readily compute the energy-momentum tensor of the full system including the linearized perturbation. We find that it assumes the standard hydrodynamic form with vanishing velocity and temperature perturbations. 11 Explicitly, where P is indeed the equilibrium pressure of the full system given by (3.6). It also follows that the shear viscosity η of the full system is given by It is to be noted that the bulk viscosity plays no role in the shear sector or in the response to a homogeneous h 13 (t) perturbation of the background metric. So even if the individual sectors have bulk viscosities all our results above remain valid. Given that S of the full system is given by (3.10) we readily obtain η/S. We may thus define the Kubo diffusion constant: (4.26) 11 Note it is a priori not obvious that even if the individual sector energy-momentum tensors are hydrodynamic, the full energy-momentum tensor also assumes a hydrodynamic form. This is particularly so because there are two independent entropy currents. Although in this specific example, the full energy-momentum tensor does indeed assume a hydrodynamic form, in the following subsection we will find counterexamples.
In Fig. 8 the shear diffusion constants D I corresponding to shear eigenmodes are compared with the overall diffusion constant D corresponding to the total shear viscosity η/S ≡ T D for various parameters. The overall result obtained from the Kubo formula is seen to be always in between the two D I 's. The left panel shows the situation for two strongly coupled systems with η i /s i = 1/4π, the right panel the one for a more weakly coupled system S 1 . The dashed lines in both panels corresponds to the case that system S 1 contributes dominantly to the pressure (n 1 > n 2 ). In this case the overall viscosity is closer to the viscosity of the dominant subsystem.
All results for the shear diffusion constants or specific viscosities decrease when the effective coupling γ 1/4 T is increased from zero. In the case of the full viscosity there is a slight nonmonotonic behavior in the crossover region between weak and strong coupling between the subsystems. At large coupling all results appear to saturate at finite values.
Solving (4.13) one in fact obtains two additional eigenmodes which are spurious. First, these are non-hydrodynamic, meaning that ω is finite as k vanishes. Second, when k vanishes, these eigenmodes correspond to spontaneous fluctuation of the effective metric components γ 13 andγ 13 without involving any fluctuation of ν,ν or any external background metric fluctuation (i.e., perturbation of g (B) µν ). Such a freaky fluctuation is possible because the perturbed energy-momentum tensors of each sector involves time-derivatives of the effective background metrics. Therefore, these make the coupling equations (4.10) dynamical in the sense that these are differential equations for γ 13 andγ 13 . 12 The spurious modes correspond to this spurious dynamics. The spurious modes are also badly behaved and are acausal (having positive imaginary parts in the dispersion relation) and this is related to the acausal behavior of first-order hydrodynamics. If we embed the hydrodynamics of each sector in kinetic theory/Israel-Stewart framework/holographic gravity, then these spurious modes disappear and are replaced by well-behaved relaxation modes. This will be one of the topics of the next section.
To summarize our findings for shear diffusion and specific viscosity: 1. The full system has two shear diffusion modes with diffusion constants D a,b such that T D a,b decrease monotonically with increasing temperature T before saturating at finite values at large T .
2. The overall specific viscosity η/S derived from the total conserved energy-momentum tensor is in between the values of T D a,b with slight nonmonotonic behavior at the phase transition.
3. When one of the systems has a dominant contribution to the total energy/pressure and a different specific viscosity, the overall specific shear viscosity is closer to that of the dominant system. 12 One may see this from (4.6) -the first-order corrections omitted in these equations involves the time derivative of γ13 andγ13. Therefore even when h13 is set to zero there exist solutions for γ13 andγ13! In this case, the two systems can just have fluctuating effective metrics without any extrinsic perturbation or change in the internal physical variables ν,ν, δT1 and δT2.

Bi-hydro sounds and their attenuations
Owing to the rotational symmetry of the thermal equilibrium state, we can consistently assume that the velocity fluctuations in both sectors are longitudinal, i.e., pointing in the same direction as the momentum k. This longitudinal alignment of the linearized velocity field defines the sound sector. Without loss of generality, we can take ν,ν and k to be in the z-direction. The consistent forms of the effective metrics are: with the non-vanishing components of δg µν and δg µν being: The four-velocity fields in the two sectors including the fluctuations thus are: We may also anticipate that the temperatures also fluctuate from their equilibrium values so that we also have δT 1 e i(kz−ωt) and δT 2 e i(kz−ωt) . (4.30) The non-vanishing components of the linearized perturbations of the individual hydrodynamic energy-momentum tensors then turn out to be: and similarly The linearized coupling equations take the form: δg µν = γ η µρ δt ρσ η σν −g + 1 2 η µρt (eq)ρσ η σν −gg αβ δg αβ +γ η ρσ δt ρσ η µν −g + 1 2 η ρσt (eq)ρσ η µν −gg αβ δg αβ , These should be utilized to eliminate δa, δã, δb, δb, χ,χ, β andβ in favour of the physical dynamical hydrodynamic variables δT 1 , δT 2 , ν andν.
Assuming that the bulk viscosities of each individual system vanishes, the hydrodynamic equations of motion in the respective effective metrics take the form: In order to find the eigenmodes, one can first solve for the effective metric fluctuations δa, δb, δã, δb, χ andχ in terms of ν,ν, δT 1 and δT 2 using the linear algebraic equations (4.33). Substituting these forms above for the effective metric fluctuations, we obtain the four dynamical equations for the four variables ν,ν, δT 1 and δT 2 which yield a determinant. Finally, the dispersion relations of the eigenmodes are obtained by requiring that this determinant vanishes as in case of the shear sector. Before considering the eigenmodes in detail, it is useful to examine the simple case of two identical perfect fluids, i.e., the case of n 1 = n 2 and η 1 = η 2 = 0 (or rather we consider only the leading order in the derivative expansion). We want to prove that in this case one of the eigenmodes propagate exactly with the speed of sound of the full system provided the thermal equilibrium solution also yields identical effective metrics, i.e., a =ã and b =b. This result is valid even if the individual subsystems are not conformal.
In the case of identical perfect fluid systems, we can also assume ν =ν and δT 1 = δT 2 , and furthermore δa = δã, δb = δb, χ =χ = 0 and β =β so that the individual energymomentum tensors and effective metrics are identical. Then this eigenmode can be obtained from ik a ν − iω δs 1 s 1 We also note that the full thermal equilibrium solution is parametrized by the temperature T . Therefore, a variation of T which preserves its relationship with the (identical) individual system temperatures given by (3.8) will produce a solution corresponding to an infinitesimal change of the full system equilibrium temperature T . Thus we can obtain a solution with δa = δã = (da(T )/dT )δT , δb = δb = (db(T )/dT )δT and δT 1 = δT 2 = (dT 1 (T )/dT )δT with δT = T 1 (T )δa + a(T )δT 1 (4.36) being satisfied. Furthermore, we can boost the full energy-momentum tensor. If the full system is boosted by an infinitesimal velocity υ in the z-direction (in background flat space), then the non-vanishing components of its energy-momentum tensor with an overall infinitesimal temperature fluctuation takes the linearized perfect fluid form: Of course if we make δT space-time dependent we also need a spacetime dependent boost υ in order that we can ensure energy-momentum conservation. The conservation of the full energy-momentum tensor in flat space yields the linearized Euler equations: It is guaranteed that the diagonal components of the fluctuations can always be mapped to a change in δT even if the systems are not identical. If we solve β andβ in terms of ν andν using the off-diagonal 03-component of the coupling equations, and then compute the off-diagonal 03-component of the full energy-momentum tensor, we can also define the υ of the full system as an appropriate linear combination of ν andν demanding the form (4.37) of the full energy-momentum tensor. This can always be done. In case of identical systems with identical energy-momentum tensors living in identical effective metrics, the procedure is simpler: eliminate β in favour of ν from the coupling equation and obtain υ in terms of ν from the computed form of the full energy-momentum tensor.
To proceed further, we thus focus on the off-diagonal component T 03 . Specifically, we observe from (4.37) that δT 03 = δT 0 3 = (E + P)υ, δT 3 0 = −(E + P)υ. With our assumptions for the effective metric and the perfect fluid forms of the energymomentum tensor, we should get Furthermore thermodynamic identities for any consistent coupling ensure that E + P = 2ab 3 ( 1 +P 1 ). Therefore it follows from (4.39), (4.40) and (4.41) that any consistent coupling equation should imply b 2 a ν + 1 a 2 β = νa = υ. (4.42) The coupling equations always ensure that conservation of the individual energy-momentum tensor in the individual effective metric leads to conservation of the full energymomentum tensor in flat space. To show that the eigenmode of the full system corresponds to the thermodynamic sound of the full system, we need to turn this argument around. We need to show that the Euler equations of the full energy-momentum tensor in flat space will lead to satisfying the individual Euler equations in individual effective metrics. Clearly, we will generically need identical systems with identical energy-momentum tensors living in identical effective metrics. Otherwise the number of conservation equations of the full system are outnumbered by the individual conservation equations. At the linearized level, we need to show that (4.38) implies (4.35).
We note that thermodynamic variation ensures that δS/S = 2δs 1 /s 1 + 3δb/b since S = 2s 1 b 3 in the case of identical systems. Similarly, δT /T = δT 1 /T 1 +δa/a since T = T 1 a. It is then easy to see that (4.38) implies (4.35) because of the two relations in (4.42) which follows from consistent coupling equations. We then conclude that for any consistent coupling between two identical systems with identical effective metric solutions at equilibrium, the thermodynamic sound will correspond to one of the eigenmodes at the leading order in the derivative expansion. In this mode, the velocity fields in the two identical systems are parallel to each other so thatν = ν.
Even for identical perfect fluid systems there is another eigenmode where δT 1 = δT 2 and ν =ν. In this mode, the velocity fields are anti-parallel to each other so thatν = −ν. Most importantly the thermodynamic relation δT = δ(T 1 a) = δ(T 2ã ) is not satisfied by the fluctuations. This mode does not travel at the speed of thermodynamic sound. When n 1 = n 2 , it turns out that neither of the two eigenmodes does; in this case the thermodynamically defined speed of sound is in between the velocities of the eigenmodes.
When the two systems are identical, and we consider the eigenmode which at leading order propagates at the speed of full system thermodynamic sound, we cannot map the first-order (identical) hydrodynamic fluctuations of the individual systems to that of a hydrodynamic form for the full system. To see this, we may repeat the steps of the above argument with χ =χ = 0 and η 1 = η 2 = 0 and find that for generic η 1 the modified form of (4.42) does not imply that we can obtain (4.35) with first-order corrections from the first-order correction of (4.38) (linearized Navier-Stokes equation in flat space).
Although in general different from the thermodynamically defined sound, the dispersion relations of the eigenmodes have the same characteristic sound-like form,  Figure 10. Sound modes and their attenuation coefficients for equal and unequal conformal systems, same κ = 1 (corresponding to η i /s i = 1/4π), with the slower mode a plotted in blue, and the faster mode b plotted in orange. The black line represents the thermodynamic speed of sound and associated attenuation coefficient from the Kubo formula. The green dashed line shows the light-cone velocities squared of the two subsystems (in the case of n 2 = 1/10 onlyṽ 2 is in plot region). In the case n 1 = n 2 the lines for c 2 a,b meet and could be continued smoothly by switching the designation; however for any n 1 = n 2 we have c b > c a at nonzero γ 1/4 T . The discontinuous behavior of the damping rates Γ a,b for n 1 = n 2 is in fact the limit of smooth curves as n 1 → n 2 from different starting values.
respective attenuation coefficients are given by attenuation coefficient of the faster mode which then coincides with the thermodynamically defined speed of sound moreover becomes independent of γ 1/4 T .
Mode a has velocity and temperature fluctuation fields with perturbative expansions ν = n 1 n 2 1 + 21 2 (n 2 − n 1 )γT 4 + O(γ 2 , γ 2 , k) ν, Above, the + sign refers to the case when the mode is propagating parallel to the momentum k and − sign refers to the case of opposite propagation. Mode b similarly is one in which For equal partial pressures, n 1 = n 2 , mode a and b haveν = ν andν = −ν, respectively, to all orders. It is instructive to compare the attenuation coefficients of these propagating modes with what would have been the hydrodynamic sound attenuation if one of these modes could have been interpreted as a sound channel hydrodynamic mode of the full system. First, assuming that each individual sector is conformal and hydrodynamic as we have done above, one can see from the expression of the conserved energy-momentum tensor of the full system that the trace of the full energy-momentum tensor does not contain any spatial or temporal derivative. This implies that the full system has vanishing bulk viscosity although it is not conformal but has a nonzero trace of the total energy-momentum tensor. Left panel: the relation between the velocity amplitudes ν andν of the sound eigenmodes displayed in Fig. 8 in the form ξ ≡ 2 π arctan(ν/ν); right panel: the corresponding fluctuation amplitude of the total entropy density, δS ≡ δS µ=0 , divided by νT 3 . Mode a and b are given in blue and orange, respectively, with full and dashed lines representing n 1 = n 2 = 1 and n 1 = 1, n 2 = 1/10. The divergence of δS/(T 3 ν) of mode a at one value of γ 1/4 T is due to a zero of ν (corresponding to |ξ a | = 1); here a velocity field is present only in subsystem 2.
Second, with the bulk viscosity vanishing, the sound dispersion relation for a hydrodynamic system in flat space is given by with c s being the speed of thermodynamic sound and the attenuation coefficient being Γ s = (2/3)(η/T S). With η given by (4.25), the latter would be the attenuation coefficient of one of the propagating modes if it could be interpreted as hydrodynamic motion in flat space. We find that none of the propagating modes attenuates in the hydrodynamic way even when one of them travels at the speed of thermodynamic sound as in the case of identical subsystems. In Fig. 10 and 11 our nonperturbative results for the speeds and attenuations of the propagating modes of two strongly coupled fluids, and weakly plus strongly coupled fluids respectively have been plotted. Comparison has been made also with the hydrodynamic sound attenuation Γ s of the full system. For equal partial pressures, n 1 = n 2 , the results for c a,b coincide at γ 1/4 T = 0 and one finite value of γ 1/4 T . The crossing at the latter point is lifted for all n 1 = n 2 such that mode b is always faster than mode a for γ 1/4 T > 0. In the limit n 1 → n 2 , the results for c a,b develop a cusp, while Γ a,b even become discontinuous. For exactly n 1 = n 2 , max(Γ a , Γ b ) is even a constant independent of γ 1/4 T and the results for the two modes could all be connected smoothly. We have however kept the mode labels corresponding to taking the limit n 1 → n 2 . Fig. 12 displaysṽ/v for the corresponding sound eigenmodes as well as the associated (adiabatic) fluctuations of the total entropy density S µ=0 . Again, the seemingly spurious discontinuities at n 1 = n 2 indeed arise when taking the limit n 1 → n 2 starting from n 1 = n 2 .
In summary, our findings for the sound sector are: 1. The thermodynamic speed of sound of the full system c s is always between the velocities of the two sound modes c a and c b , c a ≤ c s ≤ c b , and coincides exactly with one of the modes for n 1 = n 2 .
2. At temperatures above the crossover from weak to strong inter-system coupling, the velocity of the faster mode (b) quickly approaches the thermodynamically defined speed of sound.
3. Near the crossover temperature, the velocity fields of the two modes change their phase with v (orṽ) vanishing at a certain value of γ 1/4 T for mode a (or b). Mode a has out-of-phase oscillations for large γ 1/4 T with decreasing total entropy fluctuations δS µ=0 and speed slower than the thermodynamic speed of sound.
4. At high temperatures, c s and c b approach 1/ √ 3 due to emergent conformality.
5. While c s and c b can become larger than the effective lightcone speeds v,ṽ, the velocity c a of the slower mode remains smaller than v,ṽ; this mode thus lies within both effective lightcones.
6. The value of the attenuation coefficient obtained from the Kubo formula is between that of the sound modes for large γ 1/4 T .
7. While the dependence of the attenuation coefficients on γ 1/4 T is in general complicated, at temperatures sufficiently above the crossover region the slower "non-acoustic" sound mode is always the more weakly damped one.
8. The coupling studied in our setup provides no pure damping modes, i.e. the imaginary part of the speed of sound vanishes as k → 0. This reflects the fact, that this interaction is not sufficient to equilibrate the two subsystems, e.g. in a homogeneous configuration with subsystems at unequal temperatures.

Coupling a kinetic sector to a strongly coupled fluid
In order to obtain a qualitative understanding of a coupled system of weakly interacting and strongly interacting degrees of freedom, we can study the consequences of mutual effective metric coupling of a gas of massless particles (gluons) described by kinetic theory (as subsystem S 1 ) and a strongly interacting holographic gauge theory described by dual gravitational perturbations of a black hole (S 2 ). Using the fluid/gravity correspondence, we may further simplify gravitational dynamics to that of a fluid with a low value of η/s if we are interested in the long time dynamics 13 . Due to the appearance of spurious modes and associated acausalities we need to embed first-order hydrodynamics in a more complete description. Therefore, we embed the strongly coupled fluid in an Israel-Stewart framework 13 Here we are tacitly assuming that the strongly coupled sector thermalizes faster despite its coupling to the weakly coupled gluons. Our results here indicate that the correction to the relaxation dynamics of each sector can be very mild even though the hydrodynamic behavior is modified significantly by the democratic effective metric coupling.
with an extremely small relaxation time. In the future we plan to do a more complete calculation by involving the relaxation dynamics of the strongly coupled sector as described holographically via quasi-normal mode perturbations of a black brane. Following Refs. [22,23] and ignoring for simplicity effects of quantum statistics, the thermal equilibrium of the weakly coupled (and dilute) kinetic sector is described by a Maxwell-Jüttner distribution f 0 (p i ) = n 0 e p µ uµ/T 1 , where p 0 is determined using the mass-shell condition, p µ p ν g µν = 0 for massless gluons when thermal corrections to the mass are neglected. Expressed in terms of p ≡ p x2 + p y2 + p z 2 , we have p 0 = pb/a, u µ = (−a, 0, 0, 0), so that the time-dilation factor a (but not the spatial dilation factor b) drops out, giving The normalization constant n 0 can be fixed as follows. The energy-momentum tensor corresponding to a quasi-particle distribution f is [24] with p 0 = g 0µ p µ satisfying the mass-shell condition. The equilibrium energy-momentum tensor then takes our previously assumed form: where n 1 is our previously introduced (theory-dependent) parameter if n 0 = n 1 π 2 . (5.5) We will therefore set n 0 to n 1 π 2 so that we can directly use our previously obtained results to describe the equilibrium of the full system. For convenience, we use spherical coordinates for the components of the momenta so that p x = p sin θ cos φ, p y = p sin θ sin φ and p z = p cos θ. A linearized fluctuation of the quasi-particle distribution about equilibrium can be written as: For computational purposes, it is useful to split the linear term δf into two parts, each having a specific momentum k and a specific frequency ω component, according to The term δf (eq) can be defined uniquely such that it produces a perturbation δt µν (eq) in the energy-momentum tensor that is of a perfect fluid form. Only the term ∆f will then contribute to the dissipation of energy and momentum. If δg is the (self-consistent) effective metric fluctuation in the kinetic theory, then in the relaxation time approximation δf obeys the linearized Anderson-Wittig equation: with δΓ µ αβ being the linearized Levi-Civita connection obtained from δg and with p 0 also receiving a corresponding linear contribution so that the mass-shell condition is satisfied. Furthermore, in a conformal theory τ should be proportional to T −1 1 and we may parametrize where κ 1 is a constant which will be eventually identified with 4πη 1 /s 1 as before.
The relaxation time in the Israel-Stewart theory in which we are embedding the strongly coupled fluid is similarly set toτ In order to isolate the strongly coupled fluid from the relaxation dynamics, we will take λ very small so thatτ is small. Unlike the kinetic sector where τ determines the shear viscosity (this can be seen via consistent reduction to hydrodynamics), note thatτ of the Israel-Stewart theory is an independent parameter which does not affect the shear viscosity but only second-order hydrodynamics.

Branch cut in response functions of the kinetic sector
We can show that an infinite number of quasi-particle distribution fluctuations decouple from the strongly coupled sector in the sense that all perturbed observables will get contributions purely from the kinetic sector. For instance, fluctuations of the form δf = F (p)G(θ, φ)e −iωt+ik·x , with G(θ, φ) = H 1 (θ) cos(nφ) + H 2 (θ) sin(nφ) and n ≥ 3 (5.11) have vanishing fluctuations of the energy-momentum tensor If also all perturbations in the strongly coupled sector are set to zero, we can then selfconsistently assume that δg µν = δg µν = 0.
where n i = p i /p and τ is the relaxation time in the kinetic sector. Choosing without loss of generality k along the z-direction we obtain Figure 13. Analytic structure of the response function in the kinetic sector. The branch cut arising from (5.15) is given by the thick black line. The pole corresponding to the pure damping mode (5.34), which lies on the second Riemann sheet, is indicated by the cross in violet.
where we have used that T 1 a = T at equilibrium with T being the physical temperature of the full system. The above produces a cut in the response function that stretches in the lower half of the complex ω plane horizontally from −(a/b)k−i/τ (T ) to (a/b)k−i/τ (T ), see Fig. 13. Physically, the factor of a/b (the effective equilibrium lightcone velocity) reflects that the massless gluons propagate along this effective lightcone. The imaginary part turns out to receive no correction when expressed in terms of the full system temperature T .

Poles in response functions of the kinetic sector
We now consider quasi-particle distribution fluctuations which cost energy and momentum. As before, we can split the propagating modes of the full theory into shear, sound, and tensor channels. We focus on the shear and sound channels for exactly the same reason as before -the tensor channel has no hydrodynamic mode and to characterize it properly we require to embed the strongly coupled fluid into gravity which we have not done yet. As expected, we find that some of the propagating modes in both shear and sound channels are identical to the case of the conformal bi-hydrodynamic individual systems described before. This is simply because both the kinetic and Israel-Stewart sectors can be consistently truncated to conformal hydrodynamics individually. In particular, we will see that with the parametrization (5.9) of the kinetic relaxation time, we get exactly the same results as before with κ 1 identified with 4πη 1 /s 1 . This reproduction of bi-hydrodynamics provides a consistency check of our calculations.
In addition to the bi-hydrodynamic modes, there are two other non-hydrodynamic propagating modes in the full system in each shear and sound channel. These contribute poles in the response function. We find that one of these is continuously connected to the damping in the kinetic sector as we switch off the effective metric coupling. We will focus on this particularly because in case of the hydrodynamic sector, Israel-Stewart dynamics has been used simply as a tool for consistent embedding hydrodynamics and not for capturing actual relaxation dynamics. We will see that if the Israel-Stewart relaxation time is set to zero by taking λ → 0 limit, the other damping mode has a smooth limit that captures the effective metric interactions of the kinetic sector with a strongly coupled fluid. Furthermore, if we take k → 0 limit, there is no way to distinguish the shear and sound channels owing to rotational symmetry of the equilibrium. The damping coefficient of the full system will then be the same in both shear and sound channels. This also provides a consistency check of our calculations.
Let us first focus on the shear channel. The effective metric fluctuations of the two sectors will be given by (4.3) and (4.4) as before. In the kinetic sector, the local mass-shell condition p µ g µν p ν = 0 will imply that at the linearized level: δp 0 (p, θ, φ, z, t) = p β a 2 + γ 13 ab cos θ sin θ cos φ e i(kz−ωt) .
Assuming a self-consistent fluctuation of u µ = (1/a, νe i(kz−ωt) , 0, 0) as before, we also obtain: There is no fluctuation in the temperature in the shear channel. The linearized fluctuation of the quasi-particle distribution function takes the form (5.7) with k in the z-direction and Above δf (eq) arises from the fluctuation in p µ u µ since the local equilibrium distribution by definition takes the form n 1 π 2 e −pµu µ /T and T 1 does not fluctuate in the shear channel. Note that δf (eq) indeed reproduces the fluctuation in the energy-momentum tensor which takes a perfect fluid form. The linearized Anderson-Witting equation (5.8) can then be explicitly solved to obtain: In the kinetic sector, the energy-momentum tensor (5.3) after taking into account both effective metric and quasiparticle distribution fluctuations assume the linearized form with the non-vanishing components of δt µν being Comparing (5.21) with (4.6) we see that the perfect fluid parts of the energy-momentum tensor perturbation that originates in the former case from the linearized perturbation in p 0 , p 0 and δf (eq) match perfectly. The dissipative contribution in (5.21) however originates from ∆f and is given by π 13 . Using the solution (5.19) for ∆f in (5.22) we find that: In order to obtain the hydrodynamic limit we need to expand the right hand side above in τ which yields It is easy to note that the expansion in τ is essentially the derivative expansion. Substituting the above form of π 13 in (5.21) and comparing again with the hydrodynamic form (4.6), we find a perfect match with i.e., η 1 /s 1 = T 1 τ /5 and crucially κ 1 = 4πη 1 /s 1 as we have claimed.
The energy-momentum conservation equation with δt µν given by (5.21) and the metric perturbation given by (4.4) amounts to: One can check that the above reduces to the standard hydrodynamic equation (4.7) when π 13 is approximated by (5.24). We can regard (5.23) and (5.26) as the dynamical equations for π 13 and ν. It is to be noted that one can explicitly check that the conservation equation (5.26) is equivalent to the linearized version of the matching condition u µ (t µν − t µν (eq) ) = 0 which says that the projected energy-momentum tensor obtained from the full quasi-particle distribution f should agree with that obtained from f (eq) . In fact, this matching condition is necessary to ensure energy-momentum conservation. At the level of linearized shear-sector fluctuation, the matching condition reduces to Explicitly, we can check that if we use ∆t 01 = 0 with the on-shell form of ∆f given by (5.19) and the equation of motion (5.23) for π 13 to solve for the variables ν and π 13 , we find that indeed (5.26) is satisfied leading to energy and momentum conservation. Embedding the holographic conformal fluid (with 2 = 3P 2 = 3n 2 T 4 2 as before) in the Israel-Stewart framework we obtain: The linearized Israel-Stewart equation of motion ofπ 13 is: The conservation of energy-momentum tensor mirrors (5.26) of the kinetic sector and takes the form: The equations (5.29) and (5.30) are the equations of motion forπ andν. Note that once again the hydrodynamic limit is reproduced by Taylor expansion inτ aboutτ = 0.
There are four propagating modes for each k as discussed earlier. Two of these are exactly the bi-hydro shear-like eigenmodes obtained earlier with diffusion constants D a and D b . We thus reproduce our previous results.
There are additionally two relaxation eigenmodes. One of these eigenmodes is related to the Israel-Stewart relaxational mode and its damping constant becomes large for small λ and therefore can be decoupled. The corresponding propagating mode in this limit is localized mostly in the Israel-Stewart sector and involves the following combination of π 13 andπ 13 where when γT 4 is small.
The damping constant of the other relaxational mode remains finite. It is of the form with perturbative expansion in the limit λ → 0 according to This is interesting because the Anderson-Witting kinetic theory does not have on its own any non-hydrodynamic pole -the mutual metric coupling evidently causes a pole to be generated from the cut discussed above (for γT 4 → 0 it coincides with the cut). This pole is farther from the real axis than the cut when κ 1 > κ 2 , i.e., when the kinetic sector is more weakly coupled than the second sector described by pure hydrodynamics. The corresponding propagating mode involves the following combination of π 13 andπ 13 wherẽ so that it is mostly localised in the kinetic sector as expected in the limit of small γT 4 . Interestingly, when 5κ 1 = 4κ 2 all corrections to Γ 0 /T vanish so that it is exactly 4π/5κ 1 as the perturbation series (5.35) indicates. However, in this case the λ → 0 limit becomes sick because the other mode becomes unstable. This is consistent with the expectation that the non-kinetic sector should have a lower η/s as it is more strongly coupled.
Furthermore, the departure of Γ 0 /T from its decoupling limit value 4π/5κ 1 in the full calculation are found to be very small for any value of γ 1/4 T (see the left panel of Fig. 14). The damping constant Γ I of the Israel-Stewart relaxational mode is evaluated in the right panel which is indeed large for all γ 1/4 T for the small value of λ chosen. However, it turns out that one cannot take the limit λ → 0 for large γ 1/4 T , for Γ I diverges at a certain value of γ 1/4 T beyond which it turns negative, corresponding to an instability. One thus has to keep λ finite in order to decouple this mode.
Repeating the same calculation in the sound channel, we find that we indeed reproduce the bi-hydro sound sector modes and the same damping coefficient Γ 0 .
A remarkable outcome from our calculations is that non-hydrodynamic observables turn out to receive mild or no non-perturbative corrections even when the hydrodynamic sector receives large qualitative and quantitative modifications.

Conclusions
In the semi-holographic approach to the dynamics of quark-gluon plasma with its coexistence of strongly and weakly interacting sectors it has been proposed to introduce a coupling of the respective marginal operators [6,7,17]. In four dimensions, this always includes the energy-momentum tensors which can be coupled to the (effective) metric of the complement subsystem. In this paper we have determined the most general ultralocal mutual effective metric coupling which leads to a total energy-momentum tensor with respect to the flat Minkowski space of the complete system. The effective metric tensors of the subsystems encode the interactions between them; in particular they lead to state-dependent effective lightcone velocities within a subsystem that can be smaller than unity, similar to thermal masses (but different in that the latter reduce the velocity of massless particles depending on their energy).
We have then studied the consequences of mutual effective metric couplings in equilibrium and in a hydrodynamic limit of near-equilibrium situations. Assuming full thermal equilibrium, we have found an interesting phase structure, which can be separated by a first or second-order phase transition, or an analytic crossover, depending on the coupling parameters. With only the two coupling constants of inverse dimension four turned on, we obtained two distinct phases where the one at higher mutual coupling (or equivalently higher system temperature) has a larger number of degrees of freedom per volume and eventually approaches conformality, which is curiously reminiscent of the deconfinement transition in QCD.
Studying the hydrodynamic behavior of such a two-fluid system, we found two modes in both the sound channel and in the shear channel (with our detailed findings summarized already at the end of the respective subsections above). In the shear channel we found a decrease of shear diffusion constants as the mutual coupling is increased, in accordance with the fact that shear diffusion is in general weaker at strong coupling. 14 The overall shear viscosity, which is determined by the Kubo formula involving the total energy momentum tensor, shows a similar behavior, numerically intermediate between the viscosity values of the subsystems.
In the sound sector, we also found two modes, which correspond essentially to in-phase and out-of-phase density perturbations of the two subsystems. One mode always has a velocity close (or equal) to the thermodynamically defined speed of sound, while the other is slower and always below the effective lightcone velocity of the subsystems, and also more weakly attenuated, suggestive of a quasiparticle nature. In the transition region, a role reversal takes place, reminiscent of a similar phenomenon in other two-fluid systems [25]. The damping of the two modes depends in a complicated manner on the shear viscosities in the individual subsystems and their mutual interaction (or system temperature). In the long-wavelength limit, however, both attenuations vanish quadratically with the modulus of the wave vector, which means that completely homogeneous and isotropic density perturbations do not equilibrate. This is similar to the behavior found in the semi-holographic toy model of Ref. [7], where the dual gravitational theory did not permit thermalization because in this limit there are no propagating degrees of freedom (bulk gravitons). Instead one needs to turn on scalar degrees of freedom (the bulk dilaton) dual to the Lagrangian density [? ].
Finally, we also investigated the case where one subsystem is described by kinetic theory. In order to do so, we supplemented one of the subsystems with microscopic dynamics and chose to describe it in terms of a transport theory of a distribution function of particles f ( x, p, t). We observe that because the metric coupling is mediated by local stretches of space-times of the individual subsystems, all the modes with different p are affected uniformly. Because of this, the only effect the coupling between the two subsystems can have to the distribution function is to rescale the momentum variable f ( x, p, t) → f ( x, p , t) with p i = Λ ij ( x, t)p j . In consequence the coupling cannot bring the distribution to the equilibrium form, and the non-hydrodynamic modes are unaffected by the coupling. Therefore, as in the toy example studied in [7] the thermalization of the full system must rely on thermalization of the individual subsectors, in the sense that coupling a non-dissipative subsystem to a dissipative one does not allow the non-dissipative subsystem to thermalize. Considering the damping of non-hydrodynamic relaxation modes we find that they receive typically small modifications, or none, through the mutual effective metric coupling.
In order to have a conserved total energy-momentum tensor we need (for simplicity switching to a Minkowski background metric g (B) µν for now) where the Ward identities for the subsystems, (2.13) and (2.15), have been used. The terms obtained by differentiating ∆K can be easily seen (using cyclicity of the trace) to be matched by the ansatz .5) and (Indices enclosed by round parentheses are to be symmetrized.) To match terms of order 2k in ∆K we have to take = 1, . . . k, m = 0, . . . k − , and j i = k − − m. The number of coefficients γ and γ at this order (denoted by |γ k |) is therefore where it is convenient to define |κ −n | = 0 for n = 1, 2, . . .. (Note that |κ 0 | = 1, corresponding to the possibility of adding a cosmological constant, which we ignore because it does not lead to a gravitational coupling of the two sectors.) For the first few orders 2k the number of coefficients in the ansätze for ∆K and the two metric tensors read From differentiating the trace terms with powers j i at order 2k we get q(k) + q(k − 1) + . . . + q(1) different terms, while from differentiating (tt) m with m ≥ 1 we get p(k − 1) + p(k − 2) + . . . + p(0) = |κ k−1 | different terms (here p(0) = 1 corresponds to the single term where m = k). Now it turns out that 15 q(n) = p(n − 1) + p(n − 2) + . . . + p(0) = |κ n−1 | and therefore the number of relations is 2|κ k−1 | + |κ k−2 | + . . . = |γ k |. Hence, there are always as many relations as coefficients in the ansatz for the metric which can be used to determine them in terms of the (free) coefficients κ in ∆K.
It is interesting to note that the most general form of the conserved energy-momentum tensor that we obtain here is that where the interaction term is the most arbitrary symmetric polynomial of t µ ν andt µ ν , and (i) is thermodynamically consistent and (ii) such that the total entropy is the sum of the two subsystem entropies. We show in Appendix B that the latter requirements can indeed be satisfied if the interaction term is proportional to δ µ ν as we have here and is otherwise arbitrary. The combinatoric identity (A.8) along with our general construction ensures that any such arbitrary interaction term in the full conserved energy-momentum tensor satisfying the above thermodynamic requirements can be interpreted as a democratic effective metric interaction, i.e., as an interaction that can be absorbed into an appropriate mutual modification of respective effective metrics. This is important because in absence of any symmetry argument to rule out a specific interaction term, it can be generically present, and therefore indeed we should be able to obtain it via an mutual effective metric coupling.

A.1 Solutions to lowest orders
When ∆K is a quadratic expression formed of the energy-momentum tensors, there are two coefficients each in ∆K and the metric ansatz, with two equations relating them: where κ ≡ κ 010 and κ ≡ κ 10 in terms of the general multi-index coefficients introduced above, with0 denoting an infinite string of zeros; where γ ≡ γ 1|0 and γ ≡ γ 1|0 . Eq. (A.4) yields in accordance with (2.18). At the next higher order, there are 4 coefficients in ∆K, κ 20 , κ 110 , κ 020 , κ 0010 , (A. 12) and 5 coefficients in the metric tensors g µν andg µν , constrained by 5 linear relations following from (A.4), which yield

A.2 Action formulation
When the subsystems can be described by an action principle, one can formulate the democratic effective metric coupling of the two subsystems also through a joint action. Let the fundamental elementary fields of the first subsystem be denoted collectively as φ and those of the second subsystem asφ. Then the dynamics of the full system with the lowest-order effective metric coupling (A.10) can be obtained from (A.14) We note that the above action is not simply a functional of φ,φ and g µν as usually is the case but also of the two effective metrics g µν andg µν . Thus g µν andg µν appear as auxiliary fields and the interaction terms of the two subsystems represented by the last two lines of the above action merely implement the algebraic relations between the effective metrics and the subsystem energy-momentum tensors. On the other hand, if we vary with respect to φ andφ first, we evidently obtain the two subsystem dynamical equations in the respective effective metrics, since the last two lines are independent of φ andφ. Since the individual subsystem actions are diffeomorphism invariant, these automatically imply that the subsystem Ward identities ∇ µ t µ ν = 0 and∇ µt µ ν = 0 hold on-shell. We can also explicitly check that when the full action is stationary for variation of φ,φ, g µν andg µν , then its variation with respect to the background metric g (B) µν yields the full energy-momentum tensor T µν as given in Section 2. 16 B Thermodynamic consistency

B.1 General proof
We will show that any consistent effective metric coupling rule with a total conserved energymomentum tensor of the form (2.22) will imply that if the global equilibrium condition is satisfied then thermodynamic consistency is also satisfied with total entropy S(T ) = s 1 (T 1 )b 3 + s 2 (T 2 )b 3 . Note that this not only covers the simplest metric coupling rule (2.16) but the most general ansatz for ∆K in (A.2). 16 The existence of a full conserved energy-momentum tensor in any case follows from the individual subsystem Ward identities (which are easier to obtain from the above action as shown above). Typically a system has a unique energy-momentum tensor of a system up to possible local terms that are separately conserved by virtue of identities (i.e. without need of equations of motion). We can explicitly check that (A.14) does not produce such additional terms.
In order to prove this assertion, it is useful to consider the system in a static gravitational potential so that the background metric is: where φ(x) is static, i.e., not a function of time, but otherwise arbitrary. One can then make static ansätze for the effective metrics of the individual sectors, i.e. assume that The effective metric coupling equations (constructed after removing anomalous terms in the individual energy-momentum tensors) with the assumption that each sector is in thermodynamic equilibrium at temperatures T 1 and T 2 , respectively, i.e. with the forms involve no derivative of the metric so will still be algebraic (although the solutions will depend on the specific spatial point x). However, the coupling equations need to be taken in their generalized form with nontrivial background metric g (B) . We will assume that we can obtain flat-space solutions by smoothly taking φ(x) → 0.
To define an equilibrium solution, we need to relate T 1 and T 2 parametrizing each in terms of the full system temperature T . At global equilibrium, the inverse of the Euclidean time circle specifying the system temperature should be constant. The global equilibrium condition is simply then: where T 0 is a constant, parametrizing the global thermal equilibrium of the full system in the background metric (B.2). The first thing that we need to show is that the above is compatible with the conservation of the energy-momentum tensors. Using (2.13), we can check that the conservation of the individual thermal energy-momentum tensors (B.4) in the respective effective metrics (B.3) imply simply that respectively. Note that b(x) andb(x) do not feature directly in the above equations. Since dP 1 = s 1 dT 1 , 1 + P 1 = T 1 s 1 , dP 2 = s 2 dT 2 , 2 + P 2 = T 2 s 2 , the conservation equations (B.6) are equivalent to ∂ i (ln(T 1 a)) = 0, ∂ i (ln(T 2ã )) = 0, (B.7) and thus implied by the global equilibrium condition (B.5).
By construction, the effective metric couplings ensure that the total energy-momentum tensor, which can be parameterised as T µν = diag E(T (x))e 2φ(x) , P(T (x)), P(T (x)), P(T (x)) , (B.8) will be conserved in the background metric (B.2), since the individual thermal energymomentum tensors are conserved with respect to the respective effective metrics. We therefore have Equations (B.7) and (B.5) imply that and therefore Identifying E + P = T S this implies Since the above should hold for arbitrary smooth φ(x), we conclude that where the variation is taken by changing the constant parameter T 0 . Together with E +P = T S the above implies dE = T dS. (B.14) This shows that thermodynamic consistency follows from the conservation of the full energymomentum tensor as ensured by our effective metric coupling. In particular, assuming E + P = T S and the global equilibrium condition (B.5), we obtain dE = T dS from the conservation of the full energy-momentum tensor. Clearly, we can take the limit φ(x) → 0 limit to obtain the desired proof of thermodynamic consistency in flat space. We still need to show which form S takes. Since the form of the full energy-momentum tensor with one contravariant and one covariant index is such that the explicit interaction terms involving ∆K are always diagonal, they come with opposite signs for E and P. Therefore, Thus from (E + P) = T S we get The relation between the temperatures (B.5) reduces this to The above form then holds for the general consistent effective metric coupling discussed in Appendix A. Thus, we obtain a general proof of thermodynamic consistency with (B.5) and the above form of the full entropy.

B.2 Explicit check
In the following, we explicitly verify thermodynamic consistency of the equilibrium solution obtained in Sec. 3 for the simplest coupling rules (2.16).
With the results (3.6), the thermodynamic relation E + P = T S is evidently fulfilled with T = T 1 a = T 2ã and S = s 1 b 3 + s 2b 3 , when 1,2 + P 1,2 = T 1,2 s 1,2 . Here we shall check that then also S = dP dT (B.18) holds provided the two subsystems satisfy We need to evaluate Differentiating the equations for the metric factors allows us to substitute the derivatives of the parts written within curly brackets as follows: where a prime means differentiation w.r.t. T (except in the case of γ ). This leads to dP dT = d dT P 1 ab 3 + P 2ãb 3 + 1 a b 3 + 2ã b 3 − 3P 1 ab 2 b − 3P 2ãb 2b = P 1 ab 3 + ( 1 + P 1 )a b 3 + P 2ãb 3 + ( 2 + P 2 )ã b 3 = dP 1 dT 1 The two expressions within parentheses in the last line are both dT /dT = 1, which completes the proof: dP/dT = S.

C Low and high temperature behavior of the equilibrium solution for conformal subsystems
In this appendix we collect some formulae that allows one to derive analytically the behavior of the equilibrium solution for conformal subsystems (with arbitrary n 1 , n 2 ) at low and high temperatures, or, equivalently, small and large coupling γ, in the case r = −γ /γ > 1 such that solutions exist for all values of the physical temperature T . Moreover, we show how conformality in the limit of large γT 4 comes about. For small γT 4 , a power series expansion of the solutions to the set of equations (3.20) can be easily obtained. The leading terms in the metric coefficients are (in the latter the dependence on r first shows up at third order). As discussed in Section 3.3, the lightcone velocities asymptote to finite values v ∞ ,ṽ ∞ for large γT 4 , provided r > 1. These values are obtained by solving the sixth-order algebraic equations (3.28), which reduces to a quadratic equation with solution (3.27) when n 1 = n 2 .
The full, nonperturbative equation determining the lightcone velocities as a function of γT 4 is given by (3.22) and (3.23) which were obtained by solving first the quadratic equations for b 2 andb 2 that are implied by (3.20). Using At small γT 4 , all metric coefficients as well as v andṽ tend to unity, with 1 − v 2 and 1 −ṽ 2 proportional to γT 4 . As one can check easily, (C.3) confirms the first-order coefficients in (C.2).
At large γT 4 , where v andṽ approach nonvanishing values v ∞ andṽ ∞ below unity, (C.3) implies that the metric coefficients a,ã, b,b grow linearly with physical temperature T . Since the effective temperatures of the subsystems are given by T 1 = T /a and T 2 = T /ã, this means that they saturate at finite values proportional to γ −1/4 , This behavior of the metric coefficients, together with saturation of t µ ν andt µ ν , implies that at large T the coupling rules (2.16) become From the full solution we in fact find that T µ µ /T 4 ∼ γ −1/2 T −2 . Note also that T µν should be interpreted as the non-anomalous part of the full energymomentum tensor which is locally conserved by itself. (In flat background, the anomalous contribution vanishes.) Therefore, indeed in the limit γT 4 → ∞ the full system becomes conformal provided each system is conformal individually. As a corollary, the fluctuations of the full system in the limit T → ∞ with γ and γ fixed will behave as that of a conformal system.

D A new kind of second-order phase transition and its critical exponent
In this section we analyze further the second-order phase transition that occurs at a particular value of r = −γ /γ and derive the value of the critical exponent α in the specific heat of the full system, when T → T c , for the case of two conformal subsystems (3.17) where S is given by (3.18). Let us first study the simplest case of identical subsystems, n 1 = n 2 = n, where we can equate 17 v =ṽ. In place of (3.22) and (3.23) we then have the simpler relation 17 It is a priori not excluded that there are solutions v =ṽ despite having set n1 = n2. Indeed, such a spontaneous symmetry breaking happens in the region γ < 0 such that the broken phase has lower free energy. However, this region is unphysical in that the effective lightcone velocity of the subsystems can be larger than the speed of light in the physical Minkowski space. We have checked numerically that such symmetry breaking does not occur in the case γ > 0. which is plotted in Fig. 15 for the value of r where the second-order phase transition occurs and two values nearby in the crossover and in the first-order regime. Note that in this plot only the part connected to v = 1 (corresponding to γT 4 = 0) is physically realised; increasing γT 4 from zero to infinity lowers v to a finite limiting value given by the zero of [3 + v 4 − 3r(1 − v 2 ) 2 ] which we have given in (3.27).
For r > r c , T 4 is a monotonic function of v between v ∞ and 1, whereas for r < r c it has multiple extrema determined by the dT 4 /dv = 0. This equation can be written as The maximum of r(v) for 0 < v < 1 determines the critical value r c beyond which no extrema of T 4 (v) occur. This is given by dr dv = 4v v 2 + 3 5v 8 + 20v 6 − 202v 4 − 60v 2 + 45 The bi-quartic polynomial factor in the numerator has a single real root in the range 0 < v < 1. It can be given in closed form and reads The critical exponent α in the specific heat (D.1) can now be inferred from the simple relationship (3.18) between entropy and effective lightcone velocities. In the vicinity of the critical point we have, for n 1 = n 2 , As we have seen, the critical point is determined by the simultaneous vanishing of the first and second derivatives of T 4 as given by (D.2) with respect to v. Hence, up to some constant prefactor, and thus In the case of two conformal systems with n 1 = n 2 , the critical parameters can no longer be obtained in closed form. However, one can show that the critical exponent α is independent of n 2 /n 1 and only the values of r c and T c change.
In this case, one has to solve the two equations the zeros of dT /dv and dT /dṽ have to occur simultaneously in general. A critical endpoint with second-order phase transition appears when two zeros of dT /dv (or dT /dṽ) merge as r → r c from below, such that also d 2 T /dv 2 vanishes and a saddle point (in one dimension) arises. In principle, such a saddle point could have the next two higher derivatives vanish, too, which would change the critical exponent α to −4/5. However, with the one additional free parameter n 2 /n 1 there is not enough freedom for a corresponding fine-tuning. By exploring the solutions numerically, we have found that the situation analyzed above for n 1 = n 2 is indeed generic. Also for any n 2 /n 1 = 1, there is a range 1 < r < r c where there is a first-order phase transition at a finite value of γ 1/4 T , a second-order transition at r = r c , and an analytic crossover for r > r c . In fact, r c depends rather weakly on n 2 /n 1 ; it rises from its minimal value (D.6) to ≈ 1.119 for n 2 /n 1 = 0.1 or 10, and can be shown always to remain below 5 4 .