A multifluid perspective on multimessenger modelling

This brief review introduces the notion of a relativistic multifluid system -- a multi-component system with identifiable relative flows -- and outlines a set of models for scenarios relevant for different astronomical observation channels. The specific problems used to illustrate the key principles include superfluid hydrodynamics (with relevance for radio and x-ray pulsar timing and gravitational-wave searches), heat flow (connecting to the problem of neutron star cooling and associated x-ray observations) and the coupling between matter and electromagnetism (linking to explosive phenomena like gamma-ray bursts and more subtle issues like the long-term evolution of a neutron star's magnetic field). We also comment on the coupling between matter and radiation, for which the multifluid approach would seem less appropriate. The main motivation of the survey is to illustrate less familiar aspects that come into play in multifluid problems, establish the relevant"language"and provide a platform for more detailed work on these issues.


II. MULTIFLUID MODELS
The convective variational approach pioneered by Carter [14] and a number of collaborators (see [13] for a detailed review) provides a natural framework for a discussion of relativistic multifluid systems. The approach is fairly intuitive. Consider a system with a number of identifiable matter constituents, labelled by x, y, . . ., each associated with a flux (we use a, b, c, . . . to represent spacetime indices and it should be noted that the label x is not summed over when repeated). Each four velocity is normalised, g ab u a x u b x = −1 (with g ab the spacetime metric and we assume units such that c = G = 1 throughout the discussion), and the number density measured by a co-moving observer is The variational principle involves a matter Lagrangian, Λ, taken to be a relativistic invariant and hence depending only on covariant scalars formed out of the different fluxes. Effectively, this provides the equation of state for matter. In the general case, we have to consider both and n 2 xy = −g ab n a x n b y with y = x .
An unconstrained variation of Λ with respect to the fluxes n a x and the spacetime metric then gives (ignoring terms that can be written as total derivatives, which may be thought of as "surface terms" in the action) where g is the determinant of the metric. This defines the individual fluid momenta with and The covector µ x a provides the fluid momentum and its magnitude gives the chemical potential (intuitively, the energy associated with adding a particle of the x-species to the system). For a co-moving observer, we have We can now identify one of the key features of a multifluid model. Each momentum may be "tilted" relative to the corresponding matter current. This is the so-called entrainment effect [13,15,16], and it enters the model through the A xy coefficients. We will consider the relevance of this later.
Returning to the variational argument, equation (5) shows why it is necessary to constrain the action. As it stands, the variation of Λ suggests that the equations of motion would be µ x a = 0-none of the fluids carry energy or momentum, which is not a very interesting situation to consider. The actual equations of motion follow if we constrain the fluxes to be conserved [80]: ∇ a n a x = 0 .
This is typically done by changing to a Lagrangian perspective and making use of a lower-dimensional matter space [14]. The relevant steps are outlined in [13] and we will not repeat them here. Our focus will be on applying the model rather than the derivation. The constrained variation leads to where ξ a x represents the Lagrangian displacement associated with each fluid's set of worldlines, which are taken to be conserved as a consequence of (10). We have also introduced the "force" acting on the x-component, f a x , leading to the equations of motion (since the displacements are independent and arbitrary). It also follows that the stress-energy tensor (the variation with respect to the spacetime metric) takes the form where is a generalised pressure. We will make this notion more precise in the following. The equations we have written down may seem somewhat peculiar. The standard approach to relativistic fluid dynamics takes the divergence of the stress-energy tensor-and the fact that it must vanish by virtue of the Einstein equations and the Bianchi identities-as the starting point. Here we seem to instead have a set of force-balance relations represented by (12), involving the vorticity formed from each fluid momentum However, the two pictures are perfectly consistent. We have In other words, if the set of equations (12) are satisfied then ∇ a T a b = 0 follows as an identify. Intuitively, the argument in (16) is "Newton's third law" in action. Let us keep the focus on the stress-energy tensor (13). It may be helpful to, first of all, highlight the single-fluid result (simply dropping the chemical labels as there is only one matter component); In this case we have n a = nu a and µ a = µu a so the expression for the generalised pressure leads to We then have and an observer moving along with u a would measure the energy This then leads to and a comparison with the standard Gibbs relation from thermodynamics supports the interpretation of Ψ as the pressure, p, in this case. This then leads to an interpretation of (14) as a generalised Gibbs relation, appropriate for the multifluid setting and analogous to (more phenomenological) relations introduced in the context of extended irreversible thermodynamics [19].

III. LINEAR DRIFT ARGUMENT
At this point it should be clear that the multifluid model extends the familiar fluid logic and introduces features that we need to get used to. On the one hand, this is important as it allows us to describe the additional dynamical degrees of freedom of a given problem-and the coupling between them. On the other hand, the model may be too general for many situations of interest. In particular, we have to keep track of the individual fluxes, or-if we introduce a preferred observer-the different relative velocities and the associated Lorentz factors. We may also not have a good handle on the parameters we need to complete the model. In particular, the nuclear physics calculations that provide state-of-the-art neutron star equations of state tend to assume that the matter is in equilibrium. There are two aspects to this. First, we may insist on chemical and thermodynamical equilibrium, as appropriate if we want to determine a static matter configuration. If we want to model a dynamical scenario-even for a single matter component-we need to allow for the matter not being in equilibrium, e.g. due to changes in the local temperature. This is regularly done in modern numerical relativity simulations. The multifluid problem adds a different dimension, as the relative motion between the different components means that there is no "obvious" matter frame with respect to which we can develop the thermodynamics [81].
In order to explain this point, and as a step towards a simplified model, let us introduce an observer, associated with a four velocity U a such that (with the speed of light c = 1) with U a V a x = 0 and where provides the Lorentz factor associated with the different flows in the system. We then have and it is natural to introduce the number density measured by U a aŝ For the momenta we get with the relative velocity andμ x = −U a µ x a the chemical potential measured by the chosen observer. The energy measured by the observer is (recall (17)) This seems like the intuitive generalisation of (21), but if we combine the result with (14) we see that we can no longer identify −Λ with the measured energy. Instead, we need which leads tô The quadratic terms complicate the thermodynamics (essentially representing kinetic energy contributions relative to the chosen frame). In order to make progress we need to keep track of the entrainment and the relative velocities. From a formal point of view, we have everything we need to solve the problem. In particular, the prescription is Lorentz invariant and as such it is expected to prevent causality violations. However, if we want to make contact with a given model for the microphysics then we have a problem. The relative flows (obviously) impact on the energy of the system and this, in turn, affects the thermodynamics and issues like chemical equilibrium. If the chemical potentials are observer dependent-as they have to be in the general case-then we need to decide how to implement the notion of equilibrium [18]. Who measures what? This problem does not arise in standard approaches to dissipative fluid dynamics, which tend to be based on formal expansions in the deviation from equilibrium [21][22][23] (see [24][25][26] for recent progress). Our problem is different: We need to reconcile the nonlinear multifluid model with our "incomplete" description of the thermodynamics [17,27,28]. As a step towards applications, we may simplify the model by linearising in the relative flows-essentially neglecting quadratic terms in V a x [11]. This may seem like a drastic move, but it is a natural assumption because many of the situations we are interested in involve a gentle and gradual drift of one component relative to another rather than a vigorous relativistic stream. Of course, we have to be careful, because if we linearise the model we ignore the relative Lorentz factors, and we have In effect, the model is no longer Lorentz invariant and we must not extend it into the nonlinear regime. However, if we take appropriate care in this respect, the linear drift model has a number of attractive features. In particular, it is easy to see that the connection with thermodynamics becomes straightforward. As γ x ≈ 1, all linearly related observers agree on the number densities and chemical potentials. We also getε ≈ ε ≈ −Λ (in essence, we can drop the hats on all the scalar quantities) and the usual Gibbs relation applies. This is quite attractive. Given the linear drift assumption, the equations of motion (12) simplify to (after a bit of algebra and noting that the force only has three independent components-recall that n a x f x a = 0 follows from (12))-so we focus on the projection orthogonal to U a ) where For the stress-energy tensor we get (noting that the entrainment contributions cancel!) where we have used the usual multi-component Gibbs relation to define the pressure The result now looks much more familiar. If we define the relative momentum flux we get the final expression We recognise this as the stress-energy tensor for a system with a linearised heat flux [13]. Notably, the "stress" terms associated with relative flows are removed by the linearisation. In order to complete the description, we need the continuity equations. It is easy to see that these become With this we have all the equations we need to consider specific applications, notably without at this point having committed ourselves to a specific observer. We could carry on working with a general model, but this might get confusing so it is better to consider specific examples. We will do this in three steps. First we consider a superfluid system with two components, then we turn to a problem with heat flux and a non-conserved entropy current and finally we outline a model for charged systems with resistivity. These problems have all been discussed in the literature (see [11] for a closely related discussion), so we do not expect to learn anything "new" here. Rather, the point of the discussion is to highlight the common features of the different models and consider the framework required for more complex situations step by step.

IV. SUPERFLUIDS
The variational multifluid model provides a natural framework for discussing superfluid systems, both in the laboratory context and for neutron stars. The minimum requirement of a workable model is that it should represent the expected two-fluid dynamics (e.g. the second sound and the presence of oscillation modes that reflect the two degrees of freedom, see [29] for a recent discussion) and the anticipated entrainment effect (which encodes the relative "mobility" of the two components). As we will see, the variational model easily satisfies these criteria [82]. Moreover, we can go further and consider the quantisation of vorticity-required for a superfluid to exhibit bulk rotation-which is easily imposed on the canonical momentum [31] and the implications for the dynamics become quite intuitive. Alternative descriptions are, of course, available (see for example the seminal contributions by Israel [33][34][35] or the more recent discussion in [36] for the relativistic case or [29] for a recent summary of the corresponding Newtonian problem). We will only comment on a few of the relevant aspects here.
First of all, note that (12) is automatically satisfied if the momentum is given by the gradient of a scalar phase, µ x a = ∇φ x . This implies that the fluid is irrotational [83]. However, for most problems of interest one assumes an averaged vorticity such that the normal fluid equations apply. The argument comes with the caveat that the required (quantised) vortices may introduce new features (like mutual friction and vortex elasticity [37,38]), but here we will focus on the simplest model for bulk dynamics. The entrainment then plays a central role, so it is natural to start by providing an interpretation of this effect.
In order to make the example clear, let us focus on the two-fluid model for neutron star cores, with a neutron superfluid (x = n) coupled to a charge-neutral conglomerate of protons and electrons (x = p). We then have the spatial components of the proton momentum (relative to U a ) Suppose we let the observer ride along with the neutrons, in such a way that V i n = 0, then we can introduce the effective proton mass as (this argument is analogous to the Newtonian discussion in [29], see also [39]). A similar argument for the neutrons defines m * n and we see that we have In a neutron star core, the entrainment accounts for an important aspect of the strong interaction; the flow of neutrons pulls some protons along (and vice versa) [16]. The notion of entrainment is also important for the dynamics of the star's inner crust, in this case due to Bragg scattering off of the elastic lattice of nuclei [40,41]. However, the presence of the nuclear lattice brings in additional aspects which we will not go into here, see [29] for a brief summary. Also, in the general case, with a number of distinct flows (or several entrainment mechanisms), the expression for the effective mass is not as simple as (40)-we end up with an effective mass matrix-but the concept still makes sense.
Turning to the dynamics of the system, an advantage of considering superfluids is that both nuclear reactions and resistive scattering are suppressed, so we can focus on the non-dissipative problem with conserved fluxes. This is precisely the setting for which the variational approach was intended [13], so we can make direct use of the equations we have written down. Moreover, in many situations of interest the relative flow between the components is small so the linear drift assumption should apply.
When we consider the equations of motion, we need to make further decisions. We can introduce a specific observer and we may opt to work with different combinations of the equations. It is worth explaining both of these points as different strategies have been used in the literature and it is helpful to understand how they are related.
First consider the momentum equations (32). We can decide to represent the two dynamical degrees of freedom by the equations for f n a and f p a . This makes the mathematics fairly "symmetric" in that the two fluids are treated on the same footing. As an example, this strategy has been used in the context of neutron star seismology [30,42,43]. However, this approach does not exactly reflect the anticipated dynamics. It is easy to argue that the natural degrees of freedom are better described in terms of the total flux and the difference between the velocities (see [29] for a discussion and pointers to the literature). The decoupling is not complete (apart from in idealised situations) but we can opt to work with a set of equations that is closer to the physics. This could, for example, involve the sum of the two force equations-which we already know is equivalent to the vanishing divergence of the stress-energy tensor-and a weighted difference (such that the four acceleration term in (32) is removed). An example of this procedure, in the context of charged plasmas, can be found in [44]. A somewhat simpler option is to combine the equations from the stress-energy tensor with one of the fluid force equations to remove the four acceleration from the latter to get a direct expression for the evolution of the chosen velocity V a x . We will outline this approach, which is close in spirit to the more "traditional" approach to superfluid dynamics-notably championed by Gusakov and collaborators in the context of relativity [45,46]-in the following.
Before we work out the relevant details, it is useful to consider the choice of observer. First we note that, if we add the continuity equations (38) and define then we get We can clearly simplify this by letting the observer be such that V a = 0. This has the advantage that (44) reduces to the usual single-fluid continuity equation, involving only n and U a . In the two-component example we are considering at the moment, n = n n + n p represents the total baryon number density and we have This relation allows us to work with equations for U a and V a n . The specific observer choice is analogous to the Eckart frame familiar from discussions of dissipative relativistic fluid dynamics [47].
Another option is to focus on the stress-energy tensor and simplify life by making it as close to the perfect fluid as we can. This involves introducing an observer such that This choice leads to the well-known Landau-Lifshitz frame [48]. Other combinations are, of course, available. The bottom line is that we can use the observer to "simplify" some aspects of the models (possibly at the cost of complicating others...). It is also worth keeping in mind that the remaining relative velocity, V a n , refers to different observers and therefore has slightly different "meaning" in the two cases.
Suppose we commit to one of these choices, say, (46). Then we have the text-book equations of motion for a perfect fluid (as the stress-energy tensor takes the same form): and This is nice, but... we are still dealing with a two-fluid problem. Representing the second degree of freedom by the superfluid neutrons, we also need Using (48) to remove the four acceleration from this equation, we have (admittedly after some algebra...) an evolution equation for the superfluid component Recalling the effective mass from (41), we have which illustrates how entrainment impacts on the inertia of the superfluid.
It is now clear that the observer choice simplifying (48) is a mirage. We still have to deal with the second degree of freedom and the relevant momentum equation has the additional features we expect, like entrainment. The fact that the problem remains complex is also apparent when we consider the continuity equations. We can take one of the two relations to be represented by (47) (representing energy conservation), which is nice, but we need to keep in mind that the model requires us to keep track of two number densities, some combination of n, n n and n p or some function of them, like ε. In the single fluid case, we know that so (47) can be written and the conserved flux implies the energy equation (and the other way around). This does not follow in the two-fluid case where-at the linear drift level-we have ∇ a ε = µ n ∇ a n n + µ p ∇ a n p , and it follows that, if we choose to work with n and n n (say) we get This form is instructive because the second term vanishes if we impose chemical equilibrium [84], bringing us back to the single fluid result. Of course, in the case of superfluids the nuclear reactions required to establish equilibrium are suppressed so this argument does not really help us. We need to keep track of the second number density, which (in this example) is governed by (U a + V a n ) ∇ a n n + n n ∇ a (U a + V a n ) = 0 .
Although there are many issues one might want to add to the discussion, especially regarding vortices [31,38,49,50], will move on at this point. The key aim was to illustrate how the multifluid model can be adapted to account for the superfluid degree of freedom and how this introduces new microphysics parameters associated with entrainment (and which have to be provided from nuclear physics [16], see [51,52] for recent efforts in this direction). The simple fact that the two-fluid equations that result involve elements of choice is important to keep in mind.

V. THE FLOW OF HEAT
Staying at the formal level, let us turn to a problem that allows us to discuss how the model changes when the fluxes are not conserved. The simplest reasonable problem in this direction is that of heat flowing relative to matter [53][54][55][56][57]. In this case, the minimal model includes a distinct heat flow and the connection to the second law of thermodynamics (entropy must not decrease) which makes the problem dissipative. This situation can (again) be represented by two components. We will take them to be the matter (x = n) and the entropy (s) and it would usually be the case that the associated heat flux is a gentle drift so we should be able to make progress within the linear approximation.
A key aspect of the heat-flux problem is that the entropy current is not conserved. Letting s a = n a s we have (for an isolated system) in accordance with the second law of thermodynamics [85]. Given this, let us first ask a general question. What happens to the multifluid formalism if the fluxes are not conserved? This forces us to either adjust the variational argument (as explored in [17,18,28]) or proceed in a more phenomenological manner. Here, we will take the latter approach, essentially adding the main dissipation mechanism of interest and not worrying too much about the general problem. We simply note that if we have in (16), then we must adjust the force equation in such a way that [17,27] 2n and the constraint That is, the reaction rate determines the time component (in a co-moving frame) of the resistivity R x a (which, in turn, must be included in (59) since the original force term is orthogonal to n a x ). Let us now focus on a problem involving a single matter component-letting n a n = n a , for simplicity, with the corresponding chemical potential being µ n = µ-and an entropy component (from before) for which the chemical potential is the temperature, so µ s = T . At the linear drift level, we then have the (familiar) thermodynamical relation It also follows that (57) becomes where we have introduced the heat flux (in the usual way) Turning to the momentum of the thermal component, when we allow for entrainment between matter and heat-this may seem like an unusual idea, but let us go along with it and see where it takes us-we have and we may introduce the effective mass for the thermal component through This relates the entropy entrainment to the effective inertia of heat, a notion that is expected to play an important role in the relativistic heat problem [58,59]. So far, the results are straight translations of previous relations. Moreover, opting-as in the superfluid case-to work in the centre of momentum frame, we have (assuming that sT ≪ nµ, which seems likely to be the case in most situations of interest). This means that we have In order to proceed, we need an expression for the resistivity, R a x . Making use of the results from [28] we haveassuming that the reaction rates Γ x and the resistivity coefficients R xy are provided by the microphysics, for all material particles, along with the constraint on the entropy: In the present context, with only two components and Γ n = 0, this leads to We see that the R ns = R sn coefficient is required to be positive by the second law of thermodynamics. It is also useful to make contact with the standard discussion of the problem by introducing the thermal conductivity κ such that It follows that-as should have been expected-the entropy creation rate Γ s is quadratic in the relative velocity. Hence, the corresponding contribution to the momentum equations should be neglected at the linear drift level. This leaves only the resistive contribution to the entropy momentum equation. In effect, we only need to add a term to the right-hand side of (32) (where we need to keep in mind that we divided through by an overall factor of n x ) This means that we have an evolution equation for the heat flux: This model-which arises naturally in the multifluid framework-helps explain the main features of the relativistic heat problem. In the absence of heat flow the terms on the right-hand side bring out the classic result (when combined with the Tolman-Oppenheimer-Volkoff equations for hydrostatic equilibrium) that the redshifted temperature is uniform [60]. In effect, heat has weight. If we neglect the terms involving U a , we get the relation for the heat flux q a required for a derivation of the equations used to describe the gradual cooling of neutron stars [13]. In the general case, we have a telegraph-type equation-as in extended irreversible thermodynamics [19]-with the effective mass m * s regulating the thermal relaxation, which is key to keeping the model causal and stable [57]. Perhaps the main lesson of the exercise is that it paid off to take the entropy entrainment seriously.

VI. ELECTROMAGNETISM
The third problem we are going to consider, again from the two-fluid point of view, introduces charged flows and the coupling to electromagnetism. The marriage between electromagnetism and multifluid model is easy and comfortable [44], as should be expected given that both derive from an action principle [13]. When we consider the problem of charged flows in relativity, with the aim to reach beyond text-book results, a minimal requirement is to add resistivity to magnetohydrodynamics [61][62][63][64]. Essentially, we want to write down a version of Ohm's law that connects to the notion of the charge current as arising from the relative flow between two components [65,66].
Starting from the usual assumption of minimal coupling, the electromagnetic Lagrangian is built from the antisymmetric Faraday tensor; where A a is the vector potential, and the electromagnetic field couples to the matter flow through the charge current j a (via a coupling term µ 0 j a A a , where the coupling constant µ 0 is the magnetic permeability). In order for this construction to be gauge invariant, the current must be conserved. That is, if the model involves a number of charge carriers, each with a charge q x per particle, we have the constraint where Clearly, this constraint is automatically satisfied when the individual fluxes are conserved. If we account for reactions, as we may want to do for neutron stars, we must impose overall charge conservation The coupling to the vector potential impacts on the matter momentum in the fashion anticipated from text-book quantum mechanism, and we haveμ Formally, this is the only change we need to make to the fluid equations. As long as we change µ x a →μ x a the individual force equations (12) remain unchanged. Meanwhile, a variation of the action with respect to the vector potential (keeping j a fixed!), leads to the Maxwell equations completed by the symmetry relation Finally, a variation with respect to the spacetime metric leads to the electromagnetic contribution to the stress-energy tensor, which satisfies in turn, defining the Lorentz force, f b L . Basically, we may take the view that the matter stress-energy tensor still takes the form from (13) (in terms of the pure matter momenta) and that we have [86] ∇ a T ab If we insist on separating the electromagnetic contributions from the matter ones in this way, we end up with a new set of momentum equations of form Let us unpick the different contributions to (82) and (84), starting with the electromagnetic terms. In order to help intuition, it may be useful to introduce the (obviously observer dependent) electric and magnetic fields, e a and b a (although it is obvious from (84) that, if Γ x = 0 we must also keep track of the vector potential, see [67] for the starting point for such a formulation). These are defined by where ǫ abc = U d ǫ dabc is associated with a right-handed tetrad moving along with the observer (and noting that each field has three independent components as U a e a = U a b a = 0). We also have the charge current That is, we have and it is evident that-since they agree on the number densities-all observers within the linear drift family must agree on the charge density σ. This is important as it means that we can consistently impose the usual condition of charge neutrality [87] simply by setting σ = 0. This then leads to while the Lorentz force takes the form Meanwhile, for the individual momentum equations (84) we need Finally, we need the resistivity R x a , which is given by (69) In order to make the discussion more specific, let us focus on the simplest neutron star setting [88], i.e. a fourcomponent system, with neutrons (n), protons (p), electrons (e) and entropy (s). This system reduces to two fluid degrees of freedom if we allow the electrons to drift relative to the other components (which are taken to be locked [89]), thus providing the required charge current. In effect, if we take V a n = V a p = V a s = V a e then the charge current is given by (assuming local charge neutrality and using q e = −e, which should not be confused with the local electric field e a ) We also have the condition for the centre of momentum frame: where local charge neutrality leads to n p = n e . Putting things together, we have where the approximation is valid as long as n e µ e ≪ p + ε. We will assume that this is the case [90]. We now have all the ingredients we need to write down the momentum equation for the electrons (for which we ignore the possible entrainment coupling [91], leading to A ex = 0). Adapting the previous results, we get Assuming that the electron flux is conserved (Γ e = 0), defining the total resistivity and introducing the electro-chemical field [11,68] we have In general, this dynamical equation for the charge current complements the equation for total energy-momentum conservation. If the latter is taken from (83) and we work in the centre of momentum frame (as we have assumed), then T ab M notably has the perfect fluid form. This supports the logic often taken as starting point for relativistic magnetohydrodynamics in the literature, but it is evident that this assumption strictly depends on the choice of observer (or may only be approximately true [11]).
At this point, the problem still has two distinct degrees of freedom, represented by U a and J a , both governed by dynamical evolution equations. In essence, the model retains the underlying plasma properties. To simplify things, one would typically start by ignoring the right-hand side of (98) (essentially the electron inertia). This leads to a version of Ohm's law: If we further ignore the "battery terms" (the gradients in (97)) and the Hall term in (99) we are left with where we have definedR = R/e 2 n 2 e . This result should be familiar. It is obviously rewarding to see that we arrive at the expected result from the multifluid model, but the most important insights relate to the steps involved in the derivation. Depending on the degree of accuracy we require, we can undo the steps one by one and check if the assumptions are warranted for different problems of interest.
Finally, in the ideal limit-when the resistivity vanishes so we have a perfect conductor (withR = 0)-we see that the electric field vanishes in the fluid frame. This is the usual assumption of magnetohydrodynamics. Of course, we have to be mindful of the chosen frame and the assumption that the other contributions to (98) can all be neglected. Again, the usual result is most likely only approximately true.

VII. RADIATION HYDRODYNAMICS
So far we have considered three models relating to a range of observed neutron star phenomena: pulsar glitches explored through radio timing and considered as possible gravitational-wave sources (superfluid hydrodynamics), neutron star cooling observed by x-ray satellites (heat flux) and a range of magnetic field related phenomena, from explosive dynamics to the gradual long-term field evolution over thousands of years (resistive magnetohydrodynamics). These examples demonstrate the versatility of the multifluid formalism, but (at least) one interesting problem remains to be contemplated: the coupling between matter and radiation. This is a crucial issue for both supernova core collapse (where neutrinos are thought to provide the key explosion mechanism) and neutron star mergers (where the merger remnant heats up to temperatures similar to those reached in high-energy colliders). The problem of neutrino transport is extremely challenging [69], but it has one simple and intuitive limit. At high temperatures, the neutrinos are trapped by the matter and they may have a short enough mean-free path that we can meaningfully describe them as a fluid [70,71]. Without going into the fine print of the problem, it seems useful to ask to what extent we can expect to make progress with such a strategy. In essence, we want to consider if the radiation problem can be meaningfully considered from the multifluid perspective.
At a glance, the idea seems promising. Suppose we consider the simple case of a single photon-the logic is similar for neutrinos, but in that case we also have to consider lepton number conservation-then we have a radiation stress-energy tensor where h is Planck's constant, ν is the photon frequency and is null, so l a l a = 1. The vector l a provides the direction of propagation relative to the observer U a . It is easy to see that this leads to the anticipated energy We can obviously rewrite the stress-energy tensor as with and Viewing these results from a fluid perspective, it is clear that (104) reminds us of, for example, (37)-after all, it has the form of a general stress-energy tensor. The comparison is, however, misleading. For example, in the case of radiation, we do not have-at least not yet-an identifiable "drift velocity" that can be linearised. A photon moves at the speed of light regardless of the observer. There may be a meaningful concept of "collective drift", but this is not evident yet. Formally, we can always assume minimal coupling between matter (a perfect fluid, say) and radiation in such a way that This is obvious, as is the fact that we run into trouble if we try to bring the variational multifluid formalism to bear on the problem. The variational approach builds on the conservation of fluid element worldlines [13]. Because the radiation is null, this argument does not apply (in the free-streaming limit). Of course, the logic should work for neutrinos (which have a small-but nevertheless-rest mass) and one might also make progress by thinking of the radiation as associated with an effective mass (as we did in the case of heat, where the entropy is also massless, strictly speaking). Intuitively, the argument should work for systems with "trapped" radiation, so the objection may be more a formality than a real stumbling block. One might, for example, consider a picture where the interaction with matter (perhaps through some form of entrainment?) endows the radiation with an effective mass. This effective mass would then have to vanish in the free streaming limit (which inevitably becomes singular given that k a is null, but one might be able to make this workable). The real question is not if we can build a multifluid model for the radiation problem. Rather, we might want to ask if we should. In order for the effort to make sense, we would have to learn something new by going in this direction. This is where we seem to run into trouble. In a realistic setting, we are not dealing with a single photon with a well defined direction of propagation. Instead we have an energy/momentum spectrum, with a distribution typically governed by Boltzmann's equation. The formalism for this is well developed and the issues associated with it are well understood (in the context of kinetic theory). The actual challenge is computational [69,72], given the added dimensions associated with the radiation phase space. This is why practical implementations [73][74][75][76] typically involve simplifying the problem by integrating out momentum aspects, leading to a well-defined moment expansion. Truncating this expansion, we arrive at a radiation stress-energy tensor of form (104), which means that the formal fluid comparison "makes sense" [77] but in reality we do not learn much from this. The key point is that we do not have a single relative velocity that can be used to identify the "second fluid" in the problem. Instead, we have a situation where different pieces (say, frequency bins) are associated with separate momentum contributions, leading to the multi-group approach to numerical simulations. We could always think of this as a souped up multifluid problem, with different "fluids" representing different frequency bins, but there hardly seems to be any merit to this idea. It certainly does not provide us with anything that the current technology does not already consider. In essence then, the coupled problem matter and radiation does not appear to be a multifluid problem, at least not in a meaningful sense. There may be a strongly coupled limit where the idea applies, but this seems somewhat artificial (as we cannot model the transition to the free-streaming limit) and it does not really help us make progress on the problem we actually need to solve.

VIII. CONCLUDING REMARKS
This survey was written with the aim of providing a (somewhat introductory) perspective on the ongoing effort to develop a flexible formalism for relativistic multifluid systems. The connection between this (more formal) endeavour and the need for a reliable description of the physics required for multimessenger astronomy is fairly clear. A number of relevant physics scenarios require multifluid aspects to be taken into account. As illustrations, we considered three settings: i) the problem of superfluid dynamics, ii) the issue of heat flow, and iii) the combination of charge currents and resistivity leading to Ohm's law. In each case we demonstrated that the multifluid approach provides valuable insight, motivating concepts (like the effective inertia of heat) which have often been introduced phenomenologically. The simple fact that the different problems are described within a single over-arching framework is important as it makes a combination of the problems fairly straightforward [11]. For example, one can easily formulate a three-fluid model to demonstrate the thermo-electric effect or consider how the onset of superfluidity impacts on the evolution of a neutron star's magnetic field. These-and many other-interesting problems are now within reach. This is rewarding, but the job is not done. Neutron stars are complex systems, and the multifluid model is not (quite) able to cover all the relevant aspects. Most notably, many scenarios of interest involve dissipation and we also need to improve the description of the coupling between matter and radiation (like photons or neutrinos). In this survey, we only tried to explain why the latter is not a multifluid problem (at least not in the usual sense). The argument is fairly obvious and may seem rather dismissive. This was, however, not quite the intention. The aim was to argue that the fluid approach to radiation will always be limited and to suggest that the multifluid logic is only truly useful if it brings new perspective. The three models we discussed provide such insights while the consideration of the radiation problem does not (yet) do so. There is, of course, scope for progress in this direction. We may just have to dig deeper.
While the variational framework for multifluid systems is mature and covers a number of aspects with (more or less) immediate importance for astrophysical applications, a number of issues require further thought. This involves both formal aspects and applications. For example, on the formal side, it would be interesting to extend existing work on elasticity (see [13] for a summary) to account for plastic flow and viscoelasticity and consider the impact of the superfluid neutron component that is present in the inner crust of a neutron star. This formal development then has an immediate application, as one may consider the impact of the gradual superfluid decoupling (as the star cools) on the evolution of a neutron star's magnetic field. This would be naturally phrased as a three-fluid problem. Of course, we also need to consider what is calculable and what is not (given our understanding of the physics). In particular, we need to make sure that the models remain grounded in realistic microphysics. For any problem involving neutron star crust dynamics, the entrainment [40,41] is expected to play an important role but a number of related issues remain (somewhat) unsettled (see [51,52]). The context of elasticity is, of course, just one of many possible future directions. The framework we have outlined should be flexible enough to cope with pretty much anything we choose to throw at it.

Funding
The work described in this article was supported by STFC through grant number ST/R00045X/1.