text:
# Nonlinear Outcome of Gravitational Instability in Disks with Realistic Cooling

## 1 Introduction

The outer regions of accretion disks in both active galactic nuclei (AGN) and young stellar objects (YSO) are close to gravitational instability (for a review see, for AGN: Shlosman et al. 1990; YSOs: Adams & Lin 1993). Gravitational instability can be of central importance in disk evolution. In some disks, it leads to the efficient redistribution of mass and angular momentum (e.g. Larson 1984; Laughlin & Rozyczka 1996; Gammie 2001). In other disks, gravitational instability leads to fragmentation and the formation of bound objects. This may cause the truncation of circumnuclear disks (Goodman, 2003), or the formation of planets (e.g. Boss 1997, and references therein).

We will restrict attention to disks whose potential is dominated by the central object, and whose rotation curve is therefore approximately Keplerian. Gravitational instability to axisymmetric perturbations sets in when the sound speed $`c_s`$, the rotation frequency $`\mathrm{\Omega }`$, and the surface density $`\mathrm{\Sigma }`$ satisfy

$$Q\frac{c_s\mathrm{\Omega }}{\pi G\mathrm{\Sigma }}<Q_{crit}1$$
(1)
(Toomre, 1964; Goldreich & Lynden-Bell, 1965). Here $`Q_{crit}=1`$ for a “razor-thin” (two- dimensional) fluid disk model of the sort we will consider below, and $`Q_{crit}=0.676`$ for a finite-thickness isothermal disk (Goldreich & Lynden-Bell, 1965). 111For global models with radial structure, nonaxisymmetric instabilities typically set in for slightly larger values of $`Q`$ (see Boss 1998 and references therein). The instability condition (1) can be rewritten, for a disk with scale height $`Hc_s/\mathrm{\Omega }`$, around a central object of mass $`M_{}`$,

$$M_{disk}\frac{H}{r}M_{},$$
(2)
where $`M_{disk}=\pi r^2\mathrm{\Sigma }`$. For YSO disks $`H/r0.1`$ and thus a massive disk is required for instability. AGN disks are expected to be much thinner. The instability condition can be rewritten in a third, useful form if we assume that the disk is in a steady state and its evolution is controlled by internal (“viscous”) transport of angular momentum. Then the accretion rate $`\dot{M}=3\pi \alpha c_s^2\mathrm{\Sigma }/\mathrm{\Omega }`$, where $`\alpha 1`$ is the usual dimensionless viscosity of Shakura & Sunyaev (1973), and

$$\dot{M}\frac{3\alpha c_s^3}{G}=7.1\times 10^4\alpha \left(\frac{c_s}{1\mathrm{km}\mathrm{s}^1}\right)^3\mathrm{M}_{}\mathrm{yr}^1$$
(3)
implies gravitational instability (e.g. Shlosman et al. (1990)). Disks dominated by external torques (e.g. a magnetohydrodynamic \[MHD\] wind) can have higher accretion rates (but not arbitrarily higher; see Goodman 2003) while avoiding gravitational instability.

For a young, solar-mass star accreting from a disk with $`\alpha =10^2`$ at $`10^6\mathrm{M}_{}\mathrm{yr}^1`$, equation (3) implies that instability occurs where the temperature drops below $`17\mathrm{K}`$. Disks may not be this cold if the star is located in a warm molecular cloud where the ambient temperature is greater than $`17\mathrm{K}`$, or if the disk is bathed in scattered infrared light from the central star (although there is some evidence for such low temperatures in the solar nebula, e.g. Owen et al. 1999). If the vertically-averaged value of $`\alpha `$ is small and internal dissipation is confined to surface layers, as in the layered accretion model of Gammie (1996), then instability can occur at higher temperatures, although equation (2) still requires that the disk be massive.

AGN disk heating is typically dominated by illumination from a central source. The temperature then depends on the shape of the disk. If the disk is flat or shadowed, however, and transport is dominated by internal torques, one can apply equation (3). For example, in the nucleus of NGC 4258 (Miyoshi et al., 1995) the accretion rate may be as large as $`10^2\mathrm{M}_{}\mathrm{yr}^1`$ (Lasota et al., 1996; Gammie et al., 1999). Equation (3) then implies that instability sets in where $`T<10^4(\alpha /10^2)\mathrm{K}`$. If the disk is illumination-dominated then $`Q`$ fluctuates with the luminosity of the central source.

In a previous paper (Gammie, 2001), one of us investigated the effect of gravitational instability in cooling, gaseous disks in a local model. A simplified cooling function $`\mathrm{\Lambda }`$ was employed in these simulations, with a fixed cooling time $`\tau _{co}`$:

$$\mathrm{\Lambda }=\frac{U}{\tau _{co}},$$
(4)
where $`U`$ the internal energy per unit area. Disk fragmentation was observed for $`\tau _{co}\mathrm{\Omega }3`$. The purpose of this paper is to investigate gravitational instability in a local model with more realistic cooling.

Several recent numerical experiments have included cooling, as opposed to isothermal or adiabatic evolution, and we can ask whether these results are consistent with Gammie (2001). Nelson et al. (2000) studied a global two-dimensional (thin) SPH model in which the vertical density and temperature structure is calculated self-consistently and each particle radiates as a blackbody at the surface of the disk. The initial conditions at a radius corresponding to the minimum initial value of Q ($`1.5`$) for these simulations were $`\mathrm{\Sigma }_o50\mathrm{g}\mathrm{cm}^2,\mathrm{\Omega }8\times 10^{10}\mathrm{s}^1`$; the initial cooling time under these circumstances is $`\tau _{co}250\mathrm{\Omega }^1`$, so fragmentation is not expected and is not observed.

Durisen et al. (2001) consider a global three dimensional (3D) Eulerian hydrodynamics model in which the volumetric cooling rate varies with height above the midplane so as to preserve an isentropic vertical structure. The cooling time is fixed at each radius. Their cooling time $`10\mathrm{\Omega }^1`$ at all radii, so fragmentation is not expected based on the criterion of Gammie (2001). The simulations show structure formation due to gravitational instabilities but not fragmentation.

Rice et al. (2003) consider a global 3D SPH model with a cooling time that is a fixed multiple of $`\mathrm{\Omega }^1(r)`$. They find that their disk fragments when $`\tau _{co}3\mathrm{\Omega }^1`$ and $`M_{disk}=0.1M_{}`$. For a more massive disk ($`M_{disk}=0.25M_{}`$), fragmentation occurred at somewhat higher cooling times ($`\tau _{co}10\mathrm{\Omega }^1`$). This is effectively a global generalization of the local model problem considered by Gammie (2001). The fact that the results are so consistent suggests that the local, thin approximation used in Gammie (2001) and here give a reasonable approximation to a global outcome.

Mayer et al. (2002) consider a global three dimensional SPH model of a circumstellar disk. Explicit cooling is not included, but the equation of state switches from isothermal to adiabatic when gravitational instability begins to set in. This is designed to account for the inefficient cooling of dense, optically thick regions. Fragmentation is observed. Realistic cooling can have a complex influence on disk evolution, and it is not clear that switching between isothermal and adiabatic behavior “brackets” the outcomes that might be obtained when full cooling is used.

Other notable recent work, such as that by Boss (2002), includes strong radiative heating in the sense that the effective temperature of the external radiation field $`T_{irr}`$ is comparable to or larger than the disk midplane temperature $`T_c`$. In the limit that $`T_{irr}T_c`$ we recover the limit considered here and in Gammie (2001); in the limit that $`T_{irr}T_c`$ the disk is effectively isothermal.

The plan of this paper is as follows. In §2 we describe the model, with a detailed description of the cooling function given in §3. The results of numerical experiments are described in §4. Conclusions are given in §5.

## 2 Model

The model we use here is identical to that used in Gammie (2001) in every respect except that we use a more complicated cooling function. To make the description more self-contained, we summarize the basic equations of the model here. The model is local, in the sense that it considers a region of size $`L`$ where $`L/r_o1`$ and $`r_o`$ is a fiducial radius. We use a local Cartesian coordinate system $`xrr_o`$ and $`y(\varphi \mathrm{\Omega }t)r_o`$, where $`r,\varphi `$ are the usual cylindrical coordinates and $`\mathrm{\Omega }`$ is the orbital frequency at $`r_o`$. The model is also thin in the sense that matter is confined entirely to the plane of the disk.

Using the local approximation one can perform a formal expansion of the equations of motion in the small parameter $`L/r_o`$. The resulting equations of motion read, where $`𝐯`$ is the velocity, $`P`$ is the (two-dimensional) pressure, and $`\varphi `$ is the gravitational potential with the time-steady axisymmetric component removed:

$$\frac{D𝐯}{Dt}=\frac{P}{\mathrm{\Sigma }}2𝛀\widehat{𝐞}_z\times 𝐯+3\mathrm{\Omega }^2x\widehat{𝐞}_x\varphi .$$
(5)
For constant pressure and surface density, $`𝐯=\frac{3}{2}\mathrm{\Omega }x\widehat{𝐞}_y`$ is an equilibrium solution to the equations of motion. This linear shear flow is the manifestation of differential rotation in the local model.

The equation of state is

$$P=(\gamma 1)U,$$
(6)
where $`P`$ is the two-dimensional pressure and $`U`$ the two-dimensional internal energy. The two-dimensional (2D) adiabatic index $`\gamma `$ can be mapped to a 3D adiabatic index $`\mathrm{\Gamma }`$ in the low-frequency (static) limit. For a non-self-gravitating disk $`\gamma =(3\mathrm{\Gamma }1)/(\mathrm{\Gamma }+1)`$ (e.g. Goldreich, Goodman, & Narayan 1986; Ostriker, Shu, & Adams 1992). For a strongly self-gravitating disk, one can show that $`\gamma =32/\mathrm{\Gamma }`$. We adopt $`\mathrm{\Gamma }=7/5`$ throughout, which yields $`\gamma =11/7`$.

The internal energy equation is

$$\frac{U}{t}+(U𝐯)=P𝐯\mathrm{\Lambda },$$
(7)
where $`\mathrm{\Lambda }=\mathrm{\Lambda }(\mathrm{\Sigma },U,\mathrm{\Omega })`$ is the cooling function, fully described below. Notice that there is no heating term; heating is due solely to shocks. Numerically, entropy is increased by artificial viscosity in shocks.

The gravitational potential is determined by the razor-thin disk Poisson equation:

$$^2\varphi =4\pi G\mathrm{\Sigma }\delta (z).$$
(8)
For a single Fourier component of the surface density $`\mathrm{\Sigma }_𝐤`$ this has the solution

$$\varphi =\frac{2\pi G}{|𝐤|}\mathrm{\Sigma }_𝐤e^{i𝐤𝐱|kz|}.$$
(9)
A finite-thickness disk has weaker self-gravity, but this does not qualitatively change the dynamics of the disk in linear theory (Goldreich & Lynden-Bell, 1965).

We integrate the governing equations using a self-gravitating hydrodynamics code based on ZEUS (Stone & Norman, 1992). ZEUS is a time-explicit, operator-split, finite-difference method on a staggered mesh. It uses an artificial viscosity to capture shocks. Our implementation has been tested on standard linear and nonlinear problems, such as sound waves and shock tubes. We use the “shearing box” boundary conditions, described in detail by Hawley et al. (1995), and solve the Poisson equation using the Fourier transform method, modified for the shearing box boundary conditions. See Gammie (2001) for further details on boundary conditions, numerical methods and tests.

The numerical model is always integrated in a region of size $`L\times L`$ at a numerical resolution of $`N\times N`$. In linear theory the disk is most responsive at the critical wavelength $`2c_s^2/G\mathrm{\Sigma }_o`$.222The wavelength corresponding to the minimum in the dispersion relation for axisymmetric waves. We have checked the dependence of the outcome on $`L`$ and found that as long as $`L2c_s^2/G\mathrm{\Sigma }_o`$ the outcome does not depend on $`L`$. We have also checked the dependence of the outcome on $`N`$ and found that the outcome is insensitive to $`N`$, at least for the models with $`N256`$ that we use.

## 3 Cooling Function

Our cooling function is determined from a one-zone model for the vertical structure of the disk. The disk cools at a rate per unit area

$$\mathrm{\Lambda }2\sigma T_e^4,$$
(10)
which defines the effective temperature $`T_e`$. The cooling function depends on the heat content of the disk and how that content is transported from the disk interior to the surface: by radiation, convection, or perhaps some more exotic form of turbulent transport such as MHD waves. Low temperature disks are expected to be convectively unstable (e.g. Cameron 1978; Lin & Papaloizou 1980). Cassen (1993) has argued, however, that the radiative heat flux in an adiabatically-stratified disk is comparable to the heat dissipated by turbulence (in an $`\alpha `$-disk model), suggesting that convection is incapable of radically altering the vertical structure of the disk. We will consider only radiative transport.

If the disk is optically thick in the Rosseland mean sense, so that radiative transport can be treated in the diffusion approximation, then (Hubeny, 1990)

$$T_e^4=\frac{8}{3}\frac{T_c^4}{\tau }$$
(11)
where $`\tau `$ is the Rosseland mean optical depth and $`T_c`$ is the central temperature. We will assume that $`T_cT`$, where

$$T=\frac{\mu m_pc_s^2}{\gamma k_B},$$
(12)
and

$$c_s^2=\gamma (\gamma 1)\frac{U}{\mathrm{\Sigma }},$$
(13)
which follows from the equation of state and the assumption that the radiation pressure is small (we have verified that this is never seriously violated). Here $`k_B`$ is Boltzmann’s constant, $`m_p`$ is the proton mass, and $`\mu `$ is the mean mass per particle, which we have set to $`2.4`$ in models with initial temperature below the boundary between the grain-evaporation opacity and molecular opacity and $`\mu =0.6`$ in models with initial temperature above the boundary.

The optical depth is

$$\tau _0^{\mathrm{}}𝑑z\kappa (\rho _z,T_z)\rho _z$$
(14)
where $`\kappa `$ is the Rosseland mean opacity, $`\rho _z`$ and $`T_z`$ are local density and temperature, and $`z`$ is the height above the midplane. Following the usual one-zone approximation,

$$_0^{\mathrm{}}𝑑z\kappa (\rho _z,T_z)\rho _zH\kappa (\overline{\rho },\overline{T})\overline{\rho }$$
(15)
where the overbar indicates a suitable average and $`Hc_s(T)/\mathrm{\Omega }`$ is the disk scale height (we ignore the effects of self-gravity on the disk scale height, which is valid when locally $`Q1`$). Taking $`\overline{T}T`$ and $`\overline{\rho }\mathrm{\Sigma }/(2H)`$ then gives a final, closed expression for $`\mathrm{\Lambda }`$.

We have adopted the analytic approximations to the opacities provided by Bell & Lin (1994). These opacities are dominated by, in order of increasing temperature: grains with ice mantles, grains without ice mantles, molecules, H- scattering, bound-free/free-free absorption and electron scattering. The molecular opacity regime is commonly called the opacity gap; it is too hot for dust, but too cold for H- scattering to contribute much opacity. The opacity can be as much as $`4`$ orders of magnitude smaller than the $`5\mathrm{g}\mathrm{cm}^2`$ typical of the dust-dominated opacity regime. It turns out that this feature plays a significant role in the evolution of gravitationally-unstable disks.

To sum up, the cooling function is

$$\mathrm{\Lambda }(\mathrm{\Sigma },U,\mathrm{\Omega })=\frac{16}{3}\frac{\sigma T^4}{\tau }.$$
(16)
For a power-law opacity of the form $`\kappa =\kappa _0\rho ^aT^b`$, this implies that

$$\mathrm{\Lambda }\mathrm{\Sigma }^{53a/2+b}U^{4+a/2b}.$$
(17)
From this it follows that the cooling time $`\tau _cU/\mathrm{\Lambda }`$ scales as

$$\tau _c\mathrm{\Sigma }^{5+3a/2b}U^{3a/2+b}.$$
(18)
If the disk evolves quasi-adiabatically (as it does if the cooling time is long compared to the dynamical time) then $`U\mathrm{\Sigma }^\gamma `$ and

$$\tau _c\mathrm{\Sigma }^{53\gamma +(a/2)(3\gamma )+b(\gamma 1)}.$$
(19)
Table 1 gives a list of values for this scaling exponent for our nominal value of $`\gamma =11/7`$. Notice that, when ice grains or metal grains are evaporating, and in the bound-free/free-free opacity regime, cooling time decreases as surface density increases.

Our cooling function is valid in the limit of large optical depth ($`\tau 1`$). Since the disk becomes optically thin at some locations in the course of a typical run, we must modify this result so that the cooling rate does not diverge at small optical depth. A modification that produces the correct asymptotic behavior is

$$\mathrm{\Lambda }=\frac{16}{3}\sigma T^4\frac{\tau }{1+\tau ^2}.$$
(20)
This interpolates smoothly between the optically-thick and optically-thin regimes and is proportional to the (Rosseland mean) optical depth in the optically-thin limit. While it would be more physically sensible to use a Planck mean opacity in the optically-thin limit, usually the optically-thin regions contain little mass so their cooling is not energetically significant. An exception is in the opacity gap, where even high density regions become optically thin.

Our simulations begin with $`\mathrm{\Sigma }`$ and $`U`$ constant. The velocity field is perturbed from the equilibrium solution to initiate the gravitational instability. The initial velocities are $`v_x=\delta v_x`$, $`v_y=\frac{3}{2}\mathrm{\Omega }x+\delta v_y`$, where $`\delta 𝐯`$ is a Gaussian random field of amplitude $`\delta v^2/c_s^2=0.1`$. The power spectrum of perturbations is white noise ($`v_k^2k^0`$) in a band in wavenumber $`k_{crit}/4<|k|<4k_{crit}`$ surrounding the minimum $`k_{crit}=1/(\pi Q^2)`$ (with $`G=\mathrm{\Sigma }_o=\mathrm{\Omega }=1`$) in the density-wave dispersion relation. We have checked in particular cases that for $`10^3<\delta v^2/c_s^2<10`$ the outcome is qualitatively unchanged. This is expected because disk perturbations (unlike cosmological perturbations) grow exponentially and the initial conditions are soon forgotten.

Excluding the initial velocity field, the initial conditions for a spatially-uniform disk consist of three parameters: $`\mathrm{\Sigma }_o,U_o`$, and $`\mathrm{\Omega }`$. We fix $`Q=1`$, leaving two degrees of freedom. In models with simple, scale-free cooling functions such as that considered by Gammie (2001), these degrees of freedom remain and can be scaled away by setting $`G=\mathrm{\Sigma }_o=\mathrm{\Omega }=1`$. That is, there is a two-dimensional continuum of disks (with varying values of $`\mathrm{\Sigma }_o`$ and $`\mathrm{\Omega }`$, but the same value of $`Q`$) that are described by a single numerical model.

The opacity contains definite physical scales in density and temperature. The realistic cooling function considered here therefore removes our freedom to rescale the disk surface density and rotation frequency. That is, there is now a one-to-one correspondence between disks with fixed $`\mathrm{\Sigma }`$ and $`\mathrm{\Omega }`$ and our numerical models.

The choice of $`\mathrm{\Sigma }_o`$ and $`\mathrm{\Omega }`$ as labels for the parameter space is not unique. Internally in the code we fix the initial volume density (in $`\mathrm{g}\mathrm{cm}^3`$) and the initial temperature (in Kelvins). These choices are difficult to interpret, however, since they are tied to quantities that change over the course of the simulation; $`\mathrm{\Omega }`$ and the mean value of $`\mathrm{\Sigma }`$ do not.

The cooling is integrated explicitly using a first-order scheme. The timestep is modified to satisfy the Courant condition and to be less than a fixed fraction of the shortest cooling time on the grid. We have varied this fraction and shown that the results are insensitive to it, provided that it is sufficiently small.

## 4 Nonlinear Outcome

### 4.1 Standard Run

Consider the evolution of a single “standard” run, with $`\mathrm{\Sigma }_o=1.4\times 10^5\mathrm{g}\mathrm{cm}^2`$ and $`\mathrm{\Omega }=1.1\times 10^7\mathrm{sec}^1`$. This corresponds to $`T_o=1200`$ and $`\tau _{co}=9.0\times 10^4\mathrm{\Omega }^1`$. The model size is $`L=320G\mathrm{\Sigma }_o/\mathrm{\Omega }^2`$ and numerical resolution $`1024^2`$. The model initially lies at the lower edge of the opacity gap.

The evolution of the kinetic, gravitational and thermal energy per unit area ($`E_k`$, $`E_g`$ and $`E_{th}`$ respectively) normalized to $`G^2\mathrm{\Sigma }_o^3/\mathrm{\Omega }^2`$,333The natural unit that can be formed from $`G`$, $`\mathrm{\Sigma }`$ and $`\mathrm{\Omega }`$. are shown in Figure 1. After the initial phase of gravitational instability the model settles into a statistically-steady, gravito-turbulent state. It does not fragment. Cooling is balanced by shock heating. Energy for driving the shocks is extracted from the shear flow, and the mean shear flow is enforced by the boundary conditions.

The turbulent state transports angular momentum outward via hydrodynamic and gravitational shear stresses. The dimensionless gravitational shear stress is

$$\alpha _{grav}=\frac{1}{\frac{3}{2}\mathrm{\Sigma }c_s^2}_{\mathrm{}}^{\mathrm{}}𝑑z\frac{g_xg_y}{4\pi G}$$
(21)
where $`𝐠`$ is the gravitational acceleration, and the dimensionless hydrodynamic shear stress is

$$\alpha _{hyd}=\frac{\mathrm{\Sigma }v_x\delta v_y}{\frac{3}{2}\mathrm{\Sigma }c_s^2}$$
(22)
where $``$ denote a spatial average. Figure 2 shows the evolution of $`\alpha _{grav}`$ and $`\alpha _{hyd}`$ in the standard run. Averaged over the last $`230\mathrm{\Omega }^1`$ of the run, $`\alpha _{hyd}=0.0079`$, $`\alpha _{grav}=0.017`$, and so the total dimensionless shear stress is $`\alpha =0.025`$, where $``$ denote a space and time average.

The mean stability parameter $`Qc_s\mathrm{\Omega }/\pi G\mathrm{\Sigma }`$ averages $`1.86`$ over the last $`230\mathrm{\Omega }^1`$ of the run. Because the temperature and surface density vary strongly, other methods of averaging $`Q`$ will give different results.

Figure 3 shows a snapshot of the surface density at $`t=50\mathrm{\Omega }^1`$. The structure is similar to that observed in Gammie (2001), with trailing density structures. The density structures are stretched into a trailing configuration by the prevailing shear flow. Their scale is determined by the disk temperature and surface density rather than the size of the box (see Gammie 2001).

### 4.2 Varying $`\mathrm{\Sigma }_o`$ and $`\mathrm{\Omega }`$

We now turn to exploring the two-dimensional parameter space of models. First consider a series of models with the same initial central temperature, but with varying $`\tau _{co}`$. As $`\tau _{co}`$ is lowered the time-averaged gravitational potential energy per unit area $`E_g`$ increases monotonically in magnitude. The gravito-turbulent state becomes more extreme, with larger $`\alpha `$, larger perturbed velocities, and larger density contrasts. Eventually a threshold is crossed and the disk fragments.

Fragmentation is illustrated in Figure 4, which shows a snapshot from a run with $`\mathrm{\Sigma }_o=6.6\times 10^3\mathrm{g}\mathrm{cm}^2`$, $`\mathrm{\Omega }=5.4\times 10^9\mathrm{sec}^1`$. This corresponds to $`T_o=1200`$, $`\tau _{co}=0.025\mathrm{\Omega }^1`$. The run has numerical resolution $`256^2`$ and $`L=80G\mathrm{\Sigma }_o/\mathrm{\Omega }^2`$. The largest bound object in the center of the figure was formed from the collision and coalescence of several smaller bound objects. A snapshot of the optical depth at the same point in the simulation is given in Figure 5. For each snapshot, red indicates high values of the mapped variable and blue indicates low values. Much of the disk is optically thick, but most of the low density regions are optically thin in the Rosseland mean sense.

Lowering $`\tau _{co}`$ sufficiently always leads to fragmentation. We have surveyed the parameter space of $`\mathrm{\Omega }`$ and $`\mathrm{\Sigma }_o`$ to determine where the disk begins to fragment. Each model was run to $`100\mathrm{\Omega }^1`$.444In four cases we had to run the simulation longer to get converged results. Figures 6 and 7 summarize the results. Two heavy solid lines are shown on each diagram. The upper line shows the most rapidly cooling simulations that show no signs of gravitational fragmentation (nonfragmentation point). Quantitatively, we define this as the point at which the time-averaged gravitational potential energy per unit area is equal to $`3G^2\mathrm{\Sigma }_o^3/\mathrm{\Omega }^2`$.555$`3`$ is the potential energy per unit area of a wave at the critical wavelength in a $`Q=1`$ disk with $`\delta \mathrm{\Sigma }/\mathrm{\Sigma }=\sqrt{3}/\pi `$. No bound objects are observed throughout the duration of these runs. The lower line shows the most slowly cooling simulations to show definite fragmentation (fragmentation point). Quantitatively, we define this as the point at which the gravitational potential energy per unit area is equal to $`300G^2\mathrm{\Sigma }_o^3/\mathrm{\Omega }^2`$ at some point during the run.666These runs exhibit bound objects that persist for the duration of the run. Figure 6 shows the data in the $`\rho _o,T_o`$ plane, while Figure 7 shows the results in the $`\mathrm{\Sigma }_o,\mathrm{\Omega }`$ plane. The light contours are lines of constant $`\tau _{co}`$.

The transition from persistent, gravito-turbulent outcomes to fragmentation is gradual and statistical in nature. Figure 8 shows the gravitational potential energy per unit area in the transition region for a series of runs with $`T_o=1200\mathrm{K}`$. The abscissa is labeled with the initial cooling time $`\tau _{co}\mathrm{\Omega }`$. There is a gradual, approximately logarithmic increase in the magnitude of $`E_g`$ as $`\tau _{co}`$ decreases. Runs in this region exhibit the transient formation of small bound objects, which might collapse if additional physics (e.g. the effects of MHD turbulence) were included in the model. Eventually $`E_g`$ begins to increase dramatically, and we define the transition point as the beginning of this steep increase in gravitational binding energy.

Figure 9 shows the run of $`\tau _{co}\mathrm{\Omega }`$ for the fragmentation point, transition point, and nonfragmentation point as a function of $`T_o`$. It is surprising that a disk can begin to exhibit signs of gravitational collapse for $`\tau _{co}\mathrm{\Omega }`$ as large as $`10^6`$, and evade collapse for $`\tau _{co}\mathrm{\Omega }`$ as small as $`0.02`$. A naive application of the results of Gammie (2001) would suggest that fragmentation should occur for $`\tau _{co}\mathrm{\Omega }3`$. Evidently this estimate can be off by orders of magnitude, with the largest error for $`T_o10^3\mathrm{K}`$, just below the opacity gap.

The physical argument for fragmentation at short cooling times is as follows (e.g. Shlosman et al. 1990). Thermal energy is supplied to the disk via shocks. Strong shocks occur when dense clumps collide with one another; this occurs on a dynamical timescale $`\mathrm{\Omega }^1`$. If the disk cools itself more rapidly then shock heating cannot match cooling and fragmentation results. This argument is apparently contradicted by Figure 9. The resolution lies in finding an appropriate definition of cooling time. The disk loses thermal energy on the effective cooling timescale

$$\tau _c^1\frac{\mathrm{\Lambda }}{U}.$$
(23)
Figure 10 shows the run of $`\tau _c`$ at the fragmentation, transition, and non-fragmentation points. Evidently $`\tau _c`$ at transition lies between $`\mathrm{\Omega }^1`$ and $`10\mathrm{\Omega }^1`$. Figure 11 shows the run of $`\tau _{co}`$ and $`\tau _c`$ on the transition line. Just below the opacity gap they differ by as much as four orders of magnitude.

Why do $`\tau _{co}`$ and $`\tau _c`$ differ by such a large factor? The answer is related to the existence of sharp variations in opacity with temperature. Consider a disk near the lower edge of the opacity gap. Once gravitational instability sets in, fluctuations in temperature move parts of the disk into the opacity gap. There, the opacity is reduced by orders of magnitude. Since the cooling rate for an optically thick disk is proportional to $`\kappa ^1`$, the cooling time drops by a similar factor. Relatively small variations in temperature can thus produce large variations in cooling rate.

As in Gammie (2001), the result $`\tau _c\mathrm{\Omega }1`$ also implies a constraint on $`\alpha `$. Energy conservation implies that

$$\frac{3}{2}\mathrm{\Omega }W_{xy}=\mathrm{\Lambda },$$
(24)
where $`W_{xy}`$ is the total shear stress (hydrodynamic plus gravitational). Equivalently, stress by rate-of-strain is equal to the dissipation rate. Using the definition of $`\tau _c`$, this implies

$$\alpha =\left(\gamma (\gamma 1)\frac{9}{4}\mathrm{\Omega }\tau _c\right)^1.$$
(25)
Hence $`\tau _c\mathrm{\Omega }1`$ implies $`\alpha 1`$. Figure 12 shows $`\alpha `$ vs $`\tau _c`$ for a large number of runs plotted against equation (25). For small values of $`\tau _c`$ the numerical values lie below the line. These models are not in equilibrium (i.e., not in a statistically-steady gravito-turbulent state), so the time average used in equation (24) is not well defined. For larger values of $`\tau _c`$ numerical results typically (there is noise in the measurement of both $`\alpha `$ and $`\tau _c`$ because the time average is taken over a finite time interval) lie slightly above the analytic result.

The bias toward points lying slightly above the line reflects the fact that $`\alpha `$ measures the rate of energy extraction from the shear while $`\tau _c`$ measures the rate at which that energy is transformed into thermal energy. If energy is lost, perhaps to numerical averaging at the grid scale, then more energy must be extracted from the shear flow to make up the difference. Overall, however, the agreement with the analytic result is good and demonstrates good energy conservation in the code.

The relationship between $`\tau _c`$ and $`\alpha `$ is interesting but not particularly useful because $`\tau _c`$ is no more readily calculated than $`\alpha `$; it depends on a complicated moment of the surface density and temperature. Only for constant cooling time have we been able to evaluate this moment analytically.

### 4.3 Isothermal Disks

We have assumed that external illumination of the disk is negligible. This approximation is valid when the effective temperature $`T_{irr}`$ of the external irradiation is small compared to the central temperature of the disk. In the opposite limit, illumination controls the energetics of the disk and it is isothermal (if it is illuminated directly so that shadowing effects, such as those considered by Jang-Condell & Sasselov (2003) are negligible).

It is therefore worth studying the outcome of gravitational instability in an isothermal disk. The isothermal disk model has a single parameter: the initial value of $`Q`$. We ran models with varying values of $`Q`$ and with $`\delta v^2/c_s^2=0.1`$. We find that models with $`Q1.4`$ fragment.

It is likely that the mass of the fragments, etc., depends on how an isothermal disk becomes unstable. Rapid fluctuation of the external radiation field is likely to produce a different outcome than dimming on a timescale long compared to the dynamical time.

## 5 Discussion

Using numerical experiments, we have identified those disks that are likely to fragment absent external heating. Disks with effective cooling times $`\tau _c\mathrm{\Omega }^1`$ are susceptible to fragmentation. This is what one might expect based on the simple argument of Shlosman et al. (1990): if the disk cools more quickly than the self-gravitating condensations can collide with one another, then those collisions (which occur on a timescale $`\mathrm{\Omega }^1`$) cannot reheat the disk and fragmentation is inevitable. But our results are at the same time surprising.

The effective cooling time depends on the nonlinear outcome of gravitational instability. It depends on the cooling function, which in turn depends sensitively on $`\mathrm{\Sigma }`$ and $`U`$. Since $`\mathrm{\Sigma }`$ and $`U`$ vary strongly over the disk once gravitational instability has set in, it is difficult to estimate $`\tau _c`$ directly. One might be tempted to estimate $`\tau _c(\mathrm{\Sigma },\mathrm{\Omega })\tau _{co}(\mathrm{\Sigma }_o,\mathrm{\Omega },Q=1)`$, but our experiments show that this estimate can be off by as much as four orders of magnitude. The effect is particularly pronounced near sharp features in the opacity. For example, consider a model initially located just below the opacity gap with $`\tau _{co}\mathrm{\Omega }1`$. Gravitational instability creates dense regions with higher temperatures, where dust is destroyed. The result is rather like having to shed one’s blanket on a cold winter morning: the disk loses its thermal energy suddenly. Pressure support is lost and gravitational collapse ensues.

The difference between $`\tau _c`$ and $`\tau _{co}(Q=1)`$ implies that a much larger region of the disk is susceptible to fragmentation than naive estimates based on the approximation $`\tau _c\tau _{co}`$ might suggest. For example, consider an equilibrium disk model with $`Q1`$ at small $`r`$. As $`r`$ increases, $`Q`$ declines. Eventually $`Q1`$ and gravitational instability sets in. There is then a range of radii where $`Q1`$, $`\tau _c\mathrm{\Omega }1`$ and recurrent gravitational instability can transport angular momentum and prevent collapse. Generally speaking, however, the cooling time decreases with increasing radius. Eventually $`\tau _c\mathrm{\Omega }1`$ and fragmentation cannot be avoided. By lowering our estimate of $`\tau _c`$, we narrow the range of radii over which recurrent gravitational instability can occur.

The general sense of our result is that it is extremely difficult to prevent a marginally-stable, $`Q1`$, optically-thick disk from fragmenting and forming planets (in circumstellar disks) or stars (in circumstellar and circumnuclear disks). This is particularly true for disks with $`T10^3\mathrm{K}`$, whose opacity is dominated by dust grains, i.e. disks whose temperature lies within a factor of several of the opacity gap.

Our numerical model uses a number of approximations. First, our treatment is razor-thin, i.e. all the matter is in a thin slice at $`z=0`$. The effect of finite thickness on linear stability has been understood since Goldreich & Lynden-Bell (1965): it is stabilizing because gravitational attraction of neighboring columns of disk is diluted by finite thickness. The size of the effect may be judged by the fact that $`Q=0.676`$ is required for marginal stability of a finite-thickness, isothermal disk.

The behavior of a finite-thickness disk in the nonlinear regime is more difficult to predict. Shocks will evidently deposit some of their energy away from the midplane, where it can be radiated away more quickly (because the energy is deposited at smaller optical depth - see Pickett et al. 2000). Radiative diffusion parallel to the disk plane (not included here) may enhance cooling of dense, hot regions. Both these effects are destabilizing. Ultimately, however, a numerical study is required. This is numerically expensive: one must resolve the disk vertically, on the scale height $`H`$, and horizontally, at the critical wavelength $`2\pi QH`$.

Second, we have ignored magnetic fields. While there may be astrophysical situations where cool disks have such low ionization that they are unmagnetized, most disks are likely to contain dynamically important magnetic fields that give rise to a dimensionless shear stress $`\alpha 0.01`$ (e.g. Hawley et al. 1995). These fields are likely to remove spin angular momentum from partially collapsed objects, destabilizing them. Numerical experiments including both gravitational fields and magnetohydrodynamics are necessarily three dimensional (the instability of Balbus & Hawley (1991) requires $`_z0`$), and are thus numerically expensive.

Third, we have fixed $`\gamma `$ and $`\mu `$ for the duration of each simulation. This eliminates the soft spots in the equation of state associated with ionization of atomic hydrogen and dissociation of molecular hydrogen. In these locations the three dimensional $`\gamma `$ dips below $`4/3`$, which is destabilizing.

Fourth, we have treated the physics of grain destruction and formation very simply. In using the Bell & Lin (1994) opacities we implicitly assume that grains reform in cooling gas on much less than a dynamical time. It is likely that grain re-formation will take some time (e.g. Hessman 1991) and this will further reduce the disk opacity and enhance fragmentation.

Fifth, we have neglected the effects of illumination. In the limit of strong external illumination, i.e. when the effective temperature of the irradiation $`T_{irr}`$ is large compared to the disk central temperature $`T_c`$, the disk is isothermal (here $`T_c`$ is the temperature of a dense condensation). We have carried out isothermal experiments and shown that, for initial velocity perturbations with $`\delta v^2/c_s^2=0.1`$, disks with $`Q1.4`$ fragment. Weaker illumination produces a more complicated situation that we have not explored here. Illumination-dominated disks that become unstable presumably do so because the external illumination declines, and the rate at which the external illumination changes may govern the nonlinear outcome.

We conclude that disks with $`\tau _c\mathrm{\Omega }1`$ do not exist. Cooling in this case is so effective that fragmentation into condensed objects– stars, planets, or smaller accretion disks– is inevitable.

As an example application of this result, consider the model for the nucleus of NGC 1068 recently proposed by Lodato & Bertin (2003). Their model is an extended marginally-stable self-gravitating disk of the type investigated here and originally proposed by Goldreich & Lynden-Bell (1965) for galactic disks and Paczynski (1978) for accretion disks, although their disk is sufficiently massive that it modifies the rotation curve as well. Based on their Figure 3, at a typical radius of $`0.5\mathrm{pc}`$, $`\mathrm{\Sigma }_o10^4`$ and $`\mathrm{\Omega }10^9`$. According to our Figure 7 this disk is about 2 orders of magnitude too dense to avoid fragmentation. While it may be possible to avoid this conclusion by invoking strong external heating, the energy requirements are severe, as outlined in Goodman (2003). The disk proposed by Lodato & Bertin (2003) would therefore fragment into stars on a short timescale.

This work was supported by NSF grant AST 00-03091, PHY 02-05155, and NASA grant NAG 5-9180.


format: text

text:
# Bulk Viscosity in Hybrid Stars

## I Introduction

The discovery by Andersson, Friedman and Morsink of r-mode instabilities in neutron stars put rather severe limits on the highest rotation frequency of pulsars Andersson (1998); Friedman and Morsink (1998). These constraints can be incompatible with the existence of millisecond pulsars, if the instability is not suppressed by a sufficiently large viscosity. Actually, for a star composed only of neutrons and protons and for temperatures larger than roughly 1010 K, the bulk viscosity due to the modified Urca process is large enough to damp the instability Sawyer (1989a). On the other hand, when the star cools down to lower temperatures, the instability is not suppressed and the star is forced to loose angular momentum via emission of gravitational waves Lindblom et al. (1999). More recently, it has been noticed that the existence of viscous boundary layers at the interface between the fluid core and the crust can stabilize the star 111Also the presence of a magnetic field inside the compact star can damp the r-modes on a time-scale of order hours or days Rezzolla et al. (2000). if the rotation period is longer than $``$ 1.5 ms Bildsten and Ushomirsky (1999); Andersson et al. (2000). At very low temperatures, below 108 K, shear viscosity becomes large and it allows older stars to increase their angular velocity by mass accretion.

Actually, compact stars can be constituted by a larger variety of particles than just neutrons and protons. One possibility is that hyperons form at the center of the star. It has been shown that, due to non-leptonic weak reactions, bulk viscosity can be rather large for an hyperonic star Jones (2001); Lindblom and Owen (2002), which therefore can emit gravitational waves only if its temperature is $`10^{10}`$K and its frequency is larger than 10–30% of its Keplerian frequency.

The formation of quark matter inside a compact star has been discussed extensively in the literature. Bulk viscosity of non-interacting strange quark matter is very large Wang and Lu (1984); Sawyer (1989b); Madsen (1992). On the other hand, recent studies taking into account quark-quark interaction revealed the possible existence of color superconducting phases in which quarks form Cooper pairs with gaps as large as 100 MeV Rajagopal and Wilczek (2000). In that case, bulk viscosity is strongly suppressed by the large energy gaps, r-mode instabilities are not damped and pure color-superconducting quark stars seem therefore to be ruled out by the pulsar data Madsen (2000).

In the present paper we are interested in studying the viscosity and stability, respect to r-modes, of an Hybrid Star (HyS), for which no quantitative analysis has been performed so far. If the star is made of hyperonic matter and of non-interacting quarks, it is rather obvious that the viscosity of the HyS should be large, due to the large viscosity of its constituents. The real question concerns the case in which color-superconducting quark matter is considered, since its viscosity is negligible. There are in principle (at least) two sources of viscosity for HySs. One originates at the interface between the crust and the fluid interior Madsen (2000) and it is similar to the one discussed above in the case of purely hadronic stars Bildsten and Ushomirsky (1999); Andersson et al. (2000). In our work we will not discuss that possibility and we will instead concentrate on a quantitative evaluation of the bulk viscosity in a mixed phase (MP) composed either of hyperons and color-superconducting quark matter or of nucleons and non-interacting quarks. We have also investigated the effect of the existence of finite-size structures in the MP. These geometrical structures can exist if the surface tension at the interface between hadronic and quark matter is not vanishing. We will show that the bulk viscosity of the MP is in general rather large, particularly so when the surface tension is negligible. Therefore HySs are possible candidates for young millisecond pulsars.

## II Equation of State

We construct the Equation of State (EOS) of matter at high density modeling the Hadronic phase by a relativistic non-linear Walecka type model Liu et al. (2002) with the inclusion of Hyperons Glendenning (1997). Concerning quark matter, we use an MIT-bag like EOS in which quarks can pair to form a Color Flavor Locked (CFL) phase Alford et al. (1999, 2001); Alford and Reddy (2003). CFL is considered to be the energetically favored type of pairing pattern, at large densities, in the case of $`\beta `$-stable, electric-charge neutral quark matter. At lower densities, a two-flavor Color Superconducting (2SC) phase can form, with a smaller energy gap Neumann et al. (2003).

Two scenarios are possible, depending on the value of the critical density separating hadronic matter from MP. In the first case the critical density is large, therefore hyperons start being produced at a density smaller than the critical one. We will call this phase “hyperon-quark” MP Banik and Bandyopadhyay (2003). This is the first case discussed below and it corresponds to the upper panel of Fig. 1. In the second scenario, the MP starts at a density lower than the hyperonic threshold and the hyperon density is much smaller than in the previous case. We will explore in particular the situation in which hyperons are completely eaten-up by quarks and they do not appear in the “nucleon-quark” MP (see lower panel of Fig. 1). In this second scenario we will assume that in the MP the density is so low that the CFL pairing cannot form and, for simplicity, we will consider unpaired quarks in the MP, since anyway the bulk viscosity of the 2SC quark phase is similar to the viscosity of unpaired quark matter Madsen (2000). In both scenarios, we describe the composition of the MP assuming a first order transition and imposing Gibbs conditions 222The structure of the MP obtained imposing Gibbs condition is clearly different from a two-fluid model like the one adopted e.g. in Ref. Andersson et al. (2002), in which the two phases are essentially independent. Beta stability and charge neutrality are satisfied in every density region.

## III Bulk Viscosity

Bulk viscosity is the dissipative process in which a perturbation of the pressure of a fluid element is converted to heat. A small variation of the pressure of the system can be treated as a perturbation on the densities of the different species of particles bringing the system out of $`\beta `$-equilibrium. The reactions between the different particles drive the system back to an equilibrium configuration, with a delay which depends on the characteristic time scale of the interactions. Following the formalism of Ref.Lindblom and Owen (2002), we classify all the reactions either as fast or slow. Slow processes produce bulk viscosity, because their time scale is comparable with the period of the perturbation. Fast processes, instead, put constraints on the variation of the densities of the different particles. In the following we generalize the formalism of Ref. Lindblom and Owen (2002) in order to compute the bulk viscosity of MP. The equations needed in the computation are rather different in the two scenarios discussed above and we will have to deal with all the different cases separately. Moreover, if the effect of a non-vanishing surface tension is taken into account, the formalism needed to compute the bulk viscosity has to be further modified. It is well known that a finite value for the surface tension allows the formation of finite size structures whose geometrical shape is determined by the interplay of the various contribution to the energy, including the Coulomb term Heiselberg et al. (1993); Voskresensky et al. (2003). A precise estimate of surface tension $`\sigma `$ is unfortunately still lacking, and values ranging from a few MeV/fm2 to a few tens MeV/fm2 have been discussed in the literature. We can roughly divide this range of values in three windows. Values larger than $``$ 30 MeV/fm2 would not allow the formation of the mixed phase since it would not be energetically favoredAlford et al. (2001); Voskresensky et al. (2003). Our analysis is therefore restricted to smaller values. In the following we will discuss both the case in which the surface tension is so small that it can be neglected and the case in which it is non-vanishing.

### III.1 Negligible surface tension

Let us first discuss the case in which the surface tension is so small that the perturbations of the star, associated e.g. to r-modes, can break the finite-size structures present in the MP. We have estimated that this case corresponds to values of the surface tension $`\sigma `$1 MeV/fm2, since to r-modes excitations is associated an energy per baryon of a few MeV.

#### III.1.1 First scenario

We start our analysis from the first scenario discussed in Sec. II, namely from the case in which hyperons are present in the MP. We will discuss in particular a parameter set in which only $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }^{}`$ particles are produced in the MP. The only slow processes generating viscosity are non-leptonic reactions between hadrons 333We assume, as in Ref.Lindblom and Owen (2002), that all leptonic reaction rates are much smaller than those associated with non-leptonic reactions, so we do not include them in the calculation of bulk viscosity of the MP. The viscosity due to semi-leptonic reactions (modified Urca processes) becomes relevant at very high temperatures and we will include its effect when computing the stability of the star., since the large gap prevents weak reactions between quarks. As in Ref.Lindblom and Owen (2002), we consider the reactions

$`n+n`$ $`\stackrel{H_W}{}`$ $`p+\mathrm{\Sigma }^{},`$ (1)
$`n+p`$ $`\stackrel{H_W}{}`$ $`p+\mathrm{\Lambda }.`$ (2)
Following Lindblom and Owen (2002), we have not taken into account the reaction $`n+n\stackrel{H_W}{}n+\mathrm{\Lambda }`$, since the corresponding reaction rate cannot be easily estimated. In Ref.Jones (2001) it has been argued that this rate can be one order of magnitude larger that the one associated with reaction (2). In the first scenario, our results have therefore to be considered as upper limits for the viscosity, as it will be clarified in the following. Notice anyway that, concerning the damping of r-modes instabilities, to be discussed in the following, the most important region of the star corresponds to low–moderate densities, while the $`\mathrm{\Lambda }`$ is produced at larger densities as shown in Figs. (1) and (5).

Concerning fast processes, they come from the following reactions mediated by the strong interaction

$`n+\mathrm{\Lambda }`$ $`\stackrel{H_S}{}`$ $`p+\mathrm{\Sigma }^{},`$ (3)
$`\mathrm{\Lambda }+\mathrm{\Lambda }`$ $`\stackrel{H_S}{}`$ $`2(uds)_{CFL}.`$ (4)
The last process describes the only possible “melting” of hadrons into CFL phase, since the pairing forces the number of up, down and strange quarks to be equal. This process is fast because of the vanishing value of $`\sigma `$. As we will see later the “melting” process is instead forbidden if the value of $`\sigma `$ is not negligible. Concerning mechanical equilibrium, elastic scattering due to strong interactions, as well as melting processes like the one described in Eq.(4), are responsible for (rapid) momentum transfer between the two phases. The mechanical equilibrium between the two phases is reached in a time scale $`t_s`$ determined by strong interaction, $`t_s10^{23}`$ s. The full system is, instead, out of mechanical equilibrium on a time scale of the order of the period of the perturbation $`t_p10^3`$ s. Therefore, during a fluctuation the two components of the fluid remain in mutual mechanical equilibrium. The variations of the pressure in the two phases have therefore to satisfy the constraint:

$$\delta P_h=\delta P_{CFL}.$$
(5)
The variations $`\delta \rho _i`$ of the densities of the various particles and the variation $`\delta \chi `$ of the quark fraction are constrained by the following linearized equations

$`0`$ $`=`$ $`(1\chi )(\delta \rho _n+\delta \rho _p+\delta \rho _\mathrm{\Lambda }+\delta \rho _\mathrm{\Sigma })`$ (6)
$`+`$ $`\chi (\delta \rho _q)+\delta \chi (\rho _q\rho _n\rho _p\rho _\mathrm{\Lambda }\rho _\mathrm{\Sigma }),`$
$`0`$ $`=`$ $`(1\chi )(\delta \rho _p\delta \rho _\mathrm{\Sigma })\delta \chi (\rho _p\rho _\mathrm{\Sigma }),`$ (7)
$`0`$ $`=`$ $`{\displaystyle \underset{\{H\}}{}}p_H\delta \rho _Hp_q\delta \rho _q,`$ (8)
$`0`$ $`=`$ $`\beta _n\delta \rho _n+\beta _p\delta \rho _p+\beta _\mathrm{\Lambda }\delta \rho _\mathrm{\Lambda }+\beta _\mathrm{\Sigma }\delta \rho _\mathrm{\Sigma },`$ (9)
$`0`$ $`=`$ $`\alpha _{\mathrm{\Lambda }n}\delta \rho _n+\alpha _{\mathrm{\Lambda }p}\delta \rho _p+\alpha _{\mathrm{\Lambda }\mathrm{\Lambda }}\delta \rho _\mathrm{\Lambda }`$ (10)
$`+`$ $`\alpha _{\mathrm{\Lambda }\mathrm{\Sigma }}\delta \rho _\mathrm{\Sigma }3\alpha _{qq}\delta \rho _q,`$
where

$`\alpha _{ij}`$ $`=`$ $`\left(\mu _i/\rho _j\right)_{\rho _k,kj},`$ (11)
$`\beta _i`$ $`=`$ $`\alpha _{ni}+\alpha _{\mathrm{\Lambda }i}\alpha _{pi}\alpha _{\mathrm{\Sigma }i},`$ (12)
$`p_i`$ $`=`$ $`P/\rho _i.`$ (13)
Eqs. (6),(7) impose baryon number conservation and electric charge neutrality. Eq. (8) (where the sum runs on all hadrons) imposes the mechanical equilibrium defined by Eq. (5). Finally, Eqs. (9),(10) describe the equilibrium respect to the two strong processes of Eqs. (3),(4). Notice that $`\delta \chi `$ does not appear in Eqs. (8)–(10), because neither the pressure nor the chemical potentials explicitely depends on the quark volume fraction $`\chi `$. Solving the system allows to express all the $`\delta \rho _i`$ and $`\delta \chi `$ as function of $`\delta \rho _n`$ 444If the $`\mathrm{\Sigma }`$ density vanishes, the corresponding density fluctuation is identically zero and the constraint given by Eq. (3) has not to be imposed. If the $`\mathrm{\Lambda }`$ density vanishes, then $`\delta \rho _\mathrm{\Lambda }=0`$, Eq. (3) has not to be imposed, the melting process described by Eq. (4) does not exist and $`\delta \chi =0`$..

The relaxation time $`\tau `$ associated to the weak processes reads

$$\frac{1}{\tau }=\left(\frac{\mathrm{\Gamma }_\mathrm{\Lambda }}{\delta \mu }+2\frac{\mathrm{\Gamma }_\mathrm{\Sigma }}{\delta \mu }\right)\frac{\delta \mu }{\delta \rho _n}.$$
(14)
Here $`\mathrm{\Gamma }_\mathrm{\Lambda }`$ and $`\mathrm{\Gamma }_\mathrm{\Sigma }`$ are the rates of the weak interactions 555Notice that the Coulomb interaction does not play any direct role when computing the viscosity in the present scheme, since it cannot modify the value of the chemical unbalance $`\delta \mu `$.. They have been calculated in Eqs. (4.21), (4.28), (4.29) of Ref.Lindblom and Owen (2002), where it has been shown that the main dependence of the rate on the temperature is a quadratic one. Moreover, if hyperon superfluid gaps as the ones displayed in Fig. 2 do develop, they exponentially suppress the rates as it will be discussed later.

The chemical potential unbalance $`\delta \mu `$ is given by

$$\delta \mu \delta \mu _n\delta \mu _\mathrm{\Lambda }=2\delta \mu _n\delta \mu _p\delta \mu _\mathrm{\Sigma },$$
(15)
where the variations of the chemical potentials can be expressed in terms of $`\delta \rho _i`$ as $`\delta \mu _i=_j\alpha _{ij}\delta \rho _j.`$ A unique $`\delta \mu `$ appears in Eq. (15), since the variations $`\delta \mu _\mathrm{\Lambda }`$ and $`\delta \mu _\mathrm{\Sigma }`$ are constrained by the fast process given in Eq. (3). The real part of bulk viscosity, which is the relevant quantity for damping r-modes instabilities, finally reads

$$\zeta =\frac{P(\gamma _{\mathrm{}}\gamma _0)\tau }{1+(\omega \tau )^2},$$
(16)
where $`\gamma _{\mathrm{}}`$ and $`\gamma _0`$ are the “infinite” frequency adiabatic index and the “zero” frequency adiabatic index and $`\omega `$ is the angular velocity of the perturbation.

It is interesting to remark that there are two asymptotic behaviors of the viscosity as a function of the frequency $`\omega `$. In the high frequency case ($`\omega \tau 1`$), the viscosity scales as $`1/\tau `$, while in the low frequency limit the viscosity is proportional to $`\tau `$. In the parameter and temperature ranges explored in our paper it results that we always remain in the low frequency limit, as it can be seen from Fig. 3 noticing that the viscosity decreases with the temperature. In this regime the addition of another weak decay (e.g. $`n+n\stackrel{H_W}{}n+\mathrm{\Lambda }`$) decreases the relaxation time and as consequence the viscosity too. In this sense our results for the viscosity in the first scenario must be considered as upper limits.

We can observe from Fig. 3 that, in the “hyperon-quark” MP scenario discussed so far, MP bulk viscosity is comparable to bulk viscosity of purely hyperonic matter even though the density of hyperons is very small in the MP (see upper panel of Fig. 1) 666Our results in both scenarios are consistent with the general outcome of Ref.Haensel et al. (2001) stating that in the high-frequency limit the bulk viscosity is just the sum of partial bulk viscosities in different slow reaction channels.. Notice that in the small window of baryonic densities in which both $`\mathrm{\Sigma }`$ and $`\mathrm{\Lambda }`$ hyperons are present, the value of bulk viscosity of the MP and of the pure hadronic phase are essentially identical. On the other hand, in the baryonic density windows in which only one species of hyperons are present in the MP, the bulk viscosity of MP is larger than the bulk viscosity of the pure hadronic phase because only one decay channel is open ($`\mathrm{\Sigma }`$ decay of Eq. (1) or $`\mathrm{\Lambda }`$ decay of Eq. (2)) and therefore the relaxation time is larger. It is interesting to compute the bulk viscosity in this scenario also including the effect of hyperon superfluidity, as done in Ref. Lindblom and Owen (2002). If the energy gaps $`\mathrm{\Delta }_H`$ associated with the hyperons are non vanishing, the decay rates are suppressed by the factor $`e^{\mathrm{\Delta }_H/T}`$ where the gaps have a typical shape shown in Fig. 2. For low temperatures the viscosity displayes the characteristic features shown in the lower panel of Fig. 3, while these features disappear at larger temperatures or for vanishing $`\mathrm{\Delta }_H`$. Concerning the appreciable difference between the viscosity of MP and of pure hadronic matter for $`T=10^{9.5}`$K, this is due to the Fermi momentum dependence of $`\mathrm{\Delta }_H`$, which implies that the gaps are suppressed for large hyperon densities. In the MP, the density of hyperons is lower than in pure hadronic matter and therefore the effect of the gaps shows out more dramatically.

Let us briefly discuss another possible source of bulk viscosity related to the formation of a boson condensate. Since all quarks in the CFL phase are gapped, the low energy excitations are the Nambu-Goldstone bosons associated with the spontaneous symmetry breaking of the global symmetries Alford et al. (1999). In particular, a condensate of $`\pi ^{}`$ or $`K^{}`$ can appear, in the MP, as proposed in Ref. Alford and Reddy (2003). The role played by these boson condensates in the cooling of a CS has been studied in Ref. Jaikumar et al. (2002) and the condensates can play a role also in the calculation of the bulk viscosity of the MP. The weak processes involving these bosons are:

$`\pi ^{}`$ $`\stackrel{H_W}{}`$ $`e^{}+\overline{\nu }_e,`$ (17)
$`K^{}`$ $`\stackrel{H_W}{}`$ $`e^{}+\overline{\nu }_e.`$ (18)
In the absence of these bosons, the variation of the density of electrons $`\delta \rho _e`$ can be neglected at low temperature ($`T<10^{10}`$K) because, as already remarked, all leptonic reaction rates are much smaller than those associated with non-leptonic reactions. If, on the other hand, the decay channels of Eqs. (17),(18) are open, then the corresponding decay rates are of the same order of the rates of the non-leptonic processes of Eqs. (1),(2) and, therefore, a new contribution to the viscosity appears. While a detailed calculation is clearly needed, the main result of our work, namely the existence of a large viscosity for $`T10^{10}`$ K, would even be strengthened.

#### III.1.2 Second scenario

We now discuss the second scenario in which a “nucleon-quark” MP forms, hyperons are absent and CFL quark pairing cannot take place in the MP 777Notice that if CFL gaps can form at low density then the bulk viscosity of MP vanishes.. We will assume that in the pure quark matter phase CFL gaps can form so that this phase will not contribute to the viscosity. The only source of viscosity, neglecting as before semi-leptonic reactions, is therefore

$$d+u\stackrel{H_W}{}u+s.$$
(19)
Concerning fast processes, they are given by the “melting” of nucleons into unpaired quarks. The chemical equilibrium respect to these reactions must therefore be satisfied during the perturbation. The linearized equations governing the density fluctuations, analogous to Eqs. (6)-(10) of the first scenario, read

$`0`$ $`=`$ $`(1\chi )(\delta \rho _n+\delta \rho _p)+\chi (\delta \rho _u+\delta \rho _d+\delta \rho _s)/3`$ (20)
$`+`$ $`\delta \chi \left((\rho _u+\rho _d+\rho _s)/3\rho _n\rho _p\right),`$
$`0`$ $`=`$ $`(1\chi )\delta \rho _p+\chi (2\delta \rho _u\delta \rho _d\delta \rho _s)/3`$ (21)
$`+`$ $`\delta \chi \left((2\rho _u\rho _d\rho _s)/3\rho _p\right),`$
$`0`$ $`=`$ $`{\displaystyle \underset{\{H\}}{}}p_H\delta \rho _H{\displaystyle \underset{\{Q\}}{}}p_Q\delta \rho _Q,`$ (22)
$`0`$ $`=`$ $`\alpha _{pn}\delta \rho _n+\alpha _{pp}\delta \rho _p2\alpha _{uu}\delta \rho _u\alpha _{dd}\delta \rho _d,`$ (23)
$`0`$ $`=`$ $`\alpha _{nn}\delta \rho _n+\alpha _{np}\delta \rho _p\alpha _{uu}\delta \rho _u2\alpha _{dd}\delta \rho _d.`$ (24)
Eqs. (20)–(22) impose baryon number conservation, electric charge neutrality and mechanical equilibrium as in the first case. Eqs. (23),(24) describe the chemical equilibrium respect to the two fast processes of “melting”. As before we can calculate the chemical imbalance associated to the weak reaction (19) and the corresponding relaxation time. The rate of the process (19) is taken from Eqs. (5)–(7) of Ref.Sawyer (1989b)(see also footnote 888The reaction rate has been multiplied by a factor three, in agreement with Ref.Madsen (1992).). The resulting bulk viscosity is shown in Fig. 3 (lower panel). The bulk viscosity of “nucleon-quark” MP is of the same order of magnitude of bulk viscosity of strange matter Madsen (1992).

### III.2 Effects of the surface tension

Let us shortly discuss the effect of a non-vanishing surface tension $`\sigma `$ at the interface between hadronic and quark matter. We will assume that $`\sigma `$ 30 MeV/fm2, so that MP can form because it is energetically favored. Finite-size structures (drops, rods, slabs) form at different densities, to minimize the energy of the MP. In this case the perturbation of the pressure due to r-modes is too weak to induce the “melting” process which is now a very slow process (in comparison to the period of the perturbation) and it plays no role in the calculation of the bulk viscosity. The response of the finite size structures inside MP to a perturbation of the density corresponds to three possible processes: a) the formation of a drop of “new” phase; b) the merging of two structures into a single larger structure; c) the absorption of “old” phase into a structure of “new” phase. Obviously, also the reverse processes are possible. It is easy to see that large values of $`\sigma `$ suppress all these processes. In particular, in case a) the radius of a critical drop of new phase increases with $`\sigma `$, making it more difficult to produce a new drop; case b) shares some similarity with the fission problem in nuclear physics, where the process of the separation of a heavy nucleus into two lighter nuclei is suppressed for larger values of $`\sigma `$ since, during the fission (or merging) process, configurations having a large surface are produced. Finally, c) can be viewed as a special case of b), in which the absorbed hadron can be assimilated to a small drop of quark matter. Concerning the Coulomb interaction, it mainly plays a role in determining the size of the structures while it is not so important for the response of the structures to the perturbation, at least for relatively large values of the surface tension $`\sigma `$ 10 MeV/fm2. In that case, in fact, the screening due to electrons almost completely cancels the effect of the Coulomb interaction in the nucleation rate Iida and Sato (1998). In the following for simplicity we will assume that 10 MeV/fm2 $`\sigma `$ 30 MeV/fm2. In conclusion we reasonably assume that large values of $`\sigma `$ suppress all these processes, therefore the equations describing the melting process as a fast reaction, i.e. Eq. (10) (first scenario) and Eqs. (23), (24) (second scenario) are not imposed. We can compute the viscosity by requiring that the baryon number and the electric charge of the two phases are separately conserved. This implies that Eq. (6) in the first scenario and Eqs. (20), (21) in the second scenario, each of them separate into two distinct equations. In the first scenario, Eq. (7) is not modified due to the charge neutrality of the CFL quark phase. Finally, the equation of mechanical equilibrium (Eq. (8) in the first scenario and Eq. (22) in the second scenario) is still valid and it represents the only fast process connecting the two phases. Notice that since melting processes are suppressed, the only reaction which allows the system to rapidly re-equilibrate is elastic scattering between the two phases. In the following we have assumed that a residual interaction between the two phases always exists, similarly to the existence of “entrainment” in superfluid neutron matter as discussed e.g. in Ref.Andersson et al. (2002). The systems of equations read therefore:

$`0`$ $`=`$ $`(1\chi )(\delta \rho _n+\delta \rho _p+\delta \rho _\mathrm{\Lambda }+\delta \rho _\mathrm{\Sigma })`$ (25)
$``$ $`\delta \chi (\rho _n+\rho _p+\rho _\mathrm{\Lambda }+\rho _\mathrm{\Sigma }),`$
$`0`$ $`=`$ $`\chi \delta \rho _q+\delta \chi \rho _q,`$ (26)
$`0`$ $`=`$ $`(1\chi )(\delta \rho _p\delta \rho _\mathrm{\Sigma })\delta \chi (\rho _p\rho _\mathrm{\Sigma }),`$ (27)
$`0`$ $`=`$ $`{\displaystyle \underset{\{H\}}{}}p_H\delta \rho _Hp_q\delta \rho _q,`$ (28)
$`0`$ $`=`$ $`\beta _n\delta \rho _n+\beta _p\delta \rho _p+\beta _\mathrm{\Lambda }\delta \rho _\mathrm{\Lambda }+\beta _\mathrm{\Sigma }\delta \rho _\mathrm{\Sigma }`$ (29)
for the first scenario and

$`0`$ $`=`$ $`(1\chi )(\delta \rho _n+\delta \rho _p)\delta \chi (\rho _n+\rho _p),`$ (30)
$`0`$ $`=`$ $`\chi (\delta \rho _u+\delta \rho _d+\delta \rho _s)+\delta \chi (\rho _u+\rho _d+\rho _s),`$ (31)
$`0`$ $`=`$ $`(1\chi )\delta \rho _p\delta \chi \rho _p,`$ (32)
$`0`$ $`=`$ $`\chi (2\delta \rho _u\delta \rho _d\delta \rho _s)+\delta \chi (2\rho _u\rho _d\rho _s)`$ (33)
$`0`$ $`=`$ $`{\displaystyle \underset{\{H\}}{}}p_H\delta \rho _H{\displaystyle \underset{\{Q\}}{}}p_Q\delta \rho _Q`$ (34)
for the second scenario.

In Fig. 3 and 4 we show the effect of a non-vanishing surface tension on the viscosity (dotted lines). On rather general grounds, one can expect that the effect of a non-vanishing surface tension is to reduce the viscosity. Indeed, the surface tension suppresses the fast processes of “melting” which are responsible for the reduction of the chemical unbalance. In particular, the only equation connecting the two phases is now the equation corresponding to the mechanical equilibrium. Moreover, in the first scenario the number of constraints on $`\delta \rho _H`$ reduces from four to three and, in the second scenario, the constraints on $`\delta \rho _Q`$ from four to two. It is also interesting to notice that in both scenarios, near the beginning of the MP ($`\chi 0`$), baryon number conservation requires $`\delta \chi 0`$ if the surface tension is non-vanishing (see Eqs. (26) and (31)). On the other hand, in the first scenario $`\rho _\mathrm{\Lambda }=0`$ at the beginning of the MP below the $`\mathrm{\Lambda }`$ production threshold, and therefore, the constraint $`\delta \chi =0`$ is also satisfied in the absence of surface tension, as explained in footnote . Therefore the dotted lines coincide with the solid lines in both panels of Fig. 3 for $`\rho _B0.6`$ fm-3.

In the second scenario, the effect of the surface tension is more dramatic. In Fig. 6 we display the chemical unbalances corresponding to a vanishing and to a finite value of the surface tension, respectively. As already remarked $`\delta \mu /\delta \rho _n`$ is larger in presence of a surface tension. Concerning the singular behavior of the chemical unbalance near the first critical density, it stems from neglecting $`\delta \rho _e`$ in the equation of the electric charge conservation. While in general this is a safe approximation due to the slowness of the modified Urca process (see Ref. Lindblom and Owen (2002)), in this particular case this approximation implies the vanishing of the viscosity at threshold. Actually, the existence of a finite value for $`\delta \rho _e`$ implies that the viscosity is small but finite at threshold. It is not possible to apply directly the formalism of Ref. Lindblom and Owen (2002) if two independent  perturbations ($`\delta \rho _n`$ and $`\delta \rho _e`$) exist in the system. We are therefore forced to discuss separately the viscosities stemming from the two independent unbalances. Concerning the viscosity associated with the modified Urca process, we remind that it scales as $`T^6`$ and therefore it is essentially negligible below $`T10^{10}`$ K while it will be included in the calculation of the stability of the star presented in the next section. In conclusion, the result corresponding to the dotted line of Fig. 4 would not be significantly modified by taking into account a finite value of $`\delta \rho _e`$.

## IV Stability of Hybrid Stars

We can now address the problem of the stability of a rotating compact star. To compute the critical angular velocity we use the standard formalism of Refs.Lindblom et al. (1998, 1999); Lindblom and Owen (2002). We need first to integrate the viscosity on the structure of the star, which is obtained by solving Tolman-Oppenheimer-Volkov equation. In Fig. 5 we show the structure of a $`M=1.46M_{}`$ star for the two scenarios discussed above. For simplicity we computed the structure of the star assuming a negligible value for $`\sigma `$. The main effect of the presence of finite size structures is to reduce the volume occupied by the MP. For $`\sigma `$ 30 MeV/fm2, which is the limit of validity of our approach, the shrinking of the MP is rather modest Heiselberg et al. (1993). The critical angular velocity $`\mathrm{\Omega }_{crit}`$ is the one for which the imaginary part of the r-mode frequency vanishes and it is obtained by solving the equation

$$1/\tau _{GR}+1/\tau _B+1/\tau _{B(Urca)}=0.$$
(35)
Here $`\tau _{GR}`$ is the time scale for gravitational waves emission while $`\tau _B`$ and $`\tau _{B(Urca)}`$ are the time scales of the bulk viscosity produced by hadronic processes and by the modified Urca process of the nucleonsSawyer (1989a), respectively. Results for the critical angular velocity are shown in Fig. 7. In the upper panel we compare the stability of a purely hyperonic star with the stability of a HyS containing hyperon-quark MP. As it can be seen, due to the large viscosity of the MP the HyS is as stable as the hyperonic star. Let us stress again that our result indicate that the viscosity in the MP is almost independent on the hyperonic content and therefore hybrid hyperon-quark stars can be stable as long as a tiny fraction of hyperon is present. In the lower panel we compare the stability of a star made entirely of non-interacting quarks with a HyS containing a nucleon-quark MP. We have assumed that the viscosity of the pure quark matter phase in the HyS vanishes, to simulate a CFL core. As it can be seen, the small MP region, located near the edge of the star, is sufficient to damp the r-modes. For simplicity, we have assumed quarks to be unpaired in the MP, but similar results should be obtainable if a 2SC phase is present. A feature of r-modes is that they are active mostly in the outer regions of the star, and therefore the value of the bulk viscosity at a not too large densities is crucial for the stability of the star. In particular, in the first scenario the stability of the star is due to the presence of the $`\mathrm{\Sigma }`$ hyperons in the pure hadronic phase and in the MP. $`\mathrm{\Lambda }`$ hyperons, which are produced at larger densities, play a lesser role. In the second scenario, a small window of MP, present in the outer region, is sufficient to stabilize the star. Let us also remark that the star is stable at large temperatures, due to the modified Urca processes active in the crust. These processes does not exist in the case of pure quark stars, which are therefore unstable at large temperatures. Finally, when a finite value for the hadron-quark surface tension is taken into account the instability window is larger, but the main conclusion concerning the stability of a young hybrid star remains valid.

It is a pleasure to thank J.C. Miller and L. Rezzolla for very useful discussions.


text:
# 1 Introduction

## 1 Introduction

During the last few years, logarithmic conformal field theory (LCFT) has been established as a well-defined variety of conformal field theories in two dimensions. The concept was considered in its own right first by Gurarie , Since then, a large amount of work has appeared, see the reviews  and references therein. The defining feature of a LCFT is the occurence of indecomposable representations which, in turn, may lead to logarithmically diverging correlation functions. Thus, in the standard example of a LCFT a primary field $`\varphi \left(z\right)`$ of conformal weight $`h`$ has a so-called logarithmic partner field $`\psi `$ with the characteristic properties

$$\varphi \left(z\right)\varphi \left(0\right)=0,\varphi \left(z\right)\psi \left(0\right)=Az^{2h},\psi \left(z\right)\psi \left(0\right)=z^{2h}\left(B2A\mathrm{log}\left(z\right)\right).$$
(1)
To this corresponds the fact that the highest weigth state $`|h`$ associated to the primary field $`\varphi `$ is the ground state of an irreducible representation which, however, is part of a larger, indecomposable, representation created from $`|\stackrel{~}{h}`$, the state associated to $`\psi `$. The conformal weight is the eigenvalue under the action of $`L_0`$, the zero mode of the Virasoro algebra, which in such LCFTs cannot be diagonalized. Instead, we have

$$L_0|h=h|h,L_0|\stackrel{~}{h}=h|\stackrel{~}{h}+|h.$$
(2)
Thus, the two states $`|h`$ and $`|\stackrel{~}{h}`$ span a Jordan cell of rank two with respect to $`L_0`$. As can be guessed from eq. (1), there must exist a zero mode which is responsible for the vanishing of the 2-pt function of the primary field. Another characteristic fact in LCFT is the existence of at least one field, which is a perfect primary field, but whose operator product expansion (OPE) with itself produces a logarithmic field. Such fields $`\mu `$ are called pre-logarithmic fields . This is important, since in many cases, the pre-logarithmic fields arise naturally forcing us then to include the logarithmic fiels as well into the operator algebra. Note that this implies that the fusion product of two irreducible representations is not necessarily completely reducible into irreducible representations. In fact, we know today quite a few LCFTs, where precisely this is the case, such as ghost systems , WZW models at level zero or at fractional level such as $`\widehat{SU\left(2\right)}_{4/3}`$ , WZW models of supergroups such as $`GL(1,1)`$  or certain supersymmetric $`c=0`$ theories such as $`OSP\left(2n|2n\right)`$ or $`CP\left(n|n\right)`$ . Finally, many LCFTs are generated from free anticommuting fields such as the symplectic fermions . Such LCFTs have an interesting fermionic structure where logarithms may also arise in correlation functions involving spin zero anticommuting fields. This is in contrast to free bosons, which typically do not directly appear in the conformal field theory, but only in form of derivatives and exponentials of themselves.

In these notes, we generalize LCFT to the case of Jordan cells of arbitrary rank, but we will restrict ourselves to the Virasoro algebra as the chiral symmetry algebra to keep things simple. With some mild assumptions, the generic form of 1-, 2- and 3-pt functions can be given such that only the structure constants remain as free parameters. From this, the general structure of the OPE as well as some sort of selection rules that a general correlation function may be non-zero, are derived. The crucial role of zero modes, in particular in the case of a LCFT generated from fermionic fields, is emphasized. The results presented here, together with proofs and further details, can be found in . The computation of 4-pt and higher-point functions in the LCFT case is, unfortunatley, more complicated. The interested reader might consult  for some discussion on this issue.

## 2 1-, 2- and 3-pt functions

Let $`r`$ denote the rank of the Jordan cells we consider. One can show, that in LCFTs with Jordan cells with respect to (at least) the $`L_0`$ mode, the $`h=0`$ sector necessarily must carry such a Jordan cell structure. Furthermore, its rank defines the maximal possible rank of all Jordan cells. Thus, without loss of generality, we can assume that the rank of all Jordan cells is equal to $`r`$, other cases can easily be obtained by setting certain structure constants to zero. Each Jordan cell contains one proper highest weight state giving rise to one proper irreducible subrepresentation. We will label this state for a Jordan cell with conformal weight $`h`$ by $`|h;0`$. We choose a basis in the Jordan cell with states $`|h;k`$, $`k=0,\mathrm{},r1`$, such that eq. (2) is replaced by

$$L_0|h;k=h|h;k+|h;k1\text{for}k=1,\mathrm{},r1,L_0|h;0=h|h;0.$$
(3)
The corresponding fields will be denoted $`\mathrm{\Psi }_{(h;k)}`$. Although the OPE of two primary fields might produce logarithmic fields, we will further assume, that primary fields which are members of Jordan cells are proper primaries in the sense that OPEs among them only yield again primaries.

As discussed by Rohsiepe , the possible structures of indecomposable representations with respect to the Virasoro algebra are surprisingly rich. Besides the defining condition eq. (3), further conditions have to be employed to fix the structure. The simplest case is defined via the additional requirement

$$L_1|h;k=0,0k<r.$$
(4)
This condition means that all fields spanning the Jordan cell are quasi-primary. It will be our starting point in the following. Under these assumptions, as shown in , the action of the Virasoro modes receives an additional non-diagonal term. The off-diagonal action is defined via $`\widehat{\delta }_{h_i}\mathrm{\Psi }_{(h_j;k_j)}\left(z\right)=\delta _{ij}\mathrm{\Psi }_{(h_j;k_j1)}\left(z\right)`$ for $`k_j>0`$ and $`\widehat{\delta }_{h_i}\mathrm{\Psi }_{(h_j;0)}\left(z\right)=0`$. Thus,

$$L_n\mathrm{\Psi }_{(h_1;k_1)}\left(z_1\right)\mathrm{}\mathrm{\Psi }_{(h_n;k_n)}\left(z_n\right)=\underset{i}{}z_i^n\left[z_i_i+\left(n+1\right)\left(h_i+\widehat{\delta }_{h_i}\right)\right]\mathrm{\Psi }_{(h_1;k_1)}\left(z_1\right)\mathrm{}\mathrm{\Psi }_{(h_n;k_n)}\left(z_n\right)$$
(5)
for $`n`$. Only the generators $`L_1`$, $`L_0`$, and $`L_1`$ of the Möbius group admit globally valid conservation laws, which usually are expressed in terms of the so-called conformal Ward identities

$$0=\{\begin{array}{ccc}\hfill L_1G(z_1,\mathrm{}z_n)& =& _i_iG(z_1,\mathrm{}z_n),\hfill \\ \hfill L_0G(z_1,\mathrm{}z_n)& =& _i\left(z_i_i+h_i+\widehat{\delta }_{h_i}\right)G(z_1,\mathrm{}z_n),\hfill \\ \hfill L_1G(z_1,\mathrm{}z_n)& =& _i\left(z_i^2_i+2z_i\left[h_i+\widehat{\delta }_{h_i}\right]\right)G(z_1,\mathrm{}z_n),\hfill \end{array}$$
(6)
where $`G(z_1,\mathrm{}z_n)`$ denotes an arbitrary $`n`$-point function $`\mathrm{\Psi }_{(h_1;k_1)}\left(z_1\right)\mathrm{}\mathrm{\Psi }_{(h_n;k_n)}\left(z_n\right)`$ of primary fields and/or their logarithmic partner fields. Here, we already have written down the Ward identities in the form valid for proper Jordan cells in logarithmic conformal field theories. Note that these are now inhomogeneous equations. In principle, this is all one needs to compute the generic form of all $`n`$-pt functions, $`n3`$ upto structure constants. Thus, using freedom of scaling the fields, the 1-pt functions turn out to be

$$\mathrm{\Psi }_{(h;k)}=\delta _{h,0}\delta _{k,r1}.$$
(7)
The 2-pt and 3-pt functions can be written in a rather compact form by noting that derivaties of $`z^h`$ with respect to $`h`$ yields a logarithm. The structure constants depend on both, the conformal weights as well as the total level within the Jordan cells. One obtains

$$\mathrm{\Psi }_{(h;k)}\left(z\right)\mathrm{\Psi }_{(h^{},k^{})}\left(z^{}\right)=\underset{j=r1}{\overset{k+k^{}}{}}D_{(h;j)}\delta _{h,h^{}}\underset{\genfrac{}{}{0pt}{}{0ik,0i^{}k^{}}{i+i^{}=k+k^{}j}}{}\frac{1}{i!i^{}!}\left(_h\right)^i\left(_h^{}\right)^i^{}\left(zz^{}\right)^{hh^{}}$$
(8)
for the 2-pt functions, and for the 3-pt functions analogously

$`\mathrm{\Psi }_{(h_1;k_1)}\left(z_1\right)\mathrm{\Psi }_{(h_2,k_2)}\left(z_2\right)\mathrm{\Psi }_{(h_3;k_3)}\left(z_3\right)`$ $`=`$ $`{\displaystyle \underset{j=r1}{\overset{k_1+k_2+k_3}{}}}C_{(h_1,h_2,h_3;j)}{\displaystyle \underset{\genfrac{}{}{0pt}{}{0i_lk_l,l=1,2,3}{i_1+i_2+i_3=k_1+k_2+k_3j}}{}}{\displaystyle \frac{1}{i_1!i_2!i_3!}}`$
$`\times \left(_{h_1}\right)^{i_1}\left(_{h_2}\right)^{i_2}\left(_{h_3}\right)^{i_3}`$ $`\left(z_{12}\right)^{h_3h_1h_2}\left(z_{13}\right)^{h_2h_1h_3}\left(z_{23}\right)^{h_1h_2h_3}.`$ (9)
## 3 OPEs

It is now a simple matter to write down the generic form of OPEs. In essence, we have to raise one index of the 3-pt structure constants with the help of the inverse of the 2-pt structure constants, i.e. the propagators. Now, in the LCFT case, we have matrices instead, namely

$$\left(D_{h,h^{}}\right)_{k,k^{}}\delta _{h,h^{}}\mathrm{\Psi }_{(h;k)}\left(z_2\right)\mathrm{\Psi }_{(h^{},k^{})}\left(z_3\right),$$
(10)
which is an upper triangular matrix and thus invertible, and

$$\left(C_{(h_1;k1),h_2,h_3}\right)_{k_2,k_3}\mathrm{\Psi }_{(h_1;k_1)}\left(z_1\right)\mathrm{\Psi }_{(h_2,k_2)}\left(z_2\right)\mathrm{\Psi }_{(h_3;k_3)}\left(z_3\right).$$
(11)
The OPE is then given by the expression

$$\mathrm{\Psi }_{(h_1;k_1)}\left(z_1\right)\mathrm{\Psi }_{(h_2;k_2)}\left(z_2\right)=\underset{z_1z_2}{lim}\underset{(h_3;k_3)}{}\underset{k}{}\left(C_{(h_1;k_1),h_2,h_3}\right)_{k_2,k}\left(\left(D_{h_3,h_3}\right)^1\right)^{k,k_3}\mathrm{\Psi }_{(h_3;k_3)}\left(z_2\right),$$
(12)
where the limit means that we have to replace $`z_{13}`$ in the result by $`z_{23}`$ which, in fact, will cancel all dependency on $`z_3`$. In this form, the OPE does not obey a bound such as $`k_3k_1+k_2`$ for the so-called J-levels within the Jordan blocks. For example, pre-logarithmic fields are good primary fields, such that $`k_1=k_2=0`$, while there appears a term with $`k_3=1`$ on the right hand side.

A better bound is given by the zero mode content of the fields. This means the following: The basic fields of the conformal field theory might contain a certain number of zero modes $`\theta _0^{\left(\alpha \right)}`$ such that $`0\left|\theta _0^{\left(\alpha \right)}=\theta _0^{\left(\alpha \right)}\right|0=0`$. These modes will come with canonical conjugate modes $`\xi _{\left(\alpha \right)}`$, which are creators to the right as well as to the left. Thus, the zero mode content $`Z_0\left(\mathrm{\Psi }\right)`$ of a field $`\mathrm{\Psi }`$ is defined as the total number of $`\xi _{\left(\alpha \right)}`$ modes in its mode expansion, expressed in the modes of the basic fields. If the basic fields are anticommuting fermions, we will have anticommuting pairs $`\xi _{\left(\alpha \right)}^\pm `$ instead such that we can define zero mode contents $`Z_+\left(\mathrm{\Psi }\right)`$ and $`Z_{}\left(\mathrm{\Psi }\right)`$ separately with $`Z_0=Z_++Z_{}`$. Explicitly known examples of LCFTs do have realizations in fermionic free fields, and it turns out that the definition above can be extended to pre-logarithmic fields in a consistent way by assigning them fractional values $`Z_+`$ and $`Z_{}`$ such that always $`Z_0`$. In fact, a large class of LCFTs can be constructed from ordinary conformal field theories by introducing additional zero modes accompanied with a suitable deformation of the Virasoro modes, see  for details. The zero mode content does now provide a bound for OPEs, namely

$$Z_0\left(\mathrm{\Psi }_{(h_3;k_3)}\right)Z_0\left(\mathrm{\Psi }_{(h_1;k_1)}\right)+Z_0\left(\mathrm{\Psi }_{(h_2;k_2)}\right).$$
(13)
One of the best known examples for a LCFT is the $`c=2`$ ghost system, written in terms of two anticommuting spin zero fields $`\theta ^\pm \left(z\right)`$. The mode expansion reads

$$\theta ^\pm \left(z\right)=\xi ^\pm +\theta _0^\pm \mathrm{log}\left(z\right)+\underset{n0}{}\theta _n^\pm z^n,$$
(14)
where the modes $`\xi ^\pm `$ are the creator zero modes, while the modes $`\theta _0^\pm `$ are the annihilator zero modes, satisfying $`\{\xi ^\pm ,\theta _0^{}\}=1`$, $`\{\theta _n^+,\theta _m^{}\}=\frac{1}{n}\delta _{n+m,0}`$. We get back the original $`bc`$ ghost system by setting $`c\left(z\right)=\theta ^{}\left(z\right)|_{\theta _0^{}=0}`$ and $`b\left(z\right)=_z\theta ^+\left(z\right)`$. Thus, the pair $`\theta _0^{}`$ and $`\xi ^+`$ of zero modes is absent in the $`bc`$ system, and so is the logarithmic partner of the identity field, $`\mathrm{\Psi }_{\left(h=0;1\right)}\left(z\right)=\text{:}\theta ^+\theta ^{}\text{:}\left(z\right)`$ with state $`|0;1=\xi ^+\xi ^{}|0;0`$, where $`|0;0=|0`$. Thus, $`Z_0\left(\mathrm{\Psi }_{(0;0)}\right)=2`$, $`Z_\pm \left(\mathrm{\Psi }_{(0;0)}\right)=1`$, while the basic fermionic fields obey $`Z_0\left(\theta ^\pm \right)=1`$, $`Z_\pm \left(\theta ^\pm \right)=1`$ and $`Z_{}\left(\theta ^\pm \right)=0`$. This theory possesses a pre-logarithmic field $`\mu `$ of conformal weight $`h=1/8`$ with OPE $`\mu \left(z\right)\mu \left(0\right)z^{1/4}\left(\mathrm{\Psi }_{(0;1)}\left(0\right)2\mathrm{log}\left(z\right)\mathrm{\Psi }_{(0;0)}\left(0\right)\right)A+z^{1/4}\mathrm{\Psi }_{(0;0)}\left(0\right)B`$. To make everything consistent, one assigns $`Z_\pm \left(\mu \right)=1/2`$. We mention for completeness, that the excited twist field $`\sigma `$ with conformal weight $`h=3/8`$ has to be assigned the values $`Z_+\left(\sigma \right)=3/2`$ and $`Z_{}\left(\sigma \right)=1/2`$ or vice versa.

## 4 Zero mode content

Let us briefly consider a much less trivial example, the ghost system with $`c=26`$, made out of a pair of anticommuting fiels of spin $`2`$ and $`1`$, respectivelyi . In general, the $`(j,1j)`$ ghost system possesses $`2j1`$ zero modes $`b_{j1},b_{j2},\mathrm{},b_{1j}`$. The stress energy tensor reads

$$T_{bc}=j\text{:}b\left(c\right)\text{:}+\left(1j\right)\text{:}\left(b\right)c\text{:}.$$
(15)
Using a slight generalization of the deformation technique of , additional zero modes can be introduced by a modification, shown here for the $`j=2`$ case,

$$T_{\mathrm{log}}\left(z\right)=T_{bc}\left(z\right)+A\theta _1^{}b\left(z\right)+B\theta _0^{}z^1\left(z^2b\left(z\right)\right)+A\theta _1^{}z^2\left(z^4b\left(z\right)\right).$$
(16)
These additional zero modes can be thought of as modes of $`h=1`$ fields $`\theta ^\pm `$ with expansion

$$\theta ^\pm \left(z\right)=\xi _1^\pm z^2+\xi _0^\pm z+\xi _{+1}^\pm +\theta _1^\pm \frac{z^2}{2}\left(\mathrm{log}\left(z\right)\frac{3}{2}\right)+\theta _0^\pm z\left(\mathrm{log}\left(z\right)1\right)+\theta _{+1}^\pm \mathrm{log}\left(z\right)+\underset{\left|n\right|>1}{}\theta _n^\pm \frac{z^{n+1}}{1n}.$$
(17)
Again, $`b\left(z\right)=^{2j1}\theta ^+\left(z\right)`$ and $`c\left(z\right)=\theta ^{}\left(z\right)|_{\theta _{j1}^{}=\theta _{j2}^{}=\mathrm{}=\theta _{1j}^{}=0}`$ such that the $`\theta ^\pm `$ fields have twice as many zero modes as the original $`bc`$ system, and $`\{\xi _i^\pm ,\theta _i^{}\}=\pm \left(1\right)^{i+1}`$. Although the modes of the modified stress energy tensor satisfy the Virasoro algebra, they do not act consistently on the space of states, e.g. $`L_0|\xi _1^+=0`$.

However, as explained in , the doubling of the zero modes is not completely articifial, but does naturally imply that the conformal field theory (CFT) now lives on a hyperellptic Riemann surface, viewed as a double covering of the complex plane or Riemann sphere. Thus, we actually have a CFT on each of the sheets, such that the full CFT is the tensor product of the individual ones with $`T_{\mathrm{log}}=T_{\mathrm{log}}^{\left(1\right)}+T_{\mathrm{log}}^{\left(2\right)}`$ such that $`[T_{\mathrm{log}}^{\left(1\right)},T_{\mathrm{log}}^{\left(2\right)}]=0`$. In fact, this is possible and yields a consistent CFT provided we identify the zero modes on the different sheets with each other as $`\theta _i^{\left(1\right),\pm }=\theta _i^{\left(2\right),}`$ and $`\xi _i^{\left(1\right),\pm }=\xi _i^{\left(2\right),}`$ for $`i=1,0,1`$. This yields the Virasoro algebra for $`T_{\mathrm{log}}`$ with total central charge $`c=52`$ and a correct action of its modes on the space of states. The resulting theory possesses indecomposable representations despite the fact that the action of $`L_0`$ remains diagonal. The construction generalizes to other ghost systems, but is is not yet clear, how the construction works for higher ramified covering with more than two sheets.

This shows that zero modes are at the heart of LCFTs. Furthermore, the zero mode content provides strong conditions on whether correlation functions can actually be non-zero. It appears that fields $`\mathrm{\Psi }_{(h;k)}`$ forming Jordan blocks have a well defined even zero mode content $`Z_0=Z_++Z_{}`$ with $`Z_+=Z_{}`$. Fermionic fields, in turn, are characterized by $`Z_+Z_{}`$ but still satisfy $`Z_\pm `$. These fields are denoted by $`\mathrm{\Theta }_{(h;k^+,k^+)}`$. There are no examples known where such fields do form indecomposable structures, but this is only due to the fact that no LCFTs with a sufficiently high number of genuine zero modes have been explicitly examined yet. Finally, pre-logarithmic, or more generally, twist fields, have fractional zero mode contents, $`Z_\pm `$, and are denoted by $`\mu _\alpha `$. Such fields are generally believed to reside in irreducible representations. A generic correlation function is then of the form

$$G=\underset{i}{}\mathrm{\Psi }_{(h_i;k_i)}\left(z_i\right)\underset{j}{}\mathrm{\Theta }_{(h_j;k_k^+,k_j^{})}\left(w_j\right)\underset{l}{}\mu _{\alpha _l}\left(u_l\right)\underset{i}{}\mathrm{\Psi }_i\underset{j}{}\mathrm{\Theta }_j\underset{l}{}\mu _l.$$
(18)
The zero mode contents implies now that $`G=0`$ unless all three conditions

$`_+Z_0\left(G\right)={\displaystyle \underset{i}{}}Z_0\left(\mathrm{\Psi }_i\right)+{\displaystyle \underset{j}{}}Z_0\left(\mathrm{\Theta }_j\right)+{\displaystyle \underset{l}{}}Z_0\left(\mu _l\right)2\left(r1\right),`$ (19)
$`{\displaystyle \underset{j}{}}Z_+\left(\mathrm{\Theta }_j\right)={\displaystyle \underset{j}{}}Z_{}\left(\mathrm{\Theta }_j\right),`$ (20)
$`{\displaystyle \underset{l}{}}Z_+\left(\mu _l\right)\text{and}{\displaystyle \underset{l}{}}Z_{}\left(\mu _l\right)`$ (21)
are satisfied. These are very powerfull statements since they imply further that we can relax our condition that the logarithmic partners have to be quasi-primary. In fact, the action of $`L_n`$, $`n=1,0,1`$, in the Ward identities eq. (6) will yield new correlation functions, $`L_nG=_kG_k^{}`$. If there are contributions from fields failing to be quasi-primary, then these can be neglected if the resulting correlation functions $`G^{}`$ do not anylonger satisfy eqs (1921). Since $`L_n`$ act as derivations, only one field in the correlator is modified in each term. Thus, if the zero mode content of a non-quasi-primary term differs from the original zero mode content such that the balance is broken, it will not contribute to the correlation function $`G`$ since it does not affect the Ward identities. Hence, $`L_n`$ implements a BRST like structure on the complex spanned by $`Z_+`$ and $`Z_{}`$, as anticipated in .

In summary, we have provided the general structure of correlation functions and OPEs for LCFT with arbitrary high rank Jordan cells. We found strong constraints for correlation functions to be non-zero, intimately linked to the zero mode content of the involved fields. It remains an open problem, however, what the modular properties of such higher rank LCFTs are. These are only knonw in the rank two case .

Acknowledgement The research of M. F. is supported by the European Union network HPRN-CT-2002-00325 and the research of M. F. and M. K. is supported by the string theory network (SPP no. 1096), Fl 259/2-2, of the Deutsche Forschungsgemeinschaft.