text:
# Using Pulsars to Detect Massive Black Hole Binaries via Gravitational Radiation: Sagittarius A∗ and Nearby Galaxies

## 1 Introduction

Sazhin (1978) and Detweiler (1979) discussed the influence of long-wavelength (nanoHertz) gravitational radiation on the propagation of pulsar signals. Detweiler (1979) suggested a possible source of such radiation: binary massive black holes (MBHs) in distant galaxies. We have been engaged in a program to detect the stochastic background from the Universe of coalescing MBHs as well as to make estimates of the expected level. Here we consider the detection of gravitational radiation from the nearest objects.

We begin this inquiry by considering our Galactic Center (GC). There has been mounting evidence that the dark mass detected via proper motions of IR stars in the vicinity of Sagittarius A (SgrA) is a MBH (Eckart & Genzel 1997; Ghez et al. 1998; Maoz 1998; Ghez et al. 2000). Proper motion and absolute astrometry techniques lead to the identification of the compact non-thermal radio source Sgr A with the MBH (Backer & Sramek 1999; Reid et al. 1999; Menten et al. 1997). The electromagnetic emission of Sgr A may be either from the faint glow of matter being accreted on the MBH or from cooling in a small disk/jet system (see e.g. Narayan, Igumenshchev & Abramowicz 2000; Falcke 1996 ). At mm and short cm wavelengths where the effects of interstellar scattering are minimized the intrinsic source has been determined to be of order a few AU, less than 100 Schwarzschild radii of the MBH (Lo et al. 1998; Doeleman et al. 2001).

Recently, Zhao et al. (2000) reported quasi-periodic flux variations of Sgr A with a 106d period using observations at 1.3 and 2.0 cm from the VLA. The authors explored various models to account for the variation including the possibility that the periodicity is related to the orbit of a binary companion. The authors discount the binary scenario because the two holes would be easily resolved by high angular resolution VLBI imaging assuming that both were luminous. VLBI observations of Sgr A at 22-43 GHz reveal only a single source (e.g., Bower & Backer 1998; Lo et al. 1998; Doeleman et al. 2001 ). One could further argue that owing to the short lifetime for coalescence of an equal-mass binary MBH in Sgr A the system is unlikely. The residuals in proper motion measurements of Sgr A can also be used to place an orientation dependent limit on the mass of a dark companion to the Sgr A MBH. While binarity of the Sgr A MBH is an unlikely explanation for the flux variations reported by Zhao et al., we proceed in this paper to explore the detectability of the gravitational radiation from such a binary in millisecond pulsar (MSP) timing residuals.

The ratio of the hole masses now being measured in nearby galaxies (Magorrian et al. 1998; Merritt & Ferrarese 2001) to their distance is such that these objects are also candidate sources for detectable gravitational radiation if we make the binary hypothesis for them also. In this case we have no candidate period and are free to explore the limits on binary mass ratio at orbital periods where we are most sensitive, 2000 days.

In this article we first discuss in §2 the possible amplitude of perturbation of pulsar timing residuals by the conjectured Sgr A binary MBH, including a discussion of possible mass ratios. In §3 we present our recent observations of MSPs and available archival data. This is followed by a periodogram analysis that yields the best limit we can reach with current data sets. §5 discusses the possibility of detecting binary MBHs in nearby galaxies for which hole masses have been estimated. In our final section we summarize our conclusions.

## 2 Perturbation of Pulsar Timing by Gravitational Radiation from Sgr A\*

We use the superb formulation of the gravitational radiation from binary masses by Peters & Matthews (1963, hereafter PM63) to make detailed estimates of the amplitude $`h(\stackrel{}{r},t)`$ of the possible gravitational emission from Sgr A. We then follow the development by Sazhin (1978) and Detweiler (1979) of the influence of this radiation on the electromagnetic pulses emitted from a pulsar as it travels through space-time that is perturbed by $`h(\stackrel{}{r},t)`$. In short, pulse propagation through complete cycles of $`h(\stackrel{}{r},t)`$ have no net effect on the arrival time. There is only a perturbation of the arrival time by the incomplete cycles traversed at the pulsar and at Earth.

We use Eqn. 16 in PM63 to calculate the luminosity $`L_{GW}`$ of the gravitational wave (GW) from the assumed circular binary system in the GC

$$L_{\mathrm{GW}}=\frac{32}{5}\frac{m^5q^2}{a^5(1+q)^4}\frac{G^4}{c^5}=2.9\times 10^{43}\mathrm{erg}\mathrm{s}^1\frac{q^2}{(1+q)^4}$$
(1)
where $`m=2.6\times 10^6`$ $`M_{}`$ is the total mass of the systemGhez et al. (2000), $`q`$ is the mass ratio ($`q1`$), and $`a`$ is the semi-major axis. The numerical result uses Kepler’s Law $`a^3=GmP_{orb}^2/4\pi ^2`$ and the 106d orbital period of interest. For $`P_{orb}=106d`$, $`a=59\mathrm{AU}`$. Here and below we assume a circular orbit which is likely following the combined actions of dynamical friction and radiation.

From the energy density of a GW, which is $`U=c^2\dot{h}^2/32\pi G`$ (Eqn. 2 in PM63), we derive the dimensionless amplitude of the GW;

$$h=34\left(\frac{m^{1.67}G^{1.67}}{P_{orb}^{0.67}d}\right)\frac{q}{(1+q)^2}=6.3\times 10^{14}\frac{q}{(1+q)^2},$$
(2)
where $`d`$ is the distance to the emitter, 8 kpc. In this expression $`h`$ is averaged over all orientations of the observer relative to the plane of the binary orbit. Since we obtained this expression from the total energy density, which is the sum of the contributions from two polarizations, $`h_+`$ and $`h_\times `$, $`h`$ is actually the quadrature sum of these two polarizations. In order to find the dependence on inclination angle, $`i`$, we use the expression for the average power radiated per solid angle in PM63. This shows that the power radiated along the axis of the orbit is 8 times that for an edge-on view.

As discussed by Detweiler (1979) and others the dimensionless strain $`h`$ produces an apparent redshift in the pulsar frequency. A periodic source of GW then will produce a periodic shift in pulse arrival time from propagation through the gravitational radiation with an amplitude, $`\delta t`$, which is given by

$$\delta t\frac{hP_{\mathrm{GW}}}{2\pi }=22\mathrm{ns}\sqrt{(1+6\mathrm{cos}^2i+\mathrm{cos}^4i)}\frac{q}{(1+q)^2},$$
(3)
where $`P_{\mathrm{GW}}`$ represents the period of the gravitational wave, $`P_{orb}/2`$. The mass ratio factor is at most $`1/4`$. The angular factor ranges from 1 to 2.8. Therefore, $`\delta t`$ is less than 16 ns, and its average value over all solid angles is

$$\delta t46\mathrm{ns}\frac{q}{(1+q)^2},$$
(4)
which is 11 ns for $`q=1`$ and less for any other mass ratio.

Detweiler (1979) discusses the dependence of the GW signature in pulsar timing on the angle between the GW and the pulsar sightline. Pulse propagation times are perturbed by the GW owing to incomplete traversal of a cycle of the GW both at the pulsar as pulses are emitted and at the Earth upon pulse reception. In the plane wave approximation the resulting timing residual, $`\delta t`$, is

$$<\frac{\delta t}{P_{\mathrm{GW}}}>=\frac{1}{2\pi }\left[\frac{(1+\gamma )}{2}h\right](f(ct_\mathrm{r})f(ct_\mathrm{e}+\gamma l)).$$
(5)
where $`\gamma `$ is the cosine of the angle, $`\varphi `$, between the GC and the pulsar, where $`\varphi =0`$ is defined as the pulsar lying along the line of sight to the GC. $`f(t)`$ is the dimensionless phase term that comes from the fragments of the GW traversed at the emitter (e) and receiver (r) ends ($`f1`$). The times of emission and reception are $`t_\mathrm{e}`$ and $`t_\mathrm{r}=t_\mathrm{e}+l/c`$, respectively, and the factor $`\gamma l`$ is the projection of the pulsar distance $`l`$ along the GW propagation vector. Note that when $`\gamma =+1`$, there is no effect from the GW: a pulsar lying along the line of sight to the GC will experience no effect. The residual is also identically zero for $`\gamma =1`$ which describes electromagnetic waves (EMW) traveling in the opposite direction as the GW. The residual increases as the EMW and GW become perpendicular, and reaches a maximum just before they become parallel. The angle between PSR B1937+21 and Sgr A from Earth is 58$`^{}`$, and between PSR J1713+0747 and Sgr A is 37$`^{}`$. The signal from the two pulsars will therefore be diminished from our earlier estimate in Equation 4 by a factor $`(1+\gamma )/2=0.76`$ and 0.90, respectively. This factor adjusts our earlier estimate of the maximum possible effect from 16 ns to 14 ns, and our estimate of the average effect from 11 ns to 10 ns due to an equal-mass binary.

If either of the sources were much closer to the GC than the Earth, then $`h`$ in Eqn. 5 would need to be corrected for both the different wave amplitudes and the different angular factors. In our case all relevant distances (Earth to GC, J1713+0747 to GC, B1937+21 to GC) are 7-8 kpc, and are not known to better than 25% (Kaspi, Taylor, & Ryba 1994; Camilo, Thorsett, & Kulkarni 1994). Since the distances and therefore the phase factors are unknown, and the amplitudes of the two effects at the two ends are roughly equal, it is even possible that the two phase factors will be nearly the same, vastly diminishing the signal. The signal at emitter site will represent GW level $`10^4`$ years ago, probably not too different from today. The phase factor from the receiver end will produce a signature that is correlated with other pulsars, whereas the emission terms will be uncorrelated. Clearly observations with an array of pulsars at different angles and distances are critical to overcoming this “emission phase” noise.

We have used a mass ratio of $`q=1`$ to calculate the maximum signal we might expect. There are two arguments against a large mass ratio. One is based on evolutionary arguments and the other on proper motion observations of Sgr A. First, the lifetime for gravitational inspiral from a 106d orbit for $`q=1`$ is $`3\times 10^6`$ y; the lifetime increases with $`(1+q)^2/q,q1`$. After two galaxies merge the timescale for their central black holes to reach such an orbit is unknown, and could be from 30 million years to a Hubble time (e.g., Rajagopal & Romani 1995, Gould & Rix 2000 ). Toth & Ostriker (1992) rule out the possibility of a recent merger using models of disk heating via accretion of satellite galaxies. They demonstrate that no more than 4% of the mass of the galaxy could have been accreted within the last 5 billion years. Xu & Ostriker (1994) model the accumulation of a central MBH in a scenario with much less disk heating: the accretion of primordial $`10^6\mathrm{M}_{}`$ black holes. They show that a quickly accumulating MBH in the GC is actually not what we should expect; dynamical friction and gravitational radiation serve to eject massive objects from the center as well as accreting them. Galaxies such as ours, they conclude, are usually host to zero, one or two MBHs. In any case, a crude estimate of the likelihood that we happen to be living in the epoch of coalescence of a $`2.5\times 10^6\mathrm{M}_{}`$ 106d binary is $`3\times 10^6`$y/(age of Milky Way) or less than 0.1%; in other words it is unlikely that we are observing the pair of black holes at the moment of their coalescence.

The second argument uses the VLA and VLBA measurements of the proper motion of Sgr A (Backer & Sramek 1999; Reid et al. 1999) to place a limit on the mass ratio. The three reported VLBA measurements span almost exactly 7 cycles of the observed 106d quasi-periodicity. The middle point is offset in phase from the end points by 0.5. The scatter about the line connecting the VLBA measurements, about 0.5 mas, provides an estimate of the maximum separation between the observed mass and the center of mass of the system along the axis where the observational data are most sensitive. If we assume Sgr A is centered on the larger of the two masses, this places a limit on the mass ratio of $`q<(0.5\mathrm{mas})(8\mathrm{kpc})/59\mathrm{AU}=0.07`$. This limit is mainly in the EW direction where the VLBA data are most precise. Thus, the limit is strictly on $`q\mathrm{sin}(\theta )`$ where $`\theta `$ is the inclination of the orbital plane to the plane defined by the NS direction and the line of sight. Interestingly, the VLBA data nearly rule out the possibility that we are seeing the less massive member of the conjectured binary.

In conclusion, we estimate that the perturbation in pulsar timing residuals is no larger than $``$ 14 ns for the pulsars we consider here. In the following section, we present the current level of precision in our measurements, and discuss the feasibility of detecting a 14-ns sinusoidal amplitude in our timing residuals.

## 3 Observations

We have been conducting monthly observations at 0.43 GHz, 1.4 GHz and 2.4 GHz of an array of MSPs using the Arecibo Observatory 300 m telescope111The National Astronomy and Ionosphere Center Arecibo Observatory is operated by Cornell University under contract with the National Science Foundation. since December 1997. These data are used to make precision arrival time measurements for a variety of astrophysical goals. We used the Arecibo-Berkeley Pulsar Processor (ABPP), which is a multi-channel, coherent dispersion removal processor222‘coherent’ means that the dispersion is removed in the voltage domain prior to power detection. with 112-MHz total bandwidth capability.

In this paper we include only PSRs J1713+0747, B1855+09 and B1937+21 because their data sets span the longest time, 3.2 years, and they have the best timing precision. For PSRs J1713+0747 and B1855+09 we use 56-MHz bandwidth for observations at 1.4 GHz, and 112-MHz for 2.4 GHz, while for PSR B1937+21 we use 45-MHz overall bandwidth at 1.4 GHz and 55-MHz at 2.4 GHz. Calibrated total intensity profiles were formed from signals with orthogonal circular polarization. The profiles were then cross correlated with a template to measure times of arrival (TOAs) relative to the observatory atomic clock. Small errors in the observatory UTC clock, of order 1 $`\mu `$s, were corrected based on comparison of local time to transmissions from the Global Positioning System of satellites (GPS) using the Totally Accurate Clock receiver at the observatory. GPS time is then corrected to the TAI scale via publications from the Bureau International des Poids et Mesures (BIPM). The resulting TOAs are modeled for spin, astrometric, and, when relevant, binary parameters using the TEMPO software package333http://pulsar.princeton.edu/tempo. In this paper we use the residuals from the model to look for other effects.

After fitting for phase, spin period ($`P`$), period derivative ($`\dot{P}`$), right ascension, declination, and proper motion in the B1937+21 data, and additionally for the 5 Keplerian binary parameters in PSR J1713+0747 and PSR 1855+09, we have the residuals shown in Figure 1. Note that we did not fit for the Shapiro delay which has been measured in both PSRs J1713+0747 and B1855+09, but rather set the values to the best-fit values published (Kaspi, Taylor, & Ryba 1994; Camilo, Foster, & Wolszczan 1994). The weighted RMS of the day-averaged residuals are 0.35 $`\mu `$s, 0.53 $`\mu `$s, and 0.14 $`\mu `$s for J1713+0747, B1855+09, and B1937+21, respectively.

In our periodogram analysis (§4) we also use the Princeton Mark II and Mark III data published and made public by Kaspi, Taylor, & Ryba (1994, hereafter KTR94) . This was a biweekly monitoring program that spanned 8 years. KTR94 carefully removed dispersive effects from each of their TOAs on PSR B1937+21, which we have not done in the ABPP data. Due to our large bandwidth (45 MHz vs their 10 MHz) we achieve a similar level of precision with shorter data span.

PSR B1937+21 has been demonstrated by KTR94 to be unstable on time scales of about 5 years, but on the time scales over which we are interested here, 25-100 d, any instability it may have is below the noise level of our data, about 0.14 $`\mu `$s.

The timing measurements at Arecibo Observatory have a short-term precision of order 50 ns on the “best” millisecond pulsars with integrations of 3 minutes over a single hour. This suggests that with sufficient averaging and frequent observation we could detect the 14-ns amplitude perturbation discussed above. However, we have not achieved this level of precision over the full data span of our current data ( $``$ 3 y). We suspect that the discrepancy is the result of distortion of the average pulse profile by diffractive scintillation across our wide band (50-110 MHz). The pulse profile evolves with frequency over the band, and diffractive scintillation weights different parts of the band more heavily on different days. We are working on an algorithm to suppress this effect (Lommen & Backer 2001).

Nevertheless we proceed in §4 with a careful analysis of the current limits our data can place on the existence of a GC binary MBH.

## 4 Periodogram Analysis

The residual data shown in Figure 1 represent an irregularly spaced, sparsely sampled system in which we look for periodicities. As such, the method described by Cumming et al. 1999 (hereafter C99) for constructing and normalizing periodograms is ideal. C99 were searching for planets in spectroscopic velocity data from the Lick Observatory. The astronomical goal is very different, but the method is essentially identical. This periodogram approach is similar to the approach taken by Bailes, Lyne, & Shemar (1993) to look for planets around slow period pulsars.

We constructed periodograms for the residual timing data in the following way. We fit the data to the function $`A\mathrm{cos}\omega t+B\mathrm{sin}\omega t+C`$ by minimizing $`\chi ^2`$, at a range of frequencies, $`\omega =2\pi f`$, centered at the nominal GW frequency, $`f_0`$ of 1/53 d-1. The arbitrary offset ‘C’ is critical for sparsely sampled data as demonstrated by C99, even though the mean of the residuals is fit by the period polynomial. The range of frequencies we considered was dictated by the width of the periodogram feature detected by Zhao, Bower, & Goss (2001) which was approximately unity for the 2-cm data. We chose to look at a frequency range, $`\mathrm{\Delta }f=\frac{1}{2}f_0`$, which is 0.014 to 0.024 d-1 corresponding to a range of periods from 42 to 71 days. The frequency range considered in the analysis is very important as the smaller the range considered for detection, the more sensitive the measurement will be. In order to sample the possible periodicities completely, the approximate size of the steps by which we needed to sample this frequency range is $`\frac{1}{S}`$, where $`S`$ is the time span of the data set. Since we are using data sets of different spans (namely the Kaspi data which is 8 years long, and the post-upgrade data which is 3 years long) we chose a frequency spacing that would slightly oversample the Kaspi data, and 4 times oversample the post-upgrade data: stepsize $`=2\times 10^4`$ d-1.

In Figure 2 we use the following formalism, taken from C99, to plot the periodogram power $`z(\omega )`$ as a function of $`\omega `$.

$$z(\omega )=\frac{(Nm3)}{2}\frac{(\chi ^2\chi ^2(\omega ))}{\chi ^2(\omega _0)}$$
(6)
where $`N`$ is the number of degrees of freedom of the data, $`m`$ is the number of independent fit parameters, $`\chi ^2`$ is the weighted sum of the squares of the original residuals, $`\chi ^2(\omega )`$ is the weighted sum of the squares of residuals after the periodicity fit is included and $`\omega _0`$ is the frequency that gives the lowest $`\chi ^2(\omega )`$. $`Nm3`$ is the number of degrees of freedom of the periodicity fit and corresponds to a residual fit using the $`m`$ parameters in the pulsar model, and 3 additional parameters corresponding to fitting $`A`$, $`B`$, and $`C`$. Therefore $`z(\omega )`$ is the amount by which the $`\chi ^2`$ is reduced by adding the periodic signal to the data, normalized by the best fit $`\chi ^2`$. Note that $`z(\omega )`$ is normalized such that $`z(\omega )=1`$ for no sinusoidal signal present in the residuals.

We decided not to include the data from PSR B1855+09 in our analysis due to the higher RMS of its residuals compared to the other pulsars. The normalization of the statistic $`z(\omega )`$ allows us to average the three periodograms from the 3 remaining data sets to acquire the final periodogram shown in Figure 2. We also repeated the analysis using only PSR B1937+21 data and achieved very similar results.

Is there any significant detection of gravitational radiation in the periodogram shown in Figure 2? To answer this question we use Monte Carlo simulations to determine whether the peak at 60.5 d, with a value of 2.5, is a spurious effect of the noise. Following the method described in C99, we created 400 sets of simulated residuals, with no periodic signal, but with identical statistics to the real data. We asked what percentage of these 400 realizations would conspire to produce a peak between 42 and 71 days, as high or higher than 2.5. This percentage is called the “false-alarm probability”, and if it is higher than some threshold the detection is spurious. The threshold must be determined by the specifics of the problem, and can be anywhere from $`10^5`$ to 0.1. C99, for example, used a threshold of $`10^5`$ because there were systematics in their data that would imitate a signal at higher values. Creating identical statistics to the original data set proved to be challenging and also crucial to producing a meaningful false-alarm probability. In all cases we replicated the sampling of the original data. When we merely generated random numbers with a gaussian noise distribution and the same RMS as the original data, the false-alarm probability was 100%. This is due to the “redness” of the variation in the residuals, i.e., neighboring residuals have much smaller RMS with respect to each other than the overall RMS of the data set. Therefore by creating gaussian noise with a particular RMS we were creating a data set with much higher RMS between neighboring points than was present in the original set. The redness of pulsar timing residuals is commonly called “timing noise” and has been studied and discussed by a number of authors, and will not be discussed here (Cordes & Helfand 1980; Arzoumanian et al. 1994; Kaspi, Taylor, & Ryba 1994). To account for this low-frequency variation we fit for 3 additional period derivatives in each data set. This is identical to removing a 5th order polynomial from each data set. We are confident that this additional fitting would not remove any signal at a 53 day period since a 5th order polynomial has 5 zero-crossings, and our data sets are all at least 3 years long. We may, therefore, be removing variations as short as (3 yr/5) $``$ 200 days but no shorter. Additionally, to fabricate the data we merely randomly reordered the real data, i.e. randomly assigned real residuals to the wrong dates, to assure that the statistics were identical. The final false alarm probability for producing the peak shown in Figure 2 was 23%, and thus the detection is spurious. To be sure we are not washing out features in the power-spectrum present in any one data set we have performed this analysis on each data set individually and obtained very similar results.

In order to find the minimum amplitude detectable in our data we used Monte Carlo simulation in the following fashion. We simulated timing residuals with the same statistical properties as our real data using the same method described above, and injected a signal at a 53-day period at a variety of amplitudes. The amplitude that we designate as our ‘limit’ corresponds to the amplitude where 99% of the random realizations of simulated data produced a power amplitude larger than that which we actually measure. This amplitude, $`R=0.15\mu `$s, represents the minimum amplitude we can detect with our current data.

## 5 Nearby Massive Dark Objects

What is the probability that massive dark objects (MDO) in other galaxies are effecting our timing residuals? Magorrian et al. (1998) determined the masses or mass limits of 29 MDOs in nearby galaxies using HST photometry and ground-based spectroscopy of the galaxies. Several independent studies have demonstrated that these estimates are systematically high(van der Marel 1999; Wandel 1999; Ferrarese & Merritt 2000). van der Marel (1999) suggests the discrepancy is due to the assumption of velocity isotropy made by Magorrian et al. Merritt & Ferrarese (2001) improved the mass determination by using the extremely tight relationship between MDO and the velocity dispersion of the host bulge. In our analysis we use the Magorrian et al. sample of galaxies with the updated masses by Merritt & Ferrarese (2001).

We choose an advantageous orbital period for this discussion, 2000 days, where our data place the most stringent limit, and the lifetimes of the orbital systems are longer. The wavelength of the GW emitted from a 2000-day system is 1000 days, which is the length of our data set, so we must be careful that our pulsar model fitting would not remove the signal we wish to detect. A third period derivative fit or higher order would do so. Consequently, in order to obtain a meaningful limit we did not fit for the additional 3 derivatives of the period as we did for the previously described periodogram analysis, but rather used the residuals shown in Figure 1 which include only a single period derivative. By the same technique described in §4 we determined that the detectable amplitude at this period is 170 ns. In Table 2 we show the amplitude of the effect on the timing residual from each of the objects studied by Merritt & Ferrarese (2001) assuming the MDOs are equal-mass binaries with 2000-day orbits. Only those which produce a residual amplitude greater than 10 ns are tabulated. The amplitude shown is averaged over all solid angles. The amplitudes are as large as 850 ns which would be detectable in our data. However, the probability of detecting such an object during coalescence is proportional to its lifetime. Lifetime, $`\tau `$, scales with $`q`$, $`M`$, and $`P_{orb}`$ as follows.

$$\tau =8.0\times 10^4\mathrm{y}\frac{(1+q)^2}{q}\left(\frac{m}{10^9\mathrm{M}_{}}\right)^{5/3}\left(\frac{P_{orb}}{2000d}\right)^{8/3}$$
(7)
With these larger mass systems, we gain amplitude of the gravitational wave ($`m^{5/3}`$), but lose reasonable lifetime of the system by the same factor. Fortunately, lifetime also scales strongly with orbital period in our favor ($`P_{orb}^{8/3}`$), so we have much to gain by looking at longer periods. We place upper limits on q for these objects in Table 2 just as we have for the GC using the 170-ns detectable amplitude. These limits are as small as $`q0.06`$, which correspondingly give much longer lifetimes for the systems as shown in the last column. Objects with a reasonable lifetime ($``$1 Gyr) have a $`P_{orb}`$ too large to be detectable by the pulsar measurements.

## 6 Conclusion

Gravitational radiation of an equal-mass 2.5 $`\times 10^6`$ $`\mathrm{M}_{}`$ black hole binary at the GC would produce a periodicity in pulsar arrival times of order 10 ns. While this is an order of magnitude below the limits of present data presented in this paper, in the future a special observing effort might reach such a detection level. However, published proper motion measurements of Sgr A place a limit on the mass ratio of any such binary of 20:1 for low inclinations and EW orientation. In this case the gravitational radiation effects would be below detection limit in pulsar timing. A small mass ratio is also likely given consideration of the coalescence time scale and the absence of any evidence of a recent large mass accretion event in the GC.

The known MDOs in nearby galaxies, if binary MBHs with orbital periods around 2000d, would produce a larger signal, up to $``$$`\mu `$s, than that estimated for Sgr A. However the lifetimes to gravitational radiation inspiral for such binaries are shorter than the already short lifetime of Sgr A and therefore lower the probability that we are seeing them in this phase of evolution. With a number of such objects in existence, the probability increases that at least one of them is still in a ‘young’ binary state, and might be seen in pulsar timing residuals. Maintaining of precision pulsar monitoring programs with long and continuous coverage is important for the future of such detection efforts. The “Pulsar Timing Array”, our precision millisecond pulsar timing program, extends the work we have described here to include the entire ensemble of MBH-MBH systems in the universe. This ensemble produces a stochastic GW background at a level that we can detect. This measurement will place important constraints on the origin and evolution of MBHs.

We are grateful to Andrew Cumming and Yuri Levin for very useful discussions. We thank the Arecibo Observatory telescope operators for many days and nights of assistance, and Dunc Lorimer for pioneering remote observing at AO.


format: text

text:
# SPATIOTEMPORAL DYNAMICS OF PACING IN MODELS OF ANATOMICAL REENTRY

( 1 Division of Cardiology, Weill Medical College of Cornell University, New York, NY, 10021, USA
2 Centre for Condensed Matter Theory, Department of Physics, Indian Institute of Science, Bangalore 560012, India
3 Dept. of Physiology and Biophysics, Weill Graduate School of Medical Sciences of Cornell University, New York, NY, 10021
)

*Abstract-*Reentry around non-conducting ventricular scar tissue, which often causes lethal arrhythmias, is typically treated by rapid stimulation from an implantable defibrillator. However, the mechanisms of termination (success and failure) are poorly understood. To elucidate such mechanisms, we studied pacing of anatomical reentry in 1-D and 2-D excitable cardiac media models. Our results suggest that the existence of inhomogeneity in the reentry circuit is essential for pacing termination of tachycardia to be successful. Considering the role of such inhomogeneities may lead to more effective pacing algorithms.
*Keywords -*  Reentry, tachycardia, antitachycardia, pacing

I. Introduction

Trains of local electrical impulses are widely used to restore normal wave propagation in the heart during tachycardia - but this process is not always successful. Such “antitachycardia pacing” may inadvertently induce tachycardias or cause tolerated tachycardias to degenerate to more rapid and threatening tachyarrhythmias. It is known from clinical electrophysiology studies that pacing to terminate uniform ventricular tachycardia (VT) can cause tachycardia acceleration or degeneration to ventricular fibrillation (VF) . The underlying mechanisms governing the success or failure of a pacing algorithm are not yet clear. Understanding these mechanisms is essential, as a better knowledge of the processes involved in the suppression of VT through ventricular extrastimuli pacing can aid in the design of more effective therapies. The reentrant action potential may propagate around an inexcitable obstacle (“anatomical reentry”) or may reside in a region of the myocardium that is excitable in its entirety (“functional reentry”) . Although studies on the role of pacing in eliminating tachycardia have been done for functional reentry , there has been little work done for the case of anatomical reentry .

Several factors influence the ability of extrastimuli and/or rapid pacing to interact with VT. The most prominent are : (a) VT cycle length (for anatomical reentry this is determined by the length of the VT circuit around the obstacle), (b) the refractory period at the stimulation/pacing site and in the VT circuit, (c) the conduction time from the pacing site to the VT circuit, and (d) the duration of the excitable gap. As there are a number of conditions to be satisfied before reentry is successfully terminated, a single extrastimulus is rarely sufficient. Multiple stimuli are often used - where the first extrastimulus is believed to “peel back” refractoriness to allow the subsequent extrastimuli to interact with the circuit more prematurely than was possible with only a single extrastimulus.

There has been substantial work on computer modeling of pacing termination of tachyarrhythmias in a one-dimensional ring . The termination of reentry in such a geometry (which is effectively that of the reentry circuit immediately surrounding an anatomical obstacle) occurs in the following manner. Each extrastimulus splits into two branches while traveling around the reentry circuit. The retrograde branch (proceeding opposite to the direction of the reentrant wave) ultimately collides with the reentrant excitation, causing mutual annihilation. The anterograde branch (proceeding along the direction of the reentrant wave) can, depending on the timing of the stimulation, lead to either resetting, where the anterograde wave becomes the new reentrant wave, or termination of reentry, where the anterograde wave gets blocked by the refractory end of the original reentrant wave. From continuity arguments, it can be shown that there exists a range of extrastimuli phases and amplitudes that leads to successful reentry termination . Unfortunately, the argument is essentially applicable only to a 1D model - the process is crucially dependent on the fact that the pacing site is on the reentry circuit itself. In reality, however, it is unlikely that the pacing site will be so fortuitously located.

In this study, we examine the dynamics of pacing from a site some distance away from the reentry circuit. One then has to consider the passage of the extrastimulus from the pacing site to the reentry circuit. Because the reentrant circuit is the origin of rapid outwardly propagating waves, the extrastimulus will be annihilated before it reaches the circuit under most circumstances. Successful propagation to the reentry circuit will require utilising multiple extrastimuli to “peel back” refractory tissue incrementally until one successfully arrives at the circuit. Further, the anterograde branch must be blocked by the refractory tail of the reentrant wave. However, as outlined in detail in the Discussion section of our paper, this is extremely unlikely to happen in a homogeneous medium. As shown by our simulation results, the existence of inhomogeneity in the reentry circuit seems essential for pacing termination of VT.

II. Methodology

Modified Fitzhugh-Nagumo type excitable media model equations of ventricular activation were used to simulate the interaction of extrastimuli with the VT circuit around an anatomical obstacle. For the 1D model we used the Panfilov model  defined by the two equations governing the excitability $`e`$ and recovery $`g`$ variables:

$$\begin{array}{ccc}e/t\hfill & =\hfill & ^2ef(e)g,\hfill \\ g/t\hfill & =\hfill & ϵ(e,g)(keg).\hfill \end{array}$$
(1)
The function $`f(e)`$, which specifies fast processes (e.g., the initiation of the action potential) is piecewise linear: $`f(e)=C_1e`$, for $`e<e_1`$, $`f(e)=C_2e+a`$, for $`e_1ee_2`$, and $`f(e)=C_3(e1)`$, for $`e>e_2`$. The function $`ϵ(e,g)`$, which determines the dynamics of the recovery variable, is $`ϵ(e,g)=ϵ_1`$ for $`e<e_2`$, $`ϵ(e,g)=ϵ_2`$ for $`e>e_2`$, and $`ϵ(e,g)=ϵ_3`$ for $`e<e_1`$ and $`g<g_1`$. We use the physically appropriate parameter values given in Ref. . For the 2D studies, in addition to the Panfilov model we used two other two-variable modified Fitzhugh-Nagumo type models: the Aliev-Panfilov model , which simulates the restitution property of cardiac tissue, as well as the model proposed by Kogan et al , which simulates the restitution and dispersion properties of cardiac tissue.

We solve the models by using a forward-Euler integration scheme. We discretise our system on a grid of points in space with spacing $`\delta x=0.5`$ dimensionless units and use the standard three- and five-point difference stencils for the Laplacians in spatial dimensions $`d=1`$ and 2 respectively. Our spatial grid consists of a linear lattice with $`L^{}`$ points or a square lattice with $`L\times L`$ points; in our studies we have used values of $`L^{}`$ ranging from 500 to 1000 and $`128L256`$. Our time step is $`\delta t=0.022`$ dimensionless units. We define dimensioned time $`T`$ to be $`5`$ ms and $`1`$ spatial unit to be $`1`$ mm.

Our initial condition is a wave front at some point in the medium with transient conduction block on one side to permit wavefront propagation along a single direction only. The conduction block is realised by making the corresponding region refractory. This initiates a re-entrant wave which roughly corresponds to an anatomically-mediated tachycardia. Extrastimuli are then introduced from a pacing site. In the 1D model, the ends have periodic boundary conditions - such that the medium represents a ring of cardiac tissue around an inexcitable obstacle. In the 2D models, for top and bottom ends we use no-flux (Neumann) boundary conditions since the ventricles are electrically insulated from the atria (top) and the waves converge and annihilate at the region around the apex (bottom). For the sides, we use periodic boundary conditions. Two square patches of inexcitable region (whose conductivity constant is zero) are placed about the centre of the simulation domain with a narrow passage for wave propagation between them (Fig. 1). No-flux boundary conditions are used at the interface between the active medium and the inexcitable obstacle. We use this arrangement to represent anatomical obstacles (e.g., scar tissue) around which “figure-of-eight” reentry can occur. This kind of pattern has been seen previously in actual cardiac tissue (e.g., in post healed myocardial infarction tissue ).

III. Results

In the 2D simulations we looked at the effectiveness of different pacing algorithms in terminating figure-of-eight reentry. The pacing was simulated by a planar wave initiated at the base of the simulation domain. This represents a wave that has propagated a fixed distance up the ventricular wall from a local electrode at the apex - similar to a wave initiated at the tip of an inverted cone before propagating up its walls. Two different pacing patterns were applied. The first used pacing bursts at a cycle length which was a fixed percentage of the tachycardia cycle length (determined by the size of the obstacles). This corresponds to “burst pacing” used in implantable cardiac defibrillators (ICD) . The pacing cycle length as well as the number of extrastimuli in each pacing burst (4-8) were varied to look at their effects on the reentrant wave. The second was decremental ramp pacing in which the interval between successive extrastimuli was gradually reduced (corresponding to “ramp pacing” used in ICDs ).

For a homogeneous cardiac medium, both ramp and burst pacing were unsuccessful at terminating reentry (see Discussion for mechanism). We also looked at the effect of anisotropy by making the conductivities along the vertical and the horizontal axes differ in a ratio of 1:0.3, which is consistent with cardiac tissue . As with the homogeneous case, in the anisotropic case none of the pacing methods we tried were successful in terminating reentry.

However, we did achieve successful termination in the presence of inhomogeneity in the model. This was done by placing a small zone ( = 7.5 mm) of slow conduction in the narrow channel between the two non-conducting patches. The conductivity in this region was 0.05 times smaller than the rest of the tissue. A 6-stimulus pacing burst was found to successfully block the anterograde branch of the extrastimulus traveling through the narrow channel - and hence terminate the reentry.

To understand the process of inhomogeneity-mediated termination we implemented a 1D model with a zone of slow conduction whose length and conductivity constant were varied to examine their effect on the propagation of the anterograde branch of the extrastimulus. The reentrant wave was started at a point in the ring (proximal to the zone of slow conduction) chosen to be the origin ($`x=0`$) at time $`t=T_0=0`$. At time $`t=T_1`$, the first extrastimulus was given at $`x=0`$ followed by a second extrastimulus at $`t=T_2`$ (Fig. 2). The first extrastimulus is able to terminate the reentry by itself if it is applied when the region to one end of it is still refractory - leading to unidirectional propagation. This is identical to the mechanism studied previously for terminating reentry in a 1D ring . However, in this study we are interested in the effect of pacing from a site away from the reentry circuit. In that case, it is not possible for the first extrastimulus to arrive at the reentry circuit exactly at the refractory end of the reentrant wave (see the Discussion section). Therefore, we consider the case when the first extrastimulus is unable to block the anterograde conduction by itself. We have used values of $`T_1`$ for which the first extrastimulus can give rise to both the anterograde as well as the retrograde branches.

Different values of coupling interval ($`T_1T_0`$) and pacing cycle ($`T_2T_1`$) were used to find which parameters led to block of the second anterograde wave.

IV. Discussion

Because the retrograde branch of the extrastimulus always mutually annihilates the reentrant wave, the success of pacing in terminating reentry around an anatomical obstacle depends on whether or not the anterograde branch of the extrastimulus is successfully blocked. Although apparently this is often achieved in practice (i.e., ICDs are often successful at anti-tachycardia pacing), as we have seen in our 2D model simulation, this is almost impossible to observe for homogeneous tissue unless the pacing site is located on the reentry circuit. Inhomogeneities appear to be necessary for termination, as outlined in the following simple mathematical argument.

Let us consider a reentrant circuit as a 1D ring of length $`L`$ with separate entrance and exit sidebranches (Fig. 3). Further, let the pacing site be located on the entrance sidebranch at a distance $`z`$ from the circuit. We use the entrance sidebranch as the point of spatial origin ($`x=0`$) to define the location of the wave on the ring. The conduction velocity and refractory period at a location a distance $`x`$ away (in the clockwise direction) from the origin is denoted by $`c(x)`$ and $`r(x)`$, respectively.

For a homogeneous medium, $`c(x)=c`$, $`r(x)=r`$ ($`c,r`$ are constants). Therefore, the length of the region in the ring which is refractory at a given instant is $`l=cr`$. For sustained reentry to occur, there must be an area of excitable tissue between the front and the refractory tail of the reentrant wave - i.e., the excitable gap. The condition for existence of this is $`L>cr`$. For convenience, associate $`t=0`$ with the time when the reentrant wave is at $`x=0`$ (i.e., the entrance sidebranch) (Fig. 3(a)). Let us assume that an extrastimulus is given at $`t=0`$. This extrastimulus will collide with the branch of the reentrant wave propagating out through the entrance sidebranch at $`t=z/2c`$ (Fig. 3(b)). The pacing site will recover at $`t=r`$ and if an extrastimulus is given again immediately it will reach the reentry circuit at $`t=r+(z/c)`$ (Fig. 3(c)). By this time the refractory tail of the reentrant wave will be at a distance $`x=z`$ away from the entrance sidebranch and the anterograde branch of the extrastimulus will not be blocked, i.e., it is impossible for the stimulus to catch up to the refractory tail in a homogeneous medium. This results in resetting of the reentrant wave rather than its termination.

Note that if the first extrastimulus is given at a time $`t<z/c`$, it reaches the reentry circuit before the arrival of the reentrant wave. As a result, the retrograde branch collides with the oncoming reentrant wave, while the anterograde branch proceeds to become the reset reentrant wave. It appears that only under the very special circumstance that the reentrant wave reaches the entrance sidebranch exactly at the same instant that the extrastimulus reaches the circuit, is there a possibility of the two waves colliding with each other leading to termination of reentry. However, if the circuit has any width, part of the reentrant wave is likely to survive and the wave will continue. Also such fine-tuning of pacing time is almost impossible in reality. As a result, pacing termination seems all but impossible if we consider the reentry circuit to be homogeneous.

The situation changes, however, if an inhomogeneity (e.g., a zone of slow conduction) exists in the circuit (Fig. 3(d)). In this case, the above argument no longer holds as the waves travel faster or slower depending on their location in the circuit. If the pacing site is proximal to the zone of slow conduction, then the pacing site may recover faster than points on the reentry circuit. As a result, extrastimuli may arrive at the circuit from the pacing site and encounter a region that is still recovering. This leads to successful block of the extrastimulus anterograde wave, while the retrograde wave annihilates the reentrant wave, resulting in successful termination. The success of pacing will depend on the location of the inhomogeneity. If pacing is performed distal to the zone of slow conduction, termination will be harder to achieve as the extrastimulus will have a longer distance to traverse before reaching the inhomogeneity, which will correspondingly have a longer time in which to recover. This is supported by electrophysiological studies of pacing in cardiac tissue .

Similar arguments may apply for other types of inhomogeneity. For example, existence of a region having longer refractory period than normal tissue will lead to the development of patches of refractory zones in the wake of the reentrant wave. If the anterograde branch of the extrastimulus arrives at such a zone before it has fully recovered, it will be blocked. In fact, a computer modeling study of the interruption of tachycardia through pacing has been performed assuming the existence of such a region on the reentrant circuit .

By providing simulation results of 1D and 2D models in which pacing is unsuccessful in terminating anatomical reentry in the absence of any inhomogeneities and providing a simple mathematical argument in support of this conclusion, this paper outlines the potential significance of inhomogeneities in cardiac tissue in the termination of VT through pacing.

V. Conclusion

There are obvious limitations of our study - the most prominent being the use of excitable media models of ventricular activity. These models do not incorporate details of ionic currents (unlike, e.g., the Luo-Rudy model) but rather lump such details into a single model parameter. Further, we have used 1D or 2D models in our study, even though the heart’s 3D structure is important in reentry. We have also assumed the heart to be a monodomain rather than a bidomain (which has separate equations for intracellular and extracellular space) - for the low antitachycardia pacing stimulus amplitude we believe this simplification to be justified.

Despite these limitations, the results presented here may apply to the case of pacing termination of anatomical reentry in the ventricle. We have tried to develop general (i.e., model-independent) mathematical arguments about the requirement that inhomogeneities should exist in the reentrant circuit for pacing to be successful in eliminating VT.

Acknowledgment

We thank Ken Stein for helpful discussions. This work was supported by the American Heart Association (#0030028N).

References


text:
# Explicitly solvable cases of one-dimensional quantum chaos

## Abstract

We identify a set of quantum graphs with unique and precisely defined spectral properties called regular quantum graphs. Although chaotic in their classical limit with positive topological entropy, regular quantum graphs are explicitly solvable. The proof is constructive: we present exact periodic orbit expansions for individual energy levels, thus obtaining an analytical solution for the spectrum of regular quantum graphs that is complete, explicit and exact.

PACS: 05.45.Mt,03.65.Sq

A quantum graph  is a network of vertices and bonds with a quantum particle moving along its bonds. An example of a graph with five vertices and six bonds is shown in Fig. 1. The wave function and the energy levels of a quantum particle on a graph are defined by the corresponding one-dimensional Schrödinger equation. Despite the apparent simplicity of the system, quantum graphs have proven to be a rich source of physical insight. From the mathematical point of view, the spectral properties of Schrödinger operators on graphs are highly nontrivial and have been widely investigated in the mathematical literature . Quantum graphs and networks have also been used to model various phenomena in different branches of physics and chemistry for more than 30 years. The most recent physical development appeared in a series of publications , where quantum graphs were studied in the context of quantum chaos.

It is easy to see that the behavior of a particle on a quantum graph is very complex. Each time the particle encounters a vertex $`v_i`$ of the graph, it can scatter with different probabilities in the forward or backward directions along any of the bonds emanating from the vertex. As a simple physical analogue of this system one may imagine a beam of light traveling along a network of optical fibers. At every joint of the fibers the light waves scatter in such a way that the total energy flux is conserved.

As a result of the multiple scattering possibilities at the vertices, the dynamics of a classical particle on a graph is very complex and the number of possible periodic orbits traced by the particle increases exponentially with their lengths. Consequently the topological entropy of the particle is positive and since the phase space of the system is bounded, the dynamics of the particle is ergodic. The classical chaoticity notwithstanding it was shown that several important spectral characteristics of quantum graphs, such as the density of states and the spectral staircase, can be obtained exactly in terms of periodic orbit expansion series .

The “wiring” of a quantum graph, i.e. the arrangement of bonds and vertices, is called the topology of the quantum graph. For any given graph topology there exists a wide variety of possible quantum graphs. The vertices, for instance, may be realized as simple hubs that redistribute the quantum flux into various channels; or we may place $`\delta `$ function potentials of various strengths at the vertices; or there can even be potential functions on the bonds. Because of this flexibility we anticipate that the topology of a quantum graph alone does not uniquely specify the graph’s spectral properties. We conjecture that there are several different types of quantum graphs, all chaotic in their classical limit, but each exhibiting unique and precisely defined spectral characteristics.

We start a rigorous classification of quantum graphs by defining regular quantum graphs. Although generically their classical limit is chaotic, we show that their spectrum is explicitly solvable analytically, state by state, via explicit periodic orbit expansions. This result is backed up by rigorous mathematical proofs  whose basic elements are presented below. To the authors’ knowledge this is the first time that the spectrum of a quantum chaotic system is obtained exactly and explicitly. In addition regular quantum graphs show a spectral gap at small energy spacings and a cut-off at large spacings. Both the size of the gap and the location of the cut-off are computed analytically.

The regular quantum graphs are different from a class of quantum graphs described in the literature  whose level spacing distributions, to a good approximation, exhibit features of the Gaussian orthogonal (unitary) ensemble . Thus there are at least two different types of quantum graphs with distinct spectral properties.

The periodic orbits of graphs are defined as the periodic connected sequences of bonds $`B_{ij}`$. Denoting by $`k`$ the wave number (momentum) of the particle, each periodic orbit contributes

$$S_{ij}=_{B_{ij}}k(x)𝑑x,$$
(1)
to the total action of a path traced by the particle. It turns out that the information contained in the totality of all the possible classical periodic orbits often allows one to reconstruct exactly certain quantities of a purely quantum nature. For example, according to recent results , the exact periodic orbit expansion for the density of states can be written as

$$\rho (k)\underset{j=1}{\overset{\mathrm{}}{}}\delta \left(kk_j\right)=\overline{\rho }(k)+\frac{1}{\pi }Re\underset{p}{}T_p\underset{\nu =1}{\overset{\mathrm{}}{}}A_p^\nu e^{i\nu S_p(k)}.$$
(2)
Here $`\overline{\rho }(k)`$ is the average density of states, $`\nu `$ is the repetition index, $`T_p=S_p(k)/k`$, and $`S_p`$, $`A_p`$ are correspondingly the action and the weight factor of the prime periodic orbit labeled by $`p`$. We assume in what follows that the system is scaling . This means that the actions of the periodic orbits are proportional to the wave number,

$$S_p(k)=S_p^0k,$$
(3)
where $`S_p^0`$, a constant, is the reduced action. In this case $`T_p=S_p^0`$ and $`A_p`$ are $`k`$-independent constants. We define the total reduced action $`S_0`$ of the graph, $`S_0=(k𝑑x)/k`$, where the integral is extended over all of the bonds of the graph. The scaling assumption is not an artificial restriction. It occurs, for instance, in atomic physics where scaled spectroscopy is now a common experimental technique . In addition, scaling quantum systems of this kind are the analogues of certain electromagnetic ray-splitting systems, flat metal cavities partially filled with a dielectric substance such as Teflon .

The formal description of a quantum graph system proceeds as follows . On a bond $`B_{ij}`$ connecting the vertices $`v_i`$ and $`v_j`$, the wave function of a quantum particle is defined by the one-dimensional Schrödinger equation which may include potentials on the bonds . At every vertex $`v_i`$ the wave function satisfies the usual boundary conditions of continuity and flux conservation. The consistency of the boundary conditions at every vertex of the graph naturally yields the global quantization conditions that determine the momentum eigenvalues $`k_n`$. A simple and elegant method based on the scattering quantization approach was presented in , where the quantization condition is given in the form

$$det[1S(k)]=0.$$
(4)
Here $`S(k)`$ is the scattering matrix , written explicitly in terms of the connectivity matrix  of the graph. It can be shown that the modulus of the complex function (4) is a trigonometric function of the form

$$\mathrm{cos}(S_0k+\pi \gamma )\mathrm{\Phi }(k)=0,$$
(5)
where

$$\mathrm{\Phi }(k)=\underset{i}{}a_i\mathrm{cos}(S_ik+\pi \gamma _i).$$
(6)
In the scaling case the $`a_i`$ are constants and $`S_i<S_0`$ are certain combinations of the classical actions (1). In general the functions $`\gamma `$ and $`\gamma _i`$ are bounded and tend to constant values for $`k\mathrm{}`$. In the scaling case they are $`k`$-independent constants.

In order to proceed we define regular quantum graphs. A regular quantum graph fulfils the condition

$$\alpha =\underset{i}{}a_i|<1.$$
(7)
It is convenient for the following discussion to assume $`\alpha >0`$, i.e. we exclude trivial graphs with $`\alpha =0`$. They are regular quantum graphs whose spectrum can be obtained trivially. For regular quantum graphs (5) can be solved formally for the momentum eigenvalues $`k_n`$. We obtain the following implicit solution for the roots of this quasi periodic function

$$k_n=\frac{\pi }{S_0}\left[n+\mu \gamma \right]+\frac{1}{S_0}\{\begin{array}{cc}\mathrm{arccos}(\mathrm{\Phi }_n),\hfill & \text{for }n+\mu \text{ even}\hfill \\ \pi \mathrm{arccos}(\mathrm{\Phi }_n),\hfill & \text{for }n+\mu \text{ odd,}\hfill \end{array}$$
(8)
where $`\mathrm{\Phi }_n=\mathrm{\Phi }(k_n)`$ and $`\mu Z`$, a fixed integer, is to be chosen such that $`k_1`$ is the first nonnegative solution of (5). Because of (7), the boundedness of the trigonometric functions in (6) and the properties of the $`\mathrm{arccos}`$ function, the second term of (8) is bounded away from 0 and $`\pi /S_0`$ and assumes only values between $`u`$ and $`\pi /S_0u`$, where $`0<u=\mathrm{arccos}(\alpha )/S_0<\pi /2S_0`$. Thus, for regular graphs the points

$$\widehat{k}_n=\frac{\pi }{S_0}(n+\mu \gamma +1),n=1,2,\mathrm{}$$
(9)
are not roots of (4). They serve as separators between root number $`n`$ and root number $`n+1`$. In fact (8) implies even more: the existence of finite-width “root-free zones” $`R_n=(\widehat{k}_nu,\widehat{k}_n+u)`$ surrounding every separating point $`\widehat{k}_n`$, where no roots of (4) can be found. Thus, roots of (4) can only be found in the “root zones” $`Z_n=[\widehat{k}_{n1}+u,\widehat{k}_nu]`$, subsets of the root intervals $`I_n=[\widehat{k}_{n1},\widehat{k}_n]`$. Since $`S_0`$ is the largest action in (5) and (6), it can be shown  that $`k_n`$ is the only root in $`Z_n`$. Therefore, in summary, there is exactly one root $`k_n`$ inside of $`Z_nI_n`$, and this root is bounded away from the separating points $`\widehat{k}_{n1}`$ and $`\widehat{k}_n`$ by a finite amount $`u`$.

The spectral properties of regular graphs discussed above allow us to draw several important conclusions. (i) Since there is exactly one root $`k_n`$ of (4) in $`I_n`$, this proves rigorously that the number of roots of (4) smaller than $`k`$ grows like $`\overline{N}(k)S_0k/\pi `$ (Weyl’s law), while the second term of (8) represents the fluctuations around the average. (ii) The existence of the root-free zones $`R_n`$ gives rise to a spectral gap of finite width $`g=2u`$ in the nearest neighbor spacing distribution of regular quantum graphs, i.e. no level spacings smaller than $`g`$ are ever found. This, for instance, precludes the possibility of degenerate eigenvalues. (iii) The existence of the separating points (9) together with the root-free zones imply the existence of a spectral cut at $`c=2(\pi /S_0u)`$, i.e. no level spacings greater than $`c`$ are ever found. (iv) The existence of the separating points (9) and the root-free zones $`R_n`$ are the key for obtaining an explicit and exact periodic orbit expansion for every root of (4). Multiplying both sides of (2) by $`k`$ and integrating from $`\widehat{k}_{n1}`$ to $`\widehat{k}_n`$ we obtain

$`k_n=\widehat{k}_n{\displaystyle \frac{\pi }{2S_0}}{\displaystyle \frac{1}{\pi }}\mathrm{Re}{\displaystyle \underset{p}{}}{\displaystyle \underset{\nu =1}{\overset{\mathrm{}}{}}}A_p^\nu {\displaystyle \frac{e^{i\nu S_p^0\widehat{k}_n}}{\nu }}\left\{(1e^{i\nu \omega _p})\left(i\widehat{k}_n{\displaystyle \frac{1}{\nu S_p^0}}\right)+{\displaystyle \frac{i\pi }{S_0}}e^{i\nu \omega _p}\right\},`$ (10)
where $`\omega _p=\pi S_p^0/S_0`$. Since all the quantities on the right-hand side of (10) are known, this formula provides an explicit representation of the roots $`k_n`$ of the spectral equation (4) in terms of the geometric characteristics of the graph. To our knowledge, this is the first time that the energy levels of a chaotic system are expressed explicitly in terms of a periodic orbit expansion. Previously, explicit formulae for individual energy levels were known only for integrable systems. In the context of periodic orbit theory, the energy levels of integrable systems are given by the Einstein-Brillouin-Keller (EBK) formula . However, apart from a few exceptional cases  EBK quantization is only of semiclassical accuracy.

For a generic chaotic system the energy levels can only be obtained indirectly as the singularities of the periodic orbit expansion of the density of states (Gutzwiller’s formula ), an implicit method which, in general, is only of semiclassical accuracy. Formula (10), on the other hand, shows that for regular quantum graphs every quantum level can be obtained individually, explicitly and exactly in terms of classical parameters.

In order to demonstrate that the class of regular quantum graphs is not empty we present an explicit example: the one-dimensional step potential shown in Fig. 2 (a) . With

$$V_0=\lambda E$$
(11)
we turn it into a scaling system. Physically this potential is realized for a rectangular microwave cavity partially loaded with a dielectric substance . The scaling step potential is equivalent to the scaling three-vertex chain graph shown in Fig. 2 (b), with two bonds $`L_1=b`$ and $`L_2=\beta (1b)`$, $`\beta =\sqrt{1\lambda }`$. In this case the spectral equation is given by

$$det\left[1S(k)\right]=\mathrm{sin}(Lk)r\mathrm{sin}\left[\left(L_1L_2\right)k\right]=0,$$
(12)
where $`L=L_1+L_2`$, and $`r=(1\beta )/(1+\beta )`$ is the reflection coefficient at the vertex $`v_2`$ between the two bonds. The solution of (12) is

$$k_n=\frac{\pi n}{L}+(1)^n\frac{1}{L}\mathrm{arcsin}\left[r\mathrm{sin}(\left(L_1L_2\right)k_n)\right],$$
(13)
where the second term, which gives the fluctuations around the average, is bounded by

$$\left|\delta k_n\right|\frac{1}{L}\mathrm{arcsin}(r)<\frac{\pi }{2L}.$$
(14)
Therefore, for the three-vertex chain graph with $`r<1`$, the separating points (9) are given by

$$\widehat{k}_n=\frac{\pi }{L}\left(n+\frac{1}{2}\right).$$
(15)
They mark the borders of the root intervals $`I_n`$ for all possible values of the bond lengths. For a given prime periodic orbit $`p`$, such as the ones shown in Fig. 2, the coefficients $`A_p`$ are given by

$$A_p=(1)^{\chi (p)}r^{\sigma (p)}(1r^2)^{\tau (p)/2},$$
(16)
where $`r`$ is the reflection coefficient at the middle vertex and $`\sigma (p)`$ and $`\tau (p)`$ are correspondingly the number of reflections off it, and transmissions through it. Since the reflection coefficient may be positive or negative depending on whether the particle scatters from the right or from the left, the factor $`(1)^{\chi (p)}`$ is needed to keep track of how many times it appears with a minus sign. Moreover $`\chi `$ keeps track of how many times the particle scatters off the walls of the potential, since each scattering event from the walls gives rise to a minus sign in the wave function. Using the explicit form of the spectral determinant (12), one obtains the explicit form of the eigenvalues for a quantum particle in the step potential in terms of classical (Newtonian and non-Newtonian ) periodic orbits.

In order to illustrate the validity of (10) we computed $`k_1`$, $`k_{10}`$ and $`k_{100}`$ including periodic orbits that experience up to five scattering events. We chose $`b=0.3`$ (see Fig. 2) and $`\lambda =1/2`$ (see (11)). We obtain $`k_1^{(5)}=4.1161`$, $`k_{10}^{(5)}=39.2866`$ and $`k_{100}^{(5)}=394.9477`$. This can be compared with the exact $`k`$ values given by $`k_1=4.107149`$, $`k_{10}=39.305209`$ and $`k_{100}=394.964713`$. In order to illustrate convergence, we also computed $`k_1`$, $`k_{10}`$ and $`k_{100}`$ including all periodic orbits that experience up to 20 scattering events. This amounts to including more than 100,000 periodic orbits and results in $`k_1^{(20)}=4.105130`$, $`k_{10}^{(20)}=39.305212`$ and $`k_{100}^{(20)}=394.964555`$. This also indicates that the convergence of (10) is not destroyed by keeping more and longer periodic orbits.

Additional examples of regular graphs are provided by the scaling “Manhattan potentials”. These potentials are obvious generalizations of the step potential shown in Fig. 2 (a) to piecewise constant potentials where the potential heights scale with the energy. Furthermore we checked explicitly that linear chain graphs with scaling $`\delta `$ function potentials at the vertices provide more examples of regular quantum graphs. In this case the strengths of the $`\delta `$ function potentials scale linearly with the momentum.

The key physical feature of these one-dimensional quantum systems, which permits the exact periodic orbit expansion (10) for the eigenvalues, is the rigidity of their spectra. For integrable systems, spectral rigidity is due to the “geometrical rigidity” of the periodic orbits, confined to integrable tori. In the case of quantum graphs, the geometrical structure of the periodic orbits is much more complicated. The complexity of the expansion (10) compared to the EBK formula reflects the geometrical complexity of the periodic orbits.

This paper profitted from a fruitful exchange of E-mails with Uzy Smilansky. Y.D. and R.B. gratefully acknowledge financial support by NSF grants No. PHY-9900730 and PHY-9984075; Y.D. and R.J by NSF grant No. PHY-9900746.

Figure Captions

Fig. 1: A generic (quantum) graph with five vertices and six bonds.

Fig. 2: (a) Simple step potential, a basic problem in one-dimensional quantum mechanics. Also shown are examples of Newtonian (“N”) and non-Newtonian (“NN”) periodic orbits used in the periodic orbit expansion of its energy eigenvalues (see Text). (b) Three-vertex chain graph corresponding to the step potential above.


text:
# Anomalous quantum confined Stark effects in stacked InAs/GaAs self-assembled quantum dots

submitted to Physical Review Letters

## Abstract

Vertically stacked and coupled InAs/GaAs self-assembled quantum dots (SADs) are predicted to exhibit a strong non-parabolic dependence of the interband transition energy on the electric field, which is not encountered in single SAD structures nor in other types of quantum structures. Our study based on an eight-band strain-dependent $`𝐤𝐩`$ Hamiltonian indicates that this anomalous quantum confined Stark effect is caused by the three-dimensional strain field distribution which influences drastically the hole states in the stacked SAD structures.
PACS: 78.67.Hc, 73.21.La, 71.70.Fk

Zero-dimensional semiconductor structure, like InAs /GaAs self-assembled quantum dots (SADs)  have attracted considerable attention because of the new physics  of a few electron systems and potential applications in optoelectronics . Recent experiment on Stark effect spectroscopy in SADs  has demonstrated the existence of an inverted electron-hole alignment due to the presence of gallium diffusion in InAs SADs, and established a relation between the Stark shift and the vertical electron-hole separation.

The theoretical interpretation of these experimental results is based on the assumption that the applied electric field can be treated by the second-order perturbation theory, which results in a quadratical dependence of the transition energy on the applied electric field ,

$$E(F)=E(0)+pF+\beta F^2$$
(1)
where $`p`$ is the built-in dipole moment and $`\beta `$ measures the polarization of the electron and hole states, i.e., the quantum confined Stark effect. While this relation is well satisfied in many quantum systems including single SADs , and quantum well structures , we show in this work that it is not valid for vertically coupled SAD structures  where the quantum confined Stark effect deviates significantly from its quadratic dependence on the electric field. The reason for this anomalous quantum confined Stark effect is due to the three-dimensional (3D) strain field distribution in the dots and in the coupling region, which controls the localization of hole states in the respective SADs, and their sensitivity to external field. The existence of this effect is important for basic condensed matter physics because it can not be inferred from a simple superposition of the electronic properties of single SADs. It is also promising for applications in optoelectronics because interband transition energies can be significantly modulated by electric fields in quantum dot lasers and other photonic devices.

The insets of Fig. 1 show schematically a single SAD structure and a system of two vertically coupled SADs that are truncated pyramids separated by a GaAs barrier of $`1.8`$ nm, with identical base $`17.4\times 17.4nm^2`$ and individual height $`3.6`$ nm. A positive electric field $`F`$ is directed from the top to the bottom of the structures. Fig. 1 shows the calculated ground state transition energies for the single dot and for the stacked structure, as functions of electric fields. The electron and hole states of the system are obtained from the Schrödinger equation in the framework of the envelop function formalism ,

$$(𝐇_{kp}+|e|Fz)\varphi =E\varphi .$$
(2)
Here $`\varphi =(\varphi _1,\varphi _2,\mathrm{},\varphi _8)`$ is the envelop eigenvector and $`𝐇_{kp}`$ is the $`kp`$ eight-band Hamiltonian which includes the effect of strain and piezoelectricity . The Hamiltonian is discretized on a three-dimensional grid as a large sparse matrix which is solved by Lanczos algorithm. This approach has been shown to be reliable, especially in the investigation of the inverted electron-hole alignment in SAD structures .

The single dot exhibits a nearly perfect quadratical dependence on the electric field which is referred to the conventional quantum confined Stark effect and has been observed in many other types of quantum structures . The maximum transition energy occurs at $`F50`$ kV/cm because this homogeneous InAs SAD has a positive electron-hole alignment, i.e., the hole is at the bottom of the dot for zero electric field. The stacked structure shows piecewisely quasi-linear dependence, with a turning point at a smaller negative electric field $`F22`$ kV/cm. The lower panel of Fig. 1 shows that this discrepancy from the conventional quantum confined Stark effect occurs over the whole range of electric field, i.e., for $`200`$ kV/cm $`F`$ $`200`$ kV/cm. In addition, the stacked structure exhibits a Stark shift of one order of magnitude stronger than the single dot. We also notice that the maximum transition energy for both systems are roughly identical, because in the stacked structure, despite the downshift of the electron states due to tunnel coupling, the hole states are characterized by an up-shift in energy due to the strain distribution (see Fig. 2). Similar simulations have been performed for stacked structures with different sizes, different thicknesses of the coupling region, and with or without wetting layers, where the same kind of anomalous quantum confined Stark effect has been found .

Fig. 2 shows the energy levels of the stacked structure, for both electrons and holes, as a function of electric fields. All the energies are given in reference to the top of the valence bands of GaAs. For the sake of showing the fine structure of the hole spectra near zero electric field, we use an energy scale twice as small as for electrons. For comparison, the energy levels of the single dot, are shown in dash-dotted lines. It is clearly seen that, due to the coupling between the stacked structures, the double dot system exhibits much richer structure in its energy spectrum.

In the conduction band, the energy spectrum of the single dot shows the ground 1s state, the two nearly degenerate excited 2p states, and the 3d-like state with weak sensitivity to the electric field. The coupled dot system results in bonding and antibonding states originating from these states, that in the absence of electric field, are identified as $`1s^+,1s^{},2p^+,2p^{},2s^+,3d^+,\mathrm{}`$ where $`+()`$ denotes bonding (antibonding) states. Except for the ground state, all the excited states are seen to have crossings or anticrossings with other states, which reorders the states at high fields. In the valence bands of the single dot, the four hole levels ($`h_1^s`$, $`h_2^s`$, $`h_3^s`$ and $`h_4^s`$) are seen to show quadratical dependence on the electric field.

The most dramatic feature in the valence bands of the stacked structure is the quasi-linear dependence of the hole levels on the electric field for intermediate and strong field intensities, especially the ground hole state $`h_1^d`$. In addition, we could not find any bonding or antibonding hole states for these structure as the ground hole state $`h_1^d`$ and first excited state $`h_2^d`$ are localized in the bottom dot and the top one, respectively. However, at small electric fields, the hole levels show small fluctuations due to mutual anticrossings. It is also noticed that the magnitude of the Stark shift for holes in the stacked structure is significantly larger than in the single dot, which is responsible for the large shift in the transition energy shown in Fig. 1.

The main reason for the hole states to behave so different in the stacked structure than in the single dot is traced in the 3D strain field distribution. Fig. 3(a) shows the profiles and the corresponding contour plots (insets) for the hydrostatic (H) and the biaxial (B) components of the strain field, which are defined as

$`H=ϵ_{xx}+ϵ_{yy}+ϵ_{zz}`$ (3)
$`B^2=(ϵ_{xx}ϵ_{yy})^2+(ϵ_{yy}ϵ_{zz})^2+(ϵ_{zz}ϵ_{xx})^2`$ (4)
Fig. 3(b) shows the diagrams of the conduction and valence band edges which are given by

$`\mathrm{\Delta }U_c=\mathrm{\Delta }E_c+a_cH`$ (5)
$`\mathrm{\Delta }U_v=\mathrm{\Delta }E_va_vHbB/\sqrt{2}`$ (6)
where $`\mathrm{\Delta }E_c`$ ($`\mathrm{\Delta }E_v`$) is the conduction (valence) band offset, $`a_c`$ ($`a_v`$) is the conduction (valence) band hydrostatic deformation potential parameter, $`b`$ is the valence band shear deformation potential parameter. Hence, while electrons are only sensitive to hydrostatic strain, holes are also and mostly affected by biaxial strain.

In the stacked structure, the hydrostatic negative strain field resides entirely inside the dots, and is a little stronger than in the single dot, while it almost vanishes in the coupling region (lower curve in Fig. 3a). Therefore, the coupling region is seen by electrons as a conventional tunneling barrier, which results in bonding and antibonding states as shown in Fig. 2. The biaxial positive strain field is however seen very differently in the stacked structure from that in the single dot (upper curve in Fig. 3a). First, it is smaller than in the single dot, leading to hole levels with higher energies in the stacked system (see Fig. 1). Second, unlike the hydrostatic strain, the biaxial strain retains a substantial value in the coupling region in the stacked system (top inset of Fig. 3a), which noticeably reduces the barrier height in the valence bands. Third, the biaxial strain profile in the stacked structure is inverted in the upper dot compared to the lower one although symmetric with respect to a median plane between the dots.

The most important point is that the biaxial strain defines two triangular confining potentials in each dot of the stacked structure, while in the single dot the valence band edge profile exhibits a much smoother slope (see Fig. 3b). It is seen that the strain-induced triangular potential is inverted with respect to the base in the two dots: In the bottom dot, the valence band edge is higher at the base than at the top, while it is opposite in the top dot. Therefore, the favorable combination of the lower strain-induced potential and the wider base result in a weaker hole confinement in the bottom dot than in the top dot, thereby localizing the hole ground state in the lower dot. The inset of Fig. 3(b) schematically shows the potential in the lower dot of the stacked structure together with the ground hole state at different electric fields. At zero electric field, the strain-induced built-in potential variation in the dot is $`34`$ meV and the ground hole state is localized approximately one fourth of the lower dot height from the bottom base.

This situation is illustrated in the middle panel of Fig. 4(a) where we plot the probability density isosurfaces of ground electron and hole states at different electric fields. In contrast, it is seen that the ground electron state is the bonding $`1s^+`$ state extending quasi-equally in both dots and in the coupling region. Under strong positive electric field, the ground electron state is seen in the right panel of Fig. 4(a) to become a $`1s`$-like state entirely localized in the upper dot. However, the ground hole state is found to behave very differently. Because the hole state is localized inside this triangular potential, as schematically shown in the inset of Fig. 3(b), its energy level changes approximately by the same amount as external variation in the dot, i.e., the potential drop arising from the external electric field . For example, the variation of ground hole state energy between $`F=200`$ kV/cm and zero field is about $`70`$ meV, which is about the same value as the potential drop over the bottom dot caused by electric field $`F=200`$ kV/cm. Therefore it is seen in Fig. 2 that the energy of the hole state has an almost linear dependence on the electric field. This also explains the quasi-linear dependence of the ground state transition energy as shown in Fig. 1.

In Fig. 3(b), the ratio between the tunneling barrier and band offset in the valence bands is much smaller than that in the conduction bands. Consequently, it is much easier for holes to tunnel through the barrier at negative electric fields than electrons. Fig. 4(b) illustrates the transition of the ground hole state from the bottom dot to the top one occuring at a small negative electric field $`F=21.6`$ kV/cm. As the biaxial strain field has different distribution in the respective dots, the hole states exhibit different probability density profiles when localizing in different dots. In Fig. 2, this transition is seen as an anticrossing between the ground hole state and the first excited state, which also accounts for the anomalous quantum confined Stark effect.

Although the anomalous quantum confined Stark effect invalidates Eq. (1) over the whole range of electric field, it is possible to write two separate equations that describe the dependence of the transition energy on electric field for the stacked structure, except for the electric fields at which anticrossings occur.

$`E^+(F)=E(0)+p^+F+\beta ^+F^2,`$ (7)
$`E^{}(F)=E^{}+p^{}F+\beta ^{}F^2.`$ (8)
where $`E^{}`$ is a fitting parameter, $`E^+(F)`$ is for $`F>21.4`$ kV/cm, and $`E^{}(F)`$ is for $`F<21.8`$ kV/cm.

For the right branch, $`p^+/|e|=3.65`$ nm is the built-in dipole moment of the stacked structure at zero field. From the middle panel of Fig. 4(a), the ground electron and hole states are seen to be separated by a distance larger than the coupling region ($`1.8`$ nm), and much larger than the electron-hole separation in the single dot, which results in a dipole moment of more than one order of magnitude larger than that in the single dot. For the left branch, $`p^{}p^+`$ nm has a positive sign because of the inverted electron-hole alignment  in the stacked structure at negative electric fields.

The quantum confined Stark effect coefficient is defined by $`\beta =\beta _e\beta _h`$ where $`\beta _e`$ ($`\beta _h`$) can be related to oscillator strength of optical intraband transitions in conduction (valence) bands . For instance, $`\beta _e`$ is given by

$$\beta _e=\frac{\mathrm{}^2}{2m_0}\underset{n>1}{}f_{1n}/(E_nE_1)^2,$$
(9)
where $`m_0`$ is the bare electron mass, $`f_{1n}`$ is the oscillator strength for the intraband transition from the ground state to the n-th state, with polarization along $`z`$ direction. $`\beta _h`$ has a similar expression. Single SAD structures have been shown to have very weak $`z`$-polarized intraband transitions , therefore both $`\beta _e`$ and $`\beta _h`$ are small. In the stacked structure, the strengths of the $`z`$-polarized intraband transitions in the valence bands are similar to those in the single dot. i.e., $`\beta _h`$ is small, which could be observed in Fig. 2 where the hole state energies show a quasi-linear dependence on the electric field. However, in the conduction bands of the stacked structure, there are several strong intraband transitions, especially the $`1s^+1s^{}`$ transition . This results in a much larger $`\beta _e`$ than in single dots, and explains the magnitude of the quantum confined Stark effect coefficients in both branches ($`\beta ^\pm /|e|^280\text{ nm}^2/eV`$) which are six times larger than in single dots, and are responsible for the ‘bowing’ of the Stark shift in Fig. 1.

In conclusion, we have shown that vertically stacked InAs/GaAs self-assembled quantum dots exhibits, in addition to a much larger Stark shift than single dots, a strong non-parabolic dependence of interband transition energy on the electric field, not been encountered in other types of quantum structures. We have demonstrated that the 3D distribution of the biaxial strain field is mainly responsible for this anomalous quantum confined Stark effect.

###### Acknowledgements.

This work is supported by Army Research Office, grant #DAAD 10-99-10129 and National Computational Science Alliance, grant #ECS000002N.


text:
# Calogero model, deformed oscillators and the collapse

## Abstract

We discuss the behaviour of the Calogero model and the related model of $`N`$ deformed oscillators with the $`S_N`$-extended Heisenberg algebra for a special value of the constant of interaction/statistical parameter $`\nu `$. The problem with finite number of deformed oscillators is analyzed in the algebraic approach, while collective-field theory is used to investigate the large-$`N`$ limit. In this limit, the system reduces to a large number of collapsing (free) particles, for $`\nu =1/N`$.

PACS number(s): 03.65.Fd, 03.65.Sq, 05.30.Pr

The (rational) Calogero model describes $`N`$ identical particles which interact through an inverse-square two-body interaction and are subjected to a common confining harmonic force. The inverse-square potential can be regarded as a pure statistical interaction and the model maps to an ideal gas of particles obeying fractional statistics. In one dimension, the usual notion of exchange statistics is not well defined because the exchange necessarily involves scattering, but instead one can apply Haldane’s definition of exclusion statistics.

The model is completely integrable in both the classical and quantum case, and that makes it an ideal ”lab” for exploring new concepts in physics, like the already mentioned exclusion statistics. The discovery of the integrable spin-chain with the long-range interaction, the so-called Haldane-Shastry model and its connection to the trigonometric version of the Calogero model, introduced the model to condensed matter physics. The Calogero model also appears when two-dimensional QCD is formulated as a random matrix model in the large-$`N`$ limit, restricted to a singlet subspace. More precisely, pure QCD in two dimensions is equivalent to the one-dimensional unitary-matrix model, which is a special case of the Calogero model.

Recently, the interest in the model has been renewed since it was proposed that the supersymmetric extension of the model could provide a microscopic description of the extremal Reissner-Nordström black hole. Many more applications of the model have been found, and intensified the research of the model intrinsic properties.

In this letter we investigate the behaviour of the Calogero model and the related model of N deformed oscillators with the $`S_N`$-extended Heisenberg algebra. We apply the Fock-space analysis and collective-field theory. We discuss the physical picture of the problem with the negative statistical parameter $`\nu `$. To our knowledge, this is the first analysis of the Calogero model and the related algebra of deformed oscillators for the negative, but still allowed $`\nu `$ leading to the positive definite Fock space, i.e. to quadratic-integrable wave functions, for general $`N`$. The two-particle case was previously discussed in Refs..

The Calogero model is defined by the following Hamiltonian:

$$H=\frac{\mathrm{}^2}{2m}\underset{i=1}{\overset{N}{}}_i^2+\frac{m\omega ^2}{2}\underset{i=1}{\overset{N}{}}x_i^2+\frac{\mathrm{}^2\nu (\nu 1)}{2m}\underset{ij}{\overset{N}{}}\frac{1}{(x_ix_j)^2},$$
(1)
where $`\nu (\nu 1)1/4`$ is the dimensionless coupling constant and $`\omega `$ is the strength of a harmonic confinement potential. In the following we set $`\mathrm{}=m=\omega =1`$, for simplicity. The ground-state wave function is, up to normalization,

$$\psi _0(x_1,\mathrm{},x_N)=\theta (x_1,\mathrm{},x_N)\mathrm{exp}\left(\frac{1}{2}\underset{i=1}{\overset{N}{}}x_i^2\right),$$
(2)
where $`\theta (x_1,\mathrm{},x_N)=_{i<j}^N(x_ix_j)^\nu `$, with the corresponding ground-state energy $`E_0=N[1+(N1)\nu ]/2`$. The physical requirement of vanishing of the wave function $`\psi _0(x_1,\mathrm{},x_N)`$ for coincidence points implies $`\nu >0`$. We note in passing that for $`0<\nu 1/2`$ one needs to renormalize the theory, see Refs..

One can introduce the following analogs of creation and annihilation operators:

$$a_i^{}=\frac{1}{\sqrt{2}}(D_i+x_i),a_i=\frac{1}{\sqrt{2}}(D_i+x_i),$$
(3)
where

$$D_i=_i+\nu \underset{j,ji}{\overset{N}{}}\frac{1}{x_ix_j}(1K_{ij})$$

are Dunkl derivatives. The elementary generators $`K_{ij}`$ of the symmetry group $`S_N`$ exchange labels $`i`$ and $`j`$:

$`K_{ij}x_j=x_iK_{ij},K_{ij}=K_{ji},(K_{ij})^2=1,`$ (4)
$`K_{ij}K_{jl}=K_{jl}K_{il}=K_{il}K_{ij},\mathrm{for}ij,il,jl,`$ (5)
and we choose $`K_{ij}|0=|0`$. One can easily check that the commutators of the creation and annihilation operators (3) are

$$[a_i,a_j]=[a_i^{},a_j^{}]=0,[a_i,a_j^{}]=\left(1+\nu \underset{k=1}{\overset{N}{}}K_{ik}\right)\delta _{ij}\nu K_{ij}A_{ij}.$$
(6)
The algebra (6) is consistent for all values of the parameter $`\nu `$. The only restriction on the parameter $`\nu `$ comes from the representation of the algebra on the Fock space. By performing a similarity transformation of the Hamiltonian (1), we obtain the reduced Hamiltonian

$$H^{}=\theta ^1H\theta =\frac{1}{2}\underset{i=1}{\overset{N}{}}\{a_i,a_i^{}\}=\underset{i=1}{\overset{N}{}}a_i^{}a_i+E_0,$$
(7)
when restricted to symmetric functions and acting as a total number operator. The physical space of $`N`$ identical bosons is a space of totally symmetric functions (for $`N`$ identical fermions one chooses a space of totally antisymmetric functions). We can then rewrite the Hamiltonian $`H^{}`$ as a Hamiltonian for deformed oscillators:

$$H^{}=\frac{1}{2}\underset{i=1}{\overset{N}{}}(p_i^2+x_i^2),$$
(8)
where $`p_i=iD_i=i(a_ia_i^{})/\sqrt{2}`$ is the hermitian momentum operator, $`x_i=(a_i+a_i^{})/\sqrt{2}`$ and the deformed commutation relation is

$$[x_i,p_j]=i\left(1+\nu \underset{k=1}{\overset{N}{}}K_{ik}\right)\delta _{ij}i\nu K_{ij}.$$

Let us now investigate the condition under which the states in the complete Fock space of algebra (6) have positive norm. To this end, we calculate the general matrix element of the Gram matrix for one-particle excited states:

$$0|a_ia_j^{}|0=0|A_{ij}|0=\delta _{ij}(1+N\nu )\nu ,$$
(9)
where we have used the fact that $`a_i|0=0`$ and $`0|0=1`$. There is only one non-degenerate eigenvalue $`\lambda _1=1/N`$ with the corresponding eigenvector $`e_1=_{i=1}^Na_i^{}|0`$, and $`N1`$ degenerate eigenvalues $`\lambda _{N1}=1+N\nu `$ with eigenvectors $`e_{N1}=(a_1^{}a_i^{})|0,i=2,3,\mathrm{},N`$. The Gram matrix is positive definite, i. e. there are no vectors with negative norm if all eigenvalues are positive. Thus

$$1+N\nu >0\nu >\frac{1}{N}.$$

At $`\nu =1/N`$ there is a critical point where the states $`(a_i^{}a_j^{})|0`$ have null normIf we had extracted the prefactor $`\theta (x_1,\mathrm{},x_N)=_{i<j}^N(x_ix_j)^{1\nu }`$ from the wave function, we would have obtained the critical point $`\nu =1+1/N`$.. It should be mentioned that the necessary condition for positivity is $`0|a_ia_i^{}|0>0`$, i. e. $`\nu >1/(N1)`$. However, this condition is not sufficient. Namely, for $`1/(N1)<\nu <1/N`$, the states $`(a_i^{}a_j^{})|0`$ have negative norm, although the mean value $`0|a_ia_i^{}|0`$ is positive. One can show that for two-, three- and many-particle states, the same condition for positivity of eigenvalues, i. e. $`\nu >1/N`$, is required. There exists the universal critical point at $`\nu =1/N`$, where all matrix elements of the arbitrary k-multi-state matrix are equal to $`k!/N^k`$. Namely, for one-particle states, $`0|a_ia_j^{}|0=1/N`$ at the critical point $`\nu =1/N`$. For two-particle excited states, we obtain

$$0|a_{i_2}a_{i_1}a_{j_1}^{}a_{j_2}^{}|0=0|a_{i_2}A_{i_1,j_1}a_{j_2}^{}|0+0|a_{i_2}a_{j_1}^{}A_{i_1,j_2}|0=\frac{2}{N^2},$$
(10)
since $`0|A_{i,j}|0=1/N`$ and $`0|a_{i_2}A_{i_1,j_1}a_{j_2}^{}|0=1/N^2`$ at the critical point. It can be shown by induction that the general matrix element at $`\nu =1/N`$ is given by

$$0|a_{i_k}\mathrm{}a_{i_1}a_{j_1}^{}\mathrm{}a_{j_k}^{}|0=\frac{k!}{N^k},k,N2.$$
(11)
If $`\nu >1/N`$, the diagonal matrix elements satisfy the inequality

$$0|a_{i_k}\mathrm{}a_{i_1}a_{i_1}^{}\mathrm{}a_{i_k}^{}|0>\frac{k!}{N^k},$$

which ensures the positivity condition.

What is the interpretation of these results like? For $`\nu >1/N`$, we see that all k-states have positive norm. The corresponding Gram matrix of order $`N^k`$ shows $`N`$ independent $`\nu `$-deformed oscillators which can be mapped to $`N`$ bosons and vice versa. However, at the critical point there is only one state for any $`k`$ with the eigenvalue $`k!>0`$. All other eigenvalues vanish! This means that in the limit $`\nu N1`$ the system of deformed oscillators exhibits singular behaviour. There survives only one oscillator $`_{i=1}^Na_i^{}/\sqrt{N}`$ describing the motion of the centre-of-mass coordinate. All other $`N1`$ relative coordinates collapse into the same point - the centre of mass. One can easily check that $`\stackrel{~}{a}_i^{}|0>=0`$ and therefore $`\stackrel{~}{p}_i|0>=\stackrel{~}{x}_i|0>=0`$, where ”$``$” denotes operators defined in the centre-of-mass system. For example, $`\stackrel{~}{a}_i^{}=a_i^{}_ja_j^{}/N`$. This means that the relative coordinates and the relative momenta are zero at the critical point. Also, the relative energy vanishes for this critical value of the parameter $`\nu `$. This collapse in the deformed quantum-mechanical system can be regarded as a toy model for confinement.

The interval $`(1/N,0)`$ for the statistical parameter $`\nu `$ is physically acceptable for the Hamiltonian $`H^{}`$ describing $`N`$ oscillators with the $`S_N`$-extended Heisenberg algebra, but it is not allowed for the original Calogero Hamiltonian. Hence, there is no collapse of Calogero particles. The wave function of the ground state (2) contains the Jastrow factor, i. e. the product of all relative coordinates with $`\nu `$ as a power. Consequently, the wave functions diverge for the negative value of $`\nu `$, at coincidence points, although for $`\nu >1/N`$ they remain normalizable. The Hamiltonian $`H^{}`$ obtained by similarity transformation (7) leads to the deformed oscillators with $`S_N`$-extended Heisenberg algebra and is equivalent to the Calogero model only for $`\nu >0`$.

Starting from the symmetric Fock space and the algebra of $`S_N`$-symmetric observables in the Calogero model, one can construct the mapping to free Bose oscillators (and vice versa) for $`\nu >1/N`$. Then it immediately follows that the partition function is

$$𝒵=e^{\beta E_0}\underset{k=1}{\overset{N}{}}\frac{1}{1e^{\beta k}}.$$

However, all observables disappear at the collapse point, leaving only the centre-of-mass coordinate, so the partition function reduces to the partition function of a single oscillator $`𝒵_{\mathrm{crit}}=\mathrm{exp}(\beta /2)/(1\mathrm{exp}(\beta ))`$.

In order to check the consistency and further elucidate our physical picture of the behaviour of the system at the critical point, we now turn to the case of infinite number of particles. Let us consider the large$`N`$ Calogero model in the Hamiltonian collective-field approach based on the $`1/N`$ expansion. It has been shown, that in the large$`N`$ limit, the Hamiltonian $`H^{}`$ can be expressed entirely in terms of density of particles $`\rho (x)`$ and its canonical conjugate $`\pi (x)=i\delta /\delta \rho (x)`$. The Jacobian $`J`$ of the transformation from $`x_i`$ into $`\rho (x)`$ rescales the symmetric wave functional

$$\varphi (x_1,x_2,\mathrm{},x_N)=\sqrt{J}\varphi (\rho )$$

resulting in the hermitian collective-field Hamiltonian

$$H=\frac{1}{2}𝑑x\rho (x)(_x\pi )^2+\frac{1}{8}𝑑x\rho (x)\left(_x\frac{\delta \mathrm{ln}J}{\delta \rho (x)}\right)^2\frac{1}{2}𝑑x\rho (x)x^2,$$
(12)
and the Jacobian determined from the hermicity condition

$$_x\left(\rho (x)_x\frac{\delta \mathrm{ln}J}{\delta \rho (x)}\right)=(\nu 1)_x^2\rho (x)+2\nu _x\left(\rho (x)𝑑y\frac{\rho (y)}{xy}\right).$$
(13)
The first term in the collective Hamiltonian (12), quadratic in the conjugate momentum $`\pi (x)`$, represents the kinetic energy of the problem. To find the ground-state energy of our system, we assume that the corresponding collective-field configuration is static and has a vanishing momentum $`\pi `$. Therefore, the leading part of the collective-field Hamiltonian in the $`1/N`$ expansion is given by the effective potential

$$V_{\mathrm{eff}}(\rho )=\frac{1}{8}𝑑x\rho (x)\left[(\nu 1)\frac{_x\rho (x)}{\rho (x)}+2\nu 𝑑y\frac{\rho (y)}{xy}\right]^2+\frac{1}{2}𝑑x\rho (x)x^2.$$
(14)
This can be rewritten as a complete square, up to the numerical constant (the ground-state energy):

$$V_{\mathrm{eff}}(\rho )=\frac{1}{8}𝑑x\rho (x)\left[(\nu 1)\frac{_x\rho (x)}{\rho (x)}+2\nu 𝑑y\frac{\rho (y)}{xy}2x\right]^2+\frac{N}{2}[\nu (N1)+1].$$
(15)
The first term in (15) is positive semidefinite, and, therefore, its contribution to the ground-state energy vanishes if there exists a positive solution of the first-order differential Bogomol’nyi-type equation:

$$(\nu 1)_x\rho (x)+2\nu \rho (x)𝑑y\frac{\rho (y)}{xy}2x\rho (x)=0.$$
(16)
What happens at the critical value $`\nu =1/N`$? Since the collective field obeys the normalization condition $`𝑑x\rho (x)=N`$, it is legitimate to rescale $`\rho (x)=N\stackrel{~}{\rho }(x)`$. Rewriting Eq.(16) for the rescaled collective field $`\stackrel{~}{\rho }(x)`$, we obtain at the critical point, in the leading approximation in $`1/N`$:

$$_x\stackrel{~}{\rho }(x)+2\stackrel{~}{\rho }(x)𝑑y\frac{\stackrel{~}{\rho }(y)}{xy}+2x\stackrel{~}{\rho }(x)=0.$$
(17)
Having in mind the identity for the singular distribution:

$$\frac{\delta (x)}{x}+2\delta (x)P\frac{1}{x}=0,$$
(18)
where P denotes the principal distribution, it is evident that the normalizable solution of Eq.(17) is given by the $`\delta `$-function, $`\stackrel{~}{\rho }(x)=\delta (x)`$. All particles are concentrated around the origin with density given by

$$\rho (x)=N\delta (x).$$
(19)
The corresponding ground-state energy reduces to $`E_0=1/2`$, i. e. to the ground state of the single harmonic oscillator.

Let us briefly summarize the main findings of the collective-field approach. We confirm the existence of the critical point $`\nu =1/N`$ found using the exact algebraic approach valid for any $`N`$. Also, we observe the collapse of the system accompanied by the simultaneous emerging of a single free oscillator.

We point out that for $`N\mathrm{}`$, $`\nu =1/N`$, the Fock space representation of the $`S_N`$-extended Heisenberg algebra (6) reduces to the Fock space of the centre-of-mass oscillator, although the $`S_N`$-extended Heisenberg algebra becomes the ordinary Heisenberg algebra for infinitely many oscillators. Hence, for $`\nu =1/N`$, the system reduces to a large number of collapsing (free) particles. On the other hand, the limit $`\nu 0_+`$ is regular, i.e. the Fock space corresponds to the ordinary bosons for the ordinary Heisenberg algebra.

The behaviuor of the system near the critical point is a universal feature of a large class of permutation-invariant algebras, which can be written as

$$a_ia_j^{}qa_j^{}a_i=\delta _{ij}[1+\underset{l=1}{\overset{N}{}}\nu _{il}K_{il}]\nu _{ij}K_{ij},|q|1.$$

When all $`\nu _{ij}`$ tend to the critical point $`1/N`$, only a single q-deformed oscillator, representing the centre of mass, survives. It gives us a new and interesting form of matter. We hope that these ideas might be applied to problems like black-hole physics, extremely dense matter, anyons and confinement.

Acknowledgment

We thank the referee for critical consideration of the manuscript. This work was supported by the Ministry of Science and Technology of the Republic of Croatia under contracts No. 00980103 and 0119222.