Hadronic transport coefficients from the linear sigma model at finite temperature

We investigate general frameworks for calculating transport coefficients for quasiparticle theories at finite temperature. Hadronic transport coefficients are then computed using the linear sigma model (LSM). The bulk viscosity over entropy density ($\zeta/s$) is evaluated in the relaxation time approximation (RTA) and the specific shear viscosity ($\eta/s$) and static electrical conductivity ($\sigma_{el}/T$) are both obtained in the RTA and using a functional variational approach. Results are shown for different values of the scalar-isoscalar hadron vacuum mass with in-medium masses for the interacting fields. The advantages and limitations of the LSM for studies of strongly interacting matter out of equilibrium are discussed and results are compared with others in the literature.


I. INTRODUCTION
The behaviour of strongly interacting matter in extreme conditions of temperature and density is the subject of a vibrant experimental program and of numerous theoretical efforts. One of the major achievements of relativistic heavy-ion physics is the realization that an exotic phase of nuclear matter -the Quark-Gluon Plasma (QGP) [1] -has been created in experiments performed at RHIC (the Relativistic Heavy Ion Collider, at Brookhaven National Laboratory) and at the LHC (the Large Hadron Collider, at CERN).
A related collection of theoretical breakthroughs has shown that the dynamical evolution of this QGP is amenable to hydrodynamic modeling [2]. Hydrodynamics is able to interpret a large body of data that reflects the collectivity of the observed particles and can even make quantitative statements about local deviations from equilibrium. The response of a quantum system to some perturbation can be characterized in different ways, including by monitoring the time it takes for the system to relax back to the equilibrium state. The relaxation time can then be related to transport parameters, which are calculable in terms of correlation functions [3]. In spite of the existence of a general formalism to calculate transport parameters, obtaining them from QCD has remained challenging. This has motivated their extraction from analyses of heavy-ion phenomenology [2,[4][5][6] and from effective models of the strong interaction.
This work reports on studies of transport coefficients using the linear sigma model (LSM). The LSM, first proposed by Gell-Mann and Lévy [7], is a well-studied model of a simple hadronic system and is one of the most instructive and paradigmatic field theories. Hadronic physics estimates have yet to agree on details of the transport parameters that characterize the strong interaction, whether those come from analyzing and interpreting data * heffernan@physics.mcgill.ca † jeon@physics.mcgill.ca ‡ gale@physics.mcgill.ca or from attempts to calculate from more-or-less first principles. Therefore, it is instructive to make predictions using the LSM. In addition, there is a need to distinguish the effects of approximations from the consequences of model-dependent assumptions. This is one of the goals pursued herein. We adapt and develop various general theoretical techniques in order to calculate the transport coefficients of the LSM. The shear and bulk viscosity are computed and the electrical conductivity is calculated in the LSM for the first time. We discuss the effect on our results of the assumed value of the vacuum mass of the σ meson using figures that broadly span the mass of the f 0 (500) as defined by its large reported width [8]. The paper is organized as follows: Sec. II reviews and details the theoretical framework, established along the lines of work done in Refs. [9,10]. It begins with a survey of the LSM with some emphasis on the thermodynamic quantities that are used as inputs into both the relaxation time approximation and the variational method. Sec. III includes a brief overview of transport and shows the technique for calculating transport coefficients in the relaxation time approximation. To go beyond the limitations of the RTA, a general variational technique for massive theories with elastic and inelastic reactions is developed in Sec. IV. Finally, in Sec. V, we produce the results of the calculations with the different techniques and compare to the literature. Additional details of calculations are provided in the appendices.

II. THEORETICAL FRAMEWORK
We begin the development of the theoretical framework with an overview of the linear sigma model. Once this is established, we describe the treatment of the effective masses and the thermodynamics of the system.

A. The linear sigma model
The LSM is a ubiquitous relativistic field theory capable of illustrating aspects of low-energy QCD [11]. It can incorporate mean field effects and is thermodynamically arXiv:2005.12793v1 [hep-ph] 26 May 2020 consistent. It has been used extensively as an effective model of simple hadron dynamics and is tractable and well-studied. Consequently, we use it to provide parametric estimates for transport coefficients of hadronic systems, to demonstrate the effectiveness of our framework, and to highlight the quantitative effects of different approximation schemes.
The classic linear sigma model Lagrangian is In general, the bosonic field Φ has N components. When N = 4, the standard practice is to ascribe the first N − 1 components to the pion field and define the N th as the sigma field: Φ = { π, σ}. Then, the LSM model is an effective theory of soft pion dynamics owing to the isomorphism between O(4) -the symmetry of the LSM -and SU(2) L × SU(2) R -the symmetry group for two flavours of massless quarks in QCD. At low temperatures, the O(N ) symmetry is spontaneously broken to O(N −1) and a temperature-dependent σ condensate v appears: In the classic LSM, the condensate goes to zero at a critical temperature in a second-order phase transition [12]. The LSM can be made more realistic by explicitly breaking chiral symmetry with a pion vacuum mass. The Lagrangian can now be written as [13,14] where the vacuum expectation value v of the scalar field σ is determined by the symmetry-breaking term The three undetermined parameters λ, H, and f 2 are determined by the vacuum values of the pion decay constant f π and the pion and sigma masses In this subsection, all mass symbols represent vacuum masses -the reason for this clarification will soon become apparent. The vacuum pion mass is chosen to be m π = 140 MeV, the decay constant is f π = 93 MeV, and in this work the vacuum sigma mass will take one of the values m σ = {400, 600, 900} MeV. The nature of the symmetry breaking at zero temperature drastically impacts how the chiral symmetry is restored at high temperatures. We separate the LSM Lagrangian into kinetic and potential terms and expand the sigma field into a condensate and an excitation, σ → v + σ. We call the excitation σ as it is the true σ meson and the condensate v, which is the non-vanishing vacuum expectation value of the field. We perform our calculations in the isospin pion basis, which represents the physical pions. The relations between the physical pions and the Cartesian pion fields are To convert the Lagrangian to the isospin pion basis, it is simple to invert these relations. Doing so allows one to trivially rewrite the Lagrangian and read off the matrix elements. For example, in the calculation of M π a π b ;π c π d , we see that from the 4-point diagram, we get a factor −2λ. From π + π → σ → π + π in the s-channel, we get a factor 4λ 2 v 2 s−m 2 σ . Thus, the 4-point pion s-channel diagram is Including other processes with appropriate Kronecker deltas produces the full matrix element, A pole clearly arises in each channel. A consistent method of handling this pole theoretically would entail a resummation that would parametrically promote the process to higher powers of λ. However, the coupling in the LSM is larger than one and the diagrammatic expansion is not convergent. We treat the LSM as an effective theory with λ understood as a parameter adjusted to fit π-π scattering [12]. We therefore restrict the kinematics in order to bypass the singularities and the processes to tree-level. We use the Mandelstam variables in the limit s, t, u → ∞, effectively removing the 3-point interactions [9]. This results in the following matrix elements for the LSM: As we wish to study the critical dynamics at finite mass and temperature we adopt an approximate lower bound of T = 150 MeV: a common value for thermal freeze-out in studies of heavy ion collisions, which is also in the vicinity of the crossover temperature obtained in lattice calculations with baryonless QCD [15]. Higher temperatures will also be explored so that the behavior of physical quantities of interest -such as the transport coefficients -can be studied through gradual chiral symmetry restoration in the LSM.

B. Thermodynamic quantities
To ensure consistency in the calculation, it is imperative to rigorously incorporate the thermodynamics of the LSM. This will directly reflect chiral symmetry restoration and will be critical in determining the transport coefficients. The interactions in the LSM are evaluated in the mean field limit, which can in turn be absorbed in a mass redefinition. Detailed discussions of thermal effective masses exist in the literature, e.g. in [16] (based on the methods of [17][18][19]); we present a brief overview of the method here for the sake of completeness. In this study, we use classical (Boltzmann) statistics which results in simplifications of statistical factors throughout.
Mean field equations of motion are derived by taking the thermal average of the Euler-Lagrange equation with respect to a general field, which we denote ψ a = { π, σ}.
but this can be simplified further by recognizing as ψ a = 0. Therefore, the thermal average equation of motion becomes ∂U ∂ψ a = 0.
Through some further calculations, it is clear that the non-trivial solution is In order to calculate the effective masses, we solve the equation of motion for each field in the Lagrangian, producing three coupled equations that must be solved selfconsistently. These are where the thermal average of the fields is given by and solutions are shown in Fig. 1. We show results for vacuum masses of 400 and 900 MeV, the two extreme values of the range considered here. While the framework we develop is general, the evaluation of scattering matrix elements, thermal effective masses, and mean-field effects is done using the LSM. Importantly, in the rest of this work all energies and masses are thermal, i.e. the singleparticle energies are E a = p 2 a + m 2 a , where m hereon denotes the effective thermal masses discussed previously. When needed, a vacuum mass is now written as m 0,a .
As mentioned earlier, the incorporation of explicit symmetry breaking qualitatively and radically alters the symmetry restoration at high temperatures. Fig. 1 confirms that restoration now takes place over a broad crossover region. It is interesting to recall that our current understanding of the QCD transition from partonic to confined hadronic degrees of freedom is also a crossover (at zero net baryon density), albeit occurring at a lower temperature [15].
We also wish to calculate thermodynamic quantities, such as the pressure, entropy density, energy density, heat capacity, and the speed of sound. The total pressure, entropy density, energy density, and heat capacity are the sum of the contribution from each species.
The form of T µν , discussed in more detail in Sec. III A, yields where we have integrated by parts, assuming a vanishing boundary term. The energy density equation is trivial. We next provide the entropy density Finally, the expression for heat capacity is given below.
The calculation of the speed of sound is more involved and is shown in detail in [20]. However, it is important to consider that the speed of sound is a thermodynamic quantity as is the consistency condition (also detailed in [20]). Thus, these become involved in Landau matching, which we address in Appendix A. This discussion gains additional importance as we outline difficulties in performing exact calculations of the bulk viscosity in single and multi-component gases in Sec. IV.

III. THE RELAXATION TIME APPROXIMATION
In order to discuss transport coefficients, we begin with a discussion of transport with the Boltzmann equation (BE).

A. Boltzmann equation and some definitions
The Boltzmann transport equation can be written as where we have used the convenient shorthand and with particles a, b incoming and c, d outgoing. The LSM is not explicitly restricted to 2 ↔ 2 processes, but we will only consider these processes in this work as already discussed in Sec. II A. We write the symmetric energy-momentum tensor T µν as T µν = −P g µν + wu µ u ν + ∆T µν (34) with correction terms related to the dissipation properties We choose the Landau frame, where u µ is the energy transport velocity. In our convention, we use the mostly negative Minkowski metric signature (+, −, −, −) and define u such that u 2 = 1.
In this discussion, P is the pressure, w is the enthalphy density w = T s = P + , s is the entropy density, and is the energy density. We have defined the projection tensor and derivative normal to u µ as respectively. We also adopt the convention that Latin indices either label species or refer to three-components of four-vectors and are thus not affected by raising and lowering.

B. Approximate solution
The relaxation time approximation 1 is a popular approximation to a solution of the BE, and is commonly used in calculating transport coefficients [21]. However, its validity decreases as the relaxation time increases and as such it is arguably uncontrolled. To lay the foundation for a more quantitative discussion, we begin by developing this formalism in some detail before introducing an exact solution for transport coefficients in the LSM. We will then be able to precisely quantify this approximation in a consistent way. We derive the forms of the shear and bulk viscosity using the Chapman-Enskog expansion [22]. We assume a small deviation from equilibrium; we satz an that this is of the form where φ(p) quantifies deviations from equilibrium. We typically suppress the momentum dependence for clarity.
As a result of the maximally-general tensor decomposition of the form of Eq. (35), it is natural to construct φ such that it has the same decomposition where C a µν = C a p µ p ν and both C a and A a in general depend on the scalar u α p α .
We additionally require non-decrease of entropy. In the local rest frame of the fluid, the change in the entropy is given by which requires that both shear viscosity η and bulk viscosity ζ be non-negative.
Assuming that interactions are localized and are point or contact interactions, the energy-momentum tensor can be written as a sum of independent contributions. We can then return to the form T µν = T µν eq +∆T µν by writing this as For convergence, the deviation φ must be perturbatively small, i.e. |φ| 1. We now find the form of the shear and bulk viscosity in the local rest frame of the fluid. This is done by equating the two expressions for the dissipative part of the energymomentum tensor.
To calculate shear viscosity, we investigate a purely shear flow in a single direction. Without loss of generality, we choose u k = (u x (y), 0, 0). Applying this flow to both expressions for the dissipative part of T µν reduces them to and it is possible to identify shear viscosity as The same method can be used to isolate the contribution of bulk viscosity that, after some manipulations, can be written as The terms C a and A a are found by manipulations of the Boltzmann equation. Beginning with the LHS, The last line assumes that the off-equilibrium component φ a (x, p) is small. Now we move on to the RHS. Recall Eq. (39) and that, as a consequence of classical statistics, the product of the distribution functions before and after a collision are equivalent.
Keeping terms linear in φ, we return to the Boltzmann equation.
The first task is to compute the LHS. We begin with the calculation of ∂ µ f eq a where f eq a = exp(−u ν p ν /T ) and obtain We now rewrite the Boltzmann equation using some thermodynamic quantities. These are standard, but a detailed treatment can be found in [20] where many standard results are collected. Using the speed of sound, we may now make some progress in rewriting the Boltzmann equation into a more convenient form At this stage it is necessary to substitute the structure of φ and group terms in Eq. (52), which now reads To do this in a consistent way, it is necessary to consider this equation term-by-term. It is useful to define the shorthand Rewriting the φ coefficients of each term, we find that In taking the relaxation time approximation, we suppose that all particles are in equilibrium except for species a, which is out of equilibrium by a perturbatively small amount, f a = f eq a + δf a . To that order, where is the interaction frequency. The expressions and phase space for ω a are shown in detail in [20]. We define the relaxation time to be and the deviation δf a is that in Eq. (39). Both the interaction frequency and relaxation time are energydependent, but we suppress this in the notation for clarity. To find the viscosities, it is necessary to substitute the deviation φ into the Boltzmann equation. Thus, we find Therefore, we arrive at It is then possible to make the following identifications The shear viscosity is now readily calculated using Eq. (46) The inclusion of mean-field effects and ensuring thermodynamic consistency makes the evaluation of the bulk viscosity slightly more complicated than that of the other transport coefficients discussed in this work. This is addressed in Appendix A with details in [20]. The result for the bulk viscosity is With the the incorporation of mean field effects, the bulk viscosity meets the Landau matching condition. We also display the result for the electrical conductivity in the RTA [23][24][25][26] for later use: The only modification of η and of the electrical conductivity σ el due to mean fields comes from the presence of effective masses in phase space considerations. Finally, τ a is calculated numerically using Eq. (62) and inserted into these workings.

IV. VARIATIONAL METHOD
It is difficult to calculate systematic corrections to the RTA. As it is important to separate the impact of approximations from those of the transport properties of the medium itself, we turn to a method for calculating transport coefficients exactly in the limit of the linearized Boltzmann equation. We extend the variational method of [10,27,28] to massive theories and remove the small momentum transfer approximation so the technique is applicable to inelastic processes. We will follow the same notation and provide an overview of the method for completeness. For details of incorporating masses into this framework, see Appendix B.
We begin by laying out some notation and motivations of the general form. To do this, we define a collision operator that is already linear in the deviation from equilibrium when using Boltzmann statistics.
It may then be shown that the Boltzmann equation, to first order in gradients, is a linear integro-differential equation, We now move on to interpreting the LHS. In the local fluid rest frame it may be written as where q a is the conserved charge of the quantity of interest and we have separated the angular and spatial dependence into I i...j , the unique rotationally covariant tensor, and X i...j , the spatial tensor denoting the driving field: and Due to the rotational invariance of the collision operator C, the departure from equilibrium and the driving field must have the same angular form. Thus, the deviation that will solve the Boltzmann equation must be δf a (p a , x) = β 2 X i···j (x)I i···j (p a )χ a (|p a |).
In order to solve this, we additionally define an inner product as that will allow us to construct a functional, Q. We will expand this functional in a variational basis and maximize the expanded functional to determine variational coefficients and extract transport coefficients. We define a functional Q[χ i...j ] such that it is extremal when χ a (p) satisfy the linear Boltzmann equation In the above, and C is the linearized collision operator while χ a (p) is a rotationally invariant function depending only on excitation energy. χ a (p) is what we will expand in a convergent variational basis in order to calculate transport coefficients.
The source in Eq. (81) can be written in terms of an expansion basis φ m (see Appendix B) as and the collision term may be writteñ This may be assembled into the maximized functional Q max From this maximized functional, we are able to calculate the transport coefficients. Keeping in mind the angular momentum structure of Eq. (81), the expressions for transport coefficients are where details of constructing the different components are given in Appendix B. While we only calculate those applicable to our theory, the work in [10] also calculates the flavor diffusion constants; our extended framework can also calculate this for massive theories with inelastic processes.
In the bulk viscosity calculation, it is necessary to orthogonalize to the zero modes. This is because the bulk viscosity is a spin 0 quantity (a scalar), while the other transport coefficients are spin 1 (a vector) or spin 2 (a tensor). Exact zero modes are only present in the scalar quantity because when taking the dot product of the momenta, angular factors in the spin 1 or 2 modes break the degeneracy.
A zero mode is an eigenvector with a vanishing eigenvalue that presents a problem to the inversion of the matrix. It corresponds to conserved quantities in the system, so in a system with only number-conserving processes, two zero modes exist: one corresponding to the conservation of energy and the other to the conservation of total particle number. At the order we consider, our theory exactly conserves particle number, which means that the source must be absolutely orthogonal to both the energy and number conservation zero modes. This can only be accomplished by detailed accounting for a chemical potential that we do not develop, or by explicitly knowing the form of the zero mode. The exact form of the zero mode is known for a single-component gas [28] but this result is not applicable to multi-component gases and rigorous study of this is beyond the scope of this work. To avoid detailed treatment of a chemical potential, one would have to make the system particle number non-conserving. One would then need to consider 1 ↔ 2 processes or 2 ↔ 4 processes, which would mean the return of poles in the matrix elements, poorly-defined expansion to higher orders in the coupling, and/or a reevaluation of the approximations in this work. Due to a large coupling constant in the LSM, higher-order expansions are not well-defined. As a result, these considerations are beyond the scope of this work and the development of techniques to address them will be pursued elsewhere.
The details of the LSM itself are incorporated in the thermal masses, collision term, and source of the variational framework that has been developed in this section. As a result, this method remains completely general and quantitative comparisons can be made to other general methods and calculations, such as the RTA and perturbative QCD.

V. RESULTS
We present the numerical results beginning with the thermodynamics of the LSM. This reveals that the behaviour near chiral symmetry restoration contains interesting physics that can be explored in more physical models. Importantly, the features of the chiral symmetry restoration have a consequence upon the system's thermodynamics. Having verified that we are able to resolve the dynamics we expect in the thermodynamic quantities, we compute transport coefficients beginning with the vector and tensor quantities σ/T and η/s in both the relaxation time approximation and the variational method as these do not possess exact zero modes. We conclude by calculating ζ/s in the relaxation time approximation, leaving the treatment of zero modes for future work.
Integrals are evaluated using Vegas adaptive Monte Carlo [29] and numerical uncertainties are propagated through nested calculations 2 . We compare to hadron gas calculations, chiral perturbation theory, and pQCD calculations and show that the LSM demonstrates key features of these other approaches. We also use this calculation to provide insight to the possible parameters of a sigma meson, keeping in mind the characteristics of the f 0 (500) [30] and also the possible caveats associated with identifying this σ field with a physical particle [31].

A. Thermodynamic quantities
Once the effective masses shown in Fig. 1 are obtained via a self-consistent numerical optimization [32], one is able to calculate thermodynamic quantities such as energy density, pressure, and entropy density (Fig. 2); heat capacity (Fig. 3); and the speed of sound (Fig. 4). An important and clearly visible feature is the different behaviors of the thermodynamic quantities and heat capacity at different values of the vacuum sigma mass: the higher the vacuum sigma mass, the more suddenly chiral symmetry is restored, producing a peak in the heat capacity and a corresponding trough in the speed of sound. Similarly, the thermodynamic quantities have a more pronounced behavior at m 0,σ = 900 MeV than they do at 600 MeV. This behavior has a clear impact on the transport coefficients, but none more so than the bulk viscosity, which is highly sensitive to the conformality of the system as measured by the speed of sound.

B. Electrical conductivity
We now turn our attention to the transport coefficients and begin with the electrical conductivity, which quantifies the conduction properties of the medium. The DC conductivity is the real, static part of the complex conductivity tensor, σ el = Re{lim ω→0 + σ ii (ω, 0)}, and is related to the electrical field through Ohm's law: J EM = σ el E. One may use linear response theory [12] to derive a Kubo formula for the conductivity tensor, which features the electromagnetic current operator: Importantly, the emission of electromagnetic radiation is also regulated by the current-current correlator [33]. Therefore, in addition to its intrinsic interest from the point of view of transport, σ el can provide information on the ability of the hot and strongly interacting medium to emit soft photons: where # is a numerical pre-factor. Rigorous quantitative control of the electrical conductivity can thereby also provide constraints on the phenomenology of soft electromagnetic radiation.
The electrical conductivity has been the subject of previous studies with both hadronic (confined) and partonic degrees of freedom [10,[34][35][36][37][38][39][40][41][42]. The DC electrical conductivity was not previously studied in the LSM, although its vector structure is conducive to calculation using the variational method we have extended. We provide calculations here both in the relaxation time approximation and in the variational method, noting the differences between them.
Some extractions of the electrical conductivity of hadron gases are available, making this an ideal choice for further study and for validation of our method. An additional benefit comes in computational efficiency: since the conserved charge for electrical conductivity is simply the electrical charge of particle a, the source (Eq. (81)) for neutral particles is identically 0. This naturally simplifies the structure of the maximization, since these components will not contribute and do not have to be calculated in either the exact variational or RTA approaches. The σ meson contributes resistance only through interactions with π.
As current calculations of transport parameters of strongly interacting matter have not converged to a set of well-defined values, it is prudent to learn from this apparent lack of unity and to compare with various approaches and models. We begin with the electrical conductivity here. Fig. 5 contains results for the electrical conductivity over temperature plotted as a function of temperature, calculated by both the RTA and variational techniques. Some previous calculations of σ el for a hadron gas are also shown. For ease of comparison between different models and results we use a scaled temperature. Since there is no genuine phase transition occurring with increasing temperature, one must adopt an operational definition of critical temperature "T c ". Possible choices are the temperature where effective masses are minimized (see Fig. 1), or when c V /T 3 peaks (see Fig. 3). In addition, these answers would vary, depending on the choice of m 0,σ . Choosing the first criterion leads to T c = (242, 245, 259) MeV, for m 0,σ = (400, 600, 900) MeV, respectively. The second yields T c = (200, 219, 245) MeV. We will choose the first scheme and the intermediate value of the scalar-isoscalar mass. Therefore, for calculations that involve the LSM, we set T c = 245 MeV. The "critical temperature" in other approaches shown here is the one reported using each of those models.
Our results depend strongly on the value of the vacuum σ mass, as seen in Fig. 5. For comparison, we also show results obtained with the PHSD [37] and NJL [38] models. As is the case for the magnitude of the transport parameters, there is currently no consensus on the details of their temperature dependence. A hint of critical behavior is observed in the results with large σ mass, which is also seen in the results of PHSD, but that is where similarities end. This spread in theoretical results is also seen in more extensive compilations [43], which makes it difficult to single-out a preferred value of m 0,σ by comparing between theoretical results. However, the magnitude of the vast majority of results obtained in the field fall within the range spanned by the value of m 0,σ explored here.
We directly compare the calculations in the relaxation time approximation and the variational method and we find that the two are within a factor ∼ 3 of each other. This deviation between results obtained with the two  [37] and the Nambu-Jona-Lasinio model [38].
techniques using the same model is a feature seen in all of our calculations and is also observed in others [9]. We also note some difference in the parametric behaviours, particularly at low-T and near T c . This could be added to the growing body of evidence cautioning against using the RTA for precise quantitative studies of stronglyinteracting systems.

C. Shear viscosity
The shape and value of the shear viscosity to entropy density ratio for strongly-interacting matter is a topic of immense interest. Our results (Fig. 6) again reveal that relaxation time calculations of the minimum value of η/s can be as much as a factor of 3 lower than the value of a more precise calculation within the same theory. While the RTA calculation of η/s with a vacuum σ mass of 900 MeV approaches the KSS result of 1/4π [44], the same calculation in the variational method does not.
Many calculations of η/s exist in the contemporary scientific literature and we again can not show an inclusive compendium here. It is appropriate to show a direct comparison to another similar calculation in the linear sigma model [9]. As seen on Fig. 6, the numerical results reported here are close those seen in [9], with a difference increasing with decreasing temperatures. As T shrinks to the lowest values explored here, an apparent plateau in η/s is seen in [9] -and perhaps even a decreasewhereas the values of specific shear viscosity calculated here follow an almost perfect exponential increase. In addition, the figure contains results obtained with BAMPS -a relativistic Boltzmann equation solver [45], the Dual Quasi-Particle Model [38], and with the NJL model [38]. These approaches again all yield results that differ over the range of temperatures chosen here. This is a recurrent theme and is consistent with the current state of affairs in the field. Several of those approaches use different degrees of freedom but parton-hadron duality for T ∼ T c has the potential to minimize these differences. Finally, calculations of η/s exist for lower values of T [46,47]. Given that the range of validity of those decreases with increasing temperature 3 and our results have used the s → ∞ limit, showing that they agree in some range of T (they do) has questionable value. The specific shear viscosity, η/s, is calculated in the relaxation time approximation and with the functional variatio method, and compared results using BAMPS [45], the dynamical quasi-particle model (DQPM) and NJL of [38], and AdS-CFT [44]. We additionally compare to the variational technique of [9], labelled as "Alt. Var.".
In most hydrodynamic applications, quantities are often fixed by the ratios of different transport coefficients. If the system can be characterized by a single relaxation time then, once we know one transport coefficient, others can be deduced by knowing these ratios. A natural question is that, while the use of the RTA is questionable for quantitative studies, perhaps the ratios between transport coefficients calculated using the RTA are close to those obtained using the variational technique. We explicitly considered the ratio of electrical conductivity to shear viscosity, σ el /η. The value of the ratio using the RTA depends on choice of m 0,σ . Choosing the lower value more closely reproducing the ratio seen in the functional variation calculation. Even with the larger scalar-isoscalar masses, the ratios from the RTA are within ∼ 20% of those obtained using functional variations. The ratios obtained with both techniques are almost flat above T c . We will show these results situated in context and in more detail in upcoming work.

D. Bulk viscosity
We present the calculation of the linear sigma model bulk viscosity in the RTA in Fig. 7. The only computational advantage of the RTA presents itself here: zero modes of the collision matrix are not present. However, as we have established through the calculations of electrical conductivity and shear viscosity, we can produce only an estimate that exhibits the broad dynamics of the more precise calculation. The bulk viscosity to entropy density ratio for a variety of vacuum sigma masses in the relaxation time approximation. Comparisons are shown to pQCD [28].
We do not calculate bulk viscosity in the variational technique because of zero modes, as discussed in Sec. IV. We instead calculate the bulk viscosity in the relaxation time approximation as this approximation bypasses the issue. We have shown in Figs. 5, 6 evidence suggesting that the RTA is insufficiently precise for detailed studies, and discussions in the literature also suggest the RTA is insufficient for calculations of the bulk viscosity [48]. With those caveats in mind, one observes that the peak exhibited in Fig. 7 by the LSM RTA calculation with the larger σ mass is approximately a factor of 3 lower than that used in some hydrodynamics-driven phenomenological analyses [49,50]. This must be understood in context: as made clear in Ref. [51], the bulk viscosity is currently not well constrained by systematic analyses of experimental data.
As stated many times, the convergence of results for transport coefficients is currently not at hand. This is especially true for the bulk viscosity and it is clear that the proper exact treatment of bulk viscosity should be a priority for future investigations. One of the reasons for this emphasis is the special dual role enjoyed by ζ/s. On one hand, in dynamical simulations the bulk viscosity reflects the resistance of the hydrodynamic system to volumetric deformation and therefore has a direct impact on the average transverse momentum of measured hadrons [49]. On the other hand, it can also be related to the nonconformality of the underlying theory, QCD [52].

VI. CONCLUSIONS
In this work, we have presented a comprehensive general framework for the calculation of transport coefficients in massive quasiparticle theories at finite temperature with inelastic processes. Using the linear sigma model, hadronic transport coefficients were calculated and compared to results from other theories. We have produced the first calculations of the electrical conductivity in the linear sigma model using the RTA and functional variational techniques. We have shown calculations for the shear viscosity in both methods, while a calculation of the bulk viscosity was only produced in the relaxation time approximation. Reasons for this choice were provided in detail. In all cases, we observe that the RTA and the variational results can differ by a factor of ∼3. This should be interpreted as a cautionary flag for all calculations and should preface most -if not all -current theoretical attempts at a quantitative characterization of strongly-interacting matter out of equilibrium. This reinforces the need for more precise calculations using realistic models of hadron and parton dynamics and for rigorous and systematic phenomenological extractions of transport coefficients from experimental relativistic heavy-ion data. This study should also provide impetus for further phenomenological applications by influencing, and even providing, prior distributions in Bayesian analyses [6].

VII. ACKNOWLEDGMENTS
We thank S. Hauksson for many useful discussions, P. Arnold for helpful correspondence, and G. Moore for an informative discussion on bulk viscosity. We thank P. Chakraborty  This section provides a summary of the Landau matching in [20], in turn based on that in [9].
Using Eq. (59), it is simple to rewrite Eq. (56) as and extract the functions C a µν and A a . Details of this are provided in [20].
These "departure functions" characterize the shear and bulk departures from equilibrium. An important subtlety is that in the departure function decomposition of the Boltzmann equation, there is not a unique solution to A a in Eq. (A3). Infinite solutions can be generated by shifting a particular solution A(E), e.g.
where a is an arbitrary constant associated with particle conservation and b is an arbitrary constant associated with energy conservation. This degree of freedom is related to the fact that the Boltzmann equation admits summational invariants [23]. We are restricting our scope to that with no chemical potential, thus all a are 0 as there are no particle conservation considerations. Thus, a particular solution A a to the above can be related to all other solutions Thus, by considering Eq. (43) and u µ ∆T µν = 0, it is straightforward to conclude that Recalling the definitions of the single particle contributions to thermodynamic quantities, it can be seen that We now consider the variational impact from the mean field effects on Landau matching. It is of particular importance to treat δf correctly [48]. We take a small deviation from the equilibrium distribution function: The single particle energy also takes an off-equilibrium shift If the equilibrium distribution function is expressed as a function of the true energy (including off-equilibrium shifts), then and The shift in the energy density is therefore We now return to the consideration of ∆T 00 .
We obtain this final result by recalling 2T dT = dT 2 . By definition, u 2 = 1. As a result, we satz that u µ u ν is the prefactor of the T 2 term an. This follows from physical arguments, such as that this effective mass dependence does not have an impact on the pressure in the local rest frame. Thus, we generalize as follows: Recall the definitions of shear and bulk viscosity We now impose Landau matching with the new results on the effects of mean fields. This only modifies bulk viscosity ζ. As before, if a particular solution does not meet the Landau-Lifshitz matching condition, it can be made to comply by adding/subtracting a linear energy term We now use a relation to simplify the process and find the final result.
Compiling Eq. (A14, A15, A16) results in Using this expression, we may rearrange slightly and constrain the solution with the use of Landau matching.
a pa And as a result, it is possible to constrain b in Eq. (A19) We may finally substitute A a (E a ) = A a − bE a into the expression for bulk viscosity and we can conclude with an expression for bulk viscosity that meets the Landau-Lifshitz condition by construction Substituting the solution for A a from the relaxation time approximation yields Eq. (70). We now encounter the first point where the masses enter the theory: in the use of the convenient delta function expansion where we recognize w as the energy transfer. We additionally define a momentum transfer q such that p c = p a + q, p d = p b − q. To remove the integration over cosines, we must expand these delta functions and derive limits. and Implicit in each of these final delta functions is a theta function that ensures that energy is conserved; Θ(E b −w) in Eq. (B4) and Θ(E a + w) in Eq. (B6). The masses do not modify the Jacobians, so these are the same for the massless case. Performing the cosine integrals trivially using the delta function will yield limits to ensure that the delta functions are satisfied. In order to find the limits, we solve the inequality If we consider only elastic processes, we can make the further simplification that m a = m c and m b = m d . Adding the assumptions that E pa,p b ≥ 0, m pa,p b = m pc,p d ≥ 0 yields the following bounds: Our calculations also include inelastic collisions, but the limits for inelastic processes are significantly more involved and are excluded due to space constraints. However, the phase space can be thoroughly explored by taking the clear kinematic limits and checking numerically to see that the conditions from the delta functions are satisfied. We include details here to clarify what differences arise without taking the small momentum transfer approximation in a massive case. The biggest consideration is in the treatment of angular cross-terms, which are radically different. Returning to the integral, χ(p) can then be expanded in a basis, where we successfully use the same basis as that in [10], where N is the size of the basis, which has been chosen such that it converges quickly. The coefficients a m are maximized in order to maximize the functional Q, which we will denote as Q max . We must now also repeat the treatment of delta functions for cosines to account for masses. All our results for limits and angular factors recover those of [10] when m → 0.
cos θ paq = 2wE a + t + m 2 a − m 2 c 2p a q (B14) cos θ pap b = cos θ paq cos θ p b q + sin θ paq sin θ p b q cos φ (B20) cos θ pap d = cos θ paq cos θ p d q + sin θ paq sin θ p d q cos φ (B21) cos θ pcp b = cos θ pcq cos θ p b q + sin θ pcq sin θ p b q cos φ (B22) cos θ pcp d = cos θ pcq cos θ p d q + sin θ pcq sin θ p d q cos φ (B23) What remains is to properly deal with the χ factor. Knowing that the matrix elements in the LSM are constants has the potential to simplify the forms of some of the φ integral, but this is properly done at the end of the manipulations.
In not taking the small momentum transfer approximation, we must consider all of the angular factors. As a result, we require proper treatment of angular to account for the inelastic processes that take place in the linear sigma model. Taking the analogy to φ as shown before, the right hand side becomes φ a m (p a )φ b n (p b ) + φ a n (p a )φ b m (p b ) P l (cos θ pap b ). (B26) This occurs for all the cross terms and angles between particles. The result of this computation is inserted into Eq. (B11), but is not included here due to length. Care must be taken with angular cross-terms in Eq. (B11). These are explicitly I i...j (p a )I i...j (p b ) = P l (cos θ pap b ) (B27) where P l (cos θ) is the l th Legendre polynomial and refers to the spin of the underlying transport coefficient. As discussed in Sec. IV, the appropriate values are 0 for the bulk viscosity (a scalar), 1 for the conductivity (a vector), and 2 for the shear viscosity (a tensor). The electrical conductivity, at least in the case of the LSM, is a simple case as the source for neutral particles is explicitly 0 and we expect δf to be identical only for particles with identical mass and reactions. This makes it possible to easily identify components of the collision matrix that do not contribute. The shear and bulk viscosity will require appropriate consideration of different δf for π and σ mesons.
Most importantly, in shear and bulk viscosity, we must consider a cross-coupling between the sigma and the pi with appropriate consideration for their different departures from equilibrium. This means that instead of an N × N collision matrix, we instead have a 2N × 2N collision matrix. Note that it is not 4N × 4N because for the viscosities, the pions are identical. We therefore sum the three pions in the π segment of the above matrix and the πσ interactions in the πσ sections above. In the expansion in Eq. (B12), one species will have expansion coefficients a m , but the other will then have the expansion coefficient a N +m . We now turn to a discussion of the source.
The assembly of the collision matrixC and the source vectorS require some details to ensure that interactions are properly implemented. The assembly of these into block components most clearly demonstrates this process.C is assembled is by calculating each component with only one species C σ , C π and then the interactions between particle species, C σπ , and then compiling into blocks within the respective matrices.
The notation below is somewhat simplified, C π can be further expanded into C π =    C π 0 C π 0 π + C π 0 π − C π + π 0 C π + C π + π − C π − C π − π + C π − π 0    although this becomes more cumbersome to show in the full collision matrix. C σ contains all reactions with only σ and the off diagonals contain pion-sigma reactions. Similarly, the source is comprised of sub-components that are compiled into The term S π may be decomposed in the same manner as C π . Thus, Q max becomes