text:
# Neutron Stars with a Stable, Light Supersymmetric Baryon

## 1 Introduction

The interiors of neutron stars offer a unique meeting point between astrophysics on one hand and nuclear and particle physics on the other. The macroscopic properties of neutron stars, such as mass, radius, spin rate and thermal history depend on the microscopic nature of matter at very high densities. Neutrons stars are, in this sense, cosmological laboratories of physics at densities exceeding the nuclear saturation density ($`\rho _00.16`$fm$`^3`$). One aspect of matter at supernuclear densities is the possible appearance of hadronic degrees of freedom beyond nucleons. The principle candidates are hyperons, meson condensates and a deconfined quark phase (see, e.g., Glendenning 1996, Prakash et al. 1997 for reviews).

In this work we wish to illuminate another possibility, the appearance of an $`R`$-baryon. We are motivated by the observation of Spergel & Steinhardt (2000) that the observed properties of dark matter (DM) on scales of order 10-100 kpc are better described by a theory with self-interacting DM than by conventional, weakly interacting DM. Since the required cross-section to mass ratio, $`\sigma _{XX}/m_X1`$b/GeV, is characteristic of hadronic physics at low energies, a light stable neutral hadron would be a welcome feature of physics-beyond-the standard model, and such particles arise naturally in theories with a very light gluino. In some sense, this is the minimal extension of the standard model needed to produce the required DM particle, and therefore merits examination.

A supersymmetric version of the standard model contains a massless gluino ($`\stackrel{~}{g}`$), the spin-1/2 superpartner of the gluon; therefore it is endowed with a corresponding chiral symmetry called $`R`$-invariance. When supersymmetry is broken, the $`R`$-invariance may or may not break. If it does not, the gluino mass arises only from loops sensitive to the Higgs vevs (which do break $`R`$-invariance) and is naturally of order $`100`$$`\text{MeV}\text{c}^2`$. Although considerable effort has been devoted to excluding it, a very light gluino remains a possibility – mainly because of the difficult of making firm, first-principles predictions for non-perturbative phenomena in QCD. The question deserves very serious scrutiny because the existence of a light gluino would resolve several puzzles. These issues and the constraints imposed by experiment are reviewed by Farrar (1998, 2000). See also Farrar & Gabadadze (2000) and references therein.

The lightest supersymmetric “$`R`$-hadrons” are the spin-1/2 $`\stackrel{~}{g}`$$`g`$ bound state ($`R^0`$ or glueballino) and the spin-0 baryon $`uds`$$`\stackrel{~}{g}`$ bound state, the $`\mathrm{S}^0`$ . The $`\mathrm{S}^0`$ is thus a boson with baryon number unity and a strangeness charge $`1`$. The strong color hyperfine interaction unique to the $`uds`$$`\stackrel{~}{g}`$ in a flavor singlet state makes the mass of the $`\mathrm{S}^0`$ much less than that of all other $`R`$-baryons (Farrar, 1984; Buccella & Farrar, 1985), and probably less than 2 GeV$`c^2`$, which would make it stable through conservation of $`R`$parity (Farrar, 1996).

In the very early universe the $`\mathrm{S}^0`$ remains in equilibrium with nucleons via processes such as $`R^0N\mathrm{S}^0K`$, down to $`T35`$ MeV (Farrar, Spergel & Steinhardt 2000; see also in Wandelt et al. 2000). Hence for a suitable $`\mathrm{S}^0`$ mass, slightly less than the proton mass, standard freezeout estimates lead to a density ratio $`\mathrm{\Omega }_{\mathrm{S}^0}/\mathrm{\Omega }_B610`$, as required if the $`\mathrm{S}^0`$ is to be the dark matter. To posit a new baryon lighter than a nucleon is quite shocking! But it is known that the color hyperfine force is very strong in this channel, and there are innumerable examples of the inadequacy of phenomenological models of non-perturbative QCD phenomena. The possibility of such a light $`\mathrm{S}^0`$ should not be dismissed on theoretical grounds, and we must seek observational constraints. One possible astrophysical constraint is whether this new baryon can form in neutron stars, and if so, would it imply unacceptable neutron star properties?

Nucleons cannot convert directly into the $`\mathrm{S}^0`$ , even if the latter is less massive. Conservation of $`R`$parity requires that $`R`$-hadrons are produced in pairs, forbidding direct conversion of individual nucleons into $`\mathrm{S}^0`$ particles. The weak neutral current process $`nn\mathrm{S}^0\mathrm{S}^0`$ requires simultaneous creation of two units of strangeness ($`\mathrm{\Delta }S=2`$). Strangeness-changing neutral currents ($`\mathrm{\Delta }S=\pm 1,\mathrm{\Delta }Q=0`$) are highly suppressed in the standard model due to the GIM mechanism. Thus, $`nn\mathrm{S}^0\mathrm{S}^0`$ is doubly GIM suppressed; crude estimates indicate that its rate is of order the inverse Hubble time (Farrar et al., 2000). However processes such as $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{S}^0\mathrm{S}^0`$ ($`\mathrm{\Delta }S=0`$) and even $`\mathrm{\Sigma }^{}p\mathrm{S}^0\mathrm{S}^0`$ (charged current, $`\mathrm{\Delta }S=1`$) are not GIM suppressed, so that if the conditions in interiors of neutron stars allow for hyperons to appear (as is predicted by many works on this issue, see, e.g., Glendenning 1985, Weber & Weigel 1989, Knorren, Prakash & Ellis 1995, Schaffner & Mishustin 1996, Balberg & Gal 1997, Huber et al. 1998, Baldo, Burgio & Schulze 2000), these hyperons could serve as a doorway for $`\mathrm{S}^0`$ formation, which might affect observable properties of neutron stars.

Our goal in this letter is to explore whether a light $`\mathrm{S}^0`$ can form in neutron stars and assess whether the resulting global properties of these stars are compatible with observations. In § 2 we extend an equation of state for high density matter with hyperons to include the $`\mathrm{S}^0`$ . The equilibrium composition with the $`\mathrm{S}^0`$ is found in § 3, followed by a discussion of the corresponding properties of neutron stars in § 4.

## 2 $`\mathrm{S}^0`$ interactions in high density matter

In order to include $`\mathrm{S}^0`$ particles in an equation of state of high density matter, we must specify their interactions with other baryons and among themselves. Here we do so in the framework of an effective equation of state (Balberg & Gal 1997, Balberg, Lichtenstadt & Cook 1999), where the potential felt by a single, isospin singlet particle of species $`x`$ placed in matter composed of species $`y`$ at a baryon number density $`\rho _y`$ is

$$V_x(\rho _y)=a_{xy}\rho _y+c_{xy}\rho _y^\gamma +w\rho _y^\theta .$$
(1)
In equation (1) $`a`$ is negative and represents a long range, attractive component in the interaction, $`c`$ is positive, corresponding to an intermediate range repulsive component, and $`w`$ is also positive, representing a universal repulsive short range interaction, presumably arising from the quark structure of the baryons. The powers $`\gamma `$ and $`\theta `$ must be greater than unity, so that repulsion dominates at high densities.

For nucleons and hyperons, the values of the coefficients can be constrained by experimental data. Since no such data exists for the $`\mathrm{S}^0`$ , we resort to phenomenological arguments to estimate the magnitude of these terms. The attractive component is generally associated with two-pion or two-kaon exchanges. This interaction is strong because the intermediate state in the $`s`$-channel can be the same as the initial state, hence the amplitude has a small energy denominator. However the flavor singlet $`\mathrm{S}^0`$ cannot exchange $`\pi `$ or $`K`$ without becoming a much more massive flavor octet $`R`$-baryon (Buccella & Farrar, 1985). Thus the long-range attractive interaction is strongly suppressed for the $`\mathrm{S}^0`$ . A more empirical way to arrive at the same conclusion is to note that if there were a significant long-range attractive force the $`\mathrm{S}^0`$ would bind to nuclei (and to itself), and nuclei with exotic mass values have not been observed. We therefore assume that the attractive component is negligible for the $`\mathrm{S}^0`$ $``$related interactions. The intermediate repulsive component may be associated with the exchange of vector mesons, such as the $`\rho `$ and $`\eta `$. Only the flavor-singlet component of these vector-mesons can couple to the $`\mathrm{S}^0`$ , so we roughly assume that the magnitude of this term (the value of $`c`$ in Eq. ) is half of its value for the $`\mathrm{\Lambda }`$. The universal short-range interaction for the $`\mathrm{S}^0`$ is taken as identical to that of the baryons. Here we use EoS 1 of Balberg et al. (1999), where $`\gamma =\frac{4}{3}`$ and $`\theta =\frac{5}{3}`$, so that for the $`\mathrm{S}^0`$ $`a=0`$ and $`c=329.5\text{MeV}\text{fm}^4`$ for all its interactions, and $`w=223.6\text{MeV}\text{fm}^5`$. Below we refer to this choice of parameters along with an $`\mathrm{S}^0`$ mass of $`m_{\mathrm{S}_0}=900`$$`\text{MeV}\text{c}^2`$as our nominal set. Most of our conclusions prove insensitive to these choices, as will be seen below.

## 3 Composition of neutron star matter with the $`\mathrm{S}^0`$

The composition of neutron star matter is determined by chemical equilibrium for all weak-interaction processes. Specifically, these are $`B_1B_2\mathrm{}\nu _{\mathrm{}}`$, where $`B_1`$ and $`B_2`$ are baryons, $`\mathrm{}`$ is a lepton and $`\nu _{\mathrm{}}`$ is its corresponding neutrino (or antineutrino). Neutron star matter is well approximated as cold ($`T=0`$) and neutrino free ($`\mu _\nu =0`$) so that all equilibrium conditions are summarized concisely in the form:

$$\mu _{B_i}=\mu _nq_{B_i}\mu _e,$$
(2)
where $`\mu _{B_i}`$ is the chemical potential of the baryon species $`B_i`$, $`\mu _n`$ and $`\mu _e`$ are the neutron and electron chemical potentials, respectively, and $`q_{B_i}`$ is the electric charge of the baryon $`B_i`$. The equilibrium composition at any baryon number density, $`\rho _B`$, is then determined by two additional conditions: $`\mathrm{\Sigma }_ix_{B_i}=1`$ and charge neutrality, $`\mathrm{\Sigma }_ix_{B_i}q_{B_i}+\mathrm{\Sigma }_ix_{\mathrm{}}q_{\mathrm{}}=0`$, where $`x_{B_i}`$ and $`x_{\mathrm{}}`$ denote the fraction per baryon for each baryon and lepton species, respectively.

Equation (2) does not fully apply to the case of the $`\mathrm{S}^0`$ , since the process $`nn\mathrm{S}^0\mathrm{S}^0`$ is suppressed, even if energetically favorable. The $`\mathrm{S}^0`$ must maintain a direct equilibrium only with hyperons, i.e., $`2\mu _{\mathrm{S}^0}=2\mu _{Y^0}`$, and $`2\mu _{\mathrm{S}^0}=\mu _Y^{}+\mu _p`$, but $`2\mu _{\mathrm{S}^0}=2\mu _n`$ is not required (where $`Y^0`$ and $`Y^{}`$ denote a neutral and a negatively charged hyperon, respectively). The significance of this distinction is that a combination of $`\mathrm{S}^0`$ and nucleons may coexist even while $`\mu _{\mathrm{S}^0}<\mu _n`$, as long as no hyperons are present. On the other hand, a finite hyperon fraction will indirectly impose that $`\mu _{\mathrm{S}^0}=\mu _n`$, so then equation (2) applies to all species.

We emphasize that even if $`\mu _{\mathrm{S}^0}\mu _n`$, nuclear matter cannot “burn” exothermally into $`\mathrm{S}^0`$ matter. Unlike the case of “strange-quark-stars” (Witten, 1984; Alcock, Farhi & Olinto, 1986), where a deconfining front converts an entire neutron star into the (hypothetically) lower energy strange quark matter if a seed of of strange quark matter is somehow formed (Olinto, 1987), there exists no equivalent stimulant to initiate conversion of nuclear-matter into $`\mathrm{S}^0`$ -matter, regardless of density or composition.

In Fig. 3a we show for reference the equilibrium composition of neutron star matter as a function of baryon density, using the “standard” equilibrium composition of matter with nucleons and hyperons but no $`\mathrm{S}^0`$ (hereafter, NH matter) based on EoS 1 of Balberg et al. (1999), assuming that the equation of state is valid up to densities as large as $`10\rho _0`$. The results, with hyperon formation beginning at about $`2\rho _0`$, are typical of matter with hyperons (Glendenning, 1985; Weber & Weigel, 1989; Knorren et al., 1995; Schaffner & Mishustin, 1996; Balberg & Gal, 1997; Huber at al., 1998; Baldo et al., 2000).

Figure 3b depicts the equilibrium composition when the $`\mathrm{S}^0`$ is included according to the assumptions of § 2. Due to the low mass of the $`\mathrm{S}^0`$ , its chemical potential is lower than that of all hyperons, even when the fractions of the latter are infinitesimally small. As a result, all hyperons that form convert to $`\mathrm{S}^0`$ particles, and the matter consists only of nucleons and $`\mathrm{S}^0`$ . The composition is determined by the nucleons being in equilibrium with an exactly zero fraction of hyperons 111In reality, some additional nucleons would convert to hyperons through quantum and thermal ($`T0.01`$MeV in neutron star cores) fluctuations, and these hyperons would decay to $`\mathrm{S}^0`$ ’s. The actual equilibrium composition would be one where the chemical potential of the nucleons is slightly below that of zero-fraction hyperons at every density (in other words, fluctuations would deplete a thin outer layer in phase space of nucleons). Assuming that the probability of fluctuations decreases exponentially with the energy gap, the quantitative effect on the results should be small.. Even though absent from the matter, the hyperons do play a role in determining the amount of $`\mathrm{S}^0`$ particles which are created.

Since the $`\mathrm{S}^0`$ is never in chemical equilibrium with the nucleons (or the zero fraction hyperons), the equilibrium composition would be identical for any $`m_{\mathrm{S}_0}<900`$$`\text{MeV}\text{c}^2`$. In fact, for the nominal interactions, the same composition is also found for $`m_{\mathrm{S}_0}m_p`$. A qualitative change is found when $`m_{\mathrm{S}_0}`$ is larger: in Fig. 3c we show the equilibrium composition for nominal $`\mathrm{S}^0`$ interactions, but with $`m_{\mathrm{S}^0}=975`$$`\text{MeV}\text{c}^2`$. At low and intermediate densities the composition is identical to that of the $`m_{\mathrm{S}^0}=900`$$`\text{MeV}\text{c}^2`$case, but beginning at a density of 0.9 fm$`^3`$the larger mass of the $`\mathrm{S}^0`$ makes it favorable for a finite fraction of hyperons to coexist in the matter. The general equilibrium condition (Eq. ) is now satisfied for all particle species. Owing to the $`\mathrm{S}^0`$ having only repulsive interactions, the hyperons achieve chemical equilibrium with it even though it has a lower mass and no kinetic energy. Moreover, for larger densities, it becomes increasingly favorable to reduce the $`\mathrm{S}^0`$ fraction as the equilibrium composition gradually approaches that of NH matter. With similar calculations we find that for nominal interactions, if $`m_{\mathrm{S}^0}1030`$$`\text{MeV}\text{c}^2`$the $`\mathrm{S}^0`$ appears only after some hyperons already exist in the matter, and if $`m_{\mathrm{S}^0}1060`$$`\text{MeV}\text{c}^2`$$`\mathrm{S}^0`$ formation is suppressed for the entire range of densities.

## 4 Equation of State and Mass Sequences

We now turn to examine the equation of state (pressure vs. mass-energy density, hereafter EoS) of high density matter with the $`\mathrm{S}^0`$ . These are shown in Fig. 4, including the models used for Fig. 3 and also a model labeled N for nuclear matter,      Equation of state for matter of the compositions shown in Fig. 1. Matter with the $`\mathrm{S}^0`$ : nominal values (see text) - thick solid line; nominal interactions and $`m_{\mathrm{S}^0}=975`$$`\text{MeV}\text{c}^2`$\- dashed line; nominal mass and increased repulsion - dot-dashed line. Also shown are the equations of state for nucleon and hyperon matter - N+H, and nuclear matter - N.  where hyperon (and therefore also $`\mathrm{S}^0`$ ) formation is ignored.

The hierarchy of EoSs in Fig. 4 is mostly determined by the number of available species in each model. Nuclear matter naturally has the stiffest EoS, followed by matter with the $`\mathrm{S}^0`$ . Even though the $`\mathrm{S}^0`$ , being a boson, does not exert any Fermi pressure, the fact that four species of hyperons are suppressed due to its appearance makes such matter stiffer than NH matter, which is clearly the softest. An exception occurs just above the threshold density for hyperon formation, where NH matter includes only a single hyperon species (the $`\mathrm{\Sigma }^{}`$), so that the contribution of its Fermi pressure does make NH matter stiffer.

Although the mass of the $`\mathrm{S}^0`$ does not affect its contribution to the pressure, it does contribute to the mass energy-density. Hence a lower $`m_{\mathrm{S}^0}`$ yields a stiffer EoS, while a larger value yields a softer one. For $`m_{\mathrm{S}^0}=975`$$`\text{MeV}\text{c}^2`$such softening is enhanced at a density of $`1.8\times 10^{15}`$$`\text{gm}\text{cm}^3`$, when hyperons begin to accumulate; as the $`\mathrm{S}^0`$ fraction is depleted this EoS gradually approaches that of NH matter. It is noteworthy that increasing the repulsive power of the $`\mathrm{S}^0`$ does not lead to significant stiffening of the EoS. We demonstrate this with another model shown in Fig. 4 when $`m_{\mathrm{S}^0}=900`$$`\text{MeV}\text{c}^2`$but $`c=494.5\text{MeV}\text{fm}^4`$ (50% larger than the nominal value). Magnified repulsion leads to an increase of the $`\mathrm{S}^0`$ chemical potential, which in turn enhances hyperon formation. The resulting composition of this model is quite similar to that of Fig. 3c, and the EoS, while slightly stiffer than nominal until hyperons appear, is then softened. We emphasize that this is a manifestation of the ability of hyperons to exert “pressure control” (Balberg & Gal, 1997; Balberg et al., 1999), since the extent of their formation tracks the dependence of particle interactions on density.

Given an EoS, the structure of a spherical, non-rotating neutron star is found by numerically integrating the Tolman-Oppenheimer-Volkoff equations of general relativistic hydrostatic equilibrium (see, e.g., Misner, Thorne & Wheeler 1973). The resulting mass sequences (stellar mass as a function of the assumed central density) for the EoSs of Fig. 4 are given in Fig. 4. All sequences are identical up to about $`1.1M_{}`$, where the mass sequences diverge according to the hierarchy of EoSs. The maximum mass found for matter with the $`\mathrm{S}^0`$ is $`1.81M_{}`$ for the nominal case, $`1.68M_{}`$ for $`m_{\mathrm{S}^0}=975`$$`\text{MeV}\text{c}^2`$,      Nonrotating neutron star sequences for the five equations of state shown in Fig. 2.  and $`1.78M_{}`$ for $`m_{\mathrm{S}^0}=900`$$`\text{MeV}\text{c}^2`$but enhanced $`\mathrm{S}^0`$ repulsive interactions. In both latter cases, hyperons coexist with the $`\mathrm{S}^0`$ only in the largest mass stars, those of mass closer than $`0.1M_{}`$ to the maximum mass in each case. We note that for EoS 1, the maximum masses for nuclear matter and NH matter are $`2.18`$ and $`1.57M_{}`$ respectively.

At present, the primary observational constraint on high density equations of state is that the maximum mass must be at least $`1.44M_{}`$, the inferred mass of PSR 1913+16. All EoSs discussed above, with and without the $`\mathrm{S}^0`$ , are acceptable. This limit does allow to rule out very soft EoSs, such as would arise if the $`\mathrm{S}^0`$ had a low mass and an attractive component in its interactions. Such interactions would make a pair of $`\mathrm{S}^0`$ similar to the hypothetical H dibaryon - a six quark flavor singlet, boson (Jaffe 1977a, 1977b) - whose interactions are supposed to be similar to those of two $`\mathrm{\Lambda }`$ hyperons combined. The EoS for matter with the H would indeed be too soft to sustain a $`1.44M_{}`$ star, unless its mass is significantly above $`2m_\mathrm{\Lambda }`$ (Glendenning & Schaffner-Bielich, 1998). The expected strongly repulsive interactions of the $`\mathrm{S}^0`$ are thus crucial for allowing its existence in neutron stars.

## 5 Conclusions and Discussion

Our primary conclusion is that the existence of a light ($`m_{S^0}m_p`$) $`\mathrm{S}^0`$ particle is consistent with the fundamental observed properties of neutron stars, namely, their inferred masses. This conclusion applies for a wide range of models for high density matter with the $`\mathrm{S}^0`$ . In particular, the maximum mass of a static neutron star under nominal assumptions regarding the $`\mathrm{S}^0`$ is found to be about $`1.8M_{}`$. We note that since it suppresses hyperon accumulation, the appearance of the $`\mathrm{S}^0`$ yields a stiffer equation of state and higher maximum mass than obtained with nucleons and hyperons without $`\mathrm{S}^0`$ . Therefore, identification of a relatively high mass neutron star would be easier to reconcile with an $`\mathrm{S}^0`$ component than in its absence (where hyperons accumulate undisturbed). Some indication for the existence of such massive neutron stars is implied by recent analyses, although uncertainties are still significant (van Kerkwijk, 2000). We emphasize that current data on neutron star masses is not restrictive enough to decisively constrain the properties of high density matter.

The presence of $`\mathrm{S}^0`$ particles could have a significant and welcome effect on the cooling rate of neutron stars, if they do indeed replace hyperons in the cores. Hyperons can participate in the direct Urca cooling process ($`B_1B_2e\nu `$) while the $`\mathrm{S}^0`$ cannot. The threshold fraction required for hyperons to drive direct Urca cooling processes is of the order of 0.01 (Prakash et al., 1992), so neutron stars with hyperons are susceptible to rapid cooling, which seems inconsistent with observations (Schaab et al., 1996). Hyperon direct Urca cooling can be partially suppressed if the $`\mathrm{\Lambda }`$ hyperons pair to a $`{}_{}{}^{1}\mathrm{S}_{0}^{}`$ superfluid state (Schaab, Balberg & Schaffner-Bielich, 1998), but the extent of suppression is sensitive to the details of the equilibrium composition and gap energy (Page et al., 2000). On the other hand, $`\mathrm{S}^0`$ formation would prevent the hyperon direct Urca, or restrict it only to the highest mass stars, so that very moderate cooling rates are attainable. Future measurements of surface thermal emission from young neutron stars with the CHANDRA and XMM satellites may lead to stronger constraints on cooling rates; if so this will yield valuable insight into the internal structure of high density matter.

The research of GRF is supported in part by NSF-PHY-99-96173. The idea that the $`\mathrm{S}^0`$ could be a dark matter candidate, and constraints on such a scenario, were developed by GRF in collaboration with D. Spergel and P. Steinhardt; that work will be reported in greater depth elsewhere (Farrar et al., 2000); a brief description can be found in (Wandelt et al., 2000). SB and TP thank NYU for hospitality which facilitated this research.


format: text

text:
# Anatomy of giant spiral galaxies

## 1. Introduction

The progress made in last years in numerical modeling allows to study the internal structure of galaxies with great detail. This is of particular interest for stellar systems composed by several components, such as spiral galaxies. In the recent years, we produced a self-consistent, 3-D model of disk galaxy made by an exponential disk and a $`r^{1/4}`$ bulge (Pignatelli & Galletta, 1999).

The model fit simultaneously the photometric data (surface brightness, ellipticities and PAs of isophote major axis) and the kinematic data (rotation curves and velocity dispersion curves along various PAs). In such a way, the structure of the luminous matter derived from the photometric analysis is compared with the total matter deduced from the dynamical data and the presence of dark matter may be studied.

To test the model, we selected five galaxies with morphological types spanning from S0/a to Sc and absolute magnitudes from -20.6 to -22.5. They have been selected from HI observations looking for objects with a line width $`W_{20}>350`$ km/sec.

For all these galaxies we obtained images at B and I photometric wave bands using the ESO telescopes at La Silla. The used telescopes were the Danish 1.5m, the NTT+EMMI and the ESO-MPI 2.2m+EFOSC2. The spectra have been collected with the spectrographs Boller& Chivens and EFOSC2 of the 2.2m ESO-MPI telescope.

We used the model to evaluate the intrinsic properties of these galaxies, such as the disk/bulge mass ratio and the scale length of the galaxy components.

## 2. The results from the model

The galaxies can be described by superposition of different components. For each component, we separately assume that: (a) the density distribution is oblate, without strongly triaxial structures; (b) the isodensity surfaces are similar concentric spheroids; (c) the surface density profile follows a simple $`R^{1/4}`$ or exponential law; (d) the velocity distribution is locally Gaussian; (e) the velocity dispersion is isotropic $`\sigma _r=\sigma _\theta =\sigma _z`$;

The model has $`4n+1`$ free parameters, where $`n`$ is the number of adopted components: namely the luminosity $`L_{tot}`$, the effective radius $`R_e`$, the mass-luminosity ratio $`M/L`$ and the flattening $`b/a`$ of each component plus the inclination angle $`i`$ of the galaxy.

Photometry can be used to constrain all these parameters except the M/L ratio, which must be derived by kinematics. In Fig. 2 we show an example of the application of the model to the rotation curve of NGC 6925: in this case, we found that the stellar velocity can be fitted by a model without dark matter with a $`M/L=3.9`$ solar units and a total mass of $`2.710^{11}M_{}`$.

Since these galaxies have high value for the HI line width $`W_{20}`$, it has been suggested (Buson et al., 1991) that they could have in general an unusually high content of dark matter in the inner regions or, perhaps, an unusual stellar population. None of this hypothesis is confirmed by our analysis.

## References

Buson, L. M., Galletta, G., Saglia, R. P., Zeilinger, W. W., 1991, Msngr 63, 50

Pignatelli, E., & Galletta, G., 1999, A&A 349, 369


text:
# Relaxation of non-order parameter field in directed Ising systems

## Abstract

We investigate the effect of initial conditions on the dynamic exponents of the interacting monomer-monomer model with infinitely many absorbing states in one dimension. This model exhibits a directed Ising (DI) type transition from an active phase into an absorbing phase. In case of the directed percolation universality class, it has been reported that the non-order parameter as well as the order parameter exhibits critical fluctuations, relaxing algebraically to its natural value with the same scaling exponents. We numerically confirm that this is also valid for the DI universality class. We also observe continuously varying dynamic exponents with a linear dependence on the non-order parameter initial density.

Various kinds of nonequilibrium lattice models exhibiting absorbing phase transitions have been studied extensively during last few decades. Two distinct types of absorbing phase transitions have been identified in one dimension: the directed percolation (DP) and directed Ising (DI) universality class. Most models have been found to belong to the DP class, which involves typically a single absorbing state or multiple absorbing states without any symmetry. DI-type critical behavior appears in models with two equivalent absorbing states or two equivalent group of absorbing states .

Recently, the question has been addressed whether one can construct initial states that affect the entire temporal evolution of nonequilibrium systems including asymptotic dynamics. Typical examples are the systems that display an absorbing phase transition from an active phase into an absorbing phase with infinitely many absorbing (IMA) states. As in the ordinary order-disorder phase transition models, the order parameter takes a nonzero value in the active phase and vanishes in the absorbing phase. The IMA states can be characterized by a non-order parameter. In the steady state, the non-order parameter approaches algebraically to a natural value at criticality.

Nonuniversal dynamic properties have been reported in various kinds of DP-type nonequilibrium lattice models with IMA states . The order-parameter dynamic exponents vary continuously with the initial conditions characterized by the non-order parameter density. Moreover, at criticality, the non-order parameter density also exhibits critical fluctuations, relaxing algebraically to its natural value with the same dynamic exponents . The critical relaxation of the non-order parameter field in DP-type systems with IMA states was recently examined from a field theoretical (Langevin equation) approach. The evolution equation for the order parameter is non-Markovian and includes a temporal memory term due to the non-order parameter field. This non-Markovian term is known to be responsible for the nonuniversal dynamic exponents.

In this paper, we introduce an interacting monomer-monomer model with IMA states (IMA-IMA model) that belongs to the DI universality class and investigate the effects of the non-order parameter field in one dimension. Using numerical simulations, we show that the non-order parameter plays the same role in DI systems as well as DP systems.

The IMA-IMM model is an interacting monomer-monomer model with two different species of monomers, namely $`A`$ and $`B`$. Monomers $`A`$ and $`B`$ are selected with probability $`p`$ and $`1p`$ respectively. Monomers can adsorb at a randomly selected vacant site with unit probability for $`A`$ and $`1r_B`$ for $`B`$ respectively. An adsorption attempt of a monomer is rejected when both sites adjacent to a selected vacant site are occupied or at least one adjacent site is occupied with a monomer of the same species. A nearest neighbor $`AB`$ pair reacts and desorbs immediately from the lattice.

$`A+V`$ $`\stackrel{p}{}`$ $`A_s,`$ (1)
$`B+V`$ $`\stackrel{(1p)(1r_b)}{}`$ $`B_s,`$ (2)
$`A_s+B_s`$ $``$ $`\mathrm{}.`$ (3)
Here the subscript $`s`$ denotes adsorbed particles, $`V`$ denotes a vacant site where the adsorption attempt is allowed.

This model has infinitely many absorbing states. Any configuration without a nearest neighbor pair of vacant sites is one of the absorbing states. The number of absorbing states diverges exponentially with system size $`L`$ as $`2^{L/2}`$. The absorbing states can be divided into two equivalent groups; odd-site occupied and even-site occupied groups. Clearly, these two groups of absorbing states have one-to-one correspondence and the Ising ($`Z_2`$) symmetry in between. So we expect that this model exhibits the DI-type absorbing phase transition.

We locate the critical line in the $`r_bp`$ phase diagram by dynamic Monte Carlo simulations . For small values of $`r_b`$, the system is always absorbing. As $`r_b`$ increases, a window of the active phase appears and divides the absorbing phase by $`A`$-dominated and $`B`$-dominated absorbing phases. For example, as $`p`$ increases along the $`r_b=0.9`$ line, the system undergoes two continuous phase transitions from the $`B`$-dominated absorbing phase into the active phase at $`p_{c1}=0.100(3)`$ and finally into the $`A`$-dominated absorbing phase at $`p_{c2}=0.5058(7)`$. We numerically confirm that both transitions are of the DI type as expected .

The order parameter $`\rho `$ is the number density of nearest neighbor pairs of vacant sites, while we define the $`A`$ particle density as the non-order parameter $`\rho _A`$. First, we perform static Monte Carlo simulations to measure the natural density of the non-order parameter at both criticalities. We measure the $`A`$ particle density $`\rho _A(t,L)`$, averaged over $`2\times 10^35\times 10^4`$ survived samples for system size $`L=2^62^{10}`$. $`\rho _A(t,L)`$ relaxes to the natural density $`\rho _A^{nat}\rho _A(\mathrm{},\mathrm{})`$ in the thermodynamic limit. At $`p=p_{c1}`$ and $`p=p_{c2}`$, we estimate $`\rho _A^{nat}=0.101(1)`$ and $`0.458(2)`$, respectively. The distribution function for the $`A`$ particle density appears to be Gaussian, so the mean density is identical to the most probable density. We measure the most probable density in the long time limit, which turns out to be consistent with the above mean natural density within statistical errors.

Dynamic properties for the non-order parameter $`\rho _A(t,L)`$ can be extracted by studying its temporal deviation from its steady-state value (natural density) as

$$\mathrm{\Delta }\rho _A(t,L)|\rho _A(t,L)\rho _A^{nat}|.$$
(4)
As in the DP case, we assume that $`\mathrm{\Delta }\rho _A(t,L)`$ follows the same scaling behavior as the order parameter. Using the finite-scaling theory  in the steady state, $`\mathrm{\Delta }\rho _A(t,L)`$ scales at criticality as

$$\mathrm{\Delta }\rho _A(\mathrm{},L)L^{\beta /\nu _{}}.$$
(5)
One can also expect the critical short time behavior as

$$\mathrm{\Delta }\rho _A(t,\mathrm{})t^{\beta /\nu _{}},$$
(6)
and the characteristic time $`\tau _A(L)`$ scales as

$$\tau _A(L)L^{\nu _{}/\nu _{}}.$$
(7)
From static simulations, we estimate the scaling exponents at both criticalities. At $`p=p_{c1},p_{c2}`$, we estimate $`\beta /\nu _{}=0.46(2),0.46(3)`$ $`\beta /\nu _{}=0.26(2),0.25(2)`$, and $`\nu _{}/\nu _{}=1.75(5),1.90(10)`$ (see Fig. 1). As expected, these estimations involve rather large statistical and systematic errors, especially due to inaccuracy of the natural density values. However, these values agree reasonably well with the DI values , which confirm our assumption that the non-order parameter exhibits the same type of critical fluctuations as the order parameter.

Dynamic exponents for the order parameter in DP systems with IMA states are known to depend on initial conditions characterized by the non-order parameter. To investigate this nonuniversal dynamic properties in DI systems, we perform dynamic Monte Carlo simulations with various initial conditions. We start with a pair of nearest neighbor vacant sites in the absorbing background which is controlled by the $`A`$ particle density $`\rho _A^0`$.

We measure the survival probability $`P(t)`$(the probability that the system is still active at time $`t`$) and the mean number of pair of vacant sites (order parameter) $`N(t)`$ averaged over all samples. At criticality, these quantities scale algebraically in long time limit as

$$P(t)t^\delta ,N(t)t^\eta .$$
(8)
The dynamic exponents are in general functions of initial non-order parameter density; $`\delta =\delta (\rho _A^0)`$ and $`\eta =\eta (\rho _A^0)`$. With initial configurations of the natural non-order parameter density, the exponents take the ordinary DI values; $`\delta (\rho _A^{nat})=\delta _{DI}0.285`$ and $`\eta (\rho _A^{nat})=\eta _{DI}0.00`$.

At $`p=p_{c1}`$ with $`r_b=0.9`$, we estimate the values of $`\delta `$ and $`\eta `$ as $`\rho _A^0`$ varies from 0 to 0.35 (see Table I). In Fig. 2, we plot the exponent shifts from the DI values ($`\delta \delta _{DI}`$ and $`\eta \eta _{DI}`$) versus $`\rho _A^0\rho _A^{nat}`$. It shows a linear dependence of the exponent shifts on the deviation of the non-order parameter from the natural density in initial configurations. This linear dependence has been also seen in DP systems.

In summary, we investigated the IMA-IMM model in one dimension, which show the DI-type continuous phase transition from an active phase into an absorbing phase consisting of infinitely many absorbing states. We show that the non-order parameter exhibits critical fluctuations identical to the order parameter. Dynamic exponents $`\delta `$ and $`\eta `$ depends linearly on the initial non-order parameter density and coincide with the ordinary DI values only at the natural density.

This work was supported by an Inha University research grant.


text:
# Bunch Length Measurements at the ATF Damping Ring in April 2000

## 1 Introduction

We want to accurately know the energy spread and bunch length dependence on current in the ATF damping ring. One reason is to know the strength of the impedance: From the energy spread measurements we know whether or not we are above the threshold to the microwave instability, and from the energy spread and bunch length measurements we find out the extent of potential-well bunch lengthening (PWBL). Another reason for these measurements is to help in our understanding of the intra-beam scattering (IBS) effect in the ATF. The ATF as it is now, running below design energy and with the wigglers turned off, is strongly affected by IBS. To check for consistency with IBS theory of, for example, the measured vertical beam size, we need to know all dimensions of the beam, including the longitudinal one. But beyond this practical reason for studying IBS, IBS is currently a hot research topic at many accelerators around the world (see e.g. Ref. ), and the effect in actual machines is not well understood. Typically, when comparing theory with measurements fudge factors are needed to get agreement (see e.g. Ref. ). With its strong IBS effect, the ATF is an ideal machine for studying IBS, and an indispensable ingredient for this study is a knowledge of the longitudinal phase space of the beam.

The results of earlier bunch lengthening measurements in the ATF can be found in Refs. -. Measurements of current dependent effects, especially bunch length measurements using a streak camera, can be difficult to perform accurately. For example, space charge in the camera itself can lead to systematic errors in the measurement results. It is important the results be accurate and reproducible. In the measurements of both December 1998 and December 1999, by using light filters, the authors first checked that space charge in the streak camera was not significant. And then the Dec 99 authors show that their results agree with those Dec 98, i.e. on the dates of the two measurements the results were reproducible.

Since IBS is so strong in the ATF, in the Dec 99 measurements an attempt was made to estimate the impedance effect using the following method: First, from the form of the energy spread vs. current measurements it was concluded that the threshold to the microwave instability was beyond 2 mA. Then, by dividing the bunch length vs. current curve by the energy spread vs. current curve the effect of IBS was divided out, and PWBL was approximated. The assumption is that PWBL can be treated as a perturbation on top of IBS. The result was that this component of bunch lengthening was found to grow by 7-15% (depending on the rf voltage) between the currents of .5 mA and 2 mA, about a factor of 3 less than the total bunch length growth. The conclusion was that the inductive component of the impedance was small, in fact much smaller than had been concluded earlier in Ref. .

Electron machines generally run in a parameter regime where IBS is an insignificant effect, and impedance measurements and calculations have also normally been performed for machines where IBS is unimportant. To simplify the interpretation of the impedance from bunch length measurements, in April 2000 the energy spread and bunch length measurements of Dec 99 were repeated, but now with the beam on a linear (difference) coupling resonance, where the horizontal and vertical emittances were approximately equal. For this case the effect of IBS was expected to be very small. An energy spread vs. current measurement under such conditions will also allow us to more clearly see whether we reach the threshold to the microwave instability. As part of the April data taking we, in addition, repeated the earlier off-coupling measurements, in order to check the reproducibility of the earlier results. In this report we present and analyze this recent set of data, and compare it with the results of the earlier measurements, particularly those of Dec 99.

The measurements and analysis of data in this report follow essentially the same procedure as was used in Ref. . In the present report we will try to be relatively brief. The comparison of our results with IBS theory will be given in a following report. For more details about the measurement and analysis techniques presented in this report, the reader should consult Ref. .

## 2 Energy Spread Measurements

Energy spread vs. current measurements were performed on March 21-22, 2000, with the beam both on and off the coupling resonance, for rf voltages of $`V_c=150`$ kV and 300 kV. As usual the beam width is measured after extraction on a thin screen in a dispersive region. The data is fit to a Gaussian, yielding the rms energy spread $`\sigma _\delta `$ (see Fig. 1). (For these measurements the energy spread was not calibrated, but from earlier experience we expect the zero value results to be near $`5.5\times 10^4`$, which allows us to estimate the scaling factor.) In the figure the points in the plots are fit to $`y=5.5(1+ax^b)`$.

We note from Fig. 1 that, as expected, on the coupling resonance the energy spread is almost independent of current, only increasing by 3-4% by 2.5 mA. With the vertical beam size enlarged, the effect of intra-beam scattering becomes small. Off the coupling resonance the change by 2.5 mA is 35-40%. We also note in (a) and (b) that there is no evidence of the threshold to the microwave instability, whose signature would be a kink in the data. If there is no microwave instability on resonance, then almost certainly there is no microwave instability off resonance, since in the latter case the longitudinal phase space volume is increased, which tends to stabilize the beam. Finally, note that the 150 kV and 300 kV on-resonance results are statistically identical, as would be expected in the case of no microwave instability and negligible intra-beam scattering. Note also that for the off-resonance results the energy spread is slightly larger for the higher voltage. This would be expected, since a higher voltage means a smaller bunch length, which increases intra-beam scattering, which, in turn, more increases the energy spread.

The measurements were repeated for the nominal (off-resonance, $`V_c=300`$ kV) settings on April 13 and 14 (see Fig. 2), just before and after the bunch length measurements to be presented in the next section. We note that the results of the two measurements are essentially the same, and they are also essentially identical to the March results with the beam off-resonance and $`V_c=300`$ kV (Fig. 1c). This indicates that the machine conditions, as concerns intra-beam scattering (the horizontal and vertical beam sizes; the lattice), are the same during all these measurements, as well as during the off-coupling, bunch length measurements presented below.

Finally, in Fig. 3 we compare the fits to the $`V_c=150`$ and 300 kV, off-coupling measurements with those to measurements—under the same conditions—taken in Dec 99 and presented in Ref. . We note that the corresponding curves for the two dates are significantly different. For example, for $`V_c=300`$ kV and $`I=1`$ mA the new results are 8% higher than the earlier results. This suggests that (when off coupling) the machine, concerning intra-beam scattering, was significantly different for the two dates.

## 3 Bunch Length Measurements

For the streak camera measurements first, to make sure that space charge in the camera itself was not a problem, tests were made at the highest current (2.5 mA) with different light filters, and an appropriate filter was chosen. Then, the main data taking process consisted of storing a high current beam and measuring the longitudinal bunch profile $`5070`$ times at fixed time intervals, while the current naturally decreased. Each trace, along with its DCCT current reading, was automatically saved to disk. To let in more light at the lower currents the light filter was automatically switched half-way through each data taking sequence. Measurements were performed with the beam on and off the coupling resonance, and with 4 different settings of rf voltage.

As before the streak camera traces were fit to an asymmetric Gaussian, given by

$$\lambda _z=\frac{A}{\sqrt{2\pi }\sigma _0}\mathrm{exp}\left[\frac{1}{2}\frac{(z\overline{z})^2}{\sigma _0^2(1\pm t)^2}\right]z\overline{z},$$
(1)
with the convention that more negative values of $`z`$ are more to the front of the bunch. The fitting constants are $`A`$, $`\sigma _0`$, $`\overline{z}`$, and the asymmetry factor $`t`$ (a constant, platform offset is also included in the fit). Note that the full-width-at-half-maximum (FWHM) of the fit is $`z_{fwhm}=2\sqrt{2\mathrm{ln}2}\sigma _0`$ and the rms bunch length $`\sigma _z=\sigma _0\sqrt{1+(38/\pi )t^2}`$. For small $`t`$, $`\sigma _z\sigma _0`$. The skew moment, defined by $`s=(zz)^3/\sigma _z^3`$, is given by $`s4t/\sqrt{2\pi }`$, for small $`t`$. Note also that from physical considerations, we generally expect $`t>0`$, i.e. the leading edge of the distribution to be steeper than the trailing edge. Four example scans, with their fits, are shown in Fig. 4. We note that the fits are reasonably good, though we see some, what are likely anomalous, deviations from the fits in the data. We have lots of data and we will use the statistical Method of Maximum Likelihood to do the error analysis.

The results for the beam on the coupling resonance and for rf voltages $`V_c=150`$, 200, 250, 300 kV are given in Figs. 5-8, and for the beam off the coupling resonance in Figs. 9-12. Shown are the parameters of the asymmetric Gaussian fit to the measured profiles: the area (a), the rms $`\sigma _0`$ (b), the asymmetry factor $`t`$ (c), and an estimate of the relative rms error in the fit (d). Plots (a)-(c) give the zeroth, the second, and the third moments of the charge distribution. The line in (b) is a straight line fit to the bunch length, where each data point has been weighted inversely by the variance in the fit to the asymmetric Gaussian \[note that for the fit $`y=xmfit+bfit`$\]. The line in (c) is a straight line fit to $`y=xmfit`$.

In the (a) part of the figures we note that the area is roughly proportional to current—with the constant of proportionality different for the low and the high current results due to the two different light filters used—though with some scatter in the data. In (b) we note that the scatter in the measured bunch lengths is quite significant, especially at the lower currents, with an rms deviation from the straight line fit varying from .35-.80 mm. In (c) we see even more scatter in the asymmetry factor. We expect this parameter to start at zero, and to increase with increasing current, and the data does seem to support an increasingly positive asymmetry factor with current. The results are near .1/mA, and slightly decrease with voltage. We will show in a future ATF report that this parameter gives us information about the higher mode losses and the real part of the impedance. Note that the measured asymmetry factors off-resonance are about half of those on-resonance, qualitatively consistent with the bunch being longer, and the higher mode losses being less, when off-resonance. Finally, in (d) we note that the relative rms error in the fit to the asymmetric Gaussian is $``$4-10%. There is more noise than in the measurements of Dec 99, reflected in the fact that the fits then were twice as good.

### 3.1 Comparisons and Discussion

In Fig. 13 we reproduce the linear fits to the bunch length parameter for the measurements with the beam off the coupling resonance (the solid lines). Note that these fits are not valid below $`I.4`$ mA, where there is no data; in fact, the bunch lengths are expected to drop below these lines for lower currents. In (a) the results are compared with the on-coupling results; in (b) with off-coupling results of Dec 99. In all cases the bunch length curves move up as the voltage moves down, as expected. Also, in (a) the bunch length is larger off the coupling resonance than on, due to the effect of intra-beam scattering, as expected.

We expect $`\sigma _z(0)`$, for $`V_c=200`$ kV, to be 6.2 mm. For the case of on-resonance, we find the value at the origin of the curves to be 3.7%, 10.4%, 4.6%, 5.0% larger than the zero value calculation for $`V_c=150`$, 200, 250, 300 kV, respectively. For the case of off-resonance the values are 4.5%, 6.4%, 15.3%, 27.1% larger than the zero value calculations. We expect the real bunch lengths to follow curves with negative curvature, so an overshoot is not unexpected. Note that these overshoots are significantly less than the $`35\%`$ for the Dec 99 results (see Fig. 13b), which, at the time, made us suspect a scaling error in the streak camera results.

In Sec. 2 we saw that, for the on-coupling case, the energy spread was almost independent of current. For example, by $`I=2`$ mA the energy spread has grown by only 3%. Here we see that at $`I=2`$ mA the bunch length has grown by 33-39%, depending on the rf voltage. Since IBS on the coupling resonance is very weak, this growth must be almost entirely due to PWBL. Thus, it appears that, after all, PWBL is a big factor in the ATF, and the machine must have a large inductive component of impedance. Note that this conclusion is different from that reached in Dec 99, where PWBL was estimated by dividing the bunch lengthening off-resonance by the energy spread growth. There it was found, for example, that the PWBL contribution to the bunch length at 2 mA, divided by the contribution at .5 mA, ranged from 7-15% depending on voltage. Here the equivalent (on-coupling) results range from 18-26%, a marked difference. Or, if we try to compare the present on and off-coupling results using the same prescription, we find that off-coupling the results are increases of 21.4% and 6.1%, for respectively the case $`V_c=150`$ and 300 kV; on-coupling the results are 23.5% and 20.1%. These also don’t agree. In making the earlier approximate calculation of PWBL we had assumed that PWBL can be treated as a perturbation on top of IBS. That is, we assumed that (1) PWBL doesn’t affect IBS significantly and (2) the bunch length increase due to PWBL can be added on top of that due to IBS, and the result would be similar whether on or off resonance. It appears that these assumptions taken together are wrong and the interaction of IBS and PWBL is, in fact, more complicated.

In Fig. 13b we note that the new off-coupling results are very different from those of Dec 99. We saw in Sec. 2 above that the energy spread vs. current was different for the two measurement days, presumably due to a difference in the IBS effect, so we would expect differences also in bunch length, since it is also affected by IBS. In addition the bunch length (below threshold) is affected by the impedance through PWBL. To say more, quantatively, about such measurements we would need a better understanding of how IBS and PWBL interact and, in addition, a knowledge of the transverse beam sizes and optics on the days of the measurements.

## 4 Conclusions

We have performed energy spread and bunch length measurements with the beam on and off the coupling resonance. Our energy spread results show that, with the beam on-resonance, the effect of intra-beam scattering is small, and that the threshold to the micro-wave instability is beyond 2.5 mA. Our on-resonance bunch length results show that potential well bunch lengthening is large—by 2 mA the bunch has lengthened by 33-39%—indicating that the impedance has a large inductive component.

Our energy spread and bunch length measurements with the beam off-resonance give results that are very different from those of Dec 99. This suggests that the horizontal and vertical emittances and/or the lattice were different than during the earlier measurements, resulting in a different intra-beam scattering effect. These measurements should be repeated both on and off the coupling resonance, under the same conditions—verifiable by finding the same energy spread vs. current curves—to see if the bunch length measurements are reproducible. We have assumed here that the streak camera measurements are sufficiently accurate for our purposes. The results appear to be consistent, but we have no independent way of checking their accuracy. (We might suggest installing in the ATF, sometime in the future, a bunch length measuring apparatus that works in a completely different way, such as the spectrum approach of Ref. , as an independent check on the streak camera results.)

In the future, when doing any such current-dependent measurements the energy spread vs. current measurement should always be performed, as a specifier of the intra-beam scattering machine conditions, and since it is a rather quick and simple measurement. To study specifically impedance questions, measurements should be performed with the beam on the coupling resonance, where intra-beam scattering is very weak. With the beam off-resonance, bunch length and energy spread measurements are important ingredients, along with the horizontal and vertical emittances and optics, for attempting a full understanding of intra-beam scattering and its interaction with the impedance in the ATF damping ring.

## 5 Acknowledgements

One of the authors (K.B.) thanks the ATF scientists and staff for their hospitality and help during his visits to the ATF, and M. Ross for his encouragement and support to make such visits.


text:
# Application of two-parameter dynamical replica theory to retrieval dynamics of associative memory with non-monotonic neurons

## 1 Introduction

The Hopfield model has attracted interests of researchers in various fields, and enormous amount of studies, both numerical-experimental and theoretical, have been carried out on it. Of particular interest among them is the one on the model with non-monotonic units: Whereas the conventional model uses, as the output function $`f`$ of a unit, monotonic functions such as $`f(x)=\mathrm{tanh}\beta x`$ ($`\beta >0`$), the model with non-monotonic units, or the non-monotonic model for short, uses a non-monotonic function. It has been reported that the non-monotonic model has various nice properties as a model of associative memory, such as enhancement of storage capacity, enlargement of basins of attraction associated with retrieval states.

In order to rigorously argue such properties of the non-monotonic model, theoretical analyses are necessary. However, attempts to analyze the non-monotonic model are often faced with difficulty because the non-monotonic model in general does not have a Lyapunov function, which would be a powerful analytical tool for characterizing the equilibrium as well as dynamical properties of the model. Thus, applicable theories are restricted to what have been devised for analyzing the conventional model and yet are independent of the functional form of the output function $`f`$. As for equilibrium analysis, the self-consistent signal-to-noise analysis (SCSNA) has been applied to the non-monotonic model and some interesting properties, including the existence of the so-called superretrieval states, have been found. For retrieval dynamics, on the other hand, currently no exact and tractable theory has been known even for the conventional model: The path integral formalism and Gardner-Derrida-Mottishaw theory (for “zero-temperature,” or $`\beta +\mathrm{}`$ case) are the exact theories for the asynchronous (or Glauber) and synchronous (or Little) dynamics, respectively, but computation of dynamics based on each of them is prohibitively difficult. It has been generally believed that any tractable theories on retrieval dynamics necessarily incorporate some approximation. For the conventional model with the synchronous dynamics, Amari-Maginu theory (for zero-temperature case; for extension to finite-temperature case, see ) has been proposed as one of such theories, and Nishimori and Opriş have applied it to the non-monotonic case. As one of tractable theories for the conventional model with the asynchronous dynamics, Coolen and Sherrington have proposed the 2-parameter dynamical replica theory (2-DRT hereafter; it is also sometimes called the Coolen-Sherrington (CS) theory). Of course it is an approximate theory, as confirmed, for example, by Ozeki and Nishimori and Tanaka and Osawa; nevertheless, it describes retrieval dynamics of the conventional model reasonably well. As Ozeki and Nishimori have mentioned, formulation of 2-DRT does not depend on the functional form of the output function $`f`$, and therefore it is possible to apply it to the non-monotonic model. Thus, the following question naturally arises: How well does 2-DRT describe retrieval dynamics of the non-monotonic model? In this paper, we address this problem, with emphasis placed on the storage capacity, size of basins of attraction, and the superretrieval states.

## 2 Model

Let us consider a model with $`N`$ units. Each unit has a binary variable $`s_i\{1,\mathrm{\hspace{0.17em}1}\}`$, $`i=1,\mathrm{},N`$, and $`𝒔=[s_1,\mathrm{},s_N]\{1,\mathrm{\hspace{0.17em}1}\}^N`$ represents a microscopic state of the model. Each unit stochastically and asynchronously updates the value of $`s_i`$ based on the current value of the “local field”

$$h_i(𝒔)=\underset{ji}{}J_{ij}s_j,$$
(1)
where $`J_{ij}`$ is a synaptic weight from neuron $`j`$ to neuron $`i`$. The probability of state flip $`s_i:=s_i`$ is given by the following transition probability

$$w_i(𝒔)=\frac{1}{2}(1s_if(h_i(𝒔))),$$
(2)
where $`f:𝑹[1,\mathrm{\hspace{0.17em}1}]`$ is the output function.

Taking $`f(x)=\mathrm{tanh}\beta x`$ yields the conventional Hopfield model. In this paper we consider the following non-monotonic function (Fig. 1):

$$f(x)=\{\begin{array}{cc}\mathrm{sgn}(x)\hfill & (|x|\theta )\hfill \\ \mathrm{sgn}(x)\hfill & (|x|<\theta )\hfill \end{array}$$
(3)
Its functional form is the same as that treated by Nishimori and Opriş. Throughout the paper, we follow the common convention about the time scale, that it is taken in such a way that the average number per unit time (frequency) of updates each unit executes is 1.

The model memorizes $`p=\alpha N`$ binary patterns $`𝝃^\mu =[\xi _1^\mu ,\mathrm{},\xi _N^\mu ]\{1,\mathrm{\hspace{0.17em}1}\}^N`$, $`\mu =1,\mathrm{},p`$, via the Hebb rule,

$$J_{ij}=\frac{1}{N}\underset{\mu =1}{\overset{p}{}}\xi _i^\mu \xi _j^\mu (ij).$$
(4)
The quantity $`\alpha p/N`$ is called the memory rate. We consider the case where the patterns to be memorized are randomly generated, that is, each of $`\xi _i^\mu `$’s takes the value $`\pm 1`$ with probability $`1/2`$, independently of others.

For measuring how well the model retrieves a pattern $`\mu `$, the correlation, or overlap,

$$m^\mu (𝒔)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\xi _i^\mu s_i$$
(5)
is used: $`|m^\mu (𝒔)|1`$ holds by definition, and if $`m^\mu (𝒔)=1`$, then the model is in the state $`𝒔=𝝃^\mu `$, and it is regarded as perfectly retrieving the pattern $`\mu `$. When $`m^\mu (𝒔)=1`$ the model is in the state $`𝒔=𝝃^\mu `$. The state is called the reversal state, but it can also be regarded as retrieving the pattern $`\mu `$ due to the symmetry of the model.

We assume that the model is going to retrieve one single pattern, so that $`m^\mu (𝒔)`$ are of order unity for that pattern only (the condensed ansatz), and that pattern $`\mu =1`$ is nominated for retrieval, without loss of generality. Then $`mm^1(𝒔)`$ can be taken as a macroscopic variable, or order parameter, describing how well the model retrieves the nominated pattern. If the model reaches an equilibrium with $`m0`$, the model is said to successfully retrieve the pattern, and such an equilibrium is called a retrieval state. The local field $`h_i(𝒔)`$ is now rewritten as

$$h_i(𝒔)=\xi _i^1[m+z_i(𝒔)]\frac{1}{N}s_i,z_i(𝒔)\xi _i^1\underset{\mu >1}{\overset{p}{}}\xi _i^\mu \frac{1}{N}\underset{ji}{}\xi _k^\mu s_j.$$
(6)
The term $`z_i(𝒔)`$ represents the interference in the local field from non-nominated patterns $`\mu >1`$, and is called the “noise” term. Although time evolution of the microscopic state of the model is stochastic in nature, it is observed that the evolution of the order parameter $`m`$ in the course of pattern retrieval can be often seen as being governed by a certain deterministic law. To describe this retrieval dynamics is one of important problems in this field. Path integral formalism provides an exact description to the retrieval dynamics, but it requires parametrization with infinite degrees of freedom and hence practically intractable.

## 3 Two-parameter dynamical replica theory (2-DRT)

2-DRT, proposed by Coolen and Sherrington, provides a tractable, and yet reasonably well for the conventional model, description of the retrieval dynamics. It uses two parameters, $`m`$ and $`r`$, as the order parameters, the latter being defined as

$$r\frac{1}{\alpha }\underset{\mu >1}{}(m^\mu (𝒔))^2,$$
(7)
which intuitively represents the degree of interferential effect of non-nominated patterns $`\mu >1`$ onto the retrieval of pattern 1. 2-DRT derives deterministic flow equations for these two order parameters. Without any assumptions, one cannot expect that the flow equations for these two parameters are closed, and thus the time evolution of $`m`$ and $`r`$ cannot be completely determined by current values of them. In the framework of DRT in general, the following two assumptions are made:

1. Self-averaging of flow equations with respect to randomness of the system (randomly chosen patterns for the case treated in this paper).
2. Probability equipartitioning within subshells: when values of the order parameters are given, the probability distribution of the corresponding microscopic state can be regarded as being uniform over the subshell (the set of microscopic states which have the specified values of the order parameters), with regard to calculation of the flow equations.

Owing to these closing assumptions, one can derive the deterministic flow equations for $`m`$ and $`r`$:

$`{\displaystyle \frac{dm}{dt}}={\displaystyle 𝑑zD_{m,r}[z]f(m+z)}m`$
$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{dr}{dt}}={\displaystyle \frac{1}{\alpha }}{\displaystyle 𝑑zD_{m,r}[z]zf(m+z)}+1r`$
where $`D_{m,r}[z]`$ is the distribution of the noise terms. Replica calculation gives, within the replica-symmetric (RS) ansatz, the RS solution $`D_{m,r}^{\mathrm{RS}}[z]`$ for the noise distribution, which has been derived by CS as

$`D_{m,r}^{\mathrm{RS}}[z]`$ $`=`$ $`{\displaystyle \frac{e^{(\mathrm{\Delta }+z)^2/2\alpha r}}{2\sqrt{2\pi \alpha r}}}\left\{1{\displaystyle Dy\mathrm{tanh}\left[\lambda y\left(\frac{\mathrm{\Delta }}{\rho \alpha r}\right)^{1/2}+(\mathrm{\Delta }+z)\rho \frac{r_{\mathrm{AGS}}}{r}+\mu \right]}\right\}`$ (8)
$`+{\displaystyle \frac{e^{(\mathrm{\Delta }z)^2/2\alpha r}}{2\sqrt{2\pi \alpha r}}}\left\{1{\displaystyle Dy\mathrm{tanh}\left[\lambda y\left(\frac{\mathrm{\Delta }}{\rho \alpha r}\right)^{1/2}+(\mathrm{\Delta }z)\rho \frac{r_{\mathrm{AGS}}}{r}\mu \right]}\right\}`$
$$\mathrm{\Delta }=\rho \alpha (rr_{\mathrm{AGS}})$$
(9)
$$r_{\mathrm{AGS}}=\frac{\lambda ^2}{\rho ^2\alpha }$$
(10)
where $`Dy=(dy/\sqrt{2\pi })e^{y^2/2}`$ is the Gaussian measure. The parameters $`\{q,\lambda ,\rho ,\mu \}`$ are to be determined from $`m`$, $`r`$, and $`\alpha `$ by the following saddle-point equations:

$$\begin{array}{cc}r=\frac{1\rho (1q)^2}{[1\rho (1q)]^2},\hfill & \lambda =\frac{\rho \sqrt{\alpha q}}{1\rho (1q)}\hfill \\ & \\ m=Dy\mathrm{tanh}(\lambda y+\mu ),\hfill & q=Dy\mathrm{tanh}^2(\lambda y+\mu )\hfill \end{array}$$
(11)
It should be noted that the replica calculation of the noise distribution bears no relation to the dynamics of the model and the choice of the output function $`f`$; it executes averaging over the ($`m`$, $`r`$)-subshell with uniform measure, and therefore the calculation is not dynamical but configurational. This is the reason why, apart from validity of the two assumptions, 2-DRT can be straightforwardly applied to the non-monotonic model.

The latter of the two above-mentioned assumptions, the equipartitioning assumption, is a critical one because it has been known that it is not valid for the conventional model, as well as the continuous-valued linear system, or the Langevin spin system. To show this directly for the conventional model, Tanaka and Osawa have proposed a dynamics, called ($`m`$, $`r`$)-annealing. The ($`m`$, $`r`$)-annealing is defined, on the basis of the solution $`\{\rho ,\mu \}`$ of the saddle-point equations for given $`m`$ and $`r`$, as the dynamics of the conventional model with the inverse temperature $`\rho `$ (that is, it uses $`f(x)=\mathrm{tanh}\rho x`$), but an extra bias is added to the local field,

$$h_i(𝒔)=\underset{ji}{}J_{ij}s_j+b\xi _i^1,$$
(12)
where $`b\mu /\rho m`$. It has been shown that the ($`m`$, $`r`$)-annealing executes Monte Carlo sampling from a ($`m`$, $`r`$)-subshell with uniform probability (in the limit $`N\mathrm{}`$), hence realizing the equipartitioning. It is useful in investigating validity of the equipartitioning assumption, and will be utilized in this paper.

## 4 Results

### 4.1 Time evolution of order parameters

We first examined whether 2-DRT describes overall characteristics of the dynamics of the non-monotonic model. We assume the common convention that initial microscopic states of the model are given by randomly corrupting the nominated pattern; that is, the initial state $`𝒔`$ is set by the probability law $`\mathrm{Prob}[s_i]=\delta (s_i\xi _i^1)(1+m_0)/2+\delta (s_i+\xi _i^1)(1m_0)/2`$, so that the initial overlap $`m(t=0)`$ approximately equals to $`m_0`$ when $`N`$ is sufficiently large. Following this initialization procedure, the initial value of $`r`$, when $`N`$ is sufficiently large, approximately equals to $`1`$. We found that 2-DRT does describe the dynamics of the model considerably well, as shown in Fig. 2 for the case with $`\alpha =0.2`$ and $`\theta =1.4`$. 2-DRT reproduced the trajectories almost exactly when retrieval succeeded. Noticeable disagreement between simulations and 2-DRT was seen for cases where retrieval failed. Characteristic aspect of the disagreement is that the trajectories predicted by 2-DRT exhibit, in early stages, overall slowing down against the corresponding simulations. These observations are essentially the same as those found in the conventional model. The observed disagreement is due to the failure of the equipartitioning assumption of 2-DRT just as in the conventional model, as demonstrated by the following numerical experiment. The procedure of this experiment is as follows: Simulate a model with $`\theta =1.4`$ and $`\alpha =0.2`$ by setting its initial condition as $`(m,r)=(0.1,\mathrm{\hspace{0.17em}1})`$, stop it at $`t=1`$, and then execute the ($`m`$, $`r`$)-annealing for 80 units of time. The flow vectors $`(\dot{m},\dot{r})`$, evaluated from (1) the model just before executing the ($`m`$, $`r`$)-annealing, (2) the model just after it, and (3) 2-DRT, were compared, and the result is summarized in Fig. 3. This shows that the flow vectors of the model after the ($`m`$, $`r`$)-annealing and of 2-DRT are almost the same, whereas the flow vector before the ($`m`$, $`r`$)-annealing is different from these two, which means that the equipartitioning assumption does not hold in this case.

### 4.2 Capacity and basins of attraction

In this paper, we define the storage capacity $`\alpha _c`$ as the maximum value of $`\alpha `$ for which a stable macroscopic state with nonvanishing $`m`$ exists. We call an equilibrium macroscopic state with nonvanishing $`m`$ the retrieval state. Since the retrieval state may be unstable, the condition for the existence of a retrieval state will give an overestimate of the true storage capacity. On the other hand, a stable retrieval state may have a very small basin of attraction, so that we may fail to find such a retrieval state in numerical experiments even though it is stable.

Shiino and Fukai have analyzed equilibrium properties of continuous-valued continuous-time non-monotonic models using SCSNA. Although we consider the model with binary variables in this paper rather than ones with continuous values, the equilibrium condition is shown to be the same as that for the continuous-value continuous-time model owing to the current choice of the output function $`f`$ (eq. (3)), and thus SCSNA can be applied to the model treated here. We executed the SCSNA calculation on this model, and Fig. 4 shows the result for $`\alpha _e`$, the maximum of $`\alpha `$ for which a retrieval state exists, versus $`\theta `$. When $`\theta \mathrm{}`$, $`\alpha _e`$ approaches the well-known value 0.138, confirming that SCSNA is consistent with the Amit-Gutfreund-Sompolinsky (AGS) theory. As $`\theta `$ becomes smaller, $`\alpha _e`$ increases so that it reaches its maximal value $`\alpha _e=0.489`$ at $`\theta 0.7`$.

As we have already discussed, $`\alpha _e`$ gives an overestimation of the true storage capacity $`\alpha _c`$. To evaluate $`\alpha _c`$ itself, we have to take into account the dynamical aspect of the retrieval process. This discussion leads us to the idea to apply 2-DRT to determine the storage capacity $`\alpha _c`$, by observing whether or not the trajectories from arbitrary initial conditions approach a retrieval state.

Two points have to be mentioned here. First, although 2-DRT becomes exact at equilibrium for the conventional model, it is no longer so for the non-monotonic model, which means that 2-DRT may not reproduce the storage capacity. Second, since we assume the initialization procedure described above, only the retrieval states which can be reachable from initial conditions with $`r=1`$ are to be observed in the simulations. To make correspondence with this experimental setup, we estimated the storage capacity $`\alpha _c`$ by tracking 2-DRT trajectories from initial conditions with $`r=1`$. The storage capacity estimated by the simulations and by the 2-DRT trajectory tracking may be therefore an underestimate against the true storage capacity, because there might be stable retrieval states unreachable from any initial condition with $`r=1`$ (see Sect. 4.3).

Figure 5 shows the estimated storage capacity $`\alpha _c`$ by 2-DRT trajectory tracking and by simulations. It can be seen that 2-DRT trajectory tracking well reproduces the storage capacity obtained by the simulations for the whole range of $`\theta `$. Comparing it with Fig. 4 reveals that, while the agreement is good when $`\theta `$ is large, the discrepancy becomes apparent as $`\theta `$ becomes less than about $`1.5`$, showing that the storage capacity estimated by 2-DRT is considerably smaller than that estimated by SCSNA. It may be partly explained by the overestimation of SCSNA described above.

One of the main advantage of dynamical theory is that it allows us to evaluate basins of attraction, because they are essentially of dynamical nature. We used 2-DRT to evaluate the basin of attraction of the retrieval states. Since we adopt the above-mentioned initialization convention, we consider the basin of attraction as being represented in terms of $`m`$ only, that is, we regard a value of $`m_0`$ as belonging to the basin of attraction of the retrieval state if the trajectory starting at the state $`(m,r)=(m_0,\mathrm{\hspace{0.17em}1})`$ approaches the retrieval state. Evaluating the basin of attraction defines the critical initial overlap $`m_c`$, which means that initial states $`(m_0,\mathrm{\hspace{0.17em}1})`$ with $`m_0m_c`$ yield successful retrieval. Figures 6 and 7 show the critical initial overlap $`m_c`$ and the values of $`m`$ of the retrieval state, respectively, evaluated by 2-DRT and by simulations. It is readily seen that enlargement of the basin of attraction occurs as $`\theta `$ becomes small, and that 2-DRT captures this phenomenon reasonably well.

### 4.3 Superretrieval states

As a result of SCSNA analysis on non-monotonic models, Shiino and Fukai have shown that there is a phase where equilibrium states corresponding to “perfect” retrieval exist. Such states are called the superretrieval states, whose existence is one of unique features of the non-monotonic models. Here, “perfect” means that the correlation of the sign of the local field $`h_i(𝒔)`$ (not of $`𝒔`$) with the nominated pattern $`𝝃^1`$ is exactly equal to $`\pm 1`$. The correlation defined as above is called the tolerance overlap. The critical capacity $`\alpha _0`$, below which the superretrieval occurs, can be evaluated numerically by SCSNA, and is also shown in Fig. 4.

An explanation, given by Shiino and Fukai, for possibility of such states is briefly as follows: For such states $`r_{\mathrm{AGS}}+0`$ holds, which has been confirmed by numerically solving relevant self-consistent equations. Since variance of the noise term (without the “systematic” term, $`\mathrm{\Gamma }Y`$, in their terminology) is given by $`\alpha r_{\mathrm{AGS}}`$ in SCSNA, it implies that the effect of the noise completely vanishes in the states, which enables the superretrieval.

Especially interesting is whether or not 2-DRT provides an appropriate description of dynamics which is bound for superretrieval states. To see this, we examined the case where $`\theta =0.4`$ and $`\alpha =0.05`$, which, according to SCSNA analysis, is expected to have the superretrieval states. Figure 8 shows time evolution of $`r`$ evaluated by simulation and by 2-DRT, with initial condition $`m_0=0.9`$. Again, 2-DRT reproduced the simulation result fairly well. For the simulation, state transition ended (confirmed numerically) at $`t30`$, where $`(m,r)=(0.398,\mathrm{\hspace{0.17em}0.00440})`$. The tolerance overlap was evaluated to be exactly equal to $`1`$ at this state, indicating that this is the superretrieval state. Figure 9 shows the noise distribution $`D[z]`$ at the equilibrium state achieved by the simulation, and the one computed by 2-DRT for the same values of $`(m,r)`$. They are in good agreement, suggesting that 2-DRT can successfully predict the trajectories even in the case of superretrieval, provided that the system size $`N`$ is sufficiently large. As for 2-DRT, the value of $`r`$ continued to decrease at $`t10^5`$ (we stopped computation of 2-DRT at $`t=10^5`$ because rounding errors became profound beyond this point), where $`(m,r)=(0.399,\mathrm{\hspace{0.17em}0.00159})`$. We observed numerically that $`t`$-dependence of the decrease of $`r`$ can be expressed, to a good approximation, as $`r(t)1/\mathrm{ln}t`$ (Fig. 10), which strongly supports the conjecture that the trajectory obtained by 2-DRT certainly approaches $`r=0`$. As can be seen by eq. (8), variance of the noise term is, roughly speaking, given by $`\alpha r`$ in 2-DRT. Then, if the superretrieval states are described appropriately by 2-DRT, they should correspond to the states with $`r=+0`$. We therefore investigate in the following the solutions of the saddle-point equations (11) when $`r=+0`$.

By formally taking the limit $`r+0`$, the saddle-point equations (11) are reduced to the following equations

$`\rho `$ $`=`$ $`\mathrm{}`$ (13)
$`m`$ $`=`$ $`{\displaystyle Dy\mathrm{tanh}(\lambda y+\mu )}`$ (14)
$`q`$ $`=`$ $`{\displaystyle Dy\mathrm{tanh}^2(\lambda y+\mu )}`$ (15)
$`\lambda `$ $`=`$ $`{\displaystyle \frac{\sqrt{\alpha q}}{1q}},`$ (16)
and, correspondingly, the RS noise distribution $`D_{m,r=+0}^{\mathrm{RS}}[z]`$ (8) becomes

$$D_{m,r=+0}^{\mathrm{RS}}[z]=\frac{1m}{2}\delta (z\alpha )+\frac{1+m}{2}\delta (z+\alpha ).$$
(17)
The condition that the macroscopic state $`(m,r=+0)`$ is an equilibrium state, that is, $`(\dot{m},\dot{r})=(0,\mathrm{\hspace{0.17em}0})`$ holds, is thus given by

$$f(m\pm \alpha )=1.$$
(18)
For the current choice of the function $`f`$ (eq. (3)), this condition is satisfied when $`0<m\alpha <\theta <m+\alpha `$, or equivalently, $`\mathrm{max}\{\theta \alpha ,\alpha \}<m<\mathrm{min}\{\theta +\alpha ,\mathrm{\hspace{0.17em}1}\}`$. This is in consistent with the observation from the simulations that $`m`$ tends to approach $`\theta `$ when the superretrieval occurs.

The applicability of the argument above depends not only on the two assumptions made at the beginning but also on the two following points: The first one concerns the so-called freezing line, which defines the points in the $`(m,r)`$ plane where the number of microscopic states within the $`(m,r)`$-subshell changes from an exponentially large number to an exponentially small number in terms of $`N`$. The second one concerns the so-called de Almeida-Thouless (AT) line, at which the replica-symmetry-breaking (RSB) occurs and thus the RS ansatz becomes no longer valid.

The freezing line is given, under the RS ansatz, by the following equation:

$$Dy\mathrm{ln}\mathrm{cosh}[\lambda y+\mu ]\mu m\frac{\alpha }{2}\left[\mathrm{ln}[1\rho (1q)]+\frac{\rho (1q)(1\rho +3q\rho )}{[1\rho (1q)]^2}\right]+\mathrm{ln}2=0$$
(19)
The number of microscopic states within the subshell is exponentially large, and therefore taking averages over the subshell is expected to have the proper meaning, as long as the left-hand side of eq. (19) is positive. For $`r=1`$ and $`|m|<1`$ the left-hand side becomes

$$\frac{1}{2}\left[(1+m)\mathrm{ln}\frac{1+m}{2}+(1m)\mathrm{ln}\frac{1m}{2}\right]>0,$$
(20)
so that the number of microscopic states is indeed exponentially large. In the limit $`r+0`$, however, the left-hand side of eq. (19) goes to $`\mathrm{}`$, implying that $`r=0`$ is outside the freezing line. This means that there are so few microscopic states near $`r=0`$, which may in part explain the observation from the simulations that the trajectories approaching $`r=0`$ reach equilibrium before actually arriving at $`r=0`$. That $`r=0`$ is outside the freezing line also means that the argument with the formal limit $`r+0`$ eventually loses its proper meaning, because the relevant subshell average is over an exponentially small number of microscopic states. Numerical evaluation reveals, however, that the freezing line for $`r<1`$ lies very close to $`r=0`$ (Fig. 11). Moreover, the saddle-point solution of the order parameters changes smoothly as $`r`$ tends to $`+0`$: Let $`\mu _0`$, $`q_0`$ be the asymptotic values of the saddle-point solution $`\mu `$, $`q`$, respectively. Then the asymptotic form of the saddle-point solution as $`r+0`$ is given by

$`\rho `$ $`=`$ $`r^1+O(1)`$
$`\lambda `$ $`=`$ $`{\displaystyle \frac{\sqrt{\alpha q_0}}{1q_0}}+O(r)`$
$`\mu `$ $`=`$ $`\mu _0+O(r)`$
$`q`$ $`=`$ $`q_0+O(r)`$
$`\mathrm{\Delta }`$ $`=`$ $`\alpha +O(r).`$ (21)
We can therefore expect that the argument presented above on the formal limit $`r+0`$ well captures qualitative aspects of the equilibrium states achieved by the simulations.

The AT line is determined by examining stability of the RS solution. Assuming that RSB is caused by destabilization of the so-called “replicon” modes for the case $`r<1`$ just as it has been assumed for the case $`r1`$, the AT line turns out to be given by the same formula as given in  for $`r1`$:

$$\alpha =\rho ^2(\alpha +\mathrm{\Delta })^2\frac{Dy}{\mathrm{cosh}^4(\lambda y+\mu )}$$
(22)
The RS solution is valid if the right-hand side of eq. (22) is less than $`\alpha `$. We confirmed that, for the case where $`\alpha =0.05`$, the AT line lies in the region $`m>0.890`$ for $`r<1`$ (Fig. 11), which implies that the superretrieval observed in the simulations was irrelevant to RSB.

We have conducted local stability analysis of the 2-DRT stationary solutions $`r=0`$, $`\mathrm{max}\{\theta \alpha ,\alpha \}<m<\mathrm{min}\{\theta +\alpha ,\mathrm{\hspace{0.17em}1}\}`$ corresponding to the superretrieval state. The result of the analysis states that, among the superretrieval solutions, $`(m,r)=(\theta ,\mathrm{\hspace{0.17em}0})`$ is the only attractor when $`2\alpha <\theta <1`$, although it becomes unstable when $`\alpha <\theta <\mathrm{min}\{2\alpha ,\mathrm{\hspace{0.17em}1}\}`$. For the details of the local stability analysis, see Appendix. In the simulations, however, it is not for all $`(\theta ,\alpha )`$ values satisfying $`2\alpha <\theta <1`$ that the superretrieval was observed. Figure 12 shows the region where the local stability analysis of the 2-DRT predicts the stable superretrieval state, the one where the superretrieval is observed by simulations (in the sense that tolerance overlap is numerically evaluated to be 1), and the one where the superretrieval is predicted by tracking the 2-DRT trajectories starting at $`(m,r)=(m_0,\mathrm{\hspace{0.17em}1})`$. The result of 2-DRT trajectory tracking and that of simulations are in good agreement with each other. Both are within the region where the superretrieval state is locally stable, as they should be, but apparently they do not coincide. The region where the superretrieval is observed in simulations may be further restricted by the following factors:

* The superretrieval solution $`(m,r)=(\theta ,\mathrm{\hspace{0.17em}0})`$ may not be reachable from the conventional initial states with $`r=1`$, even though it is an attractor.
* The superretrieval solution $`(m,r)=(\theta ,\mathrm{\hspace{0.17em}0})`$ may be at the outside of the RS region, where the stationarity and local stability arguments, both based on the RS ansatz, are no longer valid.

A demonstration regarding the former factor is shown in Fig. 13. For the condition $`(\theta ,\alpha )=(0.4,\mathrm{\hspace{0.17em}0.15})`$ (marked by a cross in Fig. 12), for example, the superretrieval state is not observed by following time evolution by either numerical simulation or 2-DRT. Nevertheless, 2-DRT predicts that under this condition the stable superretrieval state exists at $`(m,r)=(0.4,\mathrm{\hspace{0.17em}0})`$. As shown in the figure, 2-DRT trajectory tracking shows that in this condition the superretrieval state is indeed stable, but it is not reachable from the initial states with $`r=1`$.

An interesting observation related to the latter factor is that there is a rough numerical correspondence between the region where the equilibrium superretrieval solution $`(m,r)=(\theta ,\mathrm{\hspace{0.17em}0})`$, $`\alpha <\theta <1`$, which may not be stable, satisfies the RS ansatz, and the region where SCSNA predicts superretrieval to occur, as shown in Fig. 14.

## 5 Conclusion

We have studied the question of how well 2-DRT describes retrieval dynamics of the non-monotonic model. Although there is no theoretical justification for 2-DRT to be exact either for the non-monotonic model, 2-DRT turns out to reproduce the retrieval dynamics quite well, and it gives reasonable results as for the capacity, basins of attractions, and the superretrieval states.

## Appendix A Local Stability Analysis of Superretrieval States

We first split the RS noise distribution $`D_{m,r}^{\mathrm{RS}}[z]`$ into two components, as follows:

$$D_{m,r}^{\mathrm{RS}}[z]D^{}[z]+D^+[z]$$
(23)
$$D^\pm [z]=\frac{e^{(\mathrm{\Delta }\pm z)^2/2\alpha r}}{2\sqrt{2\pi \alpha r}}w_\pm (z)$$
(24)
$$w_\pm (z)=1Dy\mathrm{tanh}\left[\lambda y\left(\frac{\mathrm{\Delta }}{\rho \alpha r}\right)^{1/2}+(\mathrm{\Delta }\pm z)\rho \frac{r_{\mathrm{AGS}}}{r}\pm \mu \right]$$
(25)
Note that $`w^\pm (z)`$ is “slowly varying” with respect to $`z`$, because

$$\left|\frac{dw^\pm (z)}{dz}\right|<\left|\rho \frac{r_{\mathrm{AGS}}}{r}\right|$$
(26)
holds and the bound $`|\rho r_{\mathrm{AGS}}/r|`$ remains finite even when $`r+0`$. We can thus regard that each component is basically a gaussian distribution centered at $`z=\mathrm{\Delta }`$ and width $`O(\sqrt{\alpha r})`$, and it has been modulated by a bounded, monotonic, and slowly varying function $`0<w_\pm (z)<2`$. From the asymptotic form of the saddle-point solution as $`r+0`$ (eq. (21)), we can expect, for small $`r`$, that the noise components $`D^\pm [z]`$ become sharply peaked around $`z=\mathrm{\Delta }\pm \alpha `$. In the limit $`r+0`$, we have

$$D^\pm [z]=\frac{1m}{2}\delta (z\alpha ),$$
(27)
so that the condition

$$\mathrm{max}\{\theta \alpha ,\alpha \}<m<\mathrm{min}\{\theta +\alpha ,\mathrm{\hspace{0.17em}1}\}$$
(28)
is obtained for the existence of equilibrium states of the form $`(m,r)=(m,\mathrm{\hspace{0.17em}0})`$, as discussed in Sect. 4.3.

In this section we analyze local stability of the equilibrium states $`(m,r=0)`$ satisfying the condition (28). Using the noise components, the time evolution equations are rewritten as

$`\dot{m}={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑zD^{}[z]f(m+z)+{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑zD^+[z]f(m+z)m`$
$`{\displaystyle \frac{1}{2}}\dot{r}={\displaystyle \frac{1}{\alpha }}\left[{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑zD^{}[z]zf(m+z)+{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑zD^+[z]zf(m+z)\right]+1r.`$ (29)
Because $`D^\pm [z]`$ are sharply peaked, as the first step of approximation we can assume that $`f(m+z)=f(m\mathrm{\Delta })`$ in the integrals with $`D^\pm [z]`$. This assumption becomes exact in the limit $`r+0`$ and when $`f(m+z)`$ is continuous around $`z=\pm \alpha `$, but for finite $`r`$ it gives an approximate result and the approximation error comes from the contribution of the tails of $`D^\pm [z]`$ where $`f(m+z)`$ changes the sign. For explanation purposes we introduce the following four regions:

$`\mathrm{I}`$ $``$ $`\{zm+z<\theta \}`$
$`\mathrm{II}`$ $``$ $`\{z\theta <m+z<0\}`$
$`\mathrm{III}`$ $``$ $`\{z0<m+z<\theta \}`$
$`\mathrm{IV}`$ $``$ $`\{z\theta <m+z\}`$ (30)
$`f(m+z)=1`$ for $`z\mathrm{I}`$ or $`\mathrm{III}`$, and $`f(m+z)=1`$ for $`z\mathrm{II}`$ or $`\mathrm{IV}`$. The equilibrium states which we are interested in correspond to the case where the peak of $`D^+[z]`$ is in the region IV and that of $`D^{}[z]`$ in the region III. In this case the time evolution equations are approximated to be

$`\dot{m}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑zD^{}[z]{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑zD^+[z]m`$
$`{\displaystyle \frac{1}{2}}\dot{r}{\displaystyle \frac{1}{\alpha }}\left[{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑zD^{}[z]z{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑zD^+[z]z\right]+1r.`$ (31)
However, direct calculation shows that the right-hand sides of these equations exactly equal to 0. This fact indicates that the time evolution near $`r=0`$ should be governed by the contribution of the tails.

The principal contribution comes from the largest one of the following three quantities:

1. Contribution of the tail of $`D^+[z]`$ in the region III:

$$I_1=2_{\mathrm{III}}𝑑zD^+[z]a(z)2_{\mathrm{}}^{z_0}𝑑zD^+[z]a(z)$$
(32)
2. Contribution of the tail of $`D^{}[z]`$ in the region IV:

$$I_2=2_{\mathrm{IV}}𝑑zD^{}[z]a(z)=2_{z_0}^{\mathrm{}}𝑑zD^{}[z]a(z)$$
(33)
3. Contribution of the tail of $`D^{}[z]`$ in the region II:

$$I_3=2_{\mathrm{II}}𝑑zD^{}[z]a(z)2_{\mathrm{}}^m𝑑zD^{}[z]a(z)$$
(34)

where $`z_0\theta m`$, and $`a(z)=1`$ or $`z`$, depending on which of $`\dot{m}`$ and $`\dot{r}`$ we are considering. We approximate each contribution by extending the integral region to $`\mathrm{}`$ or $`\mathrm{}`$. In fact this approximation does not affect the final result in the $`r+0`$ limit because it changes each quantity by a vanishingly small amount.

First let us consider the contribution to $`\dot{m}`$. Evaluation for $`I_1`$ yields

$`I_1`$ $``$ $`2{\displaystyle _{\mathrm{}}^{z_0}}𝑑zD^+[z]`$ (35)
$`=`$ $`{\displaystyle _{\mathrm{}}^{z_0}}𝑑z{\displaystyle \frac{e^{(\mathrm{\Delta }+z)^2/2\alpha r}}{\sqrt{2\pi \alpha r}}}w_+(z)`$
$``$ $`{\displaystyle _{\mathrm{}}^{z_0}}𝑑z{\displaystyle \frac{e^{(\mathrm{\Delta }+z)^2/2\alpha r}}{\sqrt{2\pi \alpha r}}}w_+(z_0)`$
$``$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}\mathrm{exp}\left[{\displaystyle \frac{(z_0+\mathrm{\Delta })^2}{2\alpha r}}\mathrm{ln}{\displaystyle \frac{z_0+\mathrm{\Delta }}{\sqrt{\alpha r}}}+\mathrm{ln}w_+(z_0)\right].`$
Similarly, for $`I_2`$ and $`I_3`$, we have

$`I_2`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}\mathrm{exp}\left[{\displaystyle \frac{(z_0\mathrm{\Delta })^2}{2\alpha r}}\mathrm{ln}{\displaystyle \frac{z_0\mathrm{\Delta }}{\sqrt{\alpha r}}}+\mathrm{ln}w_{}(z_0)\right],`$ (36)
$`I_3`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}\mathrm{exp}\left[{\displaystyle \frac{(m\mathrm{\Delta })^2}{2\alpha r}}\mathrm{ln}{\displaystyle \frac{m\mathrm{\Delta }}{\sqrt{\alpha r}}}+\mathrm{ln}w_{}(m)\right],`$ (37)
respectively. In the $`r+0`$ limit, the dominant contribution comes from the first term of the exponent for each case, so that comparison of the term is sufficient to determine which of $`I_1`$, $`I_2`$, and $`I_3`$ has the largest contribution to $`\dot{m}`$. The result of the comparison for small $`r`$ is summarized as follows:

* When $`2\mathrm{\Delta }<\theta `$, the largest contribution comes from $`I_1`$ or $`I_2`$. If $`m<\theta `$, $`I_1`$ is the largest and $`\dot{m}>0`$. Otherwise, $`I_2`$ is the largest and $`\dot{m}<0`$.
* When $`2\mathrm{\Delta }>\theta `$, the largest contribution comes from $`I_2`$ or $`I_3`$. Since both $`I_2`$ and $`I_3`$ have negative contribution, $`\dot{m}<0`$.

From this result, we can conclude that the stable superretrieval state, if it exists, should be $`(m,r)=(\theta ,\mathrm{\hspace{0.17em}0})`$, and that $`2\alpha <\theta <1`$ is a necessary condition for the existence of the stable superretrieval state.

Let us now take a closer look at the flow near the state $`(m,r)=(\theta ,\mathrm{\hspace{0.17em}0})`$. We let $`\epsilon z_0=m\theta `$, and consider time evolution of the two small quantities, $`\epsilon `$ and $`r`$. In the following arguments we assume that $`2\mathrm{\Delta }<\theta `$ holds, so that the principal contribution comes from $`I_1`$ or $`I_2`$, but not from $`I_3`$. Under this assumption, we have

$`\dot{\epsilon }`$ $``$ $`\sqrt{{\displaystyle \frac{2\alpha r}{\pi }}}e^{(\mathrm{\Delta }^2+\epsilon ^2)/2\alpha r}\overline{w}\left({\displaystyle \frac{1}{\mathrm{\Delta }}}\right)\left[\mathrm{tanh}{\displaystyle \frac{\mathrm{\Delta }\epsilon }{\alpha r}}+{\displaystyle \frac{\delta }{\overline{w}}}\right]\mathrm{cosh}{\displaystyle \frac{\mathrm{\Delta }\epsilon }{\alpha r}}`$
$`{\displaystyle \frac{1}{2}}\dot{r}`$ $``$ $`\sqrt{{\displaystyle \frac{2\alpha r}{\pi }}}e^{(\mathrm{\Delta }^2+\epsilon ^2)/2\alpha r}\overline{w}\left({\displaystyle \frac{r}{\mathrm{\Delta }^2}}\right)\left[\left({\displaystyle \frac{\mathrm{\Delta }\epsilon }{\alpha r}}{\displaystyle \frac{\delta }{\overline{w}}}\right)\mathrm{tanh}{\displaystyle \frac{\mathrm{\Delta }\epsilon }{\alpha r}}+\left({\displaystyle \frac{\delta }{\overline{w}}}{\displaystyle \frac{\mathrm{\Delta }\epsilon }{\alpha r}}1\right)\right]\mathrm{cosh}{\displaystyle \frac{\mathrm{\Delta }\epsilon }{\alpha r}},`$ (38)
where

$`\overline{w}`$ $``$ $`{\displaystyle \frac{1}{2}}\left(w_+(z_0)+w_{}(z_0)\right)`$
$`\delta `$ $``$ $`{\displaystyle \frac{1}{2}}\left(w_+(z_0)w_{}(z_0)\right).`$ (39)
Note that $`1<\delta /\overline{w}<1`$ holds.

Assuming that $`\epsilon /r`$ remains finite, we can readily see that $`\dot{r}`$ is smaller in magnitude than $`\dot{\epsilon }`$ by a factor $`r`$. Then for small enough $`r`$ the slaving principle applies and $`\epsilon `$ is expected to relax toward its equilibrium value much faster than $`r`$. This justifies the adiabatic approximation, and we can regard that the equilibrium condition for $`\epsilon `$,

$$\mathrm{tanh}\frac{\mathrm{\Delta }\epsilon }{\alpha r}+\frac{\delta }{\overline{w}}=0,$$
(40)
holds throughout the dynamics. This is indeed consistent with the assumption that $`\epsilon /r`$ remains finite. $`\dot{r}`$ is then given by

$$\frac{1}{2}\dot{r}\sqrt{\frac{2\alpha r}{\pi }}e^{(\mathrm{\Delta }^2+\epsilon ^2)/2\alpha r}\overline{w}\left(\frac{r}{\mathrm{\Delta }^2}\right)\left[1\left(\frac{\delta }{\overline{w}}\right)^2\right]^{1/2}<0.$$
(41)
This shows that the state $`(m,r)=(\theta ,\mathrm{\hspace{0.17em}0})`$ is actually a stable point of the dynamics described by 2-DRT under the condition $`2\alpha <\theta <1`$.