Bulk viscosity in the nonlinear and anharmonic regime of strange quark matter

The bulk viscosity of cold, dense three-flavor quark matter is studied as a function of temperature and the amplitude of density oscillations. The study is also extended to the case of two different types of anharmonic oscillations of density. We point several qualitative effects due to the anharmonicity, although quantitatively they appear to be relatively small. We also find that, in most regions of the parameter space, with the exception of the case of a very large amplitude of density oscillations (i.e. 10% and above), nonlinear effects and anharmonicity have a small effect on the interplay of the nonleptonic and semileptonic processes in the bulk viscosity.


Introduction
Cold, dense quark matter is one of possible states of baryonic matter formed at very high densities. While there is little doubt that such quark matter can be formed in principle, the value of the critical density needed remains unknown. This is one of the reasons why it is not settled whether such form of matter can exist inside neutron stars. Indeed, the interior region of neutron stars is the most likely place to find dense quark matter. The corresponding densities reach up to about 10 times the nuclear saturation density, while the temperatures remain moderately low, of the order of 1 MeV or less.
A recent report on a precision measurement of the mass of the binary millisecond pulsar J1614-2230 [1] seems to strongly constraint the possibility of stellar dense quark matter [2]. Nevertheless, it is hard to rule out quark matter solely based on the measurement of such a global property of a star such as its mass.. More informative probes of the state of matter at the highest densities are likely to be based on stellar characteristics determined by the transport properties of matter, various emission rates, as well as certain thermodynamic properties, which are sensitive to the spectrum of quasiparticles at the Fermi surface and which, therefore, can provide a deeper insight into the microscopic nature of dense matter.
In this paper, we study one of such transport characteristics of dense quark matter, the bulk viscosity. In general, viscosity is one of possible mechanisms responsible for damping of the so-called r-mode instabilities in compact stars [3,4,5,6,7]. The emission of gravitational waves tends to drive a differential collective motion of the stellar fluid in the form of r-modes [8,9,10,11]. Such modes in turn increase the stellar quadrupole moment and feed back to result in even stronger gravitational emission. If not damped by dissipative processes, the resonant growth of the r-modes can eventually lead to a breakup of the star.
The bulk viscosity in the normal phase of three-flavor quark matter is usually dominated by nonleptonic weak processes [3,7,12,13,14,15,16,17,18], shown in Figs. 1 (a) and (b). Under certain conditions, semileptonic processes, see Figs. 1 (c)-(f), can lead to an order of magnitude increase of the viscosity as a result of their subtle interplay with the nonleptonic processes [19,20].
When the temperature is sufficiently low or when the magnitude of the density oscillation is sufficiently large, the rates of the weak processes may get substantial nonlinear dependence on the parameter δµ i /T , where δµ i are the chemical potential imbalances that measure the departure of quark matter from β equilibrium. In the context of the nonleptonic processes in strange quark matter, such effects have already been studied [3,15,18] (see also Ref. [21]). Here we extend the analysis to include also the semileptonic processes with the nonlinear corrections in the rates.
When the nonlinear regime is realized, it may be reasonable to expect additional complications in the dynamics responsible for the energy dissipation of the hydrodynamic flow. In this paper we address one of such complications that is associated with the anharmonicity of the density oscillations driven by the r-modes. To the best of our knowledge, previously this has not been addressed in the literature. In studies of the bulk viscosity of stellar matter, one commonly assumes that the oscillations are perfectly harmonic. Of course, this is justified in the linear regime, when different harmonics do not interfere.
In general, large magnitude density oscillations can hardly remain perfectly harmonic because higher harmonics can be generated by the nonlinearity. As discussed in Refs. [22,23,24,25,26,27], the nonlinear regime in compact stars may be responsible for several qualitatively new features in the dynamics, e.g., coupling of eigenmodes, a non-sinusoidal shape of oscillations, various resonance phenomena, and the formation of shocks and turbulence. At the fundamental level, the nonlinear effects are the consequence of the equation of state of matter and the Einstein equations, e.g., see Ref. [27]. Usually, they play a profound role when the interaction energy is comparable to the energies of modes. However, the nonlinearity is present in principle even at arbitrarily small amplitudes of oscillations. Then, one of the natural questions about the nonlinear regime is: How does the associated non-sinusoidal shape of oscillations (i.e. the anharmonicity due to mode coupling) change the bulk viscosity? It is the goal of this paper to make a first attempt in exploring this issue. Our study will be limited to the case when the density fluctuations are much smaller than the equilibrium density, although they can be much larger than the temperature. In view of this limitation, we cannot address many features of dynamics that may develop in the extreme nonlinear regime. Also, instead of considering realistic nonlinear oscillations found in stellar simulations [22,23,24,25,26,27], we will simply model the shape of density oscillations by using classical solutions to an anharmonic oscillator problem. Such a toy-model study will hopefully lead to a better understanding of the problem and will trigger further developments in the future. The rest of the paper is organized as follows. In the next section, we present the general formalism for calculating the bulk viscosity in dense quark matter with multiple active weak processes and anharmonicity included. In Sec. 3, we start our study of the bulk viscosity in the harmonic regime of density oscillations. Our results qualitatively reproduce the earlier results of Refs. [3,18] and extend them to include additional effects due to semileptonic weak processes. In the same section, we also comment on the linear regime when the expression for the bulk viscosity can be obtained analytically. The anharmonic regime is studied in Sec. 4 for two different types of density oscillations, modeled by solutions to a classical anharmonic oscillator with cubic and quartic potentials, respectively. A brief discussion of our result is given in Sec. 5. General periodic solutions in the case of classical anharmonic oscillators with two different types of potentials are presented in analytical form in Appendix A.

Formalism
One can calculate the bulk viscosity ζ under conditions realized in stars by comparing the hydrodynamic relation for the energy dissipation, averaged over one period τ , with the thermodynamic relation for the mechanical work, counteracting the hydrodynamic flow, The latter is given in terms of the instantaneous pressure P and the specific volume V ≡ 1/n. In the above expressions, n 0 is the equilibrium density and δn = n − n 0 is the density deviation from the equilibrium value. (Even when the magnitude of density oscillations is not vanishingly small, we will assume that |δn| ≪ n 0 .) If the magnitude of the density oscillations δn 0 is vanishingly small, one may simulate the collective motion as a harmonic oscillation, δn(t) = δn 0 cos(ωt). However, in a resonance regime, when large density oscillations develop, nonlinear effects may start to play an important role. In this paper, we study this possibility by simulating two different types of anharmonic density oscillations. The two types correspond to oscillators with cubic and quartic terms in the potential energy, which have different symmetry properties under δn → −δn. The corresponding equations of motion read: δn + ω 2 0 δn 1 + βδn 2 = 0, (Type II).
Note that the coupling constants α and β have the dimensions of an inverse density and an inverse density squared, respectively. It is convenient, therefore, to introduce the dimensionless parameters α * ≡ αδn 0 and β * ≡ β(δn 0 ) 2 , which are given in terms of the amplitude of density oscillations δn 0 . [For asymmetric oscillations, described by Eq. (3), we assume that the amplitude is the maximum deviation from the equilibrium point.] Note that the parameters α * and β * can be either positive or negative. General periodic solutions for both types of anharmonic oscillators can be given in terms of the Jacobi elliptic functions. The corresponding solutions are presented in Appendix A. By substituting these exact solutions into Eq. (1) and making use of the result in Eq. (A.4), we derive where the constant F for each type of solution is determined in terms of α * and β * , see Eqs. (A.5) and (A.9), respectively. As is easy to check, F → 1 in the harmonic limit α * → 0 (Type I) or β * → 0 (Type II). By comparing Eqs. (2) and (5), we obtain the following expression for the bulk viscosity: When there is a departure from β equilibrium, the pressure can be given in terms of the instantaneous composition, where X e ≡ n e /n and X s ≡ n s /n is the electron and strangeness fractions, andP is the pressure in equilibrium. The susceptibility functions C 1 and C 2 were defined in Ref. [19]. By taking into account that δn is a periodic function, one finds that only the last two terms in the pressure (7) contribute to the bulk viscosity (6), The instantaneous composition of quark matter is determined by the weak processes, shown in Fig. 1. Taking all of them into account, we derive the following set of nonlinear differential equations for the electron and strangeness fractions: where δµ 1 ≡ µ s − µ d , δµ 2 ≡ µ s − µ u − µ e , and the notation for the λ-rates are the same as in Refs. [19,20], In Eqs. (9) and (10), all higher order corrections in powers of δµ i /T were taken into account, while higher order corrections in powers of δµ i /µ were neglected. This is the same approximation that was used in Ref. [18]. By calculating the semileptonic rates using the approach of Ref. [28], it is easy to check that the coefficients of the nonlinear terms are the same as in the nucleon direct Urca process [29,30]: χ 0 = 1, χ 1 = 10 17π 2 , and χ 2 = 1 17π 4 . The corresponding coefficients in the nonleptonic rates are Υ 0 = 1 and Υ 1 = 1 4π 2 [3,18,31]. The functions that describe the deviation from equilibrium, δµ i , can be equivalently rewritten in terms of the electron and strangeness fractions: δµ i = C i δn n +B i δX e +A i δX s , where coefficient functions A i , B i and C i were defined in Ref. [19]. By making use of these relations and Eqs. (9) and (10), we derive the following self-consistent set of equations for the dimensionless quantities In this study, δn is a periodic function that describes anharmonic oscillations of either Type I or Type II, see Eqs. (3) and (4). We make use of the analytical results in Appendix A and solve Eq. (14) numerically. When the solutions for ν i (with i = 1, 2) are available, one can invert the relations for δµ i in terms of δn, δX e and δX s in order to determine the deviation of the electron and strangeness fractions, where Finally, by making use of these results in Eq. (8), we can calculate the bulk viscosity,

Harmonic oscillations
In the limiting case of harmonic oscillations, the density deviations are described by δn = δn 0 sin(ω 0 t) and the constant F in Eq. (5) is equal to 1. The numerical results for the bulk viscosity as a function of δn 0 /n 0 for several fixed values of temperature are shown in Fig. 2. The linear regime corresponds to small values of δn 0 /n 0 , where the bulk viscosity saturates. It is also interesting to present the temperature dependence of the bulk viscosity. The corresponding plots for several fixed values of δn 0 /n 0 are shown in Fig. 3. As we see, with decreasing the temperature, the bulk viscosity eventually levels off. This is the outcome of reaching the nonlinear regime, and it would be absent if the linear approximation were used instead. In both Fig. 2 and Fig. 3, the results for two representative values of the period of oscillations, τ = 0.1 s and τ = 10 −3 s, are shown. (In all calculations, we used the model parameters from Ref. [19] at n = 5ρ 0 , where ρ 0 ≃ 0.15 fm −3 is the nuclear saturation density.) In most regions of the parameter space, our results qualitatively agree with earlier findings in Refs. [3,12,13,14,15,16,17,18]. The only notable difference occurs in Fig. 3 around the "semileptonic" hump (T ∼ 1 MeV). This comes from the interplay of semileptonic weak processes with the more dominant nonleptonic ones [19,20]. In most of the studies, this is neglected because of a smallness of the semileptonic rates.
Here it is appropriate to mention that, in application to stellar quark matter, the bulk viscosity may not be the only, or even the dominant mechanism responsible for damping of the r-mode instabilities. For example, at sufficiently low temperatures, the corresponding dissipative dynamics is known to be dominated by the shear viscosity [7,32]. (For several representative studies of the shear viscosity in dense quark matter see, for example, Refs. [33,34,35].) In this connection, our low-temperature results in Figs. 2 and 3 should be used only for the purpose of determining where exactly the damping by the shear viscosity takes over.   In order to better understand the role of nonlinear terms in the weak rates, see Eqs. (9) and (10), as well as their effect on the bulk viscosity, it is instructive to study the linear approximation. In this case, the expression for the viscosity can be derived analytically. As we shall see, this will be also helpful to elucidate the role of induced oscillations of δµ 2 = µ s − µ u − µ e when the semileptonic processes are formally switched off.

Harmonic oscillations: linear approximation
In the linear regime, Eqs. (9) and (10) for dimensionless functions ν i simplify down to where ϑ = 2πt/τ is the dimensionless time variable, d i = C i δn 0 /(T n 0 ) is the magnitude of the "driving force", and the other coefficient functions are The general solution to this set of equations in the steady state regime is given by where the coefficients satisfy the following set of algebraic equations: It is straightforward, although tedious to solve this set of equations. When the solution is available, the result for the bulk viscosity will follow from Eq. (20), i.e.
It can be shown that this expression (with the appropriate solutions for x 1 and x 2 ) coincides exactly with the result for the bulk viscosity, obtained in Ref. [19].

Harmonic oscillations: nonleptonic contribution in linear approximation
If one ignores the semileptonic processes (i.e. if one formally takes λ 2 = λ 3 = 0), the linearized equations (21) and (22) for ν i 's take the following form: The explicit solution to this set of equations in the steady state regime is given by and the corresponding expression for the bulk viscosity reads where we used the relation f 1 /h 1 = A 1 /A 2 to arrive at the final result. As expected, this agrees with the known result [3,13,18,19]. It is interesting to notice that the final result for the bulk viscosity receives a nonzero contribution due to the oscillation of ν 2 ≡ δµ 2 /T . Since δµ 2 ≡ µ s − µ u − µ e controls the imbalance of the rates in the semileptonic processes, shown in Figs. 1 (c) and (d), which are formally switched off in the approximation at hand, one might wonder why there should be such a contribution at all. The answer is quite simple. In absence of the semileptonic processes the electron fraction in quark matter cannot change. However, when the nonleptonic processes drive the oscillations of the strangeness composition, they inevitably induce the oscillations of δµ 2 ≡ µ s −µ u −µ e . Then, the latter contributes to the instantaneous pressure and, in turn, to the bulk viscosity. Interestingly, such a contribution due to the induced oscillation of ν 2 ≡ δµ 2 /T were ignored in all previous studies [3,12,13,14,15,16,17,18]. Fortunately, the corresponding correction is quantitatively small. The reason for its smallness seems to be rooted in the "accidental" fact that one of the susceptibility functions, B 2 , is inversely proportional to the square of the chemical potential of electrons (rather than quarks) and, thus, is considerably larger than all others [19].

Bulk viscosity in anharmonic regime
In this section we study the effect that anharmonic oscillations of the density have on the bulk viscosity of dense quark matter.

Anharmonic oscillations of Type I
Let us start by modeling the density oscillations of quark matter δn(t) by a time dependent anharmonic function of Type I, which is a solution to Eq. (3) with a fixed anharmonicity parameter α * . Before we proceed to the numerical results, it is important to notice that the corresponding oscillations are asymmetric with respect to the equilibrium point δn eq = 0. For α * < 0, the density oscillations are larger in the direction of positive δn, while for α * > 0, the oscillations are larger in the direction of negative δn, see also Fig. A1. The cases of the positive and negative anharmonicity parameters are physically equivalent, however. Indeed, they are related by the following sign reversal symmetry: α → −α and δn → −δn. Therefore, it is sufficient to study only one of them. For technical reasons, we choose α * < 0.
Typical results for the bulk viscosity as a function of anharmonicity parameter α * are shown in the left panel of Fig. 4 for the whole range of negative α * , i.e. −1 < α * < 0, for which physically meaningful periodic solutions exist. We used the following values of the period and the amplitude of density oscillations: τ = 10 −3 s and δn 0 /n 0 = 10 −3 , and plotted the results for several representative values of temperature. [It should be emphasized that the period τ is related to the "bare" frequency ω 0 by a modified relation, see Eq. (A.3).] In general, we find that the bulk viscosity decreases with increasing the degree of anharmonicity. This qualitative behavior may be understood as the result of an effective increase of the frequency of oscillations due to an admixture of higher harmonics. Quantitatively, however, the effect is rather small. Only a very large anharmonicity (α * ≈ −1) leads to a substantial decrease of the bulk viscosity.

Anharmonic oscillations of Type II
Anharmonic density oscillations of Type II are modeled by a function δn(t), which is a solution to Eq. (4) with a fixed anharmonicity parameter β * . Conceptually, this is a simpler case because the oscillations are symmetric about the equilibrium point δn eq = 0. Unlike the case of Type I oscillations, there is no reversal symmetry here. As in the previous case, however, periodic solution exist only for a range of values of the the anharmonicity parameter, β * > −1.
Numerical results for the bulk viscosity as a function of anharmonicity parameter β * are shown in the right panel of Fig. 4. The qualitative dependence of the viscosity on the parameter β * is somewhat different. While it decreases at large values of parameter β * , there is a range of small negative values of β * , where it slightly grows with increasing anharmonicity. Moreover, this feature seems to be rather general and especially pronounced in the nonlinear regime (small temperature). Just like in the case of Type I oscillations, the effects appear to be rather small.

Discussion
In this paper we studied the bulk viscosity of dense quark matter by taking into account the nonlinear dependence of the nonleptonic and semileptonic weak rates on the parameter δµ i /T , where δµ i are the chemical potential imbalances that control the departure of strange quark matter from β equilibrium. We reproduce the earlier obtained interplay of the nonleptonic and semileptonic processes, leading to an increase ("hump") of the viscosity in a narrow temperature range around 1 MeV. The nonlinear corrections have a small effect on the corresponding shape of the "hump". The reason for this is a relatively high temperature (T ∼ 1 MeV), at which the corresponding effects can occur. At such moderately high temperatures, the interplay between the two types of weak processes is substantially affected only if the nonlinearity (measured by δn 0 /n 0 ) is well above 10%.
In this study, we also found that the anharmonicity of density oscillations has an effect on the bulk viscosity, even though the effect was not large in the cases that we studied. For a strong anharmonicity, the bulk viscosity showed a substantial decrease. We also saw that different types of anharmonicity have slightly different qualitative as well as quantitative outcomes. This finding may suggest that some types of anharmonicity may be more efficient and, thus, lead to larger corrections to the bulk viscosity. In this preliminary study, we considered only a toy model to introduce the anharmonic oscillations. In the future, it may be interesting to study more realistic types of density oscillations that result from the actual nonlinear dynamics of stellar r-modes, produced by the gravitational emission.
Here we studied the bulk viscosity only in the normal phase of dense quark matter. One may wonder, however, how the results will be modified if quark matter is a color superconductor [36]. As suggested by several existing studies in the linear regime, color superconductivity can have a large effect on the bulk viscosity [20,28,37,38]. The additional effects due to the nonlinear regime may be much harder to predict. One of the complications comes from the fact that the weak rates in color superconducting matter have a strong dependence on the ratio of the superconducting energy gap and temperature, φ/T . Moreover, in the suprathermal regime, their dependence on δµ i /T will most likely be nonpolynomial. These technical difficulties can be resolved in principle, but the outcome is not obvious and the corresponding study is outside the scope of the present paper. (The apparent asymmetry between positive and negative values of αx 0 is a result of the assumption that x 0 is the maximum deviation from the equilibrium point in the positive x-direction.) It is given in terms of the Jacobi elliptic function sn (u|m) as follows: where a ± = 1 2 1 + 2α * ± (1 + 2 3 α * )(1 − 2α * ) and α * ≡ αx 0 is a dimensionless parameter that measures the maximum deviation of the solution from the harmonic regime. Note that x(t) is periodic with the period given by where K (a − /a + ) is the complete elliptic integral of the first kind. Two representative solutions are shown in Fig. A1. For anharmonic solutions of this type, the average kinetic energy is given by where the analytical expression for constant F reads Here E (a − /a + ) is the complete elliptic integral of the second kind. is described by the following equation of motion: x + ω 2 0 x 1 + βx 2 = 0, (Type II). (A.6) Using the known parametric solution to the above differential equation [39] and assuming that x 0 is the maximum deviation from the equilibrium point x = 0, we find that the periodic solution exists for β * ≡ β x 2 0 > −1. The solution is given in terms of the Jacobi elliptic function, x(t) = x 0 1 + (1/2)β * √ 1 + β * sd ω 0 (t − t 0 ) 1 + β * β * 2(1 + β * ) .
This is a periodic solution with the period equal Two representative solutions are shown in Fig. A2. For solutions of this type, the average kinetic can be given in the same form as in Eq. (A.4), but the value of the corresponding constant F is different, (A.9)