Relativistic BGK hydrodynamics

Relativistic Boltzmann equation governs the space-time evolution of the single particle phase-space distribution function of a relativistic system. Moreover, suitable moments of the Boltzmann equation are capable of describing the collective dynamics of the system. Therefore, it has been extensively used to derive equations of relativistic dissipative hydrodynamics and obtain expressions for the transport coefficients [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. The collision term in the Boltzmann equation, which describes change in the phase-space distribution due to the collisions of particles, makes it a complicated integro-differential equation. In order to circumvent this issue, several approximations have been suggested to simplify the collision term in the linearized regime [16, 17, 18, 19, 20].

Bhatnagar-Gross-Krook [16], and independently Welander [17], proposed a relaxation type model for the collision term, which is commonly known as the BGK model. This model was further simplified by Marle [18] and Anderson-Witting [19] to calculate the transport coefficients. In the non-relativistic limit, Marle’s formulation leads to the same transport coefficient as the BGK model but fails in the relativistic limit. On the other hand, the Anderson-Witting model, also known as the relaxation-time approximation (RTA), is better suited in the relativistic limit. The RTA has been employed extensively in several areas of physics with considerable success and has been widely employed in the formulation of relativistic dissipative hydrodynamics [8, 9, 10, 11, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33].

The RTA Boltzmann equation has provided remarkable insights into the causal theory of relativistic hydrodynamics as well as a simple yet meaningful picture of the collision mechanism in a non-equilibrium system. On the other hand, the BGK collision term ensures conservation of net particle four-current by construction, and is the precursor to RTA. Moreover, collision kernel in RTA involves only one timescale, the relaxation time, which governs all dissipation in the system. As a consequence, if the coefficient of shear viscosity is specified, the coefficient of bulk viscosity gets fixed automatically via the relaxation time. While the RTA has been employed extensively, a consistent formulation of relativistic dissipative hydrodynamics with the BGK collision term is relatively less explored. This may be attributed to the fact that the BGK collision kernel is ill defined for relativistic systems without a conserved net particle four-current. This has limited the use of the BGK collision kernel to the studies related to flow of particle number and/or charge [34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44].

In this article, we take the first step towards formulating a consistent framework of relativistic dissipative hydrodynamics using the BGK collision kernel. Furthermore, we propose a modified BGK collisions kernel (MBGK), which is well defined even in the absence of conserved particle four-current and is better suited for the formulation of the relativistic dissipative hydrodynamics. We find that there exists a free scalar parameter arising from the freedom of matching condition. This affects the scalar dissipation in the system, i.e., the coefficient of bulk viscosity. This feature allows us to specify the coefficient of bulk viscosity, independently from shear viscosity, as opposed to RTA. We study the effect on bulk viscosity in several different scenarios.

The conserved net particle four-current, $N^{\mu}$, and the energy-momentum tensor, $T^{\mu\nu}$, of a system can be expressed in terms of the single particle phase-space distribution function and the hydrodynamic variables as [45],

where the Lorentz invariant momentum integral measure is defined as $\mathrm{dP}=g\,d^{3}p/\left[(2\pi)^{3}E\right]$ with $g$ being the degeneracy factor and $E=\sqrt{|{\bf p}|^{2}+m^{2}}$ being the on-shell energy of the constituent particle of the medium with three-momentum ${\bf p}$ and mass $m$. Here $f\equiv f(x,p)$ and $\bar{f}\equiv\bar{f}(x,p)$ are the phase-space distribution functions for particles and anti-particles, respectively. In the above equations, $n$ is the net particle number density, $\epsilon$ is the energy density, $P_{0}$ is the equilibrium pressure, $n^{\mu}$ is the particle diffusion four-current, $\delta P$ is the correction to the isotropic pressure, and $\pi^{\mu\nu}$ is the shear stress tensor. We note that the fluid four-velocity $u^{\mu}$ has been defined in the Landau frame, $u_{\mu}T^{\mu\nu}=\epsilon u^{\nu}$. We also define $\Delta^{\mu\nu}\equiv g^{\mu\nu}-u^{\mu}u^{\nu}$ as the projection operator orthogonal to $u^{\mu}$. In this article, we will be working in a flat space-time with metric tensor defined as, $g^{\mu\nu}={\rm diag}(1,-1,-1,-1)$.

Hydrodynamic equations are essentially the equations for conservation of net particle four current, $\partial_{\mu}N^{\mu}=0$, and energy-momentum tensor, $\partial_{\mu}T^{\mu\nu}=0$. Using the expressions of $N^{\mu}$ and $T^{\mu\nu}$ from Eqs. (1) and (2), the hydrodynamic equations can be obtained as,

where we use the standard notation, $\dot{A}\equiv u_{\mu}\partial^{\mu}A$ for the co-moving derivatives, $\nabla^{\alpha}\equiv\Delta^{\alpha\beta}\partial_{\beta}$ for the space-like derivatives, $\theta=\partial_{\mu}u^{\mu}$ for the expansion scalar, and $\sigma^{\mu\nu}=\frac{1}{2}\left(\nabla^{\mu}u^{\nu}+\nabla^{\nu}u^{\mu}\right)-\frac{1}{3}\Delta^{\mu\nu}\theta$ for the velocity stress-tensor.

To express the conserved net particle four-current and the energy-momentum tensor in terms of hydrodynamic variables in Eqs. (1) and (2), we chose Landau frame to define the fluid four-velocity. Additionally, the net-number density and energy density of a non-equilibrium system needs to be defined using the so called matching conditions. We relate these non-equilibrium quantities with their equilibrium values as

where $n_{0}$ and $\epsilon_{0}$ are the equilibrium net-number density and the energy density, respectively, and, $\delta n$, $\delta\epsilon$ are the corresponding non-equilibrium corrections. For a system which is out-of-equilibrium, the distribution function can be written as $f=f_{0}+\delta f$, where $f_{0}$ is the equilibrium distribution function and $\delta f$ is the non-equilibrium correction. In the present work, we consider the equilibrium distribution function to be of the classical Maxwell-Juttner form, $f_{0}=\exp(-\beta\,u\cdot p+\alpha)$, where $\beta\equiv 1/T$ is the inverse temperature, $\alpha\equiv\mu/T$ is the ratio of chemical potential to temperature and $u\cdot p\equiv u_{\mu}p^{\mu}$. The equilibrium distribution for anti-particles is also taken to be of the Maxwell-Juttner form with $\alpha\to-\alpha$.

We can now express the equilibrium hydrodynamic quantities in terms of the equilibrium distribution function as,

Similarly, the non-equilibrium quantities can be expressed as

where $\Delta^{\mu\nu}_{\alpha\beta}\equiv\frac{1}{2}(\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}+\Delta^{\mu}_{\beta}\Delta^{\nu}_{\alpha})-\frac{1}{3}\Delta^{\mu\nu}\Delta_{\alpha\beta}$ is a traceless symmetric projection operator orthogonal to $u^{\mu}$ as well as $\Delta^{\mu\nu}$. In order to calculate the non-equilibrium quantities defined in Eqs. (10)-(14), we require the out-of-equilibrium correction to the distribution function i.e., $\delta f$ and $\delta\bar{f}$. For this purpose in the following, we consider the Boltzmann equation with BGK collision kernel, which describes the spacetime evolution of the distribution function.

The covariant Boltzmann equation, in absence of any force term or mean-field interaction term, is given by,

for a single species of particles and its antiparticles. In the above equation, $C[f,\bar{f}]$ and $\bar{C}[f,\bar{f}]$ are the collision kernels that contain the microscopic information of the scattering processes. For the formulation of relativistic hydrodynamics from the kinetic theory of unpolarized particles, the collision kernel of the Boltzmann equation must satisfy certain properties. Firstly, the collision kernel must vanish for a system in equilibrium, i.e., $C[f_{0},\bar{f}_{0}]=\bar{C}[f_{0},\bar{f}_{0}]=0$. Further, in order to satisfy the fundamental conservation equations in the microscopic interactions, the zeroth and the first moments of the collision kernel must vanish, i.e., $\int\mathrm{dP}\,C=0$ and $\int\mathrm{dP}\,p^{\mu}\,C=0$. Vanishing of the zeroth moment and the first moment of the collision kernel follows from the net particle four-current conservation and the energy-momentum conservation, respectively.

In the present work, we consider the BGK collision kernel which has the advantage that the particle four-current is conserved by construction. The relativistic Boltzmann equation with BGK collision kernel for particles can be written as [40, 41, 42, 43, 44],

and similarly for anti-particles with $f\to\bar{f}$ and $f_{0}\to\bar{f}_{0}$. Here, $\tau_{\mathrm{R}}$ is a relaxation time like parameter¹¹1A more conventional notation is the collision frequency which is defined as $\nu=1/\tau_{\mathrm{R}}$. which we assume to be the same for particles and anti-particles. It is easy to verify that the conservation of net particle four-current, defined in Eq. (1), follows from the zeroth moment of the above equations. The first moment of the above equations should lead to the conservation of the energy-momentum tensor, defined in Eq. (2). However, we find that the first moment of the Boltzmann equation, Eq. (16), leads to,

which does not vanish automatically.

In order to have energy-momentum conservation fulfilled by the Boltzmann equation with the BGK collision kernel, Eq. (16), we require that

which we identify as one matching condition. Using Eq. (6), the above constraint relation can be expressed in terms of the out-of-equilibrium contribution to energy density and net number density as

Note that, in the present work, the fluid four-velocity has been defined using Landau frame condition. In order to derive a consistent hydrodynamic description from kinetic theory, one further requires two independent matching conditions to define two hydrodynamic field variables namely, the temperature, $T(x)$ and the chemical potential, $\mu(x)$ [46]. This requirement essentially stems from the absence of a first-principle microscopic definition for nonequilibrium temperature and chemical potential. In the case of traditional RTA Boltzmann equation, $p^{\mu}\partial_{\mu}f=-\frac{(u\cdot p)}{\tau_{R}}(f-f_{0})$, the two constraints necessary for the energy-momentum and net particle four-current conservation in the Landau frame are $\epsilon=\epsilon_{0}$ and $n=n_{0}$. Moreover, these matching conditions are also sufficient to define $T(x)$ and $\mu(x)$ in the RTA case. However, in the BGK approach, the net particle four-current conservation is automatic and does not lead to any constraint. Thus the only necessary condition obtained by demanding the energy-momentum conservation is not sufficient to specify both $T(x)$ and $\mu(x)$. In fact, Eq. (18) specifies one relationship among $T(x)$ and $\mu(x)$. However, the second independent matching condition is not determined by fundamental conservation laws and remains unconstrained within this framework. In this regard, it is interesting to note that the RTA framework is recovered if this second condition is fixed as $n=n_{0}$, or equivalently, $\epsilon=\epsilon_{0}$, demonstrating that the RTA prescription is a special case of the BGK framework. We shall see later that the scalar freedom due to the unconstrained second matching condition affects the coefficient of bulk viscosity, which is the transport coefficient corresponding to scalar dissipation in the system.

So far, we have considered the case of non vanishing chemical potential. On the other hand, the equilibrium net number density, defined in Eq. (7), vanishes in the limit of zero chemical potential. This implies that the BGK collision term in Eq. (16) is ill defined in this limit. However it is important to note that the limit of vanishing chemical potential is relevant for ultra-relativistic heavy-ion collisions. Therefore, it is desirable to modify the BGK collision kernel in order to extend its regime of applicability. At this juncture, we are well equipped to propose a modification to BGK collision kernel that is well-defined for all values of chemical potential. To this end, we rewrite the condition necessary for energy-momentum conservation from BGK collision kernel, Eq. (18), in the form

Substituting the above equation in Eq. (16), we obtain Boltzmann equation for particles with a modified BGK (MBGK) collision kernel,

and similarly for anti-particles with $f\to\bar{f}$ and $f_{0}\to\bar{f}_{0}$. The advantage of the above modification is that the collision kernel conserves energy-momentum by construction and is applicable to systems even without any conserved four-current, i.e., in the limit of vanishing chemical potential. Additionally, MBGK is uniquely defined even for a system with multiple conserved charges, as opposed to the BGK collision kernel. In the case of finite chemical potential, the matching condition, Eq. (18), ensures net particle four-current conservation. It is important to note that BGK and MBGK are completely equivalent for the purpose of the derivation of hydrodynamic equations at finite chemical potential. Therefore, similar to the BGK framework, we require a second matching condition, along with the matching condition given in Eq. (19), to formulate dissipative hydrodynamics using MBGK description. In the following, we consider the MBGK Boltzmann equation, Eq. (21), to obtain non-equilibrium correction to the distribution function.

In order to obtain the non-equilibrium correction to the distribution function, we use Eq. (6) to rewrite the MBGK Boltzmann equation, Eq. (21), as

and similarly for anti-particles. The next step is to solve the above equation, order-by-order in gradients. In this work, we intend to obtain the non-equilibrium correction to the distribution function up to first-order in derivative, which we represent by $\delta f_{1}$. However, obtaining the expressions for $\delta f_{1}$ from Eq. (22) is not straightforward because it contains $\delta\epsilon$ which is defined in Eq. (11) as an integral over $\delta f$. Therefore, to solve for $\delta f_{1}$, we examine each term individually. At first, we use Eqs. (7)-(9) to substitute equilibrium hydrodynamic quantities in the conservation Eqs. (3)-(5), and obtain, up to first-order in gradients,

where,

The thermodynamic integrals are given by,

and we identify $n_{0}=I^{-}_{10}$, $\epsilon_{0}=I^{+}_{20}$ and $P_{0}=I^{+}_{21}$. Replacing $\dot{\alpha},\dot{\beta},\nabla_{\mu}\beta$ from Eq. (23), the left-hand side of Eq. (22) becomes,

where,

Similar procedure can be followed for anti-particles. We assume $\delta f_{1}$ to have the same form as in Eq. (26),

and similarly for anti-particles. In the above expression, the coefficients $B_{\Pi}$, $B_{n}$ and $B_{\pi}$ needs to be determined using Eq. (22), up to first order in derivatives. At this point, we note that some discussion is due on the choice of the structure of $\delta f_{1}$. In general, we may write $\delta f_{1}=\phi_{1}f_{0}$, with,

where $G_{\mu_{1}\cdots\mu_{\ell}}$ is some $\ell$-rank gradient [45], $A^{\langle\mu_{1}}\cdots A^{\mu_{\ell}\rangle}=\Delta^{\mu_{1}\cdots\mu_{\ell}}_{\nu_{1}\cdots\nu_{\ell}}A^{\nu_{1}}\cdots A^{\nu_{\ell}}$, and $B_{\ell}$ is the coefficient of the $\ell$-th rank gradient, which can be function of particle momenta. Substituting Eq. (30) in Eq. (22), and comparing both sides of the Boltzmann equation, one may conclude $B_{\ell}=0$ for $\ell\geq 3$, where we also have to recall the scalar nature of $\delta\epsilon$.

It should be noted here that the modified collision kernel in Eq. (22) possess a term proportional to $\delta\epsilon$ which also needs to be determined up to first order in gradient. For that purpose we substitute the above decomposition of $\delta f_{1}$ and the corresponding decomposition for anti-particles in the definition of $\delta\epsilon$ given in Eq. (11). It can be shown that the contributions from both the vector and tensor components drop out due to their orthogonality properties. Only the contribution from the scalar gradient survives and is given by

Using Eqs. (26), (29) and (31) into Eq. (22) and comparing both sides, we get

Another set of equations in terms of $\bar{A}_{\Pi}$, $\bar{A}_{n}$ and $\bar{A}_{\pi}$ can be obtained by considering MBGK equation, analogous to Eq. (22), for anti-particles. Here, coefficients $B_{n}$, $\bar{B}_{n}$, $B_{\pi}$ and $\bar{B}_{\pi}$ are easily determined but $B_{\Pi}$ and $\bar{B}_{\Pi}$ require further investigation.

To obtain their expressions, we consider $B_{\Pi}$ to be of the general form, $B_{\Pi}=\sum_{k=-\infty}^{+\infty}b_{k}\left(u\cdot p\right)^{k}$ and $\bar{B}_{\Pi}=\sum_{k=-\infty}^{+\infty}\bar{b}_{k}\left(u\cdot p\right)^{k}$. Substituting these in Eq. (32) and its corresponding equation for anti-particles, we can conclude that the only non-zero $b_{k}$ and $\bar{b}_{k}$ are the ones with $k=-1,0,1$. We obtain

Substituting Eqs. (27) and (34) in Eq. (32), we find

where we have also used the relation analogous to Eq. (32) for anti-particles. On the other hand, for $b_{0}$ and $\bar{b}_{0}$, we find two coupled equations which are identical and leads to the relation

Hence, we see that a unique solution for $b_{0}$ and $\bar{b}_{0}$ can not be obtained but they are constrained by the above relation. We need to provide one more condition, which we recognize as the second matching condition, to fix $b_{0}$ and $\bar{b}_{0}$ separately.

Nevertheless, at this stage, we can determine $\delta f_{1}$ and $\delta\bar{f}_{1}$ up to a free parameter, $b_{0}$, by using Eqs. (32)-(36) into Eq. (29), and similarly for anti-particles. We obtain,

Note that for vanishing chemical potential, we have $\alpha=\chi_{a}=0$. In this case, Eqs. (37) and (38) coincide to give

where we have used $\chi_{b}=\beta c_{s}^{2}$, with $c_{s}^{2}$ being the squared of the speed of sound, given by,

Here $z\equiv m/T$ is the ratio of particle mass to temperature.

The first-order correction to the phase-space distribution functions of the particles and anti-particles at finite $\mu$ are given by Eqs. (37) and (38). Substituting them in Eqs. (10)-(14), we obtain the first-order expressions for non-equilibrium hydrodynamic quantities as

where,

Note that the parameter $b_{0}$ appears in the expressions of $\nu$, $e$ and $\rho$. Of these, $\nu$ and $e$ vanishes for $b_{0}=-\chi_{a}$ which corresponds to the Landau matching condition and RTA collision kernel. Conductivity $\kappa$ and the coefficient of shear viscosity $\eta$ does not contain the parameter $b_{0}$, and the expressions for these two transport coefficients, given in Eq. (44), matches with those derived using RTA collision kernel [11]. Next, we analyze entropy production in the MBGK setup in order to identify dissipative transport coefficients in Eqs. (42)-(44).

To study entropy production, we start from the kinetic theory definition of entropy four-current, given by the Boltzmann’s H-theorem, for a classical system

The entropy production is determined by taking four-divergence of the above equation,

Using the MBGK Boltzmann equation, i.e., Eq. (22), and keeping terms till quadratic order in deviation-from-equilibrium, we obtain

where we have defined $\phi\equiv\delta f/f_{0}$ and $\bar{\phi}\equiv\delta\bar{f}/\bar{f}_{0}$.

Using Eqs. (37) and (38) in Eq. (47), we obtain,

where,

Note that the right-hand-side of Eq. (48) represents entropy production due to dissipation in the system. Here the shear stress tensor $\pi^{\mu\nu}$ is the tensor dissipation, the particle diffusion four-current $n^{\mu}$ is the vector dissipation and $\Pi$ is the scalar dissipation, referred to as the bulk viscous pressure²²2We can further identify that, $\left(\frac{\partial P}{\partial\epsilon}\right)_{n}=\frac{\chi_{b}}{\beta}$ and, $\left(\frac{\partial P}{\partial n}\right)_{\epsilon}=-\frac{\chi_{\alpha}}{\beta}$.. From Eq. (49), we observe that $\delta P$, $\delta\epsilon$, and $\delta n$, all contribute to the bulk viscous pressure. Comparing with the Navier-Stokes relation $\Pi=-\zeta\,\theta$, the coefficient of bulk viscosity becomes,

Demanding that Eq. (48) does not violate the second law of thermodynamics, i.e., $\partial_{\mu}S^{\mu}\geq 0$, leads to the following constraints  [46],

These three transport coefficients represent the three dissipative transport phenomena of the system related to the transport of momentum and charge. We see that out of the three transport coefficients, only $\zeta$ depends on the parameter $b_{0}$ and the second matching condition is necessary to uniquely determine $\zeta$. This is to be expected because the matching conditions are scalar conditions and should only affect the scalar dissipation in the system, i.e., bulk viscosity. In the following, we specify the second matching condition.

With the parameter, $b_{0}$ still not specified, the hydrodynamic equations obtained using the MBGK Boltzmann equation forms a class of hydrodynamic theories. A specific hydrodynamic theory is determined by a specific $b_{0}$ parameter. We can access different hydrodynamic theories by varying the $b_{0}$ parameter, which is solely controlled by the second matching condition. Thus, picking a specific second matching condition will fix $b_{0}$ and hence the hydrodynamic theory. To this end, we define a function $\mathcal{A}_{r}^{\pm}$ as [20, 47],

The second matching condition then amounts to assigning a value for a given $\mathcal{A}_{r}^{\pm}$. For instance, the RTA matching conditions can be recovered by setting $\mathcal{A}_{1}^{-}=\mathcal{A}_{2}^{+}=0$. It is apparent that the choice of a second matching condition is vast, and determination of the full list of the allowed ones is a non-trivial task that goes beyond the scope of the present work. Presently, for the second matching condition, we shall restrict our analysis to a special set $\mathcal{A}_{r}^{+}=0$. These matching conditions ensures that the homogeneous part of $\delta f$ vanishes³³3It must be noted that this is not the only class of matching conditions that guarantee the zero value of the homogeneous part. [48] and are also valid in the zero chemical potential limit. Using Eqs. (37) and (38) in our proposed matching condition $\mathcal{A}_{r}^{+}=0$, we obtain

We emphasize that the above equation along with Eq. (18) constitute the two matching conditions required to define temperature and chemical potential of the system. In order to compare with the usual Landau matching conditions, we express Eq. (6) using first order results obtained in Eqs. (41) and (42),

The above conditions reduce to the usual Landau matching condition, $n=n_{0}$ and $\epsilon=\epsilon_{0}$, for the particular choice of $b_{0}=-\chi_{a}$. In the next section, we explore the effect of different $b_{0}$ on the coefficient of bulk viscosity.

In this section, we study the effect of MBGK collision kernel on transport coefficients. In the previous Section, we found that the effect of MBGK collision kernel manifests in the parameter $b_{0}$ which affects only the scalar dissipation, namely bulk viscous pressure. On the other hand, the vector (net particle diffusion) and tensor (shear stress tensor) dissipation remain unaffected. Therefore, we study only the properties of bulk viscous coefficient in this section.

Before we proceed to quantify the effect of varying the second matching condition on the coefficient of bulk viscosity, we must establish the allowed values for the parameter $b_{0}$. To this end, we note that the second law of thermodynamics demands that the coefficient of bulk viscosity must be positive, Eq. (51). In Fig. 1, we plot $b_{0}$ vs $z$ for different values of $r$ required to define the second matching condition in Eq. (53), at zero chemical potential. The red region in Fig. 1 corresponds to the part of $b_{0}$-$z$ plane where the coefficient of bulk viscosity becomes negative. Therefore all values of $r$ for which the curves for $b_{0}$ lies in the red zone are not physical and must be discarded. The boundary of the red region corresponds to the $\zeta=0$ line and is

Note that $b_{0}^{l}$ as given in Eq. (56) is defined in terms of the ratio of thermodynamic quantities only. In order to ensure $\zeta\geq 0$, we can consider $b_{0}\geq b_{0}^{l}(\alpha,z)$ as a physical constraint on the choice of $b_{0}$. While one can use different parametrizations for $b_{0}$, this condition is more rigorous and independent of the choice of parametrization. With the current choice, we find that $b_{0}$ parameter with non-negative values of $r$ respects the requirement of the second law of thermodynamics, Eq. (51), whereas, large negative values of $r$, violate the constraint, $b_{0}\geq b_{0}^{l}$ leading to negative values of $\zeta$. In Fig. 1, we see that the curve for $b_{0}$, which corresponds to $r=-4$, passes through the physically forbidden region. The black line with $r=2$ represents the $b_{0}$ for which the MBGK reduces to the RTA, where $b_{0}$ vanishes for all $z$ .

Having determined the allowed range of $r$ and equivalently, the allowed values of $b_{0}$, we will restrict ourselves to $b_{0}$ corresponding to $r\geq 0$ values. In Fig. 2 we plot the dimensionless quantity $\zeta/\left(s_{0}\tau_{\rm R}T\right)$ for MBGK with $r=0$, and RTA ($r=2$) against $T/m$ for different values of chemical potential, where $s_{0}\equiv(\epsilon_{0}+P_{0}-\mu\,n_{0})/T$. We observe that $\zeta/\left(s_{0}\tau_{\rm R}T\right)$ is a non-monotonous function of temperature, having a maximum for each $r$ for MBGK case, similar to the behavior known from RTA [20, 47, 49]. We also note that the dependence of $\zeta/\left(s_{0}\tau_{\rm R}T\right)$ on $\alpha$ is also non-monotonous, which can be realized by observing that not only the position of the peak for $\alpha=1$ is at higher $T/m$ values than for $\alpha=0$ and $\alpha=2.5$, but the peak value is also higher for $\alpha=1$ compared to $\alpha=0$ and $\alpha=2.5$.

To better understand the effect of changing matching conditions on the behavior of the bulk viscosity for the MBGK collision kernel, we focus on the zero chemical potential limit. In this limit, we study the scaling behavior of the ratio of the coefficient of bulk viscosity to shear viscosity, $\zeta/\eta$, with conformality measure $1/3-c_{s}^{2}$. In Fig. 3, we plot the ratio $(\zeta/\eta)/(1/3-c_{s}^{2})^{2}$ as a function of $z$ for different $r$ values. We observe that this ratio saturates in both small-$z$ and large-$z$ limits indicating a squared dependence of $\zeta/\eta$ on the conformality measure, characteristic to weakly coupled systems. We also observe that in the small-$z$ limit, this ratio saturates to different values whereas in the large-$z$ limit, they all converge. In order to better understand the behavior of $\zeta/\eta$ in these regimes, we separately analyze the small-$z$ and large-$z$ limits.