format: text

text:
# The Stability of Disks in Cusped Potentials1footnote 11footnote 1Rutgers Astrophysics Preprint #252

## 1 Introduction

The absence of strong bars in a significant fraction of disk galaxies has long been thought to present a serious problem for galactic dynamicists. Many simple models of fully self-gravitating disks possess global, disruptive instabilities that form bars on a dynamical time scale (Miller, Prendergast & Quirk 1970; Hohl 1971). In an influential paper, Ostriker & Peebles (1973) argued that unbarred galaxies were stable because they contained a large fraction of mass in a dynamically hot component which was unable to participate in collective instabilities. Their paper has given rise to the impression that massive dark matter halos are required for the stability of all disk galaxies.

Toomre (1981) propounded a convincing mechanism for bar-forming modes. Not only did he provide a physical understanding of the stabilizing effect of halos for disks with gently rising rotation curves, but he also argued that disks could be stabilized by making their centers inhospitable to density waves. One way to achieve this is for the circular speed to remain high towards the center, which forces an inner Lindblad resonance (ILR) that cuts the feedback loop to the swing amplifier (e.g. Binney & Tremaine 1987 §6.3). Thus, even a fully self-gravitating disk can be stabilized against bisymmetric modes, much as Toomre’s student Zang (1976) had found for the infinite, $`V=`$const. disk (here abbreviated as the Mestel disk, following the somewhat loose designation by Binney & Tremaine).

Toomre’s prediction was almost immediately challenged by Efstathiou, Lake & Negroponte (1982), who found no evidence for stabilization by hard centers in an extensive set of numerical simulations. The contradiction between their numerical results and Toomre’s prediction was resolved by Sellwood (1989) who showed that large-amplitude disturbances in highly responsive disks could saturate the ILR, trapping particles in a large-scale bar similar to that which would have formed without the dense center. Disturbances of sufficient amplitude to trigger a bar in this non-linear manner may well have arisen from particle noise in the rather grainy simulations by Efstathiou et al., particularly when the disk is highly responsive. As a result, the stabilizing effect of a dense center has been often discounted, giving rise to the flawed, but frequently recited, argument that since maximum disks are unstable, real galaxies cannot be maximum. (A galaxy model is described as having a “maximum disk” when the dark matter contribution to the central attraction in the inner regions is negligible compared with that of the luminous disk and bulge.)

Sellwood (1985, 1999) and Sellwood & Moore (1999) report $`N`$-body simulations of counter-examples which appeared robustly bar-stable. These realistic galaxy models having cool, maximal disks with dense centers, are too complicated for their stability to be checked by linear analysis, however. The simplest models with hard centers are disks in power-law potentials, and their linear stability properties have already been calculated. For example, the global, small-amplitude modes in the Mestel disk were already determined by Zang (1976) who found, rather surprisingly, that bisymmetric ($`m=2`$) modes were easily stabilized even when the disk was very cool. Zang’s result also holds in other disks with arbitrary power-law rotation curves (Evans & Read 1998a,b). Unfortunately, all these power-law models suffer from strong one-armed instabilities instead, provoking Zang to conclude “we may merely have jumped from the frying pan into the fire.” Toomre (1981) later pointed out that $`m=1`$ instabilities can be controlled by reducing the disk surface density without changing the potential, much in the same manner as bars can be avoided in soft-centered disks by making them sub-maximal.

Our objectives in this paper are two-fold. To confirm the extraordinary bisymmetric stability of the power-law disks, and to show that the lop-sided modes can also be quelled by tapering the outer disk. Such tapering again reduces the surface density everywhere, but changing the shape of the radial surface density profile from a power law to roughly exponential increases the peak disk contribution to the central attraction, preserving a maximum disk.

We begin by demonstrating reasonable agreement between results from $`N`$-body simulations and predictions from the linear stability analyses of power-law disks. Previous work has shown that the eigenmodes predicted by linear theory can be successfully reproduced in $`N`$-body experiments. These tests were for the Kalnajs disks (Sellwood 1983), the Kuz’min-Toomre disk (Sellwood & Athanassoula 1986), and the isochrone (Earn & Sellwood 1995). However, all three tests were of models with uniform density cores and therefore finite central frequencies, and we quickly discovered that simulations of power law disks require much more care in order to reproduce the predicted modes. While we are ultimately able to claim success, our difficulties arise directly from this physical difference and illustrate interesting aspects of disk dynamics.

Bisymmetric instabilities in power-law disks are provoked by a sharp inner cut-out. Not only are we able confirm these instabilities, but we also verify that an almost full-mass disk with a sufficiently gentle inner edge is bisymmetrically stable.

We demonstrate that the lop-sided instability can be removed by tapering away the massive, extended outer disk. We have been able to adapt a power-law disk to construct a completely stable and reasonably realistic galaxy model in which the disk provides most of the central attraction in the inner parts. This model has a quasi-exponential disk and an almost flat rotation curve that arises from a combination of the massive disk, a central bulge, together with a halo having very low central density. We show, using both linear stability theory and numerical simulations, that it has no global instabilities whatsoever.

## 2 The Equilibria

Motivated by the near-flatness of galactic rotation curves, we here study a family of models which possess rotation speeds that vary as a power of the radius. We break their scale-freeness by altering their surface density profiles using inner and outer tapers without changing the gravitational field.

The circular velocity is

$$v_{\mathrm{circ}}^{}{}_{}{}^{2}=\left(\frac{R_0}{R}\right)^\beta v_0^2,$$
(1)
where $`v_0`$ is the velocity at the reference radius $`R_0`$. When $`\beta =0`$, the rotation curve is completely flat; otherwise, the rotation curve is rising if $`\beta <0`$ and falling if $`\beta >0`$. The gravitational potential in the plane of the disk is

$$\psi (R)=\{\begin{array}{cc}\frac{v_0^2}{\beta }\left(\frac{R_0}{R}\right)^\beta \hfill & \beta 0\hfill \\ v_0^2\mathrm{ln}\left(\frac{R}{R_0}\right)\hfill & \beta =0.\hfill \end{array}$$
(2)
The surface density of a full-mass disk is

$$\mathrm{\Sigma }(R)=\frac{v_0^2}{2\pi GR_0(\beta )}\left(\frac{R_0}{R}\right)^{1+\beta },$$
(3)
with

$$(\beta )=\frac{\mathrm{\Gamma }\left[\frac{1}{2}\left(1\beta \right)\right]\mathrm{\Gamma }\left[\frac{1}{2}\left(2+\beta \right)\right]}{\mathrm{\Gamma }\left[\frac{1}{2}\left(1+\beta \right)\right]\mathrm{\Gamma }\left[\frac{1}{2}\left(2\beta \right)\right]}.$$
(4)
Note that $`(\beta )=1`$ for the Mestel disk, since $`\beta =0`$.

We will need the swing-amplification parameter $`X`$, defined by Toomre (1981) for an $`m`$-armed disturbance as

$$X=\frac{2\pi R}{m\lambda _{\mathrm{crit}}}=\frac{(2\beta )(\beta )}{m}.$$
(5)
Here $`\lambda _{\mathrm{crit}}`$ is given by

$$\lambda _{\mathrm{crit}}=\frac{4\pi ^2G\mathrm{\Sigma }}{\kappa ^2},$$
(6)
which local theory predicts is the longest unstable wavelength of axisymmetric Jeans instabilities in a purely rotationally supported disk (Toomre 1964). Here, $`\kappa `$ is the radial epicyclic frequency, as usual.

The distribution functions (hereafter DFs) depend on both isolating integrals of motion: the energy per unit mass $`E`$ and the specific angular momentum $`L_z`$. The DFs are

$$f_\mathrm{s}(E,L_z)=\{\begin{array}{cc}\stackrel{~}{C}L_z^\gamma |E|^{1/\beta +\gamma /\beta \gamma /2}\hfill & \beta 0\hfill \\ \stackrel{~}{C}L_z^\gamma \mathrm{exp}\left[(\gamma +1)E/v_0^2\right],\hfill & \beta =0\hfill \end{array}$$
(7)
as given by Evans (1994) for the general case and by Toomre (1977, see also Binney & Tremaine 1987) for the Mestel disk. In both these formulae, $`\stackrel{~}{C}`$ is a normalization constant, whose value can be found in Evans (1994) or Evans & Read (1998a). The velocity anisotropy constant $`\gamma `$ prescribes the radial velocity dispersion $`\sigma _u`$ of the stars:

$$\sigma _u^2=\frac{v_0^2}{1+\gamma +2\beta }\left(\frac{R_0}{R}\right)^\beta .$$
(8)
For these full-mass, self-similar disks, the stability parameter

$$Q_\mathrm{s}=\frac{2\pi (\beta )}{3.36}\sqrt{\frac{2\beta }{1+\gamma +2\beta }},$$
(9)
is independent of radius. Local theory (Toomre 1964) tells us that the disk is everywhere stable to axisymmetric modes provided $`Q_\mathrm{s}1`$.

The DFs (7) generate the mass density implied by the gravitational potential through Poisson’s equation. A self-similar disk with no preferred scale has unusual stability properties (Goodman & Evans 1999), and in any case differs greatly from the quasi-exponential profiles of real galactic disks. We therefore taper the disks with both inner and outer cut-outs, very much akin to Zang’s (1976) pioneering investigation. The cut-out mass is still present, in the sense that it contributes to the forces experienced by the remaining stars, but it is not free to participate in any perturbations. The actual DF we use is that for the self-similar case multiplied by a taper function $`H(L_z)`$:

$$f_\mathrm{c}(E,L_z)=H(L_z)f_\mathrm{s}(E,L_z),$$
(10)
where

$$H(L_z)=\frac{L_z^{\stackrel{~}{\nu }}}{\left[L_z^{\stackrel{~}{\nu }}+\left(v_0R_0\right)^{\stackrel{~}{\nu }}\right]}\frac{L_\mathrm{c}^{\stackrel{~}{\mu }}}{\left[L_z^{\stackrel{~}{\mu }}+L_\mathrm{c}^{\stackrel{~}{\mu }}\right]}.$$
(11)
The index $`\stackrel{~}{\nu }`$ controls the inner cut-out, while $`\stackrel{~}{\mu }`$ determines the outer taper; the larger the value of either index, the more abrupt the tapering. The outer taper is centered on a characteristic angular momentum $`L_\mathrm{c}`$, or radius $`R_\mathrm{c}`$, where $`L_\mathrm{c}=R_\mathrm{c}v_{\mathrm{circ}}(R_\mathrm{c})`$. In simulations, we are free to choose whatever $`\stackrel{~}{\nu }`$ and $`\stackrel{~}{\mu }`$ we wish. However, the normal mode analysis is available only for the special values

$$\stackrel{~}{\nu }=\nu \frac{2+\beta }{2\beta },\stackrel{~}{\mu }=\mu \frac{2+\beta }{2\beta },$$
(12)
where $`\nu `$ and $`\mu `$ are integers. These preferred values greatly simplify the contour integral in the normal mode analysis (see Appendix C of Evans & Read 1998a). Note that for the important $`\beta =0`$ case, there is no difference between $`\stackrel{~}{\nu }`$ and $`\nu `$, $`\stackrel{~}{\mu }`$ and $`\mu `$.

The equilibrium and linear stability properties of these models are discussed in more detail in Evans & Read (1998a,b). For the remainder of this paper, we adopt a system of units in which $`G=v_0=R_0=1`$. Unless otherwise stated, the models in our simulations have $`\nu =4`$, $`\mu =6`$, $`R_\mathrm{c}=10`$ and $`Q_\mathrm{s}=1`$.

## 3 Numerical Techniques

### 3.1 Aliasing

The acceptable agreement reported in §§4.1 & 5 between linear theory predictions and our results was reached only after many failures. In the light of the successful experience of Earn & Sellwood (1996), we began by using the smooth-field-particle (SFP) method to determine forces. While less efficient than particle-mesh (PM) methods, particularly for larger numbers of basis functions, it does not require forces to be softened at all.

Results with first 10 – 25 members of the $`k_{\mathrm{AJ}}=7`$ set of Abel-Jacobi functions (Kalnajs 1976) were appalling. The $`m=2`$ coefficients grew in a largely incoherent manner at a rate roughly ten times higher than that predicted, with “pattern speeds” which seemed quasi-continuous over a broad swathe of values. Non-axisymmetric features grew somewhat less rapidly in simulations with larger numbers of functions, but still no single “mode” stood out.

After much experimentation with different time steps, particle numbers, etc., which yielded little in the way of improved behavior, we switched to Bessel functions as our basis set. Unlike the previous discrete basis, a Hankel transform has a continuous distribution of radial wavenumbers, $`k`$, which in practice has to be represented by a discrete set of values. While our first experiments revealed no significant improvement, we soon noticed that the average rate of growth correlated with the value of $`\delta k`$ adopted; i.e. the smaller $`\delta k`$, the lower the growth rate we observed. Furthermore, inspection of the Fourier transform of the best fitting “mode” revealed a leading signal in excess of that predicted by linear theory, which gradually diminished as $`\delta k`$ was reduced.

This behavior was a strong hint of aliasing. The transforms of somewhat tightly wrapped waves, both trailing and leading, are inadequately represented by a loose spacing in $`\delta k`$. Thus, as noise induced waves sheared towards the trailing side, the gravitational field acquired components of enhanced leading waves. This purely numerical feedback from trailing to leading was clearly responsible for the unphysically rapid growth of non-axisymmetric structure in our models. The required tighter spacing in $`k`$ led to rapidly escalating cpu time requirements. Computing significantly more than 100 functions, as our results indicated would be needed, rendered this method uneconomic. (It may perhaps be possible to reproduce the predicted behavior with a much smaller number of functions from some other basis set.)

Reverting to a PM code immediately avoided such problems, and led to much more physically reasonable behavior. However, particle number, grid size, time step and softening all needed to be substantially improved over values that seemed adequate in previous work before we were able to obtain results in reasonable agreement with the predictions.

### 3.2 Adopted method

We used a 2-D polar grid (Sellwood 1981) to determine the forces from the particles. This PM method is highly efficient but does require forces to be softened, if only slightly, in order to prevent relaxation (Rybicki 1972).

We are here concerned mostly with linear instabilities, for which the radial distribution of mass is assumed to remain unchanged. For these cases, therefore, we need not evaluate the axisymmetric part of the field from the particles, and use instead the analytic central attraction of the infinite self-similar model. In fact, as the different sectoral harmonics of the perturbing forces from small amplitude disturbances are also assumed to be decoupled in linear theory, the particles in our linearized simulations contribute to just a single azimuthal Fourier component of the force field.

As these power-law models have a large range of orbital frequencies from the center to the edge, it is advantageous to employ a range of time steps. The scheme adopted (Sellwood 1985) divides the model into a number (here 5) of radial zones, with the time step increasing by a factor 2 between zones. Forces arising from all the particles in one zone are evaluated at every step and stored separately. They are combined to find the total force from the entire disk, using linear interpolation between the nearest pair of times whenever forces are needed in the inner zones at a time not coincident with an outer step. This scheme is exactly time-reversible, therefore.

As even this short step becomes inadequate near the center, we further subdivide the time step whenever a particle moves closer to the center than a “guard radius.” For such a particle, we evaluate the force acting, which is anyway dominated by the fixed axisymmetric component this close to the center, by assuming that no others move until the next step of the innermost zone. The guard radius is always equal to $`1/6`$, within which the time step is shorter by a factor 10.

The principal result from each simulation is an estimate of the eigenfrequency of the dominant mode. We determine these frequencies using the mode fitting procedure described in detail by Sellwood & Athanassoula (1986); see also Earn & Sellwood (1995). As usual, our estimated uncertainties indicate the entire range that encompasses all credible fits to subsets of the data from that run, with several scaling factors.

### 3.3 Selection of Particles

The implementation of quiet start procedures for galaxy simulations has developed from the initial ideas outlined by Sellwood (1983), and the latest prescription is described in appendix B of Debattista & Sellwood (2000). Since the DF is expressed as a function of the two classical integrals, $`f(E,L_z)`$, we aim to create a sample of particles which has a density in this space as close as possible to $`f`$. We proceed by slicing accessible $`(E,L_z)`$ space into $`j`$ small areas in such a way that the integral of the DF over each area encloses a fraction $`1/j`$ of the total active mass. We then select (not entirely at random) a point within each of these areas to determine the $`(E,L_z)`$ values for an orbit. For each orbit selected in this way, we then choose $`l`$ radial phases to determine both the initial radial position and the azimuthal and radial velocity components of each particle. We finally place $`n`$ particles having these three phase space coordinates at equally spaced azimuthal locations. Thus, the total number of particles $`N=jln`$. While as deterministic as we could make it, this procedure requires us to chose the shapes of each of the $`j`$ areas in $`(E,L_z)`$-space, and the values of $`l`$ and $`n`$. As a result, we are able to create many different samples from the DF by varying these parameters.

It seemed possible that the behavior of the modes might depend sensitively on the representation of the inner cut-out (see §4.2). If this were true and all particles have equal masses, discreteness effects would ultimately become important in this region no matter how large the number of particles employed. We therefore experimented with packing more particles into this region while retaining the desired density distribution by assigning individual masses. In most cases, we selected particles as if there were no inner cut-out, but we also tried to bias the DF with an extra factor (a function either of $`E`$ or of $`L_z`$) in order to place a still greater fraction of the particles in the inner disk. While seemingly desirable, this strategy actually back-fired, especially for models with declining rotation curves. By packing many low-mass particles in the central region, we have fewer left to represent the massive outer part of the disk. Results for declining rotation curves were clearly better when all particles had equal mass, hinting that noise in the massive part of the disk, rather than representation of the inner edge, is the principal source of error.

## 4 Bisymmetric Modes

### 4.1 Comparison with Linear Stability Theory

Zang (1976) showed that the full-mass Mestel disk is stable to $`m=2`$ modes unless the central cut-out index $`\nu >2`$. For a model with $`\nu =4`$ and $`Q_\mathrm{s}=1`$, he predicted the eigenfrequency of the dominant $`m=2`$ mode to be $`m\mathrm{\Omega }_\mathrm{p}+is=0.879+0.127i`$; where $`s`$ is the growth rate and $`\mathrm{\Omega }_\mathrm{p}`$ the pattern speed. As simulations of an infinite disk are impossible, we add an outer taper with index $`\mu =6`$ centered at $`R_\mathrm{c}=10`$, which changes the predicted frequency very slightly to $`2\mathrm{\Omega }_p+is=0.879+0.128i`$. In order to keep the radial extent finite, we impose an additional energy truncation such that no particle’s orbit takes it outside $`R=15`$. This extra truncation removes little mass and is therefore unlikely to affect the eigenfrequency further.

Figure 1 summarizes the results from many simulations of this one model as the numerical parameters are varied. Each pair of panels shows how the measured pattern speed and growth rate vary as one numerical parameter is changed from the “standard” values given in Table 1. The horizontal dotted lines indicate the frequency predicted from linear theory.

As the number of particles is increased (panels a), the measured frequencies clearly approach the predicted value near $`N=2.5\times 10^6`$. The scatter results from making many different selections of particles from the DF for each value of $`N`$. The steady trends seen in panels (b), (c) & (d), on the other hand, result from models which all began with exactly the same sample of particles.

The data we fit to estimate the eigenfrequency are never well approximated by a single mode, and all our fits require a second “mode.” Unlike the dominant mode, however, the second wave does not have a coherent frequency throughout its growth, neither is its frequency similar from runs with different samples from the DF. Since amplified particle noise can be reasonably well approximated by a coherent wave for short time periods, our fitting procedure is never able to separate noise completely from the true mode. The scatter in panels (a) and the lack of scatter in the other panels strongly suggest that our measurements are most affected by the noise spectrum in the simulation, which is changed with every new sample drawn from the DF. We discuss particle noise in more detail in §4.3.

The trends all show that the results have yet to “converge” as each parameter is refined. Our standard parameters (Table 1) are not perfect, but seem adequate. It may seem worrisome that the trends in panels (b), (c) & (d) are not converging to the prediction. We stress, however, that these results were obtained from a single sample of particles, and the true error bars on these results should be enlarged considerably to take account of the $`510`$% scatter for this $`N`$ in panels (a).

The scatter in panels (a) seems to decrease as $`N`$ rises suggesting that the predicted values might be reliably reproduced in still larger calculations, as they should. We obviously should expect our higher quality results to cluster around the prediction, which we claim they do.

The predictions and results for several values of the inner cut-out index $`\nu `$ are shown in Figure 2. The error bars on the data points show the 1-$`\sigma `$ dispersion (not the standard error) from a number ($`4`$) samples from the DF in each case. All these simulations used the standard numerical parameters given in Table 1. We find reasonable agreement between linear theory and the simulation results, but generally far below the level Earn & Sellwood (1995), who employed one twentieth the number of particles, were able to achieve for the isochrone disk. Even with the larger $`N`$ used here, no bisymmetric mode, however strongly growing, stands out clearly from the noise in the manner observed routinely in simulations of soft-centered disks. The values we obtain have a tendency to underestimate the prediction, of $`\mathrm{\Omega }_p`$ in particular, which we find bothersome; we have been unable to eliminate this small systematic difference and believe residual noise to be the most likely culprit.

Our results for two different power-law models, one with rising and the other with falling rotation curves, are summarized in Table 2, together with the predictions from linear theory, for $`\nu =4`$ in each case. The agreement for the model with the rising rotation curve ($`\beta =0.25`$) is satisfactory, but that for the falling rotation curve model is much poorer. Further investigation seems to indicate that the number of particles needed for good agreement rises with increasing $`\beta `$. One part of the reason for this behavior was traced to the folly of filling the inner cut-out with a large fraction of low-mass particles, as noted in §3.3. But this is clearly not the full story, since the results given in the Table 2 (for this case only) are for equal mass particles. While slowly-growing modes are clearly more difficult to reproduce, §4.3 gives a further reason why the falling rotation curve case is worse.

### 4.2 Mode mechanism

The reason the modes of soft-centered disks are so much easier to reproduce than those of the power-law disks stems from the difference in reflection mechanisms between the two types of disk. In both cases, they are swing-amplified inner cavity modes (Toomre 1981) in which in-going trailing waves return as out-going leading waves providing positive feedback to the amplifier. A wave-train in a soft-centered disk, such as the isochrone, reflects off the disk center and returns as a leading wave without attenuation. The centers of the power-law disks, on the other hand, cannot support density waves for two reasons: first, frequencies become arbitrarily large as $`R0`$, which ensures an ILR at some radius for all pattern speeds, and second, the inner cut-out causes $`Q`$ to rise steeply towards the center. Feedback in this disk is through reflection off the inner cut-out which probably returns some signal in both the short-leading and long-trailing wave channels. It is also possible that some fraction of the in-coming wave is absorbed by the imperfectly shielded ILR, especially when the resonance is broadened by the rapid growth of the mode. The sharper the inner cut-out, the stronger the reflected leading wave becomes and the more vigorous the mode.

Thus, unlike for the soft-centered disks, only a fraction of the in-going trailing wave returns as out-going leading wave to provide feedback to the swing amplifier. The fraction returning as long-trailing waves also shears to short-trailing near corotation (e.g. Mark 1976), but with mild amplification only when $`Q<1.2`$. The effect of this “leakage” from the main swing-amplified feedback loop is to slow the rate of growth of the mode. The response of the disk to noise, on the other hand, is unaffected by feedback and is therefore proportionately more troublesome in the power-law disks.

### 4.3 Amplified particle noise

When the reflected leading signal is weak, as in disks with gentle cut-outs, the predicted mode is easily submerged by particle noise, but particle noise remains an unusually severe problem even in disks with sharp cut-outs. The leading component of the noise spectrum produces apparently growing non-axisymmetric features that obscure our view of the underlying true instability.

When particles are distributed randomly, the level of particle noise depends both on the number of particles per square $`\lambda _{\mathrm{crit}}`$ (e.g. Toomre & Kalnajs 1991) and the vigor of the swing amplifier. The maximum growth that the swing amplifier can produce depends on several properties of the disk: the parameters $`Q`$ and $`X`$, as well as the shear rate $`\mathrm{\Gamma }=1+\beta /2`$. We find that the expected amplification remains roughly constant as $`\beta `$ varies. (We held $`Q_\mathrm{s}`$ constant and the effects of changing both $`\mathrm{\Gamma }`$ and $`X`$ largely cancel.)

The graininess of the disk, on the other hand, does depend on the power-law index, as shown in Figure 3. When all particles have equal mass, the number of particles per square $`\lambda _{\mathrm{crit}}`$ is $`N\xi `$ where $`\xi `$ is the ratio of mass per $`\lambda _{\mathrm{crit}}^2`$ to the total active mass in the simulation. The regions of the disk where $`\xi `$ is small have high $`Q`$ (as shown in the right-hand panel) and do not therefore cause difficulties. Note that the peaks of these curves are of order unity, so we have millions of particles per $`\lambda _{\mathrm{crit}}^2`$ in our largest calculations.

A quiet start is designed to reduce the seed amplitude of non-axisymmetric features in order to be able to observe linear growth over many $`e`$-folds. By placing a number of particles at equal angular spacing on rings and restricting forces to low order sectoral harmonics, initial disturbance forces are many orders of magnitude below those that would have arisen from random initial positions. Such an arrangement is not entirely noise free, however. We provide an initial seed noise spectrum by nudging every particle away from perfect azimuthal spacing by a randomly chosen fraction of $`0.002`$ radians, or about $`0.1^{}`$. This initial noise seeds the growth both of the desired global mode (if one is present) as well as a swing-amplified response. The resulting enhanced disturbance forces further distort the regular particle arrangement, creating enhanced density fluctuations which are further amplified, and so on. Thus the quiet start is itself unstable, and non-axisymmetric amplitudes rise until the distribution of particles is essentially random. The rate at which swing-amplification drives the break up again depends on the responsiveness of the disk as well as $`N\xi `$.

In order to quantify the rate at which quiet starts disrupt, we have run simulations of our usual Mestel disk but with the inner cut-out index $`\nu =2`$. Zang (1976) concluded that the Mestel disk is well stable to modes with $`m1`$ in this case and our addition of an outer taper does not alter his finding. Swing-amplification is still present in this “stable” case, however. The amplitude of $`m=2`$ disturbances in a quiet start simulation grew in a quasi-exponential manner at an approximate rate of 0.06, but we could find no evidence for coherent modes. This simulation confirms that swing-amplification from the initial noise causes a rapid break-up of the quiet start even when the equivalent smooth disk has no instability.

Amplified particle noise is the reason that growth rates appear to increase as $`N`$ decreases in Figure 1(a). Only when $`N>10^6`$ is the break-up rate slow enough to have little effect on the growth rate of the underlying mildly growing normal mode. For fixed $`N`$, the $`\beta =0.25`$ case has the smallest number of particles per $`\lambda _{\mathrm{crit}}^2`$ (dashed curve in Figure 3), and amplified particle noise should be the greatest nuisance for this case, as we have found. It should be noted that quiet starts neither create nor alleviate this difficulty because they reduce the amplitudes of both the mode and the noise in proportion.

### 4.4 Confirmation of bisymmetric stability

The severity of the noise problem is illustrated in Figure 4 which shows snapshots from a noisy start simulation of the linearly stable model; numerical parameters were as given in Table 1 except that we placed each of the 2.5 million particles at a randomly chosen azimuth and included forces from the density sectoral harmonics 0 through 4, with only $`m=1`$ omitted. Despite there being no true global instability, after just 40 dynamical times (or one orbital period at $`R6`$), swing amplified particle noise has already produced a visible spiral pattern of sufficient amplitude to alter the density distribution and heat the inner disk.

Noise-driven activity in this lively $`Q_s=1`$ disk is so strong, even for this large number of particles, that the model later develops a strong bar. The bar does not result from a linear instability, but through saturation of the ILR, or non-linear trapping, as discussed by Sellwood (1989). The swing-amplifier is tamed by raising $`Q_\mathrm{s}`$ and a further long simulation has confirmed the global stability of a similar model with $`Q_s=1.5`$.

## 5 One-armed modes

### 5.1 Comparison with Linear Stability Theory

Linear theory (Zang 1976; Evans & Read 1998b) predicts a dominant pair of $`m=1`$ modes with closely spaced frequencies for cool power-law disks. For our centrally cut-out and outer truncated $`Q_\mathrm{s}=1`$ disks, the predicted frequencies for the dominant mode pair are given in Table 3. These frequency differences are typical. The Table also presents results from a number of separate simulations using different particle samples, but were otherwise identical for each disk and used the numerical parameters given in Table 1.

In all our simulations of the Mestel and power-law disks, the $`m=1`$ data were consistent with just a single growing mode having the indicated eigenfrequency, but we were able to find hints of the two predicted modes in one case. Separation of two such closely-spaced modes would require a long period of integration in the simulations, which cannot be continued indefinitely as the modes saturate after approximately 8 $`e`$-folds, or roughly 100 dynamical times. In the linear regime, therefore, they will have rotated relative to one another by only some 1.6 radians. It seems reasonable that we should find growth rates consistent with that of the more rapidly growing mode, but it is curious that we sometimes find a pattern speed closer to that of the more slowly growing mode. The variations within a set of experiments using the same physical model probably arise from different relative amplitudes of the two modes seeded by the different initial particle noise.

### 5.2 Discussion

The reason for the doublet of $`m=1`$ modes is that there are two possible paths from trailing to leading waves in the inner disk (Toomre 1998, private communication). This is because the reflection of an in-going trailing wave off the inner cut-out returns some short-leading and long-trailing signal. Both paths complete the feedback loop through swing amplification at corotation. As $`Q`$ increases, the distinction between the two paths diminishes, and the splitting of the doublet decreases. When $`Q1.5`$, linear theory recovers a single $`m=1`$ mode.

In order to verify that the discrepancies between predictions and simulations are due to the existence of the closely spaced doublet, we have tested the prediction of an isolated $`m=1`$ mode in a $`Q_s=1.5`$ Mestel disk. Unlike all results reported above, in this case the data from a simulation are well approximated by a single coherent mode, with little interference from noise. Furthermore, we no longer find significant variations in the measured frequency between different draws of particles when $`N=2.5\times 10^6`$. Figure 5 shows that in this clean case we are able to achieve truly convincing agreement between linear theory and simulation – the discrepancy is less than 2% for the largest $`N`$ we have used. This result suggests that in this case almost all the in-going trailing wave reflects to a leading wave, thereby allowing the mode to outgrow any interference from particle noise and to stand out quickly and cleanly. The cleaner reflection at $`m=1`$ is probably largely due to the greater spatial scale of the mode, but the absence of an ILR for $`m=1`$ waves may also be significant, since no residual damping is possible.

## 6 A Completely Stable Disk Model

### 6.1 Linear Stability Theory

Our aim here is to provide a clear-cut example of a stable, massive disk with a nearly exponential radial profile. We adopt a power-law disk with a gently falling rotation curve ($`\beta =0.05`$), but choose the inner cut-out index $`\nu =1`$ and the outer taper index $`\mu =2`$. Figure 6 shows how the pattern speed and the growth rate of the dominant one-armed mode depends on the characteristic radius of the outer taper $`R_\mathrm{c}`$ for a disk of nominal $`Q_\mathrm{s}=1.0`$. \[Note that the nominal $`Q_\mathrm{s}`$ value given in equation (9) refers to the underlying scale-free disk, whereas the actual $`Q`$ value of the tapered disk is higher.\] Figure 6 shows that when the outer taper is centered on $`R_\mathrm{c}<5`$, this disk becomes completely stable to one-armed modes. We have verified that this disk is also linearly stable to the $`m=2,\mathrm{\hspace{0.17em}3}`$ and $`4`$ modes. At least within the framework of linear theory, the strongly tapered disk is completely stable, therefore.

### 6.2 A quasi-exponential disk

We have removed material from the disk with inner and outer tapers without changing the potential. The extra central attraction is provided by unresponsive matter, which we think of as two rigid spherical components, a bulge and halo, whose density distributions are determined as follows.

The disk surface density has a logarithmic spiral decomposition:

$$\mathrm{\Sigma }(R)=\mathrm{\Sigma }_0_{\mathrm{}}^{\mathrm{}}𝑑\alpha B(\alpha )\mathrm{exp}\left[\left(i\alpha \frac{3}{2}\right)\mathrm{log}\left(\frac{R}{R_0}\right)\right],$$
(13)
where $`B(\alpha )`$ is given by

$$B(\alpha )=\frac{1}{2\pi }_0^{\mathrm{}}\frac{\mathrm{\Sigma }(R)}{\mathrm{\Sigma }_0}\mathrm{exp}\left[\left(i\alpha \frac{3}{2}\right)\mathrm{log}\left(\frac{R}{R_0}\right)\right]\frac{dR}{R}$$
(14)
For the power-law disks, the integration is analytic; the details are given in Appendix A. With this decomposition in hand, we find the spherical density $`\rho _\mathrm{d}`$ that reproduces the forces in the equatorial plane provided by the disk is

$$\rho _\mathrm{d}(r)=\frac{\mathrm{\Sigma }_0}{R_0}\left(\frac{r}{R_0}\right)^{{\scriptscriptstyle \frac{5}{2}}}_0^{\mathrm{}}𝑑\alpha (\alpha ^2+\frac{1}{4})K(\alpha ,0)B(\alpha )\mathrm{exp}\left[i\alpha \mathrm{log}\left(\frac{r}{R_0}\right)\right],$$
(15)
where $`K(\alpha ,0)`$ is the axisymmetric Kalnajs (1971) gravity function. Therefore, the spherical bulge density $`\rho _\mathrm{b}`$ and the halo density $`\rho _\mathrm{h}`$ must satisfy

$$\rho _\mathrm{b}(r)+\rho _\mathrm{h}(r)=\frac{v_0^2(1\beta )}{4\pi Gr^{2+\beta }}\rho _\mathrm{d}(r),$$
(16)
with $`\rho _\mathrm{d}`$ on the right hand side given by equation (15). There is some freedom in the partition of the rigid density between the bulge and the halo. We choose the dark halo to be a cored power-law sphere with the potential-density pair (e.g. Evans 1993)

$$\rho _\mathrm{h}(r)=\frac{v_0^2R_0^\beta (r^2(1\beta )+3r_\mathrm{h}^2)}{4\pi G(r^2+r_\mathrm{h}^2)^{2+\beta /2}},\psi _\mathrm{h}(r)=\frac{v_0^2}{\beta }\frac{R_0^\beta }{(r_\mathrm{h}^2+r^2)^{\beta /2}}.$$
(17)
To ensure the halo makes a modest contribution to the rotation curve in the inner galaxy, we chose the core radius to be large, $`r_\mathrm{h}=20`$, giving it an almost uniform density in the inner parts. At large radii, on the other hand, the halo provides almost all of the rotational support. The remainder of the cut-out density is ascribed to a bulge.

The equilibrium properties of this model are shown in Figure 7. Panel (a) shows the rotation curves of the disk, bulge and halo components. The whole model has a slowly falling rotation curve. The active disk has a quasi-exponential profile (panel b). The distribution functions given by equation (10) may be the simplest known which yield quasi-exponential disks with flat rotation curves. As the true $`Q`$ profile of the disk (panel c) is always larger than unity, the disk is axisymmetrically stable. At both small and large radii, $`Q`$ rises sharply where a small fraction of the mass remains active, while $`Q<2`$ at intermediate radii. Panel (d) shows the density of the spherical bulge and halo. As the halo makes a negligible contribution to the mass budget in the inner parts, we can regard this is a “maximum disk” model, albeit with a massive bulge.

It should be noted that although the surface density (Figure 7b) reaches at most barely more than half that of the untapered disk near $`R2`$, the tapered disk’s contribution to the circular speed (Figure 7a) peaks at some 75% of that from the full-mass disk near $`R5`$. This result should not be surprising since the taper has cut away material at larger radii which previously depressed the disk’s contribution to the central attraction.

### 6.3 Simulations

We have used our $`N`$-body code to test the linear theory prediction that this model has no global instabilities. Figure 8 shows the result from a noisy start simulation with a modest number of particles ($`N=5\times 10^4`$) and with $`m=04`$ forces calculated from the particles. All other numerical parameters were as given in Table 1. The force acting on each particle at every step included the fixed central attraction of our bulge/halo in addition to the forces computed from the particles.

The simulation quickly develops some spiral patterns, which result from swing-amplified particle noise, but there is no evidence for a coherent inner bar or any lop-sidedness. A test with ten times more particles behaved similarly, but the spirals were so weak that the disk barely heated at all, and no long-lived coherent waves were detectable at any $`m`$.

We conclude that this model is robustly stable.

## 7 Conclusions

It should no longer be doubted that real galaxies can be stabilized by dense centers, as argued by Toomre (1981). We have finally confirmed Zang’s (1976) global stability analysis for the Mestel disk, and also its generalization to other power laws by Evans & Read (1998b). This encouraging agreement between linear theory and simulations for these idealized galaxy models buttresses Sellwood’s (1985, 1999) claims of global stability for more realistic, but less analytically tractable, models with hard centers.

The lop-sided ($`m=1`$) instabilities of the power-laws disks, first discovered by Zang, can be controlled by reducing the disk surface density to $`<2/3`$ of that of a full-mass disk, without changing the potential (Toomre 1981), much in the same manner as bars can be avoided in soft-centered disks by making them sub-maximal. In this paper, we have been able to show that the $`m=1`$ mode can also be avoided by tapering the massive outer disk until its surface density declines in a quasi-exponential manner. The mechanism is no different since the surface density reaches barely more than half that of the untapered disk at any radius, but the result is a more realistic model in which the disk contribution to the central attraction in the inner parts is decreased to a smaller extent than is the mass.

Many real galaxies do appear to be stabilized by dense centers. High-resolution kinematic observations of spiral galaxies by both Rubin, Kenney & Young (1997) and by Sofue et al. (1999) have revealed high orbital speeds close to their centers. In the cases reported by Sofue et al., for example, the only galaxies with gently rising rotation curves all have rather low luminosity (a peak circular velocity $`<150\text{ km s}^1`$). Such galaxies are required to have large dark matter fractions (e.g. Broeils 1992) – presumably enough to suppress bar-forming instabilities. Almost every galaxy with a circular speed in excess of 150 km s<sup>-1</sup> has a steep inner rise in the rotation curve, and must therefore be bar-stable whatever its dark matter content.

The linear modes of these power-law disks proved to be much harder to identify than those in previous studies. The predicted bisymmetric instabilities are completely obscured by swing-amplified noise unless we employ much larger numbers of particles than are needed in soft-centered disks, and even then the underlying mode does not stand out clearly. We attribute this difference to the different reflection mechanism in the feed-back loop for the power-law disks. The mode stands out clearly from the noise only in cases where the inner reflection is total, as in soft-centered disks, and possibly also for $`m=1`$ modes in hard-centered disks (§5.2). But particle noise readily swamps the weakened modes which result when the reflection involves considerable leakage, as happens for the bisymmetric modes in hard-centered disks (§4.2).

We have been able to confirm the predicted most unstable $`m=1`$ mode in the $`Q_\mathrm{s}=1.5`$ disk with a precision approaching that obtained by Earn & Sellwood (1995) for the bar modes in the isochrone disk. Our simulations of the cooler $`Q_\mathrm{s}=1`$ disk appeared consistent with only a single mode where linear theory predicted a closely-spaced doublet; nevertheless, the frequency we obtained was always some “average” of the predicted pair.

We have shown that one-armed instabilities are tamed by tapering the outer disk. We note that the lop-sided instability has implications for the dark matter candidate proposed by Pfenniger, Combes & Martinet (1994). They suggested that flat rotation curves result from a high surface density of very cold molecular gas in a rotationally supported disk extending beyond the optical edge. This hypothesis is ruled out on stability grounds: massive extended disks with flat rotation curves are grossly unstable to $`m=1`$ modes.

We have presented a realistic galaxy model that is globally stable. It has a quasi-exponential disk in an almost flat rotation curve potential, and includes both a halo with a large core radius and a bulge. The model is stabilized against bar modes by its dense center and against one-armed modes by the exponentially decreasing surface density profile. Even though the halo makes an negligible contribution to the central attraction in the inner parts, the model is clearly globally stable as we demonstrate using both normal mode analysis and $`N`$-body simulations.

We are indebted to Alar Toomre for his comments on an early draft of this paper and detailed, patient explanations of the mode mechanisms. Useful conversations with David Earn are also acknowledged. We thank the referee (Agris Kalnajs) for a thoughtful and stimulating report. This work was begun at the Isaac Newton Institute, Cambridge and continued at the Lorentz Centre, Leiden. The hospitality of both places is acknowledged. JAS is supported by NSF grant AST 96/17088 and NASA LTSA grant NAG 5-6037. NWE thanks the Royal Society for financial support.

## Appendix A The Density Transform

In section 6.2, we exploit the decomposition into logarithmic spirals to find the three-dimensional spherical bulge density that provides the forces to compensate for the matter carved out of our disks. The integration for the logarithmic spiral density transform (14) can be carried through explicitly for the power-law disks. This is done by writing the density as an integral over velocity space of the cut-out distribution function, and then reversing the orders of integrations.

For the case of our quasi-exponential disk, the inner cut-out is $`\nu =1`$ while the outer taper is $`\mu =2`$. The transform then is

$`B(\alpha )`$ $`=`$ $`{\displaystyle \frac{1}{2(2+\beta )}}\left({\displaystyle \frac{2}{\beta }}\right)^{\stackrel{~}{\beta }}{\displaystyle \frac{B(2+1/\beta +\gamma /\beta ,\stackrel{~}{\beta })}{B(1/2+\gamma /2,\stackrel{~}{\beta })}}{\displaystyle \frac{1}{1+\eta ^{2(2+\beta )/(2\beta )}}}`$ (A1)
$`\times \left\{\eta ^{2\stackrel{~}{\beta }+(2+\beta )/(2\beta )}\mathrm{sec}\stackrel{~}{\alpha }+\eta ^{2\stackrel{~}{\beta }}\mathrm{cosec}\stackrel{~}{\alpha }2\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{\hspace{0.17em}2}\stackrel{~}{\alpha }\right\},`$
where $`\eta =r_0v_0/L_\mathrm{c}`$ and

$$\stackrel{~}{\alpha }=\frac{\pi (\frac{1}{2}\beta i\alpha )}{2+\beta },\stackrel{~}{\beta }=\frac{\beta +i\alpha \frac{1}{2}}{2\beta }.$$
(A2)
This is valid for $`\beta >0`$; slightly different expressions hold good for $`\beta 0`$.


text:
# Nonperiodic Flux to Voltage Conversion of an Arithmetic Series Array of dc SQUIDs

\[

## Abstract

A theoretical study on the voltage response function $`V`$ of a series array of dc SQUIDs is presented in which the elementary dc SQUID loops vary in size and, possibly, in orientation. Such series arrays of two-junction SQUIDs possess voltage response functions vs. external magnetic field $`𝐁`$ that differ substantially from those of corresponding regular series arrays with identical loop-areas, while maintaining a large voltage swing as well as a low noise level. Applications include the design of current amplifiers and *quantum interference filters*.

\]

Recently series arrays of dc SQUIDs have been successfully exploited as current amplifiers with wide bandwidth, large dynamic range and low noise level . Using thin-film Nb-technology, amplifiers consisting out of up to $`10^3`$ identical dc SQUID loops have been fabricated. Such serial devices are characterized by large voltage swings of several $`mV`$ and current-to-voltage transfer functions of some $`V/mA`$ so that a direct connection to a room temperature preamplifier is feasible.

At present only series arrays consisting out of *identical* dc SQUID loops have been described in the literature . The dc voltage response function $`V`$ of such regular series arrays displays a $`\mathrm{\Phi }_0`$-periodicity just as a single dc SQUID where $`\mathrm{\Phi }_0=\frac{h}{2e}`$ is the elementary flux quantum . However, this periodicity may represent a serious limitation for the modes of operation of such devices. To get a linear flux-to-voltage conversion a feedback circuit is used in most applications in which the dc SQUID acts as a null-flux detector. Furthermore special electronic devices and efficient background shielding is often required . Here theoretical studies on the voltage response of arithmetic and irregular series arrays of dc SQUID loops are presented, in which the individual loop-areas are *not* all equal. For such devices the advantages of dc SQUID series arrays are preserved but limitations due to the periodicity of the voltage response can be circumvented.

The arrays under consideration consist of $`N`$ two junction SQUID loops connected in series. The bias current $`I_b`$ is fed into the array as indicated in Fig.(1). In general the areas of the $`N`$ loops in a *generic* array differ in size and, possibly, in orientation. Let $`𝐚_n`$ be the orientated area element of the $`n^{th}`$ dc SQUID loop. The magnetic flux threading $`𝐚_n`$ is then $`\mathrm{\Phi }_n=𝐁,𝐚_n`$, where $`𝐁`$ is the total magnetic field. Taking into account all inductive effects in the dc SQUID array the total magnetic field $`𝐁=𝐁^{(p)}+𝐁^{(s)}+𝐁^{(c)}`$ is a superposition of the primary magnetic field $`𝐁^{(p)}`$ one wants to measure, the secondary magnetic field $`𝐁^{(s)}`$ induced by currents that flow in the dc SQUID array and, possibly, the magnetic compensation field $`𝐁^{(c)}`$ induced by compensation current(s) $`I_{comp}`$ flowing through a set of suitable orientated compensation coils or wires.

As a special case a generic *planar* dc SQUID array is shown in Fig.(1). The primary magnetic flux is coupled into the individual dc SQUID loops by the signal or input current $`I_{inp}`$ which flows through a common input flux bias line. In typical applications $`I_{inp}`$ is the current provided by a pick-up loop or by other signal sources. In addition the flux in the individual dc SQUID loops may be modulated by the current $`I_{comp}`$ flowing through a common compensation or control flux bias line.

Let $`M_n^{(c)}`$ and $`M_n^{(p)}`$ denote mutual inductances of the currents $`I_{comp}`$ and $`I_{inp}`$ associated with the $`n^{th}`$ loop with area-element $`𝐚_n`$ and let $`𝐋_A`$ be the inductance matrix describing all other inductive couplings in the circuit . The total magnetic flux in the $`n^{th}`$ loop is then given by $`\mathrm{\Phi }_n=M_n^{(p)}I_{inp}+M_n^{(c)}I_{comp}+𝐋_A𝐈`$, where $`𝐈=(I_1,\mathrm{},I_{2N})`$ is the vector of the $`2N`$ currents $`I_n`$ through the Josephson junctions. Since $`I_{comp}`$ should be able to compensate the primary induced magnetic fluxes $`M_n^{(p)}I_{inp}`$ in all loops simultaneously the relation $`M_n^{(c)}=\alpha M_n^{(p)}`$ should hold for all $`n`$, where $`\alpha `$ is a constant.

Within the range of validity of the RCSJ-model  the current $`I_m`$ through the $`m^{th}`$ array junction with gauge invariant phase difference $`\phi _m`$, critical current $`I_{c,m}`$, normal resistance $`R_m`$, and capacitance $`C_m`$ is given by

$$I_m=I_{c,m}\mathrm{sin}(\phi _m)+\frac{\mathrm{}}{2eR_m}_t\phi _m+\frac{\mathrm{}C_m}{2e}_t^2\phi _m.$$
(1)
For simplicity identical parameters for the junctions are supposed in the following discussion: $`I_{c,m}=I_c`$, $`R_m=R`$ and $`C_m=C`$. The dimensionless parameter $`\beta _C=\frac{2\pi I_cR^2C}{\mathrm{\Phi }_0}`$ characterizes then the capacitance effect of the junctions within this model .

In the overdamped regime, $`C=0`$, and for vanishing inductive coupling in the array, $`𝐋_A=0`$, an analytical expression for the dc voltage response function $`V`$ of the dc SQUID series can be determined. Under conditions where a constant current $`I_b>2I_c`$ is biased, and assuming for simplicity a static magnetic field $`𝐁=𝐁^{(p)}+𝐁^{(c)}`$, the output voltage $`V_n`$ of the $`n^{th}`$ dc SQUID in the array as a function of the bias current $`I_b`$ and of the applied flux $`\mathrm{\Phi }_n`$ is given by

$$V_n=I_cR\sqrt{\left(\frac{I_b}{2I_c}\right)^2\left|\mathrm{cos}(\pi \frac{\mathrm{\Phi }_n}{\mathrm{\Phi }_0})\right|^2}.$$
(2)
The voltage response $`V`$ of a series array containing $`N`$ dc SQUIDs is then given as the sum of the voltages $`V_n`$:

$$V=\underset{n=1}{\overset{N}{}}V_n.$$
(3)
Since the maximum voltage swing increases linearly with the number $`N`$ of dc SQUIDs in the array the gain $`\frac{V}{I_{inp}}`$ is proportional to $`N`$. However, the noise in the array increases only as $`N^{\frac{1}{2}}`$, because the noise voltages implied by the junctions in the individual SQUID-loops of the array are expected to add incoherently . These features make generic series arrays attractive for amplifiers and magnetometers.

Consider, as a special case, a regular array, consisting out of $`N`$ *identical* dc SQUID loops so that $`𝐚_n=𝐚_N`$ with identical mutual inductances $`M_n^{(c)}=M_N^{(c)}`$ and $`M_n^{(p)}=M_N^{(p)}`$ of the flux bias lines. As a result the flux threading each individual loop is then given by $`\mathrm{\Phi }_n=\mathrm{\Phi }_N`$. In this case the voltage response, $`V=NV_N`$, displays as a function of flux a $`\mathrm{\Phi }_0`$-periodicity.

A voltage response function $`V`$ with a much longer period may be obtained in an *arithmetic* dc SQUID series array where the orientated area elements increase in size according to the arithmetic relation

$$𝐚_n=\frac{n}{N}𝐚_N,n=1,\mathrm{},N.$$
(4)
Provided the mutual inductances are chosen such that $`M_n^{(c)}=\frac{\left|𝐚_n\right|}{\left|𝐚_N\right|}M_N^{(c)}`$ and $`M_n^{(p)}=\frac{\left|𝐚_n\right|}{\left|𝐚_N\right|}M_N^{(p)}`$ the fluxes are determined by $`\mathrm{\Phi }_n=\frac{n}{N}\mathrm{\Phi }_N`$, where $`\mathrm{\Phi }_N`$ is the flux through the largest area element $`𝐚_N`$. In this case the individual dc SQUIDs operate with different periodicity in flux. This implies that the voltage response $`V`$ of the whole arithmetic array as a function of $`\mathrm{\Phi }_N`$ is periodic with period $`N\mathrm{\Phi }_0`$. This period is $`N`$ times larger than for a regular array, consisting of $`N`$ identical dc SQUID loops with size $`𝐚_N`$.

In Fig.(2) the voltage response for an arithmetic series array with $`N=100`$ dc SQUIDs is shown. The bias current $`I_b`$ was adjusted for maximum voltage swing slightly above $`2I_c`$. In the depicted interval the voltage response shows pronounced antipeaks at $`\mathrm{\Phi }_N=0`$ and at $`\mathrm{\Phi }_N=N\mathrm{\Phi }_0`$. Between these two peaks the output voltage oscillates with a very small amplitude around an average voltage drop $`\overline{V}`$.

For a sufficiently high number of elements in the arithmetic array the flux-distribution over the dc SQUID loops is uniform between some $`\mathrm{\Phi }_{min}`$ and $`\mathrm{\Phi }_{max}`$. Let $`\kappa =\frac{\mathrm{\Phi }_{min}}{\mathrm{\Phi }_{max}}`$ be the flux spread coefficient, $`f=\pi \frac{\mathrm{\Phi }_{max}}{\mathrm{\Phi }_0}`$ the frustration of the greatest loop in the array and let $`\omega _0=\sqrt{\left(\frac{I_b}{2I_c}\right)^21}`$ be the characteristic frequency of a single dc SQUID for zero flux. If the frustration $`f`$ is limited to the periodicity interval $`\left|f\right|\pi \frac{N}{2}`$ the voltage response $`V`$ as a function of $`f`$ can be approximated by an integral:

$`V`$ $``$ $`{\displaystyle \frac{NI_cR}{\mathrm{\Phi }_{max}\mathrm{\Phi }_{min}}}{\displaystyle \underset{\mathrm{\Phi }_{min}}{\overset{\mathrm{\Phi }_{max}}{}}}\sqrt{\left({\displaystyle \frac{I_b}{2I_c}}\right)^2\left|\mathrm{cos}(\pi {\displaystyle \frac{\mathrm{\Phi }}{\mathrm{\Phi }_0}})\right|^2}𝑑\mathrm{\Phi }`$ (5)
$`=`$ $`NI_cR{\displaystyle \frac{\omega _0}{(1\kappa )}}{\displaystyle \frac{1}{f}}\left[E(f,m)E(\kappa f,m)\right],`$ (6)
where $`m=\omega _0^2`$ and $`E(f,m)`$ denotes the elliptic integral of the second kind . The voltage versus flux relation displayed in Fig.(2) shows an excellent agreement between the theoretical result Eq.(3) and its analytical approximation Eq.(5). For bias currents $`I_b=2.2I_c`$ there results an average voltage drop of $`\overline{V}=0.81I_cRN`$.

At first sight all the above results apply only in the overdamped junction regime and only for vanishing inductive couplings in the arithmetic dc SQUID array, i.e. for $`\beta _C=0`$ and $`𝐋_A=0`$. A dimensionless measure for the inductive effects in a single dc SQUID loop with loop-inductance $`L`$ is $`\beta _L=\frac{LI_c}{\mathrm{\Phi }_0}`$. Dependent on $`\beta _L`$ there exists an optimal size for any dc SQUID loop  which should coincide with the maximal loop size $`𝐚_{max}`$ in a series array. The corresponding maximal $`\beta _{L,max}`$ is then a measure of the self- and mutual-inductive couplings among the currents flowing in the dc SQUID array.

Taking into account all inductive couplings $`𝐋_A`$ in the arithmetic array our computer simulations of the full non-linear dynamics  of the $`2N`$ coupled Josephson phases $`\phi _i`$ reveal that expression Eq.(5) compares qualitatively well with the calculated voltage response functions. In particular the long periodicity of the voltage $`V`$ vs. flux $`\mathrm{\Phi }_N`$ relation is not affected by inductive effects. Qualitatively similar behavior was found for a parallel multi-junction interferometer .

As far as the irregularity is concerned, we also find that in an arithmetic array $`V`$ is very responsive to adding small random fluctuations to the size distributions of the area elements. In this case $`V`$ becomes *nonperiodic* with a pronounced antipeak only around $`\mathrm{\Phi }_N=0`$. For fixed $`\beta _C`$, but different strengths of the magnetic coupling $`\beta _{L,max}`$, Fig.(3) shows the voltage response for an arithmetic planar series array with $`N=100`$ dc SQUIDs where small random fluctuations were added to the loop size distribution Eq.(4).

For fixed $`\beta _C`$ the global minimum of $`V`$ at $`\mathrm{\Phi }_N=0`$ depends only slightly on the strength $`\beta _{L,max}`$ of the magnetic coupling whereas the voltage branch $`V\overline{V}`$ depends strongly on $`\beta _{L,max}`$, as can be seen in the inset of Fig.(3). For increasing $`\beta _{L,max}`$ the voltage swing and the maximum of the voltage transfer function $`\frac{V}{\mathrm{\Phi }_N}`$ decrease as it is the case for conventional single dc SQUIDs when $`\beta _L`$ increases . If on the other hand $`\beta _{L,max}`$ is kept constant the global minimum of $`V`$ at $`\mathrm{\Phi }_N=0`$ is an increasing function of $`\beta _C`$ as a comparison of Fig.(2) and Fig.(3) reveals. In addition Fig.(3a) shows an enhancement of the branch $`V\overline{V}`$ for $`\beta _C=0.5`$ in comparison to the underdamped case, since for single dc SQUIDs the effect of the capacitance $`\beta _C`$ is to increase the voltage near integer values of the applied flux . Although the total area $`A_{tot}=\frac{N+1}{2}\left|𝐚_N\right|`$ of an arithmetic array is smaller by a factor of $`2`$ compared to a regular array with $`A_{tot}=N\left|𝐚_N\right|`$, this effect provides that the voltage transfer function of both arrays is comparable. This suggests a higher integration density for arithmetic array circuits on chip.

In summary, arithmetic series arrays of two-junction SQUIDs possess voltage response functions with a much longer period with respect to the applied magnetic field than regular series arrays, while maintaining a comparable transfer function and a low noise level. In particular if the loop-sizes or loop-orientations are distributed randomly the voltage response becomes nonperiodic with a single pronounced antipeak only around zero applied magnetic field. These features are preserved when all capacitive and inductive effects are taken into account. Therefore arithmetic or irregular series arrays of dc SQUIDs can be used as *quantum interference filters* for various applications, including, e.g., the relatively simple and extremely sensitive measurement of the *absolute* strength of magnetic fields.

Acknowledgments: We thank R. P. Huebener, R. Kleiner and T. Träuble for useful discussions. Support by ”Forschungsschwerpunktprogramm des Landes Baden-Württemberg” is gratefully acknowledged.


text:
# Mean field theory for global binding systematics

## I Introduction

An important goal of nuclear theory is to predict nuclear binding energies. While mean field theory offers the most fundamental basis to understand nuclear structure, paradoxically it has not been as successful than other approaches in making a global fit to nuclear binding energies. The most accurate theory of nuclear binding systematics  starts from the liquid drop model, and treats shell effects perturbatively. It fits the binding energies with an RMS deviation of 0.67 MeV, taking 15 free parameters and a similar number of fixed parameters to achieve the fit. No such systematic study with so many parameters has been attempted in a purely microscopic approach, but there are a number of partial studies beginning with the pioneering work of Vautherin and Brink using an energy functional based on the Skyrme interaction. Noteworthy recent papers using Skyrme interactions are by Patyk et al.  and by Brown . Patyk et al. find that Skyrme interactions taken from the literature give RMS errors greater than 2 MeV. This level is also found for the Gogny interaction, which unlike the Skyrme is finite range. Brown recently made a new Skyrme fit to closed shell nuclei, including radii and spectroscopic properties in the fitting . He found an RMS deviation of 0.8 MeV for the 10 nuclei he considered, encouraging the hope that deviations below 1 MeV might be reached microscopically. Of course, open shell nuclei have significant correlation energies that must be included. We discuss how this might be done in Sect. 3 below. The other problem is the choice of energy density functional, and there may be reason to use other forms than the Skyrme or finite-range generalizations. This is discussed in the next section.

## II New forms for the energy functional

Some perspective on the energy functional can be obtained from the analogous problem in condensed matter physics. Correlation effects are rather mild in the many-electron problem, and the mean field approach is very successful. The energy functional analogous to the Skyrme is called the local density approximation (LDA), and it is widely used to calculate structures of many-atom systems. Its accuracy for chemical purposes is inadequate, however. For example, in a comparison of different functionals Perdew et al.  noted that the LDA had a mean absolute error of 1.4 eV in a sample of small molecules with binding energies in the range 5-20 eV*** It should be mentioned, however, that the LDA energy functional is constructed ab initio without fitting binding data.. Two refinements of the energy functional, going beyond the LDA, make a dramatic improvement in the quality of agreement. The first refinement is to include a term in the energy functional that depends on the gradient of the density. Gradient terms are already present in the Skyrme interaction, but in the electron system improvement only appears with a nonlinear functional form of the gradient term,

$$e(n,n)\frac{f(n)}{1+a(n/n)^2}$$
(1)
The mean absolute error decreases by a factor of 4, to 0.35 eV, when this term is addedAgain, no free parameters are added with this term; the functional form and dependencies on $`n`$ and $`\tau `$ are constructed to simulate the nonlocality of the exchange.. The other improvement is to allow the functional to depend on the kinetic energy density $`\tau `$ and well as on the local density $`n(r)`$. This also is a familiar feature of the Skyrme interaction, but for the electron case the $`\tau `$ is combined with the $`n`$ dependence in eq. (1). The resulting mean absolute error is reduced by almost a factor of three to give a final mean absolute error of 0.13 eV.

The nuclear problem is different from the electronic is one important respect. In the electronic problem, much of the loss of accuracy is due to the exchange potential, which is intrinsically nonlocal but must be treated in a nearly local approximation for computational reasons. In contrast, in the nuclear problem the strong interaction is short-ranged, implying that the exchange is short-range as well and thus suited to local approximations. Nevertheless, it might be that more complicated functional forms such as eq. (1) could be useful in the nuclear problem. Indeed, this kind of generalization has be examined by Fayans . He used such terms in an energy functional that was fit to 100 spherical nuclei. He obtained an RMS binding error of 1.2 MeV, a factor of two better than the Gogny or the published Skyrme functionals.

## III Correlation energy

As mentioned earlier, in open shell nuclei the correlation energy associated with nearly degenerate configurations can be of the order of several MeV, so the single-configuration mean field approximation is not accurate enough for global energy systematics.

We believe that the following correlations should be considered explicitly in the theory: $``$center-of-mass delocalization
$``$quadrupole deformations
$``$pairing.
These are the obvious correlations associated with symmetry breaking in the mean field Hamiltonian. Translational invariance is always broken in finite systems and rotational invariance is often broken as well. Pairing treated by BCS theory violates particle number conservation.

How important are these correlations to the energetics? For the center of mass, the correlation energy can be estimated in the harmonic oscillator model (but see below) as $`3\mathrm{}\omega _0/430/A^{1/3}`$ MeV. The resulting magnitude of several MeV is certainly much larger than the allowable error in a global mass theory. However, the energy varies very smoothly and one could question whether it needs to be treated explicitly as a correlation energy or whether its effect can be subsumed into the parameters of the mean-field energy functional. This might not be the case because the mean field functional determines most directly the leading terms in the liquid drop expansion, varying as $`A`$ and $`A^{1/3}`$. One should not ask it to simulate a completely different $`A`$-dependence. This argument, and some mean-field calculations to support, was given recently by Bender et al. .

The correlation energy associated with the deformations may be thought of as the energy gained by projecting the states of good angular momentum out of the deformed intrinsic state which contains many angular momenta. Its magnitude can be quite large on our accuracy scale of hundreds of keV. For example, the nucleus 20Ne has a Hartree-Fock ground state close in structure to the  SU(3) state of the harmonic oscillator. Combining probabilities of the different angular states with the energies of the states derived from the experimental spectrum, one can easily show that the energy gained by the projection is of the order of 4 MeV. A similar number was obtained recently in the projected mean-field calculations of 24Mg by ref. . Thus this correlation energy should be included to accuracy of 10% or better to achieve the desired accuracy of the mass theory.

The situation with pairing is less severe. It is certainly necessary to including pairing in a theory of masses, just to get realistic odd-even mass differences, but the need to include correlation effects beyond the Hartree-Fock-Bogoliubov theory is less clear. With certain simplified pairing interactions the pairing Hamiltonian can be solved exactly, without making the BCS approximation . The error in the BCS energy due to the number nonconservation is of the order of 0.5 MeV , Table 11.1. This might be significant in the global mass theory and we shall include it in our discussion.

There are many ways that correlation energies can be calculated. Most popular seems to be the obvious method in which the eigenstates of the symmetry are projected out of the mean-field wave function. If the energy minimization is carried out after the projection, this method is rather costly to use and probably not suitable for a global mass theory. For a global theory it is important that the method be simple computationally and also that it be systematic, applicable in principle to all possible mean field solutions. Particularly important is that it does not introduce a discontinuity when the mean field solution changes its character. In the systematic development of many-particle perturbation theory, the leading term beyond the mean-field approximation gave the correlation energy as an integral over the RPA excitation modes. In a finite system, the RPA correlation formula is given by

$$E_{corr}=\frac{1}{2}\left(\underset{i}{}\mathrm{}\omega _iTr(A)\right),$$
(2)
where $`\omega _i`$ is the (positive) frequency of the RPA phonon and $`A`$ is the $`A`$ matrix in the RPA equations. This approach was first proposed by Friedrich and Reinhard . It seems to us that this formula is well-suited to the requirements we need for a systematic mass theory. We shall first argue that the formula is adequate in principle, and then take up the issue of computational feasibility. In the next section we summarize experience with simple model Hamiltonians that show that eq. (2) is more accurate than commonly used projection techniques, or is easier to calculate. In the section following we examine specific algorithms for calculating nuclear binding energies efficiently.

## IV Experience from simple models

### A Center-of-mass localization

The Hamiltonian of two particles interacting through a quadratic potential $`V(r_{12})=m\omega _0^2r_{12}^2/2`$ is solvable exactly and also has an analytic mean-field approximate solution. One might guess that the RPA correlation energy might give the correction exactly, because one can derive the RPA by considering quadratic approximations in a path integral formulation of the problem. This is indeed the case. The RPA spectrum has two states, a zero-frequency mode and a finite frequency mode. Putting these frequencies in eq. (2), one finds that the correction to the mean-field energy $`\mathrm{}\omega _0`$ is just what is needed to give the exact energy for the total, $`\mathrm{}\omega _0/\sqrt{2}`$.

It is interesting to compare the RPA approach with other ways of dealing with correlation energies associated with broken symmetries. In the case of center-of-mass motion, a recipe that is often used is to subtract the expectation value of the center of mass operator from the mean field energy (e.g. in ref. ). With our Hamiltonian, this prescription gives

$$E_{cm}=MF\left|\frac{1}{2}\frac{(p_1+p_2)^2}{2m}\right|MF=\frac{1}{4}\mathrm{}\omega _0.$$
(3)
The total $`E_{MF}+E_{cm}=3\mathrm{}\omega _0/4`$ is not exact, although it is close to the exact energy, $`\mathrm{}\omega _0/\sqrt{2}`$. This study shows that for this first kind of correlation the RPA formula provides a better method to calculate the associated energy.

### B Deformations

When the mean field solution is deformed, a continuous symmetry is broken. As in the above example, a signature of the broken continuous symmetry is a zero-frequency RPA mode. A model to test theories of the correlation energy should thus have a corresponding continuous symmetry. We constructed a model with those properties in ref. , making a generalization of the Lipkin model. In the original Lipkin model that has been studied for 40 years, one considers many distinguishable particles each of which can be in one of two states. For that model, the ground-state correlation was studied by Parikh and Rowe . They compared various methods of treating the correlation energy, finding that the RPA formula worked best. In ref. , we extended the Lipkin model to a three-state wave function to get sufficient degrees of freedom for a continuously broken symmetry. The symmetry is imposed on the Hamiltonian by requiring it to be invariant under transformations within two of the three states. The two degenerate upper states could be thought of as the first excited states of a two-dimensional harmonic oscillator, thus allowing deformed wave functions in two dimensions. As expected, when the mean field solutions and its RPA excitations are calculated, ones sees a zero-frequency when the mean field solution is deformed. There is also another mode at finite frequency. In this case, the first mode corresponds to rotational motion perpendicular to the symmetry axis, while the second mode corresponds to the beta vibration. In the “spherical” phase, the RPA frequencies for the modes are identical and have finite frequency. Figure 1 shows the RPA frequencies as a function of $`\chi V(N1)/ϵ`$ for the particle number $`N=20`$, where $`V`$ and $`ϵ`$ are the strength of the interaction and the single-particle energy, respectively. One can clearly see the discontinuity at the critical point $`\chi =1`$.

Figure 2 compares the ground state energy as a function of $`\chi `$ obtained by several methods. The number of particle $`N`$ is set to be 20. The solid line is the exact solution obtained by numerically diagonalizing the Hamiltonian. The dashed line is the ground state energy in the Hartree-Fock approximation. It considerably deviates from the exact solution through the entire range of $`\chi `$ shown in the figure. The dot-dashed line takes into account the RPA correlation energy in addition to the HF energy. Clearly the RPA significantly improves the results. The corresponding energy as a function of the number of particles is shown in Fig. 3. One sees that the RPA correlation energy is reliable for large $`N`$, but may be problematic when there are only a few valence particles. That situation is further complicated by the pairing interaction, which may dominate the mean field solution for $`N=2`$.

### C Pairing

Pairing is the final example of a long-range correlation that significantly affects the energy. The mean field approach leads to the BCS theory, whose ground state has indefinite particle number. An early study by Bang and Krumlinde  showed that the RPA formula reproduces the exact correlation energy rather well in a schematic model. The RPA method has in fact been used in realistic models of deformed nuclei. The RPA correlation in the normal phase was studied in ref.  using the self-consistent version of RPA.

Kyotoku et al.  derived the analytic solution for a model first proposed in ref. , fermions in a space of two nondegenerate j-shells interacting with a pairing Hamiltonian. They were able to solve the model exactly, and then compare the energy with several approximate methods to calculate the correlation corrections to the BCS energy. They found that the ground state energy in the BCS + quasi-particle RPA (QRPA) coincides with the exact solution at the leading order of an expansion in 1/$`N`$, the number of particles in the system. None of the other methods obtained the correct coefficient of the leading order contribution. For example, the well-known method of Lipkin and Nogami  gave a result that is only correct in the limit of a strong pairing force (see Table I in ref. ).

In ref. , we specifically compared the RPA with the computationally attractive alternative methods, testing the behavior across shell closures and considering both even- and odd-$`N`$ systems. Taking as a test case the two-level problem with a level degeneracy of $`\mathrm{\Omega }`$ = 8 and a Fermi energy half way between the levels, the RPA was much superior to the Lipkin-Nogami method over most of the range of pairing strengths. This behavior is consistent with the result of ref. as we discussed above. However, right at the phase transition point the two methods had comparable errors of opposite sign.

We therefore also looked at a more realistic situation, varying the particle number $`N`$ rather than the interaction strength $`G`$. We consider the pairing energy in oxygen isotopes, taking the neutron 1p and 2s-1d shells as the lower and higher levels of the two-level model. The pair degeneracy $`\mathrm{\Omega }`$ thus reads $`\mathrm{\Omega }_1=3`$ and $`\mathrm{\Omega }_2=6`$, for the lower and the upper levels, respectively. The number of particle in a system is given by $`N=A82`$ for the AO nucleus. We assume that the energy difference between the two levels $`ϵ`$ is given by $`ϵ=41A^{1/3}`$ and the pairing strength $`G=23/A`$. The upper panel of Fig. 4 shows the ground state energy as a function of $`A`$. In order to match with the experimental data for the 16O nucleus, we have added a constant $`72.8`$ MeV to the Hamiltonian for all the isotopes. The exact solutions are denoted by the filled circles. The deviation from the BCS approximation (the dashed line) is around 2 MeV for even $`A`$ systems and it is around 1.2 MeV for odd $`A`$ systems. This value varies within about 0.5 MeV along the isotopes and shows relatively strong $`A`$ dependence. One can notice that the RPA approach (the dot-dashed line) reproduces quite well the exact solutions. In contrast, the Lipkin-Nogami approach (the thin solid line) is much less satisfactory and shows a different $`A`$ dependence from the exact results. The pairing gap $`\mathrm{\Delta }`$ in the BCS approximation and in the Lipkin-Nogami method is shown separately in the lower panel of Fig. 4. For the Lipkin-Nogami method, we show $`\mathrm{\Delta }+\lambda _2`$, which are to be compared with experimental data . The closed shell nucleus 16O and its neighbor nuclei 15,17O have a zero pairing gap in the BCS approximation, and the Lipkin-Nogami method does not work well for these nuclei. On the other hand, the RPA approach reproduces the correct $`A`$ dependence of the binding energy. Evidently, the RPA formula provides a better method to compute correlation energies than the Lipkin-Nogami method, especially for shell closures.

## V Is RPA computationally feasible?

We now discuss the practicality of using eq. (2) to calculate the correlation energy. In general, the RPA is computationally more demanding than the mean field minimization for the ground state by an order of magnitude or more. One must diagonalize an RPA matrix whose dimension is $`2D\times 2D`$, where $`D`$ is the number of particle-hole configurations. This number can be huge if one is interested in deformed nuclei or heavy spherical nuclei. A widely used way around this is to take a residual interaction has a separable form. Then the matrix equation to be solved has the dimension of the rank of the separable interaction; with a single term it is just the well-known algebraic dispersion relation.

Given a separable interaction, there are several efficient ways to the RPA correlation energy (2) without explicitly calculating all the roots of the dispersion relation. One method was recently proposed by Dönau et al. and also by Shimizu et al. . Instead of directly calculating the RPA correlation energy according to eq. (2), one carries out an integration of a function which depends on a free response function and its first derivative in a complex energy plane. An advantage of this method is that one can choose the integration path so that the integrand is smooth enough and thus the mesh of the numerical integration along such path can be much larger than the actual energy intervals of RPA solutions $`\omega _\alpha `$. This method is useful particularly when a separable interaction is used so that the free response function and its first derivative are analytically evaluated.

Alternatively, one can also use the Lanczos-type method proposed in ref.  to evaluate the RPA correlation energy. As we show below, this method quickly converges when the interaction is separable. The idea is to define at the outset a characteristic operator associated with each kind of correlation. That operator applied to the mean-field ground state gives an excited state, which is taken as the the first vector in a space built up by applications of the Hamiltonian to existing states. Eq. (2) is then evaluated in the restricted spaces, and the method would be computational feasible if the convergence is rapid enough.

Suppose $`A`$ and $`B`$ matrices in the RPA equation are given by

$`A_{ij}`$ $`=`$ $`ϵ_i\delta _{i,j}+\kappa f_if_j,`$ (4)
$`B_{ij}`$ $`=`$ $`\kappa f_if_j,`$ (5)
where $`i`$ and $`j`$ label particle-hole configurations and $`f_i`$ is normalized as $`_if_i^2=1`$. For such an interaction, the collective operator can be chosen as $`\psi _i=f_i`$. Notice that $`f`$ is an eigen-vector of the matrix $`B`$ with eigen-value $`\kappa `$ and also that $`B\varphi =0`$ for any vector $`\varphi `$ which is orthogonal to $`f`$. Starting from the initial vectors $`|X_1=|\psi `$ and $`|Y_1=0`$, the Lanczos manipulation  transforms the $`A`$ and $`B`$ matrices into the form of

$$𝑨^{}=\left(\begin{array}{cccc}e_1& a_1& 0& \\ a_1& e_2& a_2& \\ 0& a_2& e_3& \\ & & & \mathrm{}\end{array}\right),𝑩^{}=\left(\begin{array}{cccc}\kappa & 0& 0& \\ 0& 0& 0& \\ 0& 0& 0& \\ & & & \mathrm{}\end{array}\right).$$
(6)
The Lanczos basis for backward scattering amplitude $`|Y_i`$ remain all zero in the manipulation, and thus the transformed $`B^{}`$ matrix has only one element which is not equal to zero. This is not in general the case for non-separable interaction. Note that the $`B^{}`$ matrix measures the degree of correlation in the transformed space. Since it has a very simple form, the correlation energy evaluated in the restricted space converges very rapidly.

We have tested the method with RPA matrices given by eq. (5) with $`ϵ_i=e_0+(i1)\mathrm{\Delta }e`$. The method works very well when there is a gap in the particle-hole spectrum. For example, taking a model with $`D`$=20 particle-hole states and other parameters $`f_i=1`$, $`\mathrm{\Delta }e=e_0/D`$, and $`\kappa =0.035e_0`$, the RPA spectrum has a collective state at $`\omega _1=0.145e_0`$. The correlation energy from eq. (2) with all eigenvalues, is $`E_{corr}=0.542e_0`$. The single-state approximation starting from the state $`\psi _i=f_i/\sqrt{D}`$ is only off by 20%. The two-state approximation has an accuracy of 1%. The calculational effort to get this numbers is essentially that of three applications of the Hamiltonian to the ground-state single-particle wave function, much less than the effort to get the converged wave functions in the Hartree-Fock minimization. See Table I.

From our point of view, the problem is then to define a reliable separable interaction that can be used globally in RPA calculations. For particle-hole residual interactions, there have been proposed a few possible ways based on the self-consistency argument to construct a separable interaction . The basic argument is that the collective motions primarily result in the displacement of the nuclear surface; thus the transition potential is can be generated by displacing the self-consistency potential associated with the ground state density. In equations, the transition density has a form locally given by the gradient of the static density, and the transition potential is the corresponding gradient of the static potential:

$`\delta \rho `$ $`=`$ $`a(\mathrm{\Omega }){\displaystyle \frac{d\rho _0}{dr}}`$ (8)
$`\delta V=a(\mathrm{\Omega }){\displaystyle \frac{dU_0}{dr}},`$
where $`\mathrm{\Omega }`$ is an angle given the direction to an element of the nuclear surface. The separable residual interaction then has the shape of eq. (8) and the required magnitude to satisfied eq. (8). Thus

$$v(r_1,r_2)=\delta V(\mathrm{\Omega }_1,r_1)\delta V(\mathrm{\Omega }_2,r_2)/𝑑r\delta V\frac{d\rho }{dr}.$$
(9)
It is amusing to compare the above self-consistent definition with the microscopic particle-hole interaction proposed by Migdal . For the above transition density, his transition potential would be given by

$$\delta V(r)=V_0(f^{ex}+(f^{in}f^{ex})\rho _0(r))\delta \rho (r).$$
(10)
The two transition potentials for 208Pb are compared in Fig. 5. We generate the static potential $`U_0`$ using the velocity-independent part of the Skyrme interaction with SIII parametrisation and we use a Fermi distribution for a static density $`\rho _0(r)`$. The parameters of the Migdal interaction are given in ref. . We see both the transition potentials have a similar surface-peaked radial dependence.

The above self-consistency argument can be applied very easily to the translation mode, where $`a(\mathrm{\Omega })=\mathrm{cos}\theta `$ for translations in the $`z`$-direction. For the rotational degrees of freedom, we would use the displacement fields associated with the 5 components of the quadrupole operator, i.e. $`a_\mu (\mathrm{\Omega })=Y_{2\mu }(\mathrm{\Omega });\mu =2,1,0,1,2`$. Note that the RPA correlations are to be evaluated for spherical as well as deformed nuclei; a major point of the approach is that the theory applies to all cases.

It is not so simple to construct a global, separable pairing interaction. The commonly used forms for the pairing interaction give divergent results without a cutoff in space of states included in the BCS calculation. Here we propose a non-local surface-peaked separable form,

$$v(𝒓_1𝒓_2,𝒓_1^{}𝒓_2^{})=\delta (𝒓_1𝒓_2)\delta (𝒓_1^{}𝒓_2^{})\left[\frac{\chi }{2}\frac{dU_0}{dr_1}\frac{dU_0}{dr_1^{}}\underset{\lambda \mu }{}Y_{\lambda \mu }(\widehat{𝒓}_1)Y_{\lambda \mu }^{}(\widehat{𝒓}_1^{})\right].$$
(11)
The pairing matrix element of this interaction are to be evaluated as

$$p_1\overline{p}_1\left|v\right|p_2\overline{p}_2=\frac{\chi }{2}\underset{\lambda \mu }{}p_1\left|\frac{dU_0}{dr}Y_{\lambda \mu }(\widehat{𝒓})\right|p_1p_2\left|\frac{dU_0}{dr^{}}Y_{\lambda \mu }^{}(\widehat{𝒓^{}})\right|p_2$$
(12)
This form is inspired by recent observations that the pairing is essentially a surface phenomenon . A similar surface-peaked separable interaction was used for the particle-hole channel in ref. . We have tested the surface-peaked separable pairing interactions (11) by comparing it with the density-dependent delta interaction proposed in ref. . The problem of the cutoff is much less severe with this interaction than with a contact interaction, because the smoothness of $`dU_0/dr`$ cuts off the radial integrals. Our preliminary calculations show that the separable pairing interaction can reproduce results of the density-dependent pairing interaction reasonably well once the strength $`\chi `$ is adjusted. For the purpose of the global binding systematics, the $`A`$-dependence of the strength has yet to be sorted out.

## VI Summary

Our goal is to develop a better microscopic theory for the nuclear binding systematics. In this paper, we argued that the RPA approach provides a promising and computationally tractable way to include correlation effects in a global model, going beyond the mean-field approximation. We have shown that the HF + RPA approach works indeed well using simple Hamiltonian models for correlation associated with broken mean-field symmetries, namely the center-of-mass localization, rotation, and pairing. The RPA equation can be easily solved for a separable residual interaction. For example, the Lanczos method is quite efficient to solve the RPA equation for a separable interaction. This contrasts to other popular computational methods such as the generator coordinate method and the variation-with-projection method, which are rather complicated to apply. A work is now in progress to tackle the nuclear masses using the HF+RPA approach discussed in this paper.

## Acknowledgment

We thank P.-G. Reinhard, W. Nazarewicz, B.A. Brown, and A. Bulgac for discussions. This work was supported by the U.S. Department of Energy under Grant no. DE-FG03-00-ER41132.