text:
# Inflation in string-inspired cosmology and suppression of CMB low multipoles

## Abstract

We construct a string-inspired nonsingular cosmological scenario in which an inflaton field is driven up the potential before the graceful exit by employing a low-energy string effective action with an orbifold compactification. This sets up an initial condition for the inflaton to lead to a sufficient amount of $`e`$-foldings and to generate a nearly scale-invariant primordial density perturbation during slow-roll inflation. Our scenario provides an interesting possibility to explain a suppressed power spectrum at low multipoles due to the presence of the modulus-driven phase prior to slow-roll inflation and thus can leave a strong imprint of extra dimensions in observed CMB anisotropies.

String theory has been regarded as one of the promising candidate theories of quantum gravity over the past 20 years string . However it is unfortunately not easy to detect the signature of extra dimensions and any stringy effect in accelerators in the foreseeable future. In this sense string cosmology presumably provides the only means of testing string or M-theory concretely. It is therefore very important to investigate the viability of string theory by exploring possible cosmological implications.

From the observational side, the recent WMAP measurements showed a lack of power on the largest scales and the running of the spectral index in temperature anisotropies Spergel . This is difficult to be explained by standard slow-roll inflation, thus may have roots in esoteric Planck scale physics. One can consider an exciting possibility where the suppression of power originates from some stringy effect which may be present prior to standard inflation (see e.g., Refs. sup  for a number of attempts to explain this suppression).

String theory provides an interesting possibility to avoid the big-bang singularity due to its underlying symmetries BV . The construction of nonsingular cosmological models based on the low-energy string effective action was pioneered by Veneziano and Gasperini Veneziano:91 , whose model is called the Pre-Big Bang (PBB) scenario. According to this scenario, the universe exhibits a super-inflation when the dilaton field, which characterizes the strength of the string coupling parameter, evolves toward the strong coupling regime. The curvature singularity can be avoided by implementing the dilatonic loop and higher-derivative corrections to the tree-level action, which is followed by a post-big bang phase with a decreasing Hubble rate Bru .

While the PBB scenario is an intriguing attempt to merge string theory and cosmology, the spectrum of curvature perturbations is highly blue-tilted ($`n4`$) in its simplest form BGGMV , which is incompatible with observations Spergel . If we take into account the contribution of an axion field, the spectrum of the axion perturbation can be scale-invariant depending on the expansion rate of an internal dimension Copeland:97 . Although there remains a possibility to explain the observationally supported flat spectra if the axion plays the role of the curvaton curvaton , the PBB scenario still suffers from the problem of the requirement of an exponentially large homogenous region in order to solve the major cosmological problems Kaloper .

If standard slow-roll inflation occurs after a graceful exit, this solves the homogeneity problem above as long as we have sufficient amount of $`e`$-foldings ($`N60`$). In this paper we shall investigate a scenario in which the PBB phase is followed by slow-roll inflation with an assumption that there exists a light scalar field (inflaton) other than the dilaton*** This field can be regarded as one of the modulus fields which is always present in string theory.. Note that a similar idea was proposed in Ref. Piao  where slow-roll inflation is preceded by a contracting phase. In this work we implement the tree-level $`\alpha ^{}`$ Gauss-Bonnet correction to the string effective action in order to construct viable nonsingular solutions. We will show that the inflaton can be driven up the potential hill due to a nontrivial background dynamics prior to the graceful exit, thereby setting up sufficient initial conditions for following slow-roll inflation. This has a similarity to the inflation in loop quantum gravity quantum , but the observational signatures are different. In fact we shall see that the presence of the kinematically driven super-inflation stage before the graceful exit can lead to a stronger suppression of the power spectrum on largest scales, as favoured from observations.

We start with the following 4-dimensional action that appears for the orbifold compactification in low-energy effective string theory Antoniadis

$`S={\displaystyle d^4x\sqrt{g}\left[\frac{R}{2}\frac{1}{4}(\varphi )^2\frac{3}{4}(\sigma )^2+_{\mathrm{inf}}+_c\right]},`$ (1)
which is written in the Einstein frame. Here $`R`$ is the scalar curvature and the fields $`\varphi `$ and $`\sigma `$ correspond to the dilaton and the modulus, respectively. We shall consider a minimally coupled inflaton field, $`\chi `$, with lagrangian

$`_{\mathrm{inf}}=(1/2)(\chi )^2V(\chi ).`$ (2)
Hereafter we assume that the inflaton has a mass $`m_\chi `$ with potential $`V(\chi )=(1/2)m_\chi ^2\chi ^2`$.

The correction term $`_c`$ is given by

$`_c=(1/16)\left[\lambda e^\varphi \delta \xi (\sigma )\right]R_{\mathrm{GB}}^2,`$ (3)
where $`R_{\mathrm{GB}}^2=R^24R^{\mu \nu }R_{\mu \nu }+R^{\mu \nu \alpha \beta }R_{\mu \nu \alpha \beta }`$ is a Gauss-Bonnet term. Note that we are considering the correction up to the first order in $`\alpha ^{}`$. The coefficients, $`\lambda `$ and $`\delta `$, are determined by the inverse string tension $`\alpha ^{}`$ and the 4-dimensional trace anomaly of the $`N=2`$ sector. $`\lambda `$ is a positive constant, whereas $`\delta `$ can be positive or negative. The function, $`\xi (\sigma )`$, is written in terms of the Dedekind $`\eta `$-function and is well approximated as $`\xi (\sigma )(2\pi /3)\mathrm{cos}\mathrm{h}\sigma `$. The dilatonic higher-order loop and derivative corrections may be included as in Ref. Bru , but we do not take them into account in this work. In fact singularity-free cosmological solutions can be constructed if the modulus-dependent term dominates over the dilatonic one in Eq. (3).

With the flat Friedmann-Robertson-Walker (FRW) background described by the scale factor $`a`$, the variation of the action (1) yields the following background equations

$`12H^2=\dot{\varphi }^2+2\dot{\chi }^2+3\dot{\sigma }^2+2m_\chi ^2\chi ^296H^3\dot{f},`$ (4)
$`8(1+8H\dot{f})(\dot{H}+H^2)+4(1+8\ddot{f})H^2`$
$`+\dot{\varphi }^2+2\dot{\chi }^2+3\dot{\sigma }^22m_\chi ^2\chi ^2=0,`$ (5)
$`\ddot{\varphi }+3H\dot{\varphi }+2f_\varphi R_{\mathrm{GB}}^2=0,`$ (6)
$`\ddot{\sigma }+3H\dot{\sigma }(2/3)f_\sigma R_{\mathrm{GB}}^2=0,`$ (7)
$`\ddot{\chi }+3H\dot{\chi }+m_\chi ^2\chi =0,`$ (8)
where a dot denotes the derivative in terms of cosmic time $`t`$, $`H\dot{a}/a`$ is the Hubble parameter and $`f=(1/16)\left[\lambda e^\varphi \delta \xi (\sigma )\right]`$. The Gauss-Bonnet term is simply given as $`R_{\mathrm{GB}}^2=24H^2(\dot{H}+H^2)`$.

The universe starts out from a weak coupling regime with $`e^\varphi 1`$. In the absence of the modulus, the dynamics of the system in the weak coupling region is characterized by the dilaton-driven phase in PBB cosmology Veneziano:91 . During this phase the universe contracts as $`a(t)^{1/3}`$ in the Einstein frame. This solution does not connect to our expanding branch in a nonsingular way in the tree-level action. We can construct singularity-free bouncing solutions if higher-order loop and derivative corrections are taken into account Bru . Nevertheless it is required to fine-tune the coefficients of the correction terms. In addition the perturbations inside the Hubble radius exhibit a negative instability Cartier , which casts doubt for the validity of linear perturbation theory.

The situation is different when the modulus field is present. If the modulus dynamically controls the system rather than the dilaton in Eqs. (4)-(8), one gets the super-inflationary solution with a growing Hubble rate in the Einstein frame. In particular we can construct nonsingular solutions by taking into account the correction (3) for negative $`\delta `$. For positive $`\delta `$, the first-order $`\alpha ^{}`$ corrections do not help to lead to a successful graceful exit, as analyzed in Ref. Toporensky .

When the modulus-type $`\alpha ^{}`$ corrections dominate the system, we get the background evolution for $`t<0`$ as Kawai

$`aa_0,H=H_0/t^2,\dot{\sigma }=5/t,`$ (9)
where $`a_0`$ and $`H_0`$ ($`>0`$) are constants. When one starts in an expanding branch, the background is characterized by the super-inflationary solution (9) until the graceful exit ($`t0`$). The Universe expands slowly initially with a nearly constant scale factor. Eventually the Hubble rate reaches its peak value $`H=H_{\mathrm{max}}`$, after which the system connects to a Friedmann-like branch. In the numerical simulation of Fig. 1, $`H_{\mathrm{max}}`$ corresponds to of order 0.1. To ensure the validity of the theory, this maximum value should be less than of order unity (i.e., less than the Planck energy). In fact there is a negative instability for the perturbations if $`H`$ is larger than unity Kawai .

We include the dilaton in our simulations by generalizing the analysis of Refs. Antoniadis ; Toporensky ; Kawai . As we find in Fig. 1, it is possible to obtain singularity-free solutions even when the dilaton is present as long as the modulus $`\alpha ^{}`$ correction dominates over the dilatonic one. If the super-inflationary modulus-driven phase leads to a large number of $`e`$-foldings and also generates a nearly scale-invariant density perturbation, this can be an alternative to standard slow-roll inflation. Unfortunately the scale factor evolves very slowly during this modulus-driven phase (see Fig. 1). In addition the spectrum of the density perturbation generated in modulus-driven inflation is highly blue-tilted, as we see later.

We are interested in whether the presence of the modulus-driven phase can set the sufficient initial conditions for inflaton to drive slow-roll inflation after the graceful exit. We shall consider a natural situation where the inflaton is initially located around the potential minimum, $`\chi 0`$. If the contribution of the potential term is dropped in Eq. (8), one has $`\dot{\chi }a^3\mathrm{const}`$. Therefore $`\chi `$ grows linearly during the super-inflationary phase, whose growth continues up to the Hubble peak \[see Eq. (9)\]. In Fig. 1 we show one example of the evolution for $`\chi `$. The inflaton is driven up its potential hill during the modulus-driven stage even in the presence of the inflaton potential. The maximum value of $`\chi `$ ($`=\chi _{\mathrm{max}}`$) reached at the Hubble peak is $`\chi _{\mathrm{max}}=2.9`$ in a Planck unit in this case. This maximum value increases if the initial values of $`\dot{\chi }`$ get larger.

We wish to analyze the case in which slow-roll inflation occurs after the graceful exit. If the large-scale power spectrum observed in CMB is generated during this second stage of inflation, the mass $`m_\chi `$ is constrained to be $`m_\chi 10^6m_{\mathrm{Pl}}`$ from the COBE normalization ($`m_{\mathrm{Pl}}`$ is the Planck mass). Since this value is much smaller than $`m_{\mathrm{Pl}}`$, the energy density of the inflaton is practically negligible compared to that of the modulus around the graceful exit when $`H_{\mathrm{max}}`$ is not too much smaller than unity. Therefore the presence of the inflaton potential is not harmful for the construction of nonsingular cosmological solutions.

The radiation is expected to be created by the decays of the modulus and the dilaton after the transition to the Friedmann branch. It was shown in Ref. Barreiro:1998aj  that the dilaton can be stabilized by the potential generated by a gaugino condensation in the presence of the radiation. We shall assume that the dilaton and the modulus fall down to their potential minimums soon after the graceful exit through a similar mechanism.

The inflaton freezes with a value $`\chi \chi _{\mathrm{max}}`$ due to a large friction term until the Hubble rate drops down to of order $`m_\chi `$. Since the energy density of the radiation decreases faster than that of the field $`\chi `$, the background dynamics is eventually dominated by the inflaton. The slow-roll inflation takes place when the inflaton starts to roll down along the potential for $`Hm_\chi `$. The number of $`e`$-foldings is dependent on the field value $`\chi _{\mathrm{max}}`$. Since we numerically found that $`\chi _{\mathrm{max}}`$ can exceed 3 in a Planck unit, it is possible to have more than 60 $`e`$-foldings required to solve horizon and flatness problems. We are interested in the case where the information of the mudulus-driven phase before the transition can be imprinted on the largest scales in CMB, corresponding to the value of $`\chi _{\mathrm{max}}`$ that is slightly less than $`3m_{\mathrm{Pl}}`$.

Let us consider scalar perturbations about the FRW background with a single-field, $`\psi `$. Note that our system involves multiple scalar fields, but it is sufficient to study the perturbations of one dynamically important field in estimating the spectrum of metric perturbations. Then it is convenient to introduce a gauge-invariant variable,

$`\phi (H/\dot{\psi })\delta \psi ,`$ (10)
where $`\phi `$ is the scalar perturbation in the perturbed metric Hwang:99 . In the uniform-field gauge each Fourier mode of the perturbed Einstein equation with the correction (3) is written as Kawai

$`\ddot{}_k+(2/t)\dot{}_k\beta (k^2/a_0^2)t_k=0,`$ (11)
where $`k`$ is a comoving wavenumber and $`\beta `$ and $`a_0`$ are constants. The solution of this equation is written in terms of the Bessel functions, $`J_{\pm 1/3}`$, and the spectrum of the large-scale curvature perturbation is given by

$`𝒫_{}k^3|_k|^2/2\pi ^2k^{n_11},\mathrm{with}n_1=10/3,`$ (12)
which is a highly blue-tilted spectrum Kawai . Note that when the kinematic term dominates over the potential term the blue-tilted spectra are generic features even in generalized Einstein theories Hwang:1996hz ; Tsujikawa:01 . Around the graceful exit, $`\beta `$ can be negative when the Hubble rate grows larger than the Planck energy Kawai , which means that perturbations on the scales inside the Hubble radius exhibit a negative instability. However, we can avoid this instability as long as $`H_{\mathrm{max}}`$ is smaller than of order the Planck energy. This is contrasted with the dilaton-driven inflation in which inclusion of the $`\alpha ^{}`$ correction is accompanied by a similar negative instability even for $`H_{\mathrm{max}}m_{\mathrm{Pl}}`$ Cartier . Therefore perturbations in the modulus-driven case which crossed the Hubble radius around the graceful exit is more stable than in the dilaton-driven case.

The spectral index of curvature perturbations generated in slow-roll inflation after the graceful exit is close to scale-invariant and is given as $`n_21(1/\pi )(m_{\mathrm{Pl}}/\chi )^2`$ quantum . When $`\chi _{\mathrm{max}}3m_{\mathrm{Pl}}`$, corresponding to the $`e`$-folds $`N60`$, the CMB spectrum on cosmologically relevant scales is determined by the spectral index $`n_2`$. Meanwhile, when $`\chi _{\mathrm{max}}3m_{\mathrm{pl}}`$, the spectrum (12) is within an observational range on the scales larger than a cut-off which connects two power spectra. Then we simply model the power spectrum to be $`𝒫_{}=A_1(k/k_0)^{n_11}`$ for $`kk_c`$ and $`𝒫_{}=A_2(k/k_0)^{n_21}`$ for $`kk_c`$, where $`k_c`$ is a cut-off wavenumber. As shown in Fig. 2 it is possible to explain the loss of power at low multipoles by employing the highly blue-tilted spectrum (12), as long as $`\chi _{\mathrm{max}}`$ is slightly less than $`3m_{\mathrm{Pl}}`$. We have carried out the likelihood analysis by varying 7 cosmological & inflationary parameters using the recent WMAP Spergel  and 2dF 2dF  data sets, and found that the best-fit value of the cut-off scale corresponds to $`k_c=3.5\times 10^4\mathrm{Mpc}^1`$ (see the inset of Fig. 2). The smaller cut-off corresponding to $`\chi _{\mathrm{max}}3m_{\mathrm{Pl}}`$ is within the $`2\sigma `$ contour bound, thus not ruled out observationally. We also analyzed the smoothly connected power spectrum, $`P_{}=A(k/k_0)^{n_21}\left[1\mathrm{exp}\{(k/k_c)^{n_1n_2}\}\right]`$, around the cut-off $`k_c`$, which may be a possible case when higher-order loop corrections become important around the graceful exit. We obtained a similar likelihood value of the cut-off as shown in Fig. 2, which is also consistent with the past related works in Refs. sup ; Bridle .

In summary, we studied a nonsingular cosmological scenario based on the orbifold compactification in low-energy effective string theory and showed that the inflaton is driven up the potential hill during the modulus-driven stage before the graceful exit, thus setting up sufficient initial conditions for following slow-roll inflation. The presence of the modulus-driven phase generates a blue-tilted primordial power spectrum, which can explain the loss of power observed in CMB low multipoles. The spectral index (12) is much larger than in the case of noncommutative inflation noncom  or loop quantum gravity quantum , thus leading to a stronger suppression. Our work provides an exciting possibility to pick up the signature of extra dimensions from high-precision observations.

Acknowledgements We thank Bo Feng, Burin Gumjudpai, Roy Maartens, Parampreet Singh for useful discussions and Antony Lewis, David Parkinson for kind support in interpreting the likelihood analysis. S.T. is thankful for financial support from JSPS (No. 04942). The work of Y.P. and X.Z. is supported by the NSF of China and also by the Ministry of Sceince and Technology of China under grant No. NKBRSF G19990754.


format: text

text:
# A Hot Spot Model for Black Hole QPOs

## 1 Introduction

One of the most exciting results from the Rossi X-Ray Timing Explorer was the discovery of high frequency quasi-periodic oscillations from neutron star and black hole binaries Strohmayer et al. (1996); Strohmayer (2001); Lamb (2002). For black hole systems, these HFQPOs are observed repeatedly at more or less constant frequencies, and in a few cases with integer ratios Remillard et al. (2002). These measurements give the exciting prospect of determining a black hole’s mass and spin, as well as testing general relativity in the strong-field regime. To quantitatively interpret these observations in the context of the space-time behavior near the black hole inner-most stable circular orbit (ISCO), we have developed a ray-tracing code to model the X-ray light curve from a collection of “hot spots,” small regions of excess emission moving on geodesic orbits. The methods and basic results of this code are described in detail in Schnittman & Bertschinger (2003).

This hot spot model is motivated by the similarity between the QPO frequencies and the black hole coordinate frequencies near the ISCO Stella & Vietri (1998, 1999) as well as the suggestion of a resonance leading to integer commensurabilities between these coordinate frequencies Kluzniak & Abramowicz (2001, 2002). Stella and Vietri Stella & Vietri (1999) investigated primarily the QPO frequency pairs found in LMXBs with a neutron star accretor, but their basic methods can be applied to black hole systems as well. Perhaps the most powerful feature of this hot spot model is the facility with which it can be developed and extended to more general accretion disk geometries. For example, our ray-tracing code could be used in conjunction with a 3D MHD calculation of the accretion disk to simulate the time-dependent X-ray flux and spectrum from such a disk.

## 2 Ray Tracing in the Kerr Metric

While the details of the ray tracing algorithm are presented in full in Schnittman & Bertschinger (2003), we give here a short summary of our methods. We begin by dividing the image plane into regularly spaced “pixels” of equal solid angle in the observer’s frame, each corresponding to a single ray. Following the sample rays backward in time, we calculate the original position and direction that a photon emitted from the disk would require in order to arrive at the appropriate position in the detector. The gravitational lensing and magnification of emission from the plane of the accretion disk is performed automatically by the geodesic integration of these evenly spaced photon trajectories, so that high magnification occurs in regions where nearby points in the disk are projected to points with large separation in the image plane. To model the time-varying emission from the disk, each photon path is marked with the time delay along the path from the observer to the emission point in the disk.

To integrate the geodesic trajectories of photons or massive particles, we use a Hamiltonian formalism that takes advantage of certain conserved quantities in the dynamics: the energy, angular momentum, and mass of each particle. By eliminating the energy $`(E=p_t)`$ from the Hamiltonian, we can use the conjugate coordinate $`t`$ as the integration variable, thus reducing the problem from eight phase space dimensions to six. Another integral of the motion, Carter’s constant $`Q`$ Carter (1968), is not explicitly conserved in our equations, but rather is used as an independent check for the accuracy of the numerics. The resulting equations of motion do not contain any sign ambiguities from turning points in the orbits, as are introduced by many classical treatments of the geodesic equations in the Kerr metric.

The photon trajectories are integrated backward in time from the image plane oriented at some inclination angle $`i`$ with respect to the axis of rotation for the black hole, where $`i=0^{}`$ corresponds to a face-on view of the disk and $`i=90^{}`$ is an edge-on view. The accretion disk is confined to a finite region of latitude and is assumed to be oriented normal to the rotation axis. The photons terminate either at the event horizon or pass through the surfaces of colatitude $`\theta =\mathrm{const}`$. As trajectories pass through the disk, the photon’s position and momentum $`(x^\mu ,p_\mu )`$ are recorded for each plane intersection in order to later reconstruct an image of the disk. The results of this paper are based primarily on flat disks with simple hot spot perturbations, yet with the computational methods described above, arbitrary disk geometries and emissivity/opacity models can be simulated as well.

## 3 Hot Spot Emission

Given the map from the accretion disk to the image plane, with each photon bundle labeled with a distinct 4-momentum and time delay, we can reconstruct time-varying images of the disk based on time-dependent emission models. The simplest model we consider is a single region of isotropic, monochromatic emission following a geodesic trajectory: the “hot spot” model Stella & Vietri (1998, 1999). We treat the hot spot as a small region of the disk with additional emissivity chosen to have a Gaussian distribution in local Cartesian space. We typically take $`R_{\mathrm{spot}}=0.250.5M`$, but find the normalized light curves and QPO power spectra to be independent of spot size and shape.

By definition the hot spot will have a higher temperature or density and thus greater emissivity than the background disk, adding a small modulation to the total flux. RXTE observations find the HFQPO X-ray modulations to have typical amplitudes of 1-5% of the mean flux during the outburst Remillard et al. (2002). Assuming a Shakura-Sunyaev type disk Shakura & Sunyaev (1973) with steady-state emissivity $`g(r)r^2`$ and a similar scaling for the hot spot emission, we find that small hot spots with overbrightness of $`100\%`$ (consistent with 3D MHD calculations Hawley & Krolik (2001)) are easily capable of creating X-ray modulations on the order of 1% rms.

As the disk inclination increases, the light curve goes from nearly sinusoidal to being sharply peaked by special relativistic beaming. Thus the shape of a QPO light curve may be used to determine the disk inclination. Since current observational cannot resolve the X-ray signal over individual periods as short as 3-5 ms, instead the Fourier power spectrum can be used to identify the harmonic features of the light curve over many orbits. Disks with higher inclinations will give more power in the higher harmonic frequencies, due to the “lighthouse” effect, as the hot spot shoots a high-power beam of photons toward the observer once per orbit, approximating a periodic delta-function in time.

Figure 1 shows the quantitative dependence of harmonic power on disk inclination for a hot spot orbiting at the ISCO of a Schwarzschild black hole. Predictably, as the inclination increases, we see that both the absolute and relative amplitudes of the higher harmonics increase, almost to the limit of a periodic delta-function when $`i90^{}`$. Again, we find very little dependence of the harmonic structure on hot spot size or shape, as long as the total emission of the spot relative to the disk is constant. Here we have normalized the rms amplitudes to the background flux from the disk with a hot spot overbrightness of 100%. For a signal $`I(t)`$ with Fourier components $`a_n`$:

$$I(t)=\underset{n=0}{\overset{\mathrm{}}{}}a_n\mathrm{cos}(2\pi nt),$$
(1)
we define the rms amplitude $`a_n(\mathrm{rms})`$ in each mode $`n>0`$ as

$$a_n(\mathrm{rms})\frac{a_n}{\sqrt{2a_0}}.$$
(2)
With this normalization, the total rms can be conveniently written

$$\mathrm{rms}=\sqrt{a_n^2(\mathrm{rms})}.$$
(3)
## 4 Fitting QPO Data

The high frequency QPOs observed with RXTE can be interpreted in the context of the hot spot model by matching the model parameters to the features of the QPO power spectrum, particularly the locations, amplitudes, and widths of the peaks. This problem of fitting a power spectrum to the parameters of a multi-dimensional model is analogous to the spectacular successes that have recently been made in observing and interpreting the perturbations in the cosmological microwave background (CMB). Like the CMB power spectrum, the features of the QPO power spectrum often can be isolated and fitted to one or two of the model parameters at a time. We approach this problem in the systematic step-by-step procedure outlined below, using the example of the Type A QPOs from XTE J1550-564 Remillard et al. (2002) as case study.

First, we match the frequencies of the QPO peaks to the coordinate frequencies of specific geodesic orbits around a black hole in order to determine the black hole mass and angular momentum. The 2:3 frequency commensurability of the two major peaks (and possibly 2:3:4 in some observations) directs our attention to orbits with radial and azimuthal frequencies with a ratio of $`\nu _r:\nu _\varphi =`$1:3. The choice of the 1:3 ratio instead of 2:3 gives the appropriate relative power in the two peaks-a major peak at the fundamental azimuthal orbital frequency with additional power in the sidebands at $`\nu =\nu _\varphi \pm \nu _r`$. Conversely, a ratio of 2:3 in the coordinate frequencies would actually give a 1:3:5 ratio in the major power spectrum peaks. For any value of the black hole spin $`0a/M1`$, there exists a special radius where nearly circular orbits will have the appropriate 1:3 ratio in coordinate frequencies. Then by scaling the black hole mass, we can match the actual peak locations to the data, giving a one-dimensional degeneracy in mass-spin parameter space.

This degeneracy could be broken by identifying the low frequency QPOs with the Lens-Thirring precessional frequency at the same special radius. However, this association is still quite speculative and has difficulty explaining the large rms amplitude of the LFQPOs. For the present, we will constrain our model to a region of this degeneracy limited by the independent black hole mass determination of radial velocity measurements, giving $`8.5<M/M_{}<11.5`$ Orosz et al. (2002). For the QPO peaks at 184 Hz and 276 Hz in XTE J1550, we choose $`M=10.3M_{}`$ and $`a/M=0.5`$, giving a geodesic orbital radius of $`r_0=4.89M`$, somewhat outside the ISCO at $`R_{\mathrm{ISCO}}=4.23M`$.

Next, we match the total rms amplitudes in the two major peaks by selecting an average hot spot overbrightness of 100% and radius $`R_{\mathrm{spot}}=0.35M`$. The disk inclination is independently determined to be $`i70^{}`$ Orosz et al. (2002). As explained in the previous section, this relatively high inclination should produce significant power in many of the higher harmonics at integer multiples of the fundamental frequency $`\nu =n\nu _\varphi `$. These harmonics lie well within the sensitivity range of RXTE, yet are not detected with such high amplitudes. Instead, we find that by shearing the hot spot into an arc of finite length in azimuth, the power at higher frequencies is strongly damped, due to the countering of the lighthouse effect when the emission region is spread out over a larger portion of the orbit Schnittman & Bertschinger (2003). An arc length of $`180^{}`$ gives the appropriate amount of power around 184 Hz and 276 Hz while also successfully decreasing the power at higher frequencies.

The relative power in the two major peaks is primarily determined by the eccentricity of the geodesic orbit, as well as the arc length of the hot spot emission region. For the Type A QPOs, even the modest eccentricity of 0.11 is sufficient to produce significant power in the “side-band” at 184 Hz. Even the Type B QPOs, where the peak at 184 Hz is stronger than the fundamental peak at 276 Hz, can be explained easily in the context of the hot spot model by shearing the arc almost to a ring-shaped perturbation with arc length $`330^{}`$. For both cases, it is not obvious why there should be such a small amount of predicted power at $`\nu _\varphi +\nu _r=`$368 Hz, further advocating the necessity of a full ray-tracing algorithm, as opposed to simple analysis of the coordinate frequencies.

Up to this point, all of the ingredients of the hot spot model combine to produce a power spectrum with a series of discrete peaks describing a purely periodic light curve. In order to understand the significant broadening of the peaks (and thus the term quasi-periodic), we generalize the hot spot model by replacing a single spot orbiting indefinitely on a single trajectory by a collection of hot spots that are continually created and destroyed with random phase along their orbits. Each hot spot lifetime $`t`$ has an exponential distribution function $`f(t)\mathrm{exp}(t/T)`$ with a probability $`P=dt/T`$ of being “destroyed” during each time step $`dt`$. When one hot spot is destroyed, a new hot spot is created with the same $`r_0`$ but with random orbital phase. The effect is to spread each delta-function feature in the periodic power spectrum into a Lorentzian peak with FWHM $`1/T`$. With this multiple hot spot model, the power spectrum from XTE J1550 is closely matched by setting the typical lifetime $`T=15`$ ms, or approximately four orbits in azimuth.

The summary of these best-fit results is shown in Table 1, and the corresponding power spectrum (calculated from a Monte-Carlo simulation of the light curve over many hot spot lifetimes) is shown in Figure 2. The solid curve is the best fit of a collection of Lorentzians to the simulated data, included a small background counting noise level with a Poisson distribution. Again, the analogy to the CMB data is evident, and just like the with the CMB, each successive advance in instrumentation brings us closer to a precise determination of the model parameters. However, for the QPOs, in addition to the power spectrum, there is also significant information in the light curve phase, if we can achieve the necessary sensitivity to resolve the signal in the time domain. This phase information promises to ultimately help in distinguishing between the various proposed explanations of the cause of the black hole QPOs.

We have also investigated the effect of various distributions in orbital eccentricity, overbrightness, and the number of simultaneous hot spots at any time. Interestingly, we find that the QPO power spectrum is dependent only on the average eccentricity and overbrightness, and is not sensitive to the actual shape of the distribution. Similarly, the power spectrum is independent of the number of different hot spots in existence at any given time, sensitive only to the average lifetime of each hot spot.

One major remaining issue with the hot spot model is the preferred location of the geodesic that gives rise to 1:3 coordinate frequencies. Why should the orbital frequencies favor integer ratios, and why should the preferred ratio be 1:3 and not 1:2 or 1:4? It is possible that detailed radiation-hydrodynamic calculations with full general relativity will be required to answer this question. Perhaps the non-circular orbits can only survive along closed orbits such as these to somehow avoid destructive intersections. For now, we are forced to leave this as an open question unanswered by the geodesic hot spot model.

Yet with the computational framework developed above, these questions can be answered by modeling a whole collection of hot spots and arcs continually forming and evolving in shape and emissivity. The particular physical parameters for these hot spots can be derived from published MHD calculations such as Hawley & Krolik (2001). With the basic ray-tracing and radiation transport methods in place, it should be possible to use our code as a “post-processor” to analyze the results of these 3D simulations to simulate X-ray light curves and spectra from a realistic accretion disk.

A final piece of the black hole QPO puzzle is the spectral behavior of the source during outburst. Remillard et al. (2002) have shown there is a large range of X-ray fluxes with different relative contributions from a power-law component and a disk bolometric component. These relative and absolute fluxes seem to be correlated with the amount of power in both the LFQPOs and the HFQPOs. Future work on 3-dimensional disks and more detailed radiation transfer models should give us important insights into understanding this spectral behavior and its relation to the QPO power.

We would like to thank Ron Remillard for many helpful discussions. This work was supported by NASA grant NAG5-13306.


text:
# References

The Collatz 3n+1 Conjecture is Unprovable

Craig Alan Feinstein

2712 Willow Glen Drive, Baltimore, Maryland 21209
E-mail: cafeinst@msn.com, BS”D

Abstract: In this paper, we show that any proof of the Collatz $`3n+1`$ Conjecture must have an infinite number of lines; therefore, no formal proof is possible.

Disclaimer: This article was authored by Craig Alan Feinstein in his private capacity. No official support or endorsement by the U.S. Government is intended or should be inferred.

In 2005, the famous mathematician Freeman Dyson was asked, “What do you believe is true even though you cannot prove it?” He answered:

“Since I am a mathematician, I give a precise answer to this question. Thanks to Kurt Gödel, we know that there are true mathematical statements that cannot be proved. But I want a little more than this. I want a statement that is true, unprovable, and simple enough to be understood by people who are not mathematicians. Here it is.

“Numbers that are exact powers of two are 2, 4, 8, 16, 32, 64, 128 and so on. Numbers that are exact powers of five are 5, 25, 125, 625 and so on. Given any number such as 131072 (which happens to be a power of two), the reverse of it is 270131, with the same digits taken in the opposite order. Now my statement is: it never happens that the reverse of a power of two is a power of five.

“The digits in a big power of two seem to occur in a random way without any regular pattern. If it ever happened that the reverse of a power of two was a power of five, this would be an unlikely accident, and the chance of it happening grows rapidly smaller as the numbers grow bigger. If we assume that the digits occur at random, then the chance of the accident happening for any power of two greater than a billion is less than one in a billion. It is easy to check that it does not happen for powers of two smaller than a billion. So the chance that it ever happens at all is less than one in a billion. That is why I believe the statement is true.

“But the assumption that digits in a big power of two occur at random also implies that the statement is unprovable. Any proof of the statement would have to be based on some non-random property of the digits. The assumption of randomness means that the statement is true just because the odds are in its favor. It cannot be proved because there is no deep mathematical reason why it has to be true. (Note for experts: this argument does not work if we use powers of three instead of powers of five. In that case the statement is easy to prove because the reverse of a number divisible by three is also divisible by three. Divisibility by three happens to be a non-random property of the digits).

“It is easy to find other examples of statements that are likely to be true but unprovable. The essential trick is to find an infinite sequence of events, each of which might happen by accident, but with a small total probability for even one of them happening. Then the statement that none of the events ever happens is probably true but cannot be proved.”

In the spirit of Dyson’s observation, we shall give an example of a statement that is likely to be true and then take things one step further by presenting a formal proof that the statement is unprovable. Consider the following function:

Definition 1: Let $`T:`$ be a function such that $`T(n)=\frac{3n+1}{2}`$ if $`n`$ is odd and $`T(n)=\frac{n}{2}`$ if $`n`$ is even.

The Collatz $`3n+1`$ Conjecture states that for each $`n`$, there exists a $`k`$ such that $`T^{(k)}(n)=1`$, where $`T^{(k)}(n)`$ is the function $`T`$ iteratively applied $`k`$ times to $`n`$ . As of May 10, 2011, this conjecture has been verified for all positive integers up to about $`2^{60}`$ . Furthermore, one can give a heuristic probabilistic argument  that since every iterate of the function $`T`$ decreases on average by a multiplicative factor of about $`(\frac{3}{2})^{1/2}(\frac{1}{2})^{1/2}=`$ $`(\frac{3}{4})^{1/2}`$, all iterates will eventually converge into the infinite cycle $`\{1,2,1,2,\mathrm{}\}`$, assuming that each $`T^{(k)}`$ sufficiently mixes up $`n`$ as if each $`T^{(k)}(n)`$ (mod 2) were drawn at random from the set $`\{0,1\}`$.

However, the Collatz $`3n+1`$ Conjecture has never been formally proven. We shall prove that the Collatz $`3n+1`$ Conjecture can, in fact, never be formally proven, even though there is a lot of evidence for its truth. The underlying assumption in our argument is that any proof of a theorem can be written in a computer text-file, which is composed of bits (zeroes and ones). First, let us present a definition of “random”.

Definition 2: We shall say that vector $`𝐱\{0,1\}^k`$ is random if $`𝐱`$ cannot be specified in less than $`k`$ bits in a computer text-file .

For example, the vector of one million concatenations of the vector $`(0,1)`$ is not random, since we can specify it in less than two million bits in a computer text-file by just writing, “the vector of one million concatenations of the vector $`(0,1)`$” in the text-file. However, the vector of outcomes of one million coin-tosses has a good chance of fitting our definition of “random”, since much of the time the most compact way of specifying such a vector is to simply make a list of the results of each coin-toss, in which one million bits are necessary. We now prove three theorems.

Theorem 1: For any vector $`𝐱\{0,1\}^k`$, there exists an $`n`$ such that $`𝐱=(n,T(n),\mathrm{},T^{(k1)}(n))`$ (mod 2).

Proof: A proof of this can be found in “The $`3x+1`$ problem and its generalizations” .∎

Theorem 2: If $`k,n`$ and $`T^{(k)}(n)=1`$, then in order to prove that $`T^{(k)}(n)=1`$, it is necessary to specify the values of $`(n,T(n),\mathrm{},T^{(k1)}(n))`$ (mod 2) in the proof.

Proof: Let the vector $`(x_0(n),x_1(n),\mathrm{},x_{k1}(n))`$ equal $`(n,T(n),\mathrm{},T^{(k1)}(n))`$ (mod 2). Then notice that the formula, $`T^{(k)}(n)=\lambda _k(n)n+\rho _k(n)`$ , where

$$\lambda _k(n)=\frac{3^{x_0(n)+\mathrm{}+x_{k1}(n)}}{2^k}$$

and

$$\rho _k(n)=\underset{i=0}{\overset{k1}{}}x_i(n)\frac{3^{x_{i+1}(n)+\mathrm{}+x_{k1}(n)}}{2^{ki}},$$

is determined by the values of $`(n,T(n),\mathrm{},T^{(k1)}(n))`$ (mod 2) and there is a one-to-one correspondence between all of the possible formulas for $`T^{(k)}(n)`$ and all of the possible values of $`(n,T(n),\mathrm{},T^{(k1)}(n))`$ (mod 2); therefore, in order to prove that $`T^{(k)}(n)=1`$, it is necessary to specify the values of $`(n,T(n),\mathrm{},T^{(k1)}(n))`$ (mod 2) in the proof, since in order to prove that $`T^{(k)}(n)=1`$, it is necessary to specify the formula for $`T^{(k)}(n)`$ in the proof. ∎

Theorem 3: It is impossible to prove the Collatz $`3n+1`$ Conjecture.

Proof: Suppose that there exists a proof of the Collatz $`3n+1`$ Conjecture, and let $`L`$ be the number of bits in such a proof. Now, let $`𝐱\{0,1\}^{L+1}`$ be a random vector, as defined above. (It is not difficult to prove that at least half of all vectors in $`\{0,1\}^{L+1}`$ are random .) By Theorem 1, there exists an $`n`$ such that $`𝐱=(n,T(n),\mathrm{},T^{(L)}(n))`$ (mod 2) and $`T^{(L+1)}(n)=T^{(L)}(n)`$ (mod 2). Then $`T^{(L)}(n)>2`$, so if $`T^{(k)}(n)=1`$, then $`k>L`$. Hence, by Theorem 2 it is necessary to specify the values of $`(n,T(n),\mathrm{},T^{(L)}(n))`$ (mod 2) in order to prove that there exists a $`k`$ such that $`T^{(k)}(n)=1`$. But since $`(n,T(n),\mathrm{},T^{(L)}(n))`$ (mod 2) is a random vector, at least $`L+1`$ bits are necessary to specify $`(n,T(n),\mathrm{},T^{(L)}(n))`$ (mod 2), contradicting our assumption that the proof contains only $`L`$ bits; therefore, a formal proof of the Collatz $`3n+1`$ Conjecture cannot exist. ∎


text:
# I. Introduction

## I. Introduction

Gauge theory with $`𝒩=2`$ supersymmetry has a successful exact description of low energy effective action in terms of curve, the meromorphic differential, the period integrals and above all the prepotential . This structure is shared by the corresponding classical integrable system of particles, which may be regarded as representing universality class of the field theories of this type . In particular, the prepotential, which extracts all available information on low energy phenomena from the curve and the periods, can be extended to include infinite number of couplings as time variables  and becomes a generating functional (or $`\tau `$ function) of some kind from the point of view of the integrable hierarchy. It is therefore natural to investigate what role is played by the prepotential when the original theory is finitely as opposed to infinitesimally driven by a specific set of operators. The framework which fits this line of thoughts is deformation theory of prepotential called Whitham deformation. See  for reviews.

In more physical terms, this amounts to asking whether the notion of the prepotential be effective when the $`𝒩=2`$ theory becomes softly broken to $`𝒩=1`$  or $`𝒩=0`$. In fact there have been several such developments till now, a few of which are listed below. In , the Whitham deformation based on the formalism of  and the computation of  has been applied to $`𝒩=0`$ through the spurion formalism of . More recently, the gluino condensate prepotential has appeared through the determination of the effective superpotential for $`𝒩=1`$ super Yang-Mills with the adjoint chiral superfield . This prepotential has been identified as the matrix model free energy . In the subsequent development of , the prepotential appeared via a fermionic shift symmetry 111The developments concerning with this particular point include the work of .  associated with the system of three anomalous chiral Ward-Takahashi identities . See a recent review  for more on this last development, and  for recent discussions on the relationship of the superpotential with the Lax matrix.

In this paper, we will add one more step to this list of developments in prepotential theory: we will find a more direct route and rationale leading to the appearance of the prepotential for $`𝒩=1`$ theories. We begin with a few preliminaries. In the next section, we first give the basic ingredients of the Seiberg-Witten curve for the $`𝒩=2`$ $`SU(N)`$ pure-super Yang-Mills theory (or the spectral curve of the $`N`$ site periodic Toda chain.) Next, we briefly review Whitham deformation theory for the prepotential. It is a generic framework and we need to find a proper condition that deforms the theory and that breaks $`𝒩=2`$ supersymmetry to $`𝒩=1`$. At the same time, the superpotential responsible for this deformation must be identified. In section three, we carry out this task. We examine the mixed second derivatives. A main thrust of this paper is to observe that, in prepotential theory, the condition for a curve to degenerate or factorize is given by that the kernel of the matrix made of the mixed second derivatives of the Whitham prepotential be nontrivial:

$`\mathrm{ker}{\displaystyle \frac{^2}{a^iT_{\mathrm{}}}}\mathrm{𝟎}.`$ (1.1)
We derive the matrix model curve, identifying the tree-level superpotential with the vector belonging to this nontrivial kernel of the matrix above. In the classical limit where the scale parameter vanishes, expressions get simplified largely. We discuss this case separately in section four and show that the identification made in section three is indeed correct. We make several comments on the gluino condensate prepotential in the final section.

## II. Some Preliminaries

### A.  The Seiberg-Witten curve

Let us first collect some basic ingredients of the Seiberg-Witten curve for the $`𝒩=2`$ $`SU(N)`$ pure-super Yang-Mills theory (or the spectral curve of the $`N`$ site periodic Toda chain) and the meromorphic differential $`dS_{SW}`$ defined on this Riemann surface of genus $`g=N1`$. The curve can be represented by a hyperelliptic form:

$`Y^2`$ $`=`$ $`P_N(x)^24\mathrm{\Lambda }^{2N},`$ (2.1)
$`P_N(x)`$ $``$ $`det\left(x\mathrm{𝟏}\mathrm{\Phi }\right)={\displaystyle \underset{i=1}{\overset{N}{}}}\left(xp_i\right)=x^N{\displaystyle \underset{k=2}{\overset{N}{}}}u_kx^{Nk}`$ (2.2)
$`=`$ $`{\displaystyle \underset{k=0}{\overset{N}{}}}s_k(h_{\mathrm{}})x^{Nk},`$
where $`s_k(h_{\mathrm{}})`$ are the fundamental Schur polynomials of

$`h_{\mathrm{}}={\displaystyle \frac{1}{\mathrm{}}}tr_N\mathrm{\Phi }^{\mathrm{}}={\displaystyle \frac{1}{\mathrm{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}p_i^{\mathrm{}},\mathrm{}=2,\mathrm{}N.`$ (2.3)
The overall $`U(1)`$ can be decoupled from our consideration at this moment. In the Toda chain representation, this curve is parameterized by the spectral parameter $`z`$ by

$`P_N(x)`$ $`=`$ $`z+{\displaystyle \frac{\mathrm{\Lambda }^{2N}}{z}},`$
$`Y`$ $`=`$ $`z{\displaystyle \frac{\mathrm{\Lambda }^{2N}}{z}}.`$ (2.4)
The general variations of eq. (A.) give us the following equations:

$`P_N^{}\delta x+\delta _uP_N=Y\delta \mathrm{log}zN\left(YP_N\right)\delta \mathrm{log}\mathrm{\Lambda },`$
$`\delta Y=P_N\delta \mathrm{log}z+N\left(YP_N\right)\delta \mathrm{log}\mathrm{\Lambda }.`$ (2.5)
Here $`\delta _u`$ denotes a generic variation with respect to $`u_k`$. The distinguished meromorphic differential for the prepotential theory of the $`𝒩=2`$ $`SU(N)`$ pure-super Yang-Mills is

$`dS_{SW}`$ $`=`$ $`xd\mathrm{log}z=xt(x)dx.`$ (2.6)
$`t(x)`$ $`=`$ $`{\displaystyle \frac{P_N^{}}{\sqrt{P_N^24\mathrm{\Lambda }^{2N}}}}.`$ (2.7)
This differential possesses a double pole at $`x=\mathrm{}_\pm `$. Moduli derivatives with either $`z`$ or $`x`$ and $`\mathrm{\Lambda }`$ fixed generate the bases of the holomorphic differentials:

$$\frac{}{u_k}dS_{SW}_{z,\mathrm{\Lambda }}=\frac{x^{Nk}}{Y}dxdv_k,\frac{}{u_k}dS_{SW}_{x,\mathrm{\Lambda }}=dv_kd\left(\frac{x^{Nk+1}}{Y}\right),k=2,\mathrm{}N.$$
(2.8)
The differentials $`dv_k`$ are related by $`u_k`$ dependent numerical factors to the canonical holomorphic differentials $`d\omega _i`$ which are normalized as $`_{A_i}𝑑\omega _j=\delta _j^i`$. Moduli derivatives are in general coordinate dependent and in this paper, we take them, keeping $`z`$ fixed.

### B. Whitham Prepotential

It is possible to deform consistently both moduli of the Riemann surface and the meromorphic differential discussed in the last subsection without losing their defining properties. Namely,

$`dS_{SW}dS,{\displaystyle \frac{}{u_k}}dS_{z,\mathrm{\Lambda }}=\mathrm{holomorphic}.`$ (2.9)
This is called Whitham deformation (of Toda integrable hierarchy). It has an effect of adding higher order poles to the original double poles of the meromorphic differential $`dS_{SW}`$. The deformation is characterized by a set of the punctures and local coordinates in their neighborhood, which we denote generically by $`\xi `$. In order to accomplish this deformation, we introduce a set of meromorphic differentials $`d\mathrm{\Omega }_{\mathrm{}}`$ which satisfy

$`d\mathrm{\Omega }_{\mathrm{}}=\xi ^\mathrm{}1d\xi +\mathrm{nonsingular}\mathrm{part},\mathrm{}=1,2\mathrm{}.`$ (2.10)
This still leaves us with the possibility to add any linear combination of the holomorphic differentials $`d\omega _i`$ to the right hand side. To remove this ambiguity, we require the condition

$$_{A_i}𝑑\mathrm{\Omega }_{\mathrm{}}=0.$$
(2.11)
Construction of the Whitham prepotential begins with introducing time variables via

$`{\displaystyle \frac{dS}{T_{\mathrm{}}}}=d\mathrm{\Omega }_{\mathrm{}}.`$ (2.12)
Let

$$a^i_{A_i}𝑑S,$$
(2.13)
be the local coordinates of the moduli space. After some reasonings , one concludes

$`dS={\displaystyle \underset{i=1}{\overset{g}{}}}a^id\omega _i+{\displaystyle \underset{\mathrm{}}{}}T_{\mathrm{}}d\mathrm{\Omega }_{\mathrm{}},`$ (2.14)
and

$$T_{\mathrm{}}=res_{\xi =0}\xi ^{\mathrm{}}dS.$$
(2.15)
Both $`a^i`$ and $`T_{\mathrm{}}`$ are regarded as independent variables. In $`dS_{SW}`$, the first time alone is turned on. Invariant moduli $`h_k`$ are expressible as

$$h_k=h_k(a^i,T_{\mathrm{}}).$$
(2.16)
The prepotential $`(a^i,T_{\mathrm{}})`$ is introduced via

$`{\displaystyle \frac{}{a^i}}`$ $`=`$ $`{\displaystyle _{B^i}}𝑑S,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}1}iN1`$ (2.17)
$`{\displaystyle \frac{}{T_{\mathrm{}}}}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi i\mathrm{}}}res\xi ^{\mathrm{}}dS_{\mathrm{}}\left(h_k\right).\mathrm{}=1,2,\mathrm{}`$ (2.18)
The consistency of eq. (2.17) and that of eq. (2.18) are ensured respectively by the property that the period matrix is symmetric as well as by the Riemann identity. The right hand side of eq. (2.18), which we have denoted by $`_{\mathrm{}}\left(h_k\right)`$, is some polynomial  of the invariant moduli. Below, the local coordinates in the neighborhood of the punctures at $`x=\mathrm{}_\pm `$ are taken as

$$\xi =z^{\frac{1}{N}}.$$
(2.19)
## III. Mixed Second Prepotential Derivatives: Derivation of the Matrix Model Curve

Now we come to the point of our paper. Differentiating eqs. (2.17) and (2.18) once again, we obtain two distinct formulas for the mixed second derivatives:

$`{\displaystyle \frac{^2}{a^jT_{\mathrm{}}}}={\displaystyle _{B^j}}𝑑\mathrm{\Omega }_{\mathrm{}}={\displaystyle \frac{1}{2\pi i\mathrm{}}}res_{\xi =0}\xi ^{\mathrm{}}d\omega _j,`$ (3.1)
where $`j=1N1`$. As for $`\mathrm{}`$, it is any positive integer. Note that eq. (3.1) has been given with no reference to $`dS`$. The discussion in what follows holds, therefore, for the Whitham prepotential deformed by the arbitrary values of the time variables. As is stated in the introduction, we impose the condition that the kernel of this rectangular matrix be nontrivial:

$`\mathrm{ker}{\displaystyle \frac{^2}{a^iT_{\mathrm{}}}}\mathrm{𝟎}.`$ (3.2)
Equivalently

$$\mathrm{rank}\frac{^2}{a^iT_{\mathrm{}}}N2.$$
(3.3)
In the case in which the time variables are truncated to $`T_1,T_2\mathrm{}T_{N1}`$, the condition is of course stated as that of the vanishing determinant:

$`det{\displaystyle \frac{^2}{a^iT_{\mathrm{}}}}=0.`$ (3.4)
Eq. (3.2) tells us the existence of at least one nonvanishing column vector

$`\left(\begin{array}{c}c^1\\ \mathrm{}\\ c^{N1}\\ \mathrm{}\end{array}\right)`$ (3.9)
belonging to the kernel;

$`0={\displaystyle \underset{\mathrm{}}{}}{\displaystyle \frac{^2}{a^iT_{\mathrm{}}}}c^{\mathrm{}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{}}{\displaystyle _{B^i}}c^{\mathrm{}}𝑑\mathrm{\Omega }_{\mathrm{}}`$ (3.10)
$`=`$ $`{\displaystyle \frac{1}{2\pi i}}res_{\xi =0}\left({\displaystyle \underset{\mathrm{}}{}}{\displaystyle \frac{c_{\mathrm{}}}{\mathrm{}}}\xi ^{\mathrm{}}d\omega _i\right).`$ (3.11)
Eq. (3.10) implies the existence of the meromorphic one-form

$$d\stackrel{~}{\mathrm{\Omega }}\underset{\mathrm{}}{}c^{\mathrm{}}d\mathrm{\Omega }_{\mathrm{}}.$$
(3.12)
whose period integral over any of the $`A_i`$ and $`B^i`$ cycles vanishes. Let us argue implications of this last statement. Obviously, once this property holds, we can integrate $`d\stackrel{~}{\mathrm{\Omega }}`$ along any path ending with a point $`z`$ to define a function holomorphic except at punctures: $`f(z)=^z𝑑\stackrel{~}{\mathrm{\Omega }}`$. As for the order of the poles at the punctures, it is generically arbitrary according to the construction. This is, however, contradictory to the Weierstrass gap theorem , which says that, for $`g=N1`$ integers satisfying $`1=n_1<n_2<\mathrm{}<n_g<2g`$, such function with a pole of order $`n_j`$ does not exist. This theorem is derived from the Riemann-Roch theorem. In order to avoid the contradiction, the degeneration of the surface must take place.

Eq. (3.2) also implies the existence of at least one non-vanishing co-vector (row vector)

$$(\stackrel{~}{c}_1,\stackrel{~}{c}_2,\mathrm{},\stackrel{~}{c}_{N1})$$
(3.13)
such that

$$0=\underset{i=1}{\overset{N1}{}}\stackrel{~}{c}_i\frac{^2}{a^iT_{\mathrm{}}}=\underset{i=1}{\overset{N1}{}}\stackrel{~}{c}_i\frac{_{\mathrm{}+1}}{a^i}.$$
(3.14)
We see from eq. (2.18) that this is satisfied provided

$$\underset{i=1}{\overset{N1}{}}\stackrel{~}{c}_i\frac{h_k}{a^i}=0,k=2N.$$
(3.15)
Eq. (3.15) tells us the degeneration of the curve from the point of view of the moduli space: the invariant moduli $`h_k,k=2N`$ become functions which actually depend on less than $`N1`$ arguments. This equation is regarded as the counterpart of the statement of the vanishing discriminant .

We conclude that the original curve must get degenerated under the condition eq. (3.2) and the factorization of eq. (2.1) by the following kind  takes place:

$`Y^2(x)`$ $`=`$ $`H_{Nn}^2(x)F_{2n}(x).`$ (3.16)
$`P_N^{}(x)`$ $`=`$ $`H_{Nn}(x)R_{n1}(x)`$ (3.17)
$`t(x)`$ $`=`$ $`{\displaystyle \frac{R_{n1}(x)}{\sqrt{F_{2n}(x)}}}.`$ (3.18)
The polynomial $`R_{n1}`$ does not depend on $`N`$ . As is well-known, $`Nn`$ non-intersecting cycles vanish once this condition is satisfied and the moduli space of the resulting vacua is codimension $`Nn`$ subspace of the original $`𝒩=2`$ Coulomb branch. To smoothen out the singularities developed, let us imagine blowing up the surface and multiplying $`dv_k`$ in eq. (2.8) by a factor which cancels the poles developed and which behaves as $`x`$ at the infinities. Upon partial fractions, $`n`$ independent differentials emerge from $`dv_k,k=2,\mathrm{}N`$:

$$\frac{x^{j1}dx}{\sqrt{F_{2n}}},j=1n.$$
(3.19)
The differentials for $`j=1n1`$ are the holomorphic differentials on this reduced surface and the last one $`j=n`$ has a pole at the infinities. This last one has been included due to the blow-up process and physically implies that the overall $`U(1)`$ fails to decouple through this process.

Finally let us examine eq. (3.11). Let

$$\underset{\mathrm{}}{}\left(\frac{c_{\mathrm{}}}{\mathrm{}}\xi ^{\mathrm{}}\right)_+W_{K+1}^{}(x),$$
(3.20)
where $`(\mathrm{})_+`$ denotes the part consisting of the non-negative powers of $`\mathrm{}`$ in the Laurent expansion in $`x`$. We consider the case in which $`W_{K+1}^{}(x)`$ is a polynomial of degree $`K`$ $`(Kn)`$ in $`x`$ parameterized as

$$W_{K+1}^{}(x)\underset{j=1}{\overset{K}{}}\left(x\alpha _j\right).$$
(3.21)
We substitute eq. (3.19) for $`d\omega _i`$ in eq. (3.11) to obtain

$$0=res_{x=\mathrm{}}\left(W_{K+1}^{}(x)\frac{x^{j1}}{\sqrt{F_{2n}}}dx\right),j=1n.$$
(3.22)
Hence

$$\frac{W_{K+1}^{}(x)}{\sqrt{F_{2n}(x)}}=Q_{Kn}(x)+\underset{\mathrm{}>n}{}\frac{\beta _{\mathrm{}}}{x^{\mathrm{}}},$$
(3.23)
where $`Q_{Kn}(x)`$ is a polynomial of degree $`Kn`$. We conclude

$$F_{2n}(x)Q_{Kn}^2(x)=W_{K+1}^2(x)+f_{K1}(x),$$
(3.24)
where $`f_{K1}(x)`$ is a polynomial of degree $`K1`$. It still remains to be seen that the function $`W_{K+1}(x)`$ introduced above is in fact a tree-level superpotential. This is done in the next section. Let us also note that our discussion suggests a family of superpotentials continuously connected when the kernel is more than one-dimensional.

## IV. Classical Limit

In the classical limit, $`\mathrm{\Lambda }=0`$ and the curve gets simplified largely:

$$Y=z=\underset{\mathrm{}=1}{\overset{N}{}}(xp_{\mathrm{}}).$$
(4.1)
The original Seiberg-Witten differential becomes

$$dS_{SW}^{(class)}=xd\mathrm{log}z=\underset{i=1}{\overset{N}{}}\frac{x}{xp_i}dx.$$
(4.2)
The period integrals over the $`A_i`$ cycles just pick up the poles at $`p_i`$:

$$a_i^{(class)}=p_i.$$
(4.3)
The canonical holomorphic differentials take the following form:

$$d\omega _i^{(class)}=\frac{}{p_i}dS_{SW}^{(class)}_z=\frac{dx}{xp_i}.$$
(4.4)
Eq.(2.14) is simply

$$dS_{SW}^{(class)}=\underset{i=1}{\overset{N}{}}p^id\omega _i+Ndx.$$
(4.5)
Let us suppose that

$$z=\underset{j=1}{\overset{n}{}}(x\beta _j)^{N_j},\underset{j=1}{\overset{n}{}}N_j=N,$$
(4.6)
which means that $`N_j`$ poles of eq. (4.2) (or $`N_j`$ eigenvalues of the vev of the adjoint Higgs) coalesce to one point $`\beta _j`$ for $`j=1n`$ and $`(Nn)`$ non-intersecting $`A_i`$ cycles vanish. The canonical holomorphic differentials on this degenerate curve are

$$d\omega _j^{(class,red)}=\frac{dx}{x\beta _j},j=1n,$$
(4.7)
so that

$$a_j^{(red)}=\beta _j.$$
(4.8)
Let us now examine the condition eq. (3.22). In the classical limit, eq. (3.22) becomes

$$0=res_{x=\mathrm{}}\left(W_{K+1}^{}(x)d\omega _j^{(class,red)}\right),j=1n.$$
(4.9)
The residue is originally evaluated by the contour around infinity but by a contour deformation it becomes a residue at $`x=\beta _j`$. We see that $`\beta _j`$ must coincide with one of the roots $`\alpha _j`$ of $`W_{K+1}^{}`$. The vacuum values of the adjoint scalars are constrained to the extremum of $`W_{K+1}`$. We conclude that the function $`W_{K+1}`$ is the tree-level superpotential as is promised.

## V. Discussion

Once the condition (3.2) is imposed upon, we have a reduced curve of genus $`g=n1`$. Let us set $`K=n`$ for simplicity. In the context of the present paper, the coefficients (denoted by $`b_{\mathrm{}}`$) of the polynomial $`f_{K1}`$ are determined by the parameters $`\alpha _j`$ of the superpotential. In , these $`b_{\mathrm{}}`$ as opposed to $`u_k`$ are treated as moduli. The gluino condensate prepotential has been given by the matrix model differential whose $`b_{\mathrm{}}`$ moduli derivatives are almost holomorphic and these moduli derivatives are taken, keeping the coordinate $`x`$ fixed as opposed to $`z`$ fixed. In this setting the Whitham times have been identified with the symmetric polynomials made of $`\alpha _j`$. Note that the $`S_i`$ moduli ( the $`A_i`$ cycle integral of the matrix model differential) as opposed to $`a_i`$ moduli have a quantum mechanical origin: they vanish as $`\mathrm{\Lambda }0`$. Also $`S{\displaystyle \underset{i=1}{\overset{n}{}}}S_i`$ is non-vanishing and the cutoff must be introduced at the infinities of the curve.

## VI. Acknowledgements

The authors thank E. Date, K. Ito, A. Morozov, N. Sakai, K. Takasaki and Y. Yasui for helpful discussions and useful comments. This work is supported in part by the Grant-in-Aid for Scientific Research(14540264, 14540073) from the Ministry of Education, Science and Culture, Japan. Supports from the 21st century COE programs at Osaka City University (H.I.) and at Nagoya University (H.K) are gratefully appreciated.


text:
# Flat histogram version of the pruned and enriched Rosenbluth method

## Abstract

In this letter we present a flat histogram algorithm based on the pruned and enriched Rosenbluth method (PERM). This algorithm incorporates in a straightforward manner microcanonical reweighting techniques, leading to “flat histogram” sampling in the chosen parameter space. As an additional benefit, our algorithm is completely parameter free and, hence, easy to implement. We apply this algorithm to interacting self-avoiding walks (ISAW), the generic lattice model of polymer collapse.

Recently, there has been revived interest in flat histogram algorithms wang2001 , which strive to evenly sample configuration space with respect to a chosen parameterization, e.g. microcanonical energy. These algorithms are particular implementations of “umbrella sampling” torrie1977 , in which the configuration space is sampled according to a given probability distribution, the so-called “umbrella”. This umbrella distribution is generally chosen such that the whole configuration space of interest is accessible in one simulation. One major difficulty hereby is to find a suitable umbrella distribution.

There has also been an exciting development in stochastic growth algorithms, which are based on the Rosenbluth and Rosenbluth algorithm rosenbluth1955 . If this algorithm, which kinetically grows configurations, gets enhanced by cleverly chosen enrichment and pruning steps grassberger1997 , one obtains the pruned and enriched Rosenbluth method (PERM), a powerful algorithm for, e.g., simulation of the polymer collapse transition.

We present in this letter a new algorithm, flatPERM, which is a combination of these two types of algorithms, i.e. a flat histogram version of the pruned and enriched Rosenbluth method.

As opposed to earlier work in this direction bachmann2003 , in which an iterative scheme similar to the multicanonical algorithm berg1991  was used, we utilize the self-tuning capabilities of PERM directly. This leads to a considerable simplification of the algorithm.

While flatPERM includes umbrella sampling ideas, it is stricly speaking not a multicanonical method, as multicanonical sampling conventionally describes a particular iterative version of adaptive umbrella sampling, in which first an umbrella distribution is obtained iteratively and then a final simulation is performed with a fixed umbrella distribution.

This letter is structured as follows: we first give a pedagogical introduction to PERM, which will then allow us to introduce flatPERM as a seemingly trivial extension. As an application, we present simulations of interacting self-avoiding walks (ISAW) on the square lattice and the simple cubic lattice. We conclude with a description of further applications.

We will consider a rather abstract setting of configurations with a certain size $`n`$, which is parameterized with an additional variable $`m`$. In general, one can even consider a set of variables $`m_i`$, but for pedagogical reasons we shall in this paper restrict ourselves to the case of $`m`$ corresponding to an energy $`E=ϵm`$. Both $`n`$ and $`m`$ are assumed to have non-negative integer values. Moreover, we will need the notion of “atmosphere” of a configuration, which is the number $`a`$ of different ways to continue to grow this configuration, and is also a non-negative integer.

While it is useful to present the algorithm in such an abstract setting, it may help the reader to keep the application to interacting self-avoiding walks (ISAW) on a regular lattice in mind. In this case, the size of the configuration is the number of steps of the walk, the energy is the number of non-consecutive nearest-neighbor bonds, and the atmosphere is the number of non-occupied sites around the endpoint of the walk. If all the sites around the endpoint are occupied, then the atmosphere is zero and the walk cannot be continued.

The basis of the algorithm is the Rosenbluth and Rosenbluth algorithm, a stochastic growth algorithm in which the configurations of interest are grown from scratch. The growth is kinetic, which is to say that each growth step is selected at random from all possible growth steps. Thus, if there are $`a`$ possible ways of growth, one selects one of them with probability $`p=1/a`$. As this number generally changes during the growth process, different configurations are thus generated with different probabilities, and one needs to resort to reweighting techniques.

It is advantageous to view this algorithm as approximate counting, in which case the weights of the configurations serve as estimates of the number of configurations. To understand this point of view, imagine that we were to perform a complete enumeration of the configuration space. Doing this requires that at each growth step we would have to choose all the possible continuations and count them each with equal weight. If we now select fewer configurations, we have to change the weight of these configurations accordingly, in order to correct for missing out on some configurations. Thus, if we choose one growth step out of $`a`$ possible ones, we replace $`a`$ configurations with equal weight by one “representative” configuration with $`a`$-fold weight. In this way, the weight of each grown configuration is a direct estimate of the total number of configurations.

Let the atmosphere $`a_n=a(\phi _n)`$ be the number of distinct ways in which a configuration $`\phi _n`$ of size $`n`$ can be extended. Then, the weight associated with a configuration of size $`n`$ is the product of all the atmospheres $`a_k`$ encountered while growing this configuration, i.e.

$$W=\underset{k=0}{\overset{n1}{}}a_k.$$
(1)
After having started $`S`$ growth chains, an estimator $`C_n^{est}`$ for the total number of configurations $`C_n`$ (the “infinite-temperature” partition function, in which all configurations appear with equal weight $`1`$) is given by the mean over all generated samples $`\phi _n^{(i)}`$ of size $`n`$ with respective weights $`W_n^{(i)}`$, i.e.

$$C_n^{est}=W_n=\frac{1}{S}\underset{i}{}W_n^{(i)}.$$
(2)
Here, the mean is taken with respect to the total number of growth chains $`S`$, and not the number of configurations which actually reach size $`n`$. Configurations which get trapped before they reach size $`n`$ appear in this sum with weight zero.

The problem with the Rosenbluth and Rosenbluth algorithm is two-fold. Firstly, the weights can vary widely in magnitude, so that the mean can get to be dominated by very few samples with very large weight. Secondly, regularly occurring trapping events, i.e. generation of configurations with zero atmosphere, can lead to exponential “attrition”, i.e. exponentially strong suppression of configurations of large sizes.

To overcome both of these problems, enrichment and pruning steps have been added to this algorithm, leading to the pruned and enriched Rosenbluth method (PERM) grassberger1997 . The basic idea is that one wishes to suppress fluctuations in the weights $`W_n^{(i)}`$, as these should on average be equal to $`C_n`$. Therefore, if the weight is too large, one enriches, i.e. one makes copies of the configuration and reduces the weight accordingly. On the other hand, if the weight is too small, one prunes probabilistically, and, in case the pruning is unsuccessful, continues growing with appropriately increased weight. Note that $`S`$, and therefore the expression (2) for the estimate $`C_n^{est}`$, is not changed by either enriching or pruning steps.

We need to specify enrichment and pruning criteria as well as the actual enrichment and pruning processes. While the idea of PERM itself is straightforward, there is now a lot of freedom in the precise choice of the pruning and the enrichment steps. The “art” of making PERM perform efficiently is based to a large extent on a suitable choice of these steps. Distilling our own experience with PERM, we present here what we believe to be an efficient, and, most importantly, parameter free version.

In contrast to other expositions of PERM (e.g. hsu2003 ), we propose to prune or enrich constantly, to enable larger exploration of the configuration space (the motivation will become clear once flatPERM is introduced below). Define $`r`$ as the ratio of weight and estimated number of configurations, $`r=W_n^{(i)}/C_n^{est}`$. Then we enrich if $`r>1`$ and prune if $`r<1`$. Moreover, the actual pruning and enrichment steps are chosen such that the weights are set as closely as possible to $`C_n^{est}`$ to minimize fluctuations:

* $`r>1`$ $``$ enrichment step:
make $`c=\mathrm{min}(r,a_n)`$ distinct copies, each with weight $`\frac{1}{c}W_n^{(i)}`$ (as in nPERM hsu2003a ).
* $`r<1`$ $``$ pruning step:
continue growing with probability $`r`$ and weight $`C_n^{est}`$ (i.e. prune with probability $`1r`$).

Note that we perform pruning and enrichment after the configuration has been included in the calculation of $`C_n^{est}`$ and is used for weights during the next growth step. Initially, the estimates $`C_n^{est}`$ can of course be grossly wrong, as the algorithm knows nothing about the system it is simulating. However, even if initially “wrong” estimates are used for pruning and enrichment, in all applications considered the algorithm can be seen to converge to the true values. It is, in a sense, self-tuning.

At this point it is now straight-forward to change PERM to a thermal ensemble with inverse temperature $`\beta =1/k_BT`$ and energy $`E=ϵm`$ by multiplying the weight with the appropriate Boltzmann factor $`\mathrm{exp}(\beta E)`$, which leads to an estimate of the partition function $`Z_n(\beta )`$ of the form

$$Z_n^{est}(\beta )=W\mathrm{exp}(\beta E)_n.$$
(3)
The pruning and enrichment procedures are changed accordingly, replacing $`W`$ by $`W\mathrm{exp}(\beta E)`$ and $`C_n^{est}`$ by $`Z_n^{est}(\beta )`$, and using $`r=W_n^{(i)}\mathrm{exp}(\beta E_m^{(i)})/Z_n^{est}(\beta )`$. (It is in this setting that PERM is usually described.)

Alternatively, however, it is also possible to consider microcanonical estimators for the total number $`C_{n,m}`$ of configurations of size $`n`$ with energy $`m`$ (i.e. the “density of states”). An appropriate estimator $`C_{n,m}^{est}`$ is then given by the mean over all generated samples $`\phi _{n,m}^{(i)}`$ of size $`n`$ and energy $`m`$ with respective weights $`W_{n,m}^{(i)}`$, i.e.

$$C_{n,m}^{est}=W_{n,m}=\frac{1}{S}\underset{i}{}W_{n,m}^{(i)}.$$
(4)
Again, the mean is taken with respect to the total number of growth chains $`S`$, and not the number of configurations $`S_{n,m}`$ which actually reach a configuration of size $`n`$ and energy $`m`$. The pruning and enrichment procedures are also changed accordingly, replacing $`C_n`$ by $`C_{n,m}`$ and using $`r=W_{n,m}^{(i)}/C_{n,m}^{est}`$.

At this point it is worth pointing out that the pruning and enrichment criterion for PERM leads to a roughly constant number of samples being generated at each size $`n`$ for PERM. In fact, one can motivate the given pruning and enrichment criteria by stipulating that one wishes to have a roughly constant number of samples, as this leads to the algorithm performing an unbiased random walk in the configuration size. Correspondingly, in the flat-energy version the algorithm performs an unbiased random walk in both size and energy of the configuration, and we obtain a roughly constant number of samples for each size $`n`$ and energy $`m`$.

It is because of the fact that the number of samples is roughly constant in each histogram entry, that this algorithm can be seen as a “flat histogram” algorithm, which we consequently call flatPERM. In hindsight in becomes clear that PERM itself can be seen as a flat histogram algorithm, at it creates a roughly flat histogram in size $`n`$. Recognizing this led us to the formulation of this algorithm in the first place.

At this point we return to our discussion of the pruning and enrichment strategies. For PERM it may be more advantageous to allow for a high diffusivity along the size $`n`$. This is done by minimizing pruning and enrichment, i.e. at the cost of allowing larger weight fluctuations. For flatPERM on the other hand we need to achieve diffusive behavior also with respect to the energy variable $`m`$. Precisely this is achieved by allowing for much enrichment (and thus necessarily also pruning), as each set of configurations enriched at size $`n`$ will contribute to a range of different energies $`m`$ at size $`n+1`$.

Even though we view flatPERM as a flat histogram algorithm, we have not yet explicitly conditioned the pruning and enrichment with respect to the local number of samples $`S_{n,m}`$ created at each size $`n`$ and energy $`m`$. This can be done by multiplying $`r`$ by $`S/S_{n,m}`$, as this enhances enrichment if $`S_{n,m}<S`$, i.e. when there are too few samples, and vice versa.

One general problem with PERM is that enrichment generates large correlation between samples, so that it would be useful to replace the number $`S_{n,m}`$ of samples by a number $`S_{n,m}^{eff}`$ of effectively independent samples, thereby taking account of some autocorrelation time. This can be done heuristically by considering the average number of independent growth steps in a configuration. If the last enrichment has occurred at size $`n_{enr}`$, then a configuration of size $`n`$ and energy $`m`$ has $`n_{ind}=nn_{enr}`$ independent steps. The frequency of independent steps in a configuration is $`w_{n,m}^{(i)}=n_{ind}/n`$. We thus use

$$r=\frac{S}{S_{n,m}^{eff}}\frac{W_{n,m}^{(i)}}{C_{n,m}^{est}}\text{with}S_{n,m}^{eff}=\underset{i}{}w_{n,m}^{(i)}$$
(5)
for the pruning and enrichment criterion in our algorithm.

Now that one has a flat histogram method of PERM for the whole range of energies, one can easily modify the algorithm further to sample only a selected range by dividing $`r`$ by a “profile shape” $`f_{n,m}`$, leading to pruning whenever $`f_{n,m}`$ is close to zero. This is advantageous if one only wants to explore a restricted region of parameter space.

Even though the algorithm described thus far is free of parameters, there is a technical problem due to the fact that initial weights are grossly wrong, which can lead to overflow problems. This can easily be overcome by initially restricting the maximal size of grown configurations, e.g. by limiting the maximal size by the number of started growth chains,

$$n<cS.$$
(6)
The relevant number of growth chains at size $`n`$ is thus reduced to $`Sn/c`$. For ISAW we find that a choice of delay of $`c10`$ is sufficient.

Averages of observables $`Q`$ defined on the set of configurations can now be obtained by storing weighted sums of these observables, from which one obtains

$$Q_{n,m}^{est}=\frac{QW_{n,m}}{W_{n,m}}=\frac{_iQ_{n,m}^{(i)}W_{n,m}^{(i)}}{_iW_{n,m}^{(i)}}.$$
(7)
These can then be used for subsequent evaluations. For instance, the expectation value of $`Q`$ in the canonical ensemble at a given temperature $`\beta `$ can now be computed via

$$Q_n^{est}(\beta )=\frac{_mQ_{n,m}^{est}C_{n,m}^{est}\mathrm{exp}(\beta E_m)}{_mC_{n,m}^{est}\mathrm{exp}(\beta E_m)}.$$
(8)
We have implemented this algorithm for interacting self-avoiding walks on the square and simple cubic lattice. In both two and three dimensions, we have simulated walks up to length $`n=1024`$. Figures 1 and 2 clearly show the strength of the method. Figure 1 shows the number of configurations $`C_{n,m}`$, which vary over several hundred orders of magnitude. This range would have been inaccessible during one simulation with any canonical method. As an aside, we note that we had to rescale weights during the run to avoid overflow berg2003 . Additionally, we chose a delay of $`c=10`$ to stabilize the algorithm.

As can clearly be seen in Figure 2, the number of samples generated is roughly constant, irrespective of whether one flattens with respect to $`S_{n,m}`$ (top graph) or $`S_{n,m}^{eff}`$ (middle and bottom graphs). Only results for the square lattice are shown, as we find the same behavior for the simple cubic lattice. Using $`S_{n,m}^{eff}`$ for the flatness criterion leads moreover to an increased sampling of walks with very few and very many interactions, thus overcoming the usual difficulty of obtaining configurations with a large number of interactions. (The largest energy state gets repeatedly sampled in simulations in both dimensions.)

Once the simulations have been completed, thermodynamic quantities of interest can easily be computed. As an example, the specific heat curves $`C_n=k_B(\beta ϵ)^2\sigma ^2(m)/n`$ near the transition temperature are shown for both systems in Figure 3.

We have also tested our algorithm on the HP model dill1985 , which is a toy model for proteins. It consists of a self-avoiding walk with two types of monomers along the sites visited by the walk, which are either hydrophobic (type H) or polar (type P). One considers monomer-specific interactions, mimicking the interaction with a polar solvent such as water. The interaction strengths are thus chosen such that HH-contacts are favored, e.g. by choosing $`ϵ_{HH}=1`$ and $`ϵ_{HP}=ϵ_{PP}=0`$. The central question is to determine the density of states (and to find the ground state with lowest energy) for a given sequence of monomers. For various sequences taken from the literature we have confirmed previous density of states calculations and obtained identical ground state energies. The sequences we considered had $`n=58`$ steps ($`3`$ dimensions, $`E_{min}=44`$) and $`n=85`$ steps ($`2`$ dimensions, $`E_{min}=53`$) from hsu2003 , and $`n=80`$ steps ($`3`$ dimensions, $`E_{min}=98`$) from frauenkron1998 . We studied also a particularly difficult sequence with $`n=88`$ steps ($`3`$ dimensions, $`E_{min}=73`$) from beutler1996 , but the lowest energy we obtained was $`E=69`$. However, in this case no lower energy has been found with any other PERM implementation either, see hsu2003 .

To summarize, we have presented a new algorithm, flatPERM, which is a flat histogram version of a stochastic growth algorithm. This algorithm can in principle be applied to any statistical mechanical system for which configurations can be grown in a well-defined manner. Next to the presented applications of linear polymers, this algorithm can be applied to more complicated systems, such as lattice models of branched polymers rechnitzer2004 . Extensions to models with two energy parameters are in preparation, e.g. to the problem of adsorbing interacting polymers owczarek2004  or an extended Domb-Joyce model of polymer collapse krawczyk2004 .

The authors thank Andrew Rechnitzer and Aleks Owczarek for helpful discussions, Peter Grassberger for pointing out an incorrect formula, the referees for helpful comments, and the DFG for financial support.


text:
# Updated dispersion-theoretical analysis of the nucleon electromagnetic form factorsWork supported in part by the DFG under contract HA 3203/2-1.

## 1 Introduction

The electromagnetic form factors of the nucleon encode information about the structure of the nucleon over a wide range of scales. Depending on the kinematical conditions, they can be determined most precisely from elastic electron scattering using Rosenbluth-separation or in double polarization experiments. Only in the last decade has it become technically feasible to perform the double polarization experiments. Thus the data base has been considerably enlarged since our last dispersion analysis was performed in 1996 MMD ; HMD . This holds in particular for the neutron form factors, which can only be measured indirectly using deuterium or 3He targets. For recent reviews, see e.g. Refs. Gao ; Schmieden . At low momentum transfer, one obtains information about the various nucleon rms-radii, which are not only of interest by themselves but e.g. the proton charge radius $`r_E^p`$ has also to be known to a good accuracy to perform precision tests of QED in Lamb-shift measurements, or conversely, one can use such experiments to pin down $`r_E^p`$. For some time, low-energy electron scattering and Lamb shift determinations led to different values for $`r_E^p`$, but this discrepancy has been resolved, see Refs. Rosenfelder ; Sick:2003gm  and references therein. One can also extract vector meson–nucleon coupling constants and eventually study the transition from the non-perturbative regime of QCD to the perturbative one. Since the data span the range of four-momentum transfers from $`Q^20`$ to $`Q^231`$GeV2, the only model-independent method to analyze these data is dispersion theory. The most recent dispersion-theoretical analysis of the nucleon electromagnetic form factors dates back to almost a decade ago, including fits to space-like MMD  as well as time-like data HMD . For some recent vector-meson pole-model approaches to the nucleon electromagnetic form factors, see Refs. Lomo ; Dubni . In view of the new data and a new data base collected and discussed in Ref. FW , it seems timely to update the work of Refs. MMD ; HMD  using the data basis of FW . It is also interesting to investigate whether such a general scheme can lead to the pronounced structure around $`Q^2=0.2`$GeV2 found for all four form factors in the phenomenological analysis of Ref. FW .

## 2 Formalism

The electromagnetic structure of the nucleon is parameterized by the Dirac ($`F_1`$) and Pauli ($`F_2`$) form factors for the proton and the neutron

$$F_i^{p/n}(Q^2)=F_i^S(Q^2)\pm F_i^V(Q^2),i=1,2,$$
(1)
with $`Q^2`$ the four–momentum transfer of the virtual photon (the photon virtuality). Our conventions are such that $`Q^2>0`$ for space-like momentum transfer. We have expressed the nucleon form factors in the isospin basis ($`S=`$ isoscalar, $`V=`$ isovector) which is most appropriate for the dispersive analysis. The experimental data are usually given for the Sachs form factors, which are linear combinations of $`F_1`$ and $`F_2`$:

$`G_E^I(Q^2)`$ $`=`$ $`F_1^I(Q^2){\displaystyle \frac{Q^2}{4m^2}}F_2^I(Q^2),`$ (2)
$`G_M^I(Q^2)`$ $`=`$ $`F_1^I(Q^2)+F_2^I(Q^2),I=S,V.`$ (3)
The analysis of the nucleon electromagnetic form factors proceeds most directly through the spectral representation given by111Note that we work with unsubtracted dispersion relations. Since the normalizations of the various form factors are known, one could also work with once–subtracted dispersion relations.

$$F_i^I(Q^2)=\frac{1}{\pi }_{(\mu _0^I)^2}^{\mathrm{}}\frac{\sigma _i^I(\mu ^2)d\mu ^2}{\mu ^2+Q^2},i=1,2,I=S,V,$$
(4)
in terms of the real spectral functions $`\sigma _i^I(\mu ^2)=\mathrm{Im}F_i^I(\mu ^2)`$, and the corresponding thresholds are given by $`\mu _0^S=3M_\pi `$, $`\mu _0^V=2M_\pi `$. The spectral functions encode the pertinent physics of the nucleon form factors. In the isovector channel, the spectral function is build up by the two–pion continuum (including the $`\rho `$–resonance) as given by unitarity Hohler:1974ht  plus a series of poles, whose masses and residues are fit parameters. In the isoscalar channel, we only have poles, where the lowest two are given by the $`\omega `$ and the $`\varphi `$ mesons, respectively. We restrict the number of poles in both channels by the stability condition of Ref. SabbaStefanescu:1978sy  which requires to use the minimum number of vector meson poles necessary to fit the data. Note that this condition was also used in the earlier fits Hohler:1974ht ; MMD ; HMD . Enforcing the correct normalization conditions for all form factors, the experimental value of the neutron charge radius and the superconvergence relations from the perturbative QCD behaviour of the form factors, e.g. $`F_i(Q^2)1/(Q^2)^{i+1}`$ as $`Q^2`$ tends to infinity, reduces the number of free fit parameters considerably. For more details on the form of the spectral functions, the stability condition, and the way the various constraints are included, see Ref. MMD .

It is possible to take all pole masses from physical particles (except for one which is determined by the constraints discussed above). As a consequence, we have 2 (3) free parameters in the isovector (isoscalar) channel if we restrict the number of poles to 3 (4). Furthermore, we have one more free parameter that characterizes the onset of the leading logarithms from perturbative QCD. A difference to the earlier fits MMD ; HMD  is that the constraint from the neutron charge radius has somewhat changed. While in MMD ; HMD  an electron-neutron scattering length of $`b_{ne}=(1.308\pm 0.05)10^3`$ fm was used, the reevaluation of the data lead to $`b_{ne}=(1.33\pm 0.027\pm 0.03)10^3`$ fm for scattering off 208Pb and $`b_{ne}=(1.44\pm 0.033\pm 0.06)10^3`$ fm for 209Bi Kop . We have performed fits with both values but only show results for the Pb value which seems to be favored by the $`G_E^n`$ data.

## 3 Results

Before discussing the results, we must specify the data to which we fit. We use the data basis collected and specified by Friedrich and Walcher in Ref. FW . It consists of a total of 190 data points, about a quarter of which were not available in 1996 (mainly for the neutron). Furthermore, the data basis has also been pruned for inconsistent data and therefore does not make use of some of the data that were included in the fits of MMD ; HMD . For a detailed list of the data taken into consideration, we refer the reader to Ref. FW . In particular, the Jefferson Lab data for $`G_E^p/G_M^p`$ Jones:1999rz ; Gayou:2001qd  are treated as data for $`G_E^p`$ in FW  since $`G_M^p`$ is supposed to be known better for these virtualities.

Fitting these data with the readjusted constraints leads to the following parameters. The isoscalar masses are $`M_\omega =0.782`$ GeV, $`M_\varphi =1.019`$ GeV, $`M_S^{}=1.65`$ GeV, and $`M_{S^{\prime \prime }}=1.68`$ GeV. Note that in contrast to the earlier fits, we can not work with 3 isoscalar poles only. The isovector masses are $`M_\rho ^{}=1.05`$ GeV, $`M_{\rho ^{\prime \prime }}=1.465`$ GeV, and $`M_{\rho ^{\prime \prime \prime }}=1.70`$ GeV. Note that except from $`M_\rho ^{}`$, whose value is fixed from the various constraints, all these masses correspond to physical particles listed in the PDG tables. The corresponding residua for the isoscalar poles are: $`a_1^\omega =0.767`$, $`a_2^\omega =0.318`$, $`a_1^\varphi =0.832`$, $`a_2^\varphi =0.250`$, $`a_1^S^{}=2.09`$, $`a_2^S^{}=3.97`$, $`a_1^{S^{\prime \prime }}=2.04`$, and $`a_2^{S^{\prime \prime }}=3.76`$, where the subscript 1 (2) refers to the vector (tensor) coupling of the corresponding vector meson to the nucleon. Similarly, we have for the three isovector poles: $`a_1^\rho ^{}=0.154`$, $`a_2^\rho ^{}=0.306`$, $`a_1^{\rho ^{\prime \prime }}=1.15`$, $`a_2^{\rho ^{\prime \prime }}=4.42`$, $`a_1^{\rho ^{\prime \prime \prime }}=1.32`$, and $`a_2^{\rho ^{\prime \prime \prime }}=3.18`$. The QCD parameters (for definitions, see MMD ) are $`\gamma =2.148`$, $`\mathrm{\Lambda }^2=18.28`$GeV2 and $`Q_0^2=0.35`$GeV2. The value of $`\mathrm{\Lambda }^2`$ characterizes the onset of the leading logarithmic behavior from perturbative QCD. We note that the value for $`\mathrm{\Lambda }^2`$ chosen here is about a factor of two larger than in the earlier fits MMD ; HMD . This change is due to the strong deviation of the new Jefferson Lab data for $`G_E^p/G_M^p`$ Jones:1999rz ; Gayou:2001qd  from the asymptotic prediction of perturbative QCD.

The errors in the fit parameters quoted above are unknown. Due to the highly nonlinear nature of the problem (remember that we fit all 4 form factors simultaneously) a trustworthy error analysis is a very non-trivial problem. This is complicated by the fact that the errors of some inputs (e.g. the two-pion continuum Hohler:1974ht ) are not known. For this reason, no errors were given in the previous dispersion analyses by Höhler et al. Hohler:1976ax  and Mergell et al. MMD  (and also not in the widely used dispersion theoretical analysis of pion-nucleon scattering by Höhler and coworkers LB ). At present, it is an open problem how to assign a serious theoretical error to such type of analysis.

In the previous dispersion analyses MMD ; HMD , a strong fine tuning of the residues of the excited $`\rho `$ mesons in the isovector channel could be observed. It is well known that the experimental dipole behavior of the form factors requires two narrow structures with opposite signs close to each other. The fits in the previous analyses, however, tended to favor very large residues ($`40`$) for a small overall gain in $`\chi ^2`$. Taken at face value, these residues implied unphysically large couplings for the isovector poles to the nucleon. In order to avoid this problem, we have restricted the magnitude of the vector meson coupling constants in the present analysis. As a consequence, the residues of all poles in the isoscalar and isovector channels are of order one.

In Fig. 1, we show the resulting Sachs form factors by the solid line. The form factors are normalized to the canonical dipole fit, with the exception of the neutron electric form factor, which is not scaled. We can describe most of the data fairly well, with the exception of the Jefferson Lab data on $`G_E^p`$ in the $`Q^2`$ range from 3 to 6 GeV2 (see also the discussion in Refs. Ulfff ; Hammer:2001rh ). The $`\chi ^2`$/datum of this fit is 2.07. For comparison, we exhibit fit 1 of Ref. HMD  by the dash-dotted line. This fit is based on the 1996 data basis and describes the data clearly worse than our new fit. Also shown by the dashed line is the fit with the “phenomenological ansatz” of FW , which has 6 free parameters for each of the 4 form factors.222We have reproduced this fit from Eqs. (4-6) and the numbers given in Table 2 of the published version of Ref. FW . For $`G_M^p`$, we have used the following values for the parameters: $`a_{10}=1.0024`$, $`a_{11}=0.7554`$, $`a_{20}=1.0000a_{10}=0.0024`$, $`a_{21}=14.68`$, $`a_b=0.159`$, $`Q_b=0.325`$, $`\sigma _b=0.231`$ which differ slightly from what is given in Table 2 of the published version of Ref. FW . With the numbers in that table, one can not reproduce the curve shown in Fig. 4 of FW .  It does significantly better for $`G_E^p`$. In Ref. FW , a number of other fits with a “physically motivated ansatz” were also performed. These fits are generally comparable in quality to the one shown in Fig. 1 but give a better description of $`G_M^p`$.

In the fits of Ref. FW , a pronounced structure around $`Q^20.2`$GeV2 was found for $`G_E^n`$ (and could also be isolated in the other form factors as a small effect).

This region of momentum transfers is shown in Fig. 2 in greater detail. The structure was attributed to a contribution from a very long range pion cloud (for a discussion of this point, see also Ref. HMD2 ). Our analysis does not show the pronounced structure around $`Q^20.2`$GeV2. What this result means remains to be seen. Taking the errors of the experimental data into account, the structure is not an unambiguous feature of the data (cf. Fig. 2). Ultimately, the question of whether the structure found in Ref. FW  is real physics or not should be decided by more accurate data in the region $`Q^20.2`$GeV2. A similar plot for the proton charge form factor in the low momentum transfer region is shown in Fig. 3. The data in this range are sufficiently spread to leave room for various interpretations.

The $`\omega NN`$ and $`\varphi NN`$ coupling constants derived from the residua given above are $`g_1^{\omega NN}=21.4`$, $`g_2^{\omega NN}=0.9`$, $`g_1^{\varphi NN}=10.3`$, and $`g_2^{\varphi NN}=3.1`$. The absolute values of these coupling constants are comparable with the results of Refs. MMD ; HMD  but the small tensor couplings $`g_2^{\omega NN}`$ and $`g_2^{\varphi NN}`$ have changed sign. This change of sign indicates that the error in the extracted vector meson coupling constants is of the order of the tensor couplings.

For the electric and magnetic radii of the nucleon, we find the values $`r_E^p=0.848`$ fm, $`r_M^p=0.857`$ fm, and $`r_M^n=0.879`$ fm. As in the previous analysis MMD ; HMD , the proton charge radius is somewhat small compared to the recent precise determinations $`r_E^p=0.880(15)`$ fm Rosenfelder  and $`r_E^p=0.895(18)`$ fm Sick:2003gm  from low-momentum-transfer data and $`r_E^p=0.883(14)`$ fm Melnikov:1999xp  from Lamb shift measurements.

Note that the recent determinations of the proton charge radius by Rosenfelder Rosenfelder  and Sick Sick:2003gm  are based on analyzing low-momentum cross section data directly and include corrections from two-photon exchange. Our analysis relies on form factor data extracted in the one-photon exchange approximation where the two-photon corrections enter in the systematical error. A comparison of our value to the analyses of low-momentum transfer cross sections directly allows one to estimate systematic errors in the extraction of the radii. Of course, one might contemplate a fit to cross section data only. However, for consistency the two-photon corrections have to be applied to the full data basis (as opposed to low $`Q^2`$ data only) and this will certainly take a long time. This issue has to be addressed in the future but requires a large common effort by experimenters and theorists.

Our results for the magnetic radii of the proton and neutron are consistent with the recent values from continued fraction expansions: Kubon et al. Kubon:2001rj  extracted $`r_M^n=0.873(11)`$ fm from precise data for $`G_M^n`$, while Sick Sick:2003gm ; Sickprivate  obtained $`r_M^p=0.855(35)`$ fm from a careful analysis of the world data on elastic electron-proton scattering. To simplify the comparison, we have collected our results and the other recent radius determinations in Table 1.

## 4 Outlook

We have shown that a dispersion-theoretical analysis based on a minimal number of poles can describe most of the current world data on the nucleon electromagnetic form factors for space-like momentum transfer. While the charge radius of the proton is still somewhat small, the magnetic radii of the proton and neutron are in good agreement with the recent determinations using continued fraction expansions Sick:2003gm ; Kubon:2001rj ; Sickprivate .

The spectral functions used here can be improved. Ideally, one would like spectral functions consisting of various continuum contributions with vector mesons emerging naturally as resonances and avoid explicit pole terms altogether. In particular, the isoscalar region around masses of 1 GeV could be supplemented by explicit $`\pi \rho `$ Meissner:1997qt  and $`K\overline{K}`$ HaRM  continuum contributions. The $`\varphi `$ meson would then appear as a finite width resonance in the $`K\overline{K}`$ (and possibly $`\pi \rho `$) continua, thereby eliminating the need to include it explicitly as a pole term. It would also be desirable to replace the poles with masses in the region from 1.5 to 2.0 GeV with an appropriate continuum contribution. However, the analytical continuation of experimental data that is required to obtain the continua becomes more and more difficult the higher one gets in mass, so that the latter goal is not realistic in the near future. Another improvement is to construct a better representation of the perturbative QCD behaviour for large $`Q^2`$. This would allow to include all existing time-like data in the fits, and thus lead to a more consistent description of these fundamental quantities.

Furthermore, the proton charge radius extracted from low-momentum-transfer data and Lamb shift measurements could be included as a further constraint in the analysis. This would allow to check the consistency between the high-$`Q^2`$ data and the low-$`Q^2`$ extractions of the radii. Whether this requires the introduction of additional vector meson poles is an open question. Last but not least, a full error analysis of the extracted radii and coupling constants should be carried out. As mentioned above, this is a nontrivial task because of the highly nonlinear nature of the problem. Work along these lines is under way BHM .

## Acknowledgment

We acknowledge valuable discussions with D. Drechsel, J. Friedrich, I. Sick, and Th. Walcher. Furthermore, we thank J. Friedrich and Th. Walcher for providing us with their data basis and updated fit parameters.