text:
# MD simulations and continuum theory of partially fluidized shear granular flows

## Abstract

We carry out a detailed comparison of soft particle molecular dynamics simulations with the theory of partially fluidized shear granular flows. We verify by direct simulations a constitutive relation based on the separation of the shear stress tensor into a fluid part proportional to the strain rate tensor, and a remaining solid part. The ratio of these two components is determined by the order parameter. Based on results of the simulations we construct the “free energy” function for the order parameter. We also present the simulations of the stationary deep 2D granular flows driven by an upper wall and compare it with the continuum theory.

When the ratio of shear to normal stress in a packed granular matter exceeds a certain threshold value, the granular matter yields and a flow ensues. In the last few years there have been many experimental, numerical and theoretical studies behringer99 ; gollub98 ; gollub00 ; lemaitre ; starton ; pouliquen99 ; daerr ; Tsai02 ; grest91 ; wolf96 ; grest02 ; sparks99  that explored a broad range of granular flow conditions. While dilute granular flows can be well described by the kinetic theory of dissipative granular gases Jenkins85 , dense granular flows still present significant difficulties for theoretical description. A continuum theory of slow dense granular flows based on the so-called associated flow rule that relates the strain rate and the shear stress was proposed in Ref. savage98 . This description neglects the effects of the dry friction between the grains and works only in a fluidized state, so it cannot not describe hysteretic nature of the granular flow relevant for stick-slips gollub98 , avalanching daerr , etc. A similar model based on a Newtonian stress-strain constitutive relation with density dependent viscosity was proposed in Refs. gollub00 . In this model, the viscosity diverges at the fluidization threshold when the density approaches the random close packing density of grains.

Recently we proposed a phenomenological order parameter (OP) description of the fluidization transition AT1 . The OP specifies the ratio between solid and fluid parts of the stress tensor. The viscosity is defined as a ratio of the fluid part of the shear stress to the strain rate and remains finite at the fluidization threshold. This model yielded a good qualitative description of broad variety of phenomena occurring in granular flows.

In this Letter we report on 2D soft particle molecular dynamics simulations performed to validate and quantify our OP theory. To fit the equation for the OP and stress-strain relation we performed simulations of the granular flow in a thin Couette geometry. The obtained set of equations was used to calculate the stress and velocity profiles in a different system, a thick granular layer under non-zero gravity driven by a moving heavy upper plate.

The theory AT1  is based on a standard momentum conservation equation and the incompressibility condition applicable for slow dense flows. To close the system, we assumed that the stress tensor $`\sigma `$ is comprised of two parts, a solid part $`\sigma ^s`$, and a fluid part $`\sigma ^f`$ (taken in a purely Newtonian form)

$$\sigma _{ij}^f=p_f\delta _{ij}\mu _f\dot{\gamma }_{ij}$$
(1)
where $`p_f`$ is the isotropic “partial” fluid pressure, $`\mu _f`$ is the viscosity coefficient associated with the fluid stress tensor. We set the fluid part of off-diagonal components of the stress tensor to be proportional to the off-diagonal components of the full stress tensor with the proportionality coefficient being a function of the OP $`\rho `$,

$$\sigma _{yx}^f=q(\rho )\sigma _{yx};\sigma _{yx}^s=(1q(\rho ))\sigma _{yx}.$$
(2)
We choose a fixed range for the OP such that it is zero in a completely fluidized state and one in a completely static regime. Thus, the function $`q(\rho )`$ has the property $`q(0)=1,q(1)=0`$. In Refs.AT1  for simplicity we postulated $`q(\rho )=1\rho `$. A similar relationship can be introduced for the diagonal terms of the stress tensor

$`\sigma _{xx}^f`$ $`=`$ $`q_x(\rho )\sigma _{xx},\sigma _{yy}^f=q_y(\rho )\sigma _{yy}`$ (3)
$`\sigma _{xx}^s`$ $`=`$ $`(1q_x(\rho ))\sigma _{xx},\sigma _{yy}^s=(1q_y(\rho ))\sigma _{yy}`$ (4)
The dynamics of the OP was assumed to be relaxational in nature and controlled by the generic Ginzburg-Landau equation,

$$\frac{D\rho }{Dt}=D^2\rho F(\rho ,\delta )$$
(5)
Here $`D/Dt`$ is the material derivative, $`D`$ is the diffusion coefficient, $`F(\rho ,\delta )`$ is the derivative of the free energy density which has a quartic polynomial form to account for the bistability near the solid-fluid transition. Control parameter $`\delta `$ was taken to be a linear function of $`\varphi =\mathrm{max}|\sigma _{mn}/\sigma _{nn}|`$.

The OP $`\rho `$ which plays a pivotal role in the theory, should be associated with the ”microscopic” properties of the granular assembly. At any moment of time all contacts among the grains can be classified as either fluid-like or solid-like. A contact is considered fluid-like if two particles slide past each other or briefly collide, and is solid-like if two particles are jammed together for longer than a characteristic collision time. Here we postulate that the OP can be introduced as a ratio between space-time averaged numbers of solid contacts $`\overline{Z_s}`$ and all contacts $`\overline{Z}`$ within a sampling area,

$$\rho (y)=\overline{Z_s}/\overline{Z}.$$
(6)
where $`\xi `$ and $`\overline{\xi }`$ stand for averaging of $`\xi `$ in space and time respectively. This definition satisfies both limiting cases: when a granulate is in a solid state all contacts are stuck and $`\rho =1`$; when it is strongly agitated $`\overline{Z_s}`$ is zero and $`\overline{Z}`$ is small but finite, therefore $`\rho =0`$. Our OP is expected to be sensitive to the degree of fluidization. A small rearrangement of the force network may result in strong fluctuations of $`\rho `$, while the solid fraction and granular temperature remain virtually constant. This quantity is difficult to measure in experiments, however it can be found in soft-particle molecular dynamics.

Molecular Dynamics Simulations note1 . The grains are assumed to be non-cohesive, dry, inelastic disk-like particles. Two grains interact via normal and shear forces whenever they overlap. For the normal impact we employed *spring-dashpot* model wolf96 . This model accounts for repulsion and dissipation; the repulsive component is proportional to the degree of the overlap and the velocity dependent damping component simulates the dissipation. The model for shear force is based upon technique developed in cundall79 . The motion of a grain is obtained by integrating the Newton’s equations with forces and torques produced by interactions with neighbors and walls. For detailed discussion of the advantages and limitations of the employed model see Refs. wolf96 ; grest91 ; grest02 . The computational domain spans $`L_x\times L_y`$ area, and is periodic in horizontal direction $`x`$. The grain diameters are uniformly distributed around mean with relative width $`\mathrm{\Delta }_r`$ to avoid crystallization effects grest91 . The material parameters of grains were chosen similar to Ref.grest02 . All quantities are normalized by an appropriate combination of the average particle diameter $`d`$, mass $`m`$, and gravity $`g`$.

We studied a thin ($`50\times 10`$) granular layer sandwiched between two “rough plates” under fixed external pressure $`P_{ext}`$ and zero gravity conditions. The rough plates were simulated by two straight chains of large grains “glued” together. Two opposite forces $`𝐅_1=𝐅_2`$ were applied to the plates along the horizontal $`x`$ axis to induce shear stress in the bulk. We started with weak forces and slowly ramped them up in small increments. After that we ramped the shear forces down until the granular layer was jammed again. At every “stop” we measured all stress components, strain rate, and the OP by averaging over the whole layer and over time. These simulations were carried out at several values of the external pressure $`P_{ext}`$. Figure 1 shows $`\rho `$ as a function of the normalized shear stress $`\delta =\sigma _{yx}/P_{ext}`$. With this normalization, the results of several simulations with different $`P_{ext}`$ collapse on a single bifurcation curve. We observe a quiescent state $`\rho =1`$ until $`\delta `$ reaches a certain critical value $`\delta _10.32`$. This value differs slightly for different runs because of the finite system size and absence of self-averaging in the static regime. Above $`\delta _1`$, $`\rho `$ abruptly drops to approximately $`0.15`$. At larger $`\delta `$, the OP rapidly approaches zero. The return curve corresponding to the diminishing of the shear stress follows roughly the same path, and then continues to another (smaller) value of the shear stress, $`\delta _20.23`$. At this point the OP jumps back to one, and the granular layer returns to a jammed state.

The striking feature of this bifurcation diagram is the hysteretic behavior of the OP as a function of the shear stress. This hysteresis was anticipated in our model AT1 , and now we are in a position to describe it quantitatively. Assuming that there is an (unobserved) unstable branch of the bifurcation curve which merges with the stable branch at $`\delta \delta _1`$, we make a simple analytic fit,

$$F(\rho ,\delta )=(1\rho )\left(\rho ^22\rho _{}\rho +\rho _{}^2\mathrm{exp}[A(\delta ^2\delta _{}^2)]\right)$$
(7)
with $`\rho _{}=0.6,A=25,\delta _{}=0.26`$ (see Figure 1, dashed line) and use it in equation (5). The inset of Figure 1 depicts the density $`\nu `$ vs. the OP $`\rho `$ for the same runs. The density stays almost constant in a wide range of the OP $`0.1<\rho <1`$. This shows that unlike the particle density, our OP is a sensitive characteristic of slow dense granular flows reflecting subtle changes in the contact network and the structure of the stress distribution.

We also probed the relaxation dynamics of the OP by studying the response of the system on small perturbation in the the hysteretic region, $`\delta _2<\delta <\delta _1`$ note1 . From these simulations we find that the intrinsic time scale of the OP relaxation is rather small, $`O(1)`$. Our thin Couette flow system did not allow us to probe the local coupling of the OP since $`\rho const`$ throughout the system. In the absence of such data this coupling was modeled by the linear diffusion term in (5) with $`D=const`$.

The constitutive relation was fitted using the same Couette flow simulations. We analyzed the fluid stress $`\sigma _{ij}^f`$ and the solid stress $`\sigma _{ij}^s`$ separately during our ramp-down simulations at three different values of $`P`$. Figure 2 shows the ratios of fluid tensor components to the corresponding full tensor components as functions of $`\rho `$. We observe that data for $`\sigma _{xy}^s/\sigma _{xy}`$ and $`\sigma _{yx}^s/\sigma _{yx}`$ for different $`P`$ fall onto a single curve which is fitted by $`q(\rho )=(1\rho )^{2.5}`$ (Figure 2,a). Repeating the procedure with diagonal elements yields different scaling, see Figure 2, b. A small but noticeable difference is evident between $`\sigma _{xx}^f/\sigma _{xx}`$ and $`\sigma _{yy}^f/\sigma _{yy}`$. Detailed analysis shows that in fact fluid parts of the diagonal components of the stress tensor $`\sigma _{xx}^f`$ and $`\sigma _{yy}^f`$ are nearly identical note2 , and the difference is mainly due to the solid part of the normal stresses. Functions $`q_{x,y}(\rho )`$ approach 1 as $`\rho 0`$, but they may have different functional form to reflect the anisotropy of the solid stress tensor. In our Couette flow, $`\sigma _{xx}`$ and $`\sigma _{yy}`$ can be fitted by $`q_x(\rho )(1\rho )^{1.9}`$ and $`q_y(\rho )(1\rho ^{1.2})^{1.9}`$, respectively, see Figure 2,b. We observe that even in a partially fluidized regime, the “fluid phase” component behaves as a fluid with a well-defined isotropic “partial” pressure $`p_f`$ which is zero in a solid state and is becoming the full pressure in a completely fluidized state.

To test the stress-strain relation (1) we plot $`\sigma _{yx}^f`$ vs $`\dot{\gamma }_{yx}`$, see Figure 3. At small $`\dot{\gamma }_{yx}`$ all curves are close to the same straight line $`\sigma _{yx}^f12\dot{\gamma }_{yx}`$, i.e. the Newtonian scaling for fluid shear stress holds reasonably well. The deviations at large $`\dot{\gamma }_{yx}`$ are evidently caused by variations of density and temperature in the dilute regime as to be expected from kinetic theory of dilute granular flows Jenkins85 . The full shear stress, of course, does not vanish as $`\dot{\gamma }_{yx}0`$ (Figure 3, inset). Therefore, the standard viscosity coefficient defined as the ratio of the full shear stress and strain rate diverges at the fluidization threshold as observed in Ref. gollub00 .

We applied our theoretical description which was formulated above on the basis of MD simulations of a thin Couette flow with no gravity, to a different system, a shear granular flow in a thick granular layer under gravity driven by the upper plate which was pulled with constant speed $`V`$ (or constant force $`F`$). A similar system has been studied experimentally Refs. gollub98 ; Tsai02 . We simulated up to 20,000 particles in a periodic rectangular box under a heavy plate which was moved horizontally in $`x`$-direction.

We systematically carried out comparison between MD simulations and the continuum theory using the stress-strain relations and the specific form of the OP equation described above. We used no-flux boundary conditions for the OP and no-slip condition for the velocity at the bottom plate. The boundary condition at the top plate where the flow can be in a dilute regime, is a separate issue which we do not address here. We assumed that the shear stress component in the bulk is specified by applied force, and calculated the velocity profile using the constitutive relation (1)-(4). We found that the constitutive relations determined from the thin Couette experiment hold for this system as well. Selected results for the OP and velocity profiles are presented in Figure 4. As seen from the Figure, the vertical profiles of the OP and the horizontal velocities are reasonably well described by the theory. However, for low pressure the horizontal velocity profiles deviate from the numerical data for low pressure runs, apparently because the viscosity coefficient is no longer a constant in a dilute region near the top plate. The only fitting parameter used was the diffusion constant $`D`$ in the OP equation, which has not been determined in our simulations of the thin layer. It appears that the diffusion coefficient depends on the applied pressure and the strain rate, however more detailed numerical experiments are needed.

In conclusion, we calibrated the theory of partially granular flows AT1  on the basis of a series of 2D soft particle molecular dynamics simulations. The OP which controls the fluidization transition, was defined as the fraction of solid-like contacts among particles. Measurements of the OP, the stress tensor, and the strain rate in a thin Couette cell allowed us to quantify the constitutive relations based on the relaxational OP dynamics. We studied the flow structure of a thick surface driven granular flow under gravity and found the model predictions to be in a good quantitative agreement with soft-particle MD simulations. Our results support an intriguing interpretation for the OP description of dense and slow granular flows. The granular matter under shear stress appears to be similar to a multi-phase system with the fluid phase “immersed” in the solid phase. The fluid phase behaves as a simple Newtonian fluid for small shear rates when the density is almost constant, but exhibits shear thinning at larger shear rates when the density begins to drop. This regime can be described by the generalization of the theory note1  which includes the equations for density and granular temperature following from the kinetic theory of dilute granular gases. Our simulations were limited to 2D geometry. While we anticipate that the structure of the model should remain unchanged in 3D systems, the specific form of the fitting functions may vary. The authors are indebted to J. Gollub, P. Cvitanovic, T. Halsey, B. Beringer for useful discussions, and to J.C. Tsai for sharing his unpublished experimental data. This work was supported by the Office of the Basic Energy Sciences at the US DOE, grants W-31-109-ENG-38, and DE-FG03-95ER14516. Simulation were performed at the National Energy Research Scientific Computing Center.


format: text

text:
# Influence of temperature on the nodal properties of the longitudinal thermal conductivity of YBa2Cu3O7-x

## I Introduction

The order parameter in high-temperature superconductors (HTS) has been proved to have mainly a $`d_{x^2y^2}`$ pairing symmetry.Tsuei and Kirtley (2000); Aubin et al. (1997); Salamon et al. (1995); Ocaña and Esquinazi (2002) In particular, try-crystal phase-sensitive measurements Tsuei and Kirtley (2000) have determined the presence of a predominant $`d_{x^2y^2}`$gap symmetry up to the superconducting critical temperature $`T_c`$ over, if any, other minor components smaller than 5% of that with the $`d_{x^2y^2}`$symmetry. Thermal transport measurements have also shown the predominance of this symmetry for the gap function. In particular, when a magnetic field is rotated parallel to the CuO2 planes, the longitudinal thermal conductivity shows a fourfold oscillation which can be explained in terms of both Andreev scattering of quasiparticles by vortices (AS) and Doppler shift (DS) in the energy spectrum of the quasiparticles if one takes a $`d_{x^2y^2}`$gap into account.Aubin et al. (1997); Salamon et al. (1995); Ocaña and Esquinazi (2002) Nevertheless, the fourfold oscillation in the thermal conductivity has been resolved up to $`15`$ K. Above this temperature, no direct evidence for this type of symmetry from this kind of measurements has been published. Furthermore, deviations from the expected angular pattern within a pure $`d_{x^2y^2}`$symmetry have been attributed to the effect of pinning of vortices.Aubin et al. (1997)

The variation of the thermal conductivity as a function of angle $`\theta `$ between the heat current and the magnetic field applied parallel to the CuO2 planes depends mainly on the heat transport by quasiparticles, their interaction with the supercurrents (vortices) and the symmetry of the order parameter.Aubin et al. (1997); Salamon et al. (1995); Ocaña and Esquinazi (2002) An angular pattern showing properties of the order parameter symmetry can only be achieved when the temperature is low enough so that the quasiparticle momentum is close to the nodal directions of the gap. Otherwise, thermal activation would also induce quasiparticles at different orientations from those of the nodes and hence, the sensitivity of the probe to measure gap characteristics will be reduced. The first question we address in this paper is related to the temperature range at which the nodal characteristics of the gap are directly observable by thermal conductivity. The second question we would like to clarify in this paper is whether the thermal activation of quasiparticles with increasing temperature, treated phenomenologically within a Fermi liquid approximation for thermal transport including the $`d_{x^2y^2}`$gap function, can explain the experimental results in the whole temperature range, in particular the change of symmetry of $`\kappa _{xx}(\theta )`$ as a function of temperature. In this paper, we calculate numerically the thermal conductivity at different angles and temperatures at fixed magnetic field, assuming an Andreev reflection model for the scattering of quasiparticles by supercurrents, originally proposed by Yu et al. Yu et al. (1995); Klemm et al. (1996); Yu et al. (1996) within the two dimensional BRT expression Bardeen et al. (1959) for the thermal conductivity, and compare it to the experimental data.

Angle scans in a magnetic field applied parallel to the CuO2 planes were performed in order to measure the angular variation of the longitudinal thermal conductivity $`\kappa _{xx}(\theta )`$ in two single crystals of YBa2Cu3O7-x high-temperature superconductor. The overall results agree with the theoretical model and confirm the predominance of the $`d_{x^2y^2}`$gap up to $`T_c`$. The measurements provide also new results that improve our knowledge of the thermal transport at temperatures at which the nodal properties of the gap are not directly observable.

Following the experimental and sample details of the next section, we present in Sec. III the main experimental results. In Sec. IV we describe the used model and compare it with the experimental data. A brief summary is given in Sec. V.

## II Experimental and Sample details

In order to rule out effects concerning shape and structure characteristics of the crystal we have used two different samples of YBa2Cu3O7-x (YBCO): a twinned single crystal with dimensions (length $`\times `$ width $`\times `$ thickness) $`0.83\times 0.6\times 0.045`$mm3 and critical temperature $`T_c=93.4`$K previously studied in Refs. Taldenkov et al. (1999); Ocaña et al. (2001); Ocaña and Esquinazi (2001, 2002) and an untwinned single crystal with dimensions $`2.02\times 0.68\times 0.14`$mm3 and $`T_c=88`$K.Ocaña and Esquinazi (2001) For the measurement of the thermal conductivity, a heat current $`J`$ was applied along the longest axis of the crystal studied. In the untwinned sample, J was parallel to the $`a`$-axis and in the twinned crystal was parallel to the $`a/b`$-axes (twin planes oriented along (110)). In both cases the position of the lattice axis with respect to the crystal axis was determined using polarized light microscopy and X-ray diffraction.

The longitudinal temperature gradient $`(_xT)`$ was measured using previously calibrated chromel-constantan (type E) thermocouples Inyuskin et al. (1998) and a dc picovoltmeter. Special efforts were made in order to minimize the misalignment of the plane of rotation of the magnetic field applied perpendicular to the c-axis with the CuO2 planes of the sample. This misalignment was minimized step by step, measuring the angle dependence of $`\kappa _{xx}(\theta )`$ until a satisfactory symmetrical curve was obtained. In this way, we estimate a misalignment smaller than $`0.5^0`$. An in-situ rotation system enabled measurement of the thermal conductivity as a function of the angle $`\theta `$ defined between the applied field and the heat flow direction along $`+\widehat{x}`$, see Fig. 1. For more details on the experimental arrangement see Ref. Ocaña and Esquinazi (2002).

As pointed out by Aubin et al. Aubin et al. (1997) and observed in Refs. Taldenkov et al. (1999); Ocaña and Esquinazi (2002), the effect of the pinning of vortices plays an important role in determining the correct angle pattern in this kind of measurement. In fact, when the angle of the magnetic field is changed, a non uniform vortex distribution due to pinning forces may appear. As argued in Ref. Ocaña and Esquinazi (2002), the pinning of the Josephson-like vortices parallel to the planes is strongly affected by vortices perpendicular to the planes which may appear due to the misalignment of the crystal axes with respect to the applied magnetic field. In this situation, even hysteresis in the angular patterns of the thermal conductivity can be measured. Aubin et al. (1997) We note that pinning of vortices is influenced by the distribution of the oxygen vacancies in the sample as well as by defects and impurity centers. Thus, in order to rule out the influence of this effect in the measurements, a field-cooled procedure have to be used. This procedure consists of two steps. In the first step the sample is driven into the normal state by heating to a few Kelvin above $`T_c`$ and the angle is changed. Secondly, it is cooled down to the desired temperature at constant field.

In Fig. 2 we show the longitudinal thermal conductivity as a function of the temperature in both crystals at zero magnetic field. As argued by many authors (see for example Refs. Hirschfeld and Putikka (1996); Hirschfeld ; Kübert and Hirschfeld (1998)), the observed behavior in Fig. 2 provides qualitative information about the quality of the sample. The height of the peak observed in the temperature dependence of the thermal conductivity is related to the relative contributions between the density of impurity scattering centers and the strength of the inelastic electron-electron scattering. Crystals showing large peaks may have a small amount of impurity scattering centers zeiniepj  and/or larger quasiparticle-related inelastic contribution. The origin of the peak is explained in terms of a competition between the decrease of the inelastic scattering rate Bonn et al. (1993) and the decrease of the population of quasiparticles as temperature decreases.Hirschfeld ; Houssa and Ausloss (1995, 1996a) Taking this into account we might conclude that the untwinned crystal has a larger impurity scattering and a smaller quasiparticle-related inelastic scattering since its peak is considerably smaller than that of the twinned sample. This is supported by its reduced Tc with respect to the twinned sample (which indicates a larger density of oxygen vacancies). Since the thermal conductivity at $`T_c`$ appears to be larger for the untwinned sample, we conclude that there should be a substantial reduction of the inelastic scattering. As we show below, this makes the measurement of electronic properties related to the gap symmetry more difficult.

## III Experimental Results and Discussion

The angular patterns of the longitudinal thermal conductivity $`\kappa _{xx}`$ at different temperatures for the twinned and untwinned crystals are shown in Figs. 3 and 4 respectively. At low enough temperatures ($`T<20K`$), a magnetic field parallel to the CuO2 planes of the sample produces a fourfold oscillation in the longitudinal thermal conductivity as the field is rotated from $`\theta =90^{}`$ to $`\theta =90^{}`$. We note that at both $`\theta =90^{}`$ and $`\theta =90^{}`$ the magnetic field is perpendicular to the heat flow. At $`\theta =0^{}`$ the field is parallel to the heat flow (see Fig. 1). As noted in Refs. Aubin et al. (1997) and Ocaña and Esquinazi (2002) the angular patterns are not free of the phononic contribution to the total thermal conductivity. However, this contribution seems to be constant under variations of the magnetic field orientations for two reasons. 1) For fields parallel to the CuO2 planes, vortices are unlikely to have a normal core (Josephson vortices), making a change of the phonon attenuation with field and angle unlikely. 2) There is no experimental evidence that the phononic contribution does change substantially with field. The similarity between the temperature dependencies of $`\kappa _{xx}`$ and $`\kappa _{xy}`$ below $`T_c`$ where a large change in the quasiparticle density ocurrsOcaña et al. (2001), indicates that the phonon-electron interaction does not play a main role in the temperature dependence of $`\kappa `$. The electronic contribution is entirely responsible for the magnetic field dependence, revealed in both thermal conductivity Yu et al. (1992); Zeini et al. (1999); Krishana et al. (1999); Ocaña et al. (2001) and in microwave measurements Bonn et al. (1993), and indicates against any large phonon-electron scattering. In fact, no remarkable phonon-electron scattering have been addressed in those experiments where the measured magnetic field dependencies are larger than the oscillation amplitudes measured in this paper. Therefore, measurable contributions of the phonons via phonon-electron scattering are more unlikely to occur in the angular patterns. This characteristic makes the angular profiles of the thermal conductivity suitable to be compared with electronic models by using the quantity $`\kappa _{xx}(\theta )\kappa _{xx}(90^{})`$ where the phononic contribution is subtractedOcaña and Esquinazi (2002).

As pointed out in Ref. Ocaña and Esquinazi (2002), the variation of the thermal conductivity in Figs. 3 and 4 can be explained with a model involving a Doppler shift (DS) in the energy spectrum of the quasiparticles along with an accurate inclusion of the impurity scattering Won and Maki ; Won and Maki (2001); Hirschfeld ; Kübert and Hirschfeld (1998) and/or assuming Andreev scattering (AS) of quasiparticles by vorticesYu et al. (1995). Briefly, in the mixed state the quasiparticles are in the presence of a phase gradient produced by the superfluid flow of the vortices. Thus, as viewed from the laboratory frame they experience a Doppler shift in their energy spectrum given by the scalar product of the momentum and superfluid velocity, p$``$vs. When a quasiparticle of momentum p is moving parallel to the magnetic field, the product is zero and hence, no DS occurs. Thus, the angular characteristic of this effect is to produce an excess of quasiparticles in the direction perpendicular to the field and thereby, reducing the “local” thermal resistance at this orientation Won and Maki ; Won and Maki (2001); Hirschfeld . When the field is placed parallel to the heat current, the Doppler shift is the same for both nodal directions at $`\theta =45^{}`$ and $`\theta =45^{}`$. Therefore and since thermal resistances must be added in parallel, a simple picture in which there are only quasiparticles at the nodes would produce a fourfold oscillation in the thermal conductivity with opposite sign to that observed in the measurements. However, as pointed out in Refs. Kübert and Hirschfeld (1998); Hirschfeld ; Won and Maki (2001), the DS affects both the carrier density as well as the scattering rate of quasiparticles, and since at high enough temperatures the latter dominates, the quasiparticles move more easily parallel to the magnetic field. Therefore, the sign of the fourfold oscillation observed in the experiments is also recovered within this picture.

In the AS mechanism a phase gradient, namely, the superfluid flow surrounding the vortices may induce Andreev reflection of the quasiparticles Andreev (1964). Thus, the DS in the energy spectrum p$``$vs is implicitly taken into account in the AS picture. As viewed from the laboratory frame, when the quasiparticle energy equals the value of the gap, then the quasiparticle is transformed into a quasihole reversing its velocity and hence, decreasing its contribution to the thermal conductivity. Thus, a quasiparticle with momentum p parallel to the magnetic field does not experience a DS and hence, no Andreev reflection can take place. The field acts as a filter Matulis et al. (1994) for the quasiparticle that contributes to reduce the total temperature gradient. Note that the AS picture gives rise to a fourfold oscillation in the thermal conductivity as the magnetic field is rotated parallel to the CuO2 planes without more considerations.

Then, qualitatively both AS and DS can explain in principle the angle profiles observed at low temperatures. However as pointed out in Ref. Ocaña and Esquinazi (2002), neither AS nor DS alone can explain the magnetic field dependence of the oscillation amplitude in the whole field range $`0`$T$`B9`$T. Therefore, a more realistic picture of the thermal transport should take both into account. As argued in Ref. Vekhter and Houghton (1999), this scenario produces different regimes influenced by the predominance of either the DS or the AS in the magnetic field dependence. Thus, at constant temperature the strength of the magnetic field becomes the parameter that changes the regime. At low magnetic fields the DS dominates and the AS dominates at high fields Vekhter and Houghton (1999); Ocaña and Esquinazi (2002). In both cases, an increase of the oscillation amplitudes with increasing magnetic field strength is predicted. In Fig. 3 the oscillation amplitudes of the thermal conductivity at 3 T and 8 T for different temperatures in the twinned crystal are shown.

As the temperature increases the fourfold oscillation is no longer observable, see Figs. 3 and 4. This fact can be understood if one takes into account the thermal activation of the quasiparticles at different orientations from those of the nodes of the order parameter. In other words, if we suppose that the same processes that govern the thermal transport at low temperatures (AS and DS) are responsible for the high temperature behavior as well, we have to conclude that an increasing number of carriers appears in the direction of the heat current as the temperature is raised. As we shall see below from the numerical results using the AS model, this is, in fact, what takes place as the temperature is increased. However, although the observed angle profiles can be well described by the numerical analysis, the amplitude of the oscillations are smaller than the simulation results at $`T>30`$K. As argued for the peak of the curves in Fig. 2 and shown in the following section, this discrepancy can be solved by the inclusion of an inelastic scattering into the calculations.Hirschfeld and Putikka (1996); Bonn et al. (1993)

In Fig. 4 we show the angle dependence of the untwinned crystal at 8 T. Similar angle profiles as for the twinned sample have been found, but the amplitude of the oscillation is considerably smaller. This can be explained in terms of a larger concentration of impurity scatterers in the untwinned sample, as discussed above and in Ref. Ocaña and Esquinazi (2001). In fact, if the impurity scattering rate is sufficiently large, the relative weight of the directionality of the AS and DS mechanisms to $`\kappa _{xx}(\theta )`$ weakens. The influence of the impurity scattering has been considered in the original models by inclusion of an impurity rate.Won and Maki ; Won and Maki (2001); Hirschfeld ; Yu et al. (1995)

We note also that the symmetry of the $`\kappa _{xx}(\theta )`$curves differs slightly for the twinned and untwinned samples at the same absolute or reduced temperature. The effect of impurity scatterers can be accounted for by an impurity dependent gap parameter as done in Ref. Won and Maki . In general, however, it can be viewed as an effect related to the ratio of the modulus of the gap and the density of states of the quasiparticles at a temperature T and momentum p. Therefore, since our twinned sample has a larger critical temperature $`T_c`$ than our untwinned sample, the oxygen deficiency in the latter could be also responsible for a reduced density of states, and therefore, for a different angle patterns respect to the sample with a higher $`T_c`$.

Figure 5 shows the complete angle patterns measured in the twinned crystal at 8 T. At $`T>78`$ K no clear change in the thermal conductivity when the field is rotated parallel to the CuO2 planes is observed. This can be explained by taking into account the temperature dependence of the gap, which vanishes at $`T_c`$, and the relative increase of the inelastic scattering rate.

## IV Comparison with the two-dimensional thermal transport theory

To compare the experimental results with theory we take recently published results of the longitudinal and transverse thermal conductivityOcaña and Esquinazi (2002) into account that indicate that at high fields ($`B>2`$T), applied parallel to the planes, the main scattering mechanism for the quasiparticles appears to be the Andreev reflection by vortex supercurrents.Andreev (1964) Therefore, we use the formulation proposed by Yu et al.Yu et al. (1995) using a two dimensional model of the BRT expression for the thermal conductivity,Bardeen et al. (1959) which has been useful in the interpretation and discussion of previous results.Aubin et al. (1997); Ocaña and Esquinazi (2001, 2002) Under this model the longitudinal thermal conductivity can be written as

$$\kappa _{xx}^{el}=\frac{1}{2\pi ^2ck_BT^2\mathrm{}^2}_{p_F}^{\mathrm{}}\mathrm{d}^2p\frac{v_{gx}v_{gx}E_𝐩^2}{\mathrm{\Gamma }(𝐁,𝐩,T)}\mathrm{sech}^2\left(\frac{E_𝐩}{2k_BT}\right),$$
(1)
where $`v_{gx}`$ is the the x-axis component of the group velocity, and $`E_𝐩`$ is the quasiparticle energy. For this energy we use a free quasiparticle model $`E_𝐩^2=(𝐩^2/2m_{\mathrm{eff}}\mu )^2+\mathrm{\Delta }^2(𝐩,T)`$ where $`m_{\mathrm{eff}}`$ is the effective mass of the quasiparticle. $`\mathrm{\Gamma }(𝐁,𝐩,T)`$ may be taken as a relaxation rate given by the sum of the following scattering mechanisms acting in series: scattering of QP by impurities $`\mathrm{\Gamma }_{\mathrm{imp}}(𝐩)`$, by phonons $`\mathrm{\Gamma }_{\mathrm{ph}}(𝐩,T)`$,Houssa and Ausloss (1995, 1996a, 1996b) by quasiparticles $`\mathrm{\Gamma }_{\mathrm{qp}}(𝐁,𝐩,T)`$Hirschfeld and Putikka (1996); Quinlan et al. (1994) and AS by vortex supercurrents $`\mathrm{\Gamma }_v(𝐁,𝐩,T)`$.Yu et al. (1995) The expression for this scattering rate according to the model for the AS, proposed by Yu et al. Yu et al. (1995), is given by

$$\mathrm{\Gamma }_v(𝐁,𝐩,T)=\mathrm{\Gamma }_v^0\mathrm{exp}\left\{\frac{m^2a_\upsilon ^2\left[E_p\left|\mathrm{\Delta }(𝐩,T)\right|\right]^2}{p_F^2\mathrm{}^2\mathrm{ln}(a_\upsilon /\xi _0)\mathrm{sin}^2\psi (𝐩)}\right\},$$
(2)
where $`a_\upsilon `$ is the intervortex spacing given in this model by

$$a_\upsilon ^2=\frac{p_F^2\mathrm{}^2}{\pi \mathrm{\Delta }_0^2\gamma m_e^2}\frac{B_{c2}^{ab}}{B}.$$
(3)
We use the BCS-like parameterization for the temperature dependence of the gap amplitude

$$\mathrm{\Delta }(T)=\mathrm{\Delta }_0\mathrm{tanh}\left(2.2\sqrt{\frac{T_c}{T}1}\right),$$
(4)
and the $`d_{x^2y^2}`$-pairing symmetry in the resulting gap $`\mathrm{\Delta }(𝐩,T)`$ enters as follows

$$\mathrm{\Delta }(𝐩,T)=\mathrm{\Delta }(T)\left[\frac{\mathrm{cos}(p_xa/\mathrm{})\mathrm{cos}(p_ya/\mathrm{})}{1\mathrm{cos}(p_Fa/\mathrm{})}\right].$$
(5)
We showed recentlyOcaña and Esquinazi (2002) that at low temperatures and high fields, a quantitative agreement between this model and the oscillation amplitudes of the thermal conductivity tensor is achieved only if the intervortex spacing is increased by about five times the value defined in Eq.(3) if we use the parameters from Ref. Yu et al. (1995) ($`\mathrm{\Delta }_0=20`$ meV, Ginzburg-Landau parameter $`\kappa =100`$, anisotropy $`\gamma =4`$, $`B_{c2}^{ab}=650`$ T). Furthermore, from those fitsOcaña and Esquinazi (2002) for the twinned crystal we obtained a momentum independent impurity scattering $`\mathrm{\Gamma }_{\mathrm{imp}}^10.12`$ ps and $`(\mathrm{\Gamma }_v^0)^1(B=B_{c2})(3/2)\mathrm{\Gamma }_{\mathrm{imp}}^1`$. Assuming a total scattering rate given by the sum of the AS and impurity rates, i.e. $`\mathrm{\Gamma }(𝐁,𝐩,T)=\mathrm{\Gamma }_v(𝐁,𝐩,T)+\mathrm{\Gamma }_{\mathrm{imp}}`$, the model given by Eqs. (1) to (5) reproduces remarkably well the symmetry of the curves in Fig. 5 in the whole measured temperature range. Because in this calculation we do not take explicitly into account the inelastic scattering rate $`\mathrm{\Gamma }_{\mathrm{in}}`$ given by phonons $`\mathrm{\Gamma }_{\mathrm{ph}}`$ and by quasiparticles $`\mathrm{\Gamma }_{\mathrm{qp}}`$, the calculated oscillation amplitude above 20 K increases up to $`10`$ times that observed in the experiment.

In fig. 6 we show the results of the numerical simulation. Each of the curves in this figure has been multiplied by a factor $`f(T)`$ in order to fit the experimental oscillation amplitude (Fig. 5). We note that the symmetry of the results can be explained quite satisfactorily by this model.

This result, on one hand, supports the predominance of $`d_{x^2y^2}`$-pairing symmetry of the order parameter up to temperatures close to $`T_c`$ and, on the other hand, confirms the idea of a competition between the thermal activation of quasiparticles in the direction of the thermal current and the gap structure. We note that the $`d_{x^2y^2}`$-gap symmetry is the only gap function that can explain the whole angular patterns up to Tc. Of course, a $`s`$-gap also gives a onefold oscillation similar to the experiments above $``$ 50 K, but it does produce neither the fourfold oscillation below $``$ 15 K nor the curves between 15 K and 50 K. Thus, the set of curves in Fig. 5 can be only explained if one uses a main $`d_{x^2y^2}`$-gap into the calculations. Furthermore , the anisotropy of $`d`$-wave gap proposed for BSCCO,Ando et al. (2002) does not seem to occur in YBCO since in the latter the nodes are observed at both field orientations $`\theta `$= 45 and $`\theta `$= -45(see also Refs. Aubin et al. (1997); Yu et al. (1995); Ocaña and Esquinazi (2002)).

The factor $`f(T)`$ is related to the inelastic scattering rate $`\mathrm{\Gamma }_{\mathrm{in}}`$, which was not taken explicitly into account in the previous calculations. However, an approximation can be carried out to get roughly the temperature dependence of the total scattering rate, which could be in part associated to the temperature dependence of $`\mathrm{\Gamma }_{\mathrm{in}}`$. Since the AS mechanism is weakly temperature dependent below $`T/T_c<0.8`$, we may approximate the total scattering rate as

$$\mathrm{\Gamma }_{\mathrm{imp}}+\mathrm{\Gamma }_v+\mathrm{\Gamma }_{\mathrm{in}}(\mathrm{\Gamma }_{\mathrm{imp}}+\mathrm{\Gamma }_v)f^1(T).$$
(6)
The temperature dependence of the inelastic scattering rate is given then by the factor $`f(T)`$ in this approximation as

$$\mathrm{\Gamma }_{\mathrm{in}}f^1(T)1.$$
(7)
In Fig. 7 we plot $`f^1(T)1`$.

As expected, we recover the overall behavior of the inelastic scattering rate already described by many authorsHirschfeld and Putikka (1996); Quinlan et al. (1994); Houssa and Ausloss (1995, 1996a, 1996b). Below $`T_c`$ it decreases rapidly and becomes negligible in comparison with the impurity and AS scattering rates below $`25`$ K. Although the approximation given by (6) is too rough to obtain the true temperature dependence of the inelastic scattering rate, it is instructive to compare the obtained power dependence with results from literature. Thermal Hall angle measurements performed in the same YBCO twinned crystal show that $`\mathrm{cot}(\theta _H)=m_H/\tau _HT^4`$ down to $`20`$K ($`m_H`$ is the effective mass and $`\tau _H`$ is the Hall scattering time of the quasiparticles responsible for the Hall signal).Ocaña et al. (2001) The power law dependence obtained for $`\tau _H^1T^4`$ is similar to that obtained for the longitudinal scattering rate assuming a $`d`$wave pairing.Yu et al. (1992) Since quasiparticle-quasiparticle scattering mechanism should be the dominant temperature dependent inelastic scattering below $`T_c`$, we expect a rate proportional to the density of quasiparticles. Interestingly, within the simple two-fluid model we expect a density of quasiparticles proportional to $`(T/T_c)^4`$. On the other hand, a $`T^3`$ dependence is expected within the spin-fluctuation scattering pictureQuinlan et al. (1994).

## V Summary

In summary, we have measured the longitudinal thermal conductivity $`\kappa _{xx}`$ in two single crystals of YBCO in the presence of a planar magnetic field which was rotated parallel to the CuO2 planes, from $`T10`$ K up to a few Kelvins below $`T_c`$. Fourfold oscillations were recorded below $`20`$ K. Above this temperature the angle dependence of the longitudinal thermal conductivity changes; from the minimum at $`\theta =0^{}`$ (field parallel to the heat current) a maximum develops at high $`T`$. The observed behavior in the whole temperature range can be very well reproduced by a model involving Andreev scattering of quasiparticles by vortices and the two-dimensional BRT expression for the thermal conductivity assuming a $`d`$wave pairing. The overall results agree with the $`d_{x^2y^2}`$-pairing symmetry of the order parameter. This agreement suggests that the mechanisms that influence the behavior of the quasiparticles below $`T_c`$ and above $`10`$K are well described by the Fermi liquid theory at nearly optimal doping. The small differences found in the angle patters between the twinned and untwinned samples can be explained in terms of the different impurity concentration as well as oxygen deficiency and are accounted for by the phenomenology exposed in this work.

###### Acknowledgements.

This work was supported by the DFG under Grant DFG ES 86/4-3.


text:
# Quantum Hall fractions in rotating Bose-Einstein condensates

## Abstract

We study the Quantum Hall phases that appear in the dilute limit of rotating Bose-Einstein condensates. By exact diagonalization in a spherical geometry we obtain the ground-state and low-lying excited states of a small number of bosons as a function of the filling fraction $`\nu `$, ratio of the number of bosons to the number of vortices. We show the occurrence of the Jain principal sequence of incompressible liquids for $`\nu =1/2,2/3,3/4,4/3,5/4`$ as well as the Pfaffian state for $`\nu =1`$. The collective excitations are well described by a composite-fermion scheme.

Bose-Einstein condensates in dilute atomic gases offer a unique opportunity to investigate the physics of vortex matter when they undergo rotationMatthews99 ; Haljan01 . Indeed, recent experiments Madison00 ; Abo01  have observed the appearance of large vortex arrays at sufficient high angular velocity $`\omega `$. In addition to this phase akin to the Abrikosov lattice of type-II superconductors, there is the possibility that at larger $`\omega `$ the lattice melts Sinova02  and is replaced by a quantum Hall liquid. Consider a trap with strong confinement in the $`z`$ direction such that the system is effectively two-dimensional (2D). Then if the rotation frequency is tuned to the characteristic frequency of the harmonic confining potential in the $`xy`$ plane, the bosons feel only the Coriolis force and the system is equivalent to 2D charged bosons in a magnetic field, i.e. the conditions of the quantum Hall effect. In this regime, it has been pointed out Cooper99 ; Wilkin00  that the celebrated Laughlin wavefunction is the exact ground state for the filling fraction $`\nu =1/2`$ where $`\nu `$ is the ratio of the number of bosons to the number of vortices. Some of the excitations above this ground state are quasiparticles with fractional statistics which may eventually be probed by laser manipulations Paredes01 . Investigations by exact diagonalization have given evidence Cooper01  for even more exotic states of matterMoore91 ; Greiter92 , some involving parafermionic wavefunctions introduced in the context of the fractional quantum Hall effect for fermions ReadRezayi .

In this Letter we investigate the quantum Hall states of bosons as a function of the filling $`\nu `$ by use of exact diagonalizations in the spherical geometry Haldane83 ; Haldane85 . This allows to separate bulk from edge excitations. We show the appearance of the Bose analog of the Jain principal sequence of fractions, $`\nu =\frac{n}{n+1},\frac{p}{p1}`$. The excited states show collective modes well described by a composite fermion picture in which there is binding of one vortex per boson. We obtain evidence for the Pfaffian state Moore91 ; Greiter92  at $`\nu =1`$ by displaying its peculiar half-vortex excitations. For higher fillings, $`\nu 3/2`$, we observe some states with properties of the Read-Rezayi (RR) parafermionic states. However they show no clear tendency to convergence to the thermodynamic limit.

In the rotating frame Rosenbusch02 , the Hamiltonian describing N bosons of mass $`m`$ is given by :

$``$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{1}{2m}}(𝐩_im\omega \widehat{𝐳}\times 𝐫_i)^2+{\displaystyle \frac{1}{2}}m(\omega _0^2\omega ^2)(x_i^2+y_i^2)`$ (1)
$`+`$ $`{\displaystyle \frac{1}{2}}m\omega _z^2z_i^2+{\displaystyle \underset{i<j}{\overset{N}{}}}V(𝐫_i𝐫_j),`$
where the $`xy`$ trap frequency is $`\omega _0`$, the axial frequency is $`\omega _z`$ and the angular velocity vector is $`\omega \widehat{𝐳}`$. In the ultra-cold atomic gases the interaction takes place through s-wave scattering only and is thus given by $`V(𝐫)=(4\pi \mathrm{}^2a_s/m)\delta ^{(3)}(𝐫)`$ where $`a_s`$ is the s-wave scattering length. For $`\omega `$ close to $`\omega _0`$, the physics is that of charge-$`e`$ bosons in a magnetic field $`𝐁=(2m\omega /e)\widehat{𝐳}`$, corresponding to a magnetic length $`\mathrm{}=\sqrt{\mathrm{}/(2m\omega )}`$. There is then a 2D regime in which the boson wavefunction along the $`z`$-axis is the ground state of the harmonic oscillator and the interaction is now given by $`V^{2D}(𝐫)=g\mathrm{}^2\delta ^{(2)}(𝐫)`$ (the vector $`𝐫`$ is 2D), with $`g=\sqrt{32\pi }\mathrm{}\omega a_s/\mathrm{}_z`$, $`\mathrm{}_z=\sqrt{\mathrm{}/m\omega _z}`$ is the confinement length along $`z`$. The energy scale of the quantum Hall problem is thus set by $`g`$.

We are thus led to study the quantum Hall effect of bosons interacting through a delta potential in the lowest Landau level (LL) Sinova02-2 . To study the vortex liquids that appear as a function of the filling factor $`\nu `$, we use the spherical geometry Haldane85 ; Fano86  in which the bosons move on a sphere of radius $`R`$ in the magnetic field of a monopole $`B=\mathrm{}S/eR^2`$ at the center of the sphere, giving rise to 2S+1 cyclotron orbits in the lowest LL. In the thermodynamic limit, the filling factor $`\nu `$ is given by $`N/2S`$. However for the incompressible liquids there is in general a finite shift in the relation between the number of particles and the flux, i.e. one has generally $`2S=(1/\nu )NX`$. If we have a guess for the ground state then one can evaluate the shift and check for its validity against numerical results. For example the bosonic Laughlin state for $`\nu =1/2`$ on the sphere is realized for $`2S=2N2`$ by the wavefunction :

$$\mathrm{\Psi }_{1/2}=\underset{i<j}{}(u_iv_ju_jv_i)^2,$$
(2)
where the spinor coordinates $`(u,v)`$ are given by :

$$(u_i,v_i)=(\mathrm{cos}\theta _i/2e^{i\varphi _i/2},\mathrm{sin}\theta _i/2e^{i\varphi _i/2}).$$
(3)
This is an exact zero-energy eigenstate of the present problem Wilkin00 . We have conducted Lanczos diagonalizations for various N and flux 2S to elucidate the nature of the incompressible liquid states. States can be labelled by their total angular momentum $`L`$, contrary to the planar geometry where only the z component is conserved.

Jain sequence. The signature of incompressible states is the presence of a L=0 singlet ground state separated by a clear gap from excited states. A typical spectrum is given in fig. 1a for $`\nu =1/2`$. In the excited states we observe a well-defined collective mode which is gapped for all values of $`L`$. Clear signs of incompressible liquids are seen for the fractions $`\nu =1/2,2/3,3/4,4/3,5/4`$ : some sizes are displayed in figs. 1b-d. These fractions are the bosonic analog of the Jain sequence Jain89  of fractional quantum Hall states. They are explained by a composite particle picture in which the composite fermion (CF) is a boson bound to one vortex. Then the integer quantum Hall effect with n filled CF LLs leads to a fraction at $`\nu =\frac{n}{n+1}`$. This state is realized on the sphere for $`2S=(n+1/n)N(n+1)`$ which is exactly what we observe. The collective mode is then an exciton-like mode obtained by promoting one CF from the highest occupied LL to the next LL. On the sphere the maximal $`L`$ is then given by $`L_{max}=N/n+n1`$ in complete agreement with our results. We estimate the gap to these neutral excitations by finite-size scaling Ambrumenil89  : $`0.09`$g ($`\nu =1/2`$), 0.05g (2/3), 0.04g (3/4). The series of values for (N, 2S) is aliased Ambrumenil89  with the sequence of fractions at $`\nu =\frac{p}{p1}`$. Indeed, if (N, 2S) matches the fraction $`n/n+1`$ then it also matches fraction $`p/p1`$ for p=N/n. Hence the same data set points to the presence of the fractions $`\nu =4/3`$ and $`5/4`$ (we do not have enough points to provide gap estimates).

The two-particle correlation function $`g(r)`$ is displayed in fig. 2 for some states in the CF sequence. For $`\nu =1/2`$ it is essentially free of finite-size effects and shows the strong correlation hole characteristic of a Laughlin ground state. The state with $`\nu =2/3`$ no longer vanishes at the origin and hence has a nonzero ground state energy.

The CF sequence also include the fractions $`\nu =3/2`$ and $`\nu =2`$. For these values we find families of incompressible states in the (N, 2S) plane but they show no sign of convergence toward the thermodynamic limit, neither in the ground state energies nor in the gap values. For these fractions, we have candidates possibly originating from the Read-Rezayi parafermionic wavefunctions (see below).

For fillings less than 1/2, the spectrum has many zero-energy eigenstates that are the quasiholes of the Laughlin state eq.(2). This is a special property of the delta function interaction, the ”quasielectron” being gapped. This obscures the appearance of fractions less than 1/2. If we change the interaction from pure delta by adding a pseudopotential $`V_2`$ Haldane85  in the next allowed partial wave for bosons, i.e. $`\mathrm{}=2`$, then we find other states from the hierarchy, the strongest being $`\nu =2/5`$.

Pfaffian state. The filling $`\nu =1`$ corresponds to the absence of magnetic field acting upon the composite fermions. Previous studies Canright89  are indicative that pairing of the CFs takes place instead of a Fermi surface. An appealing wavefunction describing this phenomenon is the so-called Pfaffian state. On the sphere it can be written as :

$$\mathrm{\Psi }_{\nu =1}=\mathrm{Pf}\{\frac{1}{u_iv_ju_jv_i}\}\underset{i<j}{}(u_iv_ju_jv_i),$$
(4)
where Pf stands for the Pfaffian which is the antisymmetrized product of pair wavefunctions deGennesBook  (a fermionic version of this state is a good candidate to describe the enigmatic $`\nu =5/2`$ quantum Hall state). Calculations of overlaps between the model Pfaffian wavefunction Eq.(4)and the exact ground state suggest that it describes the physics of bosons at $`\nu =1`$ in toroidal and disk geometry Cooper99 ; Wilkin00 ; Cooper01 . The bosonic Pfaffian state is realized on the sphere for 2S=N-2 for all N even. Our calculations lead to incompressible states at these special values for N=4,6,8,10,12. There is a clear gap that extrapolates smoothly to $`0.05`$g. It appears on fig. 4a for N=12. This state has charged excitations that are different from those of a Laughlin fluid. If we add or remove one flux quantum, then two quasiparticles are created, leading to a set of low-lying states with an alternate even-odd character : see fig. 4b. This is consistent with the spectrum for two identical particles with repulsive interactions. This is observed for all accessible sizes. We consider this as a proof that the physics is different from that of the CF sequence and is the hallmark of the Pfaffian state Greiter92 .

The correlation function is shown in fig. 5. It has now a hump at the origin possibly due to the pairing of the CFs. The energy of this state is lowered by a correlation dip (instead of a hole) which appears at some characteristic radius of order one magnetic length. For larger separation, $`g(r)`$ approaches 1 but with a characteristic length scale which is definitely larger than that occurring in fig. 2. The correlation length of Pfaffian is larger than for the Jain-like fractions.

Read-Rezayi states. For larger fillings it has been suggested by Cooper et al. Cooper01  that fractions occur at $`\nu =k/2`$ and are well described by the Read-Rezayi parafermionic wavefunctions. These functions involve clustering of $`k`$ particles and are a generalization of the Pfaffian which corresponds to simple pairing, i.e. $`k=2`$. On the sphere they are realized for $`2S=\frac{2}{k}N2`$. There are possible candidates at $`\nu =3/2`$ for N=6,9,12,15, $`\nu =2`$ for N=8,12,16 and the fraction $`\nu =5/2`$ may be realized for N=15 and 20 (but the gap is very small). Contrary to the samples belonging to the CF sequence (they have different fluxes since the shifts are different between the CF sequence and the RR states) we find no sign of convergence towards a thermodynamic limit. The gap is non-monotonous as a function of the size for $`\nu =3/2`$, and for $`\nu =2`$ finite-size effects are very large, preventing any sensible extrapolation. One possibility is that these states have very large correlation lengths and are not accommodated on our largest spheres. This is consistent with the fact that the correlation function shows very strong oscillations and no hint of incompressibility.

Evidence for clustering of more particles comes from the correlation function where we see the same phenomenon as in the Pfaffian case. The hump at the origin is even more pronounced and surrounded by a correlation hole. The hump also increases with the filling albeit we cannot make quantitative statements.

If we increase the number of bosons at fixed flux, then the spectrum becomes rotor-like : the levels lie on parabolas described by effective Hamiltonian $`\frac{1}{2I}𝐋^2`$ and correlation effects disappear.

We have shown the appearance of the Jain principal sequence of quantum hall fractions in rotating Bose-Einstein condensates. The composite fermion picture gives a successful account of the observed fractions as well as their collective mode excitations. The Pfaffian state is realized at $`\nu =1`$ as seen from the special matching of flux and number of bosons as well as its half-flux quasiparticle. The gaps we estimate from our diagonalizations are all of the order of $`\mathrm{}\omega a_s/\mathrm{}_z`$.

###### Acknowledgements.

We thank Yvan Castin and Jean Dalibard for numerous discussions.


text:
# Stability of accretion disk around rotating black holes: a pseudo-general-relativistic fluid dynamical study

## 1 Introduction

For over last three decades, fluid dynamical studies of the accretion disk around black holes have been extensively studied. Shakura & Sunyaev (1973) initiated this discussion considering very simplistic but effective model of accretion disk. They chose the Newtonian gravitational potential. Novikov & Thorne (1973) and Page & Thorne (1974) studied a few aspects of the accretion disk in a full relativistic treatment. Paczyński & Wiita (1980) proposed one pseudo-potential which could approximately describe relativistic properties of the accretion disk around non-rotating black holes. Using that pseudo-potential, Abramowicz & Zurek (1981) studied some transonic properties of accretion flows. Later on, with time, number of studies of accretion disk have been carried out with modified disk structures (e.g. Abramowicz et al. 1988; Abramowicz & Kato 1989; Chakrabarti 1990, 1996a; Narayan & Yi 1994; Narayan et al. 1997, 1998; Kato et al. 1998; Mukhopadhyay 2002a) for non-rotating black holes. Also the studies of accretion disk in Kerr geometry have been explored by several authors (Sponholz & Molteni 1994; Chakrabarti 1996b, 1996c; Abramowicz et al. 1996; Peitz & Appl 1997; Gammie & Popham 1998; Popham & Gammie 1998; Miwa et al. 1998; Lu & Yuan 1998; Manmoto 2000). Those studies have been made either in a full general relativistic or in a pseudo-Newtonian approach.

The study on stability of accretion disk is one of the important criteria in this context. The possibility of steady and stable disk formation by incoming matter towards a black hole is allowed only for certain sets of initial parameters. Abramowicz & Zurek (1981) studied the effect of angular momentum on the accretion and corresponding stability of the transonic nature of the infalling matter onto the black hole. After that, Chakrabarti (1989, 1990) discussed about the formation of sonic points and shock waves in accretion disks. He showed that the formation of stable sonic points and shock in the disk are possible for a particular range of physical parameters. Recently, Mukhopadhyay & Chakrabarti (2001) have analysed the stability of accretion disk in presence of nucleosynthesis. As the temperature of the disk is very high (at least of the order of $`10^9`$K), nuclear burning can take place and generate or/and absorb energy in a large scale in the disk. Because of this generation/absorption of nuclear energy, fundamental properties of the disk, like sonic locations, temperature, Mach number, shock location (if any), etc., may be affected. Thus the stability gets influenced by the nuclear burning. If the rate of generation/absorption of nuclear energy is large, the disk structure may become unstable. In their discussion, Mukhopadhyay & Chakrabarti (2001) have also shown that the sonic locations disappear for a particular shell of accretion disk where the nuclear energy is significant compared to the mechanical energy of the accreting matter.

In all the above mentioned works, when the stability was analysed for the accretion disk, the central black hole was chosen to be non-rotating. As it is well known that no object is static in the Universe, before making any serious conclusions about the inner properties of the accretion disk, the consideration of rotation of a black hole is essential. Though a few works have been done for Kerr black holes, no such comparative study has been explored so far with different Kerr parameters. Thus automatically the following questions arise, does the rotation of black hole play an important role in the accretion properties of a disk? Does the consideration of rotation of black hole affect the disk structure for a particular parameter regime which is known to be stable for the non-rotating case? Here, we would like to address all these questions one by one.

Recently, we have proposed a new pseudo-Newtonian potential (Mukhopadhyay 2002b, hereinafter Paper-I) to describe the accretion disk around Kerr black holes. Following Paper-I, we will describe our accretion disk in a pseudo-Newtonian manner. The potential, that we are using, can describe approximately all the essential general relativistic properties of the accretion disk. It can reproduce the locations of marginally stable and bound orbit exactly or almost exactly for different Kerr parameters (see Paper-I), which are the pure general relativistic properties. Therefore, we call this study as a pseudo-general-relativistic approach. As our main interest is to study the inner region of the accretion disk, we will concentrate upon the sub-Keplerian flow, where the effect of rotation of black hole is significant. In Paper-I, we already raised several questions, which are important to address in physics of accretion disk. Here we will show, how useful is the pseudo-potential prescribed in Paper-I to study the global behaviour of the accretion disk. In full general relativity, the basic equations of accretion disk are very complicated and tedious to handle. In Early, there were no suitably good pseudo-potential for a rotating black hole. It seems that Paper-I has come up with the answer and one of our aim in this paper is to show its applicability. As the location of horizon, marginally bound and marginally stable orbits change for different Kerr parameters, any inner property of the disk is very much related to the rotation of black hole. We will follow Mukhopadhyay & Chakrabarti (2001) to analyze the stability of the disk. We will study, whether or not the sonic location(s), shock formation(s) (if any), etc., gets affected by applying rotation to the black hole. In an accretion disk, if the sonic point looses entropy or disappears by the rotational effect of black hole, we will understand that the sonic point as well as the disk is unstable because stable sonic point is necessary for the accretion onto black hole. In the next section, we will present the basic equations needed to describe the accretion disk. In §3, we will discuss about the parameter space of accretion disk and how does the rotation affect it. Subsequently, in §4, we will describe the fluid dynamical results. Finally, in §5, we will sum up all the results and make the overall conclusions.

## 2 Basic Equations

Here, throughout in our calculations, we express the radial coordinate in unit of $`GM/c^2`$, where $`M`$ is mass of the black hole, $`G`$ is the gravitational constant and $`c`$ is the speed of light. We also express the velocity in unit of speed of light and the angular momentum in unit of $`GM/c`$. The equations to be solved are given below as:

$$\frac{d}{dx}(x\mathrm{\Sigma }v)=0,$$
(1)
$$v\frac{dv}{dx}+\frac{1}{\rho }\frac{dP}{dx}\frac{\lambda ^2}{x^3}+F(x)=0$$
(2)
where, $`\mathrm{\Sigma }`$ is the vertically integrated density and $`F(x)`$ is the gravitational pseudo-Newtonian force given in Paper-I as

$$F(x)=\frac{(x^22a\sqrt{x}+a^2)^2}{x^3(\sqrt{x}(x2)+a)^2},$$
(3)
where $`a`$ indicates the specific angular momentum of the black hole (Kerr parameter). We do not consider any kind of energy dissipation in the accretion disk, as we want to check, how the black hole rotation solely can affect the disk properties. Thus the angular momentum of accreting fluid ($`\lambda `$) remains constant throughout a particular flow. Following Matsumoto et al. (1984), we can calculate the vertically integrated density as

$$\mathrm{\Sigma }=I_n\rho _eh(x),$$
(4)
where, $`\rho _e`$ is the density at equatorial plane, $`h(x)`$ is the half-thickness of the disk and $`I_n=\frac{(2^nn!)^2}{(2n+1)!}`$, where $`n`$ is the polytropic index. From the vertical equilibrium assumption, the half-thickness can be written as

$$h(x)=c_sx^{1/2}F^{1/2},$$
(5)
where $`c_s`$ is the speed of sound. We also consider the equation of state to be $`\gamma P/\rho =c_s^2`$, where $`\gamma `$ is the gas constant.

Now, combining (1) and (2) we get

$$\frac{dv}{dx}=\frac{\frac{\lambda ^2}{x^3}F\left(x\right)+\frac{c_s^2}{\gamma +1}\left(\frac{3}{x}\frac{1}{F}\frac{dF}{dx}\right)}{v\frac{2c_s^2}{\left(\gamma +1\right)v}}.$$
(6)
Far away from the black hole, $`v<c_s`$ and close to it, $`v>c_s`$, thus there is an intermediate location where the denominator of (6) must vanish. Therefore, to have a smooth solution at that location, the numerator has to be zero. This location is called the sonic point or critical point ($`x_c`$). The existence of this sonic location plays an important role in the accretion phenomena. For an accretion disk around black hole, sonic point must exist. From the global analysis of sonic point, one can understand the stability of physical parameter region which we will discuss in the next section.

As it is described above that, at $`x=x_c`$, $`\frac{dv}{dx}=\frac{0}{0}`$, using l’Hospital’s rule and after some algebra, we can get the velocity gradient of accreting matter at the sonic location as

$$\frac{dv}{dx}|_c=\frac{+\sqrt{^24𝒜𝒞}}{2𝒜},$$
(7)
where,

$`𝒜`$ $`=`$ $`1+{\displaystyle \frac{2c_{sc}^2}{(\gamma +1)v_c^2}}+{\displaystyle \frac{4c_{sc}^2(\gamma 1)}{v_c^2(\gamma +1)^2}},`$
$``$ $`=`$ $`{\displaystyle \frac{4c_{sc}^2(\gamma 1)}{v_c(\gamma +1)^2}}\left({\displaystyle \frac{3}{x_c}}{\displaystyle \frac{1}{F_c}}{\displaystyle \frac{dF}{dx}}|_c\right),`$
$`𝒞`$ $`=`$ $`{\displaystyle \frac{c_{sc}^2}{\gamma +1}}\left[\left({\displaystyle \frac{1}{F_c}}{\displaystyle \frac{dF}{dx}}|_c\right)^2{\displaystyle \frac{1}{F_c}}{\displaystyle \frac{d^2F}{dx^2}}|_c{\displaystyle \frac{3}{x_c^2}}\right]{\displaystyle \frac{3\lambda ^2}{x_c^4}}{\displaystyle \frac{dF}{dx}}|_c`$
$``$ $`{\displaystyle \frac{c_{sc}^2(\gamma 1)}{(\gamma +1)^2}}\left[{\displaystyle \frac{3}{x_c}}{\displaystyle \frac{1}{F_c}}{\displaystyle \frac{dF}{dx}}|_c\right]^2.`$
Thus, we have to integrate (6) and (7) with an appropriate boundary condition to get the fluid properties in accretion disk. From (6), we can easily find out the Mach number at the sonic point as

$$M_c=\frac{v_c}{c_{sc}}=\sqrt{\frac{2}{\gamma +1}}$$
(9)
and the corresponding sound speed as

$$c_{sc}=\sqrt{(\gamma +1)\left(\frac{\lambda ^2}{x_c^3}F_c\right)\left(\frac{1}{F_c}\frac{dF}{dx}|_c\frac{3}{x_c}\right)^1}.$$
(10)
Now integrating (1) and (2), we can write down the energy and entropy of the flow at sonic point as

$$E_c=\frac{2\gamma }{(\gamma 1)}\left(\frac{\frac{\lambda ^2}{x_c^3}F_c}{\frac{1}{F_c}\frac{dF_c}{dx}|_c\frac{3}{x_c}}\right)+V_c+\frac{\lambda ^2}{2x_c^2},$$
(11)
and

$$\dot{}_c=(\gamma K)^n\dot{M}=x_c^{3/2}F_c^{1/2}(\gamma +1)^{q/2}\left(\frac{\frac{\lambda ^2}{x_c^3}F_c}{\frac{1}{F_c}\frac{dF}{dx}|_c\frac{3}{x_c}}\right)^{\frac{\gamma }{\gamma 1}},$$
(12)
where $`K`$ is the gas constant that is basically the entropy of the system, $`V_c=(F𝑑x)|_c`$ and $`q=\frac{(\gamma +1)}{2(\gamma 1)}`$. Actually $`\dot{}`$ itself is not entropy but carries the information of it. For a non-dissipative system, $`\dot{}`$ is conserved for a particular $`\dot{M}`$, unless the shock forms there. We have to supply the sonic energy $`E_c`$ as the boundary condition for a particular flow. Then from (11), we can find out the sonic location $`x_c`$. Therefore, knowing $`x_c`$, one can easily find out the fluid velocity and sound speed at the sonic point from (9) and (10) for a particular accretion flow. These have to be supplied as further boundary conditions of the flow.

Another important issue is the formation of shock in the discussion of accretion disk. Here, whenever we mention the shock, will mean the Rankine-Hugoniot shock (Landau & Lifshitz 1959). If we generalize the conditions to form a shock in accretion disk given by Chakrabarti (1989), to the case of a rotating black hole, we get

$$\frac{1}{2}M_+^2c_{s+}^2+nc_{s+}^2=\frac{1}{2}M_{}^2c_s^2+nc_s^2,$$
(13)
$$\frac{c_{s+}^\nu }{\dot{}_+}\left(\frac{2\gamma }{3\gamma 1}+\gamma M_+^2\right)=\frac{c_s^\nu }{\dot{}_{}}\left(\frac{2\gamma }{3\gamma 1}+\gamma M_{}^2\right),$$
(14)
$$\dot{}_+>\dot{}_{},$$
(15)
where

$$\dot{}=Mc_s^{2(n+1)}\frac{x_s^{3/2}}{\sqrt{F(x_s)}}.$$
(16)
Here, subscripts ’$``$’ and ’$`+`$’ indicate the quantities just before and after the shock respectively and $`x_s`$ indicates the shock location. $`M`$ denotes the Mach number of the matter and $`\nu =\frac{3\gamma 1}{\gamma 1}`$. From (13)-(16), it is very clear that except the entropy expression, all remain unchanged with respect to the case of a non-rotating black hole. Also from (14) and (16), we get the shock invariant quantity as

$$C=\frac{\left(\frac{2}{M_+}+(3\gamma 1)M_+\right)^2}{M_+^2(\gamma 1)+2}=\frac{\left(\frac{2}{M_{}}+(3\gamma 1)M_{}\right)^2}{M_{}^2(\gamma 1)+2}$$
(17)
which remains unchanged with respect to a non-rotating case. If all the conditions (13)-(15) and (17) are simultaneously satisfied by the matter, shock will only form in an accretion disk.

## 3 Analysis of the parameter space

One of our aim is to check, how does the rotation affect the disk parameter region known for Schwarzschild black hole. Therefore, we have to analyse globally the accretion disk properties around rotating black holes. We also have to check, how does the rotation of black hole affect the sonic location, structure as well as the stability of disk. In Fig. 1a, we show the variation of disk entropy as a function of sonic location. The intersections of all the curves by the horizontal line (which is a constant entropy line) indicate the sonic points of the accretion disk for that particular entropy and rotation of the black hole. It is clearly seen that, at a particular entropy, if the rotation of black hole increases, marginally bound ($`x_b`$) and stable ($`x_s`$) orbits as well as sonic points shift to a more inner region and the possibility to have all four sonic points in the disk outside the horizon increases. As an example, for $`a=0.998`$, the inner edge of accretion disk enlarges in such a manner that the fourth sonic point in the disk appears outside the horizon. On the other hand, for retrograde orbits (counter rotating cases), $`x_b`$ and $`x_s`$ move to greater radii and all the sonic points come close to each other for a particular Kerr parameter. Similar features are reflected from Fig. 1b, where the sonic energy is plotted as a function of sonic location. Here also, the intersections of the horizontal line (which indicates a constant energy line) by all the curves indicate the sonic points of the accretion disk for that particular energy and rotation of the black hole. For both the Fig. 1a and 1b, sonic points with negative slope of the curve indicate the locations of ’saddle-type’ sonic point and positive slopes indicate the ’centre-type’ sonic point. Thus the rotation of black hole plays an important role to the formation and location of sonic points which are related to the structure of accretion disk.

In Fig. 2, with the consideration of rotation of black hole, we show how the disk region having three sonic points is affected. Figure 2a,b,c show the variation of sonic entropy as a function of sonic point, when the angular momentum of accreting matter is chosen as parameter for three different rotations of the black hole as $`a=0.5,0.5,0`$ respectively. If we join the extrema of the curves, we get the bounded parameter region of the disk where the three sonic points exist. As the angular momentum of accreting matter reduces, the possibility of forming three sonic points reduces. This is very well-understood physically; as the angular momentum of accreting matter reduces, disk tends to a Bondi-Flow like structure, that has single saddle-type sonic point. The solid curves of all three figures indicate the locus of sonic entropy with a single extremum point. Only the flows of angular momentum greater than that of the solid curve have three sonic points, otherwise for a particular flow, we have only one sonic point. In Fig. 2d, we show the projection of all the curves in Fig. 2a,b,c to the $`\dot{}_{cm}\lambda ^2`$ plane, where $`\dot{}_{cm}`$ indicates the extremum entropy at the sonic points. Figure 2d indicates the valid physical parameter region of the accretion disk that has three sonic points in the flow. It should be noted that, for a particular Kerr parameter, the accreting flow can have three sonic locations for different pairs of ($`\lambda ,\dot{}_{cm}`$). However, for a particular $`\lambda `$, if $`\dot{}_{cm}`$ decreases the stability of the sonic points as well as that of the disk diminishes too. Thus, in a section of the diagram for particular $`a`$, all points are not equally probable set of physical parameters in the formation of three sonic points. It is clear that, as the co-rotation of black hole increases, the parameter region enlarges, where three sonic points may exist. Therefore, the disk with three sonic points becomes more probable at high co-rotation. On the other hand, this region reduces for the counter-rotating flow. Similar features are noticed from Fig. 3, where the variation of sonic energies are depicted. In a similar manner, Fig. 3d indicates that the parameter space with three sonic points enlarges if the co-rotation of black hole increases, where $`E_{cm}`$ indicates the extremum energy of the different curves in Fig. 3a,b,c. The physical reasons are the following. Higher co-rotation results in a shift of $`x_b`$ and $`x_s`$ to a more inner edge of the accretion disk and an enlargement of the stable disk region (as we know that the stable circular orbit is possible in the disk upto the radius $`x_s`$). On the other hand, higher counter-rotation results in an outward shift of $`x_b`$ and $`x_s`$ and the size of the stable disk region reduces.

Figure 4 indicates the variation of energy with entropy at sonic locations for different Kerr parameters. ’I’ indicates the branch of inner ’X-type’ sonic point that is different for different Kerr parameters, while ’O’ indicates the outer ’X-type’ sonic point branch. Along the branch O, all the curves (solid, dotted and dashed) lie on top of each other. Therefore it is clear that the rotation of black hole does not affect the outer ’X-type’ sonic point. Figures 4a,b show that with the increase of Kerr parameter (co-rotation), inner sonic points of the disk shift towards the lower entropy region and the disk becomes unstable. But the outer sonic point branch remains unchanged for all values of rotation. We know that, if there is a possibility of matter in the outer sonic point branch of lower entropy to jump to the inner sonic point branch of higher entropy, the shock can form in accretion disk. As, with the increase of Kerr parameter, the inner sonic point branch shifts to the lower entropy region, the stability of both the shock and the disk is at stake. For example, in Fig. 4a, if we join the points A, B and A, D, we will find two sections of the curve namely, ABC (for solid curve) and ADE (for dashed curve) which are the parameter regions where the shock may form in accretion disk for $`a=0`$ and $`a=0.5`$ respectively. If we draw the lines of constant energy (which are horizontal in Fig. 4a), connecting the outer sonic point branch (AC for $`a=0`$ and AE for $`a=0.5`$) to the inner one (BC for $`a=0`$ and DE for $`a=0.5`$), section ABC will contain more number of lines with respect to that in section ADE (as the section ABC is bigger than ADE). On the other hand, a constant energy line indicates the possibility of shock (as it shows the possibility of increment of entropy in the flow, keeping the energy constant). Clearly, for a co-rotating black hole, the possibility of formation of shock reduces as well as the shock itself becomes unstable. The physical reason is that, with the increase of Kerr parameter, the angular momentum of the system increases, which helps the disk to maintain a high azimuthal speed of matter upto a very close proximity to the black hole, as a result the radial matter speed can overcome the corresponding sound speed at an inner radius only. On the other hand, the inner edge of the accretion disk is comparatively less stable as the entropy decreases at lower radii. As the inner sonic points form at an inner radii, with the increase of $`a`$, the disk tends to an unstable situation that decreases the value of entropy. Figure 4b shows that the possibility of shock disappears completely when $`a=0.998`$ as there is no transition possible from the outer to inner sonic branch that can increase entropy of the system. Here, the inner edge is extremely unstable in such a manner that its entropy is even lower than that of outer sonic point branch. As the inner edge of the disk spreads out for $`a=0.998`$, there are two branches of the centre-type sonic.

In Fig. 4c, the situation is opposite. As the counter-rotation of black hole increases in magnitude, the inner sonic point branch shifts towards a higher entropy region. Thus, the disk as well as the shock become more stable. Still, the outer sonic point branch is not affected due to the rotation of black hole and the branches for different Kerr parameters are superimposed. In Fig. 4c, the sections ABC (for solid curve) and ADE (for dashed curve) are for $`a=0`$ and $`a=0.5`$ respectively, show the regions where the shock may form in the disk. The physical reason to have a more stable inner sonic point branch for the counter-rotating case is as follows. Higher the counter rotation implies that the net angular momentum of the system reduces, as a result, the matter falls faster at larger radii and it becomes supersonic at a comparatively outer edge of the disk that is more stable. In a similar manner, as explained for Fig. 4a, the greater size of the section ADE than that of ABC implies the greater possibility of shock. However, for a very high counter-rotation of the black hole (see Fig. 4d), the entire system shifts to the outer side that is more stable, and the inner sonic point branch itself tends to merge to the branch of outer sonic point. Since the net angular momentum of the system becomes very small for higher counter-rotation, the corresponding centrifugal pressure onto the matter becomes insignificant, as a result the possibility of shock formation diminishes again. Therefore we can conclude that, if the angular momentum of black hole increases or decreases significantly, the shock wave in accretion disk becomes unstable and the disk itself looses stability.

## 4 Fluid Properties of Accretion Disk

We will discuss now the effect of rotation of black hole on the disk fluid. In the previous section, we have already confirmed that the rotation significantly affects the global parameter space of accretion disk. Our intention is, to incorporate the rotation of black hole and to see the change of fluid properties known for the non-rotating case.

Figure 5a shows the variation of Mach number for a particular set of physical parameter. If the black hole is chosen to be non-rotating, shock forms in the disk. When the rotation is considered keeping other parameters unchanged, the solution changes and shock disappears. If the black hole co-rotates, the inner-edge of disk becomes unstable and the matter does not find any physical path to fall steadily into the black hole. On the other hand, if the black hole counter-rotates, there is a single sonic point in the disk and the matter attains a supersonic speed only close to the black hole and falls into it. Thus, for the non-rotating and counter-rotating cases, there are smooth solutions of matter passing through the inner sonic point, while for the co-rotating cases, it is not so. The physical reasons behind these phenomena can be explained as follows. When the black hole co-rotates, the angular momentum of system increases, the radial speed of matter may not be able to overcome the centrifugal barrier to fall into the black hole for that parameter set. If the black hole counter-rotates, the angular momentum of system reduces, the centrifugal barrier smears out and matter falls steadily into the black hole. It can be mentioned here that, for the other choice of physical parameters (e.g. reducing/increasing the angular momentum of accreting matter for co/counter-rotation, changing the sonic locations etc.) the shock may appear again, even for the rotating black holes (a few such examples are depicted in Fig. 6). Figure 5b shows the corresponding density profiles in unit of $`\frac{c^6}{G^3M^2}`$ for a particular accretion rate. As the matter slows down abruptly, the density of disk fluid jumps up at the shock location for a non-rotating case. For the counter-rotating black holes, the angular momentum of the system is less, thus, to overcome the centrifugal barrier, matter does not need to attain a high radial speed away from the black hole. In this situation, the radial matter speed is less and the corresponding density of accreting fluid is high compared to a non-rotating case. We also show the density profiles for co-rotating, unstable cases. In Fig. 5c, we compare the corresponding temperature profiles in unit of $`10^9`$ ($`T_9`$). The variation of temperature is similar to that of density away from the black hole but is opposite when close to the horizon (as the density varies inversely with the temperature at a small $`x`$ but not so at a large $`x`$, see Eqn. (1)). As the higher counter-rotation Mach number lowers down, the temperature of accretion disk becomes higher.

The variation of shock invariant quantity, $`C`$ is shown in Fig. 5d for the non-rotating case. The intersections between the solution of the inner and outer sonic point branches indicate the shock locations. Since the inner shock location is unstable, the only shock forms at an outer radius, $`x=37.77`$.

Figure 6 gives examples where the shock forms in an accretion disk around rotating black hole. The variation of Mach numbers for different Kerr parameters are shown in Fig. 6a. It reflects that, for a small rotation of black hole ($`a=0.1`$), matter adjusts the sonic locations in such a manner that the outer ($`x_o`$) and inner ($`x_i`$) sonic locations shift in a more outer and more inner region respectively with respect to a non-rotating case, and the shock forms in disk at $`x=17.18`$. For the counter-rotating case ($`a=0.1`$), $`x_o`$ may remain unchanged, but $`x_i`$ has to be shifted outside to form a shock at $`x=17.25`$ (that keeps the shock location (almost) unchanged). As for the counter rotating cases, marginally stable ($`x_s`$) and marginally bound ($`x_b`$) orbits shift away from the black hole with respect to a co-rotating one, $`x_i`$ also shifts outside to stabilize the disk. If the co-rotation of black hole is higher, say $`a=0.5`$, to form a stable shock, not only the inner sonic location has to be smaller, the disk angular momentum ($`\lambda `$) should also be reduced. High rotation of a black hole results in the location of horizon ($`x_+`$) as well as $`x_b`$ and $`x_s`$ to shift inwards, and the matter attains a supersonic speed at a point in more inside of the inner edge of the disk to fall into the black hole. Thus all the phenomena, like sonic transitions, shock formations, etc., shift inside. Moreover, for the stable solution, the angular momentum of the system can not be very high, otherwise matter will be unable to overcome the centrifugal barrier. Thus, for the high rotation of a black hole, like $`a=0.5`$, $`\lambda `$ has to be reduced to form a shock (at $`x=9.51`$) and stabilize the disk. However, for $`a=0.5`$ and onwards, there is no sub-Keplerian flow for which a shock may form. For the higher counter-rotation of a black hole, the angular momentum of system reduces. Since a higher centrifugal barrier is essential to form a shock, at the higher counter-rotation, shock completely disappears in the accretion disk. Figures 6b and 6c show the corresponding variations of density and temperature respectively. In Figs. 6d,e,f, we show the corresponding variation of shock invariant quantity $`C`$ for the inner and outer sonic point branch, when $`a=0.1,0.5,0.1`$ respectively. Again it becomes clear that, for the faster rotation of a black hole, the shock location shifts inside the disk.

In Fig. 7, we show the fluid properties for $`a=\pm 0.998`$. At a very high rotation of black hole, the possibility to form a shock in accretion disk totally disappears. For the co-rotating case, the disk cannot be stable unless the angular momentum of accreting matter is very low (as we have already mentioned that the high angular momentum of the system makes the disk unstable). Therefore, there is no scope to form a shock, as to its formation, matter needs a significant centrifugal barrier as well as a high angular momentum of the accreting fluid. Here we choose, $`\lambda =1.8`$ and the matter falls almost freely. For the counter-rotating case, the matter angular momentum need not be low to stabilize the disk. However, in the sub-Keplerian accretion disk, $`\lambda `$ cannot attain a high value to increase the angular momentum as well as the centrifugal barrier of a system to form a shock. Thus, shock also disappears for a high counter-rotating case. Figure 7a shows that $`x_b`$ and $`x_s`$ shift outwards for the retrograde orbit and the matter attains a supersonic speed at an outer radius than the case of a direct orbit of the disk. However, far away from black hole, the radial velocity is less (as it needs to overcome a little centrifugal barrier) than that of co-rotating case. Figure 7b shows the variation of corresponding density. Less velocity implies a high rate of piling up of matter at a particular radius that produces a high density in the disk and thus, the density behaves in an opposite manner to Mach number. Figure 7c shows the variation of virial temperature in accretion disk. The temperature is low for the co-rotation and high for the counter-rotation. The physical reason behind it is as follows. For a co-rotation, the net speed as well as the kinetic energy of matter is less, producing less temperature. On the other hand, for a counter-rotating case, the situation is opposite. Although we show the virial temperature of the disk, in reality several cooling processes may take place that can reduce the disk temperature upto a factor of two orders.

## 5 Discussion and Conclusion

Here, we have studied the stability of accretion disk around rotating black holes in a pseudo-Newtonian approach. Following Paper-I, we have incorporated the general relativistic effects (according to Kerr geometry) to the accretion disk in such a manner that all the essential relativistic properties, like, the locations of the marginally bound and stable orbits, black hole horizon and mechanical energy of the accreting matter are reproduced in the disk for various values of Kerr parameters. Thus we can say, our method is a ’pseudo-general-relativistic’. As our interest is to check how the rotation of a black hole solely affects the disk structure, we have chosen the inviscid fluid in accretion disk. Most of the earlier studies on structure and stability analysis of the accretion disk have been done around a non-rotating black hole (Schwarzschild geometry). As no cosmic object is static, it is very important to consider the effect of rotation of the black hole, particularly for the discussion of an inner edge of accretion disk. We have found that when the rotation of a black hole is incorporated, the valid, known disk parameter region of a non-rotating black hole is dramatically affected. As a result, the location of various sonic points and shock (if any) get shifted severely and influence the disk structure. Therefore, we can conclude that in the study of an accretion disk, the rotation should be considered. Any related conclusion depends on the rotation parameter of a black hole. So, we can say that the known physical parameter regime (Chakrabarti 1989, 1990, 1996a) to form a shock, sonic locations etc. in the disk around a non-rotating black hole is not global. The mentioned disk properties and phenomena not only get modified and/or shifted in location at the disk, but also sometimes completely disappear. In our study, we have chosen the Kerr parameter from $`a=0`$ to the case of an extremely high rotation as $`a=0.998`$ and thus, we now have a clear global picture of the accretion disk.

Though, in the literature, there are a few works on the accretion disk around Kerr black holes (e.g. Chakrabarti 1996b,c; Gammie & Popham 1998; Popham & Gammie 1998), those do not concentrate on the comparative study of a rotating and non-rotating black hole and the instability of disk induced by the rotation of black hole. Chakrabarti (1996b,c) mainly concentrated on fast rotating black holes and showed solution topology in full Kerr geometry. He does not show, how the solution varies with the change of rotation of the black hole, keeping unchanged other physical parameters. Thus, it is not clear from those works, how the Kerr parameter affects the disk solution solely. On the other hand, Gammie & Popham (1998) and Popham & Gammie (1998) do not discuss the effect of rotation on the sonic points globally, which is shown here in Figs. 1-4. However, the fluid dynamical results achieved by them tally with the present work. For example, they (Gammie & Popham 1998; Popham & Gammie 1998) showed, the effect of non-zero Kerr parameter is dominant close to the black hole, that is also reflected from Fig. 4 of the present work, which shows the outer sonic point branches are unchanged, while the inner ones are shifted. They also showed, the (inner) sonic points move to smaller (larger) radii as co-rotation (counter-rotation) increases in magnitude. In our work, similar features come out from Figs. 1,2,3 and 6. Thus, we can conclude once again that, the pseudo-potential prescribed in Paper-I can reproduce the general relativistic properties in accretion disk. Unlike the works by Gammie & Popham (1998) and Popham & Gammie (1998), here we have chosen inviscid flow to understand the sole effect of rotation of the black hole on the fluid and corresponding disk structure. Finally we can conclude that the effect is significant.

Here, we have carried out the stability analysis of accretion disk by choosing the black hole to be rotating. Prasanna & Mukhopadhyay (2002) have worked on the stability of accretion disk around rotating compact objects (but not for black holes) in a perturbative approach. They incorporated the rotational effect of central object in a disk indirectly, by the inclusion of Coriolis acceleration term. Here, the compact object is black hole and its rotational effect is brought from the Kerr metric only. As we have a very self-consistent pseudo-potential (Paper-I), it has been now easy to study the accretion phenomena around a rotating black hole.

From the global analysis of sonic points, we have established that the disk becomes unstable for a higher co-rotation at a particular angular momentum of the accreting matter. On the other hand, to form a shock, accretion disk must have a stable inner sonic point. Thus for higher $`a`$ (Kerr parameter), shock is unstable, as the inner region of disk itself is not stable. In that regime, shock may disappear by any disturbance created on the accreting matter and thus the matter may not get a steady inner sonic point to form a stable inner disk. For the counter-rotating cases, the angular momentum of the system reduces and the matter falls to the black hole more rapidly. As there is no significant centrifugal barrier to slow down the matter, the possibility of shock reduces again. For a very high counter-rotation, shock disappears completely. Thus we can conclude that the parameter region where the shock is expected to form (Chakrabarti 1989, 1990, 1996a) for non-rotating black holes, gets affected for the rotating ones. For a positive $`a`$, shock might have been formed for a lesser value of the angular momentum of accreting matter. On the other hand, for a negative $`a`$, to form a shock, the angular momentum of accreting matter should have a larger value than that of non-rotating cases. As we are studying the sub-Keplerian accretion flow, the angular momentum of matter cannot be unlimitedly high (e.g., for a non-rotating black hole, $`\lambda \stackrel{<}{}3.6742`$ at last stable orbit). Thus, to form a shock, matter cannot attain the angular momentum beyond a certain limit. Also for a highly co-rotating black hole, to stabilize the disk, $`\lambda `$ has to be reduced by some physical process, otherwise the radial speed could not overcome the centrifugal pressure to form a steady inner edge of the disk. Therefore, the disk would not be transonic and could be unstable in this situation. Because of a reduced $`\lambda `$, centrifugal force as well as shock strength reduces, as a result, the possibility of a shock is also reduced for both the co-rotating and counter-rotating black holes. However, at an intermediate co-rotation, the parameter region of accretion disk enlarges having three sonic points and there is a possibility of a stable shock formation. On the other hand, for a counter-rotation, the valid parameter region with three sonic points at the disk shortens and thus the possibility of shock reduces.

Now the question may arise, how would a flow behave in the region where sonic points disappear because of the higher rotation of black hole? It can be thought as, matter would try to pass through same sonic point as of non-rotating (or slowly rotating) case. Because of the absence of sonic point, the searching procedure for the true sonic location should give rise to an unstable solution. Only the time-dependent simulation of accretion disk can tell whether this idea is correct or not. Overall we can say that the rotation of a black hole affects the sonic locations of an accretion disk which are directly related to the formation of stable disk structure. Therefore, the stability of an accretion disk is strongly related to the rotation of central black hole.

It has also been confirmed that the potential given in Paper-I is very well applicable for the global solution and corresponding fluid dynamical study of the accretion disk around a rotating black hole. The set of a full general relativistic disk equations is very complicated for a Kerr black hole, particularly, when the accreting fluid is chosen to be viscous. Thus, in the study of a relativistic accretion disk, this pseudo-potential is very effective to avoid the complexities of full general relativistic fluid equations. The next step should be the study of a viscous accretion disk using this general pseudo-potential.

I am grateful to my friend Subharthi Ray (of IF-UFF) for careful reading the manuscript and amending the grammatical errors.