text:
# Remnants from Gamma-Ray Bursts

## 1 Introduction

If a $`\gamma `$-ray burst (GRB) originate within a galactic disk then the large deposition of energy will result in a blast wave whose initial phase produces the afterglow. The late phase of the blast wave evolution would result, as noted by Chevalier (1974) in the context of supernova remnants (SNRs), in a cool expanding H i shell. The shell will remain distinct from its surrounding until it has slowed down to a velocity of $`10\mathrm{km}\mathrm{s}^1`$. This would happen within $`2.3E_{51}^{0.32}10^6\mathrm{yr}`$ where $`E_{51}`$ is the initial energy in units of $`10^{51}\mathrm{erg}`$ (Loeb & Perna, 1998). The current rate of GRBs is one per $`10^7\mathrm{yr}`$ per galaxy (Schmidt, 1999). This leads us to expect a few remnants per galaxy at any given time. However SN are $`10^5`$ times more frequent and it would be difficult to distinguish between a SNR and a GRB remnant. Today there is mounting evidence that some GRBs are beamed (Halpern et al., 1999; Harrison et al., 1999; Kulkarni et al., 1999; Sari, 1999; Kuulkers et al., 2000; Salmonson, 2000). Beamed GRBs will illuminate only a fraction $`f_b`$ of the sky, and their rate should be higher by a factor of $`f_b^1`$. With $`f_b0.01`$ we would expect a hundred GRB remnants per galaxy. At an early stage the morphology of a “jetted” GRB remnant would be very different from a spherical explosion. This can be used to identify those remnants. We study this phase here.

Both beamed GRBs and SN deposit a comparable ($`10^{51}\mathrm{erg}`$) energy into the ISM. In both cases the evolution is expected to be similar since both are in the Sedov (Sedov, 1959) regime where all the energy is the initial explosion energy and all the mass is in the surrounding ISM. A key distinguishing feature unique to GRB remnants could be their beamed nature which we expect would lead to a distinct double shell morphology at intermediate times (the late time behavior of the GRB remnant is expected to be spherical in any case). In order to establish how many H i shells are GRB remnants we must find out the expected morphology of GRB remnants and how long they remain non-spherical, and distinguishable from SNRs.

Expanding H i shells have been found in many spiral galaxies (Tenorio-Tagle & Bodenheimer, 1988). The interior of these shells is relatively empty and their current expansion velocity is in the order of tens of Km s<sup>-1</sup>. Models for the origin of these H i supershels involve a large number of spatially correlated supernova events (Heiles, 1979, 1984), infall of massive gas clouds on to the galactic plane (Tenorio-Tagle & Bodenheimer, 1988) and flaring of radio lobes formed by jets ejected from the galactic nucleus during an active phase (Gopal-Krishna & Irwin, 2000). It has been previously noted (Efremov et al., 1998; Loeb & Perna, 1998) that a subset of these H i supershells may be the late signature left by GRBs on the interstellar medium (ISM). Establishing how many of the H i shells are GRB remnants would make it possible to directly estimate the local rate of GRBs, determine $`ϵ`$, the efficiency of converting the explosion energy into $`\gamma `$-rays, and the beaming factor $`f_b`$ (Loeb & Perna, 1998).

We model the intermediate evolution of a beamed GRB by two blobs of dense material moving into the ISM in opposite directions and we follow numerically their hydrodynamical evolution. Our results can be rescaled to fit a variety of initial energies. We find that at a time of $`5000\mathrm{yr}`$, when the ratio $`Mc^2/E_0`$ between the accumulated mass, $`M`$, and the initial GRB energy,$`E_0`$, reaches $`9.6\times 10^4`$ the remnant becomes spherical, similar in shape to a SNR. This value is independent of the exact initial details of our model such as the opening angle, the velocity and the morphology. We compare our results to the H i shells DEM L 316 (Williams et al., 1997) previously classified as two colliding SNRs. As the accumulated mass there is much larger it is most likely not a GRB remnant.

We describe our model and the numerical methods in section 2 . In section 3 we describe our results. Discussion and summary including a comparison with DEM L 316 are given in section 4.

## 2 The Simulation

### 2.1 Model

We study the intermediate evolution of the relativistic ejecta that caused a GRB. We assume that this matter was ejected in the form of two ultra relativistic blobs moving in opposite directions with a bulk Lorenz factor $`\mathrm{\Gamma }100`$. The emission from internal collisions in the blobs during an early stage comprises the GRB. Late external shocks caused by collisions with circumstellar matter produce the afterglow. The matter slows down during this interaction and its bulk Lorenz factor $`\mathrm{\Gamma }`$, decreases. The ejecta stays collimated only until $`\mathrm{\Gamma }`$ drops below $`1/\theta _0`$, at approximately $`3.4(E_{51,\mathrm{iso}}/\rho _\mathrm{I})^{1/3}(\theta _0/0.1)^{8/3}\mathrm{hr}`$ after the GRB, when it has accumulated $`10^5E_{51,\mathrm{iso}}(\theta _0/0.1)^2\mathrm{M}_{}`$ of ISM mass (Sari et al., 1999) where $`\theta _0`$ is the initial angular width, $`E_{51,\mathrm{iso}}`$ the isotropic energy in units of $`10^{51}\mathrm{erg}`$ and $`\rho _\mathrm{I}`$ is the surrounding ISM density in units of $`10^{24}\mathrm{gm}\mathrm{cm}^3`$. Note that in the rest of the paper $`E_{51}`$ denotes the actual energy in units of $`10^{51}\mathrm{erg}`$. At this time the matter will start to expand sideways causing, for an adiabatic evolution, an exponential slowing down (Rhoads, 1997). The ejecta continues to expand sideways at an almost constant radial distance from the source $`R_00.3E_{51}^{1/3}\rho _\mathrm{I}^{1/3}\mathrm{pc}`$ until it becomes non-relativistic. At this stage, we begin our simulation.

Without a detailed numerical modeling of the relativistic phase of the ejecta we have only an approximate description of how to construct the initial conditions. We expect the angular width of the ejecta to be $`1\mathrm{rad}`$. Additionally we are constrained by Newtonian energy conservation which yields a relation between $`R_0`$, $`\mathrm{\Gamma }`$ and $`E_0`$:

$$R_00.3E_{51}^{1/3}\rho _\mathrm{I}^{1/3}(v_0/c)^{2/3}\mathrm{pc}.$$
(1)
The exact shape and energy distribution of the ejecta are uncertain. This may influence the numerical coefficient in Eq. (1). However we show that the intermediate and late evolution of the ejecta are insensitive to these uncertainties in the initial conditions. This is due to the fact that since we are in the Sedov regime, the mass is dominated by the “external” ISM gas which washes out any variations in the initial conditions of the ejecta. This same mechanism will smear out any non-spherical features in a SN explosion, producing spherical SNRs.

We realize these initial conditions by two identical blobs moving in opposing directions into the ISM. The computational space is a cylinder in which the blobs move along the symmetry axis. Both the blobs and the ISM are modeled by a cold $`\gamma =5/3`$ ideal gas. The blobs are denser than the ISM. We consider various initial densities, angular widths and shapes of the blobs in order to make sure that these are indeed unimportant in determining the final morphology. In all runs $`v_0c/3`$. To evolve these initial conditions we use a hydrodynamic code, neglecting radiative effects.

We note that we have also used the post Newtonian version of our code with the same initial blob velocity of $`c/3`$. We found that the blobs slow down rapidly and the results were similar to the Newtonian results. Within a few timesteps we have a Newtonian system. In essence this is equivalent to slightly changing the initial shape of the blobs which, as noted before, is insignificant due to the mass dominance of the ISM.

### 2.2 The Numerical Method

The code we use is based on the Newtonian version of the smooth particle hydrodynamics (SPH) code introduced in Ayal et al. (2001). The code was adapted for the specific problem at hand. The three computationally important features of the problem are the negligible role of gravity, the large amount of stationary gas representing the ISM, and the symmetry. The fact that gravity is unimportant has made the equations very simple and greatly increased the speed of the code. Implementing single particle time steps allowed us to put in a large volume of stationary gas without making a big impact on the computational time. This way most of the computational effort is invested in the “important” interacting particles. Typically only 0.1% of the particles actually move during each timestep.

The initial conditions are set up so that the blobs are on the $`z`$ axis with a velocity along the $`z`$ axis. With these initial conditions there is a rotational symmetry about the $`z`$ axis and a reflection symmetry about the $`xy`$ plane. We implement the reflection symmetry exactly and the rotation symmetry only approximately. The implementation, similar to Libersky et al. (1993), consists of considering only $`1/8`$ (one quadrant) of the computational volume and adding reflecting boundary conditions on the three inner boundaries defined by the $`xy`$, $`xz`$, and $`yz`$ planes. We implement no outer boundary condition since the surrounding gas is very cold and there is almost no pressure so the outer particles gain only a negligible velocity throughout the course of the simulation. The quadrant we evolve consists of 62,800 SPH particles of which only $`10`$ represent the initial blob, again stressing the total mass dominance of the ISM.

The reflecting boundary conditions are implemented using pseudo-particles. At the beginning of each timesteps, all particles which intersect one of the boundaries (In SPH each particle has a finite size called the smoothing length) are reflected about this boundary and added to the simulation as additional pseudo-particles. After this is done for all reflecting boundaries, we calculate all the time derivatives in the usual manner treating all particles in an equal way. When all time derivatives are calculated we delete all the pseudo-particles. Only the “real” particles are then advanced in time. This simple algorithm allows us to implement reflecting boundary conditions using only a small number (typically 10% of the total number of particles for the three boundaries we use) of additional particles.

The use of SPH requires adding some artificial viscosity in order to resolve shocks. We use the standard artificial viscosity (e.g. Monaghan, 1992; Benz, 1990) consisting of a term analogous to bulk viscosity and a Von Neuman-Richtmyer artificial viscosity term. For the time integration we used a second order Runge-Kutta integrator with an adaptive stepsize control.

## 3 Results

We choose the initial velocity to be in the range $`0.22c`$ to $`0.33c`$, making relativistic effects small. Equation (1) leaves the freedom of choosing two out of the three parameters $`E_0`$, $`R_0`$ and $`\rho _\mathrm{I}`$, the initial energy, distance and ISM density respectively. In presenting the results we choose $`E_0`$ and $`\rho _0`$. To parameterize the evolution of the remnant we utilize the fact that mass scales linearly with the initial energy and we define the dimensionless parameter $`\mu =Mc^2/E_0`$ where $`M`$ is the accumulated shell mass. We define $`M`$ as all mass with density above $`2\rho _0`$ corresponding to the shocked material of the shell and to the accumulated ISM mass. All subsequent results are presented as functions of $`\mu `$. Our simulation begins approximately 2 years after the burst, when $`\mu 18`$. In Figure 1, we show the linear scaling of time with $`\mu `$. The scaling relation is $`t0.053\mu (E_{51}/\rho _\mathrm{I})^{1/3}\mathrm{yr}`$ which is close to the Sedov-Taylor (Sedov, 1959) blast wave result of $`tR^{5/2}\mu ^{5/6}`$ where $`R`$ is the radius of the blast wave.

A bow shock forms as each blob collides with the ISM. The shock also propagates in the direction perpendicular to the blob’s velocity and over time, backwards. The expected morphology of the remnant will therefore be of two expanding shells which will eventually join, producing yet another shock. At late times the shells will merge and become a single spherical shell. We made several runs, differing in the initial density, shape, initial velocity, angular width and the numerical coefficient in Eq. (1) of the ejecta as summarized in Table 1. In five runs the initial ejecta was shaped like a disk and in one case it was shaped like a sphere. As we show in Fig. 2, for $`\mu 5\times 10^3`$ the $`z`$ positions of the shock and its maximal radius in the $`xy`$ plane, $`r_{xy}`$ are indeed unaffected by these variations in initial conditions. The $`\mu `$ at which the shells touch depends on the exact initial conditions, but all other events (such as the $`\mu `$ at which the remnant becomes spherical) are unaffected. Therefore in the following discussions we will consider only one of the runs.

The mass inside the shell, defined as all mass with density below $`\rho _0/2`$ also evolves almost linearly with $`\mu `$ (Fig. 3) as $`0.06(E_0/c^2)\mu `$. In Fig. 4, we show the density contours as a function of $`\mu `$. The two shells touch and a shock forms along the equatorial plane between them at $`t50260(E_{51}/\rho _\mathrm{I})^{1/3}\mathrm{yr}`$ when $`\mu =15\times 10^3`$ depending on the initial conditions. The maximal $`z`$ position of the shock and $`r_{xy}`$ can be fitted with a power law as shown in Fig. 5. The scaling relation for $`r_{xy}`$ is $`r_{xy}=0.07\mu ^{0.45}(E_{51}/\rho _\mathrm{I})^{1/3}\mathrm{pc}t^{0.45}`$. The effective radius of the shell, $`(r_{xy}^2z)^{1/3}`$, scales as $`t^{0.4}`$, exactly the result expected for a Sedov-Taylor blast wave. The ratio $`z_{max}/r_{xy}`$ decreases in time staying always between 1 and 2 (see Fig 6). Extrapolating we see that this ratio reaches a value of 1 at $`\mu 9.6\times 10^4`$. At this time the shock has a spherical shape with $`z=r_{xy}12(E_{51}/\rho _\mathrm{I})^{1/3}\mathrm{pc}`$. Even then the shock will not be completely spherically symmetric as there would still be a ring of shocks around the “equator” where the shells have collided.

During the whole run, total energy is conserved to within 1%. In Fig. 7 we show $`E_z`$ the total kinetic energy in the $`z`$ direction, $`E_{xy}`$ the total kinetic energy in the $`xy`$ plane and $`E_i`$ the total internal energy in the simulated volume. As the blobs interact with the ISM $`E_i`$ and $`E_{xy}`$ increase at the expense of $`E_z`$. $`E_{xy}`$ increases until $`\mu 5\times 10^3`$ and then it remains constant at $`0.22E_0`$. $`E_i`$ increases by a factor of 200 to a value of $`0.72E_0`$, most of the increase occurs before $`\mu 2500`$. In the final configuration ($`\mu 2\times 10^4`$) 72% of the energy is in internal energy, 22% is in kinetic energy in the $`xy`$ plane and only 6% remains in the kinetic energy in the $`z`$ direction.

Figures 8 and 9 depict the images of the remnant as a function of time and angles of inclination. We show images due to bremsstrahlung emission and synchrotron emission. The images are constructed assuming that all the gas is optically thin in the relevant frequencies. The bremsstrahlung luminosity (Fig. 8) was calculated assuming that the volume emissivity is proportional to $`\rho ^2\epsilon ^{1/2}`$ (Lang, 1980). In calculating the synchrotron emissivity (Fig. 9) we assumed that both the magnetic field energy density and the energy of the relativistic electrons are proportional to the internal energy density of the gas with the proportionality factors $`ϵ_B`$ and $`ϵ_e`$ respectively. We further assume that the relativistic electron number density is a power law in energy. Under these assumptions the volume emissivity is proportional to $`\rho ^2\epsilon ^2`$ (e.g. Shu, 1991). In the late images there are two bright circles at the lines where the colliding blobs forming a hot shocked region.

## 4 Discussion

We follow the hydrodynamic evolution of two blobs colliding with the ISM. This scenario is a model for the intermediate time behavior of matter ejected from a central engine producing a GRB. Our results can be rescaled to fit a variety of initial energies and ISM densities. The two shells touch and a shock forms along the equatorial plane between them at $`t50160(E_{51}/\rho _\mathrm{I})^{1/3}\mathrm{yr}`$ when $`\mu =15\times 10^3`$. We show that the late time remnant is insensitive to the exact initial morphology, angular width and density of the ejecta. Although initially the remnant may be highly non-spherical, the ratio between its height and radius will approach unity and it will eventually become spherical in shape after a time of $`5\times 10^3E_{51}^{1/3}\rho _\mathrm{I}^{1/3}\mathrm{yr}`$ when $`\mu 9.6\times 10^4`$. After this time it will be difficult to distinguish a GRB remnant from a SNR using the morphology. The expected number of non-spherical GRB remnants is, therefore, $`5\times 10^4f_b^1E_{51}^{1/3}\rho _\mathrm{I}^{1/3}`$ per galaxy (using the observed present GRB rate (Schmidt, 1999) of $`10^7\mathrm{yr}^1\mathrm{gal}^1`$).

The results of Gaensler (1998) show a tendency for the bilateral axis of the non-spherical SNRs to be aligned with the galactic plane. This presents strong evidence in favor of an extrinsic model for the origin of non-spherical SNRs. Using our results, we can propose an alternative explanation to the origin of some of the highly non-spherical SNRs and H i shells. Instead of assuming an extrinsic model, namely spherical energy deposition into a non isotropic medium we propose an intrinsic model: non spherical energy deposition into an isotropic medium. Our GRB remnant model can explain non-spherical SNRs with a two shell morphology provided that $`\mu 9.6\times 10^4`$. One extreme and well studied case, the non-spherical SNR DEM L 316 (Williams et al., 1997), has a distinct double shell morphology. The external model clearly fails here. However the ratio $`\mu `$ measured for DEM L 316 is $`7.1\times 10^5`$ and a a GRB remnant would already be spherical at this stage, suggesting that DEM L 316 is not a GRB remnant. This conclusion could be revised if for some reason the initial blobs stay confined for a much longer period (so that eq. 1 is not satisfied) or if the solution becomes radiative before becoming spherical and the radiative cooling slow down the evolution.

We thank the anonymous referee for his helpful comments. This research was supported by a ISRAEL-US BSF grant.


format: text

text:
# Phase Transition of a Heisenberg Spin-Glass Model in Three Dimensions

## Abstract

We study the phase transition of the $`\pm J`$ Heisenberg model in three dimensions. Using a dynamical simulation method that removes a drift of the system, the existence of the spin-glass (SG) phase at low temperatures is suggested. The transition temperature is estimated to be $`T_{\mathrm{SG}}0.18J`$ from both equilibrium and off-equilibrium Monte-Carlo simulations. Our result contradicts the chirality mechanism of the phase transition reported recently by Kawamura which claims that it is not the spins but the chiralities of the spins that are ordered in Heisenberg SG systems.

It has been believed that the spin-glass (SG) phase is realized in three dimensions (3D) for the Ising model  but not for the XY and Heisenberg models. Thus, the SG phases observed in experiments were suggested to be realized due to magnetic anisotropies. However, it remains puzzling that no sign of Heisenberg-Ising crossover has been detected, which is expected if the observed finite temperature phase transition is due to a weak magnetic anisotropy. Recently, Kawamura and coworkers have proposed a chirality mechanism to solve the puzzle. They considered the chirality described by neighboring three spins. It has either positive or negative value, just like the Ising spin. Kawamura calculated the chirality autocorrelation function $`C_\chi (t)`$ and the spin autocorrelation function $`C_S(t)`$ of 3D Heisenberg SG models with and without a weak random magnetic anisotropy $`D`$ and found that, even for $`D=0`$, $`C_\chi (t)`$ exhibits a pronounced aging effect reminiscent of the one observed in the mean-field model, while $`C_S(t)`$ exhibits a similar aging effect only when $`D0`$. He claimed that, in the 3D Heisenberg SG model, a chiral-glass (CG) phase transition occurs at a finite temperature $`T_{\mathrm{CG}}0`$, but the SG phase is absent. He argued that, in real SG magnets, the spin and the chirality are mixed due to a weak magnetic anisotropy, and the CG transition is revealed as anomalies in experimentally accessible quantities. According to this interpretation, the SG phase transition never occurs in Heisenberg-like systems, and what was observed in experiments is nothing but the CG phase transition. This mechanism is quite interesting, because it calls for re-consideration of the SG phase transition from both theoretical and experimental points of view, because most of SG magnets are rather Heisenberg-like.

In this Letter, we demonstrate that the apparent difference between $`C_S(t)`$ and $`C_\chi (t)`$ comes from a drift of the system and suggest that, contrary to the chirality mechanism, the SG phase is realized in the Heisenberg SG model even when the anisotropy is absent. We estimate the SG transition temperature of the $`\pm J`$ Heisenberg model as $`T_{\mathrm{SG}}0.18J`$ from both equilibrium and off-equilibrium Monte-Carlo simulations. We believe that what was revealed by Kawamura and coworkers is nothing but the ordering of the spins themselves. The most important point revealed here is that the SG phase will be much more stable than what were believed so far and re-examination is also necessary for different models.

We start with the $`\pm J`$ Heisenberg model on a simple cubic lattice of $`L\times L\times (L+1)(N)`$ described by

$`H={\displaystyle \underset{ij}{}}J_{ij}𝑺_i𝑺_j,`$ (1)
where $`𝑺_i`$ is the Heisenberg spin of $`|𝑺_i|=1`$ and $`ij`$ runs over all nearest-neighbor pairs. The exchange interaction (bond) $`J_{ij}`$ takes on either $`+J`$ or $`J`$ with the same probability of 1/2.

First we consider the dynamical properties of the system calculating the spin autocorrelation function. When one considers the dynamics of the Heisenberg model with $`D=0`$, one should pay a special attention to a drift of the system, i.e., the global rotation of the system. Although this drift does not affect the static spin correlation and becomes negligible in large systems, some data of finite systems inevitably suffer from this effect. The spin autocorrelation $`C_S(t)`$ diminishes by this drift, whereas the chirality autocorrelation $`C_\chi (t)`$ does not, because the chiralities are invariant under the global rotation of the system. Therefore the difference in the behavior between $`C_S(t)`$ and $`C_\chi (t)`$ might come from this drift.

In order to consider the dynamical properties of an infinitely large system, we remove this effect applying the uniform rotation to the system so that

$`S={\displaystyle \underset{i}{}}|R(t)𝑺_i(t+t_w)𝑺_i(t_w)|^2`$ (2)
becomes minimum, where $`R(t)`$ is the rotation matrix, and $`\{𝑺_i(t_w)\}`$ and $`\{𝑺_i(t+t_w)\}`$ are the spin configurations at times $`t_w`$ and $`t+t_w`$, respectively. We can successively determine $`R(t)`$, because $`R(t)=E`$ at $`t=0`$ and it changes step by step as $`t`$ goes by, where $`E`$ is the unit matrix. Hereafter, we call the system described by $`\{𝑺_i(t)\}`$ the drifting system and that by $`\{\stackrel{~}{𝑺}_i(t+t_w)R(t)𝑺_i(t+t_w)\}`$ the fixed system. The spin autocorrelation function $`C_S(t_w,t+t_w)`$ of the fixed system is defined as

$`C_S(t_w,t+t_w)={\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}[\stackrel{~}{𝑺}_i(t+t_w)\stackrel{~}{𝑺}_i(t_w)],`$ (3)
where $`\mathrm{}`$ and \[$`\mathrm{}`$\] mean the thermal average and the bond distribution average, respectively. The chirality autocorrelation function is also defined as

$`C_\chi (t_w,t+t_w)={\displaystyle \frac{1}{3N}}{\displaystyle \underset{i\mu }{}}[\chi _{i\mu }(t+t_w)\chi _{i\mu }(t_w)],`$ (4)
where $`\chi _{i\mu }`$ is the chirality at the $`i`$th site and in the $`\mu `$th direction defined by $`\chi _{i\mu }=𝑺_{i+\widehat{𝒆}_\mu }(𝑺_i\times 𝑺_{i\widehat{𝒆}_\mu })`$ with $`\widehat{𝒆}_\mu (\mu =x,y,z)`$ being a unit lattice vector along the $`\mu `$ axis. Note again that the chirality $`\chi _{i\mu }`$ is invariant under the global spin rotation, and then the chirality autocorrelation functions of both drifting and fixed systems are completely the same and described by Eq. (4). Before applying this method to the present model, we test it in the ferromagnetic Heisenberg model, results of which are shown in Fig. 1. It is found that the spin autocorrelation function of the fixed system with a large $`t_w`$ approaches to the equilibrium value of the square of the order parameter, $`(𝑴/N)^2`$, with $`𝑴`$ being the total magnetization, whereas that of the drifting system rapidly decays with increasing $`t`$. This result indicates that the effect of the drift can be removed properly by the present method.

We perform the same simulation that was done by Kawamura using the standard single spin-flip heat-bath Monte-Carlo(MC) method. That is, starting with a random initial spin configuration, the system is quenched to a working temperature, and after waiting $`t_w`$ MC steps per spin (MCS) the autocorrelation functions are measured up to about $`2\times 10^5`$ MCS. A sample average is taken over about 32 independent bond distributions. The results of $`C_S(t_w,t+t_w)`$ for different $`t_w`$ are presented in Fig. 2 as functions of $`t`$. Contrary to that reported by Kawamura, $`C_S(t_w,t+t_w)`$ exhibits a quite similar aging effect that was found for $`C_\chi (t_w,t+t_w)`$ (see Fig. 1 in ref.), which indicates that the spins and the chiralities possess similar dynamical properties. To examine this point in more detail, we plot in Fig. 3 the ratio $`T(t_w,t+t_w)C_\chi /C_S`$ for a fixed $`t_w=10^5`$ at different temperatures. At all the temperatures, $`T(t_w,t+t_w)`$ slowly decreases as $`t`$ goes by, revealing that $`C_S`$ decays more slowly than $`C_\chi `$. This result clearly indicates that, if the ordering of the chiralities is realized, then the same is true for the spin. We estimate values of the Edward-Anderson order parameters of the spin and the chirality, $`q_{\mathrm{SG}}`$ and $`q_{\mathrm{CG}}`$, using the method proposed by Parisi et al. and used by Kawamura. That is, these values are extracted by fitting the data of $`C_S`$ and $`C_\chi `$ for $`t_w=3\times 10^5`$ to the power-law form of

$`C_{S,\chi }(t_w,t+t_w)q_{\mathrm{SG},\mathrm{CG}}+{\displaystyle \frac{C}{t^\lambda }},`$ (5)
in the time ranges of $`30t3000`$ and $`50t5000`$ . The obtained $`q_{\mathrm{SG}}`$ and $`q_{\mathrm{CG}}`$ are plotted as functions of temperature in Fig. 4. Both quantities have positive, non-vanishing values at low temperatures and seem to vanish at almost the same temperature of $`T0.18J`$. This result suggests that, in fact, the phase transitions of the spin and the chirality exist at $`T0`$ and they occur at the same temperature. It is not easy, however, to estimate the definite value of the transition temperature as well as values of the order parameter exponents, because the data of $`q_{\mathrm{SG}}`$ around $`T=0.18J`$ depend strongly on the time range used for extracting them.

Next we perform a detailed MC simulation in the thermal equilibrium to examine the occurrence of the SG phase transition itself. We calculate the SG susceptibility $`\chi _{\mathrm{SG}}`$ defined by

$`\chi _{SG}={\displaystyle \frac{1}{N}}{\displaystyle \underset{ij}{}}[𝑺_i𝑺_j^2],`$ (6)
using the exchange MC algorithm. The sizes of the lattice studied here are $`L=519`$. Equilibration is checked by monitoring the stability of the results against at least two-times longer runs. The numbers of the samples are 320 for $`L11`$, 96 for $`L=15`$, and 48 for $`L=19`$. The results of $`\chi _{\mathrm{SG}}`$ for different $`L`$ are plotted in Fig. 5. At low temperatures, $`\chi _{\mathrm{SG}}`$ exhibits a strong size dependence implying the occurrence of the phase transition at $`T0`$. If the lower critical dimension $`d_l`$ is less than the lattice dimension, $`d_l<3`$, and the phase transition really occurs at $`T=T_{\mathrm{SG}}`$, the data for different $`L`$ will be scaled as

$`\chi _{SG}=L^{2\eta }F(L^{1/\nu }(TT_{\mathrm{SG}})),`$ (7)
where $`\nu `$ is the exponent of the correlation length and $`\eta `$ is the exponent which describes the decay of the correlation function at $`T=T_{\mathrm{SG}}`$. The scaling plots obtained by assuming $`T_{\mathrm{SG}}0`$ and $`T_{\mathrm{SG}}=0`$ are shown in Fig. 6. It should be evident the scaling with $`T_{\mathrm{SG}}0`$ works much better than that with $`T_{\mathrm{SG}}=0`$, even if the data for the smallest size $`L=5`$ are ignored in the latter. The phase transition temperature and the values of the critical exponent estimated here are $`T_{\mathrm{SG}}/J=0.18\pm 0.01`$, $`\nu =0.97\pm 0.05`$ and $`\eta =0.1\pm 0.1`$. Note that, in the $`\pm J`$ model with the random magnetic anisotropy $`D0`$, the transition temperature was estimated as $`T_{\mathrm{SG}}0.30J`$ for $`D=0.1J`$, and $`T_{\mathrm{SG}}0.25J`$ for $`D=0.05J`$. Thus our value of $`T_{\mathrm{SG}}0.18J`$ is consistent with the one obtained by extrapolation of those values. The value of $`\nu 1.0`$ for the correlation length also is not strongly different from that of $`\nu 1.2`$ estimated in the case of $`D0`$.

In summary, all the results presented here seem to support the view that the SG phase is realized in the $`\pm J`$ Heisenberg model and the phase transition to this phase occurs at $`T_{\mathrm{SG}}0.18J`$. This view is compatible with our recent studies of the stiffness of the system, which suggest that, contrary to previous studies, the stiffness exponent $`\theta _\mathrm{S}`$ of the spin has the positive value of $`\theta _\mathrm{S}0.8`$ at $`T=0`$ and changes its sign at $`T0.19J`$. Thus we believe that the SG phase is realized in this model even when the anisotropy is absent. A major role of the anisotropy will be neither to bring the phase transition nor to mix the spin and the chirality, but to fix the system and lead anomaly of the susceptibility in a finite system. One might think that the present suggestion is quite surprising, because it disproves the common belief. It should be noted again that the occurrence of the ordered phase of the chiralities was already suggested by Kawamura et al. and what is suggested here is that the spins themselves also order. We think the origin of the occurrence of the phase transition in the Heisenberg SG model is the ordering of the spins and that the ordering of the chiralities happens to be accompanied by it. However, further studies are necessary to establish this concept as well as to reveal the nature of the phase transition. We hope that the present paper stimulates studies of SG systems with the continuous spin symmetry.

The authors would like to thank Professor K. Sasaki, Dr. T. Nakamura, Professor S. Miyake, and Professor H. Takayama for their valuable discussions.


text:
# Low energy scattering and photoproduction of 𝜂-mesons on deuterons.

## 1 Introduction

The production of $`\eta `$ mesons and their collisions with nuclei have been studied experimentally and theoretically with increasing interest during the last years. To a large extent this is motivated by the fundamental problem of charge–symmetry breaking of strong interaction . Another relevant questions concerns the nature of $`S_{11}(1535)`$–resonance  and the possible formation of $`\eta `$–nucleus quasi–bound states . It is worth mentioning that according to  the mean free path of $`\eta `$ mesons in a nuclear medium is about 2 fm, i.e., less than the size of a typical nucleus. A necessary condition for the existence of $`\eta `$-nuclei, hence, appears to be satisfied.

For the calculation of these states various model treatments were employed, among them the optical potential method , the Green’s function method , the modified multiple scattering theory  and few–body calculations . The predictions concerning the possibility of $`\eta `$–mesic nucleus formation are very diverse. One obvious reason for such a diversity is the poor knowledge of the $`\eta N`$ forces. Another reason comes from the differences among the employed approximations some of which might be faulty in view of the resonant character of the $`\eta N`$ dynamics and the delicacy of the quasi–bound state problem.

In the present paper we treat the $`\eta `$–deuteron system on the basis of the exact few–body Alt–Grassberger–Sandhas (AGS) equations. The Faddeev–type coupling of these equations guarantees uniqueness of their solutions. Moreover, as equations for the elastic and rearrangement operators they are well–defined in momentum space, providing thus the desired scattering amplitudes in a most direct and technically reliable manner. The advantage of working with coupled equations involving the elastic and rearrangement operators of the final state, is not only suggested by questions of uniqueness, but also by the relevance of rescattering effects which were found to give a significant contribution to the corresponding amplitude .

## 2 AGS formalism

In terms of the AGS transition operator $`U_{11}`$ the $`\eta d`$ elastic scattering amplitude is represented as

$$f(𝐩_1^{},𝐩_1;z)=(2\pi )^2M_1𝐩_1^{};\psi _d|U_{11}(z)|𝐩_1;\psi _d$$
(1)
with the on–energy–shell conditions $`|𝐩_1^{}|=|𝐩_1|`$ and $`z=p_1^2/2M_1+E_d`$ with $`E_d`$ being the deuteron energy. Here the subscript 1 labels the $`\eta (NN)`$ partition and the $`\eta `$–deuteron channel. The transition operator $`U_{11}`$ obeys the system of AGS equations

$$U_{\beta \alpha }(z)=(1\delta _{\beta \alpha })G_0^1(z)+\underset{\gamma =1}{\overset{3}{}}(1\delta _{\beta \gamma })T_\gamma (z)G_0(z)U_{\gamma \alpha }(z),$$
(2)
with $`G_0(z)`$ being the free resolvent (Green’s operator) of the three particles involved. This set of equations couples all $`3\times 3`$ elastic and rearrangement operators $`U_{\alpha \alpha }`$ and $`U_{\beta \alpha }`$. Here each of the subscripts runs through the values 1, 2 and 3, indicating the two-fragment partitions (1,23), (2,31) and (3,12) respectively. Therefore, $`U_{11}`$ describes the elastic transition $`1(23)1(23)`$, while $`U_{21}`$ represents the rearrangement process $`1(23)2(13)`$. These subscripts are also used in the complementary notation to label the two–body T–operator $`t_\alpha (z)`$ of the ($`\beta \gamma `$) pair, for instance $`t_1(z)=t_{NN}(z)`$. It should be noticed, however, that it is not this genuine two–body operator which enters the AGS equations, but the operator

$$T_\alpha (z)=t_\alpha (z𝐪_\alpha ^2/2M_\alpha ),$$
(3)
which is to be understood as the two–body operator embedded in the three–body space, with the relative kinetic energy operator $`𝐪_\alpha ^2/2M_\alpha `$ of particle $`\alpha `$ being subtracted from the total energy variable z. Considered in momentum space, (3) thus reads

$$𝐩_{}^{}{}_{\alpha }{}^{},𝐪_{}^{}{}_{\alpha }{}^{}|T_\alpha (z)|𝐩_\alpha ,𝐪_\alpha =\delta (𝐪_{}^{}{}_{\alpha }{}^{}𝐪_\alpha )𝐩_{}^{}{}_{\alpha }{}^{}|t_\alpha (zq_\alpha ^2/2M_\alpha )|𝐩_\alpha ,$$
(4)
where $`𝐩_\alpha `$ and $`𝐪_\alpha `$ are the Jacobi momenta of the pair $`\alpha `$ and the spectator $`\alpha `$ respectively, and $`M_\alpha `$ is the corresponding reduced mass.

For both $`T_{\eta N}`$ and $`T_{NN}`$ we used one–term separable forms

$$T_\alpha (z)=|\chi _\alpha \tau _\alpha (z)\chi _\alpha |.$$
(5)
For the $`NN`$ subsystem Eq. (5) implies that the asymptotic wave function is related to the form-factor $`|\chi _1`$ according to

$$|𝐩_1;\psi _d=G_0(z)|\chi _1|𝐩_1.$$
(6)
Due to (5) and (6) the scattering amplitude (1) can be rewritten as

$$f(𝐩_1^{},𝐩_1;z)=(2\pi )^2M_1𝐩_1^{}|X_{11}(z)|𝐩_1,$$
(7)
where the operators $`X_{\beta \alpha }`$, defined as

$$X_{\beta \alpha }(z)=\chi _\beta |G_0(z)U_{\beta \alpha }(z)G_0(z)|\chi _\alpha ,$$

obey the system of equations

$$X_{\beta \alpha }(z)=Z_{\beta \alpha }(z)+\underset{\gamma =1}{\overset{3}{}}Z_{\beta \gamma }(z)\tau _\gamma \left(z\frac{p_\gamma ^2}{2M_\gamma }\right)X_{\gamma \alpha }(z)$$
(8)
with

$$Z_{\beta \alpha }(z)=(1\delta _{\beta \alpha })\chi _\beta |G_0(z)|\chi _\alpha .$$

The identity of the nucleons implies that $`X_{31}=X_{21}`$, $`\tau _3=\tau _2`$, and $`Z_{31}=Z_{21}`$, which reduces the system (8) to two coupled equations

$$\begin{array}{ccc}\hfill X_{11}(z)& =& 2Z_{12}(z)\tau _2\left(z\frac{p_2^2}{2M_2}\right)X_{21}(z),\hfill \\ & & \\ \hfill X_{21}(z)& =& Z_{21}(z)+Z_{21}(z)\tau _1\left(z\frac{p_1^2}{2M_1}\right)X_{11}(z)+Z_{23}(z)\tau _2\left(z\frac{p_2^2}{2M_2}\right)X_{21}(z).\hfill \end{array}$$
(9)
Eventually, after making the $`S`$–wave projection of the matrix elements $`𝐩_\beta ^{}|X_{\beta \alpha }|𝐩_\alpha `$ and $`𝐩_\beta ^{}|Z_{\beta \alpha }|𝐩_\alpha `$, we end up with one–dimensional integral equations which can be solved numerically.

## 3 Two–body T–matrices

The $`S`$-wave nucleon-nucleon separable potential is adopted from Ref. with its parameters slightly modified to be consistent with more recent $`NN`$ data (see Ref.). The $`\eta `$-nucleon $`T`$-matrix is taken in the form

$$t_{\eta N}(p^{},p;z)=(p_{}^{}{}_{}{}^{2}+\alpha ^2)^1\frac{\lambda }{(zE_0+i\mathrm{\Gamma }/2)}(p^2+\alpha ^2)^1$$
(10)
consisting of two vertex functions and the $`S_{11}`$-propagator in between . It corresponds to the process $`\eta NS_{11}\eta N`$ which at low energies is dominant. The range parameter $`\alpha =3.316`$ fm-1 was determined in Ref. , while $`E_0`$ and $`\mathrm{\Gamma }`$ are the parameters of the $`S_{11}`$ resonance ,

$$E_0=1535\mathrm{MeV}(m_N+m_\eta ),\mathrm{\Gamma }=150\mathrm{MeV}.$$

The strength parameter $`\lambda `$ is chosen to reproduce the $`\eta `$-nucleon scattering length $`a_{\eta N}`$,

$$\lambda =\frac{\alpha ^4(E_0i\mathrm{\Gamma }/2)}{(2\pi )^2\mu _{\eta N}}a_{\eta N}.$$
(11)
the imaginary part of which accounts for the flux losses into the $`\pi N`$ channel.

The two–body scattering length $`a_{\eta N}`$ is not accurately known. Different analyses  provided for $`a_{\eta N}`$ the values in the range

$$0.27\mathrm{fm}\mathrm{Re}a_{\eta N}0.98\mathrm{fm},0.19\mathrm{fm}\mathrm{Im}a_{\eta N}0.37\mathrm{fm}.$$
(12)
Recently, however, most of the authors agreed that Im $`a_{\eta N}`$ is around $`0.3`$ fm. But for Re $`a_{\eta N}`$ the estimates are still very different (compare, for example, Refs.  and ). Most of the results presented in this paper are, therefore, obtained using Im $`a_{\eta N}=0.3`$ fm and several values of Re $`a_{\eta N}`$ within the above interval.

We solved Eqs. (9) for $`\eta d`$ collision energies varying from zero ($`\eta d`$–threshold, $`z=E_d`$) up to $`22`$ MeV by replacing the integrals by Gaussian sums. As is well known (see, for example, ), the kernels of these equations, when expressed in momentum representation, have moving logarithmic singularities for $`z>0`$. In the numerical procedure, we handle it with the method suggested in Ref.. The main idea of this method consists in interpolating the unknown solutions (in the area covering the singular points) by certain polynomials and subsequent analytic integration of the singular part of the kernels.

## 4 Scattering results

The results of our calculations are presented in Figs. 24 and in Table 1. In Fig. 2 the energy dependence of the $`\eta d`$ phase-shifts for five different choices of Re $`a_{\eta N}`$ is shown. The curves correspond (starting from the lowest one) to Re $`a_{\eta N}=0.55`$ fm, $`0.65`$ fm, $`0.725`$ fm, $`0.75`$ fm, and $`0.85`$ fm. The larger this value, the stronger is the $`\eta `$N attraction.

The change in the character of these curves, hence, reflects the growth of the attractive force between the $`\eta `$ meson and the nucleon. The lower three curves for Re $`\delta _{\eta d}`$ corresponding to the smaller values of Re $`a_{\eta N}`$ start from zero, the two curves corresponding to the strong attraction start from $`\pi `$. According to Levinson’s theorem, the phase shift at threshold energy is equal to the number of bound states $`n`$ times $`\pi `$. We found that the transition from the lower family of the curves to the upper one happens at the critical value Re $`a_{\eta N}=0.733`$ fm. Therefore, the $`\eta `$N force, which generates Im $`a_{\eta N}=0.3`$ fm and Re $`a_{\eta N}>0.733`$ fm, is sufficiently attractive to bind $`\eta `$ inside the deuteron.

For three other choices of Im $`a_{\eta N}`$ within the uncertainty interval, namely, $`0.20`$ fm, $`0.25`$ fm, and $`0.35`$ fm, the corresponding critical values of Re $`a_{\eta N}`$ turned out to be $`0.788`$ fm, $`0.761`$ fm, and $`0.698`$ fm. In the complex $`a_{\eta N}`$–plane (see Fig. 2) the corresponding points form a curve separating the uncertainty area (given by formulae (12), dashed rectangle) into two parts. If $`a_{\eta N}`$ is to the right of this curve, the strength of $`\eta `$N attraction is sufficient for $`\eta d`$ bound state formation.

In Fig. 4 we present the result of our calculations of the total cross-section (integrated over the angles) for elastic $`\eta `$-deuteron scattering as a function of collision energy. The five curves correspond (starting from the lowest one) to Re $`a_{\eta N}=0.55`$ fm, $`0.65`$ fm, $`0.725`$ fm, $`0.75`$ fm, and $`0.85`$ fm. The dashed line indicates the deuteron break–up threshold. The peaks in the energy dependence of the total elastic cross–section indicate that a resonance appears in the $`\eta d`$–system. Of course, not every maximum of the cross–section is a resonance, but we plotted the Argand plots also (Fig. 4) for the $`\eta d`$ elastic scattering amplitude in the energy interval from $`0`$ to $`22`$ MeV. The five curves correspond (from right to left) to Re $`a_{\eta N}=0.55`$ fm, $`0.65`$ fm, $`0.725`$ fm, $`0.75`$ fm, and $`0.85`$ fm. When the energy increases the corresponding points move anticlockwise. So that the Argand plots prove that the maxima we found are resonances. Their positions and widths for various choices of Re $`a_{\eta N}`$ are given in Table 1. It should be noted that, while the resonance energy is determined in our calculations exactly (as the maximum of the function $`\mathrm{sin}^2\mathrm{Re}\delta _{\eta d}`$), the corresponding width is obtained by fitting the cross-section with a Breit-Wigner curve. Therefore, the values of $`\mathrm{\Gamma }_{\eta d}`$ given in Table 1 should be considered only as rough estimates.

The resonant behavior of $`\eta d`$ elastic scattering should be seen in various processes involving $`\eta d`$-system in their final states, such as $`\gamma d\eta d`$ and $`np\eta d`$. Indeed, the corresponding amplitudes $`\psi _{\mathrm{out}}|𝒪|\psi _{\mathrm{in}}`$ involve the $`\eta d`$ wave function $`\psi _{\mathrm{out}}`$ which, in the vicinity of the resonance, strongly depends on the total energy. Its resonant growth at short distances may enhance the transition probability. To check this suggestion we performed the calculations of the $`\gamma d\eta d`$ reaction.

## 5 Photoproduction of $`\eta `$-mesons

Theoretical analysis of $`(\gamma ,\eta )`$-reactions on nuclei is hampered by the three major problems: the unknown off-shell behavior of the two-body $`\gamma N\eta N`$ amplitude, inaccuracies in the description of the nuclear target as a many-body system, and rescattering effects in the final state. The simplest is, of course, the process of coherent $`\eta `$ photoproduction on deuteron. There are many theoretical studies devoted to ($`\gamma ,\eta `$) reactions on deuteron. Early attempts to go beyond a simple impulse approximation led to very different conclusions  as do more recent approaches based on the effective two-body formulations . Moreover, the experimental cross-section  of the reaction

$$\gamma +d\eta +d$$
(13)
in the near-threshold region is far above these theoretical predictions. Therefore, a reliable description of $`\eta `$ photoproduction on deuteron on the basis of exact equations is desirable.

## 6 Formalism

To consider reaction (13), we employ the exact Alt-Grassberger-Sandhas formalism modified to include the electromagnetic interaction. A photon can be introduced into this formalism by considering the $`\eta N`$ and $`\gamma N`$ states as two different channels of the same system. This means that we should replace the T-operator $`t_{\eta N}`$ by the $`2\times 2`$ matrix. It is clear, that such replacements of the kernels of integral equations (2) lead to the corresponding solutions having a similar matrix form

$$t_{\eta N}\left(\begin{array}{cc}t^{\gamma \gamma }& t^{\gamma \eta }\\ t^{\eta \gamma }& t^{\eta \eta }\end{array}\right)U_{\alpha \beta }\left(\begin{array}{cc}W_{\alpha \beta }^{\gamma \gamma }& W_{\alpha \beta }^{\gamma \eta }\\ W_{\alpha \beta }^{\eta \gamma }& W_{\alpha \beta }^{\eta \eta }\end{array}\right).$$
(14)
Here $`t^{\gamma \gamma }`$ describes the Compton scattering, $`t^{\eta \gamma }`$ the photoproduction process, and $`t^{\eta \eta }`$ the elastic $`\eta N`$ scattering111Another method of treatment electromagnetic process of the type (13) in the frame of AGS equations was used in Ref...

It is technically more convenient to consider the reaction of $`\eta `$-photoabsorption, which is inverse to reaction (13). Then the photoproduction cross-section can be obtained by applying the detailed balance principle. We are interested in the coherent process, therefore we need the amplitude $`W_{11}^{\gamma \eta }`$ obeying the equation

$`\left(\begin{array}{cc}W_{11}^{\gamma \gamma }& W_{11}^{\gamma \eta }\\ W_{11}^{\eta \gamma }& W_{11}^{\eta \eta }\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}T_2^{\gamma \gamma }& T_2^{\gamma \eta }\\ T_2^{\eta \gamma }& T_2^{\eta \eta }\end{array}\right)G_0\left(\begin{array}{cc}W_{21}^{\gamma \gamma }& W_{21}^{\gamma \eta }\\ W_{21}^{\eta \gamma }& W_{21}^{\eta \eta }\end{array}\right)`$ (26)
$`+\left(\begin{array}{cc}T_3^{\gamma \gamma }& T_3^{\gamma \eta }\\ T_3^{\eta \gamma }& T_3^{\eta \eta }\end{array}\right)G_0\left(\begin{array}{cc}W_{31}^{\gamma \gamma }& W_{31}^{\gamma \eta }\\ W_{31}^{\eta \gamma }& W_{31}^{\eta \eta }\end{array}\right).`$
In the first order on electromagnetic interaction:

$$W_{11}^{\gamma \eta }T_2^{\gamma \eta }G_0(z)W_{21}^{\eta \eta }+T_3^{\gamma \eta }G_0(z)W_{31}^{\eta \eta }$$
(27)
and we can see, that in above approximation the transition operator $`W_{11}^{\gamma \eta }`$ can be found not from an equation, but from the expression (27). It is physically clear since in the first order on electromagnetic interaction the only contributions are from the $`\eta `$–meson rescattering which is described by the $`W_{21}^{\eta \eta }`$ and $`W_{31}^{\eta \eta }`$ operators.

It was experimentally proven  that at low energies the reaction $`\gamma N\eta N`$ mainly goes via formation of the $`S_{11}`$-resonance, which means that $`t^{\gamma \eta }`$ in the near–threshold region can be written in a separable form similar to (5). To construct such separable T–matrix, we used the results of Ref.  where $`t^{\gamma \eta }`$ was obtained as an element of a multi–channel T–matrix which simultaneously describes data for the processes

$`\pi +N\pi +N,\pi +N\eta +N,`$
$`\gamma +N\pi +N,\gamma +N\eta +N`$
on the energy shell in the $`S_{11}`$–channel. For our calculations, we extended this T–matrix off the energy shell,

$$t_{off}^{\gamma \eta }(p^{},p;E)=\frac{\kappa ^2+E^2}{\kappa ^2+p_{}^{}{}_{}{}^{2}}t_{on}^{\gamma \eta }(E)\frac{\alpha ^2+2\mu E}{\alpha ^2+p^2},$$
(28)
using the Yamaguchi form–factors which become unit on the energy shell. Here $`\kappa `$ stands for some unknown parameter. It is known, that $`t^{\gamma \eta }`$ is different for neutron and proton, we assumed that they have the same functional form and differ by a constant factor $`t_n^{\gamma \eta }=At_p^{\gamma \eta }.`$ A multipole analysis  gives for this factor the following estimate: $`A=0.84\pm 0.15.`$

## 7 Photoproduction results

As was expected, our calculations revealed very strong final state interaction in the reaction (13). A comparison of the corresponding cross–sections obtained by solving the AGS equations and by using the Impulse Approximation (IA) is given in Fig. 6, where the IA-results are multiplied by 10. Besides the fact that the IA-curve is generally an order of magnitude lower, it does not show a resonant enhancement which is clearly seen when all the rescattering and re-arrangement processes are taken into account. In this connection, it should be noted that experimental data  for the reaction (13), given in Figs. 68 show a pronounced enhancement of the differential cross–section at low energies.

In order to examine a dependence of our calculations on the choice of the parameters of the T–matrices $`t^{\eta \eta }`$ and $`t^{\gamma \eta }`$, we did variations of $`A=t_n^{\gamma \eta }/t_p^{\gamma \eta }`$, Re $`a_{\eta N}`$ and $`\kappa `$ within the corresponding uncertainty intervals. One of the most important parameters of the theory is the ratio of the photoproduction amplitudes for neutron and proton ($`A`$). Six curves corresponding to different choices of $`A`$ are depicted in Fig. 6. These curves were calculated with $`a_{\eta N}=(0.75,0.30)`$ fm and $`\kappa =\alpha =3.316`$ fm-1. The experimental data are taken from Ref. .

In Fig. 8, we present the result of our calculations for five different choices of Re $`a_{\eta N}`$, namely, $`0.55`$ fm, $`0.65`$ fm, $`0.725`$ fm, $`0.75`$ fm, and $`0.85`$ fm. This sequence of Re $`a_{\eta N}`$ corresponds to the upward sequence of the curves in the near–threshold region. The parameters are: $`A=0.75`$ and $`\kappa =\alpha =3.316`$ fm-1.

In Fig. 8 the results of our calculations with three choices of the range parameter $`\kappa `$ (2 fm-1, 3 fm-1, and 5 fm-1) of the electromagnetic vertex $`\gamma NS_{11}`$ are given. These curves correspond to $`a_{\eta N}=(0.75,0.30)`$ fm and $`A=0.75`$. A comparison of the curves depicted in Figs. 68 and 8 with the corresponding experimental data shows that no agreement with the data can be reached unless the ratio $`A`$ is greater than $`0.80`$.

Under all variations of the parameters, however, the resonant peak remains about 15 MeV to the left of the experimental peak. Since this peak is due to the resonant final state interaction between the $`\eta `$ meson and deuteron, we may expect that it can be shifted to the right by introducing a repulsion into the $`\eta N`$ interaction. To introduce an $`\eta N`$ repulsion which preserves the separable form of the corresponding T–matrix, we used the method suggested in Ref.  where a separable nucleon–nucleon T–matrix includes an energy dependent factor, $`b(E)=\mathrm{tanh}(1E/E_c).`$ This factor causes the NN phase–shift to change sign at the energy $`E_c=0.816`$ fm-1, which is equivalent to presence of an NN repulsion. Since the purpose of our numerical experiment was to check if an $`\eta N`$ repulsion could shift the peak to the right and there is no information about such repulsion, we used the same function $`b(E)`$ and did variations of $`E_c`$, namely $`E_c/3`$, $`E_c/2`$, $`E_c`$, $`2E_c`$, and $`3E_c`$. The corresponding curves are shown in Fig. 9 where the larger $`E_c`$ the lower is the curve. All the curves depicted in Fig. 9 were calculated with $`a_{\eta N}=(0.75,0.30)`$ fm, $`A=0.85`$ and $`\kappa =\alpha =3.316`$ fm-1.

Therefore, comparison of our calculations with the experimental data suggests that $`A>0.80`$, Re $`a_{\eta N}>0.75`$ fm, and the $`\eta N`$ interaction is likely to be repulsive at short distances. Our calculations are highly sensitive to the ratio of the photoproduction amplitudes for neutron and proton ($`A`$) and depends rather weakly on the parameter $`\kappa `$.

Authors would like to thank Division for Scientific Affair of NATO for support (grant CRG LG 970110) and DFG-RFBR for financial assistance (grant 436 RUS 113/425/1).


text:
# On the Use of Optimized Monte Carlo Methods for Studying Spin Glasses

## Abstract

We start from recently published numerical data by Hatano and Gubernatis  to discuss properties of convergence to equilibrium of optimized Monte Carlo methods (bivariate multi canonical and parallel tempering). We show that these data are not thermalized, and they lead to an erroneous physical picture. We shed some light on why the bivariate multi canonical Monte Carlo method can fail.

One of the main problems of numerical results originated from large scale numerical simulations is that checking them is a task that is frequently of the order of magnitude of checking a real experiment: only repeating the full simulation, that demands availability of computer time and codes, allows a full check of the results.

Here we will use as a starting point the work of reference  to discuss a few points both about optimized Monte Carlo algorithms and about the behavior of $`3D`$ Edwards-Anderson (EA) spin glasses in the low $`T`$ phase. We will start by showing that the numerical results reported in reference , as far as the low $`T`$ values are concerned, are wrong: they are not equilibrium averages over the Boltzmann probability. Because of that the physical conclusions reached in the paper, supporting a trivial behavior of the broken phase of $`3D`$ spin glasses, are wrong. On the contrary recent numerical simulations  support, in this respect, a behavior of the system consistent with the Replica Symmetry Breaking (RSB) picture . We will also shed some light on why the optimized Monte Carlo method used in  can fail.

In the following we will first analyze our numerical data obtained by the Parallel Tempering Monte Carlo method , focusing on the analysis needed to establish that thermal equilibrium has been reached : we will use a large number of severe criteria that ensure that thermalization has been reached. After showing that the results of  are not correct in the low $`T`$ region we will discuss some preliminary simulations done using the same method used in , a bivariate version of the Multi-Canonical Monte Carlo , and we will point out a series of reasons for which a non careful implementation of this strategy can fail.

Let us start from our numerical data obtained through parallel tempering***For sake of a complete reliability and without fear of appearing over cautious we have chosen to rewrite all our codes in a double blind pattern, with two different sets of programmers, using different programming languages and different random number generators: they always give statistically compatible results.. We have simulated a $`3D`$ Edwards-Anderson spin glass, with binary random quenched couplings, linear size $`L=8`$ (the largest size used in ), down to $`T=0.50.5T_c`$: let us note that in our simulations for the same $`T`$ values we are able to thermalize reliably lattices up to $`L=16`$, and that we just discuss here results about the $`L=8`$ lattice, where we are completely confident about thermalization, only because this is the largest lattice studied in . We use a minimum value of the temperature $`T_{\text{min}}=0.5`$, a number of temperatures $`N_T=49`$ and a constant temperature step $`\delta T=\frac{1}{30}`$. The measured correlation times are always smooth functions of $`T`$ and no anomalies are detected.

Our data at high $`T`$ turn out to be statistically compatible with the ones of : in the high $`T`$ region there are no problems.

In figure 1 we plot the value of the Binder parameter,

$$B(t)\frac{1}{2}\left(3\frac{\overline{q^4(t)}}{\overline{q^2(t)}^2}\right),$$
(1)
averaged over logarithmic time windows, as a function of time at $`T=0.5`$ (close to $`0.5T_c`$). Averaging over logarithmic windows is the safe approach to check convergence in time. We first average over the last half of the total time extent of the run: this is the last point on the right of the plot. We subdivide in the same way the other half of the data, and the second point on the right is the average over the second half of this time span: we continue in this way till the origin of our Monte Carlo run. With a straight line we plot the asymptotic data from  as extracted from figure $`7`$ in the paper (since we were estimating by hand we have been generous on the statistical error): here there is no time dependence, we only plot with a straight line the asymptotic value. The discrepancy of our data and the data of  is very large and statistically very significant: definitely not an accident.

In figure 2 we plot the $`T=0.6`$ data from the same run, always for the Binder parameter averaged over logarithmic time windows: here $`T`$ is higher, and one could feel safer about thermalization, but again there is a clear and significant discrepancy among our data and the ones of . The dramatic stability of our data for $`B(t)`$ at low $`T`$ is already a very good indicator of a high level of thermalization. The results are stable at least during the last eight subdivisions of our two million step runs, i.e. at least from times going from $`10^4`$ to $`210^6`$.

In order to be sure we are not trapped in some metastable situation we have to check standard criteria about convergence, that in the case of optimized dynamics can be quite difficult to assert . Let us note for example that in recent numerical simulations  a careful discussion shows that weaker criteria can be sufficient to guarantee thermalization, making in this way possible to simulate more disorder sample with the same amount of computer time (since one needs less thermal sweeps per sample). Here, since thermalization is the main issue, we will check all of the most stringent criteria.

First of all we have checked the acceptance rates of the tempering sweeps in temperature: a bad choice of the $`T`$ values can make the swap of the temperature value too rare. In our case the rates are very high, of the order of $`.7`$ in all the temperature range: our parallel tempering scheme is performing very well.

Secondly we have checked, as customary, if all configurations (we have, as we said, $`49`$ of them) have spent a similar amount of time in each one of the $`49`$ allowed $`T`$ values. This criterion is important, since the first one could not be sufficient: spin configurations could be spending time swapping among neighboring $`T`$ values locally, but never leave the high or the low $`T`$ region. Our permanence histograms are very good: because of the large time extent of the runs all configurations have visited all regions of the $`T`$ phase space, and the permanence histograms are very flat. Again, this is a powerful test of thermalization.

The last point we have checked is the symmetry under the exchange $`qq`$ of the $`P_J(q)`$ for the individual samples. Since the overall flip of all spins is supposed to be a very slow mode of the dynamics, once we have good statistics on this mode we expect to have reached all the relevant regions of the phase space. Again, the symmetry is excellent for all individual samples (even for the more complex samples where the $`P_J(q)`$ has a non-trivial structure).

We consider this body of evidence as clear: our data are thermalized, the numerical data hint evidence in favor of the RSB picture (as confirmed by the data of , where even at very low $`T`$ values one sees that $`P(0)`$ does not depend on $`L`$) and the method used in  did not allow a proper thermalization.

In order to get a better understanding of the situation, and some hints about the reason of the failure of  we have implemented a code for rerunning their bivariate multi-canonical simulations.

Our simulations closely follows the description given in the Appendix of reference  and by Hatano himself . The analysis of few samples of sizes $`L=4`$, $`6`$, $`8`$ has been sufficient in order to understand where the thermalization problems may come from. Unless differently specified we have always used $`10^6`$ Monte Carlo Sweeps (MCS) for thermalizing and $`10^7`$ MCS for taking measurements in each multi-canonical cycle. The same number of MCS has been used by the authors of  only for $`L=10`$  (less iterations have been used for smaller lattice sizes).

The most delicate point during the thermalization process is the role played by the entropic barriers during the multi-canonical simulation. In a model which undergoes a first order transition the slowing down of the simulation at the critical point is essentially due to the presence of a huge energetic barrier between the two free energy minima. In this case the multi-canonical simulation works fine , and it rapidly converges towards a regime where every energy is sampled with the right probability, i.e. uniformly. Problems may arise when the multi-canonical method is applied to spin glasses or in general to models where entropic barriers play a central role. To this respect the study of its performances in models with only entropic barriers (e.g. backgammon model ) would be illuminating.

Let us focus now specifically on the $`3D`$ EA model, and see how the estimated density of states (DoS), $`D(e,q)`$, converges to the exact one. In particular we are interested in the histogram $`h(e,q)`$ which counts the number of times, during a multi-canonical cycle, the system is in a macroscopic state $`(e,q)`$ with energy $`e`$ and overlap $`q`$. Thermalization is achieved when $`h(e,q)`$ is flat and much larger than $`1`$ for all the physically allowed pairs $`(e,q)`$. Starting from a flat DoS, the region where $`h(e,q)1`$ broadens with the number of multi-canonical cycles and eventually reaches the boundaries of the allowed domain, $`e[e_0,e_0]q[1,1]`$, where $`e_0`$ is the ground-states energy (see the first two snapshots in figure 4, that we will discuss in better detail later on). In order to describe quantitatively the histogram evolution we plot in figure 3 the fraction of the $`(e,q)`$ space where $`h(e,q)0`$, that is the fraction of macroscopic $`(e,q)`$ configurations visited by the system during a multi-canonical cycle. We expect this fraction to increase more or less linearly during the first multi-canonical cycles and then to reach a plateau when simulation is thermalized (see figure 3.a, where things look good). For all the $`L=4`$ samples simulated we have observed this correct behavior. On the contrary for the $`L=6`$ samples, problems arise. At first, if the number of MCS is not large enough the simulation does not converge at all. In figure 3.b we show the results for the same sample shown in figure 3.a, with the only difference that $`10^6`$ MCS were used instead of $`10^7`$: here thermalization problems are evident, since in some situations the system simply gets trapped in a very small region of the phase space. In different samples we have found analogous problems also when using $`10^7`$ MCS (see figures 3.c and 3.d). With $`10^6`$ MCS the parallel tempering method is able to thermalize samples up to $`L=8`$ for temperatures down to $`T=0.3`$ (for example at the lowest $`T`$ value the Binder parameter thermalizes in $`10^6`$ MCS): the bivariate multi-canonical method does not seem to be very efficient for spin glasses.

In figure 4 we show the histogram evolution for sample #2 (the same used in figure 3.c). The four snapshots correspond to the black dots in figure 3.c and clearly show that the system, after reaching an apparently thermalized state with a flat and broad $`h(e,q)`$, instead of keeping it for all subsequent times, gets trapped in very small regions of the $`(e,q)`$ space (the third and fourth snapshots in figure 4).

How can we explain this behavior? During the first multi-canonical cycles the dynamics of the system in the $`(e,q)`$ space is diffusive in character, while when approaching the boundaries of the $`eq`$ plane (especially the energy ones) the system often gets trapped for very long times. The end of the diffusive behavior near to the ground states can be easily explained in terms of accessibility, that is the probability of decreasing the energy when the system is in a $`(e,q)`$ configuration and it makes a random move to a neighbor configuration. For not too low energies the accessibility is high: in this case a random walk in the configuration space corresponds to a random walk in the $`(e,q)`$ space, which is a projection of the previous one. On the contrary for energies close to the one of the ground states the accessibility is very low, due to the presence of a large number of higher local minima. For example if the system is at the bottom of a valley in the space of microscopic configurations, in order to further decrease its energy (a little step in the macroscopic $`(e,q)`$ space) it may need a long time, the time to find a deeper valley. The dynamics turns out to be strongly constrained for energies close to the boundaries.

Having in mind that the dynamics becomes slower and slower close to the energy boundaries, one can easily explain the peaks in figure 4. The system firstly relaxes in a uniform way on a large part of the $`(e,q)`$ space, the more accessible one. Still many allowed $`(e,q)`$ values are unvisited (because of the low accessibility), their DoS estimation becomes very small and their corresponding weights, $`W(e,q)=1/D(e,q)`$, huge. When the system reaches one of this configurations it can not leave it until the end of the multi-canonical cycle, when $`W(e,q)`$ will be updated again.

In order to improve the convergence we have also tried to start with a DoS estimated from the one of a thermalized $`L=4`$ sample. The convergence seems to be faster, however the problems giving rise to the peak structure in the histogram remain unaltered.

Given that the thermalization task appears to be very hard, one should at least try to use all thermalization checks available. For example the one based on the symmetry of the overlap distribution for every sample, $`P_J(q)`$ should always be carefully checked: this analysis is lacking in .

In figure 5 we show the overlap distribution $`P_J(q)`$ for the single $`L=6`$ sample considered in figure 3.d at a low temperature $`T=0.3`$ (these data come from a further parallel tempering simulation, pushed to lower $`T`$ values). In figure 5.a we show the $`P(q)`$ measured with a parallel tempering simulation. Its very accurate symmetry is a strong evidence of complete thermalization. In the next $`3`$ plots (b,c and d) we show with continuous lines the $`P(q)`$ measured with the multi-canonical method (the chosen times correspond to the dots in figure 3.d). We always superimpose the thermalized $`P(q)`$ for comparison. It is clear that, in the best case (see figure 5.b), the multi-canonical method is not able to give results as good as the parallel tempering does: in the worst cases it just gives a completely wrong $`P_J(q)`$, with a single or a double peak. The system may very easily get stuck somewhere, and in these cases the estimated $`P(q)`$ would look much narrower than the correct one (see figure 5.c and figure 5.d): measurements taken in such a biased situation hint for a fake evidence in favor of a single peak $`P(q)`$, and consequently of the droplet scenario.

As a last piece of evidence we consider the samples where the bivariate multi canonical has been well behaved: the scaling of the visited fraction of the $`(e,q)`$ phase space (for well thermalized samples) reported in figure 6 supports the picture of a diffusion-like evolution of the histogram. The area of support of the histogram grows more or less linearly with the number of multi-canonical cycles (the best exponent estimate is 0.9). Moreover, the time for reaching the plateau (equilibration time) grows with $`\tau L^{3.37}N^{1.12}`$, which seems to be very close to the theoretical lower bound ($`\tau N`$). However this result would hold only if the number of MCS per multi-canonical cycle necessary for a proper thermalization is independent from the system size $`N`$. As we have already seen this is not true. Indeed, using the same $`10^7`$ MCS per multi-canonical cycle, the fraction of well thermalized samples we have obtained is 100% for $`L=4`$, around 40% for $`L=6`$ and 0% for $`L=8`$. Because the requested number of MCS per multi-canonical cycle grows with $`N`$ (apparently very fast), our conclusion is that $`\tau `$ grows much faster than $`N`$ (simple arguments by Berg  suggest at least as $`N^2`$).

Concluding, we have seen how difficult it is to bring a bivariate multi-canonical simulation of spin glasses to equilibrium and, consequently, one possible reason of the failure of  to thermalize for $`L=8`$ (we have checked the failure of thermalization with independent parallel tempering simulations). When we say that the simulation is not thermalized we mean that we can not use the resulting DoSNote that the DoS estimation actually used in the measurements in  is $`D(e,q)h(e,q)`$ and so it is strongly affected by non-uniformities in the histogram. in order to estimate the observables averages at all the temperature. In particular, as long as the simulation does not visit many times the ground-states, we cannot believe to have enough information on the ground-states structure. However it may perfectly be that, after a certain number of multi-canonical cycles, the estimated DoS gives good averages at higher temperatures, which do not change if new low energy states are reached. We believe this is the case in , where data at not too low temperatures are perfectly compatible with the ones obtained in previous work and fit the RSB scenario.

We thank N. Hatano for an useful correspondence regarding the bivariate method.


text:
# Isotope effects in underdoped cuprate superconductors: a quantum critical phenomenon

## Abstract

We show that the unusual doping dependence of the isotope effects on transition temperature and zero temperature in - plane penetration depth naturally follows from the doping driven 3D-2D crossover, the 2D quantum superconductor to insulator transition (QSI) in the underdoped limit and the change of the relative doping concentration upon isotope substitution. Close to the QSI transition both, the isotope coefficient of transition temperature and penetration depth approach the coefficient of the relative dopant concentration, and its divergence sets the scale. These predictions are fully consistent with the experimental data and imply that close to the underdoped limit the unusual isotope effect on transition temperature and penetration depth uncovers critical phenomena associated with the quantum superconductor to insulator transition in two dimensions.

PACS numbers: 74.20.Mn, 74.25Dw, 5.40.-a

The observation of an unusual isotope effect in underdoped cuprate superconductors on transition temperature  and zero temperature penetration depth  poses a challenge to the understanding of high temperature superconductivity. In this letter we show that this behavior naturally follows from the doping driven 3D-2D crossover, the 2D superconductor to insulator (QSI) transition in the underdoped limit and the change of the relative doping concentration upon isotope substitution. Thus, the unusual isotope effect on superconducting properties of underdoped cuprates is derived to be a quantum critical phenomenon driven by variation of the dopant concentration.

Consider the empirical phase diagram of La2-xSrxCuO4 depicted in Fig.1. It shows that after passing the so called underdoped limit $`\left(x=x_u0.05\right)`$, where the quantum insulator to superconductor (QSI) transition occurs , $`T_c`$ rises and reaches a maximum value $`T_c^m`$ at $`x_m0.15`$. With further increase of $`x`$, $`T_c`$ decreases and finally vanishes in the overdoped limit $`x_o0.3`$. This phase transition line $`T_c\left(x\right)`$, separating the superconducting from the normal conducting phase, appears to be a generic property of cuprate superconductors. In La2-xSrxCuO4, HgBa2CuO4+x  and Bi2Sr2CuO6+x both, the underdoped and overdoped limits, corresponding to critical endpoints, are experimentally accessible. In other cuprates, including Bi2Sr2CaCu2O8+x , YBa2Cu3O7-x  and Y$`{}_{1x}{}^{}\mathrm{Pr}_{x}^{}`$Ba2Cu3O7-δ , only the underdoped and optimally doped regimes appear to be accessible. As shown in Fig.1 for La2-xSrxCuO4 the effective mass anisotropy, measured in terms of $`\gamma =\sqrt{M_c/M_{ab}}`$ increases drastically by approaching the underdoped limit. This property, also observed in HgBa2CuO4+δ  and YBa2Cu3O7-x , appears to be generic and reveals the crossover from three (3D) to two dimensional (2D) behavior.

In our considerations the starting point is the critical endpoint of the phase transition line $`T_c\left(x\right)`$, where at $`T=0`$ and $`x=x_u`$, driven by variation of the dopant concentration $`x`$, the QSI transition occurs. Close to this critical endpoint the bulk superconductors correspond to a stack of independent superconducting slabs of thickness $`d_s`$. Moreover, close to such a critical point, low energy properties depend only on the spatial dimensionality of the system, the number of components of the order parameter and the range of the interaction. The theory of quantum critical phenomena predicts that the transition temperature and zero temperature in - plane penetration depth scale as

$$T_c=a\delta ^{z\overline{\nu }},\frac{1}{\lambda _{ab}^2\left(0\right)}=b\delta ^{z\overline{\nu }}.$$
(1)
$`z`$ is the dynamic critical exponent, $`\overline{\nu }`$ the exponent of the diverging length, $`a`$ and $`b`$ nonuniversal critical amplitudes, and

$$\delta =xx_u$$
(2)
measures the distance from the quantum critical point at $`x_u`$, where in cuprate superconductors the QSI transition occurs. For a complex scalar order parameter and in D=2 the critical amplitudes $`a`$ and $`b`$ and the slab thickness $`d_s`$ are not independent but related by the universal relation

$$\underset{\delta 0}{lim}\frac{T{}_{c}{}^{}\lambda _{ab}^{2}(0)}{d_s}=\frac{a}{bd_s}=\frac{1}{\overline{Q}_2}\left(\frac{\mathrm{\Phi }_0^2}{16\pi ^3k_B}\right),$$
(3)
where $`\overline{Q}_2`$ is a universal number. Although the experimental data are rather sparse close to the QSI transition, the overall picture turns out to be highly suggestive and provides consistent evidence for the QSI transition in D=2 at the critical endpoint in the underdoped limit with critical exponents $`z1`$ and $`\overline{\nu }1`$. This estimate is close to theoretical predictions , from which $`z=1`$ is expected for a bosonic system with long-range Coulomb interactions independent of dimensionality and $`\overline{\nu }11.03`$ in $`D=2`$. In Fig.1 it is seen that $`\overline{\nu }1`$ is also consistent with the doping dependence of the effective mass anisotropy $`\gamma =\sqrt{M_c/M_{ab}}\delta ^{\overline{\nu }}`$.

We are now prepared to explore the implications of the QSI transition on the isotope effect. From the definition of the isotope coefficient

$$\beta _{T_c}=\frac{m}{T_c}\frac{dT_c}{dm}$$
(4)
and Eq. (1) we obtain the scaling expression

$$\beta _{T_c}=\beta _a+\beta _\delta ,$$
(5)
where

$$\beta _a=\frac{m}{a}\frac{da}{dm},\beta _\delta =\frac{z\overline{\nu }}{\delta }\overline{\beta }_\delta ,\overline{\beta }_\delta =m\frac{d\delta }{dm}.$$
(6)
Since in the doping regime of interest, isotope substitution lowers the transition temperature , while the dopant concentration $`x`$ remains nearly unchanged , there is a positive shift of the underdoped limit $`x_u`$, and $`\overline{\beta }_\delta `$ reduces to

$$\overline{\beta }_\delta =m\frac{d\delta }{dm}m\frac{dx_u}{dm}>0.$$
(7)
Thus, by approaching the QSI transition $`\beta _\delta `$ diverges as $`\beta _\delta \delta ^1`$ \[Eq.(4)\] and provided that $`\beta _a`$ remains finite, $`\beta _{T_c}`$ is predicted to tend to $`\beta _\delta `$ so that

$$\frac{1}{\beta _{T_c}}\frac{1}{\beta _\delta }=\overline{r}\left(\frac{T_c}{T_c^m}\right)^{1/z\overline{\nu }},\overline{r}=\frac{T_c^m}{z\overline{\nu }\overline{\beta }_\delta a^{1/z\overline{\nu }}}.$$
(8)
Here we expressed $`\delta `$ in terms of $`T_c`$ \[Eq.(1)\] and rescaled $`T_c`$ by $`T_c^m`$, the transition temperature at optimum doping, to reduce variations of $`T_c`$ between different cuprates .

To verify this prediction we show $`1/\beta _{_{T_c}}`$ versus $`T_c/T_c^m`$ for La1.85Sr0.15Cu1-xO4 , YBa2-xLaxCu3O7  and Y$`{}_{1x}{}^{}\mathrm{Pr}_{x}^{}`$Ba2Cu3O7  in Fig.2. As predicted, by approaching the underdoped limit, corresponding to $`T_c/T_c^m=0`$, the data collapses on a single curve, consistent with a straight line pointing to the expected value $`z\overline{\nu }1`$ , yielding the estimate $`\overline{r}6`$ and confirming that $`\beta _a`$ is finite and $`\overline{\beta }_\delta >0`$. Thus, the OSI transition in $`D=2`$ accounts naturally for the unusual doping dependence of $`\beta _{T_c}`$ and identifies the isotope induced change of the underdoped limit as the main mechanism. Notice that the empirical estimate for $`\overline{r}`$ implies that the maximum transition temperature is fixed by $`T_c^m6a\overline{\beta }_\delta `$.

Another important element brought by the QSI transition in D=2 is the universal relation given in Eq. (3). It implies with Eq.(5) that the isotope coefficients of $`T_c`$, $`1/\lambda _{ab}^2`$ , critical amplitudes $`a`$ and $`b`$ and slab thickness $`d_s`$ are related by

$$\beta _{T_c}=\beta _{1/\lambda _{ab}^2}+\beta _{d_s}=\beta _a+\beta _\delta ,\beta _a=\beta _b+\beta _{d_s},$$
(9)
where $`\beta _F=\frac{m}{F}\frac{dF}{dm}`$ and $`F=T_c,`$ $`1/\lambda _{ab}^2,d_s,a`$ and $`b`$. Noticing that $`\beta _a`$ was confirmed to be bounded, this is also true for $`\beta _b`$ and $`\beta _{d_s}`$. Since $`\beta _{T_c}`$ diverges as $`\beta _{T_c}=\beta _\delta =\left(T_c/T_c^m\right)^{1/z\overline{\nu }}/\overline{r}`$ \[Eq.(8\], $`\beta _{1/\lambda _{ab}^2}`$ is predicted to approach

$$\beta _{1/\lambda _{ab}^2}=\beta _{T_c}=\beta _\delta =\frac{1}{\overline{r}}\left(\frac{T_c}{T_c^m}\right)^{1/z\overline{\nu }}$$
(10)
close to the QSI transition. Although the experimental data for $`\beta _{1/\lambda _{ab}^2}`$ and $`\beta _{T_c}`$ on identical samples are rather sparse, the results shown in Fig.3 for the oxygen isotope effect in La2-xSrxCu1-xO4  clearly confirm this prediction. Indeed as the underdoped limit is approached the data points tend to the solid line, marking $`1/\beta _{T_c}=1/\beta _{1/\lambda _{ab}^2}`$.

In conclusion, we have shown that the strong doping dependence of the isotope effects on transition temperature and zero temperature in - plane penetration depth naturally follows from the doping driven 3D-2D crossover and the 2D quantum superconductor to insulator transition in the underdoped limit. As the quantum superconductor to insulator transition is approached, the isotope coefficient of transition temperature $`\beta _{T_c}`$ and penetration depth $`\beta _{1/\lambda _{ab}^2}`$ tend to the coefficient of the relative dopant concentration $`\beta _\delta `$, and its divergence sets the scale. These predictions are fully consistent with the experimental data and imply that close to the underdoped limit the unusual isotope effects on transition temperature and penetration depth uncovers critical phenomena, associated with the quantum superconductor to insulator transition in two dimensions. Moreover, since lattice distortions are the primary consequence of isotope substitution it becomes clear that the isotope effects on transition temperature and zero temperature penetration depth probe the effect of lattice degrees of freedom on the condensate. However, until now most of the proposed theories of high temperature superconductivity in the cuprates have considered charge and spin fluctuations in a doped two dimensional (2D) antiferromagnet and disregarded any coupling to the lattice degrees of freedom.

The authors are grateful to K.A. Müller, S. Roos and G.M. Zhao for very useful comments and suggestions on the subject matter. This work was partially supported by the Swiss National Science Foundation.


text:
# References

A coherent approach to Spacetime Foam

R. Garattinia

a Facoltà di Ingegneria, Università di Bergamo, Viale Marconi 5, 24044 Dalmine (Bergamo), Italy.

E-mail: Garattini@mi.infn.it

A coherent superposition of $`N`$ Schwarzschild wormholes is proposed as a model for spacetime foam. Following the subtraction procedure for manifolds with boundaries, we calculate by variational methods the Casimir energy. A proposal for an alternative foamy model formed by $`N`$ Schwarzschild-Anti-de Sitter wormholes is here considered. Finally, a conjecture about the foam evolution is proposed.

Spacetime Foam; Black Holes; Quantum Gravity.

The term Spacetime Foam was used for the first time by J.A. Wheeler to indicate that spacetime may be subjected to quantum fluctuations in topology and metric at the Planck scale . In a series of papers we have proposed a model for such a foamy space made by $`N`$ coherent Schwarzschild wormholes . The relevance of this proposal is based on a Casimir energy computation showing that quantum fluctuations of the gravitational metric shift the minimum of the effective energy from flat space ( the classical minimum of the energy) to a multi-wormhole configuration. Indeed the final energy contribution to one loop is

$$\mathrm{\Delta }_{N_w}E\left(M\right)N_w^2\frac{V}{64\pi ^2}\frac{\mathrm{\Lambda }^4}{e},$$
(1)
where $`V`$ is the volume of the system, $`\mathrm{\Lambda }`$ is the U.V. cut-off and $`N_w`$ is the wormholes number . This expression shows that a non-trivial vacuum of the multi-wormhole type is favoured with respect to flat space. It is important to remark that it is the $`N`$ \- coherent superposition of wormholes that it is privileged with respect to flat space and not the single wormhole, because the single wormhole energy contribution has an imaginary contribution in its spectrum: a clear sign of an instability. Nevertheless the presence of an unstable mode is necessary to have transition from one vacuum (the false one) to the other one (the true vacuum)  . Three consequences of this multiply connected spacetime are:

1. the event horizon area of a black hole is quantized and by means of the Bekenstein-Hawking relation , also the entropy of a black hole is quantized. In particular for a Schwarzschild black hole

$$M=\frac{\sqrt{N}}{2l_p}\sqrt{\frac{\mathrm{ln}2}{\pi }},$$
(2)
namely the black hole mass is quantized. Here $`l_p`$ is the Planck length in natural units.
2. A cosmological constant is induced by vacuum fluctuations as shown by $`\left(\text{1}\right)`$ whose value is

$$\mathrm{\Lambda }_c=\frac{\mathrm{\Lambda }^4l_p^2}{N_w8e\pi }.$$
(3)
When the area-entropy relation is applied to the de Sitter geometry, we obtain

$$\frac{3\pi }{\mathrm{ln}2l_p^2N_w}=\mathrm{\Lambda }_c.$$
(4)
3. Combining $`\left(\text{3}\right)`$ and $`\left(\text{4}\right)`$, one gets

$$\mathrm{\Lambda }_c=\frac{\mathrm{\Lambda }^4l_p^2}{N_w8e\pi }=\frac{3\pi }{\mathrm{ln}2l_p^2N_w}.$$
(5)
This means that we have found a constraint on the U.V. cut-off

$$\mathrm{\Lambda }^4=\frac{24e\pi ^2}{\mathrm{ln}2l_p^4}.$$
(6)

It is interesting to note that $`\left(\text{1}\right)`$ can be obtained at least for $`N_w=1`$ even for Schwarzschild-Anti-de Sitter wormholes (S-AdS) , whose line element is $`ds^2=`$

$$f\left(r\right)dt^2+f\left(r\right)^1dr^2+r^2d\mathrm{\Omega }^2,$$
(7)
where

$$f\left(r\right)=\left(1\frac{2MG}{r}+\frac{r^2}{b^2}\right).$$
(8)
$`\mathrm{\Lambda }_{AdS}`$ is the negative cosmological constant and $`b=\sqrt{\frac{3}{\mathrm{\Lambda }_{AdS}}}`$. To consider a large $`N_w`$ approach to spacetime foam even with S-AdS wormholes, one has to consider the following rescaling

$$\{\begin{array}{c}R_\pm R_\pm /N_w\\ l_p^2N_wl_p^2\\ \mathrm{\Lambda }_{AdS}\mathrm{\Lambda }_{AdS}/N_w^2\end{array},$$
(9)
where $`R_\pm `$ are the boundaries related to the single wormhole. This rescaling is a consequence of the boundary reduction related to the coherent superposition of wormholes wave functionals leading to the stabilization of the system This means that a selection rule has to emerge to compare the quantity

$$\mathrm{\Gamma }_{\mathrm{N}\mathrm{S}\mathrm{AdS}\mathrm{holes}}$$

$$=\frac{P_{\mathrm{N}\mathrm{S}\mathrm{AdS}\mathrm{holes}}}{P_{\mathrm{AdS}}}\frac{P_{\mathrm{foam}}}{P_{\mathrm{AdS}}}$$
(10)
with

$$\mathrm{\Gamma }_{\mathrm{N}\mathrm{S}\mathrm{holes}}$$

$$=\frac{P_{\mathrm{N}\mathrm{S}\mathrm{holes}}}{P_{\mathrm{flat}}}\frac{P_{\mathrm{foam}}}{P_{\mathrm{flat}}}.$$
(11)
In both cases, we find a non-vanishing probability that a non-trivial vacuum has to be considered. Moreover the rescaling in $`\left(\text{9}\right)`$ leaves the potential $`\left(\text{8}\right)`$ invariant and when $`N_w`$ is very large $`\mathrm{\Lambda }_{AdS}0`$. This seems to suggest that not only $`\left(\text{10}\right)`$ and $`\left(\text{11}\right)`$ have to be compared, but it is likely that a hierarchical mechanism of the type

$$\begin{array}{c}\frac{P_{\mathrm{N}\mathrm{S}\mathrm{AdS}\mathrm{holes}}}{P_{\mathrm{AdS}}}\frac{P_{\mathrm{foam}}}{P_{\mathrm{AdS}}}\\ _{\mathrm{\Lambda }_{AdS}0}\\ FlatSpace\\ \\ \frac{P_{\mathrm{N}\mathrm{S}\mathrm{holes}}}{P_{\mathrm{flat}}}\frac{P_{\mathrm{foam}}}{P_{\mathrm{flat}}}\end{array},$$
(12)
when $`N_w\mathrm{}`$. Although this picture has to be examined in detail and $`\left(\text{12}\right)`$ is only at the conjecture level, an important question comes into play: why a spacetime formed by wormholes has to be preferred with respect to flat spacetime, when the last one is the one we observe. The answer that at this stage can only be conjectured is that if we consider the following expectation value on the foam state

$$\frac{\mathrm{\Psi }_F\left|\widehat{g}_{ij}\right|\mathrm{\Psi }_F}{\mathrm{\Psi }_F|\mathrm{\Psi }_F},$$
(13)
when the number of wormholes is large enough, i.e. the scale is sufficiently large, we should have to obtain

$$\frac{\mathrm{\Psi }_F\left|\widehat{g}_{ij}\right|\mathrm{\Psi }_F}{\mathrm{\Psi }_F|\mathrm{\Psi }_F}\eta _{ij},$$
(14)
where $`\eta _{ij}`$ is the flat space metric. This is a test that this foamy model has to pass if phenomenological aspects have to be considered.

Acknowledgments

I would like to thank Prof. R. Bonifacio who has given to me the opportunity of participating to the Conference.