text:
# Picture of the low-dimensional structure in chaotic dripping faucets

## Abstract

Chaotic dynamics of the dripping faucet was investigated both experimentally and theoretically. We measured continuous change in drop position and velocity using a high-speed camera. Continuous trajectories of a low-dimensional chaotic attractor were reconstructed from these data, which was not previously obtained but predicted in our fluid dynamic simulation. From the simulation, we further obtained an approximate potential function with only two variables, the drop mass and its position of the center of mass. The potential landscape helps one to understand intuitively how the dripping dynamics can exhibit low-dimensional chaos.

It is well known that the beat of a dripping water faucet is not always regular and exhibits complex behaviour including chaos . Chaotic dripping was originally suggested by Rössler  and experimentally confirmed by Shaw and his collaborators . Since then, many experimental and theoretical studies have established the dripping faucet as a sort of paradigm for chaotic systems . Most previous studies have involved measuring the time interval $`T_n`$ between successive drips, because the dripping time is easily measured using a drop-counter apparatus . The time intervals are then plotted in pairs $`(T_n,T_{n+1})`$ for each $`n`$ to give a return map. Because the return maps typically appear low dimensional ($``$ one-dimensional), the behavior is often described by a simple dynamical model composed of a variable mass and a spring . In this mass-spring model, a mass point, whose mass increases linearly with time at a given flow rate $`Q`$, oscillates with a fixed value of the spring constant $`k`$; and a part $`\mathrm{\Delta }m`$ of the total mass $`m`$ is removed when the spring extension exceeds a threshold, which describes the breakup of a drop. Although the model exhibits chaotic return maps similar to those obtained experimentally, its empirical nature means that it does not provide a unified explanation for the complex behavior of the real dripping faucets, typically seen in bifurcation diagrams (plotting of $`T_n`$ vs. the control parameter $`Q`$). Thus, the connection between the low-dimensional dynamical system and a presumably infinite-dimensional fluid dynamical system remains elusive.

Recently, we carried out fluid dynamic computations (FDC) based on a new algorithm involving Lagrangian description and succeeded in reproducing not only the time-dependent shapes of the drops, but also various characteristics of nonlinear dynamical systems, such as period-doubling, intermittency, hysteresis, observed in experiments of dripping faucets . A detailed analysis of our FDC made it possible to improve the mass-spring model by using more realistic approximations including taking into account the mass dependence of the spring constant, and the necking process . This improved mass-spring model, described only by three variables: mass $`m`$, its position $`z_\mathrm{G}`$ and the velocity $`\dot{z}_\mathrm{G}`$, closely described the dynamical behavior exhibited by the FDC. In particular, qualitatively good agreement among bifurcation diagrams (especially the repeating structure which reflects the oscillation frequency of each drop), obtained by (i) the experiment , (ii) the FDC  and (iii) the improved mass-spring model , over a wide range of the control parameter $`Q`$ provided a theoretical model of the basic low-dimensional structure inherent in the system, which is presented here.

An experimental chaotic attractor is shown in Fig. 1. In addition to the discrete-time variable $`T_n`$, the mass $`m(t)`$ of the liquid suspended by the nozzle, the position of the center of mass $`z_\mathrm{G}(t)`$ and velocity $`\dot{z}_\mathrm{G}`$ were also measured. As seen in Fig. 1(b), $`z_\mathrm{G}(t)`$ oscillates because the surface tension works as a restoring force. Similar oscillation was observed in both our FDC and improved mass-spring model . Note that the frequency of the oscillation during each time interval $`T_n`$ is determined by the flow rate. At $`Q=0.24`$ g/s, $`z_\mathrm{G}`$ oscillates $`6\pm 1/2`$ times before the breaking up and the phase of the oscillation in the plane $`(z_\mathrm{G},\dot{z}_\mathrm{G})`$ fluctuates chaotically at breakup. A return map of $`\{z_{\mathrm{peak},n}\}`$ for every fourth peak of $`z_\mathrm{G}(t)`$ is presented in Fig. 1(c). The return map of $`T_n`$ (Fig. 1(a)) looks like an entangled string and thus the map function $`T_{n+1}=G(T_n)`$ cannot be single-valued. In contrast, the map function $`z_{\mathrm{peak},n+1}=F(z_{\mathrm{peak},n})`$ is approximately single-valued, although broadened by ’noise’. Low-dimensional multi-valued map functions of $`T_n`$ have often been obtained in both experiments , and numerical simulations . For the improved mass-spring model , the map function $`m_{\mathrm{r},n+1}=H(m_{\mathrm{r},n})`$ ($`m_{\mathrm{r},n}`$ being the remnant mass just after the $`n`$-th breakup), as well as the map function $`F(z_{\mathrm{peak},n})`$ is single-valued, while $`G(T_n)`$ is generally multi-valued. This is because the mapping of the time interval $`T_n`$ is not topologically conjugate to that of the dynamical variable $`m_{\mathrm{r},n}`$. The FDC results also showed that $`H(m_{\mathrm{r},n})`$ is approximately single-valued even when $`G(T_n)`$ is multi-valued. Thus, $`T_n`$ does not directly specify the state of the system despite its easy measurement.

Measurement of the continuous-time variables $`\{z_\mathrm{G}(t),\dot{z}_\mathrm{G}(t),m(t)\}`$ made it possible to visualize the trajectory of the chaotic attractor in a continuous state space for the first time. The projection of the attractor in the plane $`(z_\mathrm{G},\dot{z}_\mathrm{G})`$ corresponding to the chaotic motion in Fig. 1 is shown in Fig. 2(top). On average, $`z_\mathrm{G}(t)`$ oscillated six times during each time interval $`T_n`$ (see Figs. 1(b) and 2(top)). The cross section of the attractor in Fig. 1(a) looks nearly one-dimensional and explicitly shows that the drop formation is well characterized by a few state variables. Figure 2(b) is a blow-up of the region S in which trajectories starting with slightly different remnant masses separate rapidly from each other. The transition from the oscillating process to the ’necking’ process (i.e., stretching of the attractor) occurs in this region S . In Fig. 2(bottom) the FDC results are presented including profiles of the evolving drops (Fig. 2(d)), which show qualitatively good agreement with the experimental observations. The cross section of the theoretical attractor (Fig. 2(c)) also exhibits the low-dimensionality. An attractor in the space $`(z_\mathrm{G},\dot{z}_\mathrm{G},m)`$ with a similar spiral structure has already been obtained using the improved mass-spring model . These results strongly support the use of the three state variables: $`z_\mathrm{G}`$, $`\dot{z}_\mathrm{G}`$ and $`m`$ to describe the drop motion.

The motion of the drop is subjected to gravitational force and surface tension. Since the surface energy depends on the shape of the drop, the total potential energy $`U`$ (i.e., gravitational plus surface) should be a function of the many degrees of freedom of the liquid. Our FDC show, however, that approximating $`U`$ as a function of two variables, $`m`$ and $`z_\mathrm{G}`$, yields a conceptually clear picture for the basic low-dimensional structure of the system. Cross sections of $`U/m`$, the potential energy per unit mass, at several fixed values of $`m`$ for the chaotic attractor shown in Fig. 2(bottom) are presented in Fig. 3(a). Poincaré cross sections at the same values of $`m`$ as Fig. 3(a) are given in Fig. 3(b). Although the cross sections in Fig. 3(a) look two-folded, if they are approximated as single-valued functions, then a sheet of the potential surface $`U(m,z_\mathrm{G})/m`$ is obtained as shown in Fig. 3(c). The surface is characterized by a U-shaped valley and a ridge which converge as $`m`$ increases and have totally merged when $`m=m_{\mathrm{crit}}`$, the maximum mass of the static stable shape. We have also carried out a numerical experiment based on FDC with $`Q=0`$ by fixing the value of $`m`$. Various artificial shapes of the drops were used as initial conditions, and it was found that the point $`(m,z_\mathrm{G})`$ rapidly converged onto a well-defined potential surface $`U(m,z_\mathrm{G})/m`$ similar to that shown in Fig. 3(c) .

The potential surface for the closed system (i.e., the pendant drop for $`Q=0`$) is useful for understanding the dynamics of the open system ($`Q>0`$). The minimum and a maximum of $`U(m,z_\mathrm{G})/m`$ for the $`Q=0`$ system correspond to a stable fixed point (static equilibrium shape) and a saddle point (unstable equilibrium shape, see Fig. 3(c)), respectively. These are indicated by red points in the top two panels of Fig. 3(b). When $`m=m_{\mathrm{crit}}`$ ($`=0.144`$ g), the stable fixed point disappears, and the necking process begins leading to the pinching off of the drop. When $`m<m_{\mathrm{crit}}`$, the trajectories of the drops oscillate in the valley of $`U(m,z_\mathrm{G})/m`$, which corresponds to the spiral structure of the attractor in Fig. 2. The flow of the $`Q=0`$ system as depicted in Fig. 3(b), also suggests that the state point of the drop causes the spiral orbits in the phase space $`(z_\mathrm{G},\dot{z}_\mathrm{G},m)`$. The chaotic dynamics, especially the stretching and folding of the attractor, can be interpreted as follows. The stable fixed point approaches the saddle point as $`m`$ increases. The state point thus visits the vicinity of the saddle point after the spiral motion, which causes the instability of the orbit, or equivalently, the stretching of the attractor (as in Fig. 2(b)). The Poincaré section of the chaotic attractor (Fig. 3(b)), which rotates around the stable fixed point with increasing $`m`$, is gradually crushed and asymptotes to a single line through the necking process.

The above analysis leads to the motion of the dripping faucet being described in terms of the approximate potential function $`U`$ as

$$\ddot{z}_\mathrm{G}=\frac{U(m,z_\mathrm{G})/m}{z_\mathrm{G}}\frac{(\gamma +Q)}{m}\dot{z}_\mathrm{G},\dot{m}=Q,$$
(1)
where $`\gamma `$ is the damping constant. This equation describes the motion of a particle with unit mass in 2-dimensional space $`(m,z_\mathrm{G})`$, where the particle is forced to have a constant velocity $`\dot{m}=Q`$ in the direction of the $`m`$-axis. Our improved mass-spring model corresponds to a local quadratic approximation for the potential near the bottom $`z_\mathrm{G}=z_0`$ of the valley in which the $`m`$-dependence of the spring constant $`k(m)^2U(m,z_0)/z_0^2`$ has been included . The present description is thus a natural extension of the improved mass-spring model. As we have seen, however, the picture of the global landscape of the potential is much clearer than the improved mass-spring model for intuitively understanding the dynamics of the system. If the flow rate $`Q`$ is small enough so that the initial oscillation after the breakup is damped, the particle (equivalently the state point of the drop) goes along the bottom of the valley as $`m`$ increases. At $`m=m_{\mathrm{crit}}`$, where the valley merges into the ridge, the particle’s motion loses its stability (start of the necking process) and rapidly approaches the breakup point (on the broken red line in Fig. 3(c)). Since the instability always begins at the same point $`m=m_{\mathrm{crit}}`$ in this case, the size of each drop is uniform. If, on the other hand, the oscillation of the drop affects the start of necking, the dripping motion exhibits a variety of flow rate-dependent periodic and chaotic patterns. Note that the particle generally gets over the ridge before $`m`$ reaches $`m_{\mathrm{crit}}`$ when it is oscillating. At a flow rate resulting in a chaotic motion, two trajectories starting with slightly different $`m`$ values (different remnant masses) are initially close to each other. As $`m`$ increases, however, one trajectory may get over the ridge at a certain $`m`$ value, while the other remains in the valley owing to a small difference in $`\dot{z}_\mathrm{G}`$ (two red lines in Fig. 3(c)).

Similar potential landscapes are expected for other subjects involving drop formation like ink-jet printing  and atomization . The potential surface describing the motion in the attractor was confirmed by FDC to have essentially the same topology even for high viscosity liquids and for microfluidic systems (drop size $``$ pl). Hence, the present picture of a particle moving along a low-dimensional potential surface is relevant to modeling drop formation in general.

We are grateful to Professor W. S. Price and Professor A. L. Salas-Brito for careful reading of the manuscript.


format: text

text:
# Chandra observations of the galaxy cluster A478: the interaction of hot gas and radio plasma in the core, and an improved determination of the Compton y-parameter

## 1 Introduction

A478 is an X-ray-luminous cluster at a redshift of 0.088 (Zabludoff, Huchra & Geller 1990). Previous X-ray observations with Ginga, Einstein, and ROSAT suggested that it contains one of the largest central cooling flows with estimates up to $``$ 1000 M yr-1 (Johnstone et al. 1992; Allen et al. 1993; White et al. 1994; Allen 2000). These observations also indicated excess X-ray absorption above the Galactic value. The X-ray morphology of A478 implies it is a very relaxed cluster (e.g., Allen et al. 1993), although it is elongated in the Northeast-Southwest direction (e.g., White et al. 1994). This elongation is along the same direction as the cluster galaxies (Bahcall & Sargent 1977), as well as the light distribution of the cD galaxy (White et al. 1994). A478 also has a marginal CO detection of 4.5$`\pm `$2.2 $`\times `$ 109 M (Edge 2001). A strong SZE decrement is detected in this cluster, making it an excellent target to measure the Hubble constant (Myers et al. 1997, M97 hereafter; Mason et al. 2001, M01 hereafter).

In this Letter, we present initial results from a recent Chandra 42 ks ACIS observation of A478, with emphasis on the central X-ray gas/radio source interaction and an improved result for the predicted Compton y-parameter using the temperature and density profiles measured with Chandra. A detailed analysis of the cooling flow will be presented in another paper. Throughout this paper we assume H0 = 70 km s-1 Mpc-1, $`\mathrm{\Omega }`$M=0.3, and $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.7 unless specified. At a redshift of 0.088, the luminosity distance to A478 is 403 Mpc and 1′′=1.65 kpc.

## 2 Chandra observation & data reduction

A478 was observed for 42.4 ksec on January 29, 2001 by Chandra with the Advanced CCD Imaging Spectrometer (ACIS). The observation used CCDs I2-I3 and S1-S4 with the cluster centered on S3. The data were telemetered in Faint mode and ASCA grades 1, 5 and 7 were excluded. Known bad columns, hot pixels, and CCD node boundaries also were excluded. We investigated the background light curves from the back-illuminated (BI) S1 CCD and from the front-illuminated (FI) CCDs, using the standard software1. No background flares were found in this observation. The streak events in the S4 CCD were removed by CIAO DESTREAK software.

In our analysis we used period D blank field background data by M. Markevitch 111http://cxc.harvard.edu/contrib/maxim/bg/index.html. The particle background level was measured in the hard band (7 - 11 keV) for the S1 CCD (far from the cluster) and in PHA channels 2500-3000 ADU for all CCDs. We found that the particle background level was 10% lower than that of the period D background data. This is within the uncertainties for the background files. Moreover, the observation was taken at the time that S3 background rate was the lowest 222http://cxc.harvard.edu/ccw/proceedings/presentations/m\_markevitch/pg002.html. Thus, we decreased the background normalization by 10% to fit the particle background level at the time of the A478 observation as recommended in Markevitch et al. (2002). This correction has only a small effect on the best-fit temperature (e.g., 10% for the outermost annulus in $`\mathrm{\S }`$4). Note that although the correction should apply only to the particle component of the total background and not to the cosmic x-ray background or CXB, the particle component is dominant at energies greater than $``$ 3 keV which are most important for measuring the gas temperature. This small rescaling of the CXB component has little effect on the spectral fitting (e.g., $``$ 5% for the outermost annulus). We also corrected for the readout effect in the S3 CCD caused by the bright cluster core (see Markevitch et al. 2000 for details). Two corrections were made to the ACIS low energy quantum efficiency (QE). The first uses CXC (Chandra X-ray Center) contributed software to correct for the QE degradation 333http://cxc.harvard.edu/cal/Links/Acis/acis/Cal\_prods/qeDeg/index.html (Plucinsky et al. 2002). The second corrects the QE by an empirical factor of 0.93 below 1.8 keV in the FI CCDs to improve the cross-calibration with the BI CCDs 444http://asc.harvard.edu/cal/Links/Acis/acis/Cal\_prods/qe/12\_01\_00/ (also see Markevitch & Vikhlinin 2001 for details). The tools calcrmf and calcarf by Vikhlinin were used to generate response files.

The calibration files used correspond to CALDB 2.15 from the CXC. The uncertainties quoted in this paper are 90% confidence intervals unless specified. The solar abundance table by Anders & Grevesse (1989) and a cluster redshift of 0.088 were used in the spectral fits and computations.

## 3 Central region of A478

On large scales from 50 kpc to several Mpc, A478 appears quite relaxed, although the X-ray morphology is elliptical rather than circular (Fig. 1a). The superior angular resolution of Chandra allows us to examine structures on scales of several kpc at the center of A478. In the central 20 kpc (Fig. 1b), the cluster shows significant structures. The X-ray emission appears double-peaked (separated by $``$ 8 kpc) with the two peaks lying east and west of the cD galaxy nucleus. In contrast to the double peaks, there are depressions North and South of the cD nucleus. To better quantitatively show the depressions, we measured the surface brightness profiles in four directions from the nucleus within $``$ 15 kpc (Fig. 2). The surface brightness in the North and South sectors is 20% - 30% lower than that to the West and East. Beyond $``$ 20 kpc, no significant difference between the surface brightness of the sectors is seen. Since the X-ray depressions seen in A478 are similar to the X-ray cavities observed in other clusters (e.g., Hydra A, McNamara et al. 2000) where central radio sources are certainly responsible, we investigated the radio activity of the central source in A478.

A 1.4 GHz VLA image of A478, obtained in 630 s with the A array on March 1, 1998, is shown as isointensity contours in Fig. 1b. The central radio source is composed of a bright point-like source at the position of the cD nucleus, with weak extensions (lobes) towards the Northeast and Southwest. The measured 1.4 GHz total flux is about 30 mJy, while the lobes contribute $``$ 35% of this. The radio lobes are quite small: $``$ 4 - 5 kpc extension for the Southwest one and $``$ 3 - 4 kpc extension for the Northeast one (the VLA angular resolution at 1.4 GHz is 1.5′′, or 2.5 kpc). They lie roughly along the direction of the X-ray depressions, although the X-ray depressions are much larger than the detected radio lobes. The two X-ray peaks are located along the East and West sides of the radio source. The correspondence of radio lobes and X-ray structure suggests that the X-ray structure was created by the interaction of the hot gas with the radio lobes. Similar X-ray structures (X-ray cavities) have been found at the centers of a number of cooling flow clusters (Hydra A - McNamara et al. 2000; Perseus cluster - Boehringer et al. 1993, Fabian et al. 2000; A2052 - Blanton et al. 2001; A2597 - McNamara et al. 2001; A4059 \- Heinz et al. 2002). Generally, they are thought to form when the hot gas is displaced by the expansion of radio lobes or rising buoyant bubbles (e.g., Churazov et al. 2001; Soker, Blanton & Sarazin 2002). For consistency, hereafter we also call the X-ray structure in A478 cavities.

No X-ray point source is detected at the position of the nucleus of the cD galaxy, where a radio point source is detected (Fig. 1b). Assuming a $`\mathrm{\Gamma }`$=1.7 power law modified by an absorption of 2.6$`\times `$1021 cm-2 (the best-fit, see $`\mathrm{\S }`$4) and a cluster redshift of 0.088, the 3 $`\sigma `$ upper limit on the point source luminosity is 3$`\times `$1041 ergs s-1 (0.5 - 10 keV). If the central source is as highly obscured as the one in Hydra A (McNamara et al. 2000), the upper limit will be higher (e.g., 1.1$`\times `$1042 ergs s-1 for NH=2.6$`\times `$1022 cm-2, which is comparable to the source detected in Hydra A).

Spectra of the two X-ray peaks were extracted from 10′′ (major axis) $`\times 4^{\prime \prime }`$ (minor axis) ellipses centered on them. We fit each with a MEKAL model and found consistent temperatures and abundances: 3.2$`\pm `$0.2 keV and 0.60$`{}_{0.17}{}^{}{}_{}{}^{+0.20}`$ for the eastern peak, 3.1$`\pm `$0.2 keV and 0.65$`{}_{0.18}{}^{}{}_{}{}^{+0.21}`$ for the western peak. A hardness ratio map (2 - 5 keV / 0.5 - 2 keV) for the central 60 kpc$`\times `$60 kpc region is shown in Fig. 3. The coolest region contains the two peaks and regions between them. The overall temperature increases outwards from the core. There is no indication at shock heated regions around the X-ray cavities. The spectrum of the central 20 kpc (radius) region was fit with a MEKAL + MKCFLOW model to search for any cooling flow component. The inclusion of a cooling flow model (MKCFLOW) did not significantly improve the fit, if we allow the gas to cool to zero. However, if we use a lower temperature limit of $``$ 1 keV, the fit improves significantly. Thus, while there is a clear signature of gas cooling towards the center (see $`\mathrm{\S }`$4), the spectral fits do not require gas at very low temperatures.

## 4 Analysis of the temperature, absorption and heavy element abundance profiles

We first compare the integrated spectrum of A478 from the entire S2, S3 and S4 CCDs with previous results. Data in the 0.7 - 8 keV energy band are used. The results are NH = (0.259$`\pm `$0.003)$`\times `$1022 cm-2, T = 7.18$`\pm `$0.11 keV, and abundance 0.37$`\pm `$0.02 solar. From ASCA GIS 1 - 9 keV data, White (2000) obtained a temperature of 6.6$`\pm `$0.4 keV and an abundance of 0.31$`\pm `$0.05 with a fixed absorption of 3$`\times `$1021 cm-2. If we fix the absorption at their value, we find a consistent temperature T = 6.62$`\pm `$0.09 keV and an abundance of 0.37$`\pm `$0.02.

Although in X-rays A478 is somewhat elliptical, it is useful to derive radial properties of the temperature, absorption and heavy element abundance. For the profiles, we centered all annuli at R.A. = 04h13m25s.3, decl. = 10$`^{}`$ 27$`^{}`$ 53$`^{\prime \prime }`$and required each annulus to contain 6000 - 8000 source counts. Each annulus was fit by a MEKAL model, with the absorption, temperature and abundance as free parameters. The reduced $`\chi ^2`$ ranged from 0.8 to 1.4 and most were around 1.0 (for 125 - 245 degrees of freedom), which is acceptable considering the large number of source counts and the current ACIS calibration. The derived temperature profile is shown in Fig. 4. The temperature falls from $``$ 8.5 keV at large radii to $``$ 3 keV at the center. Between radii of $``$ 110 kpc and $``$ 500 kpc, gas temperatures are consistent with a constant temperature of $``$ 8.5 keV. To better constrain the absorption and abundance, we used larger radial bins and obtained the results shown in Fig. 5. The abundance increases from $``$ 0.3 solar at several hundred kpc to $``$ 0.6 solar at the center. The absorption also increases towards the center. The spectra in the regions outside $``$ 500 kpc (in CCDs S2, S4, I2, and I3) were fitted simultaneously with a fixed absorption of 0.24$`\times `$1022 cm-2 (extrapolated from the absorption profile shown in Fig. 5). The best-fit temperature is 6.45$`{}_{0.75}{}^{}{}_{}{}^{+0.80}`$ keV and will be taken as the “background” temperature in the deprojection analysis ($`\mathrm{\S }`$5).

## 5 Deprojection analysis

Although the X-ray morphology of A478 is elliptical rather than circular, the ellipticity is small ($``$ 0.2) and a spherical approximation serves as a simple way to do the deprojection and other analyses. Thus, based on the 2D temperature profile obtained in $`\mathrm{\S }`$4, we performed a deprojection analysis to obtain the 3D temperature profile and abundance profile (see e.g., David et al. 2001 for details of the spectral deprojection). Nine radial bins were chosen to do the spectral deprojection. These were the same ones used above to derive the absorption and abundance profiles. To account simply for the slight decrease of absorption with radius, we used two absorption values, 0.28$`\times `$1022 cm-2 for annuli 1-3 and 0.25$`\times `$1022 cm-2 for annuli 4-9 (Fig. 5), derived by fitting a single absorption for these regions ($`\chi ^2`$ = 8.1/7). The deprojection result is shown in Fig. 4. For regions within a radius of $``$ 110 kpc, we also fit the 3D temperatures with a power law given by T(r) = T0(r/10kpc)α. The best fit is T0 = 2.79$`\pm `$0.09 keV and $`\alpha `$ = 0.412$`\pm `$0.017 (1$`\sigma `$ errors). Deprojected abundances are shown in Fig. 5.

The surface brightness profile of A478 is shown in Fig. 4. There is at least one inflection point at $``$ 40 kpc. If the surface brightness profile is fit by a single $`\beta `$-model, the central excess is quite significant and the fit is poor. Thus, we use a deprojection technique to determine the electron density profile rather than fitting with a multivariate function (see a recent example in David et al. 2001). This technique converts the observed surface brightness to the electron density starting in the largest annulus and then determines the density at progressively smaller radii, after removing the projected emission from larger radii. The fraction of cluster emission in all radial bins arising from gas outside the outermost annulus ($``$ 5 radius) is corrected using the power law fit to both Chandra and PSPC surface brightness profiles in regions outside the outermost annulus. The deprojected temperature and abundance profiles were combined to produce an emissivity profile (with uncertainties) used in the conversion. We performed 400 Monte Carlo simulations to estimate the errors.

The deprojected electron density profile is shown in Fig. 4, as well as the pressure profile. The density profile gradually steepens from n$`{}_{\mathrm{e}}{}^{}`$ r-0.25 at 10 kpc to n$`{}_{\mathrm{e}}{}^{}`$ r-1.65 at 250 kpc, and there is no flat density core. The isobaric cooling time is 4$`\times `$108 yr in the innermost 10 kpc and less than a Hubble time ($``$ 1010 yr) within the central 140 kpc. The pressure reaches a maximum at r $``$ 15 kpc, at roughly the outer edge of the X-ray cavities.

## 6 Discussion

### 6.1 Central X-ray gas / radio source interaction

As shown in $`\mathrm{\S }`$3, there are two X-ray cavities in A478. The alignment of the radio lobes and X-ray cavities implies that the X-ray gas was swept out of the cavities by expanding radio plasma. We suggest that the X-ray double peaks result from a “donut” configuration for the hot gas around the central nucleus, with the plane of the “donut” perpendicular to the direction of the radio lobes, while the central parts of the “donut” were displaced by the expanding radio plasma. The absence of shocked regions in the A478 core (Fig. 3) is consistent with the radio plasma expanding subsonically, which is similar to what was found in other clusters (e.g., Hydra A and Perseus). If we assume the X-ray cavities are shaped like circular cones, the minimum energy required to sweep the gas out of the current cavities is pV $``$ 3 $`\times `$ 1058 ergs. Churazov et al. (2001) showed that a large and strongly underdense bubble will rise with a velocity comparable to the Keplerian velocity. Therefore, as a first order approximation, we take the Keplerian velocity at 10 kpc to represent the rising velocity of the bubble ($``$ 460 km/s, the mass is estimated from the deprojected temperature and density under the assumption of hydrostatic equilibrium). Thus, the time to form the current cavities is $``$ 3 $`\times `$ 107 yr. Therefore, the minimum average power from the central source to create the X-ray cavities is $``$ 3 $`\times `$ 1043 ergs s-1, which is at least 10 times larger than the current total radio power from the central source. Since the total radio power is simply a lower limit to the kinetic power of the jet (e.g., Perseus — Fabian et al. 2002), more radio information is required to conclude whether the current kinetic power of the jet could create the X-ray cavities. However, the observed radio lobes are considerably smaller than the X-ray cavities, which may imply a fading of the radio source within $``$ 107 — 108 yr.

Five other clusters also contain X-ray cavities. Although each is very complicated, in Fig. 6 we compare the cavity size and the 1.4 GHz luminosity of the central radio source for all of them. Although the kinetic power of the radio source is generally unknown, Heinz et al. (2002) argued that the upper limit on the current kinetic power was proportional to the current 1.4 GHz luminosity. In A478 the 1.4 GHz luminosity of the central source is the smallest in this sample (less than 7% of any others) and the radio lobes are the smallest. Thus, A478 serves as an important example for understanding X-ray gas / radio plasma interaction in the center of cooling-flow clusters. The size of the X-ray cavities in A478 are $``$ 3/4 the size of those in Perseus (Fig. 6), but their 1.4 GHz luminosities differ by $``$ 20. In Fig. 6, if we use the three clusters where the central radio sources are powerful enough to create the X-ray cavities (Hydra A, A2052 and Perseus) to define a line, A478 and A4059 have much lower 1.4 GHz luminosities ($``$ 10%) than the values predicted by that line. Heinz et al. (2002) suggested that the central radio source in A4059 had faded by about an order of magnitude on a timescale less than 108 yr, since the current kinetic power they estimated is too small to create the X-ray cavities. If we adopt their method, A478 is very similar to A4059: in both clusters, the estimated current kinetic power is too small to create the X-ray cavities. Thus, A478 and A4059 may be in a different “phase” of their active nuclei cycles than the other sources.

### 6.2 An improved estimation of the Compton y-parameter

To date, X-ray observations have been combined with SZE measurements to estimate H0 for over 20 galaxy clusters (e.g., M97; Hughes & Birkinshaw 1998; Reese et al. 2000; Mauskopf et al. 2000). However, these analyses generally used only single temperatures (or simple temperature profiles) and density profiles derived by simple $`\beta `$-models. These can introduce significant errors in some cases, such as for cooling flow clusters where a single $`\beta `$-model usually fails to describe the central region, or clusters with significant temperature gradients (Hughes & Birkinshaw 1998). It is ideal and more accurate to apply directly the observed density and temperature profiles to models. In this paper, for the first time, we apply the deprojected density and temperature profiles and use direct integration to derive the Compton y-parameter

$$y=\frac{k\sigma _T}{m_ec^2}\underset{\mathrm{}}{\overset{\mathrm{}}{}}T_en_e𝑑\zeta .$$
(1)
In this section, we use H0 = 100 h km s-1 Mpc-1 and assume spherical symmetry of the cluster X-ray emission. In $`\mathrm{\S }`$5, we derived the deprojected density profile within the inner 5 (Fig. 4). For the outer regions, we combine Chandra and ROSAT PSPC images and fit the surface brightness profile with a $`\beta `$-model. The derived parameters are $`\beta `$ = 0.68$`\pm `$0.01 and core radius r0 = 1.02$`\pm `$0.07 (1 $`\sigma `$ errors). The density profile beyond 5 was obtained from this $`\beta `$-model. Combining the deprojected density profile with the deprojected temperature profiles (ne(r) and Te(r) in $`\mathrm{\S }`$5), the Compton y-parameters can be computed for different lines of sight. The uncertainties are estimated from Monte Carlo simulations. M97 and M01 measured the Compton y-parameter based on their Owens Valley 5.5 m data at 32 GHz. To compare the predicted y-parameter from our X-ray observations with their measured value, we correct for specifics of the 5.5 m telescope (the OVRO beam and the switching pattern; see eq. 32 - 37 of M97 for details). After those corrections, we obtain a predicted y-parameter specific for the 5.5 m telescope: (6.16$`\pm `$0.06) $`\times `$ 10-5 (ypred, 1 $`\sigma `$ error). The systematic errors from X-ray observations can further contribute up to 4% error (mainly from the uncertainties of temperatures). This value has a much smaller statistical uncertainty than the estimates in M97 and M01 (4.2$`\pm `$1.0 $`\times `$ 10-5 and 6.05$`\pm `$0.54 $`\times `$ 10-5). Moreover, our derivation of the Compton y-parameter is more accurate than the simple methods adopted in M97 and M01, which were limited by the X-ray data at that time.

We can use this result to obtain an estimate of H0. M01 obtained the observed y-parameter: (7.77$`\pm `$0.58) $`\times `$ 10-5 (1 $`\sigma `$ error). However, after we examined the relativistic correction they used, which depends on the measured temperature, we changed the relativistic correction factor to 1.026 following Itoh, Kohyama & Nozawa (1998). Thus, the observed y-parameter should be (7.71$`\pm `$0.57) $`\times `$ 10-5 (yobs, 1 $`\sigma `$ error). In principle, h = (ypred / yobs)2. Including the uncertainty from intrinsic CMB anisotropies (see details in M01), we obtain a Hubble constant H = 64$`{}_{18}{}^{}{}_{}{}^{+32}`$ km s-1 Mpc-1 (1 $`\sigma `$ error). Intrinsic CMB anisotropies and the uncertainty of the radio observation are the main sources of the statistical uncertainties. Besides X-ray calibration uncertainties mentioned above, other systematic uncertainties include calibration uncertainties in the radio observations (3% - M01), uncertainty from peculiar motion (2% for 500 km s-1 peculiar motion), clumping of X-ray gas ($``$ -15%), undetected radio sources ($``$ 10%), and unknown geometry ($``$ 10% based on the analysis in Hughes & Birkinshaw, 1998). A478 only serves as an example. As more SZE clusters are observed by Chandra and XMM-Newton, and as better radio observation are available, a systematic analysis will provide a better constraint on H0.

The results presented here are made possible by the successful effort of the entire *Chandra* team to build, launch, and operate the observatory. We acknowledge helpful discussions with L. David, D. Harris, M. Markevitch, and A. Vikhlinin. We thank the Smithsonian Institution for support. We thank the referee for comments. This study was also supported by NASA contract NAS8-38248.


text:
# Energies of the ground state and first excited 0⁺ state in an exactly solvable pairing model

## I Introduction

The random-phase approximation (RPA)  is one of the most popular method in the theoretical microscopic study of nuclear structure. It includes the correlations beyond the mean field models such as the Hartree-Fock (HF) approximation with a phenomenological interaction, e.g. the Skyrme interaction. Being also a computationally simple method, the RPA serves as a powerful tool in treating all the excited states in nuclei, which are beyond the reach of the full diagonalization within a shell-model basis. The quasiparticle RPA (QRPA), in which the quasiparticle correlations and excitations are considered, has shown its prominent role in treating the open-shell nuclei, where the superfluid pairing correlations are important.

Recent developments in nuclear structure including $`\beta `$-decay and double-beta ($`\beta \beta `$) -decay physics, the study of various types of giant resonances, and the prospect of using radioactive beams to study nuclei far from the line of $`\beta `$-stability have sparked off a renewed interest in the efficiency of the RPA compared to other microscopic calculations . More refined and exact treatments of the pairing problem in nuclei have also been proposed in .

Recently the accuracy of the RPA in describing the binding energies has been tested using the HF+RPA calculations within schematic exactly solvable models  as well as for nuclei throughout the $`sd`$-shell and the lower $`pf`$-shell . The tests using exactly solvable Lipkin models  has shown that the HF+RPA (QRPA) calculations yield a very good ground state energy except for the region of the pairing interaction around the point where the pairing gap collapses. Mean while, the tests using more realistic shell-model Hamiltonian have pointed out that the binding energy predicted by RPA is generally, but not always, satisfactory . A number of suggestions have been made to improve the reliability of the RPA. They can be classified into two groups. The first group includes the approaches to improve the treatment of pairing correlations. Among them are the Hartree-Fock-Bogolyubov (HFB) + QRPA, the number projection , the variational method using boson expansion , and the Lipkin-Nogami (LN) method . The LN method is an approximation to remove the fluctuations due to the violation of particle number conservation of the BCS approximation, which leads to the collapses of the pairing gap at a certain critical value of the pairing interaction. The second group includes various treatments of the ground-state correlations beyond RPA, which occur due to the violation of the Pauli principle when the quasiboson approximation is used within the RPA or QRPA. In a way similar to the particle number violation within the BCS approximation, the violation of the Pauli principle by treating the fermion pairs as bosons within the quasiboson approximation leads to the collapse of the RPA at a critical value of the interaction parameter. Various approaches have been proposed to renormalize the RPA (QRPA) to remove this inconsistency .

Since the RPA, first of all, is a theory of excited states, an improvement of the accuracy of the RPA should give, as a first step, a better description for both of the energies of the ground state and the first excited state simultaneously. In the present paper, we are going to propose a recipe for such improvement using existing approximation schemes, namely the LN method for the pairing gap and the renormalized RPA (renormalized QRPA) to treat the ground state correlations beyond RPA. In order to compare the result with the exact solution we limit ourselves to the study of the ground state and the first excited $`0^+`$ state within an well-known exactly solvable two-level model, which was introduced for the first time in Ref. , and widely used in the study of pairing correlations.

The paper is organized as follows. The pairing Hamiltonian applied to the two-level model is discussed in Sec. I. The BCS approximation, the LN method, and the results obtained within these approximations for the pairing gap in the two-level model under consideration are summarized in Sec. II. The RPA, QRPA and their renormalized versions are presented in Sec. III, where the results for the ground-state energy and energy of the first excited $`0^+`$ state are analysed. Also discussed in the same section are results obtained following a recipe which combines the LN-method and the renormalized QRPA (renormalized RPA). The last section summarizes the paper, where conclusions are drawn.

## II Pairing Hamiltonian

We consider the well-known pairing Hamiltonian of nucleons interacting via a pairing force between the time-conjugate orbitals with angular-momentum quantum number $`j`$

$$H=\underset{j}{}ϵ_jN_jG\underset{jj^{}}{}\sqrt{\mathrm{\Omega }_j}\sqrt{\mathrm{\Omega }_j^{}}A_j^{}A_j^{},$$
(1)
where $`ϵ_j`$ are the single-particle energies, and $`G`$ is the pairing constant. The nucleon number operator for the $`j`$ shell is

$$N_j=\underset{m}{}a_{jm}^{}a_{jm}.$$
(2)
The pairing operators $`A_j^{}`$ and $`A_j`$ are given by

$$A_j^{}=\frac{1}{\sqrt{\mathrm{\Omega }_j}}\underset{m>0}{}a_{jm}^{}a_{j\stackrel{~}{m}}^{},A_j=(A_j^{})^{},$$
(3)
where the tilde denotes the time-reversal operation, e.g. $`a_{j\stackrel{~}{m}}^{}=(1)^{jm}a_{jm}^{}`$. These operators satisfy the following commutation relations:

$$[A_j,A_j^{}^{}]=\delta _{jj^{}}\left(1\frac{N_j}{\mathrm{\Omega }j}\right),$$
(4)
$$[N_j,A_j^{}^{}]=2\delta _{jj^{}}A_j^{},[N_j,A_j^{}]=2\delta _{jj^{}}A_j.$$
(5)
The exact solutions of this pairing Hamiltonian have been found in Ref.  and known as the Richarson’s solution. In the present paper, we consider only a simple schematic model. It has $`N`$ particles occupying two levels with the same shell degeneracy $`\mathrm{\Omega }=j+1/2`$ so that $`\mathrm{\Omega }=N/2`$. The upper level simulates the degenerates states with energy $`ϵ/2`$ and magnetic quantum numbers $`m`$ (0$`<m\mathrm{\Omega }`$), while the lower level is for the states with energy $`ϵ/2`$ and the magnetic quantum numbers $`m`$. The distance between two levels is, therefore, $`ϵ`$, which will take the value equal to 1 MeV in all calculations in this paper. The exact solution in this case is easy to be found using the SU(2) algebra, which the operators $`J_+=\sqrt{\mathrm{\Omega }_j}A_j^{}`$, $`J_{}=(J_+)^{}`$, and $`J_0=(N_j\mathrm{\Omega }_j)/2`$ obey. This leads to the diagonalization of the Hamiltonian (1), whose matrix elements are

$`J,M_1^{},M_2^{}|H|J,M_1,M_2=ϵ(M_1M_2)\delta _{M_1,M_1^{}}\delta _{M_2,M_2^{}}`$

$`G[2J(J+1)M_1(M_11)M_2(M_21)]\delta _{M_1,M_1^{}}\delta _{M_2,M_2^{}}`$

$$G\sqrt{J(J+1)M_1(M_11)}\sqrt{J(J+1)M_2(M_2+1)}\delta _{M_1,M_1^{}+1}\delta _{M_2,M_2^{}1}$$
(6)
$`G\sqrt{J(J+1)M_1(M_1+1)}\sqrt{J(J+1)M_2(M_21)}\delta _{M_1,M_1^{}1}\delta _{M_2,M_2^{}+1},`$

with $`J=\mathrm{\Omega }/2`$ and $`JM_iJ`$ ($`i=`$ 1, 2). Among the obtained eigenenergies $`_i`$, the lowest one, $`_0`$, is the exact ground-state energy, while the exact energy $`\omega _{\mathrm{ex}}`$ of the first excited $`0^+`$ state is $`\omega _{\mathrm{ex}}=_1_0.`$

It is convenient to study a fermion system with superfluid pairing using the Bogolyubov transformation from particle operators $`a_{jm}^{}`$ and $`a_{jm}`$ to the quasiparticle ones, $`\alpha _{jm}^{}`$ and $`\alpha _{jm}`$. The quasiparticle representation of the Hamiltonian (1) is given in , which we quote here again for the convenience in further discussions:

$`H=a+{\displaystyle \underset{j}{}}b_j𝒩_j+{\displaystyle \underset{j}{}}c_j(𝒜_j^{}+𝒜_j)+{\displaystyle \underset{jj^{}}{}}d_{jj^{}}𝒜_j^{}𝒜_j^{}`$

$$+\underset{jj^{}}{}g_j(j^{})(𝒜_j^{}^{}𝒩_j+𝒩_j𝒜_j^{})+\underset{jj^{}}{}h_{jj^{}}(𝒜_j^{}𝒜_j^{}^{}+𝒜_j^{}𝒜_j)+\underset{jj^{}}{}q_{jj^{}}𝒩_j𝒩_j^{}.$$
(7)
Here $`𝒩_j`$ is the operator of the quasiparticle number on the $`j`$-shell, while $`𝒜_j^{}`$ and $`𝒜_j`$ are the creation and destruction operators of a pair of time-conjugate quasiparticles:

$$𝒩_j=\underset{m}{}\alpha _{jm}^{}\alpha _{jm},𝒜_j^{}=\frac{1}{\sqrt{\mathrm{\Omega }_j}}\underset{m>0}{}\alpha _{jm}^{}\alpha _{j\stackrel{~}{m}}^{},𝒜_j=(𝒜_j^{})^{}.$$
(8)
Their commutation relations are similar to those for nucleon operators in Eqs. (4) and (5), namely:

$$[𝒜_j,𝒜_j^{}^{}]=\delta _{jj^{}}(1\frac{𝒩_j}{\mathrm{\Omega }j}),$$
(9)
$$[𝒩_j,𝒜_j^{}^{}]=2\delta _{jj^{}}𝒜_j^{},[𝒩_j,𝒜_j^{}]=2\delta _{jj^{}}𝒜_j.$$
(10)
The coefficients in Eq. (7) are

$$a=2\underset{j}{}\mathrm{\Omega }_jϵ_jv_j^2G(\underset{j}{}\mathrm{\Omega }_ju_jv_j)^2G\underset{j}{}\mathrm{\Omega }_jv_j^4,$$
(11)
$$b_j=ϵ_j(u_j^2v_j^2)+2Gu_jv_j\underset{j^{}}{}\mathrm{\Omega }_j^{}u_j^{}v_j^{}+Gv_j^4,$$
(12)
$$c_j=2\sqrt{\mathrm{\Omega }_j}ϵ_ju_jv_jG\sqrt{\mathrm{\Omega }_j}(u_j^2v_j^2)\underset{j^{}}{}\sqrt{\mathrm{\Omega }_j^{}}u_j^{}v_j^{}2G\sqrt{\mathrm{\Omega }_j}u_jv_j^3,$$
(13)
$$d_{jj^{}}=G\sqrt{\mathrm{\Omega }_j\mathrm{\Omega }_j^{}}(u_j^2u_j^{}^2+v_j^2v_j^{}^2)=d_{j^{}j},$$
(14)
$$g_j(j^{})=Gu_jv_j\sqrt{\mathrm{\Omega }_j^{}}(u_j^{}^2v_j^{}^2),$$
(15)
$$h_{jj^{}}=\frac{G}{2}\sqrt{\mathrm{\Omega }_j\mathrm{\Omega }_j^{}}(u_j^2v_j^{}^2+v_j^2u_j^{}^2)=h_{j^{}j},$$
(16)
$$q_{jj^{}}=Gu_jv_ju_j^{}v_j^{}=q_{j^{}j},$$
(17)
where $`u_j`$ and $`v_j`$ are the coefficients of the Bogolyubov transformation. Hereafter we will refer to the terms at the right-hand side (RHS) of Eq. (7), which contain the coefficients $`a_j`$ , $`b_j`$, etc. as the a-term, b-term, etc., respectively.

## III Gap equations

### A BCS approximation

The BCS equation is usually obtained using the variational procedure to get the minimum of the average value of $`H\lambda \widehat{N}`$ ($`\lambda `$ is the chemical potential, $`\widehat{N}=_jN_j`$ is the particle-number operator) over the BCS ground state taken as the quasiparticle vacuum $`|0_\alpha `$, i.e. $`\alpha _{jm}|0_\alpha =`$ 0, where $`|0_\alpha =_{j,m>0}(u_j+v_ja_{jm}^{}a_{j\stackrel{~}{m}}^{})|0`$ with $`a_{jm}|0=`$ 0. Within the BCS approximation, only the a-term in Eq. (7) contributes, which leads to the well-known BCS equation to determine the gap $`\mathrm{\Delta }`$ and chemical potential $`\lambda `$:

$$\mathrm{\Delta }=G\underset{j}{}\mathrm{\Omega }_ju_jv_j,N=\underset{j}{}\mathrm{\Omega }_j\left(1\frac{ϵ_j^{}\lambda }{E_j}\right),$$
(18)
where the single-particle energy is $`ϵ_j^{}=ϵ_j`$ if the self-energy term $`Gv_j^2`$ is neglected, or $`ϵ_j^{}=ϵ_jGv_j^2`$ if the self-energy term is included. The quasiparticle energy is $`E_j=\sqrt{(ϵ_j^{}\lambda )^2+\mathrm{\Delta }^2}`$. The $`u_j`$ and $`v_j`$ coefficients are given as

$$u_j^2=\frac{1}{2}\left(1+\frac{ϵ_j^{}\lambda }{E_j}\right),v_j^2=\frac{1}{2}\left(1\frac{ϵ_j^{}\lambda }{E_j}\right).$$
(19)
The a-term (Eq. (11)) is actually the ground state energy within the BCS approximation, since this is the only term that remains in the average over the quasiparticle vacuum $`|0_\alpha `$, where the second term can be now replaced with $`\mathrm{\Delta }^2/G`$ using the BCS equation (18).

In the present two-level model, using the BCS equations (18) and the property $`u_i^2+v_i^2=`$ 1 with $`i=`$ 1 (lower level) or 2 (upper level), it is easy to see that the quasiparticle energy $`E`$ and the chemical potential $`\lambda `$ are state-independent, namely

$$E_1=E_2=E=G\mathrm{\Omega },\lambda =\frac{G}{2}.$$
(20)
The gap $`\mathrm{\Delta }`$ and the $`u`$ and $`v`$ coefficients are

$$\mathrm{\Delta }=G\mathrm{\Omega }\sqrt{1\left[\frac{G_{\mathrm{cr}}}{G}\right]^2},$$
(21)
$$u_1^2=v_2^2=\frac{1}{2}\left(1\frac{\stackrel{~}{ϵ}}{2G\mathrm{\Omega }}\right)=u^2,v_1^2=u_2^2=\frac{1}{2}\left(1+\frac{\stackrel{~}{ϵ}}{2G\mathrm{\Omega }}\right)=v^2,$$
(22)
where

$$\stackrel{~}{ϵ}=\{\begin{array}{cc}ϵ\hfill & \text{neglecting the self-energy term ,}\hfill \\ & \\ 2\mathrm{\Omega }ϵ/(2\mathrm{\Omega }1)\hfill & \text{including the self-energy term ,}\hfill \end{array}$$
(23)
and

$$G_{\mathrm{cr}}=\frac{\stackrel{~}{ϵ}}{2\mathrm{\Omega }}.$$
(24)
The BCS ground state energy in this model becomes

$`E_{\mathrm{BCS}}=\mathrm{\Omega }ϵ(12v^2)\mathrm{\Delta }^2/GG\mathrm{\Omega }(u^4+v^4)`$

$$=\{\begin{array}{cc}ϵ^2/2G\mathrm{\Delta }^2(11/2\mathrm{\Omega })/GG\mathrm{\Omega }\hfill & \text{neglecting the self-energy term ,}\hfill \\ & \\ ϵ^2/(2G(11/2\mathrm{\Omega }))\mathrm{\Delta }^2(11/2\mathrm{\Omega })/GG\mathrm{\Omega }\hfill & \text{self-energy term included .}\hfill \end{array}$$
(25)
This shows that the effect of the self-energy term on the ground-state energy becomes negligible at a large particle number when $`1/2\mathrm{\Omega }`$ 1. The major drawback of the BCS method is that the BCS wave function is not an eigenstate of the particle number operator $`\widehat{N}`$. Therefore the particle-number fluctuations $`\mathrm{\Delta }N^2=_\alpha 0|\widehat{N}^2|0_\alpha _\alpha |0\widehat{N}|0_\alpha ^2=4_j\mathrm{\Omega }_j(u_jv_j)^2=2\mathrm{\Delta }^2/G^2\mathrm{\Omega }`$  cause the inaccuracy of this method. Using the exact commutation relation (9), and $`{}_{\alpha }{}^{}0|𝒩_j|0_{\alpha }^{}=`$ 0, we see that

$$_\alpha 0|[𝒜_j,𝒜_j^{}^{}]|0_\alpha =\delta _{jj^{}},_\alpha 0|[𝒩_j,𝒜_j^{}^{}]|0_\alpha =0,$$
(26)
which imply that, within the BCS approximation, the quasiparticle-pair operators $`𝒜_j^{}`$ and $`𝒜_j`$ behave like bosons (the so-called Cooper pairs). Such violation of the Pauli principle between quasiparticles has the same origin as the quasiboson approximation used in the RPA, which will be discussed in the next section. This leads to the collapse of the pairing gap at the critical value $`G_{\mathrm{cr}}`$, below which the BCS equations (18) yields the imaginary solution. In the present two-level model, $`G_{\mathrm{cr}}`$ is given by Eq. (24). Such kind of critical behavior inspired a speculation of the existence of a phase transition from the normal phase, where the gap is absent, to the superfluid phase with a nonzero gap. In the absence of the pairing gap, only the sum over the hole ($`h`$) state remains in expression for the ground-state energy (11), as $`v_h=`$ 1 and $`v_p=`$ 0. This gives the HF ground-state energy within the present two-level model as

$$E_{\mathrm{HF}}=\mathrm{HF}|H|\mathrm{HF}=2\underset{j_h}{}\mathrm{\Omega }_{j_h}ϵ_{j_h}G\underset{j_h}{}\mathrm{\Omega }_{j_h}=\mathrm{\Omega }(ϵ+G),$$
(27)
where $`|\mathrm{HF}=_ja_{jm}^{}|0`$ is the HF ground state. However, the superfluid-normal phase transition in a system with a finite particle number is spurious as it does not exist in the exact calculations  as well as in the method using particle-number projection , where the gap is finite at all finite values of $`G`$ . Since carrying out the particle-number projection in calculations using realistic spectra is numerically complicate, a simple approximate number projection has been proposed, which is known as the Lipkin-Nogami (LN) method  and summarized below for the present two-level model.

### B Lipkin-Nogami (LN) method

The LN method has gained a great popularity as it provides a simple and computationally easy way to go beyond the pairing mean-field. Within this method, the particle-number fluctuations are removed by adding the term $`\lambda _2\widehat{N}^2`$, and carrying out the variational procedure over the average value of $`H\lambda \widehat{N}\lambda _2\widehat{N}^2`$ in the quasiparticle vacuum $`|0_\alpha `$. Details of this method is given in Ref. . For the present two-level model, it gives

$$\mathrm{\Delta }_{\mathrm{LN}}=G\mathrm{\Omega }\sqrt{1\frac{\overline{ϵ}^2}{4G^2\mathrm{\Omega }^2}},$$
(28)
$$u^2=\frac{1}{2}\left(1\frac{\overline{ϵ}}{2G\mathrm{\Omega }}\right),v^2=\frac{1}{2}\left(1+\frac{\overline{ϵ}}{2G\mathrm{\Omega }}\right),$$
(29)
where

$$\overline{ϵ}=\frac{2G\mathrm{\Omega }ϵ}{2G\mathrm{\Omega }+\alpha },\alpha =4\lambda _2G.$$
(30)
Setting the factor $`\alpha =`$ 0 recovers the BCS equation without the self-energy term, while setting $`\lambda _2=`$ 0 brings back to the BCS equation including this term. The factor $`\alpha `$ is found substituting $`(uv)^2`$ in the equation for $`\lambda _2`$:

$$\frac{4\lambda _2}{G}=\frac{\mathrm{\Omega }2(uv)^2}{2(2\mathrm{\Omega }1)(uv)^2}.$$
(31)
This leads to the following cubic equation for $`\alpha `$:

$$\alpha (2\mathrm{\Omega }1)[(\alpha +2G\mathrm{\Omega })^2ϵ^2]2G\mathrm{\Omega }ϵ^2=0.$$
(32)
The ground state energy is given by

$$E_{\mathrm{LN}}=\mathrm{\Omega }(u^2v^2)\frac{\mathrm{\Delta }_{\mathrm{LN}}^2}{G}G\mathrm{\Omega }[12(uv)^2]2(\alpha +G)\mathrm{\Omega }(uv)^2.$$
(33)
The gaps obtained within the BCS approximation and LN method using Eqs. (21) and (28) for $`\mathrm{\Omega }=`$ 4 and 8 are plotted against the interaction parameter $`G\mathrm{\Omega }/(2ϵ)`$ in Fig. 1. With decreasing the interaction, the BCS gap $`\mathrm{\Delta }`$ decreases and collapses at $`G_{\mathrm{cr}}`$, whose value for the case when the self-energy term is neglected is smaller than that obtained including the self-energy term. The LN gap $`\mathrm{\Delta }_{\mathrm{LN}}`$, on the contrary, decreases monotonously with decreasing $`G`$ until $`G=0`$, where $`\mathrm{\Delta }_{\mathrm{LN}}`$ vanishes, showing no signature of the superfluid-normal phase transition. Hence, by removing the particle number fluctuations, the LN method erases completely the odd behavior of the BCS gap characterized by this phase transition. The system remains in the superfluid phase at all nonzero values of the pairing interaction. It is interesting to see that, due to the suppression of the self-energy term $`Gv^2`$ by $`4\lambda _2v^2`$ in the LN method, the BCS gap without the self-energy correction is closer to the LN result than the BCS gap including this correction. By comparing Fig. 1 (a) and (b), we also see that the difference between two BCS versions (with and without the self-energy term) decreases significantly when the particle number increases.

## IV QRPA (RPA), its renormalization, and combination with the LN method

The RPA ground state includes $`2p2h`$-, $`4p4h`$-, $`6p6h`$-, etc. excitations on top of the HF ground state. The RPA including the pairing correlations within the quasiparticle representation is usually referred to as the QRPA. The correlations in the QRPA ground state, therefore, are much richer than the $`2p2h`$-, $`4p4h`$-, etc. correlations in the “diffuse” (quasiparticle) ground state $`|0_\alpha `$ created by scattering of particle pairs within the BCS approximation. The LN method takes into account some $`2p2h`$-correlations beyond the BCS approximation by using the $`\lambda _2\widehat{N}^2`$ term. However, this method still uses the same BCS ground state $`|0_\alpha `$ since it approximately expresses the expectation value of the Hamiltonian with respect to the projected state in terms of that with respect to the BCS ground state . The discussion below is conducted within the QRPA, where the RPA is obtained as the limit when $`\mathrm{\Delta }=`$ 0. Consequently, we summarize also the main features of the renormalized QRPA, whose zero-pairing limit is the renormalized RPA. In the present two-level model, the QRPA works at $`GG_{\mathrm{cr}}`$, while the RPA is applied at $`G<G_{\mathrm{cr}}`$.

### A QRPA

The standard QRPA operators, called phonon operators, have the following form in the present two-level model

$$Q_\nu ^{}=\underset{j}{}(X_j^{(\nu )}𝒜_j^{}Y_j^{(\nu )}𝒜_j),Q_\nu =(Q_\nu ^{})^{},$$
(34)
The QRPA ground state $`|0_\mathrm{Q}`$ is defined as the vacuum for the phonon operator, i.e. $`Q_\nu |0_\mathrm{Q}=0=_\mathrm{Q}0|Q_\nu ^{}`$. The $`0^+`$ excited state $`|0^+`$ is obtained by acting $`Q_\nu ^{}`$ on this vacuum, i.e. $`|0^+=Q_\nu ^{}|0_\mathrm{Q}`$. The excitation energy $`\omega _\nu `$ of the state $`|0^+`$, and the amplitudes $`X_j^{(\nu )}`$ and $`Y_j^{(\nu )}`$ are found respectively as the eigenenergy and the components of the eigenvector of the QRPA equation, which is derived from the following equation of motion for the Hamiltonian (7):

$${}_{Q}{}^{}0|[\delta Q,[H,Q_\nu ^{}]]|0_{\mathrm{Q}}^{}=\omega _\nu {}_{Q}{}^{}0|[\delta Q,Q_\nu ^{}]|0_{\mathrm{Q}}^{}.$$
(35)
In the standard way of derivation of the QRPA equations, the BCS equation is solved first. Then the a- and b-terms in the Hamiltonian (7) is replaced with the BCS result, which is $`H_{\mathrm{BCS}}=_jE_j𝒩_j`$. Using the exact commutation relations (9), we see that, among the remaining terms of (7), which do not contribute in the BCS, the d-, h-, and q-terms start to contribute within the QRPA. The c-term and g-term do not contribute since, in the commutation with the phonon operators (34), the former gives a number, while the latter leads to the terms $`𝒜_j^{}𝒜_j^{}^{}`$, $`𝒜_j^{}𝒜_j^{}`$, and $`𝒩_j(1𝒩_j^{}/\mathrm{\Omega }_j^{})`$, which are left out by linearizing the equation of motion according to (35). Moreover, in order to obtain a set of QRPA equations, linear with respect to the $`X_j^{(\nu )}`$ and $`Y_j^{(\nu )}`$ amplitudes, another approximation called the quasiboson approximation is made, which implies that the following approximate commutation relation holds

$$[𝒜_j,𝒜_j^{}^{}]=\delta _{jj^{}},$$
(36)
instead of Eq. (9). The definition of phonon operators (34) and the quasiboson approximation (36) lead to the well-known normalization of the QRPA $`X_j^{(\nu )}`$ and $`Y_j^{(\nu )}`$ amplitudes

$$\underset{j}{}[X_j^{(\nu )}X_j^{(\nu ^{})}Y_j^{(\nu )}Y_j^{(\nu ^{})}]=\delta _{\nu \nu ^{}}$$
(37)
so that the phonon operators are bosons, i.e.:

$$[Q_\nu ,Q_\nu ^{}^{}]=\delta _{\nu \nu ^{}}$$
(38)
The quasiboson approximation (36) shows that the quasiparticle-pair operators $`𝒜_j^{}`$ and $`𝒜_j`$ behave like boson operators when interacting with each other. The effect of Pauli principle represented by the last term at the RHS of (9) is just ignored. The set of QRPA equations obtained in this way is written in the matrix form as

$$\left(\genfrac{}{}{0pt}{}{AB}{BA}\right)\left(\genfrac{}{}{0pt}{}{X}{Y}\right)=\omega \left(\genfrac{}{}{0pt}{}{X}{Y}\right).$$
(39)
The explicit form of the matrices $`A`$ and $`B`$ depend on the approximation. Below we compare the results obtained within the boson and fermion formalisms.

#### 1 Boson formalism

The boson formalism is based on two following assumptions:

(a) It considers $`𝒜_j^{}`$ and $`𝒜_j^{}`$ as ideal bosons $`𝐛_j^{}`$ and $`𝐛_j`$, respectively, according to the quasiboson approximation (36). The $`j`$-shell quasiparticle number operator $`𝒩_j`$ is then mapped onto a boson pair as

$$𝒩_j=2𝐛_j^{}𝐛_j.$$
(40)
This mapping preserves the commutators (10). The Hamiltonian (7) can be then fully expressed in terms of the boson operators $`𝐛_j^{}`$ and $`𝐛_j`$.

(b) In deriving the QRPA equations according to (35), the q-term of the Hamiltonian (7) is neglected. This term is a special case of the so-called scattering term in the general Hamiltonian with a two-body residual interaction. For instance, when the residual interaction is separable, the q-term involves a sum of products of two scattering quasiparticle pairs $`B_{jj^{}}^{}(LM)B_{j_1j_1^{}}(L\stackrel{~}{M})`$, where $`B_{jj^{}}(LM)=_{mm^{}}jmj^{}m^{}|LM\alpha _{jm}^{}\alpha _{j^{}\stackrel{~}{m^{}}}`$ is the scattering quasiparticle-pair operator. The latter is equal to $`𝒩_j/\sqrt{2j+1}`$ when $`L=0`$. The scattering term is usually omitted in most of numerical calculations within QRPA for realistic nuclei in literature.

The boson mapping of the phonon operator (34) becomes

$$Q_\nu ^{}\underset{j}{}(X_j^{(\nu )}𝐛_j^{}Y_j^{(\nu )}𝐛_j).$$
(41)
The QRPA matrices $`A_{jj^{}}`$ and $`B_{jj^{}}`$ have the simple form:

$$A_{jj^{}}=2E_j\delta _{jj^{}}+d_{jj^{}},B_{jj^{}}=2h_{jj^{}},$$
(42)
which, in the present two-level model, become

$$A_{11}=A_{22}=G\mathrm{\Omega }+\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }},$$
(43)
$$A_{12}=A_{21}=\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }},$$
(44)
$$B_{11}=B_{22}=\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }},$$
(45)
$$B_{12}=B_{21}=G\mathrm{\Omega }\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }}.$$
(46)
The QRPA equations give one positive solution equal to

$$\omega _{\mathrm{QRPA}}^{(b)}=2\mathrm{\Delta },$$
(47)
while the spurious state associated with the non-conservation of particle number is located exactly at zero energy. We see that the energy of the first excited $`0^+`$ state above the phonon ground state $`|0_\mathrm{Q}`$ is just twice the pairing gap, the same for the lowest two-quasiparticle excitation above the quasiparticle ground state $`|0_\alpha `$ of a system with an even particle number.

In the normal-fluid phase ($`\mathrm{\Delta }`$=0), the quasiparticle operator $`\alpha _{jm}`$ becomes the particle-($`p`$) creation or hole-($`h`$) destruction operator depending on whether the level is located above or below the Fermi level, namely

$$\alpha _{jm}^{}=\{\begin{array}{cc}a_{jm}^{}\hfill & \text{if }ϵ_j>\lambda \hfill \\ a_{jm}\hfill & \text{if }ϵ_j<\lambda \hfill \end{array}$$
(48)
Therefore, the following boson mapping for operators $`A_j^{}`$ and $`A_j`$ holds

$$A_j^{}\{\begin{array}{cc}𝐛_j^{}\hfill & \text{if }ϵ_j>\lambda \hfill \\ 𝐛_j\hfill & \text{if }ϵ_j<\lambda \hfill \end{array},A_j\{\begin{array}{cc}𝐛_j\hfill & \text{if }ϵ_j>\lambda \hfill \\ 𝐛_j^{}\hfill & \text{if }ϵ_j<\lambda \hfill \end{array},$$
(49)
The boson mapping for the number operator $`N_j`$ is given as

$$N_j=\{\begin{array}{cc}2𝐛_j^{}𝐛_j\hfill & \text{if }ϵ_j>\lambda \hfill \\ 2(\mathrm{\Omega }_j𝐛_j^{}𝐛_j)\hfill & \text{if }ϵ_j<\lambda \hfill \end{array},$$
(50)
which preserves the commutation relations (5) and the particle number on the j-shell equal to $`2\mathrm{\Omega }_j`$. The RPA matrices for the two-level model are

$$A_{11}=A_{22}=2(ϵG\mathrm{\Omega }),A_{12}=A_{21}=0.$$
(51)
$$B_{11}=B_{22}=0,B_{12}=B_{21}=2G\mathrm{\Omega }.$$
(52)
They lead to the RPA phonon energy

$$\omega _{\mathrm{RPA}}^{(b)}=2ϵ\sqrt{1\frac{G}{G_{\mathrm{cr}}^{(b)}}},G_{\mathrm{cr}}^{(b)}=\frac{ϵ}{2\mathrm{\Omega }}.$$
(53)
The critical value of $`G_{\mathrm{cr}}^{(b)}`$ where the RPA collapses is the same collapsing point of the BCS obtained without the self-energy term (Eqs. (23) and (24)). Therefore, in order to have the superfluid regime start at the same critical point, the gap $`\mathrm{\Delta }`$ in Eqs. (43) – (47) should be calculated neglecting the self-energy term in the BCS equation (21). This gives

$$\mathrm{\Delta }^{(b)}=G\mathrm{\Omega }\sqrt{1\left(\frac{ϵ}{2G\mathrm{\Omega }}\right)^2}.$$
(54)
The phonon energies (47) and (53) are exactly those obtained for the first time in Ref.  using the space-variable technique and the gap $`\mathrm{\Delta }^{(b)}`$ (54).

#### 2 Fermion formalism

The fermion formalism does not use directly the quasiboson approximation (36). Instead, it employs the exact commutation relations (9) and (10) to rearrange the results of calculating the commutators $`[H,𝒜_j^{}]`$ and $`[H,𝒜_j]`$ into the normal order. Then, in the process of linearizing the equation of motion (35) the following “average” quasiboson approximation is used

$$_\mathrm{Q}0|[𝒜_j,𝒜_j^{}^{}]|0_\mathrm{Q}=\delta _{jj^{}},$$
(55)
assuming that the quasiparticle occupation number $`n_j`$ in the RPA ground state is zero, i.e.

$$n_j\frac{{}_{\mathrm{Q}}{}^{}0|𝒩_j|0_{\mathrm{Q}}^{}}{2\mathrm{\Omega }_j}=0.$$
(56)
The QRPA matrices are obtained in this way as

$$A_{jj^{}}=2(E_j+2q_{jj})\delta _{jj^{}}+d_{jj^{}},B_{jj^{}}=2(1\frac{1}{\mathrm{\Omega }_j}\delta _{jj^{}})h_{jj^{}}.$$
(57)
Their explicit form in the two-level model is given as

$$A_{11}=A_{22}=G\mathrm{\Omega }\frac{\mathrm{\Delta }^2}{G\mathrm{\Omega }^2}+\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }},$$
(58)
$$A_{12}=A_{21}=\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }},$$
(59)
$$B_{11}=B_{22}=\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }^2}+\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }},$$
(60)
$$B_{12}=B_{21}=G\mathrm{\Omega }\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }}.$$
(61)
The BCS equation (21) for the gap $`\mathrm{\Delta }`$ includes the self-energy term. The $`4q_{jj}\delta _{jj^{}}`$-term in the expression for $`A_{jj}`$ (Eqs. (57) and (58)) appears due to the use of the exact commutation relation (10) to calculate the commutator between the q-term of the Hamiltonian (7) and $`𝒜_j^{}`$ (or $`𝒜_j`$) as follows:

$$\underset{j_1j_1^{}}{}q_{j_1j_1^{}}[𝒩_{j_1}𝒩_{j_1^{}},𝒜_j^{}]=4q_{jj}𝒜_j^{}+4\underset{j^{}}{}q_{jj^{}}𝒜_j^{}𝒩_j^{}.$$
(62)
The last term at the RHS of (62) does not contribute in the linearization within the QRPA. The first term leads to the above-mentioned 4$`q_{jj}\delta _{jj^{}}`$-term. It is worthwhile to notice that, if one rearranges the quasiparticle operators in the q-term of (7) to the normal order as

$$\underset{jj^{}}{}q_{jj^{}}𝒩_j𝒩_j^{}=\underset{j}{}q_{jj}𝒩_j\underset{jmj^{}m^{}}{}q_{jj^{}}\alpha _{jm}^{}\alpha _{j^{}m^{}}^{}\alpha _{jm}\alpha _{j^{}m^{}},$$
(63)
and then drops the last term at the RHS (although there is no rigorous justification for doing so), the remaining term in the commutation with $`𝒜_j^{}`$ gives

$$\underset{j_1j_1^{}}{}q_{j_1j_1^{}}[𝒩_{j_1}𝒩_{j_1^{}},𝒜_j^{}]\underset{j^{}}{}q_{j^{}j^{}}[𝒩_j^{},𝒜_j^{}]=2q_{jj}𝒜_j^{},$$
(64)
i.e. just the half of the first term at the RHS of Eq. (62). This explains the difference between Eq. (58) above and Eq. (38) of Ref.  by a factor of 2 in the denominator of the second term at their RHS.

The energy of the first $`0^+`$ excited state $`\omega `$ obtained solving the QRPA equation with the matrices (58) – (61) has the form

$$\omega _1=2\mathrm{\Delta }\sqrt{\left(1\frac{3}{4\mathrm{\Omega }}\right)\left(1\frac{\mathrm{\Delta }^2}{4G^2\mathrm{\Omega }^3}\right)}.$$
(65)
The energy of the spurious state in this case is imaginary. If one uses the approximation (64) instead of (62), the positive solution is given by Ref.  as:

$$\omega _{01}=2\mathrm{\Delta }\sqrt{1\frac{1}{2\mathrm{\Omega }}},$$
(66)
while the energy of the spurious state is exactly zero. Neglecting the q-term in (7) yields the energy of the first $`0^+`$ state as

$$\omega _0=2\mathrm{\Delta }\sqrt{\left(1\frac{1}{4\mathrm{\Omega }}\right)\left(1+\frac{\mathrm{\Delta }^2}{4G^2\mathrm{\Omega }^3}\right)}.$$
(67)
There are two features of the fermion formalism, which are worth noticing here. The fact that in the process of linearizing the equation of motion, the exact commutation relations (9) and (10) have been used to calculate the commutators $`[H,𝒜_j^{}]`$ and $`[H,𝒜_j]`$ does not reduce matrices (57) by setting $`q_{jj^{}}=`$ 0 to those given within the boson formalism (42). The difference in $`B_{jj}`$ still remains. Another feature is that, when the q-term in (7) is omitted, the spurious mode in the fermion formalism is shifted to a positive value equal to

$$\omega _0^s=\mathrm{\Delta }\sqrt{\frac{1}{\mathrm{\Omega }}\left(1\frac{\mathrm{\Delta }^2}{4G^2\mathrm{\Omega }^3}\right)}.$$
(68)
However it is much smaller than the pairing gap $`\mathrm{\Delta }`$ especially at large $`\mathrm{\Omega }`$, therefore, remains well-isolated. Comparing Eqs. (65) – (67) with the boson energy $`\omega _{\mathrm{QRPA}}^{(b)}`$ (Eq. (47)), it is easy to see that

$$\omega _1<\omega _{01}<\omega _0<\omega _{\mathrm{QRPA}}^{(b)}=2\mathrm{\Delta }.$$
(69)
The ground-state energy is now calculated following Refs.  as

$$E_{\mathrm{QRPA}}=E_{\mathrm{BCS}}+\frac{1}{2}[\omega (A_{11}+A_{22})].$$
(70)
In the limit of zero gap, the $`pp`$-RPA equation is obtained from the QRPA equation discussed above putting $`\mathrm{\Delta }=`$ 0, $`v=v_1=u_2=`$1, $`u=u_1=v_2=0`$. The solution of this RPA equation is decoupled into the addition and removal modes, which have been discussed in detail in Ref. . In the present two-level model, these two sets of equations can be written in one matrix equation as

$$\left(\genfrac{}{}{0pt}{}{A_{11}B_{12}}{B_{12}A_{22}}\right)\left(\genfrac{}{}{0pt}{}{R_p^{(\tau )}}{R_h^{(\tau )}}\right)=\left(\genfrac{}{}{0pt}{}{10}{01}\right)\left(\genfrac{}{}{0pt}{}{R_p^{(\tau )}}{R_h^{(\tau )}}\right)\omega _\tau ,(\tau =1,2),$$
(71)
where

$$A_{11}=ϵ+2\lambda +2GG\mathrm{\Omega },A_{22}=ϵ2\lambda G\mathrm{\Omega },B_{12}=G\mathrm{\Omega },$$
(72)
and $`R_p^{(1)}=X_\mathrm{a}`$ , $`R_h^{(1)}=Y_\mathrm{a}`$ , $`R_p^{(2)}=Y_\mathrm{r}`$ , $`R_h^{(2)}=X_\mathrm{r}`$ , $`\omega _1=\omega _\mathrm{a}`$ , and $`\omega _2=\omega _\mathrm{r}`$ . The energies $`\omega _\mathrm{a}`$ and $`\omega _\mathrm{r}`$ are found as

$$\omega _\mathrm{a}=G2\lambda \pm \sqrt{(ϵ+G)(ϵ+G2G\mathrm{\Omega })}.$$
(73)
$$\omega _\mathrm{r}=G+2\lambda \pm \sqrt{(ϵ+G)(ϵ+G2G\mathrm{\Omega })},$$
(74)
where only the sign $`+`$ in front of the square roots corresponds to the positive values for these energies. Here we still keep the factor $`\pm 2\lambda `$, which is useful to connect the RPA solutions with the QRPA one at the critical point $`G=G_{\mathrm{cr}}`$, where 2$`\lambda =G_{\mathrm{cr}}`$ according to Eq. (20).

In order make a comparison with the boson formalism, where there is only one boson state, we apply the sum-rule method, representing the RPA phonon operator as

$$Q^{}=Q_\mathrm{a}^{}+Q_\mathrm{r}^{}.$$
(75)
The RPA ground-state wave function $`|\mathrm{RPA}`$ can be written as a direct product of the ground-state wave functions of the (orthogonal) additional and removal modes

$$|\mathrm{RPA}=|0_\mathrm{a}|0_\mathrm{r}.$$
(76)
Using the Thouless theorem for the energy-weighted sum rule $`𝒮_1`$  with respect to the phonon operator $`Q^{}`$ (75)

$$𝒮_1\underset{i}{}\omega _{\mathrm{RPA}}^{(i)}|\nu _i|Q^{}|\mathrm{RPA}|^2=\frac{1}{2}\mathrm{HF}|[Q^{},[H,Q^{}]]|\mathrm{HF},$$
(77)
it is easy to see that the LHS of Eq. (77) is equal $`\omega _{\mathrm{RPA}}`$ since $`|\nu _i=Q^{}|\mathrm{RPA}`$ ($`i=`$ 1), while the RHS is equal to

$`{\displaystyle \frac{1}{2}}\left[\mathrm{HF}|[Q_\mathrm{a}^{},[H,Q_\mathrm{a}^{}]]|\mathrm{HF}+\mathrm{HF}|[Q_\mathrm{r}^{},[H,Q_\mathrm{r}^{}]]|\mathrm{HF}\right]=`$

$$\omega _\mathrm{a}|\nu _\mathrm{a}|Q_\mathrm{a}^{}|\mathrm{RPA}|^2+\omega _\mathrm{r}|\nu _\mathrm{r}|Q_\mathrm{r}^{}|\mathrm{RPA}|^2=\omega _\mathrm{a}|\nu _\mathrm{a}|\nu _\mathrm{a}|^2+\omega _\mathrm{r}|\nu _\mathrm{r}|\nu _\mathrm{r}|^2=\omega _\mathrm{a}+\omega _\mathrm{r}.$$
(78)
(The crossing terms $`\mathrm{HF}|[Q_\mathrm{i}^{},[H,Q_\mathrm{i}^{}^{}]]|\mathrm{HF}`$ ($`ii^{}`$) vanish as can be easily checked using $`i|i^{}=\delta _{ii^{}}`$). Therefore

$$\omega _{\mathrm{RPA}}=\omega _\mathrm{a}+\omega _\mathrm{r}.$$
(79)
The ground-state energy is then given by

$$E_{\mathrm{RPA}}=E_{\mathrm{HF}}+\frac{1}{2}(\omega _{\mathrm{RPA}}A_{11}A_{22})=\sqrt{(ϵ+G)(ϵ+G2G\mathrm{\Omega })}(ϵG\mathrm{\Omega }+G),$$
(80)
which is exactly the same expression obtained previously in Eq. (45) of Ref. .

The energies of the ground-state and the first $`0^+`$ state obtained within the fermion formalism of RPA and QRPA for $`\mathrm{\Omega }=`$ 4 are compared with the exact energies in Fig. 2. The figure shows that, in the superfluid regime, except for the region close to the critical point $`G=G_{\mathrm{cr}}`$ where the QRPA collapses, the approximation (64) fits the exact result for the ground-state energy better than (62). For the energy of the first $`0^+`$ state, both of the QRPA versions, which include the q-term of the Hamiltonian (7) give practically the same values. However, they are obviously smaller compared to the exact solution. This discrepancy increases with increasing the interaction. From Eqs. (65) – (67) we found that, in the limit of infinite $`G`$ the solution $`\omega _0`$ (q-term neglected) becomes $`2\mathrm{\Delta }\sqrt{11/4\mathrm{\Omega }}3.88\mathrm{\Delta }`$, while the solutions $`\omega _1`$ based on the approximation (62), and $`\omega _{10}`$ based on (64 are equal to $`0.96\omega _0`$ and $`0.93\omega _0`$, respectively. The only QRPA approximation that fits well the exact results for both of the ground-state and excited-state energies the one, which neglects the q-term of the Hamiltonian (7). The results within this approximation practically coincide with the exact ones at large $`G`$. The RPA energy $`\omega _{\mathrm{RPA}}=\omega _\mathrm{a}+\omega _\mathrm{r}`$ of the first $`0^+`$ state exhibits the well-known behavior. It decreases with increasing $`G`$ and collapses at the same critical point $`G=G_{\mathrm{cr}}`$, from which the the normal-fluid phase ceases to exist, and the superfluid phase begins. Meanwhile, the exact solution for the first excited state has only a bending in this region, showing no signature of such phase transition. For the ground-state energy the exact result shows a completely smooth curve, while the critical point is clearly seen in the approximations.

As shown in Fig. 3 for $`\mathrm{\Omega }=`$ 4, the QRPA energy $`\omega _0`$ obtained without the q-term of (7), and the RPA energy $`\omega _{\mathrm{RPA}}`$ also match well the solutions of the boson formalism, especially after shifting $`G_{\mathrm{cr}}=G_{\mathrm{cr}}^{(b)}`$ in the latter to the value of $`ϵ/(2\mathrm{\Omega }1)`$ used in the fermion formalism due to the inclusion of the self-energy term (Eq. (24)).

### B Renormalized QRPA

The collapse of the BCS approximation and QRPA (RPA) has the same origin of neglecting the Pauli principle between quasiparticle pairs operators in the BCS approximation (26) and the quasiboson approximation strictly (36) or in average (55). The LN method approximately corrects this inconsistency within the BCS approximation. For the QRPA this is done by the renormalized QRPA.

The essence of the renormalized QRPA is to replace the quasiboson approximation in the form of Eqs. (36) or (55) with the average value of the commutator

$$\stackrel{~}{0}|[𝒜_j,𝒜_j^{}^{}]|\stackrel{~}{0}=D_j\delta _{jj^{}},D_j=12n_j,$$
(81)
in a new ground state $`|\stackrel{~}{0}`$, where the correlations beyond the QRPA due to the fermion structure of the quasiparticle pairs $`𝒜_j^{}`$ and $`𝒜_j`$ are taken into account, namely

$$n_j=\frac{1}{2\mathrm{\Omega }_j}\stackrel{~}{0}|𝒩_j|\stackrel{~}{0}0,$$
(82)
instead of the assumption (56). The phonon operators are renormalized as

$$𝒬_\nu =\underset{j}{}\frac{1}{\sqrt{D_j}}(𝒳_j^{(\nu )}𝒜_j^{}𝒴_j^{(\nu )}𝒜_j),𝒬_\nu =(𝒬_\nu ^{})^{},$$
(83)
so that the condition for phonons to be bosons within the correlated ground state $`|\stackrel{~}{0}`$

$$\stackrel{~}{0}|[𝒬_\nu ,𝒬_\nu ^{}^{}]|\stackrel{~}{0}=\delta _{\nu \nu ^{}}$$
(84)
leads to the same normalization condition for the amplitudes $`𝒳_j^{(\nu )}`$ and $`𝒴_j^{(\nu )}`$ as that of the QRPA, i.e. $`_j(𝒳_j^{(\nu )}𝒳_j^{(\nu ^{})}𝒴_j^{(\nu )}𝒴_j^{(\nu ^{})})=\delta _{\nu \nu ^{}}.`$ The factor $`D_j`$ is calculated according to the approximation in Ref.  as

$$D_j=\frac{1}{1+(𝒴_j^{(\nu )})^2/\mathrm{\Omega }_j}.$$
(85)
The renormalized-QRPA matrices $`A_{jj^{}}`$ and $`B_{jj^{}}`$ are given as

$$A_{jj^{}}=2(E_j+2q_{jj})\delta _{jj^{}}+4\underset{j\mathrm{"}}{}\mathrm{\Omega }_{j^{\prime \prime }}q_{j^{}j^{\prime \prime }}(1D_{j^{\prime \prime }})+D_jd_{jj^{}},$$
(86)
$$B_{jj^{}}=2(D_j\frac{1}{\mathrm{\Omega }_j}\delta _{jj^{}})h_{jj^{}}.$$
(87)
Since only the omission of the q-term in (7) within the QRPA reproduces well the exact solution in the present two-level model at large values of $`G`$, we discuss below only the case when $`q_{jj^{}}=`$0. The renormalized QRPA matrices in this case read

$$A_{11}=A_{22}=G\mathrm{\Omega }(2D)+\frac{D\mathrm{\Delta }^2}{2G\mathrm{\Omega }^2},A_{12}=A_{21}=\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }}(43D),$$
(88)
$$B_{11}=B_{22}=\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }^2}+\frac{D\mathrm{\Delta }^2}{2G\mathrm{\Omega }},B_{12}=B_{21}=D(G\mathrm{\Omega }\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }}).$$
(89)
The energy of the first excited $`0^+`$ state is found as

$$\omega _{\mathrm{RQRPA}}=2\mathrm{\Delta }\sqrt{\left[\left(D\frac{1}{4\mathrm{\Omega }}\right)+\frac{G^2\mathrm{\Omega }^2(1D)}{\mathrm{\Delta }^2}\right]\left(1+\frac{\mathrm{\Delta }^2}{4G^2\mathrm{\Omega }^3}\right)},$$
(90)
By setting $`D_j=`$ 1 ($`n_j=`$ 0) in Eqs. (88) – (90), the QRPA limit in Eqs. (58) – (61) (without the second term at the RHS of (58) since the q-term is neglected), and $`\omega _0`$ given by Eq. (67) are recovered . The ground-state energy is calculated using Eq. (70), using the energy $`\omega _{\mathrm{RQRPA}}`$ given by Eq. (90) and the matrices $`A_{11}=A_{22}`$ given by Eq. (88).

The renormalized RPA matrices ($`G<G_{\mathrm{cr}}`$) are given as

$$A_{11}=ϵ+2GDG\mathrm{\Omega },A_{22}=ϵDG\mathrm{\Omega },A_{12}=A_{21}=0,$$
(91)
$$B_{11}=B_{22}=0,B_{12}=B_{21}=DG\mathrm{\Omega }.$$
(92)
The positive renormalized RPA phonon energies are found as

$$\stackrel{~}{\omega }_\mathrm{a}=G2\lambda +\sqrt{(ϵ+G)(ϵ+G2DG\mathrm{\Omega })}.$$
(93)
$$\stackrel{~}{\omega }_\mathrm{r}=G+2\lambda +\sqrt{(ϵ+G)(ϵ+G2DG\mathrm{\Omega })},$$
(94)
$$\omega _{\mathrm{RRPA}}=\stackrel{~}{\omega }_\mathrm{a}+\stackrel{~}{\omega }_\mathrm{r}=2\sqrt{(ϵ+G)(ϵ+G2DG\mathrm{\Omega })},$$
(95)
and the ground-state energy is

$$E_{\mathrm{RRPA}}=\sqrt{(ϵ+G)(ϵ+G2DG\mathrm{\Omega })}(ϵ+GDG\mathrm{\Omega }).$$
(96)
As has been shown in Ref. , the presence of the factor $`D`$ in the renormalized RPA matrices reduces the interaction in such a way that the critical point where the RPA collapses is completely washed out. The energy $`\stackrel{~}{\omega }_\mathrm{a}`$, $`\stackrel{~}{\omega }_\mathrm{r}`$, and consequently, $`\omega _{\mathrm{RRPA}}`$ are always real. As for the renormalized QRPA, the presence of the pairing gap $`\mathrm{\Delta }`$ makes it still collapse if $`\mathrm{\Delta }`$ is calculated within the BCS (Eq. (21), but it is no longer the case if the LN pairing gap $`\mathrm{\Delta }_{\mathrm{LN}}`$ (28) is used.

### C LN-method + renormalized RPA (renormalized QRPA)

We have seen in the preceding sections that the LN-method allows us to avoid the phase transition of the BCS approximation, the QRPA without the q-term in (7) provides us with the best fit of the exact results for the energies of the ground state and the first excited $`0^+`$ state, while the renormalized RPA is known to smear out the collapse of the RPA due to the violation of the Pauli principle within the quasiboson approximation . Therefore, in order to avoid the phase transition and to give at the same time a good description for both of the energies of the ground state as well as the fist excited $`0^{}`$ state, we propose here the following recipe.

i) The QRPA equations (without the q-term in (7)) are solved using the pairing gap $`\mathrm{\Delta }_{\mathrm{LN}}`$ (28) found in the LN method.

ii) The ground state energy is calculated using Eq. (70), in which the LN pairing gap $`\mathrm{\Delta }_{\mathrm{LN}}`$ is used instead of the BCS gap $`\mathrm{\Delta }`$ (21) and the renormalized QRPA matrices given by Eq. (88) are used instead of Eq. (58).

iii) The first excited $`0^+`$ state is presented as a mixed state of those obtained within the renormalized RPA and renormalized QRPA (in (i)). Applying the sum-rule method and following a derivation similar to Eqs. (76) – (78), we find that the energy $`\omega _{0^+}`$ of this mixed state is the sum of the renormalized RPA and renormalized QRPA energies, namely

$$\omega _{0^+}=\omega _{\mathrm{RRPA}}+\omega _{\mathrm{RQRPA}},$$
(97)
where $`\omega _{\mathrm{RRPA}}`$ and $`\omega _{\mathrm{QRPA}}`$ are given by Eqs. (95) and (90), respectively, with $`\mathrm{\Delta }_{\mathrm{LN}}`$ used in the latter. The results obtained using this recipe will be referred to as the combined results.

The combined results for ground state energy and the energy of the first excited $`0^+`$ state are presented as thin solid lines in Fig. 4 at several values of $`\mathrm{\Omega }`$ in comparison with the exact results and those of RPA, QRPA and their renormalized versions. The result obtained within the renormalized RPA and the combined result (i) for the energy of the first excited $`0^+`$ state show that the superfluid-normal phase transition is completely washed out. For the ground-state energy, the combined result (ii) is closer to the exact one compared to the LN result in the weak-coupling region ($`GG_{\mathrm{cr}}`$) and in the region around $`G_{\mathrm{cr}}`$. It coincides with the LN and exact results in the strong-coupling region ($`GG_{\mathrm{cr}}`$) (See the thin lines in Fig. 4 (a) – (c)). The combined result (iii) agrees reasonable well with the exact results (See the thin lines in Fig. 4 (d) – (f)), especially for larger $`\mathrm{\Omega }`$. Hence, this numerical test confirms the validity of the assumption (iii) that the exact excited $`0^+`$ state may be considered as a mixed state where both of the normal and superfluid phases coexist at $`G`$ 0. The weak-coupling region is dominated by the normal-fluid phase, where the energy $`\omega _{0^+}`$ of the first excited $`0^+`$ state decreases monotonously with increasing the pairing interaction $`G`$. The strong-coupling region is dominated by the superfluid phase, in which $`\omega _{0^+}`$ increases with increasing $`G`$.

## V Conclusions

In this work, a schematic model of two symmetric levels interacting via a pairing force has been used to test several well-known approximations by calculating the ground-state energy and the energy $`\omega _{0^+}`$ of the first excited $`0^+`$ state. Results obtained within the BCS approximation, LN method, boson and fermion formalisms for RPA and QRPA, renormalized RPA, and renormalized QRPA have been analyzed in details. It is shown that the common version of the QRPA, which neglects the scattering term (the q-term of the Hamiltonian (7)) and the boson formalism give the closest results to the exact ones for both of the ground-state energy and $`\omega _{0^+}`$ at the values of the pairing interaction parameter $`GG_{\mathrm{cr}}`$. Meanwhile the QRPA version used in Ref.  as well as the QRPA using the exact commutation relations between quasiparticle-pair operators to treat the scattering term do not describe well the energy of the first excited $`0^+`$ state at large $`G`$.

A recipe has been proposed which combines the results obtained within the LN method for the pairing gap $`\mathrm{\Delta }_{\mathrm{LN}}`$, the renormalized RPA, and the renormalized QRPA neglecting the last term (q-term) of the Hamiltonian (7). The combined results agree reasonably well with the exact ones for both of the ground-state energy and $`\omega _{0^+}`$, showing no signature of a sharp superfluid-normal phase transition. The agreement is better at a larger particle number. The results suggest that the exact excited $`0^+`$ state can be decomposed approximately into two components, which correspond to the normal and superfluid regimes, respectively. The weak-coupling region is dominated by the normal-fluid phase, while the strong-coupling region - by the superfluid phase. Since the proposed scheme is based on rather simple, well-known, and numerically accessible approximations, its future extension toward an application in realistic nuclei may be useful.

For a more self-consistent approach, one can derive the set of RPA equations using the Hamiltonian (7) with the coefficients $`u_j`$, $`v_j`$, $`X_j^{(\nu )}`$ and $`Y_j^{(\nu )}`$ left as variational parameters when minimizing the average energies over the ground state and the first excited $`0^+`$ state. This may serve as a goal of the future study.

## Acknowledgments

Thanks are due to V. Zelevinsky for reading the manuscript, valuable comments and suggestions.


text:
# Self-gravitating discs: what can we learn from the dynamics of maser spots ?

## 1. The inverse problem

Finding the spatial distribution of matter $`\rho (\stackrel{}{r})`$ that reproduces a given rotation curve $`\mathrm{\Omega }_{\mathrm{obs}.}(\stackrel{}{r})`$ is a common inverse problem in Astrophysics. It can be addressed at the scale of AGN gaseous discs if one interprets the rotation curve of water masers detected in the core of some objects (Greenhill, 2002) in terms of a pure gravitational attraction. Efforts must therefore be made in developing reliable theoretical and/or numerical techniques to perform the inversion (e.g. Sibgatullin, Garcia & Manko, 2002). In the present situation, this must account for long-range interacting components (black hole, stellar cluster, bars, galactic bludge, etc.) as well as for objects populating the parsec-scale (disc, clouds, compact objects, torus, etc..) where masers are seen. Uncertainties are large, limiting our predictions and interpretations. We discuss here an example of inversion by considering a very simplified, two component system, made of a black hole and a disc, with an application to NGC 4258 and NGC 1068.

## 2. Assumptions, equations and method

At the parsec-scale, the standard theory predicts that all discs are (strongly) self-gravitating (Collin & Huré, 2000) and so, they might evolve towards a state with a unity $`Q`$-Toomre parameter. Accounting for disc mass effects is therefore necessary to investigate these outer regions. But a reliable modeling of such discs is not permitted yet, given difficulties of various origins (nature of the flow, turbulent transport, numerical resolution, instabilities, boundary conditions), and some assumptions are necessary. The simplest way to tackle the question is a “minimal approach” which assumes that the velocity field has a pure gravitational origin. For a black hole with mass $`M_{\mathrm{BH}}`$ and a disc only, the most basic equation set to solve is

$$\{\begin{array}{cc}& \mathrm{\Omega }_{\mathrm{obs}.}^2=\frac{GM_{\mathrm{BH}}}{r^3}\frac{g_R^{\mathrm{disc}}}{R}\hfill \\ & .𝐠^{\mathrm{disc}}=4\pi G\rho ,\hfill \end{array}$$

where $`𝐠^{\mathrm{disc}}`$ is the gravitational acceleration due to the disc, $`R`$ is the cylindrical radius, $`r`$ is the spherical radius and $`\rho `$ is the unknown density, a tri-dimensional function. In other words, if we are able to compute the self-gravitating field for discs with various size and mass distributions, we can select among possible solutions those matching observations. Note that, in this kind of approach, no physical model for the disc imposes the $`\rho `$function.

Here, the Poisson equation is solved numerically assuming azimuthal invariance for a given $`\rho (\stackrel{}{r})`$ from the integral form

$$𝐠^{\mathrm{disc}}(\stackrel{}{r})=G_{\mathrm{vol}}\frac{\rho (\stackrel{}{r})(\stackrel{}{r}\stackrel{}{r}^{})d^3r}{|\stackrel{}{r}\stackrel{}{r}^{}|^3},$$
(1)
by the method described in Huré (2002a). We stress that the full regularization of the volume integral requires three successive integrations, otherwise diverging fields (and subsequently infinite velocities) can result at the disc edges (e.g. Mestel 1963). This can be critical if the maser spots rotate near the disc outer edge. The present inversion method works with three inputs: the black hole mass $`M_{\mathrm{BH}}`$, the mass density in the disc $`\rho (R,z)`$, and its outer radius $`R_{\mathrm{out}}`$. As several solutions can lead, at least locally, to similar rotation curves, we impose that the total surface density $`\mathrm{\Sigma }=\rho 𝑑z`$ is of the form $`\mathrm{\Sigma }R^s`$, with a smooth connection to the ambient medium at the disc edges (see Huré 2002b for details). The mass density inside the disc obeys a parabolic law with the altitude, and the disc aspect ratio $`H/R`$ is held constant. These additional restrictions however do not remove the degeneracies. Once a model fits observations, inputs become outputs.

## 3. NGC 4258: the possible ambiguity

NGC 4258 is the most famous example of Seyfert-2 active galaxy showing a thin outer disc with masers in quasi-perfect keplerian motion (Myoshi et al. 1995). A priori, this sustains the simple idea that most of the mass is effectively gathered at the center, implying a black hole mass $`3.9\times 10^7`$ M. If this is correct, the inversion method can furnish an upper limit for the disc mass: as figure 1 shows, the mass of the masing disc at $`0.28`$ pc (or $`8`$ mas) must be less than $`10^5`$ M (that is a mass ratio $`q2.5\times 10^3`$), otherwise a deviation to the keplerian rotation results. But the possibility that a massive disc is present around a lighter black hole must not be ruled out. Actually, the rotation curve concerns a short spatial range and various $`R`$-laws can mimic a keplerian profile. Figure 1b shows a case where the rotation curve is still well reproduced with a black hole $`25\%`$ less massive and a disc with a mass $`8\times 10^6`$ M (or $`q0.26`$) and with an outer edge at $`0.11`$ pc (in this case, the masing clouds rotate “outside” the disc). In this example, we have deliberately selected a disc with $`Qcst`$ (for the reason mentioned before), but the upper limit derived in this way is almost not affected by this specific choice.

## 4. NGC 1068

In contrast, the rotation curve of the outer disc in NGC 1068 appears in noticeable departure with respect to the keplerian case. Greenhill et al. (1996) found $`\mathrm{\Omega }_{\mathrm{obs}.}R^{1.31}`$ in the range $`0.651.1`$ pc, and mentioned that the disc self-gravity could be the possible source of this non-keplerian behavior. Other mechanisms can be invoked, like radiation pressure effects (Pier & Krolik 1992) or gravitational attraction by a stellar cluster (Kumar 1999), but the disc self-gravity appears as a natural explanation. Applying the inversion method gives several families of solutions. The sensitivity to $`s`$ is important. Two possible solutions characterized by a mean angular rotation index of $`1.31`$ are shown in Fig. 2. The parameters deduced for the disc and central black hole are reported in Table 1. We find that the disc mass is $`10^7`$ M, comparable to that of the black hole. We stress that surface densities that do not vary approximately as $`R^1`$ do not match observations: with $`s>1`$, the global rotation curve is flatter (and steeper with $`s<1`$). The important point is that the value inferred for the $`Q`$parameter is of the order of unity, meaning that the outer disc could be in a marginally stable state. Moreover, the temperature and the density in the outer disc seem appropriate for excitation of water molecules (Huré 2002b).

## 5. Conclusions

We have discussed possible constraints (and uncertainties) on the mass of the central black hole and on structure of the parsec-scale region in NGC 4258 and NGC 1068, by a method involving a simple axi-symmetrical disc model and a black hole. The subsequent challenges will then be i) to understand, in terms of physical processes, why the disc has the inferred properties, ii) to check the the minimal approach, and iii) to discuss the implications of these solutions for the inner disc (see Lodato & Bertin in these proceedings).

### Acknowledgments.

I am grateful to C. Boisson, S. Collin, S. Courty, F. Hersant, F. Le Petit, D. Richard. and J.-P. Zahn for fruitful discussions.

## References

Collin S., Huré J.M., 2001, A&A, 372, 50

Greenhill L.J., Proceedings of IAU Symposium 206: Cosmic Masers (Mangaratiba, Brazil; March 2001)

Greenhill, L. J. et al., 1996, ApJ, 472, L21

Huré J.M., 2002a, submitted

Huré J.M., 2002b, A&A Let., in press

Kumar P., 1999, ApJ, 519, 599

Mestel L., 1963, ApJ, 126, 553

Miyoshi, M. et al. 1995, Nature, 373, 127

Pier, E. A., Krolik J. H. , 1992, ApJL, 399, L23-L26.

Sibgatullin, N.R., Garcia, G., Manko, V.S., 2002, astro-ph/0209200