text:
# Analytical estimate of effective charges at saturation in Poisson-Boltzmann cell models

## I Introduction

In colloidal suspensions, the charge ratio between the poly-ions (colloids) and the counter-ions or salt ions may be larger than 104. This considerable asymmetry invalidates most of standard liquid state theories. The concept of charge renormalization however allows to significantly simplify the description of such systems: the electrostatic coupling between the poly-ions and oppositely charged micro-species induces a strong accumulation (or “condensation”) near the poly-ion surface. The idea is then to consider a colloid plus its “captive” micro-ions as an entity carrying an effective (or renormalized) charge $`Z_{\text{eff}}`$, a priori much smaller than the bare charge $`Z_{\text{bare}}`$. Consequently, except in the immediate vicinity of the colloids where linearization schemes fail, the interactions in the suspension are well described by Debye-Hückel-like linearized theories, provided that the bare charge is replaced by the renormalized one $`Z_{\text{eff}}`$. Schematically, the linear effects of screening by an electrolyte induce a dressing of the bare Coulomb potential \[$`Z_{\text{bare}}/rZ_{\text{bare}}\mathrm{exp}(\kappa r)/r`$ where $`\kappa ^1`$ is the Debye length, see below\], while the non-linear effects of screening imply the identification $`Z_{\text{bare}}\mathrm{exp}(\kappa r)/rZ_{\text{eff}}\mathrm{exp}(\kappa ^{}r)/r`$ (with possibly $`\kappa ^{}\kappa `$ ).

Several reviews discussing the notion of charge renormalization have appeared recently . In the colloid science field this concept has been introduced by Alexander et al.  in the context of the Poisson-Boltzmann (PB) cell model, but had been widely accepted since the fifties in the field of linear polyelectrolytes . The definition of an effective charge from the far field potential created by an isolated macro-ion in an electrolyte is unambiguous , at least for the simple cases of spherical or infinitely long cylindrical macro-ions we shall consider here, even if the case of less symmetric poly-ions deserves more attention (see e.g. ). Within a cell model, introduced to replace the complicated many body problem of colloids in solution by a simpler one-particle system , the definition of an effective charge is more elusive, but the celebrated proposal made by Alexander et al. amounts to find the optimal linearized PB potential matching the non-linear one at the cell boundary. The cell approach, which validity has been assessed by several studies , appears to provide a reasonable description of solutions containing mono-valent counter-ions . The effective charge is obtained by integrating the charge density deduced from the linearized potential over the region accessible to micro-ions, or equivalently from Gauss theorem at the colloid’s surface. Note that a renormalization of $`Z`$ also implies a renormalization of the screening constant $`\kappa `$. For low poly-ion bare charges $`Z_{\text{bare}}`$, linearizing PB equation is a valid approximation so that $`Z_{\text{eff}}Z_{\text{bare}}`$, whereas in the opposite limit of high bare charges, $`Z_{\text{eff}}`$ was found numerically to saturate to a value independent of $`Z_{\text{bare}}`$: $`Z_{\text{eff}}Z_{\text{sat}}`$. The saturation value $`Z_{\text{sat}}`$ depends on the geometry of the colloid, the temperature or the quantity of added-salt. Unfortunately, no analytical prediction is available for these dependences.

In this letter, we propose an approximate analytical expression for the effective charge of spherical or rod-like macro-ions, as a function of macro-ion density or salt content. Our approach starts with an estimate of the saturation value of the effective charge in the infinite dilution limit, as deduced from recent analytical results . In the case of finite colloidal dilution, the osmotic equilibrium of the suspension with a salt reservoir is modeled using a Donnan equilibrium approximation. This allows to derive a simple polynomial equation of degree 4 fulfilled by the saturation value $`Z_{\text{sat}}`$ of the effective charge $`Z_{\text{eff}}`$, in the case of a symmetric 1:1 electrolyte. In the limiting case of no added salt, the above equation is easily solved analytically and provides results in quantitative agreement with the saturation effective charges following Alexander’s prescription.

## II General framework and method

We consider first the situation of an isolated macro-ion of given surface charge density in a electrolyte of bulk density $`n_0`$ (no confinement). The solvent is considered as a medium of uniform dielectric (CGS) permittivity $`\epsilon `$. Within Poisson-Boltzmann theory, the electrostatic potential, when assumed to vanish far from macro-ion obeys the equation

$$^2\varphi =\kappa _0^2\mathrm{sinh}\varphi ,$$
(1)
where the screening factor $`\kappa _0`$ is defined as $`\kappa _0^2=8\pi \mathrm{}_Bn_0`$ and the Bjerrum length quantifies the strength of electrostatic coupling: $`\mathrm{}_B=e^2/(\epsilon kT)`$ ($`e>0`$ denotes the elementary charge and $`kT`$ is the thermal energy). A complete asymptotic solution of Eq. (1) has been obtained recently by Shkel et al.  in spherical and cylindrical geometry. From the far-field behaviour of the corresponding solutions, one obtains after some algebra the effective charges at saturation

$`Z_{\text{sat}}{\displaystyle \frac{\mathrm{}_B}{a}}=4\kappa _0a+6+𝒪\left({\displaystyle \frac{1}{\kappa _0a}}\right)\text{for spheres}`$ (2)
$`\lambda _{\text{sat}}\mathrm{}_B=2\kappa _0a+{\displaystyle \frac{3}{2}}+𝒪\left({\displaystyle \frac{1}{\kappa _0a}}\right)\text{for cylinders}.`$ (3)
In these equations, $`a`$ denotes the radius of the macro-ion under consideration, $`Ze`$ the total charge in the case of spheres and $`\lambda e`$ the line charge density in the case of cylinders (with infinite length). Expressions (2) and (3) are the exact expansions of the saturation charges in the limit of large $`\kappa _0a`$ but are in practice accurate as soon as $`\kappa _0a>1`$ (not shown).

Our goal is to translate these relations into expressions that would approximate the effective charges in confined geometry, where the macro-ion is enclosed in a cell . To this aim, we first find an approximation for the relevant screening factor $`\kappa _{}`$ before inserting it into (2) and (3) making the substitution $`\kappa _{}\kappa _0`$. The implicit assumption is that the mean salt density in the cell is related to the effective charge in a similar manner as in the infinite dilution limit.

In confined geometry, PB equation still takes the form (1), where $`\kappa _0=(8\pi \mathrm{}_Bn_0)^{1/2}`$ should now be considered as the inverse screening length in a (neutral) salt reservoir in osmotic equilibrium with the solution, through a membrane permeable to micro-species ($`n_0`$ is thus now the salt density in the reservoir). It may be shown that the relevant screening factor in the cell is related to the micro-ions density at the cell boundary :

$$\kappa _{}^2=4\pi \mathrm{}_B\left[\rho ^+(R_{WS})+\rho ^{}(R_{WS})\right].$$
(4)
In this equation $`R_{WS}`$ is the radius of the Wigner-Seitz (WS) confining cell. We have recently proposed an efficient prescription to compute $`\kappa _{}`$ without solving the complicated non-linear problem . It is however impossible to deduce an analytical expression of $`\kappa _{}`$ from this approach and we resort to the following approximation. We assume that the micro-ion densities are slowly varying in the WS cell so that the mean densities $`n^\pm `$ provide a reasonable estimation of the boundary densities $`\rho ^\pm (R_{WS})`$. We thus write

$$\kappa _{}^2=4\pi \mathrm{}_B(n^++n^{})=4\pi \mathrm{}_B(2n^++Z_{\text{sat}}\rho ),$$
(5)
where the last equality follows from the electroneutrality constraint ($`\rho `$ is the density of colloids, assumed positively charged without loss of generality). Note that at the level of a linear theory, the effective charge $`Z_{\text{sat}}`$ and not the bare one $`Z_{\text{bare}}`$ enters this expression.

We now need to relate the mean densities $`n^\pm `$ to $`n_0`$, the concentration in the reservoir (the so called Donnan effect, see e.g. ). Chemical equilibrium imposes $`\rho ^+\rho ^{}=n_0^2`$ at any point in the cell. We again assume this relation to hold for the mean densities, so that $`n^+n^{}=n_0^2`$. This leads to

$$\kappa _{}^4=\kappa _0^4+(4\pi \mathrm{}_BZ_{\text{sat}}\rho )^2.$$
(6)
Up to now the reasoning is quite general and independent of the geometry. For spheres with radius $`a`$ and packing fraction $`\eta =4\pi \rho a^3/3=(a/R_{WS})^3`$, and for rods with packing fraction $`\eta =(a/R_{WS})^2`$, we obtain

$$(\kappa _{}a)^4=(\kappa _0a)^4+(3\eta Z_{\text{sat}}\mathrm{}_B/a)^2\text{(spheres)}\text{and}(\kappa _{}a)^4=(\kappa _0a)^4+(4\eta \lambda _{\text{sat}}\mathrm{}_B)^2\text{(cylinders)}$$
(7)
Supplementing these equations with the $`\kappa `$ dependence of $`Z_{\text{sat}}`$ obtained in Eq. (2) for spheres and (3) for cylinders, leads to an equation of degree 4 satisfied by $`\kappa _{}`$ (or equivalently $`Z_{\text{sat}}`$).

This equation can be solved analytically. We give the general solution in appendix A. However, without added salt (formally $`\kappa _0=0`$), the solutions of these equations take a particularly simple form:

$`Z_{\text{sat}}{\displaystyle \frac{\mathrm{}_B}{a}}=6+24\eta +12\sqrt{2\eta +4\eta ^2}`$ (8)
$`\lambda _{\text{sat}}\mathrm{}_B={\displaystyle \frac{3}{2}}+8\eta +2\sqrt{6\eta +16\eta ^2}`$ (9)
We compare on Figs. (1) and (2), the results obtained following this route with the full PB estimate following Alexander’s prescription. The latter involves a numerical resolution of the PB equation. Both the no-salt and the finite ionic strength case are considered. The agreement is seen to be quite good in view of the minimum number of ingredients involved in the present approach.

## III Conclusion

In this contribution, we propose a simple approximation scheme to compute the saturation value of the effective charge in concentrated suspensions of highly charged colloidal particles. Our estimate starts with an asymptotic expression of the effective charge in the infinite dilution limit, obtained from recent analytical results . The case of finite colloidal dilution is described using a Donnan equilibrium approximation for the osmotic equilibrium of the suspension with a salt reservoir. This calculation leads to an analytical estimate of $`Z_{\text{sat}}`$, the saturation value of the effective charge, as a function of the density of colloids and salt concentration. Our starting points \[Eqs. (2) and (3)\] neglect contributions of order $`(\kappa _0a)^1`$ and become inaccurate for $`\kappa a<1`$. This implies that our effective charges at finite density becomes less reliable in the salt free case for small packing fractions $`\eta `$. However, this is seen to occur for very small $`\eta `$ only (see e.g. Fig. 1).

These results are obtained at the level of the mean-field Poisson-Boltzmann theory. We note however that the existence of a saturation value of the effective charge, $`Z_{\text{sat}}`$, independent of the bare charge, is indeed confirmed in more refined approaches in the colloidal limit $`a\mathrm{}_B`$, with $`a`$ the colloid size and $`\mathrm{}_B`$ the Bjerrum length (see e.g. Groot  using the primitive model). These results show moreover that PB theory becomes successful in the afore-mentioned colloidal limit, $`a\mathrm{}_B`$ . Eventually the saturation picture within PB theory is in quantitative agreement with experimental data for the osmotic pressure .

It is eventually instructive to reconsider the results reported by Alexander et al. in their original paper . At a packing fraction $`\eta =0.125`$, they find numerically a saturation value for $`Z_{\text{eff}}`$ of the order of 15, in $`a/\mathrm{}_B`$ units (spherical colloid). On the other hand our expression (8) for the same $`\eta `$ gives $`Z_{\text{sat}}\mathrm{}_B/a=9+3\sqrt{5}15.7`$, i.e. very close to the value found in . Note that this number of $`15`$ has subsequently been often quoted in the literature as a “standard” value of the effective charge. The saturation value of the effective charge however crucially depends on the colloid volume fraction, and the full density dependence of $`Z_{\text{sat}}`$ should be taken into account in finite concentration cases.

Acknowledgments: We would like to thank M. Deserno and H.H. von Grünberg for useful discussions.

## A Analytical expression for the effective charge

We give here the explicit expression for the effective charge as a function of the colloid packing fraction $`\eta `$ and salt concentration in the reservoir, $`n_0=\kappa _0^2/8\pi \mathrm{}_B`$.

Inserting Eqs. (2) or (3) into Eq. (7), one gets the quartic equation for $`X=\kappa _{}a`$ :

$$X^4=X_0^4+\eta ^2(\alpha +\beta X)^2$$
(A1)
with $`X_0=\kappa _0a`$. Values for $`\alpha `$,$`\beta `$ are respectively $`\alpha =12`$ and $`\beta =18`$ for spheres, and $`\alpha =8`$ and $`\beta =6`$ for cylinders.

Let us introduce the following quantities :

$`\mathrm{\Lambda }=`$ $`12X_0^412\beta ^2\eta ^2+\alpha ^4\eta ^4`$ (A2)
$`\mathrm{\Gamma }=`$ $`108\alpha ^2\beta ^2\eta ^42\alpha ^6\eta ^672\alpha ^2\eta ^2(X_0^4+\beta ^2\eta ^2)`$ (A3)
$`\mathrm{\Phi }=`$ $`72X_0^4\alpha ^2\eta ^2+36\alpha ^2\beta ^2\eta ^42\alpha ^6\eta ^6`$ (A4)
$`\mathrm{\Psi }=`$ $`{\displaystyle \frac{1}{32^{1/3}}}\left(\mathrm{\Phi }+\sqrt{\mathrm{\Gamma }^24\mathrm{\Lambda }^3}\right)^{1/3}`$ (A5)
$`\mathrm{\Theta }=`$ $`{\displaystyle \frac{\mathrm{\Lambda }}{9\mathrm{\Psi }}}+\mathrm{\Psi }`$ (A6)
The solution $`X=\kappa _{}a`$ is then found as

$`\kappa _{}a=`$ $`{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{2\alpha ^2\eta ^2}{3}}+\mathrm{\Theta }}+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{4\alpha ^2\eta ^2}{3}}\mathrm{\Theta }+{\displaystyle \frac{4\alpha \beta \eta ^2}{\sqrt{\frac{2\alpha ^2\eta ^2}{3}+\mathrm{\Theta }}}}\right)^{1/2}`$ (A7)
The effective charge is then obtained by replacing this value for $`\kappa _{}`$ in Eq. (2) for spherical macro ions, and in Eq. (3) for rod-like macro ions.

In the no-added salt case, these expressions reduce to Eqs. (8) and (9), as indicated in the text.


format: text

text:
# The oscillation and stability of differentially rotating spherical shells: The normal-mode problem

## 1 Background and motivation

This paper describes a dynamical problem that is relevant to many areas of astrophysics. We consider the oscillations of a simple system exhibiting differential rotation, and discuss issues concerning the role of corotation points and the emergence of dynamical instabilities. This is known to to be a technically challenging problem (analogous to that of Couette flow in standard fluid dynamics) that remains to be understood in detail. In addition to being mathematically interesting, the problem is of potentially great astrophysical relevance since the oscillations of compact objects such as neutron stars or strange stars are a promising source of detectable gravitational waves. Of particular interest are oscillations that can grow via either a dynamical instability or a secular instability to the emission of gravitational waves, such as the Chandrasekhar-Friedman Schutz (CFS) instability (Friedman & Schutz, 1978a, b). The most promising instability scenarios concern hot newly born neutron stars, which one would expect to possess some differential rotation. This expectation is borne out by the latest studies of rotational core collapse, which indicate that the remnant will indeed emerge from the collapse rotating differentially (Zwerger & Müller, 1997; Dimmelmeier, Font & Müller, 2002). Subsequent accretion of supernova remnant material (Watts & Andersson, 2002), or material from a companion (Fujimoto, 1993), may also drive differential rotation. Differential rotation may also be maintained by the oscillations themselves. Recent work on r-modes in neutron stars suggests that the modes can drive a uniformly rotating star into differential rotation via non-linear effects (Rezzolla et al, 2000; Levin & Ushomirsky, 2001; Lindblom, Tohline & Vallisneri, 2001). A differentially rotating massive neutron star may also be generated following a binary neutron star merger. The maximum mass of a differentially rotating merger remnant can exceed significantly the maximum mass of a uniformly rotating remnant, remaining stable over many dynamical timescales (Shibata & Uryu, 2000; Baumgarte, Shapiro & Shibata, 2000).

It is clear that differential rotation is an important factor in modelling these systems. Nevertheless, the vast majority of secular instability studies to date have focused on uniformly rotating neutron stars (exceptions include studies by Imamura et al (1995), Yoshida et al (2002), and Karino & Eriguchi (2002)). The situation is somewhat different for the dynamical bar mode instability (see for example Shibata, Baumgarte & Shapiro (2000) and New & Shapiro (2001)). There have been a number of large scale numerical simulations aimed at understanding the bar mode instability. These studies provide significant insights into the dynamical stability of differentially rotating systems. Yet it is clear that a better theoretical understanding of this problem is desirable. This need is well illustrated by the recent discovery of dynamical instabilities that become active at very low values of $`T/|W|`$ (the ratio of rotational kinetic energy to gravitational potential energy) (Centrella et al, 2001; Shibata, Karino & Eriguchi, 2002).

There are two issues associated with differential rotation that make this problem particularly interesting. Firstly, differential rotation may introduce new dynamical instabilities (Balbinski, 1985; Luyten, 1990a). Secondly, the dynamical equations are complicated by the presence of corotation points. These are points where the pattern speed of a mode matches the local angular velocity, and at these points the governing equations are formally singular. This gives rise to a continuous spectrum and could in principle permit other corotating solutions. These two important issues have however attracted little attention in the literature. In particular, it is not at all clear what the dynamical role of the continuous spectrum may be.

Karino, Yoshida & Eriguchi (2001) studied r-mode oscillations in rapidly and differentially rotating Newtonian polytropes for two particular rotation laws. They found that the r-mode frequencies shifted towards the corotation region as the amount of differential rotation increased and at high differential rotation the codes failed, apparently due to the presence of corotation points. Solutions with corotation points were not, however, investigated. Two subsequent studies of a differentially rotating shell for the same two rotation laws also focused on modes without corotation points (Rezzolla & Yoshida, 2001; Abramowicz, Rezzolla & Yoshida, 2002). The authors claimed that for these rotation laws the lowest frequency r-modes did not develop corotation points, even at high differential rotation.

For papers that deal more fully with dynamical instabilities and corotation points we must go back a number of years (or consider the related problem of differentially rotating disks for which there is a vast literature). In the 1980’s Balbinski carried out a number of studies on $`z`$-independent modes of differentially rotating cylinders (for cylindrical coordinates $`\varpi ,\varphi ,z`$) that addressed both issues (Balbinski, 1984a, b, 1985). The first paper reported an analytic solution for the continuous spectrum for one particular rotation law. Balbinski then solved the initial value problem for simple initial data, and found that the physical perturbation associated with the continuous spectrum died away on a dynamical timescale. The second paper examined the secular stability of the continuous spectrum to the emission of gravitational waves. The final paper investigated the behaviour of modes that had no corotation points at low differential rotation. As the degree of differential rotation increased, some of the modes were seen to enter corotation and in due course merge to give rise to a dynamical instability. In the early 1990’s, Luyten generalised the work in Balbinski (1985) to other rotation laws (Luyten, 1990a). In two subsequent papers he solved for the same $`z`$-independent modes in self-gravitating homogeneous masses rather than on cylinders (Luyten, 1990b, 1991). Although Luyten’s method could not track modes with corotation points, modes appeared to enter and leave the corotation region as the degree of differential rotation increased. Luyten too reported finding an instability at high differential rotation.

There is clearly a need to analyse modes with a $`z`$-dependence such as the toroidal r-modes. As a first step towards this goal, we present a detailed analysis of the modes of a differentially rotating spherical shell. We hope to improve our understanding of the technical challenge posed by the presence of corotation points, and develop techniques that may later be generalised to the full problem of differentially rotating stars. By investigating dynamical instabilities we may also be able to rule out certain rotation laws, or set a threshold for the maximum degree of differential rotation that could be sustained. We select four rotation laws from the literature in order to compare our results to previous work, and one new law. We note however that in all but one case these laws were selected for mathematical convenience rather than physical reasons. One would clearly prefer to use a rotation law for a young neutron star that is grounded in physics, emerging for example from simulations of core collapse (Zwerger & Müller, 1997; Dimmelmeier, Font & Müller, 2002). We hope to do this in future work.

## 2 Formulating the problem

We consider an incompressible fluid confined to move in a thin spherical shell. We take the shell radius to be $`R`$. Its thickness is taken to be so small that we do not have to account for dynamics in the radial direction. On these assumptions, and taking the perturbation to depend on exp($`im\phi `$), the $`\theta `$-component of the Euler equation simplifies to

$$_t\delta v_\theta +im\mathrm{\Omega }\delta v_\theta 2\mathrm{\Omega }\mathrm{cos}\theta \delta v_\phi =\frac{1}{\rho R}_\theta \delta P,$$
(1)
and the $`\phi `$-component to

$$_t\delta v_\phi +im\mathrm{\Omega }\delta v_\phi +\stackrel{~}{\mathrm{\Omega }}\delta v_\theta =\frac{im}{\rho R\mathrm{sin}\theta }\delta P,$$
(2)
where the equilibrium vorticity is

$$\stackrel{~}{\mathrm{\Omega }}=2\mathrm{\Omega }\mathrm{cos}\theta +\mathrm{sin}\theta _\theta \mathrm{\Omega }.$$
(3)
The continuity equation for an incompressible fluid is

$$_\theta (\mathrm{sin}\theta \delta v_\theta )+im\delta v_\phi =0.$$
(4)
This tells us that only toroidal modes are permitted in this problem. Defining the stream function (in the standard way) as

$$U=\delta v_\theta \mathrm{sin}\theta ,$$
(5)
we can combine equations (1) and (2) to give the vorticity equation

$$(i_tm\mathrm{\Omega })_\theta ^2U+\frac{m_\theta \stackrel{~}{\mathrm{\Omega }}}{\mathrm{sin}\theta }U=0,$$
(6)
where

$$_\theta ^2=_\theta ^2+\frac{\mathrm{cos}\theta }{\mathrm{sin}\theta }_\theta \frac{m^2}{\mathrm{sin}^2\theta }$$
(7)
is the Laplacian on the unit sphere. Changing variables to $`x=\mathrm{cos}\theta `$ the vorticity equation becomes,

$$(i_tm\mathrm{\Omega })\left[(1x^2)_x^22x_x\frac{m^2}{1x^2}\right]Um_x\stackrel{~}{\mathrm{\Omega }}U=0.$$
(8)
In terms of $`x`$ the equilibrium vorticity is

$$\stackrel{~}{\mathrm{\Omega }}=2x\mathrm{\Omega }(1x^2)_x\mathrm{\Omega }.$$
(9)
Equation (8) is singular at the poles ($`x=\pm 1`$). Performing a Frobenius expansion around these points,

$$U=(1\pm x)^\alpha \underset{k=0}{\overset{\mathrm{}}{}}a_k(1\pm x)^k$$
(10)
we find that the regular solution has $`\alpha =m/2`$. Factoring out this behaviour and defining a new variable $`u`$ by

$$U=u(1x^2)^{m/2},$$
(11)
we find that $`u`$ obeys the equation

$`(m\mathrm{\Omega }i_t)\left[(1x^2)_x^22(m+1)x_xm(m+1)\right]u`$
$`+m_x\stackrel{~}{\mathrm{\Omega }}u=0.`$ (12)
We will use both equations (8) and (2) in this paper.

## 3 Uniform rotation: The r-modes

Before we proceed to discuss differential rotation it is useful to recall the standard result for uniform rotation. After all, we are interested in understanding how the mode spectrum changes as differential rotation is introduced. The oscillations of a uniformly rotating incompressible shell have been studied by a number of authors (Haurwitz, 1940; Stewartson & Rickard, 1969; Levin & Ushomirsky, 2001), and are sometimes referred to as Rossby-Haurwitz waves. Rossby-Haurwitz waves are however the limit on the shell of the r-modes of stellar perturbation theory, and we will use the latter term. For uniform rotation we have

$$_\theta \stackrel{~}{\mathrm{\Omega }}=2\mathrm{sin}\theta \mathrm{\Omega }.$$
(13)
Assuming a mode solution with time dependence $``$ exp($`i\sigma t`$) equation (6) becomes

$$(\sigma m\mathrm{\Omega })_\theta ^2U2m\mathrm{\Omega }U=0.$$
(14)
Assuming that $`\sigma m\mathrm{\Omega }`$ we have

$$_\theta ^2U=\frac{2m\mathrm{\Omega }}{\sigma m\mathrm{\Omega }}U.$$
(15)
Expanding $`U`$ in spherical harmonics

$$U=\underset{l}{}C_l^mY_l^m,$$
(16)
then by Legendre’s equation we must have

$$_\theta ^2U=\underset{l}{}l(l+1)C_l^mY_l^m.$$
(17)
Thus we must have

$$\underset{l}{}l(l+1)C_l^mY_l^m=\frac{2m\mathrm{\Omega }}{\sigma m\mathrm{\Omega }}\underset{l}{}C_l^mY_l^m,$$
(18)
from which we deduce that the mode solutions are such that (i) only a single $`C_l^m`$ is non-vanishing and (ii) the corresponding frequency is

$$\sigma =m\mathrm{\Omega }\frac{2m\mathrm{\Omega }}{l(l+1)}.$$
(19)
These are the familiar r-modes.

There is, however, another possible solution, the geostrophic solution. Consider what happens if we set $`\sigma =m\mathrm{\Omega }`$. It is clear from equation (14) that the only possible solution in this case is $`U=0`$, which corresponds to a zero velocity perturbation. Consider however what this means for the displacement $`\stackrel{}{\xi }`$. In the case of uniform rotation,

$$\delta v_\theta =i(m\mathrm{\Omega }\sigma )\xi _\theta ,$$
(20)
and

$$\delta v_\phi =i(m\mathrm{\Omega }\sigma )\xi _\phi .$$
(21)
It is clear that, if $`\sigma =m\mathrm{\Omega }`$, a vanishing velocity perturbation is compatible with a non-zero displacement. This solution, trivial in the case of uniform rotation, must be considered when the rotation is differential as it is this solution that gives rise to the continuous spectrum. Note the similarity with the standard r-modes, which are a trivial solution of the perturbation equations for a non-rotating star.

## 4 The normal mode problem

Let us now consider the case of a differentially rotating shell. Assuming a mode solution with time dependence $``$ exp($`i\sigma t`$) equation (2) becomes

$`(m\mathrm{\Omega }\sigma )\left[(1x^2)u^{\prime \prime }2(m+1)xu^{}m(m+1)u\right]`$
$`+m\stackrel{~}{\mathrm{\Omega }}^{}u=0,`$ (22)
where the primes indicate derivatives with respect to $`x`$. At the poles, we require that the solution be regular. In addition, symmetry allows us to require that $`u`$ be either even or odd under reflection through the equator. This implies that we can set boundary conditions on the equator that either $`u`$ or $`u^{}`$ be zero there. These boundary conditions allow us to solve the equation on the range $`0x1`$ rather than on the full range $`1x1`$.

If $`\sigma m\mathrm{\Omega }(x)`$ for all points on the shell we have a regular eigenvalue problem. If, however, $`\sigma =m\mathrm{\Omega }(x)`$ at a point $`x=x_c`$, then $`x_c`$ is a corotation point and equation (4) is formally singular. Corotation points are points where the pattern speed of a perturbation is equal to the local angular velocity. Note that the term corotation point is a little misleading as $`x=x_c`$ really represents a line of latitude on the shell. Equation (4) has three classes of solution:

* Discrete real-frequency modes outside corotation, i.e. for which $`\sigma m\mathrm{\Omega }(x)`$ for $`|x|1`$.
* Discrete complex frequency modes for which $`\sigma =\sigma _r+i\sigma _i`$. Dynamically unstable (growing) modes correspond to $`\sigma _i>0`$.
* Corotating solutions for which $`\sigma =m\mathrm{\Omega }(x)`$ at a corotation point $`x_c`$, where $`|x_c|1`$.

The first two kinds of solutions are easy to determine numerically for any given rotation law. In contrast, the solutions with corotation points present a real challenge. Nonetheless, these solutions are of great interest, not least because of their association with the continuous spectrum. For this reason much of the following discussion will concern solutions of the third kind. In the following sections we shall refer to solutions with corotation points as being “in corotation”. Modes without corotation points that subsequently develop them as the amount of differential rotation is increased will be referred to as “crossing the corotation boundary”. It should be noted, however, that solutions in corotation are in fact only corotating at one particular latitude.

## 5 A necessary stability criterion

Before discussing detailed solutions for particular rotation laws it is useful to check whether unstable mode solutions are likely to exist. Such solutions would typically correspond to dynamical instabilities, and hence they are of great physical significance.

For inviscid parallel shear flow it is possible to derive a necessary condition for stability that depends on the form of the rotation law, the Rayleigh inflexion point criterion (see for example Drazin & Reid (1982)). Here we derive a similar condition for modes of a differentially rotating shell. We start with equation (8), assuming a mode solution with complex $`\sigma `$,

$$(m\mathrm{\Omega }\sigma )\left[\frac{d}{dx}[(1x^2)U^{}]\frac{m^2U}{1x^2}\right]+m\stackrel{~}{\mathrm{\Omega }}^{}U=0.$$
(23)
Multiply by the complex conjugate $`U^{}`$ and integrate to get

$`\left[(1x^2)U^{}U^{}\right]_0^1{\displaystyle _0^1}(1x^2)|U^{}|^2𝑑x`$ $``$
$`m^2{\displaystyle _0^1}{\displaystyle \frac{|U^2|}{1x^2}}𝑑x+m{\displaystyle _0^1}{\displaystyle \frac{\stackrel{~}{\mathrm{\Omega }}^{}|U^2|}{m\mathrm{\Omega }\sigma }}𝑑x`$ $`=`$ $`0.`$ (24)
The first term is clearly zero due to the boundary conditions. The imaginary part of equation (5) is

$$m\sigma _i_0^1\frac{\stackrel{~}{\mathrm{\Omega }}^{}|U|^2}{|m\mathrm{\Omega }\sigma |^2}𝑑x=0.$$
(25)
For dynamically unstable modes $`\sigma _i0`$ and equation (25) can only hold if $`\stackrel{~}{\mathrm{\Omega }}^{}`$ changes sign at some point on the interval. This is a necessary (although not sufficient) condition for instability.

## 6 Corotating solutions

We begin our study of corotating solutions by considering the case when there is a corotation point in the domain $`0<x_c<1`$ (the symmetry of the problem means that there will be a second corotation point in the lower hemisphere, at $`x=x_c`$). We will consider the special cases of $`x_c=1,0`$ (corotation point at the poles or the equator) later. Employing the method of Frobenius to investigate the solutions in the vicinity of the corotation point, we let

$$u=\underset{n=0}{\overset{\mathrm{}}{}}a_n(xx_c)^{n+\gamma }.$$
(26)
For all of the rotation laws examined in this paper $`m\mathrm{\Omega }\sigma `$ contains only a single factor of $`xx_c`$, also $`\mathrm{\Omega }`$ is an analytic function of $`x`$. Assume also (for now) that $`\stackrel{~}{\mathrm{\Omega }}^{}`$ does not have a zero at $`x=x_c`$. Substituting into equation (4) to obtain the indicial equation, we find that $`\gamma =0,1`$, so we will have one regular and one singular solution. The general solution is

$$u=u_{reg}+u_{sing},$$
(27)
where

$$u_{reg}=\underset{n=0}{\overset{\mathrm{}}{}}a_n(xx_c)^{n+1},$$
(28)
and

$$u_{sing}=\underset{n=0}{\overset{\mathrm{}}{}}\left[a_n(xx_c)^{n+1}\mathrm{ln}(xx_c)+c_n(xx_c)^n\right].$$
(29)
Note that although $`u_{sing}`$ is finite and continuous at $`x_c`$ its derivatives are singular.

### 6.1 Possible regular solutions

Let us begin by asking whether the problem admits any solutions that are regular at corotation (i.e. contain $`u_{reg}`$ only). First of all note that if $`u`$ is regular in corotation then so is $`U`$. We now apply the method that was used to derive a necessary condition for instability in section 5.

We begin with equation (23), multiply by $`U^{}`$ and as before integrate. This time, however, the frequency is taken to be real and we are going to treat the two domains $`0xx_c`$ and $`x_cx1`$ separately. On the first domain

$`\left[(1x^2)U^{}U^{}\right]_0^{x_c}{\displaystyle _0^{x_c}}(1x^2)|U^{}|^2𝑑x`$ $``$
$`m^2{\displaystyle _0^{x_c}}{\displaystyle \frac{|U^2|}{1x^2}}𝑑x+m{\displaystyle _0^{x_c}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Omega }}^{}|U^2|}{m\mathrm{\Omega }\sigma }}𝑑x`$ $`=`$ $`0.`$ (30)
The first term vanishes due to boundary conditions and we are left with

$`{\displaystyle _0^{x_c}}\left[(1x^2)|U^{}|^2+{\displaystyle \frac{m^2|U^2|}{1x^2}}\right]𝑑x`$ $`+`$
$`m{\displaystyle _0^{x_c}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Omega }}^{}|U^2|}{m\mathrm{\Omega }\sigma }}𝑑x`$ $`=`$ $`0.`$ (31)
It is clear that we must have

$$m_0^{x_c}\frac{\stackrel{~}{\mathrm{\Omega }}^{}|U^2|}{m\mathrm{\Omega }\sigma }𝑑x>0.$$
(32)
(Note that $`U/(m\mathrm{\Omega }\sigma )`$ is well behaved at $`x_c`$ if $`uu_{reg}`$.) Treating the second domain in the same way we can derive a second condition, which is that we must have

$$m_{x_c}^1\frac{\stackrel{~}{\mathrm{\Omega }}^{}|U^2|}{m\mathrm{\Omega }\sigma }𝑑x>0.$$
(33)
Now $`m\mathrm{\Omega }\sigma `$ clearly changes sign at the corotation point. So unless $`\stackrel{~}{\mathrm{\Omega }}^{}`$ also changes sign somewhere on the shell it will be impossible to satisfy both conditions. A rotation law will therefore only admit purely regular solutions if $`\stackrel{~}{\mathrm{\Omega }}^{}=0`$ somewhere on the shell.

There is one other class of regular solution that might be encountered. The governing equation (4) will remain non-singular if, for a particular frequency, the zero of $`m\mathrm{\Omega }\sigma `$ coincides with that of $`\stackrel{~}{\mathrm{\Omega }}^{}`$. This will in fact always be the case for some frequency in the continuous spectrum range if $`\stackrel{~}{\mathrm{\Omega }}^{}`$ has a zero on the shell. Unlike the regular Frobenius solution $`u_{reg}`$, however, solutions in this class need not have a zero in the eigenfunction at the corotation point.

### 6.2 General (singular) solutions

The coefficients of both the regular and singular solutions (equations (28) and (29)) are found by substituting the solutions into equation (4). We find $`c_0`$ to be related to $`a_0`$,

$$c_0=\frac{\mathrm{\Omega }^{}(x_c)}{\stackrel{~}{\mathrm{\Omega }}^{}(x_c)}(x_c^21)a_0.$$
(34)
Hence $`a_0`$ and $`c_1`$ (equation 29) are the two arbitrary constants in the general power series solution.

The Frobenius expansion gives us an approximation to the general solution close to the corotation point. This solution is useful in a numerical search for solutions. However, the logarithm term causes a problem. In particular, it gives rise to a branch point at $`x_c`$ and we would in principle need to integrate around this point in the complex plane. An alternative is to seek an equivalent well-behaved solution for $`x<x_c`$ and match the two solutions at the corotation point.

Following Balbinski (1985) we make a change of variable

$$z=\frac{x_c}{x}.$$
(35)
The corotation point is now at $`z=1`$. Equation (4) becomes

$`(m\mathrm{\Omega }\sigma )[{\displaystyle \frac{(z^2x_c^2)}{x_c^2}}(z^2u^{\prime \prime }+2zu^{})+2(m+1)zu^{}`$
$`m(m+1)u]+mu\stackrel{~}{\mathrm{\Omega }}^{}=0.`$ (36)
As before, we look for a power series solution in the vicinity of the corotation point, letting

$$u=\underset{n=0}{\overset{\mathrm{}}{}}b_n(z1)^{n+\gamma }.$$
(37)
Once again the indicial equation yields $`\gamma =0,1`$, and the general solution is $`u=u_{reg}+u_{sing}`$, where

$$u_{reg}=\underset{n=0}{\overset{\mathrm{}}{}}b_n(z1)^{n+1},$$
(38)
and

$$u_{sing}=\underset{n=0}{\overset{\mathrm{}}{}}\left[b_n(z1)^{n+1}\mathrm{ln}(z1)+d_n(z1)^n\right].$$
(39)
The free parameters in the general solution are $`b_0`$ and $`d_1`$. In the region where $`x<x_c`$, $`z>1`$. So by writing the equation in terms of $`z`$ we have found a singular solution whose logarithmic terms are real in the region where $`x<x_c`$.

When solving for corotating solutions we follow Balbinski (1985) and use the two different solutions in the regions where the logarithms have positive argument, and match at the corotation point. Denoting by $`u_R`$ the solution in the region where $`x>x_c`$, and by $`u_L`$ the solution in the region where $`x<x_c`$ (and $`z>1`$), we have

$`u_R`$ $``$ $`c_0+(a_0+c_1)(xx_c)`$ (40)
$`+[a_0+a_1(xx_c)](xx_c)\mathrm{ln}|xx_c|,`$
$`u_L`$ $``$ $`d_0+(b_0+d_1)(z1)`$ (41)
$`+[b_0+b_1(z1)](z1)\mathrm{ln}|z1|`$
$``$ $`d_0+(b_0+d_1){\displaystyle \frac{(x_cx)}{x}}`$
$`+\left[b_0+b_1{\displaystyle \frac{(x_cx)}{x}}\right]{\displaystyle \frac{(x_cx)}{x}}\mathrm{ln}\left|{\displaystyle \frac{x_cx}{x}}\right|.`$
In the limit as the solutions approach the corotation point, all of the terms tend to zero apart from the constant terms. For $`u`$ to be continuous we must have

$$c_0=d_0.$$
(42)
Now consider the first derivatives.

$`{\displaystyle \frac{du_R}{dx}}`$ $``$ $`2a_0+c_1+(3a_1+2c_2)(xx_c)`$ (43)
$`+a_0\mathrm{ln}|xx_c|+2a_1(xx_c)\mathrm{ln}|xx_c|,`$
$`{\displaystyle \frac{du_L}{dx}}`$ $`=`$ $`{\displaystyle \frac{z^2}{x_c}}{\displaystyle \frac{du_R}{dz}}`$ (44)
$``$ $`{\displaystyle \frac{x_c}{x^2}}[2b_0+d_1+(3b_1+2d_2){\displaystyle \frac{(x_cx)}{x}}`$
$`+b_0\mathrm{ln}\left|{\displaystyle \frac{x_cx}{x}}\right|+2b_1{\displaystyle \frac{(x_cx)}{x}}\mathrm{ln}\left|{\displaystyle \frac{x_cx}{x}}\right|].`$
In the limit as the derivatives approach the corotation point all of the terms tend to zero apart from the constant terms and the plain logarithmic terms. We get

$`\underset{xx_c^+}{lim}{\displaystyle \frac{du_R}{dx}}`$ $`=`$ $`2a_0+c_1+a_0\underset{xx_c}{lim}\mathrm{ln}|xx_c|,`$ (45)
$`\underset{xx_c^{}}{lim}{\displaystyle \frac{du_L}{dx}}`$ $`=`$ $`{\displaystyle \frac{1}{x_c}}[2b_0+d_1+b_0\underset{xx_c}{lim}\mathrm{ln}|xx_c|`$ (46)
$`b_0\mathrm{ln}x_c].`$
Our insistence on the continuity of the function (equation (42)) means that

$$b_0=x_ca_0.$$
(47)
Recall that $`a_0`$ and $`c_0`$ are related by equation (34), and so the coefficients of the singular $`\mathrm{ln}|xx_c|`$ term are equal at the corotation point. But in addition to this logarithmic term there is also a finite step in the first derivative. By permitting a step in the first derivative there will be a solution that meets the boundary conditions for any frequency within the range of the continuous spectrum. In fact, for each frequency there will be two solutions, one of each symmetry, with a particular step in the first derivative, that meet the boundary conditions (note that for a given frequency the step size for the solutions of opposite symmetry need not be the same).

The step in the first derivative varies across the range of continuous spectrum frequencies and it can, for certain frequencies, be zero. If the step in the first derivative vanishes then

$$2a_0+c_1=\frac{1}{x_c}(2b_0+d_1b_0\mathrm{ln}x_c).$$
(48)
Using equation (47) we can rewrite equation (48) as

$$d_1=x_c(c_1+a_0\mathrm{ln}x_c).$$
(49)
Removing the step in the first derivative means that there will be no delta function in the second derivative, merely a simple pole. Throughout the rest of this paper we will refer to these solutions as zero-step solutions. Balbinski (1985) suggested that such solutions were modes in the corotation region. Using the criterion of a vanishing step in the first derivative he was able to track the propagation of two such solutions into corotation as the amount of differential rotation was increased. These solutions subsequently merged to give birth to a dynamically unstable mode. Whilst it is not clear to us that these solutions are modes in the traditional sense, they are clearly of interest if they can give rise to dynamical instabilities. They are also notable because at the corotation point the Wronskian of the two solutions being matched ($`u_R`$ and $`u_L`$) is zero.

We have seen in this section that corotating solutions will in general be singular at the corotation point. At first sight this may seem unphysical. In solving the normal mode problem, however, it is implicitly understood that our ultimate goal is the evolution of initial data. The physical time-varying solution is obtained by convolving the solutions to the normal mode problem with the initial data and integrating over the frequency range (an inverse Laplace transform). The resulting integral solution will not be singular, and the continuous spectrum will indeed have physical relevance (Van Kampen, 1955). In fact consideration of the initial value problem leads us to suspect that the zero-step solutions might exhibit behaviour that is physically distinct from the rest of the continuous spectrum. We are currently investigating this issue (Watts et al, in preparation).

## 7 Solutions on the boundary of the corotation region

As discussed in Section 3, the r-modes are the only non-trivial solutions in the case of uniform rotation. An interesting question is whether these modes enter the corotation region as the degree of differential rotation is increased. In order to do this they must somehow develop corotation points. For the simple rotation laws we consider in this paper there are two possibilities. Either the corotation point appears first at the poles, or at the equator. In order for the mode to enter the corotation region, there must exist a solution with a corotation point at the corresponding endpoint. We will consider each possibility in turn.

### 7.1 Corotation point at the pole

First consider (monotonic) rotation laws for which the angular velocity is lowest at the pole. If modes are to enter corotation there must be mode solutions with corotation points at the pole. We can try to rule out such solutions by considering the corresponding behaviour at the pole. This problem is easiest if we work with $`U`$, i.e consider equation (23).

If the corotation point is at the pole we can search for a power series expansion of the form

$$U=\underset{n}{}a_n(1x^2)^{n+\gamma }.$$
(50)
If $`\stackrel{~}{\mathrm{\Omega }}^{}(1)0`$ then the indicial equation is

$$\gamma ^2=\frac{m^2}{4}+\frac{\stackrel{~}{\mathrm{\Omega }}^{}(1)}{2\mathrm{\Omega }^{}(1)}.$$
(51)
We now need to determine what values of $`\gamma `$ are permitted by the boundary conditions. It is clear from consideration of the physical variables that we require $`\delta v_\theta `$ and $`\delta v_\phi `$ to be zero at the pole. From the continuity equation we can see that the second condition corresponds to requiring $`dU/d\theta =0`$. Now $`U=\delta v_\theta \mathrm{sin}\theta `$ and

$$\frac{dU}{d\theta }=(1x^2)^{\frac{1}{2}}\frac{dU}{dx}.$$
(52)
We therefore have two conditions to meet at the pole: $`U=0`$ and $`dU/dx`$ less singular than $`(1x^2)^{1/2}`$. These conditions immediately rule out the possibility of $`\gamma `$ being imaginary. For real $`\gamma `$ the conditions exclude the possibility of $`\gamma `$ being negative, moreover we must have $`\gamma >1/2`$. Whether or not this condition can be met will depend on the rotation law and the degree of differential rotation.

### 7.2 Corotation point at the equator

Consider now rotation laws for which the rotation rate is lowest at the equator. If modes are to cross the corotation boundary there must be mode solutions with corotation points at the equator. Consider equation (4). If the corotation point is at the equator the symmetry of the rotation law means that $`m\mathrm{\Omega }\sigma `$ will contain a factor $`x^2`$. As before we assume that $`\stackrel{~}{\mathrm{\Omega }}^{}(0)0`$ and look for power series expansions of the form

$$u=\underset{n}{}a_nx^{n+\gamma }.$$
(53)
The indicial equation yields

$$\gamma =\frac{1}{2}\left[1\pm \left(1\frac{8\stackrel{~}{\mathrm{\Omega }}^{}(0)}{\mathrm{\Omega }^{\prime \prime }(0)}\right)^{\frac{1}{2}}\right].$$
(54)
If $`\gamma `$ is real we have two possible power series, one for each value of $`\gamma `$ (which we call $`\gamma _+`$ and $`\gamma _{}`$). $`\gamma _+1/2`$ and $`\gamma _{}1/2`$. The $`\gamma _{}`$ series has a singular derivative, which would correspond to a singular $`\delta v_\phi `$. If $`\gamma `$ were imaginary we would have $`\gamma =1/2\pm i\stackrel{~}{\gamma }/2`$, so once again the derivative would be singular.

Suppose that we permit singular solutions but restrict the degree of singularity by requiring that the physical variables $`\delta v_\theta `$ and $`\delta v_\phi `$ be square integrable. This is a reasonable restriction as it enables us to define an energy measure. For this restriction to hold $`u^{}`$, and hence $`\delta v_\phi `$, must be less singular than $`x^{1/2}`$. The only permissible power series solution is therefore that with index $`\gamma _+`$.

If the maximum value of $`\gamma _+`$ is less than 1 then $`u^{}(0)`$ is singular for even this solution. This would preclude the existence of even mode solutions with equatorial corotation points (for even modes we expect $`u^{}(0)=0`$). In other words, it appears even r-modes cannot cross into corotation for certain rotation laws. The development of equatorial corotation points in odd modes for these laws remains however a possibility.

## 8 Application to specific rotation laws

The majority of the differential rotation laws examined in the literature were introduced because of their mathematical simplicity rather than being motivated by the anticipated physics. This is obviously fine as long as one is mainly interested in the properties of typical solutions, but it is clear that one ought to be somewhat cautious before drawing astrophysical conclusions from the relatively small set of rotation laws that have been considered to date. Ideally one would like to use a rotation law for a young neutron star that is grounded in physics, emerging for example from simulations of core collapse (Zwerger & Müller, 1997; Dimmelmeier, Font & Müller, 2002). At this stage, however, we are trying to deal with the technical challenges posed by the differential rotation problem. It therefore makes sense to consider simple rotation laws expressed in terms of a few parameters. Since we also wish to compare our results to previous work we have selected four rotation laws discussed in other papers. We have also considered a fifth law that exhibits two phenomena not seen for the other laws.

The first two rotation laws, the j-constant and v-constant laws, have been used extensively in the literature (see for example Hachisu (1986)). In particular, they were used in an examination of stable non-corotating modes of a differentially rotating shell (Rezzolla & Yoshida, 2001). The j-constant and v-constant laws are examples of rotation laws where the angular velocity is greatest at the pole.

The third rotation law was discussed in a paper on solar oscillations (Wolff, 1998). This is one of the few laws encountered in the literature that is grounded in observations. Wolff used this law to examine oscillations of a differentially rotating shell. This model caught our attention because Wolff’s results suggested a large number of corotating solutions that, as differential rotation increased, combined with the more familiar r-modes and disappeared, potentially giving rise to a dynamical instability. However, Wolff did not discuss either corotation points or dynamical instabilities. This law is also of interest as it is an example of a rotation law where the angular velocity is greatest at the equator.

Our fourth rotation law is a simplified version of the Wolff law and is a particular instance of a law discussed by both Balbinski (1985) and Luyten (1990a). This law is interesting because (as we will see) it exhibits both dynamical instabilities and regular solutions in corotation.

The fifth rotation law has angular velocity greatest at the pole, but unlike the j-constant and v-constant laws, exhibits dynamical instabilities and modes with equatorial corotation points.

All five laws obey the Solberg criterion for dynamical stability to axisymmetric perturbations (Tassoul, 1978). This criterion, which generalises the Rayleigh condition for convective stability in an inviscid, incompressible fluid, states that we must have

$$\frac{d}{d\varpi }(\mathrm{\Omega }^2\varpi ^4)>0,$$
(55)
where $`\varpi `$, the cylindrical radius, is given by $`\varpi =\mathrm{sin}\theta =(1x^2)^{1/2}`$ (we have set the radius of the shell equal to one).

During violent events such as core collapse or neutron star neutron star merger it is conceivable that a star might achieve an intermediate rotation state that does not obey the Solberg criterion. The analysis conducted in the next section for the shell will only be valid for a Solberg unstable law if the timescale for growth of the Solberg instability is low compared to the timescales associated with the non-axisymmetric oscillations.

### 8.1 Method

In our analysis of each rotation law we took the following steps (in our view these would also be the natural steps in the analysis of other cases):

1. Check whether the law obeys the Solberg criterion
2. Check whether the law has an inflexion point in the vorticity ($`\stackrel{~}{\mathrm{\Omega }}^{}=0`$). This will determine whether the law can, in principle, have dynamical instabilities or regular modes in corotation.
3. Solve for real-frequency modes outside corotation. This is the easiest calculation to do, and it will indicate whether modes are approaching the corotation region.
4. Solve for modes at the edge(s) of the continuous spectrum. This will provide information concerning whether modes cross into the corotation region. Furthermore, such solutions may (as we discuss later) indicate the onset of dynamical instability.
5. Search for regular corotating modes
6. Search for corotating modes with zero step in the first derivative at the corotation point. These are solutions of the kind considered by Balbinski (1985).
7. Search for dynamically unstable modes developing at points where

* Stable real-frequency modes breach the instability criterion
    * Stable real-frequency modes merge
    * Zero-step corotating solutions merge

The notion that dynamical instabilities develop when modes merge is supported by empirical evidence from other instability calculations. The dynamical bar mode instability that limits the angular velocity of rotating stars sets in at points where two non-radial modes merge. The radial instability that sets the maximum mass of a neutron star also sets in through the merger of two modes. Theoretical support for this idea is given by Schutz (1980), who notes that whenever instability sets in along a sequence the Hermitian eigenfrequency equation will have a double root. In addition, we have already discussed how Balbinski (1985) found dynamical instabilities developing at the point where two zero-step corotating solutions merged.

The numerical integrations necessary to find the modes were carried out using a variable order variable step Adams method. A solution at one value of differential rotation was then used as an initial estimate for the solution at a slightly increased or decreased value. In this way we were able to track modes through the desired range of differential rotation.

Modes without corotation points (real and complex frequencies) were located by picking a trial frequency, integrating from pole to equator and then iterating using the Newton-Raphson method until the equatorial boundary condition appropriate to the mode symmetry was met. A similar routine was used to search for solutions on the corotation boundary, although in this case because the location of the corotation point was fixed there was a fixed relationship between frequency and amount of differential rotation.

Corotating solutions were found by splitting the domain of integration into two domains: $`x_c<x1`$ and $`0x<x_c`$. As the first domain has singularities at both ends we integrated away from both singularities and evaluated the Wronskian at a matching point in the middle of the domain. On the second domain we integrated away from the singularity at $`x=x_c`$ and evaluated the function and its derivative at the equator. The desired solutions were those for which both the Wronskian on the first domain and either the function or its derivative at the equator (depending on the symmetry of the solution) vanished simultaneously. For both regular and general solutions in corotation the power series expansions of Section 6 were used to approximate the solution at $`x=x_c\pm \delta `$ (where we used $`\delta 10^5`$ in the numerical calculations). When searching for regular solutions in corotation for which $`\stackrel{~}{\mathrm{\Omega }}^{}(x_c)0`$, we fix the arbitrary constant $`a_0`$ in the power series expansion and vary the frequency. For the general singular solutions the situation is a little more complicated. In the power series expansion for $`x=x_c+\delta `$ there are two free parameters, $`a_0`$ and $`c_1`$. By fixing these and imposing the matching conditions at $`x=x_c`$ we completely determine the two free parameters in the power series solution at $`x=x_c\delta `$. In actuality the two desired functions become zero only for particular values of both the frequency and the ratio $`c_1/a_0`$ (we use $`c_1/a_0`$ rather than $`a_0/c_1`$ as our variable because there are solutions where $`c_1=0`$). We therefore fix $`a_0`$ and vary both the frequency and this ratio, using a modified version of the Powell hybrid method to locate zeroes of the two functions. Fortunately the frequencies and ratios of solutions turn out to be smooth functions of the parameters that determine the amount of differential rotation.

### 8.2 Case 1: the j-constant rotation law

The j-constant law for a shell takes the form

$$\mathrm{\Omega }=\frac{\mathrm{\Omega }_cA^2}{1+A^2x^2}.$$
(56)
The parameter $`A`$ controls the amount of differential rotation<sup>1</sup><sup>1</sup>1In the limit $`A0`$ the j-constant law tends to zero rotation. We do not discuss this limit further in this paper and hence assume $`A>0`$., uniform rotation ($`\mathrm{\Omega }=\mathrm{\Omega }_c`$) being achieved in the limit as $`A\mathrm{}`$. The rotation rate for this law is lowest at the equator. We get

$$\stackrel{~}{\mathrm{\Omega }}^{}=\frac{2\mathrm{\Omega }_cA^4(1+A^2+3x^2)}{(1+A^2x^2)^3}.$$
(57)
Since $`\stackrel{~}{\mathrm{\Omega }}^{}`$ is always positive for real $`A`$ there will be no dynamically unstable modes or regular corotating modes for this rotation law.

Figure 1 shows the frequencies of a selection of modes and corotating zero-step solutions for $`m=1`$ and $`m=2`$, respectively. Note that for any value of $`A`$ there is a singular solution at every frequency within the corotation band. In the Figures for this and the other rotation laws, however, we only show the zero-step solutions.

The frequencies derived for modes outside corotation confirm the results of Rezzolla & Yoshida (2001). As the amount of differential rotation increases, the frequencies of the r-modes decrease and approach corotation. The code crashes as $`\sigma m\mathrm{\Omega }`$ at the equator, crashing at larger values of $`A`$ for modes that at uniform rotation have larger $`l`$. Modes that have very large $`l`$ in uniform rotation are not shown, but follow the general trend illustrated in the figures. Whether modes continue to exist very close to the corotation boundary as the degree of differential rotation is increased still further is difficult to say. We have however tracked modes very close to the corotation boundary and our results suggest that modes outside corotation cease to exist below a certain value of $`A`$. Disappearance of modes as a system parameter changes is not unexpected. Several examples of eigenvalues vanishing when parameters change are given in Kato (1995).

The question of whether modes can enter the corotation region by developing equatorial corotation points was discussed in section 7.2. For this law the maximum value of $`\gamma _+`$ is 1. This means that, analytically, we were able to rule out any even solutions with equatorial corotation points for any value of $`A`$, and any odd solutions for $`A>0.354`$. This was confirmed by our numerical results. For $`A<0.354`$ we were not able to rule out analytically the possibility of odd r-modes crossing the boundary. Our numerical results suggest however that there are no odd r-mode solutions existing below $`A0.36`$. For this rotation law, therefore, we find no evidence that r-modes cross into corotation.

### 8.3 Case 2: The v-constant rotation law

Another law that has been examined in the literature is the v-constant rotation law. There are however two different laws with this name in the literature. The first, used by Hachisu (1986) and Karino & Eriguchi (2002), takes the form

$$\mathrm{\Omega }=\frac{\mathrm{\Omega }_cA}{(A^2+1x^2)^{1/2}}.$$
(58)
The second, used by Eriguchi & Müller (1985), Karino, Yoshida & Eriguchi (2001), and Rezzolla & Yoshida (2001), takes the form

$$\mathrm{\Omega }=\frac{\mathrm{\Omega }_cA}{A+(1x^2)^{1/2}}.$$
(59)
In both cases uniform rotation ($`\mathrm{\Omega }=\mathrm{\Omega }_c`$) is reached in the limit $`A\mathrm{}`$. For the first law it is not possible to have $`\stackrel{~}{\mathrm{\Omega }}^{}=0`$, so we can rule out the possibility of dynamically unstable modes or stable modes within corotation. For the second, $`\stackrel{~}{\mathrm{\Omega }}^{}`$ can be zero for $`A>1`$ or $`A<1/2`$, so we cannot a priori rule out the possibility of unstable modes. No instabilities were found, however, by our numerical analysis. In fact the modes and continuous spectrum for the two v-constant laws behave qualitatively identically to the modes and continuous spectrum of the j-constant law. For this reason we have not included figures of the results.

### 8.4 Case 3: Wolff’s solar rotation law

Wolff’s solar rotation law (Wolff, 1998) takes the form

$$\mathrm{\Omega }=2\pi (454.860.4\beta x^271.4\beta x^4).$$
(60)
The parameter $`\beta `$ measures the amount of differential rotation, uniform rotation corresponding to $`\beta =0`$. Wolff considers the range $`0<\beta <1`$ but we will consider an upper limit of $`\beta =2`$ as there are interesting phenomena in this range. The rotation rate for this law is lowest at the pole. We get

$$\stackrel{~}{\mathrm{\Omega }}^{}=2\pi (909.6+120.8\beta +132\beta x^22142\beta x^4).$$
(61)
In this case $`\stackrel{~}{\mathrm{\Omega }}^{}`$ has a zero on the domain $`0x1`$ for $`\beta >0.4815`$. Both dynamical instabilities and regular solutions in corotation are therefore possible.

Figure 2 shows the frequencies of a selection of modes and zero-step corotating solutions for $`m=1`$ and $`m=2`$. As the amount of differential rotation increases, the frequency of the each r-mode shifts and in due course all modes encounter the corotation boundary. Modes of low $`l`$ appear to approach the boundary directly, modes of higher $`l`$ seem to approach tangentially. Searching for mode solutions on the corotation boundary, as described in section 7.1, we find valid solutions at the points where the $`l=(2,3),m=1`$ and $`l=(2,3,4),m=2`$ modes encounter the boundary, suggesting that these modes can develop corotation points<sup>2</sup><sup>2</sup>2Figure 2 appears to suggest that the $`l=4,m=1`$ mode crosses the boundary. There is however no physically valid solution on the boundary. Resolving this region in more detail it is possible to see the $`l=4`$ mode and the corresponding zero-step corotating solution approaching the boundary tangentially but not crossing.. As $`m`$ increases so does the number of valid solutions on the corotation boundary. For $`m=3`$ and $`m=4`$ there are valid solutions at the point where the four lowest $`l`$ modes encounter the corotation boundary. Modes of larger $`l`$, which appear to approach tangentially, do not seem to cross the boundary.

A numerical search for modes that are regular in corotation but for which the corotation point does not coincide with the zero of $`\stackrel{~}{\mathrm{\Omega }^{}}`$ proved fruitless. For the range of $`\beta `$ examined there were no valid solutions on the domain $`x_c<x<1`$ (see the discussion in Section 8.1).

What about the possibility of regular solutions, with corotation points at the point where $`\stackrel{~}{\mathrm{\Omega }}^{}=0`$? Such solutions could only exist at one frequency for a given $`\beta `$, this frequency being a smooth function of $`\beta `$. If we represent the solution as a sum over spherical harmonics, $`U=C_l^mY_l^m`$, it is possible to derive a recursion relation for the coefficients $`C_l^m`$. A numerical search for these modes yields nothing, however, and a numerical examination of the coefficients $`C_l^m`$ shows why. The coefficients are divergent, with $`C_l^m\mathrm{}`$ as $`l\mathrm{}`$. There are therefore no regular modes in corotation for this law.

Continuous spectrum solutions with zero-step in the first derivative are present for this rotation law. Such solutions are generated at the point where real-frequency modes appear to enter corotation (i.e. where there are valid solutions on the corotation boundary). There are also many such solutions present at low differential rotation which, in contrast, seem to approach both edges of the corotation region tangentially.

Dynamical instabilities are found to develop at the point where modes encounter the corotation boundary at values of $`\beta `$ for which $`\stackrel{~}{\mathrm{\Omega }}^{}`$ has a zero on the shell (in accordance with the instability criterion derived in Section 5). Modes that encounter the corotation boundary below the threshold value of $`\beta `$ do not exhibit instabilities. There are five unstable modes, corresponding to $`l=(2,3),m=1;l=(2,3)m=2`$ and $`l=3,m=3`$. We find no instabilities for modes of higher $`l`$ or $`m`$. Frequencies and growth times for these unstable modes are shown in figure 3. The behaviour of the $`l=m=3`$ mode is particularly interesting as it indicates that the mode becomes stable again for high differential rotation. The $`l=3,m=2`$ mode exhibits similar behaviour at higher $`\beta `$.

We can see from equations (25) and (61) why we should not expect to find any unstable modes at high $`\beta `$. As $`\beta `$ increases, the negative part of $`\stackrel{~}{\mathrm{\Omega }}^{}`$ will come to dominate the integrand in equation (25). Beyond a certain threshold value it will no longer be possible to satisfy the integral equation. The precise value of the threshold will depend on $`U`$ but all modes will reach this point for some $`\beta `$.

The results reported in this paper agree with those reported by Wolff for real-frequency modes outside corotation. However, they differ dramatically for solutions with corotation points. In addition, we find dynamical instabilities not mentioned by Wolff. The discrepancy between the results likely originates with Wolff’s perturbation method and with his choice of basis eigenfunctions. He uses not only the standard r-mode eigenfunctions discussed in section 3 but a set of toroidal geostrophic eigenfunctions of single $`l,m`$. But the geostrophic solutions form a degenerate set in uniform rotation and cannot be classed as eigenfunctions. As such they cannot be used as a basis for perturbation.

We can compare the results from our investigation with measurements of the degree of differential rotation at various depths in the Sun (Howe et al, 2000). The degree of differential rotation is highest in the surface layers and the convective envelope, where values of $`\mathrm{\Omega }(0)/2\pi =450`$ nHz and $`\mathrm{\Omega }(\mathrm{cos}(\pi /6))/2\pi =370`$ nHz are reported. This would correspond to having $`\beta 0.9`$. At this level in our model we see dynamically unstable modes, the fastest having a growth time of just a few years.

Our model is very simplistic. It does not, for example, include viscosity, which would be expected to shift the points of marginal stability to a higher degree of differential rotation. However, in the outer layers of the Sun the degree of differential rotation is almost independent of the radius, suggesting that a shell model may not be a bad approximation. In this case a dynamical instability of the type found in this model could be involved in limiting the degree of differential rotation observed in the Sun. See also Dzhalilov, Staude & Oraevsky (2002) for a discussion of various Solar properties that might be explained by the action of unstable r-modes with growth times of $`101000`$ years. Whether dynamical instabilities of the type found in this paper could play a similar role is an open question.

### 8.5 Case 4: A simple quadratic rotation law

The fourth rotation law we have considered takes the form

$$\mathrm{\Omega }=\mathrm{\Omega }_c(1\beta x^2).$$
(62)
Here the parameter $`\beta `$ measures the degree of differential rotation and we have examined the range $`0\beta 1`$. Uniform rotation ($`\mathrm{\Omega }=\mathrm{\Omega }_c`$) is recovered if $`\beta =0`$. As with Wolff’s law the rotation rate for this law is greatest at the equator. We get

$$\stackrel{~}{\mathrm{\Omega }}^{}=2\mathrm{\Omega }_c(1+\beta 6\beta x^2).$$
(63)
From this we see that $`\stackrel{~}{\mathrm{\Omega }}^{}`$ has a zero on the domain if $`\beta 0.2`$. Both dynamical instabilities and regular solutions in corotation are therefore possible. Figure 4 shows the frequencies of a selection of modes and zero-step corotating solutions for $`m=1`$ and $`m=2`$.

Just as for Wolff’s law, modes of low $`l`$ approach the corotation boundary directly as differential rotation is increased. Modes of higher $`l`$ tend to approach tangentially. And in accordance with this there are valid solutions on the corotation boundary for the $`l=(2,3),m=1`$ and $`l=(2,3),m=2`$ modes.

The $`l=3`$ mode is a little different however, as it passes straight into corotation as a regular mode. To see why this happens consider the following argument. For this rotation law we can write

$$m\mathrm{\Omega }\sigma =m\beta \mathrm{\Omega }_c\left[x^2\frac{(m\mathrm{\Omega }_c\sigma )}{m\mathrm{\Omega }_c\beta }\right],$$
(64)
and

$$\stackrel{~}{\mathrm{\Omega }}^{}=12\beta \mathrm{\Omega }_c\left[x^2\frac{1+\beta }{6\beta }\right].$$
(65)
As discussed in section 6.1, the zeroes of these two functions may coincide, making the governing equation regular and giving rise to regular solutions in corotation. This will occur if

$$\sigma =\frac{m\mathrm{\Omega }_c(5\beta )}{6}.$$
(66)
In this case the governing equation (23) becomes

$$\left[(1x^2)U^{\prime \prime }2xU^{}\frac{m^2U}{1x^2}\right]+6U=0.$$
(67)
This equation is simply Legendre’s equation for $`l=3`$. Thus there is a pure $`l=3`$ mode that corotates but is completely regular. Just as for Wolff’s law, however, there are no regular corotating solutions of the type given in equation (28).

Only the $`l=2,m=1`$ and $`l=2,m=2`$ modes become dynamically unstable at the point where they encounter the corotation boundary. These are the only modes that hit the boundary after passing the threshold value of $`\beta `$ for which $`\stackrel{~}{\mathrm{\Omega }}^{}`$ first has a zero on the shell. Frequencies and growth times for both modes are shown in figure 5.

### 8.6 Case 5: A simple trigonometric rotation law

The final rotation law we have considered is

$$\mathrm{\Omega }=\mathrm{\Omega }_c\left(1\beta \mathrm{cos}x\right)$$
(68)
The parameter $`\beta `$ measures the degree of differential rotation and we have examined the range $`0\beta 1`$. Uniform rotation ($`\mathrm{\Omega }=\mathrm{\Omega }_c`$) is recovered if $`\beta =0`$. The rotation rate for this law is greatest at the pole, as it is for the j-constant and v-constant laws. What is different about this law, however, is that $`\stackrel{~}{\mathrm{\Omega }}^{}=0`$ for $`\beta >2/3`$, which could permit dynamical instabilities. Figure 6 shows mode frequencies for $`m=2`$ for a reduced range of the differential rotation parameter $`\beta `$. The figure shows interesting phenomena that we do not see for any of the other rotation laws. Modes crossing the corotation boundary for this law would have to develop equatorial corotation points, as discussed in section 7.2. For this law, the maximum value of $`\gamma _+`$ is 2. Modes of odd symmetry can in principle cross the boundary if $`\beta >16/25`$, modes of even symmetry if $`\beta >2/3`$.

Let us start by considering the $`l=3`$ mode. As the amount of differential rotation is increased, the frequency approaches the corotation boundary. The mode crosses the corotation boundary above the instability threshold value of $`\beta `$, and above this point continues to exist as a dynamically unstable mode. An odd zero-step solution also develops at this point, but this is not shown in Figure 6, as the frequency is very close to the real part of the frequency of the dynamically unstable mode.

The situation for the even modes is slightly different. Although the modes could in principle cross the corotation boundary, we find that the $`l=2,m=2`$ and $`l=4,m=2`$ modes merge outside corotation, becoming dynamically unstable. Such a route to instability is not seen for the other rotation laws examined. With reference to Figure 6, note that we were unable to track the $`l=4`$ mode for $`0.644<\beta <0.678`$ because the frequency was too close to the corotation boundary. We are however confident that it is this mode that is merging with the $`l=2`$ mode, for three reasons. Firstly, the merging mode has even symmetry (as we would expect). Secondly, all even modes with $`l>4`$ approach the corotation boundary at lower values of $`\beta `$. The $`l=4`$ mode is the only other even mode apart from $`l=2`$ that appears to persist to high $`\beta `$. Thirdly, the merging mode is not emerging from corotation, as the only valid solution on the corotation boundary for $`m=2`$ occurs at the point where the $`l=3`$ mode crosses. The most reasonable conclusion is that the $`l=4`$ mode continues to exist just outside corotation for $`0.644<\beta <0.678`$ before moving away from the boundary to merge with the $`l=2`$ mode at higher $`\beta `$.

## 9 Summary

The oscillations of neutron stars are thought to be a promising source of gravitational waves. In order to develop accurate models of gravitational wave emission we need to develop a full understanding of the dynamics of these stars. One factor that has been neglected in many previous studies is differential rotation. Neutron stars are thought to be born with differential rotation, and there are several plausible mechanisms that could sustain differential rotation beyond birth. Differential rotation introduces the possibility of dynamical shear instabilities (Papaloizou & Pringle, 1984, 1985). In addition, the perturbation equations of differentially rotating inviscid systems are complicated by the presence of singularities at points where the pattern speed of a mode equals the local angular velocity of the background flow. These singularities mean that differentially rotating inviscid systems possess a continuous spectrum in addition to discrete modes. To begin to understand this problem we have in this paper analysed the normal mode problem for a toy model, a differentially rotating shell. Although the model is simplistic, it exhibits many of the key features of the more complex three-dimensional problem.

In the uniform rotation limit the shell problem admits the standard r-modes as solutions. As the degree of differential rotation increases the frequencies of these modes gradually approach the lower edge of the expanding corotation band ($`m\mathrm{\Omega }_{min}<\sigma <m\mathrm{\Omega }_{max}`$). Modes that cross the corotation boundary are found to be rare. Indeed, we can rule out completely the crossing of even modes for certain rotation laws. For other rotation laws and symmetries we find crossings only for modes that have low $`l`$ in the uniform rotation limit. All other modes approach the boundary tangentially and appear to cease to exist when the degree of differential rotation exceeds a threshold value. Modes that have higher $`l`$ in the uniform rotation limit cease to exist at lower degrees of differential rotation.

For frequencies in the corotation band the shell problem admits a continuous spectrum of corotating solutions with discontinuous first derivatives. Although such solutions may appear at first sight unphysical, they will have physical significance when we consider the initial value problem. In general, corotating solutions have both a logarithmic singularity and a finite step in the first derivative at the corotation point. For certain special frequencies, however, the step in the first derivative vanishes. In differentially rotating cylinders these zero-step solutions can merge to give rise to dynamical instabilities (Balbinski, 1985), but we observe no such behaviour in our problem for the rotation laws examined. We cannot rule out the possibility, however, for other rotation laws.

The character of the modes that cross into corotation varies. At all points where modes cross the corotation boundary we observe the development of zero-step corotating solutions (except in one unique case, discussed below). At some points we also observe the development of dynamically unstable modes. A necessary (but not sufficient) condition for dynamical instability is that $`\stackrel{~}{\mathrm{\Omega }}^{}`$ be zero somewhere on the shell. For the rotation laws examined this condition is met only when the degree of differential rotation exceeds a threshold value. All modes crossing the corotation boundary above this threshold become dynamically unstable. We also observe, for one rotation law, the development of a dynamical instability via the merger of two stable modes outside corotation. Only in one exceptional case for the simplest rotation law did we observe a mode crossing the corotation boundary and continuing to exist as a regular stable mode. Such regular corotating modes are unlikely to exist for more realistic rotation laws.

Before moving on to the three-dimensional problem, we intend to investigate two related problems for the simple shell model. The first factor of interest is the effect of viscosity, an essential component of any realistic model. If we consider the Navier-Stokes equations instead of the Euler equations the vorticity equation for the stream function becomes

$`(\sigma m\mathrm{\Omega })_\theta ^2U+{\displaystyle \frac{m_\theta \stackrel{~}{\mathrm{\Omega }}}{\mathrm{sin}\theta }}U{\displaystyle \frac{i\nu }{\rho R^2}}[_\theta ^2(_\theta ^2U)`$ $`+`$
$`{\displaystyle \frac{m^2}{\mathrm{sin}^4\theta }}(\mathrm{sin}\theta \mathrm{cos}\theta _\theta U(\mathrm{cos}^2\theta +m^2)U)]`$ $`=`$ $`0,`$ (69)
where $`\nu `$ is the coefficient of shear viscosity. Equation (9) is not singular, as the coefficient of the leading derivative is never zero. So in the viscous problem there is no continuous spectrum, simply sets of discrete modes. In the zero viscosity limit, however, the solutions to the viscous problem should reduce to the modes and continuous spectrum of the inviscid problem. At first sight it is difficult to see how discrete well-behaved modes could reduce to the singular continuous spectrum solutions that we derived for the inviscid problem in Section 6.2. Boyd (1981) resolves this difficulty by noting that any physical continuous spectrum perturbation is made up of an integral over frequency of the singular inviscid solutions. When solving the viscous problem, the modes are discrete for even the tiniest viscosity and no integration over frequency is required. Solution of the viscous problem should therefore shed light on the physical manifestation of the continuous spectrum.

We are also pursuing the initial value problem. It is clear from our analysis of the normal mode problem that the time dependence of the continuous spectrum will only be understood by solving the initial value problem. In fact the singular continuous spectrum solutions found in this paper only acquire physical meaning when we consider the evolution of initial data. We are particularly interested to see whether the zero-step solutions appear physically distinct in time evolutions. That the zero-step solutions might exhibit different time-dependence is suggested by the fact that these solutions have zero Wronskian, which will contribute an additional pole to the integrand of the Laplace Transform inversion integral (Watts et al, in preparation). We also hope that an analysis of the initial value problem will enable us to develop necessary and sufficient conditions for instability.

## 10 Acknowledgements

We would like to thank Luciano Rezzolla, Shin’ichirou Yoshida and Emanuele Berti for helpful comments. We also acknowledge support from the EU Programme ’Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137). AW is supported by a PPARC postgraduate studentship and NA acknowledges support from the Leverhulme Trust in the form of a prize fellowship and PPARC grant PPA/G/1998/00606.


text:
# An information-theoretic analysis of Grover’s algorithm

## 1 Introduction

Grover ,  discovered a quantum algorithm for identifying a target element in an unstructured search universe of $`N`$ items in approximately $`\pi /4\sqrt{N}`$ queries to a quantum oracle. For classical search using a classical oracle, the search complexity is clearly of order $`N/2`$ queries since on average half of the items must be searched. It has been proven that this square-root speed-up is the best attainable performance gain by any quantum algorithm. In work preceding Grover’s, Bennett et al.  had shown that no quantum algorithm can solve the search problem in fewer than $`O(\sqrt{N})`$ queries. Following Grover’s work, Boyer et al.  showed that Grover’s algorithm is optimal asymptotically, and that square-root speed-up cannot be improved even if one allows, e.g., a 50% probability of error. Zalka  strengthened these results to show that Grover’s algorithm is optimal exactly (not only asymptotically). In this correspondence we present an information-theoretic analysis of Grover’s algorithm and show the optimality of Grover’s algorithm from a different point of view.

## 2 A general framework for quantum search

We consider the following general framework for quantum search algorithms. We let $`X`$ denote the state of the target and $`Y`$ the output of the search algorithm. We assume that $`X`$ is uniformly distributed over the integers $`0`$ through $`N1`$. $`Y`$ is also a random variable distributed over the same set of integers. The event $`Y=X`$ signifies that the algorithm correctly identifies the target. The probability of error for the algorithm is defined as $`P_e=P(YX)`$.

The state of the target is given by the density matrix

$$\rho _T=\underset{x=0}{\overset{N1}{}}(1/N)|xx|,$$
(1)
where $`\{|x\}`$ is an orthonormal set. We assume that this state is accessible to the search algorithm only through calls to an oracle whose exact specification will be given later. The algorithm output $`Y`$ is obtained by a measurement performed on the state of the quantum computer at the end of the algorithm. We shall denote the state of the computer at time $`k=0,1,\mathrm{}`$ by the density matrix $`\rho _C(k)`$. We assume that the computation begins at time 0 with the state of the computer given by an initial state $`\rho _C(0)`$ independent of the target state. The computer state evolves to a state of the form

$$\rho _C(k)=\underset{x=0}{\overset{N1}{}}(1/N)\rho _x(k)$$
(2)
at time $`k`$, under the control of the algorithm. Here, $`\rho _x(k)`$ is the state of the computer at time $`k`$, conditional on the target value being $`x`$. The joint state of the target and the computer at time $`k`$ is given by

$$\rho _{TC}(k)=\underset{x=0}{\overset{N1}{}}(1/N)|xx|\rho _x(k).$$
(3)
The target state (1) and the computer state (2) can be obtained as partial traces of this joint state.

We assume that the search algorithm consists of the application of a sequence of unitary operators on the joint state. Each operator takes one time unit to complete. The computation starts at time 0 and terminates at a predetermined time $`K`$, when a measurement is taken on $`\rho _C(K)`$ and $`Y`$ is obtained. In accordance with these assumptions, we shall assume that the time index $`k`$ is an integer in the range 0 to $`K`$, unless otherwise specified.

There are two types of unitary operators that may be applied to the joint state by a search algorithm: oracle and non-oracle. A non-oracle operator is of the form $`IU`$ and acts on the joint state as

$$\rho _{TC}(k+1)=(IU)\rho _{TC}(k)(IU)^{}=\underset{x}{}(1/N)|xx|U\rho _x(k)U^{}.$$
(4)
Under such an operation the computer state is transformed as

$$\rho _C(k+1)=U\rho _C(k)U^{}.$$
(5)
Thus, non-oracle operators act on the conditional states $`\rho _x(k)`$ uniformly; $`\rho _x(k+1)=U\rho _x(k)U^{}`$. Only oracle operators have the capability of acting on conditional states non-uniformly.

An oracle operator is of the form $`_x|xx|O_x`$ and takes the joint state $`\rho _{TC}(k)`$ to

$$\rho _{TC}(k+1)=\underset{x}{}(1/N)|xx|O_x\rho _x(k)O_x^{}.$$
(6)
The action on the computer state is

$$\rho _C(k+1)=\underset{x}{}(1/N)O_x\rho _x(k)O_x^{}.$$
(7)
All operators, involving an oracle or not, preserve the entropy of the joint state $`\rho _{TC}(k)`$. The von Neumann entropy of the joint state remains fixed at $`S[\rho _{TC}(k)]=\mathrm{log}N`$ throughout the algorithm. Non-oracle operators preserve also the entropy of the computer state; the action (5) is reversible, hence $`S[\rho _C(k+1)]=S[\rho _C(k)]`$. Oracle action on the computer state (7), however, does not preserve entropy; $`S[\rho _C(k+1)]S[\rho _C(k)]`$, in general.

Progress towards identifying the target is made only by oracle calls that have the capability of transferring information from the target state to the computer state. We illustrate this information transfer in the next section.

## 3 Grover’s algorithm

Grover’s algorithm can be described within the above framework as follows. The initial state of the quantum computer is set to

$$\rho _C(0)=|ss|$$
(8)
where

$$|s=\underset{x=0}{\overset{N1}{}}(1/\sqrt{N})|x.$$
(9)
Since the initial state is pure, the conditional states $`\rho _x(k)`$ will also be pure for all $`k1`$.

Grover’s algorithm uses two operators: an oracle operator with

$$O_x=I2|xx|,$$
(10)
and a non-oracle operator (called ‘inversion about the mean’) given by $`IU_s`$ where

$$U_s=2|ss|I.$$
(11)
Both operators are Hermitian.

Grover’s algorithm interlaces oracle calls with inversion-about-the-mean operations. So, it is convenient to combine these two operations in a single operation, called Grover iteration, by defining $`G_x=U_sO_x`$. The Grover iteration takes the joint state $`\rho _{TC}(k)`$ to

$$\rho _{TC}(k+1)=\underset{x}{}(1/N)|xx|G_x\rho _x(k)G_x^{}$$
(12)
In writing this, we assumed, for notational simplicity, that $`G_x`$ takes one time unit to complete, although it consists of the succession of two unit-time operators.

Grover’s algorithm consists of $`K=(\pi /4)\sqrt{N}`$ successive applications of Grover’s iteration beginning with the initial state (8), followed by a measurement on $`\rho _C(K)`$ to obtain $`Y`$. The algorithm works because the operator $`G_x`$ can be interpreted as a rotation of the $`x`$$`s`$ plane by an angle $`\theta =\mathrm{arccos}(12/N)2/\sqrt{N}`$ radians. So, in $`K`$ iterations, the initial vector $`|s`$, which is almost orthogonal to $`|x`$, is brought into alignment with $`|x`$.

Grover’s algorithm lends itself to exact calculation of the eigenvalues of $`\rho _C(k)`$, hence to computation of its entropy. The eigenvalues of $`\rho _C(k)`$ are

$$\lambda _1(k)=\mathrm{cos}^2(\theta k)$$
(13)
of multiplicity 1, and

$$\lambda _2(k)=\frac{\mathrm{sin}^2(\theta k)}{N1}$$
(14)
of multiplicity $`N1`$. The entropy of $`\rho _C(k)`$ is given by

$$S(\rho _C(k))=\lambda _1(k)\mathrm{log}\lambda _1(k)(N1)\lambda _2(k)\mathrm{log}\lambda _2(k)$$
(15)
and is plotted in Fig. 1 for $`N=2^{20}`$. (Throughout the paper, the unit of entropy is bits and $`\mathrm{log}`$ denotes base 2 logarithm.) The entropy $`S(\rho _c(k))`$ has period $`\pi /\theta (\pi /2)\sqrt{N}`$.

Our main result is the following lower bound on time-complexity.

###### Proposition 1

Any quantum search algorithm that uses the oracle calls $`\{O_x\}`$ as defined by (10) must call the oracle at least

$$K\left(\frac{1P_e}{2\pi }+\frac{1}{\pi \mathrm{log}N}\right)\sqrt{N}$$
(16)
times to achieve a probability of error $`P_e`$.

For the proof we first derive an information-theoretic inequality. For any quantum search algorithm of the type described in section 2, we have by Fano’s inequality,

$$H(Y|X)(P_e)+P_e\mathrm{log}(N1)(P_e)+P_e\mathrm{log}(N),$$
(17)
where for any $`0u1`$

$$(u)=\delta \mathrm{log}\delta (1\delta )\mathrm{log}(1\delta ).$$
(18)
On the other hand,

$`H(X|Y)`$ $`=`$ $`H(X)I(X;Y)`$ (19)
$`=`$ $`\mathrm{log}NI(X;Y)`$
$``$ $`\mathrm{log}NS(\rho _C(K))`$
where in the last line we used Holevo’s bound \[6, p. 531\].

Let $`\mu _k`$ be the largest eigenvalue (sup-norm) of $`\rho _C(k)`$. We observe that $`\mu _k`$ begins at time 0 with the value 1 and evolves to the final value $`\mu _K`$ at the termination of the algorithm. We have

$$S(\rho _C(K))\mu _K\mathrm{log}\mu _K(1\mu _K)\mathrm{log}[(1\mu _K)/(N1)](\mu _K)+(1\mu _K)\mathrm{log}N.$$
(20)
since the entropy is maximized, for a fixed $`\mu _K`$, by setting the remaining $`N1`$ eigenvalues equal to $`(1\mu _K)/(N1)`$. Combining (19) and (20),

$$\mu _K\mathrm{log}NP_e\mathrm{log}N+(P_e)+(\mu _K)P_e\mathrm{log}N+2$$
(21)
Now, let

$$\mathrm{\Delta }=sup\{|\mu _{k+1}\mu _k|:k=0,1,\mathrm{},K1\}.$$
(22)
This is the maximum change in the sup norm of $`\rho _C(k)`$ per algorithmic step. Clearly,

$$K\frac{1\mu _K}{\mathrm{\Delta }}.$$

Using the inequality (21), we obtain

$$K\frac{1P_e+2/\mathrm{log}N}{\mathrm{\Delta }}.$$
(23)
Thus, any upper bound on $`\mathrm{\Delta }`$ yields a lower bound on $`K`$. The proof will be completed by proving

###### Lemma 1

$`\mathrm{\Delta }2\pi /\sqrt{N}`$.

We know that operators that do not involve oracle calls do not change the eigenvalues, hence the sup norm, of $`\rho _C(k)`$. So, we should only be interested in bounding the perturbation of the eigenvalues of $`\rho _C(k)`$ as a result of an oracle call. We confine our analysis to the oracle operator (10) that the Grover algorithm uses.

For purposes of this analysis, we shall consider a continuous-time representation for the operator $`O_x`$ so that we may break the action of $`O_x`$ into infinitesimal time steps. So, we define the Hamiltonian

$$H_x=\pi |xx|$$
(24)
and an associated evolution operator

$$O_x(\tau )=e^{i\tau H_x}=I+(e^{i\pi \tau }1)|xx|.$$

The operator $`O_x`$ is related to $`O_x(\tau )`$ by $`O_x=O_x(1)`$.

We extend the definition of conditional density to continuous time by

$$\rho _x(k_0+\tau )=O_x(\tau )\rho _x(k_0)O_x(\tau )^{}$$
(25)
for $`0\tau 1`$. The computer state in continuous-time is defined as

$$\rho _C(t)=\underset{x}{}(1/N)\rho _x(t).$$
(26)
Let $`\{\lambda _n(t),u_n(t)\}`$, $`n=1,\mathrm{},N`$, be the eigenvalues and associated normalized eigenvectors of $`\rho _C(t)`$. Thus,

$$\rho _C(t)|u_n(t)=\lambda _n(t)|u_n(t),u_n(t)|\rho _C(t)=\lambda _n(t)u_n(t)|,u_n(t)|u_m(t)=\delta _{n,m}.$$
(27)
Since $`\rho _C(t)`$ evolves continuously, so do $`\lambda _n(t)`$ and $`u_n(t)`$ for each $`n`$.

Now let $`(\lambda (t),u(t))`$ be any one of these eigenvalue-eigenvector pairs. By a general result from linear algebra (see, e.g., Theorem 6.9.8 of Stoer and Bulirsch \[7, p. 389\] and the discussion on p. 391 of the same book),

$$\frac{d\lambda (t)}{dt}=u(t)|\frac{d\rho _C(t)}{dt}|u(t).$$
(28)
To see this, we differentiate the two sides of the identity $`\lambda (t)=u(t)|\rho _C(t)|u(t)`$, to obtain

$`{\displaystyle \frac{d\lambda (t)}{dt}}`$ $`=`$ $`u^{}(t)|\rho _C(t)|u(t)+u(t)|{\displaystyle \frac{d\rho _C(t)}{dt}}|u(t)+u(t)|\rho _C(t)|u^{}(t)`$
$`=`$ $`u(t)|{\displaystyle \frac{d\rho _C(t)}{dt}}|u(t)+\lambda (t)[u^{}(t)|u(t)+u(t)|u^{}(t)]`$
$`=`$ $`u(t)|{\displaystyle \frac{d\rho _C(t)}{dt}}|u(t)+\lambda (t){\displaystyle \frac{d}{dt}}u(t)|u(t)`$
$`=`$ $`u(t)|{\displaystyle \frac{d\rho _C(t)}{dt}}|u(t)`$
where the last line follows since $`u(t)|u(t)1`$.

Differentiating (26), we obtain

$$\frac{d\rho _C(t)}{dt}=\underset{x}{}(i/N)[H_x,\rho _x(t)]$$
(29)
where $`[,]`$ is the commutation operator. Substituting this into (28), we obtain

$`\left|{\displaystyle \frac{d\lambda (t)}{dt}}\right|`$ $`=`$ $`\left|u(t)|{\displaystyle \frac{i}{N}}{\displaystyle \underset{x}{}}[H_x,\rho _x(t)]|u(t)\right|`$
$``$ $`{\displaystyle \frac{2}{N}}\left|{\displaystyle \underset{x}{}}u(t)|H_x\rho _x(t)|u(t)\right|`$
$`\stackrel{(a)}{}`$ $`{\displaystyle \frac{2}{N}}\sqrt{{\displaystyle \underset{x}{}}u(t)|H_x^2|u(t)}\sqrt{{\displaystyle \underset{x}{}}u(t)|\rho _x^2(t)|u(t)}`$
$`\stackrel{(b)}{=}`$ $`{\displaystyle \frac{2}{N}}\sqrt{{\displaystyle \underset{x}{}}\pi ^2\left|u(t)|x\right|^2}\sqrt{Nu(t)|\rho _C(t)|u(t)}`$
$`=`$ $`{\displaystyle \frac{2\pi }{\sqrt{N}}}\mathrm{\hspace{0.33em}1}\sqrt{\lambda (t)}`$
$``$ $`{\displaystyle \frac{2\pi }{\sqrt{N}}}`$
where $`(a)`$ is the Cauchy-Schwarz inequality, $`(b)`$ is due to (i) $`\rho _x^2(t)=\rho _x(t)`$ as it is a pure state, and (ii) the definition (26). Thus,

$$|\lambda (k_0+1)\lambda (k_0)|=\left|_{k_0}^{k_0+1}\frac{d\lambda (t)}{dt}𝑑t\right|2\pi /\sqrt{N}.$$
(30)
Since this bound is true for any eigenvalue, the change in the sup norm of $`\rho _C(t)`$ is also bounded by $`2\pi /\sqrt{N}`$.

### Discussion

The bound (16) captures the $`\sqrt{N}`$ complexity of Grover’s search algorithm. As mentioned in the Introduction, lower-bounds on Grover’s algorithm have been known before; and, in fact, the present bound is not as tight as some of these earlier ones. The significance of the present bound is that it is largely based on information-theoretic concepts. Also worth noting is that the probability of error $`P_e`$ appears explicitly in (16), unlike other bounds known to us.