Multiﬂuid Modelling of Relativistic Radiation Hydrodynamics

: The formulation of a universal theory for bulk viscosity and heat conduction represents a theoretical challenge for our understanding of relativistic ﬂuid dynamics. Recently, it was shown that the multiﬂuid variational approach championed by Carter and collaborators has the potential to be a general and natural framework to derive (hyperbolic) hydrodynamic equations for relativistic dissipative systems. Furthermore, it also allows keeping direct contact with non-equilibrium thermodynamics, providing a clear microscopic interpretation of the elements of the theory. To provide an example of its universal applicability, in this paper we derive the fundamental equations of the radiation hydrodynamics directly in the context of Carter’s multiﬂuid theory. This operation unveils a novel set of thermodynamic constraints that must be respected by any microscopic model. Then, we prove that the radiation hydrodynamics becomes a multiﬂuid model for bulk viscosity or heat conduction in some appropriate physical limits.


Introduction
The hydrodynamic modelling of dissipative systems should guarantee the stability of the homogeneous perfect-fluid states under perturbations. Furthermore, a realistic model for dissipation should be presented in a form which enables an unambiguous contact with microphysics: this requirement, although not strictly necessary from the mathematical point of view, allows for a clear implementation of microscopic inputs into the macroscopic description. For the case of relativistic fluids the problem of finding a theory which fulfils both the requirements is still considered unsolved [1].
The natural-looking relativistic generalization of the Navier-Stokes equations [2], which maintains direct contact with the common notions of viscosity and heat conduction, was shown to admit runaway solutions when the homogeneous perfect-fluid states are perturbed [3,4]. This was shown to be a consequence of the fact that its equations do not admit a well posed initial-value problem [5]. Hence, this hydrodynamic description of relativistic viscous fluids is not suitable for computational applications. On the other hand, the second-order theory of viscous fluids of Israel and Stewart [6] introduces some phenomenological coefficients which have a clear microscopic interpretation only in the ideal relativistic gas limit. Moreover, there are cases in which the model of Israel and Stewart underestimates the number of non-equilibrium degrees of freedom, so that it cannot be considered a universal approach to model relativistic dissipative fluids [7].
The multifluid formalism of Carter [8] may be the solution to this long-standing problem. Partially arising from an action principle, Carter's variational approach leads in a natural way to a well-posed initial value problem governed by hyperbolic equations [1]. Thus, its mathematical structure has all the required properties to give rise to a causal and stable theory if the appropriate equation of state is assumed [9].
Since the multifluid concept was explicitly developed to describe relativistic conducting media [10], Carter's formalism was successfully used as the natural scheme for modelling superfluidity in a covariant framework [11][12][13]. In this context, the correspondence with microphysics was completely established [14][15][16][17]. The formalism was applied to the study of the structure of superfluid neutron stars [18][19][20] and to the formulation of a relativistic theory of vortex dynamics [21][22][23] and represents a fundamental tool in relativistic modelling of pulsar glitches [24][25][26]. In addition, it was recently proposed that a multifluid approach might find interesting application in cosmological models [27,28].
For the case of relativistic dissipation, a general multifluid model was proposed by Carter [29]. This model satisfies the conditions of stability and causality (that coincide with the ones of the Israel and Stewart [6] formulation) for small deviations from equilibrium [30]. Therefore, apart from superfluidity, relativistic multifluids were also studied in the context of heat conduction [31][32][33]. These models, however, still lack of a clear connection with microphysics and are thus not fundamentally preferable to the one of Israel and Stewart [6].
Recently, by using arguments of non-equilibrium thermodynamics, it was shown that any bulk-viscous fluid can always be described as a Carter's multifluid if an appropriate choice of thermodynamic variables is adopted [7]. Therefore, at least for the case of dissipation due to bulk viscosity, this result represents a formal justification of the universality of the multifluid formalism and provides a technique for connecting the hydrodynamic model with thermodynamics and kinetic theory.
In this paper our aim is to provide further insight into the connection of the multifluid theory with microphysics and its universal applicability: we already discussed the link with the equilibrium thermodynamics of a superfluid [17] and the link with kinetic theory for non-conducting bulk-viscous fluids [7]. Here, we add another piece to the global picture by studying how to model a perfect fluid interacting with a radiation fluid within the Carter multifluid framework. Due to its simplicity and wide applicability in astrophysical contexts, this system was widely studied in the literature [34] and the understanding of its properties can be considered satisfactory at every level: statistical mechanics [35], kinetic theory [36], thermodynamics [37] and hydrodynamics [38]. Therefore, Carter's theory (which may provide a universal hydrodynamic framework) should be able to capture the essential physics of this system. In this sense, investigating the properties of a fluid interacting with radiation in this multifluid framework represents a fundamental test for its descriptive power.
In addition, it is well known that radiation hydrodynamics admits a diffusion-type limit which reduces the theory to a conventional model for heat conduction [39,40]. This implies that the matter-radiation fluid may be the first realistic heat-conducting fluid to be rigorously described in Carter's framework, giving us precious insights about the microscopic origin of the so-called entrainment coupling (a non-dissipative coupling between the species in a multifluid) and the correct implementation of the dissipation coefficients.
Throughout the paper we adopt the spacetime signature (−, +, +, +) and work in natural units c = G = k B = 1. Moreover, the symbol γ is used to label the quantities related to the photon fluid and should not be interpreted as a spacetime index.

Multifluid Hydrodynamics
We briefly review the basic ideas of the multifluid formalism. The general theory was formulated in [8,10], see also [1] for a review. The formalism was also extended to incorporate shear viscosity [29] and elasticity [41] but these effects will not be considered in the present work.

Non-Dissipative Evolution of Relativistic Multifluids
The variational approach of Carter builds on the assumption that one can identify a set of four-currents n ν x describing different flows in the system. Given the scalars n 2 xy := −n ν x n yν , an equation of state for the fluid must be provided in terms of a Lagrangian density The condition x ≤ y is imposed to avoid repeated arguments, since n 2 xy = n 2 yx . Introducing the bulk coefficients and the anomalous coefficients we can define the conjugate momenta of the currents as The anomalous coefficients (whose presence in a multifluid is the norm) incorporate the entrainment effect, a non-dissipative coupling between the currents. Historically, the importance of entrainment was first recognized in the context of superfluid mixtures [42,43], but it is a general feature of Carter's variational approach.
In the literature it is common to find an alternative procedure of differentiating the Lagrangian density which includes also the terms with x > y and treats n 2 xy and n 2 yx as independent variables, see e.g., [8,22]. We discuss the connection with the present approach in Appendix A.
A non-dissipative hydrodynamic model can be obtained by considering an action of the form where R is the scalar curvature and √ −g is the square root of the absolute value of the determinant of the metric. The domain of the action is set imposing that the currents are conserved, both on-shell and off-shell. To make sure that this is indeed satisfied, the variations of the currents are taken in the Taub form [44] δn ν where the vector field ξ ν x describes an arbitrary infinitesimal displacement of the world-lines of the fluid elements of the species x. The Euler-Lagrange equations are obtained imposing the stationarity of the action with respect to arbitrary infinitesimal variations δg νρ and displacements ξ ν x . The first one produces Einstein's equations, where the energy momentum tensor has the form The scalar Ψ can be interpreted as a generalised thermodynamic pressure and is given by Ignoring the boundary terms which do not contribute to the equations of motion, the variation of the action produced by the displacements ξ ν x has the form where can be interpreted as the force per unit volume acting to the species x. The condition δI = 0 for any independent choice of ξ ν x produces the Euler-Lagrange equations Equations (7), (9) and (14) constitute a system which arises from a well posed action principle and, therefore, are given in the form of an initial value problem [1]. Please note that from (7) and (10), However, taking the divergence of (9), one immediately has the energy-momentum conservation Therefore even in the case in which the forces f x ν were not zero, their sum must vanish, which is Newton's third law. The formalism is easily extended to the case in which there are some currents that are locked to each other. For example, assume that the species x and the species y interact, and are coupled on time-scales much shorter than those we are interested in. In this case, we can take them to be at rest with respect to each other. The motion is still adiabatic, but now it is subject to the geometrical constraint This is implemented by imposing the world-line displacements ξ ν x and ξ ν y to satisfy the constraint Thus, from the variation (12) we find that f x ν and f y ν do not need to vanish separately, but the total force density does, Again, this condition is Newton's third law.

Including Dissipation
The second law of thermodynamics is not automatically provided by the action principle, so that the dissipative terms of the theory have to be supplied in some other way and inserted by hand into the equations of motion. However, the study of the adiabatic regimes presented in the previous subsection can be used as a guideline for a consistent inclusion of these additional terms.
The current conservation (7) in a dissipative regime may not hold, but chemical-type transfusions may be allowed. Therefore, we need to replace Equation (7) with where r x describes the rate (per unit volume and time) of production of the species x. The second law of thermodynamics is implemented by considering an additional current s ν (interpreted as the entropy current) whose production rate must satisfy the constraint The energy-momentum tensor is assumed to maintain the form (10) and we require it to still satisfy Einstein's equations. Now, its four-divergence takes the form where represent the dissipative generalization of the Lagrangian forces f x ν . Comparison with (9) tells that Newton's third law is still valid, but the terms R x ν do not need to vanish separately. The quantities r x and R x ν incorporate dissipation in the theory and have to be modelled according to microphysical arguments. Now, consider Equation (24) for x = s, where we adopted the notation Θ ν := µ s ν for the conjugate momentum to the entropy current. Contracting with the four-velocity and using Equation (25) we find where we introduced the quantity Thus, only the coefficients r x and R x ν for x = s need to be computed from microphysical calculations. Then, r s and R s ν are obtained through the identities (25) and (28). Finally, there is a subtlety we need to remark on. We have introduced the forces R x ν as dissipative contributions, but from (28) we see that this is not strictly necessary. In fact, one may in principle design them in such a way that is guaranteed by construction, without requiring that the forces themselves vanish. Under this condition, no entropy can be produced and the theory is still non-dissipative. We do not consider this possibility in the following. However, we will be forced to come back to discuss this point in greater detail in Section 6.3.

Heat Conduction and Bulk Viscosity
Relativistic models for heat conduction and bulk viscosity naturally arise as particular cases of the general multifluid theory. In this section, we briefly summarize some results that were obtained up to now.

Heat Conduction
Consider a fluid comprised of indistinguishable particles of a single type, whose number four-current n ν is conserved, In the presence of heat conduction, the entropy current s ν is generally not aligned with the particle flux. Therefore, we consider a minimal two-fluid model with two independent currents, n ν and s ν . The Lagrangian density takes the form The conjugate momenta µ ν and Θ ν , to particle and entropy current respectively, are thus we only need to determine R n ν . It is possible to further reduce the number of unknowns by means of geometrical arguments. In the simplest model, proposed by Carter [10], it is assumed that R n ν is a function of the currents n ν and s ν only (i.e., not of their derivatives). Then the force assumes the form where is a transport coefficient to be determined from kinetic theory. Equation (45) can be derived from the fact that by isotropy, R n ν is a linear combination of s ν and n ν and, from the first equation of (43), needs to be orthogonal to u ν n . The positivity of α is ensured by the equation where we have assumed that Θ E and Θ s are positive. In fact, they reduce to the usual notion of temperature at equilibrium, so this is equivalent to assume that the system is sufficiently close to thermodynamic equilibrium (out of equilibrium a rigorous definition of temperature does not exist and only on equilibrium states Θ E and Θ s both coincide with the thermodynamic temperature). Equation (45) models the force R n ν as a viscous friction between the particle current and the entropy current. Clearly, the effect of such a friction is to drive the system towards a state in which the entropy and the particles flow together. Alternative models for R n ν were proposed, which include terms involving also the derivatives of the hydrodynamic quantities [32]. As a result, in this case the force also has a component which is orthogonal to both n ν and s ν . For small deviation from equilibrium the two models coincide and both reduce to the one of [6]. In this section, we adopt the model of [10] for its simplicity, but the possible existence of terms which contain derivatives cannot be ruled out in principle. We will come back to this point in Section 6.3.
In both cases, this system of equations was shown to have the structure of a relativistic Cattaneo equation [31,32,45]. It is given in a form which is naturally hyperbolic, and therefore compatible with causality, and it becomes a good model for the second sound for high frequency perturbations [38].
When perturbations are slow, i.e., evolve on timescales that are longer than the characteristic relaxation time-scale [31] the conventional Navier-Stokes model for heat conduction is recovered. In the Navier-Stokes limit of Carter's model the thermal conductivity coefficient is given by so that Formula (47), for the entropy production, acquires the more familiar form

Bulk Viscosity
Bulk viscosity arises from the fact that the fluid has internal degrees of freedom which go out of equilibrium due to expansion and contraction of the volume elements in the hydrodynamic evolution. These degrees of freedom can always be modelled as additional currents n ν A , A = 1, ..., l − 1, which are locked to the conserved particle current n ν , provided that the volume element is locally isotropic (which implies the absence of shear viscosity and heat conduction) in the particle rest-frame [7]. Hence, a bulk-viscous fluid can always be modelled as a multifluid whose currents n ν , s ν and n ν A are all subject to the geometrical constraint (18). In this case, the Lagrangian density Λ reduces to where U is the internal energy of the fluid. Its differential is where we use the Einstein summation convention for the chemical index A = 1, ..., l − 1. It is easy to prove that where u ν is the (unique) four-velocity of this non-conducting multifluid. Moreover, the pressure and the energy-momentum tensor take the familiar perfect fluid forms The symbols −A A are adopted for the chemical potentials of the species A because the currents n ν A exist only as a parametrization of the out-of-equilibrium states (see section II-B of [7]). In local thermodynamic equilibrium, according to the minimum energy principle [46], they have the value which minimizes the energy at fixed n and s, and this gives rise to the condition For this reason the A A can be interpreted as generalised chemical affinities, justifying the adopted notation.
The collinearity condition (18), valid for all the currents, simplifies the equations of motion considerably. The independent rates which need to be provided by microphysics are the coefficients r A , which near equilibrium, can be expanded to the linear order in the affinities, The (l − 1) × (l − 1) matrix Ξ AB is symmetric as a result of Onsager's principle. The remaining equations of motion, which are needed to completely specify the hydrodynamic evolution, are given by the particle and energy-momentum conservation Equations (56) and (57) can be combined, giving the formula for the entropy production: which implies that Ξ AB is definite non-negative (strictly positive if the ergodic assumption is made, see e.g., [47]). The model we have presented is constructed in a form that is naturally hyperbolic, and therefore it is compatible with the basic requirement necessary for causality and stability. It was shown in [7] that when l = 2 the above model reduces, for small deviations from equilibrium, to the bulk viscosity prescription derived by Israel and Stewart [6], which is known to be (conditionally) causal and stable.
When the hydrodynamic evolution is slow enough compared to the microscopic equilibration timescales, the model reduces to a relativistic Navier-Stokes description of bulk viscosity (see sections II-D and VII in [7]), with a bulk viscosity coefficient given by Here the matrix Ξ AB is the inverse of Ξ AB , while x s = s/n is the entropy per particle and is the equilibrium fraction of the effective chemical species labelled by A. It is, finally, important to remark a subtlety about the dissipation in a multifluid context. As can be seen from the foregoing discussion, in a multifluid approach, heat conduction and bulk viscosity are not implemented directly as small corrections to the stress-energy tensor but they are modelled in a non-perturbative way by introducing further (non-equilibrium) degrees of freedom in the theory. The immediate consequence is that heat conduction (i.e., the flow of energy in the matter's rest-frame) and bulk viscosity (i.e., the non-equilibrium correction to the pressure) are, in a generic multifluid, interconnected (influencing each other at every order, higher than the first [6]) and cannot be completely separated. For this reason, in the present paper, we have introduced the purely heat-conducting fluid and the purely bulk-viscous fluid separately, while in principle a generic multifluid will contain both the processes.

Radiation Hydrodynamics
We show that the equations of the radiation hydrodynamics in the M 1 closure scheme [48][49][50] can be conveniently obtained directly from Carter's multifluid formalism. This alternative derivation provides considerable thermodynamic and geometrical insight.
Our study will be specifically devoted to photon radiation, so that we use the label γ to indicate the quantities related to the radiation fluid. We remark, however, that the following discussion holds in principle also for any kind of radiation which does not carry any conserved charge (for example, it would apply also to the case of a real scalar boson or a Majorana fermion). In a thermodynamic perspective, this condition corresponds to the requirement that the exchange of radiation between matter elements is a pure phenomenon of heat transfer [51] and not a chemical transfusion.
Although Carter's formalism could be used to model also neutrino radiation (which will be studied in detail in future work), such a system would not admit a chemical-equilibrium limit that is a model of heat conduction. Instead, it would become a charge-conducting fluid, where the transported charge is the lepton number. Therefore, the presence of a conserved charge associated with the radiation fluid is the reason neutrino radiation is physically different from the photon case.

The Hydrodynamic Model
We consider systems containing two particle currents, n ν and γ ν (to avoid confusion, the labels n and γ will always be used as chemical labels and never as space-time indices). The first is assumed to be an exactly conserved current, and represents the flow of the matter component of the multifluid. The second, namely γ ν , is the current density associated with photons, whose number is not conserved (it can change in absorption and emission processes), We impose Boltzmann's molecular chaos ansatz [35], namely that the statistical correlations between matter and radiation can be neglected. This allows defining two separate entropy currents s ν n and s ν γ associated with the matter and the radiation, whose sum gives the total entropy current [52]: Please note that the second law requires that but the two entropies do not need to grow separately. To simplify the system we impose that the heat conduction parameters of matter and radiation vanish, namely Finally, we assume that the interactions between matter and the radiation have the form of local collision processes, which occur for sufficiently short times that the statistical average of the interaction term of the microscopic Hamiltonian can be neglected. This allows us to decompose the Lagrangian density into a matter and a radiation part, for which we will adopt a simple separability prescription where ρ is a pure function of n = √ −n ν n ν and s n = √ −s nν s ν n , while ε is a pure function of γ = √ −γ ν γ ν and s γ = −s γν s ν γ . By comparison with (51) we interpret ρ as the internal energy of the matter fluid measured in its own rest-frame, that is identified by the four-velocity u ν n = n ν /n. Therefore, its differential takes the form dρ = Θ n ds n + µdn, where Θ n and µ are the temperature and the chemical potential of the fluid. Analogously, ε is the internal energy of the radiation fluid, measured in the frame defined by u ν γ = γ ν /γ, and its differential has the form where Θ γ is the temperature of the radiation fluid. In the above equation we introduced the notation µ γ = −A γ to recall that in thermodynamic equilibrium the chemical potential of the radiation fluid µ γ must vanish. In particular, if the emission/absorption process is interpreted as a chemical-type reaction between the matter and radiation, the affinity A γ associated with the above reaction is minus the chemical potential of photons. Please note that the reaction (69) is possible only because the radiation does not carry any conserved charge. For neutrino radiation this is no longer the case, due to the conservation of the lepton number. As a result, the neutrino chemical potential does not vanish in chemical equilibrium [53] and the present discussion does not apply.
It is possible to show that the conjugate momenta to the currents s ν n , n ν , s ν γ and γ ν are, respectively, The generalised thermodynamic pressure given in (11) splits into where is the pressure of the matter fluid, while is the pressure of the radiation fluid. Thus, the variational principle presented in Section 2.1 leads to a completely decoupled energy-momentum tensor where are respectively the energy-momentum tensor of the matter and of the radiation fluid. The expressions in (75) indicate that matter and photons are described as two perfect fluids, so that the stress-energy tensor of the radiation fluid is isotropic in the radiation rest-frame. This is exactly the M 1 closure scheme described by Sadowski et al. [49]. Please note that up to this point, we have not made any assumption about the equation of state of the radiation fluid.

The Dissipative Terms
We have seen that the model is constructed with four currents, but the locking constraints (65) imply that there are only two independent four-velocities. Combining this with Equation (25) we find that there is only one independent four-force which needs to be provided by microphysics, which is Therefore, Equation (24) are given by where we have used the conservation of the matter current to remove the term µ ν ∇ ρ n ρ in the first equation. The system above may seem unfamiliar at a first sight, but with a little algebra it can be shown that it is equivalent to Therefore, the covector G ν , representing the dissipative force of the theory, corresponds to the radiation four-force density and we have finally recovered all the basic elements of the radiation hydrodynamics [34].
The equations of motion given in the natural multifluid form (77) provide an immediate insight into the thermodynamic interpretation of G ν . Let us make the orthogonal decomposition Compared to the first equation of (77), and considering that the first two terms are orthogonal to u ν n , we obtain that The first equation implies that Q, the projection of G ν parallel to u ν n , can be interpreted as the heat exchanged or produced by the matter as a result of the interaction with the radiation fluid. This equation also shows us that the rate r s n does not need to be provided by microphysics, because it must coincide with Q/Θ n . From the second equation we see that f ν can be seen as the part of the radiation force which tends to accelerate the fluid element. More directly, this can be seen by projecting the first equation of (78) orthogonally to u ν n : This is nothing but Newton's second law for a matter fluid element, the inertia of which is the enthalpy and which is subject to the action of a pressure force and of the radiation force f ν .
It is interesting to remark that from the thermodynamic point of view, G ν is the local version (per unit space-time volume) of the heat four-vector acting on the matter fluid, in agreement with the covariant definition proposed by [51]. Thus, Q and f ν can be rigorously identified respectively with the heat and the friction (per unit space-time volume) experienced by the matter element. Now, let us turn our attention to the second equation of (77). If we contract it with u ν γ and invoke the decomposition (79) we obtain where we have introduced the Lorentz factor Equation (82) implies that once G ν and r γ are provided by microphysics, r s γ is automatically constrained. Thus, our analysis indicates that out of three reaction-type rates r s n , r s γ and r γ , only one needs to be given as an external input. One of them, say r γ , should be provided by microphysics: the independent degrees of freedom of the model are 10 (i.e., the two scalars s n , s γ and the eight components n ν , γ ν ) but there are 5 conservation laws given by Equations (16) and (61). Hence, there is room for 5 equations of motion. When the functional dependence of a dissipative term is provided by means of microphysics, the relation which defines it (i.e., Equation (21) or (24)) becomes an equation of motion. Therefore, we need 5 independent microphysical inputs to close the system, namely the four components of G ν and the scalar r γ .
There is also a more physical argument to justify why microphysics should provide both the force and r γ : G ν represents the energy-momentum exchange per unit time between the matter and the radiation fluid, but the same exchange may be originated by scattering processes (which preserve the number of radiation particles) and absorption/emission processes (which modify the number of photons). This immediately tells us that the knowledge of the force G ν is not sufficient to constrain r γ .
In the standard approach which is used in the literature (see e.g., [49]) there is no need to provide r γ and all the knowledge about the interaction processes between matter and radiation is incorporated into G ν . This apparent contradiction with the multifluid approach disappears if the radiation fluid is modeled as an ideal ultrarelativistic gas. In fact, under this condition, the relation holds not only as an equation of state, but also as a kinematic identity (i.e., it is valid also out of thermodynamic equilibrium). This implies that if there are two different thermodynamic states (γ , s γ ) and (γ , then they will have also the same pressure (this is commonly enclosed in the statement that the second viscosity coefficient of an ultrarelativistic ideal gas is always identically zero [7,54]). The mathematical implication is that if in our model the degrees of freedom of the radiation fluid are 5 (i.e., γ ν and s γ ), the energy-momentum tensor is degenerate and only 4 independent degrees of freedom have to be specified (i.e., ε and u ν γ ). Invoking the expression (86) for R νρ , the 10 degrees of freedom can be reduced to 9. Therefore, using the equation of motion in the form (78) together with the matter-particles conservation (61), it is possible to obtain a closed system of 9 equations, in agreement with the standard approach.
We remark, however, that in doing so one is implicitly making the assumption that also G ν has the same degeneracy, in particular that it is not affected by deviations of A γ from zero. The conditions under which this assumption is verified will be discussed in Section 4.4. For now we will keep our analysis general, maintaining the general multifluid formulation based on 10 degrees of freedom.

Thermodynamic Analysis of the Dissipative Terms
One of the biggest advantages of working in the multifluid framework (with 10 degrees of freedom) is that it keeps direct contact with the thermodynamics of the system. In this subsection we show how it can be used to derive useful thermodynamic relations which remain hidden in the standard formulation (based on 9 degrees of freedom) discussed in the previous subsection. These relations can be particularly useful in those situations in which the radiation can have a finite chemical potential for a long time (like in scattering-dominated materials) in which the standard approach may be inapplicable if G ν strongly depends on A γ .
To provide a clear comparison with the existing literature it is convenient to work in the rest frame of the matter fluid. To do this we define a tetrad (i.e., an orthonormal basis of the tangent space) e a = e ν a ∂ ν which is comoving with the matter-fluid element, namely e 0 = u n . We use this tetrad to decompose the radiation stress-energy tensor as which are the radiation energy density, the radiation flux and the radiation pressure tensor (or radiative stress) in the matter rest-frame [34]. By comparison with (86), we find that where we have introduced the three-velocity v j = u j γ /u 0 γ . Please note that the pressure tensorP jk can be written entirely in terms ofε and F j : correcting a typo in Equation (34) of Sadowski et al. [49], in accordance with [55], we obtain where F j is the reduced radiative flux and z is the Eddington factor [48], Equation (89) contains both the essence and the limitations of the closure scheme:ε, F j andP jk are respectively the zeroth, the first and the second moment of the radiation specific intensity [34] and we are closing the system by assuming that the last can be uniquely written in terms of the first two. This has also the natural implication that in the reference frame of the matter element, the isotropy is broken only along the direction identified by F j . For this reason, if we use the tetrad e a to decompose the radiation four-force, see Equation (79), it is legitimate to assume that The coefficient χ can be interpreted as the total opacity, a parameter that sets the attenuation rate of the radiation flux in terms of the flux itself. If, now, we promote the three-vector F j to a space-like four-vector through the construction F ν := F j e ν j , we have the orthogonal decomposition Therefore, using simple geometrical assumptions, the 5 independent dissipative terms of the theory were reduced to 3: Q, χ and r γ .
It is now possible to make a thermodynamic study of these terms near equilibrium. We can rewrite Equation (64) using (80) and (82), obtaining the entropy production where we have introduced the relative speed ∆ and the scalar F, defined through the relations Equation (94) shows that the entropy production is the sum of three contributions. Since only the total is constrained to be non-negative, in principle far from equilibrium r γ , χ and Q can have arbitrary sign, provided that all the contributions compensate each other giving r s ≥ 0. However, if we limit ourselves to near-equilibrium situations it is possible to obtain stronger constraints.
Fist of all, we assume that the 3 dissipative terms Q, χ and r γ are functions only of the local thermodynamic state of the multifluid. This implies that in principle they are functions of 5 independent thermodynamic variables (2 identifying the thermodynamic state of the matter, 2 identifying the thermodynamic state of the radiation, 1 identifying the relative motion). In particular, we decide to work with the 5 state variables that turn out to constitute a convenient choice since it is easy to check that the local thermodynamic equilibrium state is given by This makes A γ , F and Θ γ − Θ n the natural variables that can be used to parametrise a displacement of the system from local thermodynamic equilibrium. If we expand the dissipative terms to the linear order in these variables we obtain Since G ν and r γ vanish at equilibrium, there are no zeroth order terms in the above expansions for r γ and Q. Moreover, there are no contributions at the first order coming from F due to the symmetry of the coefficients under a transformation F j −→ −F j . The 7 expansion coefficients are all functions of the thermodynamic properties of matter only, i.e., n and Θ n . Onsager's principle imposes the reciprocal relation (see Appendix C.1 for the proof) After rewriting the entropy production Equation (94) by using the expansion (98), keeping only the second order in the displacement from equilibrium (which also implies Γ nγ ≈ 1) and imposing the positivity for all small deviation from equilibrium, we obtain the conditions and, making use of the reciprocal relation (99), The coefficients χ γ and χ T appear only in higher order terms in the equation for the entropy production (94); for this reason they can be neglected in the present study.
The constraints (99), (100) and (101) hold independently from the details of the matter-radiation interaction, so they can be used to check the thermodynamic consistency of any model of radiation hydrodynamics.

Application: Deriving the Four-Force G ν from Thermodynamic Arguments
The radiation four-force is usually (see e.g., [38,40,56]) computed from the kinetic theory of radiation assuming that i the thermal coefficients obey the Kirchhoff law, ii the scattering is isotropic and coherent, iii the opacities have a grey-body form.
The line of reasoning which leads to an expression for G ν starting from the foregoing assumptions is briefly sketched in Appendix B.
As a first application of our thermodynamic study, we now show that the same form of G ν can be derived directly in a hydrodynamic framework if one requires that i every dissipative process contributes additively to the transport coefficients and the thermodynamic constraints presented in the previous subsection hold separately for every microscopic contribution, ii the degeneracy assumption we presented in Section 4.2, which allows reducing the degrees of freedom of the model from 10 to 9, is fulfilled also by G ν .
Let us assume that this is the case, namely that the coefficients χ and Q can be written as functions of 4 independent state variables only, instead of the 5 in (96). We retain the variables n and Θ n because they identify the state of the matter fluid. From (89) we know that the flux F j and the energy densityε can be used to identify the radiation energy-momentum tensor completely, which in turn constitutes the reduced degree of freedom of the model. It follows that the natural choice of variables now is (n, Θ n , F,ε).
As in the previous subsection, we can perform a linear expansion in the deviations from equilibrium. As we saw, the deviations of χ are of higher order in the model (in Equation (36) F ν is already a fist-order term), therefore we will not analyse them and we will focus our attention on Q. Recalling that the symmetries of the problem impose that the linear corrections in F must vanish, the expansion contains only one term, namely The functionε eq is the equilibrium value of the energy density, which by comparison with (97), iŝ The assumption that the radiation fluid is an ultrarelativistic ideal gas implieŝ whereB(Θ) is the frequency integrated equilibrium intensity, which can be written aŝ The coefficient a R is a constant factor which depends on the type of radiation [40]. The expansion (103) is a particular case of (98). It is possible to relate k with k γ and k T considering that to first orderε This expression was obtained by computing the first-order expansion coefficients ofε from the equation of state (84), together with the fact that (73) defines the Legendre transformation Plugging this expansion into (103) we find which by comparison with the general formula (98), gives the relations We recall that since we are making a linear study, the densities γ and s γ can be identified with those in equilibrium, which for an ulrarelativistic ideal gas satisfy the relation where the specific entropy b R of the radiation gas is a constant (b R ≈ 3.6 for a Bose gas and b R ≈ 4.2 for a Majorana Fermi gas). Therefore, we have proven that the multifluid formulation of the radiation hydrodynamics, which is a theory with 10 degrees of freedom, reduces (for small deviations from equilibrium) to the standard formulation based on 9 degrees of freedom if and only if To complete the reduction to the standard theory and to see the implications of (112) we divide the dissipative processes at the origin of G ν and r γ into three different categories. We call elastic scatterings (e) those processes which conserve the number and the total energy of the radiation particles which are involved (measured in the fluid rest-frame). These processes give no contribution to r γ and Q. The inelastic scattering processes (I) are those in which only the number of radiation particles is conserved. These do not give any contribution to r γ . Finally we have the absorption processes (A), which do not conserve the radiation particle number. Please note that the absorption processes include also the emission processes. In fact, absorption and emission processes are the time reversed of each other and are mediated by the same matrix element, which implies that in thermodynamic equilibrium they must obey the detailed balance. As a result, it is necessary to consider their joint action at a thermodynamic level and they must not be separated. Using our assumption (i), we can, thus, split the dissipative coefficients according to the different contributions as This separation will also result into a kinetic subdivision of the expansion coefficients given in (98) and (103).
From a purely thermodynamic point of view, the constraints (100) and (101) and the reciprocal relation (99) hold in principle only for the total coefficients, not for the separate contributions.
For example we have the constraint χ A + χ I + χ e ≥ 0, but this does not necessarily imply the separate non-negativity of all the three parts (because only the total appears in (94)). Our assumption (i), on the other hand, consists of requiring that all the thermodynamic constraints hold separately, so that in our example which is in agreement with their interpretation as opacities (and, therefore, as inverses of mean-free-paths). Now, if we turn our attention to Onsager's relation (99), for the case of the inelastic scattering we find which using Equation (110), gives k I = 0 .
We have verified that in order for (103) to hold, one should assume that all the scattering processes are elastic, giving Q I = χ I = 0. This is consistent with the assumption of coherent scattering invoked by [38,40,56].
There is a final constraint we can impose on the kinetic coefficients which arises directly from the assumption (103). Let us consider a situation in which all the radiation particles are in their equilibrium distribution, but there is an excess of one particle of momentum q (measured in the rest frame of the matter). This condition can be modelled as a state of the system in whicĥ Ignoring the scattering processes, this particle will have a life-time τ A before being absorbed. The absorption process can be modelled as the action of the radiation four-force for a time τ A , giving the conditions Simplifying the energy of the particle with the aid of (103) and (92) we obtain k A τ A = 1 and Therefore, the multifluid approach is equivalent to the one presented by [38,40,56] if the radiation four-force can be put into the form which is what we wanted to prove.

Deriving the Reaction Rate From Thermodynamic Arguments
If we assume G ν to be given by Equation (120), we are automatically assigning a value to Ξ γT , as a result of the Onsager relation (99). On the other hand, Ξ γγ remains undetermined, therefore we are not constraining r γ completely. In Appendix B, however, we show that the same microscopic assumptions which lead to (120) can be invoked to prove (assuming a non-relativistic relative speed between the matter and the radiation fluid) that where This expression for the rate was adopted in the literature to model systems in which photon-conserving processes are dominant (see e.g., [57]) and can be used to derive a formula for Ξ γγ .
Before doing this, however, as a second application of our formalism, we will prove that (121) is consistent with the Onsager principle. More specifically, we will show that if one assumes that the rate can be written in the generic form and G ν is given by (120), then the Onsager relation (99) demands To do this, we expand γ near equilibrium, and compute the partial derivatives of γ (in equilibrium) from the ideal gas equation of state, obtaining where c R is a constant coefficient (c R ≈ 1.37 for a Bose gas and c R ≈ 0.91 for a Majorana Fermi gas). Please note that the second relation can be easily obtained taking the derivative with respect to Θ γ of the (equilibrium) relation (111). Employing the expansion (125), Equation (123) can, therefore, be rewritten as Comparing with (98) we obtain The Onsager relation (99), then, implies which can be compared with (110). Recalling that k = χ A , we finally obtain Υ = χ A , which is what we wanted to prove. In conclusion, we have shown that (if one negects the compton scattering) the model for radiation hydrodynamics adopted by Sadowski et al. [57] can be translated into the multifluid framework by imposing and that this choice respects all the thermodynamic constraints we derived in Section 4.3. In fact, the only non-trivial inequality which is left to check is (101), which in our case reduces to and is satisfied by both the Bose gas (b R c R ≈ 4.93) and the Majorana Fermi gas (b R c R ≈ 3.83).

Radiation as a Source of Bulk Viscosity
In Section 3 we anticipated that relativistic models for heat conduction and bulk viscosity can be naturally obtained as particular cases of the general multifluid theory. Hence, also any hydrodynamic model which is formulated in a multifluid framework can be interpreted as a heat-conducting or a bulk-viscous fluid whenever it arises from a Lagrangian density of the form (32) or (51) and there is only one strictly conserved current n ν . It is clear that the model for radiation hydrodynamics we presented fails to satisfy the first condition and therefore does not admit a straightforward interpretation as a heat-conducting or as a bulk-viscous fluid. It is possible, however, to impose further constraints, besides those given in Equation (65), to recover the canonical models for dissipation given in Section 3. In this subsection we focus on the possibility of transforming the presence of radiation into a contribution to bulk viscosity.
First, to recover a model for pure bulk viscosity the system should be non-conducting. To implement this physical requirement, we impose the constraint which implies u ν n = u ν γ . This can be obtained at a dynamic level taking the limit χ −→ +∞ and it corresponds to the infinitely optically thick regime in which Under this condition the radiation fluid is completely advected by the matter fluid (i.e., there is no net conduction of photons in the matter frame). Since now matter and radiation have the same rest-frame, then it is possible to define the conglomerate internal energy which by comparison with (75), is the total energy density measured in the common rest-frame. By comparison with (66) we see that the Lagrangian density has the form (51): we have constructed a bulk-viscous fluid. An alternative way of seeing the emergence of bulk viscosity is to start from Equations (67) and (68) and to write dU = Θ n ds n + µdn which can be recast as Making the identifications the differential of the internal energy (136) has the form (52). The equilibrium conditions (97) tell us that the displacement of the matter-radiation multifluid from local thermodynamic equilibrium is given by a non-vanishing value of the generalised affinities A γ and A s n . These affinities are associated with the densities γ and s n , which are in turn interpreted as generalised reaction coordinates. This is in agreement with the equilibrium condition (55) for bulk-viscous fluids. The pressure and energy-momentum tensor, which for the matter-radiation multifluid are equal to (71) and (74), coincide with the ones calculated according to the prescription (54), namely This is due to the fact that the variationally defined energy-momentum tensor (10) is invariant under changes of chemical basis [8].
We can also map the dissipative expansion coefficients introduced in (98) into the reaction matrix Ξ AB presented in Equation (56). With the aid of Equation (80), we can rewrite (98) in the form which is the defining relation for the 2 × 2 matrix Ξ AB . The Onsager reciprocal relation (99) ensures the symmetry condition Ξ AB = Ξ BA , while the thermodynamic constraints (100) and (101) imply the non-negativity of Ξ AB . We can also compute the bulk viscosity coefficient (59) associated with the radiation processes (the details of the calculations are reported in Appendix C.2), that turns out to be (140) The above formula is the general expression for the bulk viscosity of the multifluid for arbitrary values of the kinetic coefficients k T , Ξ γT and Ξ γγ . If we impose the validity of the condition (112), which combined with the Onsager relation (99) implies together with the consequent expression for the four-force (120), we find that the bulk viscosity coefficient simplifies to The coefficient Ξ γγ simplifies and does not play any role in the final expression for the bulk viscosity. This stems from the fact that the condition (112) is imposed to guarantee that the hydrodynamic evolution is decoupled from the chemical evolution of the degree of freedom γ. The term in the brackets in Equation (142) strongly depends on the equation of state of the matter fluid and essentially nothing can be said a priori if ρ(n, s n ) is not specified. In Appendix D we compute ζ explicitly for the case of a non-degenerate gas in the ultra-relativistic and non-relativistic limits.
Equation (142) can be used to replace the multifluid with a Navier-Stokes (or Israel-Stweart) bulk-viscous fluid (the final equations governing these fluids can be taken to be the ones presented in section IV-C or VII of [7]) using the matter+radiation equilibrium equation of state and ζ as the prescription for the bulk-viscosity coefficient.

Radiation Hydrodynamics as a Model for Heat Conduction
To obtain a model for heat conduction we need to consider an opposite situation with respect to the bulk-viscous case. We need to assume that the radiation-particle production rate and the exchange of energy in the rest-frame of the matter are faster than the hydrodynamic time-scale. In this way they are always in equilibrium and do not contribute to the entropy production. This, compared with (94), implies (assuming that the relative speed is non-relativistic) This condition can be formally achieved by sending at least two out of the three coefficients Ξ γγ , Ξ γT and k T to infinity. In the case in which the radiation four-force has the form (120), the limit Ξ γγ −→ +∞ can be safely imposed (even in those cases in which this is not a rigorous assumption) because the value of Ξ γγ does not have any influence on the hydrodynamic evolution of the model based on 9 degrees of freedom (n ν , Θ n ,ε, F ν ), which are usually the variables of physical interest.
For the other coefficients there is the complication that the constraint (119) implies that if k −→ +∞, then also χ will diverge. However, in this case we would have also the locking condition (132) and the multifluid would simply reduce to a single perfect fluid. Since we need to keep χ finite to enable heat conduction, in the following we will simply assume that (143) holds as a mathematical constraint, without specifying under which physical conditions this constraint is respected.

The Origin of the Entrainment
An easy way of studying the implications of the constraint (143) is to analyse its effect on the Lagrangian density. From (5) and (70) we have that the variation of Λ under arbitrary variations of the currents (at constant metric components) is We can use (63) to rewrite this variation in the equivalent form where analogously with the notation (137), we have introduced the covectors Now we impose that the variations of the currents, which in principle may be all independent, satisfy the constraints (65): we can write s ν n = x sn n ν γ ν = y s ν γ = y(s ν − x sn n ν ) and a variation of the first expression reads δs ν n = n ν δx sn + x sn δn ν .
The variation of Λ becomes where we have introduced the momentum Finally, we impose the constraints (143) on the non-varied state (i.e., the original state in which the system is, before we make the variation), but not on the varied state. The meaning of this procedure is that we are assuming a reference state which is solution of Equation (143), but the variations δs ν , δn ν , δx sn and δy are completely arbitrary (the reason why we leave them arbitrary will be clarified in the next subsection). It is easy to see that for non-relativistic relative speeds (143) is equivalent to the chemical-type equilibrium conditions so that all the contributions arising from the arbitrary variations δx sn and δy vanish and (149) reduces to δΛ = Θ ν δs ν +μ ν δn ν .
This proves that the variation of the Lagrangian density is indistinguishable from the one of a heat-conducting fluid, provided that we interpret the momenta Θ ν andμ ν as the conjugate momenta respectively to the entropy and particle current (cfr. with Section 3.1). These momenta, written in terms of the currents s ν and n ν , are respectively which lead to the identification of the bulk coefficients and of the anomalous coefficient A sn , which encodes the entrainment phenomenon, see Equation (5), It is useful, now, to summarize what we have obtained so far. We started with the model (66) which was built considering 4 currents. No entrainment was assumed in this model, namely the conjugate momenta to all the currents were collinear with the respective currents. Then we have reduced the dynamical degrees of freedom to 2 independent currents invoking 2 collinearity constraints (65) and 2 chemical-type equilibrium conditions (151). We have found that the reduced theory, described using only 2 independent four-currents, reproduces a fluid with entrainment (i.e., a fluid in which the conjugate momenta are linear combinations of both the currents, see Section 2.1).
This gives us a deep insight on the fundamental nature of the entrainment. In fact, the anomalous coefficient (155) is proportional to x sn , which represents the entropy per-particle carried by the matter fluid (see (147)), a quantity that comoves with n ν . Thus, we see that the entrainment between two currents may arise also in a theory which originally does not admit it: an effective entrainment coupling emerges whenever a fraction of the constituents of one current is forced to comove with the second current and the processes which tend to alter this fraction are in equilibrium.
As a final remark, we note that for the case of the superfluid Helium, this mechanism for the emergence of entrainment is at the origin of the equivalence between the Tisza-Landau two-fluid model and the multifluid model of Carter and Khalatnikov [15]. In the Landau model it is assumed that the Helium current can be split into two non-conserved and non-entrained currents, one of which (the so-called normal current) is locked with the entropy current. On the other hand, in model of Carter and Khalatnikov this splitting is not explicitly done, but its existence is reflected into a non-vanishing entrainment between the entropy and the total particle current. The steps of the proof of the equivalence between these two approaches are analogous to the calculations performed in this subsection [15].

Energy-Momentum Tensor
The energy-momentum tensor introduced in Section 2.1 can be equivalently defined as where the variation is performed keeping constant the components of the Hodge duals of the currents In Section 4 the energy-momentum tensor was computed treating the currents s ν n and γ ν as independent variables, therefore in the calculation of the derivative (156) their Hodge duals where held fixed, imposing the condition On the other hand, if we want to eliminate these currents from the set of possible degrees of freedom and work with a theory in which the only two fundamental currents are n ν and s ν , we have to impose the constraints (143) also in the varied state. This implies that x sn has to be considered a function of the three fundamental scalars of the model, giving The same argument should in principle hold also for y; however we know from microphysics that in equilibrium y = 1/b R , therefore the condition δy = 0 is left unchanged. However, in deriving the formula of the variation (152) no assumption on the variation of x sn and y was made (they were completely arbitrary). This means that the energy-momentum tensor we obtain from the formula (156) is the same both in the original model with 4 currents, imposing the constraints only at a dynamical level, and in the reduced model with 2 currents, in which the constraints hold also off-shell (and therefore remain valid when the variation is performed, see section 4 of [15]). The implication is that the pressure (71) and energy-momentum tensor (74) can be equivalently rewritten from (152) in the canonical forms (11) and (10): This can also be checked with direct calculations. Thus, we have proven that the energy-momentum tensor of the radiation hydrodynamics assumes the form of the energy-momentum tensor of a heat-conducting fluid, Equation (36), provided that we interpret the entrained momenta (153) as the canonical conjugate momenta to the entropy current and to matter-particle current respectively.
Let us perform the Eckart-frame decomposition (42). We note that the Eckart four-velocity identifies the matter rest-frame. However, the decomposition of the radiation stress-energy tensor has already been performed in Section 4.3. Comparing (42) with (89) we have the straightforward identifications so we see that Eckart's heat-flux is simply the radiation flux. In Appendix C.3 we show that this identification (which we have obtained from the comparison of the energy-momentum tensors) is also consistent with the Eckart-frame decomposition (39) of the entropy current (63). Finally, we can study the second order term in (42), whose coefficient D, introduced in (41), is given by where we have used the second equation of (154). It is easy to prove that the pressure tensorP jk introduced in Section 4.3 can be written aŝ which provides, in the case of radiation hydrodynamics, an immediate microscopic interpretation of the phenomenological second order term appearing in (42) as the anisotropic contribution to the radiation pressure tensor.

The Hydrodynamic Equations
The equations of motion (77) describe the dissipative interaction between the currents n ν , s ν n and the currents γ ν , s ν γ . To complete the construction of the model for heat conduction we need to recast these equations into the form (43), which describes the dissipative interaction between n ν and s ν .
Before doing this, however, there is an important remark to make, which was pointed out in [8]. Let us consider a generic multifluid and define the forces They might be thought to constitute a (1, 1) tensor in the chemical species index x. In fact, if we change the fundamental currents of our theory through a change of chemical basis, i.e., a transformatioñ where the coefficients N y x are some constants, then the conjugate momenta transform according to the contravariant law and the forces (164) will consequently have a mixed transformation law. Since, from (24), we have that the dissipative forces represent only the diagonal part of the tensor R y x ν . As a result, after a change of basis the new forcesR x ν will not be linear combinations of the old forces R x ν only, but the summation will involve also off diagonal terms R y x ν with x = y. This has two remarkable consequences. The first consequence is that if we impose R x ν = 0 for all x in the non-dissipative limit, in principle this will cease to hold if we change the chemical basis (R x ν = 0), due to the presence of mixed terms R y x ν = 0, for x = y. This shows that in general, there is no way to guarantee that the forces R x ν vanish in a non-dissipative theory, unless one has a microscopic argument to support the choice of a preferred chemical basis in which this occurs. This is related to the problem we presented at the end of Section 2.2, namely the fact that one can impose the vanishing dissipation condition (30) even in a context in which R x ν = 0. The second consequence is that since the terms R y x ν contain derivatives (both in space and time), if one imposes that the forces R x ν depend only on the value of the hydrodynamic fields in the point and not on their derivatives, this will be no longer true when we change chemical basis. Therefore, as we have already pointed our in Section 3.1, one cannot impose that the forces do not depend on the derivatives of the hydrodynamic fields without a microscopic justification. Now, we can study the heat-conductive limit of radiation hydrodynamics, knowing that both these complications may arise (in fact, we are going to perform substantially a change of chemical basis).
Let us take the second equation of (77) and use the second constraint of (143) to remove the terms proportional to A γ ν : From the second equation of (43), we know that Recalling the first definition of (146), we can use (168) to rewrite R s ν in the form We can further simplify this expression by invoking the decomposition (93), which recalling the second equation of (146), gives us the final formula The equation for the entropy production (94), on the other hand, reduces in our case to Therefore, we see that if χ = 0 the theory is non-dissipative, even if R s ν = 0. Furthermore we note that having imposed that χF ν depends only on the value of the hydrodynamic variables in the point, the force R s ν contains terms which are linear in the derivatives. We also note that Equation (171) describes a force which is not necessarily a linear combination of n ν and s ν , but contains a component which is orthogonal to both. This is in agreement with the discussion in [58].

Heat Conductivity Coefficient
We conclude the section by calculating the coefficient of thermal conductivity. This is conveniently done by comparing the formulas for the entropy production (50) and (172). Imposing the equality of the two, we find the condition (neglecting the Lorentz factors) However, from (A40) and (A41) we obtain which plugged into (173), gives It is well known that the radiation hydrodynamics have a diffusion-type limit which makes it analogous to a phenomenon of heat conduction. Shapiro [39] and Farris [40] have computed explicitly the corresponding coefficient κ, obtaining the formula (175). We have generalized this result, showing the complete formal analogy between the two systems in the framework of the multifluid formalism.

Limitations of the Model
We conclude with a few comments about the limitations of our model. We have seen that assuming the Lagrangian density (66) as a starting point for Carter's approach implies that the matter and the photons are both described as two perfect fluids. The fact that the matter can be modelled as a perfect fluid is justified if the collisions between matter-particles are faster than the hydrodynamic time-scale. However, the same argument cannot be applied to the radiation fluid, whose particles typically do not interact with each other and, therefore, an H-theorem for the radiation gas alone does not exist. This implies that the closure scheme cannot, in general, be justified using thermodynamic or kinetic arguments. Indeed, it is well-known that the local properties of the radiation stress-energy tensor depend on the global structure of the radiation field (in particular on the disposition of its sources) and, as a consequence, the M 1 closure scheme fails the multiple-source shadow test [49].
Given the fact that it is not possible to justify the closure scheme as a physical limit (and therefore it is not guaranteed to provide an accurate description of reality), we can understand its appearance in the multifluid theory as a byproduct of applying the principles of information theory in the context of Carter's formalism. In fact, one is required to provide a limited set of macroscopic parameters (the radiation particle total current γ ν and the rest-frame energy density ε, both appearing in the Lagrangian density (66)) and all the remaining properties of the radiation field need to be written in terms of this limited (local) information. Then, following the philosophy of information theory, well summarised by Jaynes [59], we have to assume that the system is in the microstate that maximizes the entropy (or, equivalently, minimizes the information) compatibly with the values of the macroscopic parameters which are known. Thus, denoting the microscopic single-particle-state occupation numbers by N(p), the Shannon entropy s γ , which for an ideal gas is the opposite of Boltzmann's H-function, is given in the radiation rest-frame by where j = −1 for Bosonic radiation and +1 for Fermionic radiation, g is the spin degeneracy and h p is Planck's constant. The particle and energy density are Hence, the most probable state must be obtained imposing where α and β are two Lagrange multipliers. This operation gives the Bose-Einstein/Fermi-Dirac occupation which justifies the interpretation of s γ , Θ γ = β −1 and −A γ = αβ −1 as respectively the thermodynamic entropy, temperature and chemical potential of the radiation gas.
To summarize, the multifluid formalism forces us to assume that the currents and the Lagrangian density are the only information we are given (it is our macroscopic knowledge about the system), and this leads to a fluid-model of the radiation gas.
Apart from failing the multiple-source shadow test, this approach has also clear limitations when the opacity has a strong dependence on the frequency. In fact, under this condition, the expression (120) for the radiation four-force ceases to hold and a hydrodynamic approach may be inapplicable. In this situation one should follow the evolution of each radiation frequency separately, requiring a kinetic theory approach.
Furthermore, even in the case in which there was only one relevant radiation frequency, which in principle may still allow the use of a hydrodynamic treatment, the M 1 closure scheme of Sadowski et al. [49] would be inconsistent with the given information. This was shown in the context of information theory by Minerbo [60], who worked in the rest frame of matter and considered a monochromatic radiation flux with a given frequency (measured in the matter's frame). By choosing the energyε and the components radiation flux F j as basic information about the system, he found a maximum entropy distribution which produces an energy-momentum tensor which obeys to a different closure scheme (see also [48] for a comparison between the different closure schemes).

Conclusions
We have studied how radiation hydrodynamics can be modelled in the context of Carter's multifluid formalism. The radiation stress-energy tensor was found to obey to the M 1 closure scheme presented by Sadowski et al. [49] and the hydrodynamic equations were shown to be equivalent to those which are often employed in the literature [34]. Moreover, we connected the hydrodynamic theory with non-equilibrium thermodynamics and performed an Onsager analysis of the dissipative terms of the model.
As an immediate application, we showed that the grey-body radiation four-force [38,40,56] is the only thermodynamically consistent expression for the force between the matter and the radiation fluid which can be used in a model with 9 independent degrees of freedom.
The multifluid formalism, therefore, perfectly captures and describes the physics of radiation hydrodynamics in detail, offering novel insight into a well understood subject.
In the second part of the paper, we reinterpreted radiation hydrodynamics as a particular case of relativistic dissipation and we used this reformulation to gain new understanding in the latter. In the infinitely optically thick limit, the interaction between the matter and the radiation fluid was shown to become a source for bulk viscosity. This is in accordance with the more general result that any locally isotropic fluid is a Carter bulk-viscous multifluid [7].
In the opposite limit, in which the radiation fluid is assumed to be in chemical equilibrium with respect to particle production processes and to have a rest-frame temperature equal to the one of the matter fluid, the multifluid reduces to a model for heat conduction. We found that the entrainment between the entropy and the matter current arises naturally from the original splitting of the entropy current into a matter part and a radiation part which were not entrained. When we used the condition of equal temperature to reduce the number of degrees of freedom of the model, the momentum associated with the entropy of matter was naturally divided into two parts which were redistributed between the currents, generating an effective entrainment in a theory that originally (because of the assumed form (66) of the Lagrangian density) was entrainment-free.
The equations of motion of the resulting heat-conducting fluid have, in the parabolic limit, the well known diffusive form of the radiative transport [39,40] and we verified the correspondence of the respective transport coefficients. In the hyperbolic regime, however, beyond first order in the deviations from equilibrium, the model was shown to be different from all the proposed universal models for heat conduction [6,10,32]. In fact, the structure of the hydrodynamic equations (determined by the expression of the four-force acting between the currents) preserves many details of the physics of the radiation hydrodynamics and no universal available theory for the heat conduction is so general to be able to encompass all these details. This paper constitutes a further step forward towards the global unification of the relativistic hydrodynamics, showing with a concrete example how multifluids and dissipative single-fluids can arise as two different mathematical descriptions of the same theory. In particular, our study might constitute the beginning of the construction of a bridge between the hydrodynamic models employed in simulations of super-novae (where a multifluid approach is usually adopted) and those employed in simulations of neutron-star mergers (where a single-fluid approach is preferred). Funding: This research received no external funding. whose differential, recalling (3) and (4), is This representation of the equation of state gives Λ as a function of the upper triangle x ≤ y of the matrix n 2 xy and for this reason it may be called triangular formulation. An other approach consists of considering Λ as a function of the whole square matrix (and for this reason we can call it square formulation) trough the equation Λ(n 2 AA , n 2 AB , n 2 BA , n 2 BB ) := Λ n 2 AA , where in the right-hand side we are using the functional dependence of Λ presented in Equation (A1). Equations (A1) and (A3) describe the same physical quantity, in fact in any real state n 2 AB = n 2 BA . However in the square formulation n 2 AB and n 2 BA are treated in the equation of state as independent variables. This allows writing the differential of Λ in the more compact form where the coefficients K xy form the symmetric 2 × 2 matrix From (A4) it is immediate to prove (5), in fact one can easily see that leading to an explicit expression for the momenta in the matrix form which is equivalent to Equation (5). The distinction between the triangular and the square formulation is never explicitly discussed in the literature and the two are used interchangeably according to convenience. However it is important to keep the distinction clear in mind, because in the square formulation while in the triangle formulation There is no contradiction between the two, because in the first case the derivative is performed keeping n 2 BA constant, while in the second case it is performed along the curve n 2 AB = n 2 BA , producing a double variation of Λ.

Appendix C. Calculations
In this appendix we report in detail the calculation which were omitted from the main text.

Appendix C.1. Onsager Symmetry of the Dissipative Coefficients
We consider a homogeneous matter-radiation multifluid prepared in an initial state such that both the components are at rest. Under this condition we can combine Equations (67) and (68) to write the differential of the total entropy as Imposing the conservation of energy and particle number, we obtain the differential Therefore we are able to identify the thermally fluctuating variables and their conjugates Onsager principle states that if we write the evolution of the variables y A in the forṁ then L AB is symmetric. From (62) and (80) we find that in homogeneous configurations and in the absence of a relative flowρ = Qγ = r γ , which recalling (98), gives the formulaṡ It is easy to rewrite this system in the form (A28): which produces the reciprocal relation k γ = Θ n Ξ γT .

Appendix D. Radiation-Mediated Bulk Viscosity of Non-Degenerate Ideal Gases
Let us consider Equation (142) and let us assume that the matter fluid is a non-degenerate ideal gas. Since all the coefficients in the formula for ζ are computed in equilibrium, in this appendix we will impose Θ n = Θ γ = Θ and s γ = b R γ. Therefore we can write where w = 3/2 if the matter-fluid is a non-relativistic gas and w = 3 if it is an ultra-relativistic gas. Now, we immediately see that in the latter case thus the adiabatic curves are given by vΘ 3 = const.
Considering that we find ∂x eq γ ∂v x s = 0 (A49) and, therefore, This result is in agreement with the well-known fact that in ultra-relativistic ideal gases the bulk viscosity is identically zero [2,7].
Let us focus on the case w = 3/2 (non-relativistic matter-gas). Starting from the obvious relation we obtain, recalling (A45), ∂Θ ∂v x s = − 2Θ 3v This can be used to show that ∂x eq γ ∂v x s which plugged into (142), gives It is useful to rewrite this formula using more standard notation. To do this, we introduce the pressure ratio where the second identity follows from the ideal gas assumption P n = nΘ. Plugging it into (A52) we recover the well-known formula [34] ∂ log Θ ∂ log v x s = − 1 + 4α P 3/2 + 12α P . (A56) Finally, our expression for the bulk viscosity becomes We see that the second fraction tends to suppress ζ as α P −→ 0 or α P −→ +∞. Intuitively, this is due to the fact that since the bulk viscosity is due to the dissipative processes which tend to equilibrate the temperatures of matter and radiation (which would depart from each other during a fast expansion), it becomes important only if both the species give a relevant contribution to the overall stress-energy tensor.