text:
# Shear flow and Kelvin-Helmholtz instability in superfluids

\[

## Abstract

The first realization of instabilities in the shear flow between two superfluids is examined. The interface separating the A and B phases of superfluid 3He is magnetically stabilized. With uniform rotation we create a state with discontinuous tangential velocities at the interface, supported by the difference in quantized vorticity in the two phases. This state remains stable and nondissipative to high relative velocities, but finally undergoes an instability when an interfacial mode is excited and some vortices cross the phase boundary. The measured properties of the instability are consistent with a modified Kelvin-Helmholtz theory.

\]

Instabilities in the shear flow between two layers of fluids  belong to a class of interfacial hydrodynamics which is attributed to many natural phenomena. Examples are wave generation by wind blowing over water , the flapping of a sail or flag in the wind , and even flow in granular beds . In the hydrodynamics of inviscid and incompressible fluids the transition from calm to wavy interfaces is known as the Kelvin-Helmholtz (KH) instability . Since Lord Kelvin’s treatise in 1871, difficulties have plagued its description in ordinary fluids, which are viscous and dissipative. They also display a shear-flow instability, but its correspondence with that in the ideal limit is not straightforward. The tangential velocity discontinuity in the shear-flow instability is created by a vortex sheet. In a viscous fluid a planar vortex sheet is not a stable equilibrium state and not a solution of the hydrodynamic equations .

Superfluids provide a close variation of the ideal inviscid limit considered by Lord Kelvin and thus an environment where the KH theory can be tested. The initial state is a non-dissipative vortex sheet – the interface between two superfluids brought into a state of relative shear flow. So far the only experimentally accessible case where this can be studied in stationary conditions, is the interface between 3He-A and 3He-B , where the order parameter changes symmetry and magnitude, but is continuous on the scale of the superfluid coherence length $`\xi 10`$nm. We discuss an experiment, where the two phases slide with respect to each other in a rotating cryostat: 3He-A performs solid-body-like rotation while 3He-B is in the vortex-free state and thus stationary in the laboratory frame. While increasing the rotation velocity $`\mathrm{\Omega }`$, we record the events when the AB phase boundary becomes unstable – when some circulation from the A-phase crosses the AB interface and vortex lines are introduced into the initially vortex-free B phase. On increasing the rotation further, the instability occurs repeatedly. Such a succession of instability events can be understood as a spin-up of 3He-B by rotating 3He-A.

Our experimental setup is shown in Fig. 1. The AB boundary is forced against a magnetic barrier in a smooth-walled quartz container, by cooling the sample below $`T_{\mathrm{AB}}`$ at constant pressure in a rotating refrigerator. The number of vortices in both phases is independently determined from the simultaneously measured nuclear magnetic resonance (NMR) spectra . By tuning the barrier field or the temperature, the state of the sample can be changed from all A phase to all B phase or to a two-phase configuration. The evolution of the quantized vorticity as a function of $`\mathrm{\Omega }`$ is then observed to be radically different when the AB interface is present.

The quasi-isotropic 3He-B supports singly quantized vortices with a core size comparable to $`\xi `$ . In the anisotropic 3He-A we form vortex lines  with continuous skyrmion topology: Inside its central part the order-parameter amplitude remains constant but the axis of the orbital anisotropy covers a solid angle of $`4\pi `$. Such a “soft-core” structure carries continuous vorticity with two circulation quanta and is three orders of magnitude larger than the core of the B-phase vortex. Thus converting an A-phase vortex into a B-phase vortex requires large concentration of the flow energy.

The large difference in core radii is the origin for the much lower rotation at which vortices start forming in A phase at the outer sample circumference, compared to B phase . If the sample consists of only B phase then the critical velocity for forming the first vortex line in the setup of Fig. 1 is $`v_{\mathrm{cB}}>7`$mm/s , while the critical velocity $`v_{\mathrm{cA}}`$ is a factor of 20 smaller , independently of the presence of the AB boundary.

Vortex lines are thus easily created at low rotation in the A-phase section of the sample, while no vortices are detected in the B phase section. This is the ideal non-dissipative initial state where the two superfluid phases slide along each other without friction. When $`\mathrm{\Omega }`$ is increased further, sudden bursts of vortex lines are observed in the B-phase section, as shown in Fig. 2. The onset $`\mathrm{\Omega }_\mathrm{c}`$ depends on temperature and on the current in the barrier magnet. Overall, the measured characteristics of the bursts fit a shear flow instability, which provides a mechanism for the circulation to cross the AB interface.

3He-A and 3He-B are states of the same order-parameter manifold. One of the conditions on their phase boundary is that the phase of the order parameter has to be continuous : If the circulation is not continuing across the interface ( Fig. 3A) then the existing A-phase vortex lines have to bend and form a vortex layer on the AB interface. This is the only hydrodynamically stable state, with vortex-free flow seen by the B-phase spectrometer and a large cluster of vortex lines detected by the A-phase spectrometer. The build up of a sheet-like vortex layer means that the A-phase circulation cannot easily penetrate into the vortex-free B-phase and the AB interface remains stable up to the measured critical velocity $`u_c=\mathrm{\Omega }_\mathrm{c}R2`$$`4\mathrm{mm}/\mathrm{s}<v_{\mathrm{cB}}`$.

At $`\mathrm{\Omega }>\mathrm{\Omega }_\mathrm{c}`$, some vortex lines have broken through the AB phase boundary. Thus there exists also a stable configuration in which vortices from the A phase continue into the B phase, where they form a cluster in the center (Fig. 3B). The likely topology of the intersection is illustrated in Fig. 3C. Thus, although the AB interface is not directly monitored by the two NMR spectrometers, the hydrodynamically stable states below and above $`\mathrm{\Omega }_\mathrm{c}`$ can only be the configurations in Fig. 3. By comparing B-phase vortex-line creation in the presence and absence of the AB interface, we conclude that the events in Fig. 2 originate from the AB phase boundary.

In the classical KH instability the interface between fluids with densities $`\rho _1`$ and $`\rho _2`$ becomes destabilized by inertial effects, which are normally balanced by gravity $`g`$ and surface tension $`\sigma `$. The relative velocity $`|v_2v_1|`$ acts as a drive. When it reaches a critical value  given by

$$\frac{\rho _1\rho _2}{\rho _1+\rho _2}(v_2v_1)^2=2\sqrt{\sigma F},$$
(1)
where $`F=g(\rho _1\rho _2)`$ is the gravitational restoring force, waves with wave vector $`k=\sqrt{F/\sigma }`$ are created on the interface.

This approach can be generalized for superfluids in terms of two-fluid hydrodynamics . For superfluid 3He under our experimental conditions it is safe to assume that the normal fractions are always in solid-body rotation. Instead of gravity, the restoring force is now produced by the magnetic barrier $`H(z)`$, owing to the difference in the susceptibilities $`\chi _\mathrm{A}`$ and $`\chi _\mathrm{B}(T,H)`$: $`F=(1/2)(\chi _\mathrm{A}\chi _\mathrm{B}(T,H))(H^2)`$ at $`H_{\mathrm{AB}}(T)`$. The restoring force can thus be calculated as a function of temperature and current in the barrier solenoid. The motion of the AB interface with respect to the normal component is subject to a finite damping . Since the initial state is non-dissipative, damping does not explicitly appear in the stability condition, but modifies the build-up rate of the interface perturbations when the instability develops. The onset of the instability can be derived from the dynamics of small amplitude perturbations, as Eq. (1) is generally derived in textbooks , or from the thermodynamics when perturbations of the interface lead to a negative free energy in the rotating frame. The mode for which the interface first becomes unstable has the same wave vector $`k=\sqrt{F/\sigma }`$ as before at a drive given by

$$\frac{1}{2}\rho _{\mathrm{sA}}(v_{\mathrm{sA}}v_\mathrm{n})^2+\frac{1}{2}\rho _{\mathrm{sB}}(v_{\mathrm{sB}}v_\mathrm{n})^2=\sqrt{\sigma F},$$
(2)
where $`\rho _\mathrm{s}`$, $`\rho _\mathrm{n}`$, and $`v_\mathrm{s}`$, $`v_\mathrm{n}`$ are the densities and velocities of superfluid and normal components . In Fig. 4 we plot the prediction for the instability for three magnetic barrier profiles. With no fitting parameters the agreement with experiment is surprisingly good.

From Fig. 2 it is seen that of order $`\delta N10`$ vortex lines enter the B phase in one instability: in fact, at $`0.77T_\mathrm{c}`$ (with $`u_c=0.39`$ cm/s at 29.0 bar) the measured average for more than 100 events is $`\delta N11`$. To compare with Eq. (2), $`\delta N`$ corresponds to the number of circulation quanta $`\kappa =h/2m_3`$, which fit in one corrugation of the interface mode of size $`\lambda _\mathrm{c}/2=\pi /k`$. In practice, the counterflow $`|v_{\mathrm{sA}}v_\mathrm{n}|`$ in A-phase is limited by the small critical velocity $`v_{\mathrm{cA}}`$ , and we may set $`v_{\mathrm{sA}}v_\mathrm{n}0`$ in Eq. (2). In solid-body rotation there are then $`N\pi R^2\mathrm{\Omega }_\mathrm{c}/\kappa `$ vortex lines in the A phase which all flare out into the lateral sample boundary at the AB interface (Fig. 3A). Measured along the perimeter of the sample there are $`\mathrm{\Omega }_\mathrm{c}R/\kappa `$ circulation quanta per unit length and thus in one corrugation $`\delta N(\pi u_\mathrm{c})/(k\kappa )`$. From Eq. (2) this is seen to be $`\delta N(2\pi \sigma )/(\kappa u_\mathrm{c}\rho _{\mathrm{sB}})`$, which at $`0.77T_\mathrm{c}`$ gives 9 vortex lines in agreement with the measured number.

Eq. (2) has more interesting properties. As the normal components of both phases are in solid-body rotation and thus at rest in the rotating frame, it is the “superfluid winds” – the flow of the superfluid component, $`|v_\mathrm{s}v_\mathrm{n}|`$, on each side of the interface – which produce the instability. It takes place even if the two superfluids have the same densities and velocities. In this sense it resembles the flapping flag instability discussed by Rayleigh : The role of the flagpole, which fixes the reference frame, falls here on the normal component, which moves with the rotating container.

It is interesting to compare Eqs. (1) and (2) in the low temperature limit, when $`\rho _\mathrm{n}0`$ and $`\rho _{\mathrm{sA}}\rho _{\mathrm{sB}}\rho `$. We find that Eq. (2) gives a critical velocity which is by $`\sqrt{2}`$ smaller than that from Eq. (1). In fact, at $`T=0`$ the valid condition is Eq. (1), since $`\rho _\mathrm{n}`$ is exactly zero and no normal component is left to provide a reference frame for the now undamped interface . The ideal KH condition is then restored, as the only remaining drive is the velocity difference between the superfluid components. Not surprisingly, the two cases $`\rho _\mathrm{n}=0`$ and $`\rho _\mathrm{n}0`$ give different onsets for the instability: For a finite-size system the result would depend on the observation time – the time one waits for the interface to be coupled to the laboratory frame.

The shear-flow instability in superfluids thus provides an explicit bridge between the ideal inviscid and the viscous hydrodynamics. For two-fluid hydrodynamics it may have important consequences. It represents a new surface mechanism for vortex generation which might be applicable also to other interfaces, of which the normal state – superfluid boundary is of major interest. This interface occurs, for example, when the bulk superflow instability limit is reached in rotating 3He-B at the sharpest surface spike on the lateral sample boundary . A bubble of normal liquid is then created around this spike. Since the surrounding superflow moves at a high relative velocity, the bubble’s interface may undergo a shear-flow instability and promote the generation of a vortex line.

A related case has been investigated in bulk 3He-B which is in vortex-free rotation while it is irradiated with slow neutrons . The absorption process of one neutron creates locally a bubble of $`100\mu `$m size which is overheated up to the normal state. The bubble then quench cools back to the temperature of the surrounding superflow, as the normal–superfluid interface of the bubble contracts at high velocity. The measurements prove directly that vortices are formed in this rapid non-equilibrium transition and that the process is consistent with a volume effect known as the cosmological Kibble-Zurek (KZ) mechanism . It has been suggested as the origin for cosmic-string formation when rapid non-equilibrium phase transitions swept across the Early Universe. The 3He-B neutron experiment has been examined by solving analytically and numerically the time-dependent Ginzburg-Landau equation . This study shows that in the non-equilibrium transition of the neutron bubble there exists also a competing surface mechanism, when the rapidly moving normal–superfluid interface becomes unstable and vortex rings form around the bubble in a corrugation instability. This is direct theoretical evidence for the existence of the shear-flow instability also at the normal–superfluid interface, between the surrounding superflow and the stationary normal liquid within the bubble. A similar situation is the normal–superfluid boundary in a gaseous rotating Bose-Einstein condensate. Here the shear-flow instability also leads to a deformation of the surface of the condensate cloud by ripplons and to vortex-line formation .

The shear-flow instability thus appears to be a universal surface mechanism for defect formation in superfluids and as intrinsic as the KZ mechanism in the bulk liquid. We are left to wonder whether or not these two are the only possible processes in the clean limit at the lowest temperatures, to which other observations of vortex formation can be reduced.


format: text

text:
# Ultra-High Energy Cosmic Rays: The Annihilation of Super-Heavy Relics

## 1 Introduction

The detection of particles with energy in excess of $`10^{20}`$ eV is one of the unsolved mysteries of modern astrophysics and might represent one of the many signals of the existence of physics well beyond the standard model of particle physics. In  it was proposed the appealing possibility that this mystery might be related to another problem, the origin of dark matter. In fact the dark matter might be made of SH relics of the big bang  with lifetimes possibly exceeding the age of the universe. The rare decays of these SH particles (with masses in excess of $`10^{12}`$ GeV) naturally generates UHECRs from the top, that is as end-product of their decay. In order to saturate the dark matter content of the universe and at the same time explain the fluxes of UHECRs, the lifetime of the SH relics should be of order $`10^{22}`$ years, which is difficult to achieve for such massive particles unless some discrete conserved symmetry is introduced that prevents the natural decay of the relics on much shorter time scales.

In this paper we investigate an alternative possibility, consisting in the annihilation of SH relics rather than in their decay. The simple assumption that dark matter is a thermal relic limits the maximum mass of the dark matter particle to less than a few hundred TeV. However, it has been recently pointed out that dark particles might have never experienced local chemical equilibrium during the evolution of the universe, and that their mass may be in the range $`10^{12}`$ to $`10^{19}`$ GeV, much larger than the mass of thermal WIMPs .

In this paper we calculate the expected flux of UHECRs from the smooth dark matter content of our Galaxy, and the flux due to clumps of dark matter in the halo.

## 2 SH Dark matter annihilation and UHECRs

In calculating the UHE cosmic ray flux from a smooth SH dark matter distribution in the galactic halo, we assume a halo density spherically symmetric about the galactic center, in the form of a NFW profile  $`n_X(d)=\frac{N_0}{d(d+d_s)^2}`$. We will use $`d_s=3d_{}=24`$kpc in our numerical estimates, where $`d_{}8`$kpc is the distance of the solar system from the galactic center. The dimensionless parameter $`N_0`$ may be found by requiring that the total mass of the Galaxy is $`2\times 10^{12}M_{}`$. For simplicity, we will assume that each annihilation produces two quark initiated jets, each of energy $`M_X`$, while the decay of a SH particle produces two jets, each of energy $`M_X/2`$. The energy spectrum of observed UHE cosmic ray events from annihilation is determined by the energy distribution of the particles in each jet, which we take to equal the fragmentation function in the modified leading logarithmic approximation (MLLA) , which was also employed in .

We also use the MLLA limiting spectrum in the results of Fig. 1. Calculating the resulting UHE cosmic ray flux in the annihilation model, and comparing it to the similar calculation in the decay scenario, we obtain the results shown in Fig. 1 (the data points are from AGASA). In order to provide enough events to explain the observed UHE cosmic rays, $`\sigma _Av`$ has to be in the range $`10^{25}\text{cm}^2`$ to $`10^{27}\text{cm}^2`$, which is well in excess of the unitarity bound to the $`l`$-wave reaction cross section . The unitarity bound essentially states that the annihilation cross section must be smaller than $`M_X^2`$. However, as emphasized by Hui , there are several ways to evade the bound. The annihilation cross section may be larger if there are fundamental length scales in the problem larger than $`M_X^1`$.

A related issue is the typical energy of the annihilation products. In this paper we assume that annihilation produces two jets, each with energy approximately $`M_X`$. It is easy to imagine that with the finite-size effects discussed above, there is the possibility that annihilation will produce many soft particles, rather than essentially two particles each of energy $`M_X`$. An example that suits our needs is the annihilation of a monopole-antimonopole pair, as discussed at length in . We regard the requisite size of the annihilation cross section to be the most unattractive feature of the annihilation scenario.

This problem becomes less severe when the clumped component of the dark matter distribution in the Galactic halo is taken into account. The spatial and mass distribution of the clumps is taken from numerical N-body simulations (see  and ). The presence of cuspy density profiles in the clumps makes their contribution to the diffuse flux of UHECRs dominant over the NFW smooth component. The ratio of the clumped to smooth contributions for isothermal and NFW profiles for the clumps and for different values of $`M_X`$ is plotted in Fig. 2.

In the case of isothermal clumps, where the density profile of a single clump scales like $`1/r^2`$, the main contribution to the annihilation comes from the inner regions of the clumps, so that it is crucial to define a minimum radius $`R_{min}`$, determined by the annihilation efficiency. In  it is explained in detail the effect of two different definitions for $`R_{min}`$.

Another peculiar prediction of the annihilation of SH relics in clumps of dark matter in the halo consists of a peculiar pattern of small scale anisotropies. Several multiple events within degree scales would be observed if annihilation in the clumps generates the UHECRs. A detailed simulation was carried out in . The model also predicts a large scale anisotropy in the direction of the galactic center, which is an important test to be carried out by the future experiments, that will be able to see the direction of the galactic center.

## 3 Conclusions

Our calculations can be summarized in the following points: 1) the annihilation scenario for the origin of UHECRs requires cross sections that may be larger that the unitarity bound. Although not a killing factor, for the reasons explained above, this is certainly an unappealing feature of the model; 2) the annihilations in clumps of galactic dark matter dominate over the contribution of a smooth NFW dark matter profile, requiring correspondingly smaller values for the annihilation cross section; 3) the pattern of anisotropies is such that both a large scale (dipole) anisotropy and a degree scale anisotropy is to be expected if UHECRs are generated in the annihilations of SH relics in galactic SH dark matter.


text:
# Radiative cooling in SPH hydrodynamical simulations of X-ray clusters

## 1. Introduction

Hydrodynamical simulations have been widely used to predict for different theoretical frameworks the time evolution of the temperature distribution of X-ray clusters. With increasing availability of computational power numerical hydrodynamical simulations have been attempted to model the effects of radiative cooling on the gas in the formation and evolution of cluster galaxies (Katz & White 1993; Anninos & Norman 1996; Yoshikawa, Jing & Suto 2000; Pearce et al. 2000; Lewis et al. 2000). In this contribution I will show the preliminary results that have been obtained from a series of hydrodynamical SPH simulations of galaxy clusters. The simulations have different numerical resolution and include the effects on the gas component of radiative losses, star formation and energy feedback from SN. Final profiles are compared in order to assess the effects of numerical resolution, or different star formation prescriptions, on the cluster X-ray variables.

## 2. Simulations

In a previous paper (Valdarnini et al. 1999, hereafter VGB) a large set of hydrodynamical simulations was used to study the global X-ray cluster morphology and its evolution. The simulations were run using a TREESPH code with no gas cooling or heating. I refer to VGB for a detailed description of the simulations. In Valdarnini (2001) I have used the same cluster sample to study the effects of including in the simulations additional physics such as gas cooling and star formation. I will report here the results for the cluster with label $`\mathrm{\Lambda }`$CDM$`00`$ in VGB. This is the most massive cluster ($`M_v1.510^{15}M_{}`$) extracted from a cosmological $`\mathrm{\Lambda }`$CDM N-body simulation with size $`L=200h^1Mpc`$, matter density $`\mathrm{\Omega }_m=0.3`$ and Hubble constant $`H_0=70Kmsec^1Mpc^1`$. In order to check the effects of radiative cooling for this cluster a set of TREESPH simulations was performed, with initial conditions provided by the cosmological simulation. The effects of radiative cooling are modelled in these simulations by adding to the SPH thermal energy equation an energy-sink term represented by the radiative cooling function. The simulation parameters are: number of gas particles, $`N_g22,500`$; gas particle softening parameter $`\epsilon _g30Kpc`$; gas particle mass, $`m_g10^{10}M_{}`$; number of dark particles, $`N_d45,000`$ (including those in an external shell surrounding the cluster, cf. VGB).

Allowing the gas to cool radiatively will produce dense clumps of gas at low temperatures ($`10^4{}_{}{}^{}K`$). In these regions the gas will be thermally unstable and will likely meet the physical conditions to form stars. In TREESPH simulations star formation (SF) processes have been implemented using different approaches . According to Katz, Weinberg & Hernquist (1996, hereafter KWH) a gas particle is in a star forming region if the flow is convergent and the local sound crossing time is larger than the dynamical time (i.e. the fluid is Jeans unstable). In a simplified version Navarro & White (1993 , hereafter NW) assume as a sufficient condition that a gas particle must be in a convergent flow and its density exceeds a threshold ($`\rho _g>\rho _{g,c}=710^{26}grcm^3`$). If a gas particle meets these criteria then it is selected as an eligible particle to form stars. The local star formation rate (SFR) obeys the equation

$$d\rho _g/dt=c_{}\rho _g/\tau _g=d\rho _{}/dt,$$
(1)
where $`\rho _g`$ is the gas density, $`\rho _{}`$ is the star density , $`c_{}`$ is a characteristic dimensionless efficiency parameter, $`\tau _g`$ is the local collapse time and it is the maximum of the local cooling time $`\tau _c`$ and the dynamical time $`\tau _d`$. At each step the probability that a gas particle will form stars in a time step $`\mathrm{\Delta }t`$ is compared with a uniform random number. If the test is successful then a mass fraction $`\epsilon _{}`$ of the gas mass is converted into a new collisionless particle. This star particle has the position, velocity and gravitational softening of the original gas particle. Typical assumed values are $`\epsilon _{}=1/3`$ and $`c_{}=0.1`$ (KWH). For the local SFR, NW adopt Eq. 1 with $`c_{}=1`$, $`\tau _g=\tau _d`$ and $`\epsilon _{}=1/2`$ when a gas particle can convert part of its mass into a star particle. The numerical tests have been performed following the NW prescriptions for selecting gas particles which can form stars. Once a star particle is created it can release energy into the surrounding gas through supernova (SN) explosions. All the stars with mass above $`8M_{}`$ end as a SN, leaving a $`1.4M_{}`$ remnant. The SN energy ($`10^{51}erg`$) is released gradually into the gas according to the lifetime of stars of different masses.

## 3. Results

Three simulations were performed with the same numerical parameters, but using different SF methods or parameters. This in order to test the final dependence of the cluster X-ray variables on the assumed star formation model for the cooled gas. In simulation I cold gas particles have been converted into stars according to the NW method. In the other two runs (II and III) the KWH prescription is adopted for converting gas particles into stars, but with different values of the star formation efficiency parameter $`c_{}`$$`0.1`$ and $`1.`$. As is shown in Fig. 1 the most important differences are found between the KWH simulations with different $`c_{}`$. The differences are dramatic in the final X-ray luminosities, which differ by a factor $`40`$. The source of this discrepancy relies in the different gas density profiles, which have substantial differences in the cluster core regions for $`r\stackrel{<}{}100Kpc`$. These differences are localized at the cluster center, beyond $`r100Kpc`$ all the profiles converge, as it is shown in the plots of Fig. 1. The temperature profiles have a peak value of $`10^8{}_{}{}^{}K`$ at $`r100Kpc`$ and thereafter decline outward by a factor $`2`$ out to $`r_v`$. Below $`100Kpc`$ the profiles instead show large differences. Compared to the NW run the simulation with $`c_{}=0.1`$ has gas temperatures which decline inwards by a factor $`10`$ from $`100Kpc`$ down to $`r10Kpc`$. These radial decays follow because of the less efficient conversion of the cooled gas into stars compared to the NW run. There is a remarkable agreement for the ratio of cluster mass locked into stars to the gas mass, which is $`10\%`$ at $`r_v`$ for all the runs considered. In order to assess the effects of numerical resolution upon final results simulations I, II and III have been run again but with a number of particles increased by a factor $`3`$. These simulations will be referred as IH, IIH and IIIH, respectively. The simulation results are shown in Fig. 1. For simulations IH there are not appreciable differences in the radial profiles. The profiles of simulation IIH are instead different from those of run II at $`r\stackrel{<}{}50Kpc`$. The strong drop in $`T(r)`$ has been removed and the gas density profile is much closer to the NW one. Simulation IIIH yields final profiles very similar to the ones of the parent simulation. The bottom left panel of Fig. 1 shows that high resolution runs have final X-ray luminosities which can differ within a factor $`2`$ from the parent simulations.

To summarize, the above results demonstrate that simulations I and III have an adequate numerical resolution to reliably predict X-ray cluster properties, such as the X-ray luminosity. For simulation II ( KWH with $`c_{}=0.1`$ ) there are large differences at the cluster core between the final profiles when the numerical resolution is increased.

## References

Anninos, P. & Norman, M.L. 1996, ApJ, 459, 12

Katz, N. & White, S.D.M. 1993, ApJ, 412, 455

Katz, N., Weinberg, D.H. & Hernquist, L. 1996, ApJS, 105, 19

Lewis, G.F., Babul, A., Katz, N., Quinn, T., Hernquist, L. & Weinberg, D.H. 2000, ApJ 536, 623

Navarro, J. & White, S.D.M. 1993, MNRAS, 265, 271 (NW)

Navarro, J., Frenk C.S. & White, S.D.M. 1995, MNRAS, 275, 720

Pearce, F.R., Thomas, P.A., Couchman, H.M.P. & Edge, A.C. 2000, MNRAS 317, 1029

Valdarnini, R., Ghizzardi, S. & Bonometto, S. 1999, New Astr, 4, 71 (VGB)

Valdarnini, R. 2001, ApJ, in press

Yoshikawa, K., Jing, Y.P. & Suto, Y. 2000, ApJ, 535, 593 (YJS)


text:
# Fresnel laws at curved dielectric interfaces of microresonators

\[

## Abstract

We discuss curvature corrections to Fresnel’s laws for the reflection and transmission of light at a non-planar refractive-index boundary. The reflection coefficients are obtained from the resonances of a dielectric disk within a sequential-reflection model. The Goos-Hänchen effect for curved light fronts at a planar interface can be adapted to provide a qualitative and quantitative extension of the ray model which explains the observed deviations from Fresnel’s laws.

\]

The fabrication of lasing microresonators  and opto-mechanical microdevices  has generated a surge of interest in the confinement and propagation of light in small dielectric structures. Some understanding has been achieved from the ray optics of these systems, complemented by Fresnel’s laws of refraction and reflection at the interfaces, e.g., in order to identify and describe the relevant resonator modes . Fresnel’s laws give the probability of reflection and refraction of plane electromagnetic waves at planar interfaces of media with different refractive index $`n`$. Microresonators, however, often are so small that the curvature of their boundary cannot be neglected.

In this paper we investigate, in the limit of large wavenumbers, the corrections to the Fresnel coefficients which appear due to the curvature of the dielectric interface. The reflection coefficients are obtained via a sequential-reflection model  from the resonance widths of a microresonator, which are analytically accessible for large wavenumbers. The deviations from Fresnel’s laws are most noticeable around the critical angle for total internal reflection, $`\chi _c=\mathrm{arcsin}(1/n)`$ (where the refractive index of the surrounding medium is set to unity) and amount in a systematic reduction of the reflection probability.

The reduction of the reflection probability is conventionally related to tunneling escape at the curved interface. In view of the previous success of the ray model, which often is desired to be retained for its simplicity, we provide an alternative qualitative and quantitative explanation of the deviations by incorporating into this model the Goos-Hänchen effect . This effect results from the interference of rays in a beam with slight variations of the angle of incidence and consists of a shift of the effective plane of reflection. At a curved interface, it turns out that the reflection probability is then reduced because the angle of incidence at the effective interface is smaller than at the physical interface.

There is evidence obtained in the context of quantum-mechanical scattering problems  that incorporating the Goos-Hänchen effect is equivalent to a semiclassical approximation. Also the effect has been used in Ref.  to explain the decreased spacing of resonances observed in experiments with dielectric spheres, in terms of an effective optical size of the cavity which is larger than its physical size. Our work can be seen as complementary to this previous study, because we are concerned with the resonance widths, not only the resonance energies.

Curvature corrections to Fresnel coefficients have been addressed in the past in a number of works, by applying various techniques; see for example Refs. . The works closest in spirit to the present paper are those which employ the complex ray method, e.g., to describe light rays approaching a disk from outside .

Although we restrict our discussion to circular interfaces, the results for the reflection coefficients should be applicable to microresonators of general shape as soon as locally the curvature can be approximated as a constant.

According to Fresnel’s laws, a plane electromagnetic wave incident on a planar dielectric interface with angle of incidence $`\chi `$ is reflected with the polarization-dependent coefficients

$`R^{\mathrm{TM}}={\displaystyle \frac{\mathrm{sin}^2(\chi \eta )}{\mathrm{sin}^2(\chi +\eta )}},R^{\mathrm{TE}}={\displaystyle \frac{\mathrm{tan}^2(\chi \eta )}{\mathrm{tan}^2(\chi +\eta )}},`$ (1)
where TM (TE) signifies transverse polarization of the magnetic (electric) field at the interface and $`\eta =\mathrm{arcsin}(n\mathrm{sin}\chi )`$ is the direction of the refracted beam (according to Snell’s law).

Let us compare the Fresnel coefficients with the reflection coefficients at a curved interface with radius of curvature $`r_c`$. Their angular dependence can be conveniently obtained from the energies and widths of resonance states in a two-dimensional circular disk of radius $`r_c`$. In this geometry the two possible polarization directions decouple and angular momentum (quantum number $`m`$) is conserved. We introduce polar coordinates $`r`$ and $`\varphi `$ and denote the (complex) wavenumber by $`k`$. We will concentrate on the case close to geometric optics $`\mathrm{Re}kr_c1`$.

The resonance states are obtained by matching the wave field $`J_m(nkr)e^{im\varphi }`$ inside the disk (with the Bessel function $`J`$) at $`r=r_c`$ to the wave field $`H_m^{(1)}(kr)e^{im\varphi }`$ outside the disk (with the Hankel function $`H^{(1)}`$), where the matching conditions follow from Maxwell’s equations:

$$J_m(nkr_c)H_{m}^{(1)}{}_{}{}^{}(kr_c)=nJ_m^{}(nkr_c)H_m^{(1)}(kr_c),$$
(2)
for TM polarization, and

$$nJ_m(nkr_c)H_{m}^{(1)}{}_{}{}^{}(kr_c)=J_m^{}(nkr_c)H_m^{(1)}(kr_c),$$
(3)
in the TE case (primes denote derivatives). Given a complex solution $`k`$, the angle of incidence is obtained from the real part by comparing the angular momentum in the ray picture and in the wave picture,

$$\mathrm{sin}\chi =\frac{m}{n\mathrm{Re}kr_c},$$
(4)
while the reflection probability

$$R=\mathrm{exp}(4n\mathrm{Im}kr_c\mathrm{cos}\chi )$$
(5)
follows from the imaginary part of $`k`$ because it determines the escape rate, which in turn can be related to $`R`$ by the sequential-reflection model .

Because the discrete set of resonance energies obtained from Eqs. (2) and (3) only is meaningful for the disk, let us first derive analytical expressions for the resonance width as a function of a continuous resonance energy which smoothly interpolate between these solutions. It is interesting to note  that one cannot simply expand Eqs. (2), (3) in $`\mathrm{Im}k`$ when $`k`$ is not close to an exact solution. Moreover, for TE polarization $`\mathrm{Im}k`$ would diverge at the Brewster angle if it is calculated by inserting the Fresnel coefficient (1) into Eq. (5). In order to achieve a more accurate expansion we separate out the problematic term and cast Eqs. (2) and (3) into the form

$$\frac{J_m^{}(nkr_c)}{J_m(nkr_c)}=(kr_c),$$
(6)
with

$`^{\mathrm{TE}}(x)=n{\displaystyle \frac{H_{m}^{(1)}{}_{}{}^{}(x)}{H_m^{(1)}(x)}},^{\mathrm{TM}}(x)`$ $`=`$ $`n^2^{\mathrm{TE}}(x),`$ (7)
depending on the polarization. In both cases, $`(kr_c)`$ is a slowly varying complex function of its argument, and the argument can be taken real because $`\mathrm{Re}k|\mathrm{Im}k|`$. The logarithmic derivative of Bessel functions, however, is a rapidly fluctuating function, and its dependence on $`\mathrm{Im}k`$ has to be worked out carefully. This can be achieved by approximation by tangents ,

$`{\displaystyle \frac{J_m^{}(nkr_c)}{J_m(nkr_c)}}=\mathrm{tan}\alpha \mathrm{cos}\chi ,`$ (8)
$`\alpha =m\mathrm{cot}\chi +m\chi m{\displaystyle \frac{\pi }{2}}{\displaystyle \frac{\pi }{4}}+in\mathrm{Im}kr_c\mathrm{cos}\chi ,`$ (9)
where $`\chi `$ is given as a function of $`\mathrm{Re}k`$ by Eq. (4). We expanded $`\alpha `$ linearly in $`\mathrm{Im}k`$ and neglected terms of order $`(\mathrm{Re}kr_c)^1`$. Equation (6) now can be solved for $`\alpha `$, without any further approximation. ¿From the imaginary part one deduces

$$\mathrm{Im}kr_c=\frac{1}{n\mathrm{cos}\chi }\mathrm{Im}\mathrm{arctan}\frac{}{\mathrm{cos}\chi }.$$
(10)
Although this can already be taken as the final result, we further may insert the uniform approximation

$$^{\mathrm{TE}}=in\mathrm{cos}\eta \left[1+\frac{1}{\mathrm{sin}^2\eta }\left(\frac{K_{2/3}(z)}{K_{1/3}(z)}1\right)\right],$$
(11)
and similarly for $`^{\mathrm{TM}}=^{\mathrm{TE}}/n^2`$, with the modified Bessel function $`K`$, the angle of refraction $`\eta `$ (which is a complex number for $`\chi >\chi _c`$), and $`z=i\mathrm{Re}kr_c\mathrm{cos}^3\eta /(3\mathrm{sin}^2\eta )`$. In Fig. 1 we illustrate that Eqs. (10), (11) agree very well with the exact solutions of Eqs. (2) and (3), even close to the Brewster angle for TE polarization, and interpolates smoothly in between.

The angular dependence of the reflection coefficients can now be obtained by combining Eqs. (5) and (10), giving the final result

$$R=\left|\frac{\mathrm{cos}\chi +i}{\mathrm{cos}\chi i}\right|^2.$$
(12)
Figures 2 and 3 show the result for the two values $`\mathrm{Re}kr_c=50`$ and $`\mathrm{Re}kr_c=150`$, respectively. Deviations from Fresnel’s laws (1) are most visible around the critical angle $`\chi _c=\mathrm{arcsin}(1/n)`$ where the reflection coefficients increase rapidly as the regime of total internal reflection is approached. The correction consists not only in a broadening of the transition interval, but most notably also in a shift of this transition region towards higher angles of incidence, resulting in a systematic reduction of the reflection coefficient.

The deviations from Fresnel’s laws in Figs. 2 and 3 increase as $`\mathrm{Re}kr_c`$ is reduced, that is, the more noticeable the curvature of the interface is. On the other hand, in the zero-wavelength limit $`\mathrm{Re}kr_c\mathrm{}`$ of geometric optics any interface appears planar, and Fresnel’s laws should apply without modification. Indeed it can be seen that they are recovered from Eq. (12) when the approximation by tangents is also applied to the Hankel functions in $``$, resulting in $`^{\mathrm{TM}}=in^1\mathrm{cos}\eta `$, $`^{\mathrm{TE}}=in\mathrm{cos}\eta `$. The deviations close to the critical angle are directly related to the break-down of this approximation when the argument of the Hankel functions becomes smaller than the index. As we will discuss now, the curvature corrections to Fresnel’s coefficients can be obtained within a minimal extension of the ray picture when the Goos-Hänchen effect is taken into account (the result obtained is given by the dashed curves in Figs. 2 and 3).

The Goos-Hänchen effect  refers to the displacement of the reflected beam when the incident beam consists of rays with slight variations of the angle of incidence, and arises because each ray experiences a slightly different phase shift when it is reflected. As is illustrated in Fig. 4a, the lateral shift along the interface also can be interpreted as resulting from a displacement of the effective plane of reflection (a signature of this displacement is the increased Wigner delay time which recently has been measured at metallic gratings .)

Since at a planar interface the reflection law is not affected by parallel displacement and fulfilled for the mean angles of incidence, the reflection coefficients are not affected by the Goos-Hänchen effect—the only consequence of the slight variation of angles is that their angular dependence is smeared out. However, the situation changes at a curved interface  as shown in Fig. 4b. The intersection of the incident and the laterally shifted reflected ray defines an effective boundary of radius $`r_c^{}>r_c`$. We now can assume that the ray is specularly reflected at the effective boundary, resulting in a smaller effective angle of incidence $`\chi ^{}<\chi `$, and evaluate the Fresnel reflection coefficients (1) at this smaller angle $`\chi ^{}`$. Since the reflection probability is then reduced this qualitatively explains the observed deviations from Fresnel’s laws in Figs. 2 and 3.

For a quantitative comparison we need the distribution of angles of incidence $`P(\chi )`$, which can be related to the radial width $`r_c(n\mathrm{Re}kr_c)^{2/3}`$ of the caustic in the radial wave field $`J_m(nkr)`$. (The same distribution of angles will also be used to smooth out the reflection coefficients as is appropriate even for the planar interface.) The variation of angles of incidence arises from the intrinsic curvature of the beam wave front and also directly from the curvature of the interface—note that both mechanisms should contribute equally to the Goos-Hänchen effect. We found that for our purposes $`P(\chi )`$ is sufficiently well approximated by that of a Gaussian beam of half width $`\sigma =r_c(n\mathrm{Re}kr_c)^{2/3}r_c`$. It is good to observe that this Gaussian beam geometry does not put us into conflict with the finite disk size while at the same time leaving us in the paraxial regime $`n\mathrm{Re}k\sigma =(n\mathrm{Re}kr_c)^{1/3}1`$. For TM polarization the effective radius of curvature $`r_c^{}`$ can then be calculated by applying the formulas of Ref.  (which are lengthy expressions and hence not given here). The result of this Goos-Hänchen-effect-based approach is presented as the dashed curves in the upper panels of Figs. 2 and 3. There is good agreement with the exact results obtained within the sequential-reflection model. For TE polarization and the chosen refractive index, the analytic results in the literature becomes applicable only for $`\mathrm{Re}kr_c1000`$. Surprisingly (see however Ref. ), in the current situation nice agreement is found by simply assuming that the Goos-Hänchen shift is identical to the TM case, as is indicated by the dashed curves in the lower panels of Figs. 2 and 3. (For $`\mathrm{Re}kr_c1000`$, however, it is appropriate to work with the correct TE formulas.)

In conclusion, we investigated the reflection coefficients at a curved refractive index boundary by relating them to resonances in a circular dielectric disk and derived analytic expressions valid in the limit of large wavenumbers. The deviations from Fresnel’s laws can be explained within geometric optics by incorporating the Goos-Hänchen effect. In this work we concentrated on the wave field confined by reflection inside the microcavity. The Goos-Hänchen shift also affects the wave field outside the cavity , because the modified angle of incidence results in a change of the angle of refraction. It is desirable to investigate the implications on the remarkable emission directionality of non-circular devices , which sometimes departs substantially from what is expected from geometric optics .

We gratefully acknowledge helpful discussions with M. Berry, S. W. Kim, J. U. Nöckel, K. Richter, R. Schubert, and J. Wiersig.


text:
# Spectrum of a duality-twisted Ising quantum chain

## References


text:
# Chiral Properties of Quenched and Full QCD

## Abstract

We extend a technique for the chiral extrapolation of hadron masses calculated with dynamical fermions to those generated by quenched simulations. The method ensures the correct leading and next-to-leading non-analytic behaviour for either QCD or quenched QCD in the chiral limit, as well as the correct large quark mass behaviour. We find that the primary difference between quenched and dynamical baryon masses can be described by the meson loops which give rise to the different leading and next-to-leading non-analytic behaviour.

preprint:  ADP-01-49/T481, LTH-528

Modern computing facilities, combined with innovations in improved actions for lattice QCD, mean that it is now possible to perform accurate quenched QCD (QQCD) simulations at quite low quark masses Bowler:1999ae ; Kanaya:1998sd ; Bernard:2001av ; Zanotti:2001yb . For simulations with dynamical fermions (full QCD) the situation is much more difficult, but there are initial results at quark masses as low as 30 MeV Bernard:2001av ; Aoki:1999ff . The latter development has inspired studies of chiral extrapolation aimed at using the full QCD data over a range of masses to reliably extract the physical hadron mass.

In general, effective field theories, such as chiral perturbation theory, lead to divergent or asymptotic expansions Dyson:1952tj ; LeGuillou:1990nq . While this raises doubts about the direct application of chiral perturbation theory to lattice data, studies of the mass dependence of hadron properties in QCD-inspired models magmom ; Hackett-Jones:2000js ; Detmold:2001jb , as well as the exactly soluble Euler-Heisenberg problem Dunne:2001ip , suggest that one can develop surprisingly accurate extrapolation formulas, provided one builds in the correct behaviour in both the small and large mass limits. For the nucleon ($`N`$) and delta ($`\mathrm{\Delta }`$) masses (and by extension all other baryons), Leinweber et al. Leinweber:1999ig  have suggested an extrapolation method which ensures both the exact low mass limit of chiral perturbation theory (technically its leading (LNA) and next-to-leading non-analytic (NLNA) behaviour) and the heavy quark limit of heavy quark effective theory (HQET). The transition between the chiral and heavy quark regimes is characterised by a mass scale $`\mathrm{\Lambda }`$, related to the inverse of the size of the pion cloud source. The rapid, non-analytic variation of hadron properties, characteristic of chiral perturbation theory, is rapidly suppressed once the pion Compton wavelength is smaller than this size (i.e. $`m_\pi >\mathrm{\Lambda }`$).

It is straightforward to extend the method of Ref. Leinweber:1999ig  to QQCD. One simply includes all the Goldstone loops (including both $`\pi `$ and $`\eta ^{}`$) which give rise to the LNA and NLNA behaviour of quenched chiral perturbation theory (Q$`\chi `$PT) Bernard:1992mk ; Labrenz:1996jy . In principle, the parameters of the chiral Lagrangian are dependent on the number of dynamical fermion flavours, $`N_f`$, but are independent of the masses of the quarks. This is a celebrated feature of partially quenched chiral perturbation theory Sharpe:2000bc ; Sharpe:2001fh . The extent of the $`N_f`$ dependence is not well known, and a precise determination awaits dynamical fermion simulations with light dynamical quarks of varying number.

Phenomenological investigations Hackett-Jones:2000js ; Leinweber:1993hj  of the role of the pion cloud in hadronic charge radii indicate that results consistent with experiment can be obtained by adding full-QCD chiral corrections to the results of quenched simulationsLeinweber:1991dv  at moderate to heavy quark masses. This suggests that the size of the pion source is not changed dramatically in going from the quenched approximation to full QCD, motivating the use of a common scale, $`\mathrm{\Lambda }`$, in QQCD and QCD. We proceed to fit quenched and dynamical lattice data assuming negligible $`N_f`$ dependence in the chiral parameters and $`\mathrm{\Lambda }`$. The extent to which these assumptions hold can only be determined via further dynamical fermion calculations in the light quark regime.

By incorporating the chiral loops which give rise to the LNA and NLNA behaviour in QCD and QQCD we find a remarkable agreement between the fit parameters of each simulation. This supports our hypothesis that the behaviour of the source of the pion cloud within baryonic systems behaves much the same in both QQCD and QCD. The differences between quenched and full QCD is primarily described by the differences in the associated chiral loops. Since the chiral corrections are expected to be larger for $`N`$ and $`\mathrm{\Delta }`$ than for others, this suggests that a similar technique may be applicable to all baryons. This would enable quenched simulations to play a more valuable role, together with new experimental information from JLab and elsewhere, in deepening our understanding of baryon spectroscopy.

With regard to the properties of the $`N`$ and $`\mathrm{\Delta }`$ we find a spectacular difference in QQCD. Whereas the extrapolation of the $`N`$ mass is essentially linear in the quark mass, the $`\mathrm{\Delta }`$ exhibits some upward curvature in the quenched chiral limit. As a result, the $`\mathrm{\Delta }`$ mass in QQCD is expected to be of the order 300–400 MeV above its mass in full QCD. The success of the extrapolation scheme also lends confidence to the interpretation of the $`\mathrm{\Delta }N`$ mass splitting as receiving a contribution of order 50 MeV from pion loops in full QCD and up to 250 MeV in QQCD tocome . The residual splitting in full QCD would then be naturally ascribed to some shorter range mechanism, such as the traditional one-gluon-exchange Isgur:1999jv .

The method for extrapolating baryon masses proposed by Leinweber et al. Leinweber:1999ig  is to fit the lattice data with the form:

$$M_B=\alpha _B+\beta _Bm_\pi ^2+\mathrm{\Sigma }_B(m_\pi ,\mathrm{\Lambda }),$$
(1)
where $`\mathrm{\Sigma }_B`$ is the sum of those pion loop induced self-energies which give rise to the LNA and NLNA behaviour of the mass, $`M_B`$. In the case of the $`N`$ this is the sum of the processes $`NN\pi N`$ and $`N\mathrm{\Delta }\pi N`$ , while for the $`\mathrm{\Delta }`$ it involves $`\mathrm{\Delta }\mathrm{\Delta }\pi \mathrm{\Delta }`$ and $`\mathrm{\Delta }N\pi \mathrm{\Delta }`$. In the heavy baryon limit, these four contributions ($`BB^{}\pi B`$) can be summarised as:

$$\sigma _{BB^{}}^\pi =\frac{3}{16\pi ^2f_\pi ^2}G_{BB^{}}_0^{\mathrm{}}𝑑k\frac{k^4u_{BB^{}}^2(k)}{\omega (k)(\omega _{BB^{}}+\omega (k))},$$
(2)
where $`\omega (k)=\sqrt{k^2+m_\pi ^2}`$ and $`\omega _{BB^{}}=(M_B^{}M_B)`$, and the constants $`G_{BB^{}}`$ are standard SU(6) couplings Leinweber:1999ig .

The factor $`u(k)`$, which acts as an ultraviolet regulator, may be interpreted physically as the Fourier transform of the source of the pion field. Whatever choice is made, the form of these meson loop contributions guarantees the exact LNA and NLNA structure of chiral perturbation theory ($`\chi `$PT). Furthermore, such a form factor causes the self-energies to decrease as $`1/m_\pi ^2`$ for $`m_\pi >>\mathrm{\Lambda }`$. One commonly uses a dipole, $`u(k)=(\mathrm{\Lambda }^2\mu ^2)^2/(\mathrm{\Lambda }^2+k^2)^2`$, with $`\mu `$ the physical pion mass.

Quenched $`\chi `$PT is a low energy effective theory for quenched QCD Bernard:1992mk ; Labrenz:1996jy , analogous to $`\chi `$PT for full QCD Bernard:1995dp . Sea quark loops are formally removed from QCD by including a set of degenerate, bosonic quarks. These bosonic fields have the effect of cancelling the fermion determinant in the functional integration over the quark fields. This gives a Lagrangian field theory which is equivalent to the quenched approximation simulated on the lattice. The low energy effective theory is then constructed using the symmetry groups of this Lagrangian.

A study of the chiral structure of baryon masses within the quenched approximation has been carried out by Labrenz and Sharpe Labrenz:1996jy . The essential differences from full QCD are: a) in the quenched theory the chiral coefficients differ from their standard values and b) new non-analytic structure is also introduced. The leading order form of the baryon mass expansion about $`m_\pi =0`$ is

$`M_B`$ $`=`$ $`M_B^{(0)}+c_1^Bm_\pi +c_2^Bm_\pi ^2+c_3^Bm_\pi ^3`$ (3)
$`+c_4^Bm_\pi ^4+c_{4L}^Bm_\pi ^4\mathrm{log}m_\pi +\mathrm{}`$
where the coefficients of terms non-analytic in the quark mass are model-independent. Throughout we use couplings as given in Ref. Labrenz:1996jy . In addition, we have included an octet–decuplet mass splitting to explicitly give a value for $`c_{4L}^B`$ tocome . We also stress that the term in $`m_\pi `$ is absent in full QCD — such a term being unique to the quenched case.

In fitting quenched data we wish to replicate the analysis for full QCD while incorporating the known chiral structure of the quenched theory. The meson-loop, self-energy corrections to baryon masses can be described in the same form as for full QCD. The effect of quenching appears as a redefinition of the couplings in the loop diagrams in order that they yield exactly the same LNA and NLNA structure as given by Q$`\chi `$PT. For example, the analytic expressions for the pion-cloud corrections to the masses of the $`N`$ and $`\mathrm{\Delta }`$ have the same form as the full QCD integrals (Eq. 2) with redefined quenched couplings. We refer to Ref. tocome  for details. Assuming a weak $`N_f`$ dependence of the chiral parameters, we describe the quenched self-energies using the same tree level values of $`D=0.76`$ and $`F=0.50`$ as in full QCD.

In addition to the usual pion loop contributions, QQCD contains loop diagrams involving the flavour singlet $`\eta ^{}`$ which also give rise to important non-analytic structure. Within full QCD such loops do not play a role in the chiral expansion because the $`\eta ^{}`$ remains massive in the chiral limit. On the other hand, in the quenched approximation the $`\eta ^{}`$ is also a Goldstone boson Bernard:1992mk ; Sharpe:1992ft  and the $`\eta ^{}`$ propagator is exactly the same as that of the pion.

As a consequence there are two new chiral loop contributions unique to the quenched theory. The first of these, $`\stackrel{~}{\sigma }_B^{\eta ^{}(1)}`$, corresponds to a single hairpin diagram such as that indicated in Fig. 1(a).

This diagram is the source of the term proportional to $`m_\pi ^3`$ (involving the couplings $`\gamma `$ and $`\gamma ^{}`$ Labrenz:1996jy ) in the chiral expansion Eq. (3). The structure of this diagram is exactly the same as the pion loop contribution where the internal baryon is degenerate with the external state. The second of these new $`\eta ^{}`$ loop diagrams, $`\stackrel{~}{\sigma }_B^{\eta ^{}(2)}`$, arises from the double hairpin vertex as pictured in Fig. 1(b). This contribution is particularly interesting because it involves two Goldstone boson propagators and is therefore the source of the non-analytic term linear in $`m_\pi `$.

The total meson loop contribution to the baryon self-energy within the quenched approximation is given by the sum of these four diagrams:

$$\stackrel{~}{\mathrm{\Sigma }}_B=\stackrel{~}{\sigma }_{BB}^\pi +\stackrel{~}{\sigma }_{BB^{}}^\pi +\stackrel{~}{\sigma }_B^{\eta ^{}(1)}+\stackrel{~}{\sigma }_B^{\eta ^{}(2)}.$$
(4)
As the resultant pion couplings in QQCD are quite a bit smaller than the corresponding full QCD couplings, $`\stackrel{~}{\mathrm{\Sigma }}_B`$ is smaller in magnitude than $`\mathrm{\Sigma }_B`$. We display the individual contributions to the $`\mathrm{\Delta }`$ mass for both quenched and full QCD in Fig. 2. It is notable that $`\stackrel{~}{\sigma }_{\mathrm{\Delta }N}^\pi `$ and $`\stackrel{~}{\sigma }_\mathrm{\Delta }^{\eta ^{}(2)}`$ are repulsive, with $`\stackrel{~}{\sigma }_\mathrm{\Delta }^{\eta ^{}(1)}`$ vanishing, so that at light quark masses the total quenched chiral loop contribution to the $`\mathrm{\Delta }`$ mass is repulsive, whereas it is attractive in full QCD.

It is now straightforward to fit the quenched lattice data with the form:

$$\stackrel{~}{M}_B=\stackrel{~}{\alpha }_B+\stackrel{~}{\beta }_Bm_\pi ^2+\stackrel{~}{\mathrm{\Sigma }}_B(m_\pi ,\mathrm{\Lambda }).$$
(5)
Once again the linear part describes how the mass of the pion-cloud source varies with quark mass. This form includes the expected behaviour of HQET where the $`\pi `$ and $`\eta ^{}`$ loop contributions are suppressed. Since, as discussed earlier, the meson–baryon vertices are characterised by the source distribution, which is argued to be similar in quenched and full QCD, we take all vertices to have the same momentum dependence — i.e. a common dipole mass, $`\mathrm{\Lambda }`$. With this parameter fixed there are just two free parameters, $`\stackrel{~}{\alpha }`$ and $`\stackrel{~}{\beta }`$, to fit the quenched data for each baryon.

As described in Ref. Leinweber:2001ac  we replace the continuum integral over the intermediate pion momentum by a discrete sum over the pion momenta available on the lattice, thus encapsulating finite lattice spacing and volume artifacts.

The lattice data used for this analysis comes from a recent paper by Bernard et al. Bernard:2001av . These results are obtained using an improved Kogut-Susskind quark action, which is known to give good scaling properties Bernard:1999xx . Unlike the standard Wilson fermion action, masses determined at finite lattice spacing are excellent estimates of the continuum limit results. The physical scale of both full and quenched data sets has been set via a variant of the Sommer scale Bernard:2001av . This procedure, based on the static-quark potential where chiral corrections are negligible tocome , provides a self-consistent determination of the scale for both simulations.

Fits of the the form for full QCD, Eq. 1, and quenched QCD, Eq. 5, are shown in Fig. 3. In fitting to data we choose a dipole mass of $`\mathrm{\Lambda }=0.8\mathrm{GeV}`$, which has been optimised to highlight the remarkably similar behaviour of the pion-cloud source in both quenched and dynamical simulations. Phenomenologically, this agrees with quite general expectations that it should be somewhat smaller than that for the axial form factor Thomas:2001kw ; Guichon:1983zk ; Thomas:1989tv . The parameters obtained from our fits are shown in Table 1.

We highlight the fact that the best fit parameters, $`\alpha `$ and $`\beta `$, obtained from both the quenched and full simulations, agree within errors. This suggests that the quark mass behaviour of the pion-cloud source is quite similar in full and quenched QCD, lending support to the hypothesis made in this analysis. This leads one to conclude that the dominant effects of quenching can be attributed to the first order meson loop corrections.

We have investigated the quark mass dependence of the $`N`$ and $`\mathrm{\Delta }`$ masses within the quenched approximation. The leading chiral behaviour of hadron masses in quenched QCD is known to differ from the full theory. This knowledge has been used to guide the construction of a functional form which both reproduces this correct chiral structure, is consistent with current lattice simulations and encompasses the HQET properties. The success of this method in the quenched case further verifies the importance of including meson-induced self-energies when extrapolating lattice results.

We find that, although the quenched approximation gives rise to more singular behaviour in the chiral limit, these contributions are quickly suppressed with increasing quark mass. In the nucleon, the effects of quenching reduce the amount of curvature expected as lighter quark masses are simulated. In contrast, for the $`\mathrm{\Delta }`$ we find some upward curvature of the mass in QQCD as the quark mass approaches zero. In addition, the $`\mathrm{\Delta }`$$`N`$ mass splitting increases to around 400 MeV at the physical point. As a consequence of this behaviour, the $`\mathrm{\Delta }`$ mass in the quenched approximation is expected to differ from the physical mass by approximately $`25\%`$.

Our calculations suggest that the one-loop meson graphs which generate the leading and next-to-leading non-analytic behaviour are the primary difference between baryon masses in quenched and full QCD. Thus, rather than quenched lattice QCD being regarded as an uncontrolled approximation, we are able to make a quantitative estimate of errors over a range of quark mass. It is vital to investigate the assumptions made with further dynamical fermion simulations at low quark masses to test the extent to which the presented results hold. Nevertheless, this discovery represents a remarkable step forward in relating lattice QCD to observed hadronic properties.

We would like to thank S. Sharpe for numerous enlightening discussions concerning Q$`\chi `$PT as well as W. Detmold, M. Oettel, A. Williams and J. Zanotti for helpful conversations. This work was supported by the Australian Research Council and the University of Adelaide.


text:
# A characteristic frequency of two mutually interacting gas bubbles in an acoustic field

## Abstract

Transition frequencies of two spherical gas bubbles interacting in an acoustic field are discussed theoretically. In the present study, transition frequency is defined as the frequency of external sound for which the phase difference between a bubble’s pulsation and the external sound is $`\pi /2`$. It is shown by a linear theory that a bubble interacting with a neighboring bubble has three (or fewer) transition frequencies but only two natural frequencies. This result means that the bubble has a characteristic frequency besides the natural frequencies.

It is known that radiative interaction between gas bubbles varies their natural (or resonance) frequencies, around which the bubbles indicate resonance response. This variation in the natural frequencies has been the subject of many studies . In a case where two spherical gas bubbles of different radii interact acoustically, a bubble has two (or one) natural frequencies (see, e.g., ).

The aim of this Letter is to show that a gas bubble interacting with a neighboring bubble has an alternative characteristic frequency; we now call this a transition frequency. In the present work, a transition frequency is defined as the frequency of an external sound (the driving frequency) for which the phase difference between a bubble’s pulsation and the external sound becomes $`\pi /2`$. It is well known that, in a single-bubble case, the phase difference between the bubble and an external sound becomes $`\pi /2`$ when the driving frequency is equal to the bubble’s natural frequency . However, it is shown in this Letter that, in a double-bubble case, the number of transition frequencies differs from the number of natural frequencies. This result means that the bubbles in the latter case have a characteristic frequency that differs from the natural frequency.

First, we derive a mathematical model to use for the present investigation. This model describes the pulsation of two coupled bubbles. Since we are interested only in the qualitative natures of the transition frequencies, we use a classical theoretical model, called a coupled-oscillator or self-consistent model , although models with greater accuracy have been proposed .

A gas bubble immersed in a liquid pulsates when a sound wave is applied. The sound pressure at the bubble position drives the pulsation. When other bubbles (named “bubble 2” $``$ “bubble $`N`$”, where $`N`$ is the total number of bubbles) exist near the bubble (bubble 1), the sound waves scattered by the other bubbles also drive the pulsation of bubble 1. Namely, the driving pressure acting on bubble 1, $`p_{\mathrm{d}\mathrm{\hspace{0.17em}1}}`$, is expressed as

$$p_{\mathrm{d}\mathrm{\hspace{0.17em}1}}=p_{\mathrm{ex}}+\underset{j=2}{\overset{N}{}}p_{\mathrm{s}\mathrm{\hspace{0.17em}1}j},$$
(1)
where $`p_{\mathrm{ex}}`$ and $`p_{\mathrm{s}\mathrm{\hspace{0.17em}1}j}`$ are the sound pressures of the external sound field and the scattered wave emitted by bubble $`j`$, respectively, at the position of bubble 1. When the wavelength of the external sound wave is much larger than the bubbles’ radii and the distances between the bubbles, $`p_{\mathrm{ex}}`$ can be considered spatially uniform. By assuming that the sphericity of the bubbles is thoroughly maintained and that the liquid surrounding the bubbles is incompressible, the scattered pressure can be estimated with (see, e.g., )

$$p_{\mathrm{s}\mathrm{\hspace{0.17em}1}j}\frac{\rho }{r_{1j}}\frac{d}{dt}(R_j^2\dot{R}_j),$$
(2)
where $`\rho `$ is the density of the liquid, $`r_{1j}`$ is the distance between the centers of bubble 1 and bubble $`j`$, $`R_j`$ is the radius of bubble $`j`$, and the over dot denotes the time derivative.

The radial oscillation of a spherical bubble in an incompressible viscous liquid is described by the RPNNP equation :

$$R_1\ddot{R}_1+\frac{3}{2}\dot{R}_1^2\frac{1}{\rho }p_{w\mathrm{\hspace{0.17em}1}}=\frac{1}{\rho }p_{\mathrm{d}\mathrm{\hspace{0.17em}1}},$$
(3)
$$p_{w\mathrm{\hspace{0.17em}1}}=\left(P_0+\frac{2\sigma }{R_{10}}\right)\left(\frac{R_{10}}{R_1}\right)^{3\kappa }\frac{2\sigma }{R_1}\frac{4\mu \dot{R}_1}{R_1}P_0,$$

where $`P_0`$ is the static pressure, $`\sigma `$ is the surface-tension coefficient at the bubble surface, $`\mu `$ is the viscosity of the liquid, $`\kappa `$ is the polytropic exponent of the gas inside the bubbles, and $`R_{10}`$ is the equilibrium radius of bubble 1. Substituting Eqs. (1) and (2) into Eq. (3) yields

$$R_1\ddot{R}_1+\frac{3}{2}\dot{R}_1^2\frac{1}{\rho }p_{w\mathrm{\hspace{0.17em}1}}=\frac{1}{\rho }[p_{\mathrm{ex}}+\underset{j=2}{\overset{N}{}}\frac{\rho }{r_{1j}}\frac{d}{dt}(R_j^2\dot{R}_j)].$$

When $`N=2`$, this equation is reduced to

$$R_1\ddot{R}_1+\frac{3}{2}\dot{R}_1^2\frac{1}{\rho }p_{w\mathrm{\hspace{0.17em}1}}=\frac{1}{\rho }[p_{\mathrm{ex}}+\frac{\rho }{D}\frac{d}{dt}(R_2^2\dot{R}_2)],$$
(4)
where $`D=r_{12}`$ ($`=r_{21}`$). Exchanging subscripts 1 and 2 yields the model equation for bubble 2:

$$R_2\ddot{R}_2+\frac{3}{2}\dot{R}_2^2\frac{1}{\rho }p_{w\mathrm{\hspace{0.17em}2}}=\frac{1}{\rho }[p_{\mathrm{ex}}+\frac{\rho }{D}\frac{d}{dt}(R_1^2\dot{R}_1)],$$
(5)
$$p_{w\mathrm{\hspace{0.17em}2}}=\left(P_0+\frac{2\sigma }{R_{20}}\right)\left(\frac{R_{20}}{R_2}\right)^{3\kappa }\frac{2\sigma }{R_2}\frac{4\mu \dot{R}_2}{R_2}P_0.$$

By assuming that the time-dependent bubble radii can be represented as $`R_1=R_{10}+e_1`$, $`R_2=R_{20}+e_2`$, and $`\left|e_1\right|R_{10}`$, $`\left|e_2\right|R_{20}`$, Eqs. (4) and (5) are reduced to the following linear formulae:

$$\ddot{e}_1+\omega _{10}^2e_1+\delta _1\dot{e}_1=\frac{p_{\mathrm{ex}}}{\rho R_{10}}\frac{R_{20}^2}{R_{10}D}\ddot{e}_2,$$
(6)
$$\ddot{e}_2+\omega _{20}^2e_2+\delta _2\dot{e}_2=\frac{p_{\mathrm{ex}}}{\rho R_{20}}\frac{R_{10}^2}{R_{20}D}\ddot{e}_1,$$
(7)
where

$$\omega _{j0}=\sqrt{\frac{1}{\rho R_{j0}^2}\left[3\kappa P_0+(3\kappa 1)\frac{2\sigma }{R_{j0}}\right]}\mathrm{for}j=1,2,$$

are the partial natural (angular) frequencies of the bubbles, and

$$\delta _j=\frac{4\mu }{\rho R_{j0}^2}\mathrm{for}j=1,2,$$

are the damping coefficients. (In general, the damping coefficients are determined by the sum of viscous, radiation, and thermal damping ; we however use the above setting to simplify the following discussion.)

The system of equations (6) and (7) corresponds to that derived by Shima when $`\sigma =0`$, $`\delta _1=0`$, and $`\delta _2=0`$ , and to that by Zabolotskaya when $`\sigma =0`$ . This kind of system of differential equations, called the coupled-oscillator model or self-consistent model, has been used repeatedly to analyze acoustic properties of coupled bubbles , and has been proved to be a useful model that can provide a qualitatively accurate result despite its simple form.

Let us analyze in detail the transition frequencies of the linear coupled system. By assuming that the external sound pressure at the bubble position is written in the form of $`p_{\mathrm{ex}}=P_a\mathrm{sin}\omega t`$, a harmonic steady-state solution of the linear coupled system is given as

$$e_1=K_1\mathrm{sin}(\omega t\varphi _1),$$
(8)
where

$$K_1=\frac{P_a}{R_{10}\rho }\sqrt{A_1^2+B_1^2},$$
(9)
$$\varphi _1=\mathrm{tan}^1\left(\frac{B_1}{A_1}\right),$$

with

$$A_1=\frac{H_1F+M_2G}{F^2+G^2},$$
(10)
$$B_1=\frac{H_1GM_2F}{F^2+G^2},$$
(11)
$$F=L_1L_2\frac{R_{10}R_{20}}{D^2}\omega ^4M_1M_2,$$

$$G=L_1M_2+L_2M_1,H_1=L_2+\frac{R_{20}}{D}\omega ^2,$$

$$L_1=(\omega _{10}^2\omega ^2),L_2=(\omega _{20}^2\omega ^2),$$

$$M_1=\delta _1\omega ,M_2=\delta _2\omega .$$

Here the solution for only bubble 1 is shown. The phase difference $`\varphi _1`$ becomes $`\pi /2`$ when

$$A_1=0.$$
(12)
It should be noted that a case in which both $`A_1`$ and $`B_1`$ become zero does not exist. From Eqs. (10) and (11), one obtains

$$A_1^2+B_1^2=\frac{H_1^2+M_2^2}{F^2+G^2}\left(\frac{P_a}{R_{10}\rho }\right)^2.$$

The numerator of this equation always has a nonzero value, since $`M_2>0`$; this result denies the existence of a case where both $`A_1=0`$ and $`B_1=0`$ are true. Also, it should be noted that $`F^2+G^2`$ appearing in the denominator of Eq. (10) always has a nonzero value. When $`G=0`$, $`F`$ is reduced to

$$F=\frac{M_2}{M_1}L_1^2\frac{R_{10}R_{20}}{D^2}\omega ^4M_1M_2.$$

This has a nonzero, negative value because $`M_2L_1^2/M_10`$, $`R_{10}R_{20}\omega ^4/D^2>0`$, and $`M_1M_2>0`$. This result means that no case exists where both $`F=0`$ and $`G=0`$ are true. As a consequence, Eq. (12) is reduced to

$$H_1F+M_2G=0.$$
(13)
In the following, we analyze this equation.

When the damping terms in Eq. (13) are negligible (but exist), one can easily obtain the solution for this equation. By assuming that $`M_10`$ and $`M_20`$, one obtains

$$H_1F0.$$

This equation can be rewritten into two independent equations:

$$FL_1L_2\frac{R_{10}R_{20}}{D^2}\omega ^4=0$$
(14)
and

$$H_1=L_2+\frac{R_{20}}{D}\omega ^2=0.$$
(15)
Equation (14) is the same as the theoretical formula given in Refs.  to derive the natural frequencies of a double-bubble system; namely, the transition frequencies given by this equation correspond to the natural frequencies. This equation predicts the existence of two natural frequencies,

$$\omega _{1\pm }^2=\frac{\omega _{10}^2+\omega _{20}^2\pm \sqrt{\left(\omega _{10}^2\omega _{20}^2\right)^2+4\omega _{10}^2\omega _{20}^2\frac{R_{10}R_{20}}{D^2}}}{2\left(1\frac{R_{10}R_{20}}{D^2}\right)},$$
(16)
and is symmetric; namely, it exchanges subscripts 1 and 2 (or 10 and 20) in this equation to reproduce the same equation. This means that the two bubbles have the same two natural frequencies. One of the natural frequencies converges to $`\omega _{10}`$ and the other to $`\omega _{20}`$ for $`D\mathrm{}`$, and the higher natural frequency rises and the lower one falls as $`D`$ decreases .

The solution of Eq. (15) given for the first time is

$$\omega _1^2=\frac{\omega _{20}^2}{1R_{20}/D}.$$
(17)
This converges to $`\omega _{20}^2`$ for $`D\mathrm{}`$, and increases as $`D`$ decreases. In contrast to Eq. (14), Eq. (15) is asymmetric ($`H_1H_2`$); namely, this serves to break the symmetry of the natural frequency mentioned above. The transition frequency given by Eq. (17) is not the natural frequency because it cannot be given by the natural frequency analysis performed in Refs. . (Even the other models used in, e.g.,  give only two natural frequencies.)

The results given above show that, when the damping effect is negligible and the radii of bubbles are not equal, the bubbles have three asymmetric transition frequencies; one of these, for $`D\mathrm{}`$, converges to the partial natural frequency of a corresponding bubble, while the remaining two converge to that of a neighboring bubble. (In the following, the former transition frequency is called “fundamental transition frequency” (FTF), and the latter two are called “sub transition frequencies” (STFs).) One of the STFs always increases as bubbles approach each other. The other STF decreases (increases) and the FTF increases (decreases) when the partial natural frequency of the bubble is higher (lower) than that of the neighboring bubble.

Figure 1 shows the transition frequencies for $`\mu 0`$ as a function of $`l=D/(R_{10}+R_{20})`$. Other parameters are set to $`\rho =1000`$ kg/m3, $`P_0=1`$ atm, $`\sigma =0.0728`$ N/m, and $`\kappa =1.4`$. Equations (16) and (17) are used in plotting those graphs. The radius of bubble 1 is fixed to 10 $`\mu `$m and that of bubble 2 varies from 5 $`\mu `$m to 10 $`\mu `$m. As discussed above, three transition frequencies changing with $`D`$ appear in each graph, except for the case of $`R_{10}=R_{20}`$, where only one decreasing transition frequency appears. (When the two bubbles have the same radii and pulsate in phase, Eq. (6) is reduced to

$$\left(1+\frac{R_{10}}{D}\right)\ddot{e}_1+\omega _{10}^2e_1+\delta _1\dot{e}_1=\frac{p_{\mathrm{ex}}}{\rho R_{10}}.$$

This equation has only one transition frequency  of

$$\omega _1^2=\frac{\omega _{10}^2}{1+R_{10}/D},$$
(18)
which converges to $`\omega _{10}^2`$ for $`D\mathrm{}`$ and decreases with decreasing $`D`$.) It is interesting to point out that the larger bubble, whose partial natural frequency is lower than that of the smaller bubble, has the highest transition frequency among all the bubbles.

Now we present numerical solutions of Eq. (13) to examine the influences of viscosity on transition frequencies. Figures 2$``$4 show the results obtained by using $`\mu =1.137\times 10^3`$ kg/(m s), which corresponds to the viscosity of water at room temperature. From those figures, we can observe that, as the viscous effect grows strong, i.e., the bubbles’ radii become smaller, the STFs vanish gradually from the large-distance region (and sometimes from the small-distance region), and only the FTF remains. In the case of $`R_{10}=1`$ $`\mu `$m, the STFs of a larger bubble disappear, and, when the bubbles are of similar size to each other, the FTF and the higher STF of a smaller bubble vanish in the small-distance region. In the case of $`R_{10}=0.1`$ $`\mu `$m, it is difficult to distinguish the STFs from the FTF since only a smooth curve, decreasing almost monotonically with decreasing $`l`$, appears. (The transition frequency of the smaller bubble remaining in the small-distance region may be the lower one of the STFs.) Those results show that the influence of the mutual interaction of the bubbles on the transition frequencies depends on the viscosity of the surrounding material, and that this interaction weakens as the viscous effect strengthens, i.e., the threshold of the distance for the appearance of the STFs is shortened.

Lastly, we briefly discuss the pulsation amplitudes of the bubbles. Figure 5 shows $`\alpha _j(P_a/R_{j0}\rho )A_j`$ and $`\beta _j(P_a/R_{j0}\rho )B_j`$ for $`R_1=10`$ $`\mu `$m, $`R_2=8`$ $`\mu `$m, $`l=3`$, and $`\mu =1.137\times 10^3`$ kg/(m s) as functions of $`\omega /\omega _{10}`$, where we set $`P_a=0.1P_0`$. The present setting is equivalent to that used for one of the graphs shown in Fig. 2 except for $`l`$. In each graph shown in Fig. 5, we observe three reversals of the sign of $`\alpha _j`$ and only two resonance responses, as expected. At the transition frequency that does not correspond to the natural frequency (the highest one of bubble 1 and the second highest one of bubble 2), no resonance response is obtained. This result confirms that the characteristic frequency given by Eq. (15) is not the natural frequency.

In summary, it was predicted in this Letter that a gas bubble interacting acoustically with a neighboring bubble has three transition frequencies that make the phase difference between the bubble’s pulsation and an external sound be $`\pi /2`$, while readymade theories predict only two natural frequencies. This present result shows that, in a double-bubble case, the phase reversal of a bubble (e.g., from in-phase to out-of-phase with an external sound) can take place not only at the bubble’s natural frequencies but also at some other frequency. The present theory for transition frequencies may be useful for understanding the reversal of the sign of the secondary Bjerknes force, which has so far been interpreted using only the natural frequencies . (This subject is discussed in a separate paper of ours .) Furthermore, the present results might also affect understandings of some related subjects such as acoustic localization  and superresonances , previous investigations of which have employed systems containing only identical bubbles. The present theory for a double-bubble system will be verified  by the direct numerical simulation technique proposed in Ref. , and will be extended to a theory for an $`N`$-bubble system , where $`N`$ denotes an arbitrary positive integer.


text:
# Statistical spectroscopic calculation of expectation values and spin-cutoff factors

## Abstract

Recently we proposed a formalism, based in nuclear statistical spectroscopy, for efficient computation of nuclear level density, or densities of states, through a sum of partitioned binomial functions (SUPARB). In this Letter we extend the formalism to the calculation of locally averaged expectation values, with specific application to spin-cutoff factors and the angular momentum dependence of the nuclear density of states.

In a recent Letter  we proposed a computationally efficient method to compute nuclear level densities, based in the ideas of nuclear statistical spectroscopy . While we refer interested readers to our original Letter for details, we briefly summarize here. Consider the density of states of a Hamiltonian $`\widehat{H}`$,

$$\rho (E)=\mathrm{tr}\delta (E\widehat{H}).$$
(1)
It is useful to partition the finite model space into subspaces, labeled by Greek letters $`\alpha ,\beta `$, etc., each with an associated projection operator $`P_\alpha `$. Then one can define partial or configuration densities:

$$\rho _\alpha (E)=\mathrm{tr}P_\alpha \delta (E\widehat{H}).$$
(2)
The total density is the sum of the partial densities. (NB:We always include $`2J+1`$ degeneracies and so are formally considering state densities.)

If the subspaces are single-particle configurations, e.g., $`(0d_{5/2})^4`$, $`(0d_{5/2})^2(1s_{1/2})^2`$, etc, then the partial, or configuration, moments up to fourth order for any system with $`Z`$ protons and $`N`$ neutrons in the valence space can be computed directly from the one+two-body matrix elements of $`\widehat{H}`$ . For any partition (configuration) $`\alpha `$ let $`d_\alpha =\mathrm{tr}P_\alpha `$ be the dimension of the subspace $`\alpha `$, and define the configuration average over the subspace to be $`\mathrm{}_\alpha d_{\alpha }^{}{}_{}{}^{1}\mathrm{tr}P_\alpha (\mathrm{})`$. Then $`\overline{E}(\alpha )=\widehat{H}_\alpha `$ is the (configuration) centroid, $`\gamma (\alpha )=(\widehat{H}\overline{E}_\alpha )^2_\alpha ^{1/2}`$ the configuration width, $`m_3(\alpha )=(\widehat{H}\overline{E}_\alpha )^3_\alpha /\gamma ^3(\alpha )`$ the scaled (dimensionless) third configuration moment, and $`m_4(\alpha )=(\widehat{H}\overline{E}_\alpha )^4_\alpha /\gamma ^4(\alpha )`$ the scaled fourth configuration moment.

With these moments in hand, we model the partial densities as binomial distributions, following a recent suggestion of Zuker . Starting with the binomial expansion of $`(1+\lambda )^N`$ and representing $`\left(\genfrac{}{}{0pt}{}{N}{k}\right)`$ with gamma functions, one can derives a continuous distribution,

$$\rho (E_x)=\lambda ^{E_x/ϵ}\frac{\mathrm{\Gamma }(E_{max}/ϵ+1)}{\mathrm{\Gamma }(E_x/ϵ+1)\mathrm{\Gamma }((E_{max}E_x)/ϵ+1)}$$
(3)
where $`E_x`$ is the excitation energy and $`E_{max}=\gamma (1+\lambda )\sqrt{N/\lambda }`$ is the maximum excitation energy in the binomial distribution. The binomial is appealing because one can easily compute the scaled third and fourth central moments: $`m_3=\frac{1\lambda }{\sqrt{N\lambda }}`$ and $`m_4=3\frac{4\lambda }{N}+\frac{1}{N\lambda },`$ allowing one to control the shape of the binomial through $`N`$ and $`\lambda `$.

To compute the density of states, we take the following steps: (1) We compute the configuration moments up to 3rd or 4th order. (2) We model the partial density for each configuration as a binomial. The binomial parameters $`N`$ and $`\lambda `$, as well as the overall energy scale and centroid, are fitted to the configuration moments. (3) The partial densities are summed to yield the total density of states. Because of these ingredients, we will refer to our approach as SUPARB (SUm of PARtitioned Binomials) state densities. We have shown that our method models well the density of states, when compared against exact calculations, and that one needs third configuration moments and occasionally, but not always, fourth configuration moments to achieve accurate results . Incidentally, none of this is special to the atomic nucleus. One could easily apply it to atomic electrons as well.

We now turn to the computation of expectation values of some general operators, $`\widehat{𝒪}`$, and define the energy-dependent locally averaged expectation value (LEV) as

$$\widehat{𝒪}(E)=\frac{\mathrm{tr}\widehat{𝒪}\delta (E\widehat{H})}{\mathrm{tr}\delta (E\widehat{H})}.$$
(4)
We also have configuration LEVs:

$$\widehat{𝒪}(E)_\alpha =\frac{\mathrm{tr}P_\alpha \widehat{𝒪}\delta (E\widehat{H})}{\mathrm{tr}P_\alpha \delta (E\widehat{H})}.$$
(5)
Our strategy for computing the total $`\widehat{𝒪}(E)`$ is the same as for the total density of states. In addition to the configuration moments, we compute the weighted averages $`\widehat{𝒪}_\alpha `$, $`\widehat{𝒪}(\widehat{H}\overline{E}_\alpha )_\alpha `$, and $`\widehat{𝒪}(\widehat{H}\overline{E_\alpha })^2_\alpha `$. These can be computed directly from one+two-body matrix elements similar to those for the moments .

The only tricky point is if $`\widehat{𝒪}`$ and $`\widehat{H}`$ do not commute; where does one insert the projection operator $`P_\alpha `$ in eqn. (5)? Because $`_\alpha P_\alpha =1`$ by the completeness of the projection operators, we conclude consistency is the only requirement, and take $`\mathrm{tr}P_\alpha \widehat{𝒪}\widehat{H}^k`$, $`k=0,1,2`$.

Assume within any subspace a quadratic energy dependence, that is,

$$\widehat{𝒪}(E)=𝒪_0+𝒪_1(E\overline{E})/\gamma +𝒪_2(E\overline{E})^2/\gamma ^2$$
(6)
(dropping for the moment the subspace label $`\alpha `$). Then

$`\widehat{𝒪}=𝒪_0+𝒪_2,`$ (7)
$`\widehat{𝒪}(\widehat{H}\overline{E})=\gamma (𝒪_1+𝒪_2m_3),`$ (8)
$`\widehat{𝒪}(\widehat{H}\overline{E})^2=\gamma ^2(𝒪_0+𝒪_1m_3+𝒪_2m_4).`$ (9)
Solving for coefficients,

$`𝒪_0={\displaystyle \frac{\widehat{𝒪}(m_4m_3^2)+m_3\widehat{𝒪}(\widehat{H}\overline{E})/\gamma \widehat{𝒪}(\widehat{H}\overline{E})^2/\gamma ^2}{m_4m_3^21}},`$ (10)
$`𝒪_1={\displaystyle \frac{(m_41)\widehat{𝒪}(\widehat{H}\overline{E})/\gamma +m_3\left(\widehat{𝒪}(\widehat{H}\overline{E})^2/\gamma ^2\widehat{𝒪}\right)}{m_4m_3^21}},`$ (11)
$`𝒪_2={\displaystyle \frac{\widehat{𝒪}(\widehat{H}\overline{E})^2/\gamma ^2m_3\widehat{𝒪}(\widehat{H}\overline{E})/\gamma \widehat{𝒪}}{m_4m_3^21}}`$ (12)
If we limit ourselves to only a linear dependence (assume $`𝒪_2=0`$), then

$`𝒪_0=\widehat{𝒪},`$ (13)
$`𝒪_1=\widehat{𝒪}(\widehat{H}\overline{E})/\gamma .`$ (14)
Note: in terms of the binomial parameters $`N,\lambda `$, one finds that $`m_4m_3^21=2(1N^1)`$.

The generalization to the full SUPARB case, summing over partitions, is easy:

$$\widehat{𝒪}(E)=\frac{\underset{\alpha }{}\widehat{𝒪}(E)_\alpha \rho _\alpha (E)}{\rho (E)}$$
(15)
In Figures 1-3 we illustrate our method, comparing to “exact” shell model calculations in full $`0\mathrm{}\omega `$ spaces. For $`sd`$-shell nuclides we used the Wildenthal USD interaction , and compared against direct diagonalization. Fig. 1 shows the LEV of $`\stackrel{}{Q}\stackrel{}{Q}`$ and $`S^2`$ for 20Ne, while Fig. 2 shows the LEV of $`J^2`$ (which will be important for the spin-cutoff factor below) for 22Na, 23Mg, and 32S.

Fig. 3 compares the LEV of $`J^2`$ for 48Cr and 54Fe in a full $`0\mathrm{}\omega `$ $`pf`$-shell calculation; here the “exact” calculation was through Monte Carlo sampling of path integrals. To avoid the well-known sign problem  we fitted a schematic multipole-multipole interaction to the $`T=1`$ matrix elements of the FPD6 interaction of Richter et al. . Clearly the full SUPARB calculation, with a quadratic dependence in each partition, works very well, and the quadratic is a significant improvement over the linear approximation. One could go to a cubic dependence, but computation of the necessary moments, while possible, would be very time consuming and there is not much room for improvement.

Now we apply our formalism to spin-cutoff factors for state densities. Define the $`J`$-dependent density as

$$\rho (E,J)=\mathrm{tr}\left(\delta (E\widehat{H})\delta (J(J+1)\widehat{J}^2)\right).$$
(16)
Traditionally this is factorized to

$$\rho (E,J)=\rho (E)(2J+1)\mathrm{\Omega }(J,E),$$
(17)
and one assumes a weighted Gaussian form for $`\mathrm{\Omega }`$:

$$\mathrm{\Omega }(J,E)=(2J+1)\frac{1}{4\sigma (E)\sqrt{2\pi }}\mathrm{exp}\left(\frac{(J+\frac{1}{2})^2}{2\sigma (E)^2}\right)$$
(18)
We have normalized $`𝑑J\mathrm{\Omega }(J,E)=1`$. Here $`\sigma (E)`$ is the “spin-cutoff” factor and is dependent on the energy. The spin-cutoff factor has been considered before in the context of nuclear statistical spectroscopy , but the moments were computed by a random sampling of representative vectors and the level densities were approximated by Hermite polynomials (which do not guarantee nonnegative densities).

Given the form (18) one finds that

$$\widehat{J}^2=3\sigma ^2\frac{1}{4}$$
(19)
(The factor of 3 is because we must include the $`2J+1`$ degeneracy in our traces.) Therefore, the problem of describing $`\rho (E,J)`$ with SUPARB reduces to computing $`\widehat{J}^2`$ as a function of energy, which were given in figures 2 and 3.

In figure 4 we plot the exact and SUPARB $`J`$-projected level densities (without the $`2J+1`$ degeneracy) of 32S, for $`J=0,1,2`$ and 8. The results are very good. Incidentally, one could compute the spin-cutoff factors and thus the $`J`$-projected level densities using Monte Carlo evaluation of path integrals just as easily, although this has not yet been done.

In summary, we have extended our previous SUPARB technique for state densities to expectation values of operators, including application to spin-cutoff factors. Our results look very good. Other possible application include estimating the contamination of spurious states in cross-shell calculations, and computation of total strengths and energy-weighted sum rules for transitions.

This work was performed under the auspices of the Louisiana Board of Regents, contract number LEQSF(1999-02)-RD-A-06; and under the auspices of the U.S. Department of Energy through the University of California, Lawrence Livermore National Laboratory, under contract No. W-7405-Eng-48.