text:
# 1 The induced magnetization 𝑚(ϵ)=𝑀_𝑧/𝑀₀ as a function of the normalized applied electric field ϵ=𝐸_𝑥/𝐸₀, according to Eq.(10). The inset shows a blow up of the paramagnetic region of 𝑚(ϵ) (exhibiting its linear superconducting magnetoelectric effect regime).

NOVEL EFFECTS IN HIGH-TC GRANULAR SUPERCONDUCTORS

PREDICTED BY A MODEL OF 3D JOSEPHSON JUNCTION ARRAYS

SERGEI SERGEENKOV

Bogoliubov Laboratory of Theoretical Physics,

Joint Institute for Nuclear Research, Dubna, Russia

1. INTRODUCTION

Despite the fact that Josephson Junction Arrays (JAA) have been actively studied for decades, they continue to contribute to the variety of intriguing and peculiar phenomena (both fundamental and important for potential applications) providing at the same time a useful tool for testing new theoretical ideas. To give just a few recent examples, it is sufficient to mention paramagnetic Meissner effect (PME)\[1-5\] as well as the recently introduced thermophase  and piezophase  effects suggesting, respectively, a direct influence of a thermal gradient and an applied stress on phase difference between the adjacent grains. At the same time, an artificially prepared islands of superconducting grains , well-described by the various models of JJAs, proved useful in studying the charging effects in these systems, ranging from Coulomb blockade of Cooper pair tunneling and Bloch oscillations  to propagation of quantum ballistic vortices .

The present paper reviews some of the recently suggested novel effects \[14-16\] which should manifest themselves either in weak-links-bearing superconductors or in an artificially prepared JJAs. In Section 2 we consider the appearance of an electric-field induced magnetization (analog of the so-called magnetoelectric effect) in a model granular superconductor. The dual effect, that is an appearance of a magnetic field induced electric polarization is studied in Section 3. Finally, in Section 4 we discuss a possibility of two other interesting effects which are expected to occur in a granular material under mechanical loading. Specifically, we predict the existence of stress induced paramagnetic moment in zero applied magnetic field (Josephson piezomagnetism) and study its influence on a low-field magnetization (leading to a mechanically induced PME).

2. MAGNETOELECTRIC EFFECT

To account for the unusual behavior of the critical current under the influence of a strong enough applied electric field observed in high-$`T_c`$ superconductos (HTCS) \[17-20\], a possibility of the superconducting analog of the so-called magnetoelectric effect (MEE) in JJAs has been recently suggested . This effect is similar to (but physically different from) the MEE seen in magnetic ferroelectrics (like antiferromagnetic $`BiFeO_3`$) . It has also been discussed in the context of some exotic nonmagnetic normal metal conductors, with the symmetry of a mirror isomer , and pyroelectric superconductors (where the supercurrent passing through a metal of a polar symmetry is assumed to be accompanied by spin polarization of the carriers ). The effect discussed here entails the electric field generation of a magnetic moment (paramagnetic in small to diamagnetic for large electric fields) due to superconducting currents that circulate between grains. The basic physical reason for the appearance of this effect is that the applied electric field induces a magnetization $`M`$, that changes its sign ($`MM`$) under a time-parity transformation. In the case of weakly-coupled superconducting grains, the phenomenological reason for the time-parity violation comes from the fact that the total free energy $`F`$ for the material exhibiting both magnetic and electric properties contains the term $`\alpha _{mn}f_{mn}(\stackrel{}{E},\stackrel{}{H})`$, with the coefficients $`\alpha _{mn}0`$ (here $`\{m,n\}=x,y,z`$). In the standard MEE, the function $`f_{mn}E_mH_n`$, which leads to the corresponding linear effect, which is either an electric-field induced magnetization $`M_m(\stackrel{}{E})F/H_m=\alpha _{mn}E_n`$ (in zero magnetic field) or a magnetic-field induced polarization $`P_m(\stackrel{}{H})F/E_m=\alpha _{mn}H_n`$ (in zero electric field). Since $`H`$ (and $`M`$) changes its sign under time-parity transformation while $`E`$ (and $`P`$) remains unchanged, an electric field induced magnetization will break time-parity symmetry even in zero magnetic field applied. As we show below, in our case the symmetry breaking term in $`F`$, represented by a nonzero coefficient $`\alpha _{mn}`$, has a more general nonlinear form for $`f_{mn}`$. We recall that the standard linear MEE can appear when an external electric field $`\stackrel{}{E}`$ interacts with an inner magnetic field $`\stackrel{}{h}_{DM}`$ of the Dzyaloshinskii-Moriya (DM) type . The DM interaction leads to a term, apart from the standard isotropic term in the Heisenberg Hamiltonian, of the form $`H_{DM}=_{i,j}\stackrel{}{D}_{i,j}(\stackrel{}{S}_i\stackrel{}{S}_j)`$, where $`\stackrel{}{S}_i`$ is a Heisenberg spin and the constant vector $`\stackrel{}{D}_{i,j}`$ arises from the spin-orbit coupling. An analogous situation occurs in our case, as we describe below. To see how we can get a nonzero $`\alpha `$, or equivalently a DM type interaction in a granular superconductor, we model a HTCS ceramic sample by a random three-dimensional (3D) overdamped Josephson junction array. This model has proven to be useful in describing the metastable magnetic properties of HTCS . In thermodynamic equilibrium, this model has a Boltzmann factor with a random 3D-XY model Hamiltonian. Specifically, the general form of the Hamiltonian (describing both DC and AC effects) reads

$$(t)=\underset{i,j}{}J(r_{i,j})[1\mathrm{cos}\varphi _{i,j}(t)]$$
(1)
Here $`\{i\}=\stackrel{}{r}_i`$ is a 3D lattice vector; $`J(r_{i,j})`$ is the Josephson coupling energy, with $`\stackrel{}{r}_{i,j}=\stackrel{}{r}_i\stackrel{}{r}_j`$ the separation between the grains; the gauge invariant phase difference is defined as

$$\varphi _{i,j}(t)=\varphi _{i,j}(0)A_{i,j}(t),$$
(2)
where $`\varphi _{i,j}(0)=\varphi _i\varphi _j`$ with $`\varphi _i`$ being the phase of the superconducting order parameter; $`A_{i,j}(t)`$ is (time-dependent, in general) frustration parameter, defined as

$$A_{i,j}(t)=\frac{2\pi }{\mathrm{\Phi }_0}_i^j\stackrel{}{A}(\stackrel{}{r},t)𝑑\stackrel{}{l},$$
(3)
with $`\stackrel{}{A}(\stackrel{}{r},t)`$ the (space-time dependent) electromagnetic vector potential which involves both external fields and the electric and magnetic possible self-field effects (see below); $`\mathrm{\Phi }_0=h/2e`$ is the quantum of flux, with $`h`$ Planck’s constant, and $`e`$ the electronic charge. Expanding the cosine term, and using trigonometric identities we can explicitly rewrite the above Hamiltonian as $`=_{i,j}J[\mathrm{cos}(A_{i,j})\stackrel{}{S}_i\stackrel{}{S}_j\mathrm{sin}(A_{i,j})\widehat{k}\stackrel{}{S}_i\stackrel{}{S}_j]`$, where the two-component XY spin vector is defined as $`\stackrel{}{S}_i(\mathrm{cos}\varphi _i,\mathrm{sin}\varphi _i)`$, and $`\widehat{k}`$ is a unit vector along the z-axis . We see that the second term in this Hamiltonian has the same form as in the DM contribution, and thus we can surmise that the time parity will be broken by applying an external field (electric or magnetic) to the granular system. To bring to the fore this possibility, we show below in a simple but yet nontrivial model that this is indeed the case.

There are different types of $`A_{i,j}`$ randomness that can be considered . For simplicity, in the present paper, we consider a long-range interaction between grains (assuming $`J(r_{i,j})=J`$) and model the true short-range behavior of a ceramics sample through the randomness in the position of the superconducting grains in the array (using the exponential distribution law $`f_r(r_{i,j})`$, see below). Here we restrict our consideration to the case of an external electric field only but it can be shown that the scenario suggested will also carry through when applying an external magnetic field (which will induce another time-parity breaking phenomenon in the granular material; namely, magnetic field induced electric polarizability, see the next Section). Besides, in what follows we also ignore the role of Coulomb interaction effects assuming that the grain’s charging energy $`E_CJ`$ (where $`E_C=e^2/2C`$, with C the capacitance of the junction).

In the case of a granular material, we show here that the corresponding time-parity breaking DM internal field can be related to the electric field induced magnetic moment produced by the circulating Josephson currents between the grains. As is known , a constant electric field $`\stackrel{}{E}`$ applied to a single Josephson junction (JJ) causes a time evolution of the phase difference. In this particular case Eq.(2) reads

$$\varphi _{i,j}(t)=\varphi _{i,j}(0)+\frac{2e}{\mathrm{}}\stackrel{}{E}\stackrel{}{r}_{i,j}t$$
(4)
The resulting AC superconducting current in the junction is

$$I_{i,j}^s(t)=\frac{2eJ}{h}\mathrm{sin}\varphi _{i,j}(t)$$
(5)
If, in addition to the external electric field $`\stackrel{}{E}`$, the network of superconducting grains is under the influence of an applied magnetic field $`\stackrel{}{H}`$, the frustration parameter $`A_{i,j}(t)`$ in Eq.(3) takes the following form

$$A_{i,j}(t)=\frac{\pi }{\mathrm{\Phi }_0}(\stackrel{}{H}\stackrel{}{R}_{i,j})\stackrel{}{r}_{i,j}\frac{2\pi }{\mathrm{\Phi }_0}\stackrel{}{E}\stackrel{}{r}_{i,j}t$$
(6)
Here, $`\stackrel{}{R}_{i,j}=(\stackrel{}{r}_i+\stackrel{}{r}_j)/2`$, and we have used the conventional relationship between the vector potential $`\stackrel{}{A}`$ and a constant magnetic field $`\stackrel{}{H}=rot\stackrel{}{A}`$ (with $`\stackrel{}{H}/t=0`$), as well as a homogeneous electric field $`\stackrel{}{E}=\stackrel{}{A}/t`$ (with $`rot\stackrel{}{E}=0`$). In the type II HTCS the magnetic self-field effects for the array as a whole are expected to be negligible . The grains themselves are in fact larger than the London penetration depth and we must then have that the corresponding Josephson penetration length must be much larger than the grain size (since the self-induced magnetic fields can in principle be quite pronounced for large-size junctions even in zero applied magnetic fields ). Specifically, this is justified for short junctions with the size $`d\lambda _J`$, where $`\lambda _J=\sqrt{\mathrm{\Phi }_0/4\pi \mu _0j_c\lambda _L}`$ is the Josephson penetration length with $`\lambda _L`$ being the grain London penetration depth and $`j_c`$ its Josephson critical current density. In particular, since in HTCS $`\lambda _L150nm`$, the above condition will be fulfilled for $`d1\mu m`$ and $`j_c10^4A/m^2`$ which are the typical parameters for HTCS ceramics . Likewise, to ensure the uniformity of the applied electric field, we also assume that $`d\lambda _E`$, where $`\lambda _E`$ is an effective electric field penetration depth .

When the AC supercurrent $`I_{i,j}^s(t)`$ (defined by Eqs.(2), (5) and (6)) circulates around a set of grains, that form a random area plaquette, it induces a random AC magnetic moment $`\stackrel{}{\mu }_s(t)`$ of the Josephson network

$$\stackrel{}{\mu }_s(t)\left[\frac{}{\stackrel{}{H}}\right]_{\stackrel{}{H}=0}=\underset{i,j}{}I_{i,j}^s(t)(\stackrel{}{r}_{i,j}\stackrel{}{R}_{i,j})$$
(7)
Notice that in the MEE-like effect discussed here for a granular superconductor, the electric-field induced magnetic moment in the system is still present in zero applied magnetic field due to the phase coherent currents between the weakly-coupled superconducting grains.

To consider the essence of the superconducting analog of MEE, we assume that in a zero electric field the phase difference between the adjacent grains $`\varphi _{i,j}(0)=0`$ which corresponds to a fully coherent state of the array. In this particular case, the electric-field induced averaged magnetization reads

$$\stackrel{}{M}_s(\stackrel{}{E})\overline{\stackrel{}{\mu }_s(t)}=\frac{1}{\tau }\underset{0}{\overset{\tau }{}}𝑑t\underset{0}{\overset{\mathrm{}}{}}𝑑\stackrel{}{r}_{i,j}𝑑\stackrel{}{R}_{i,j}f(\stackrel{}{r}_{i,j},\stackrel{}{R}_{i,j})\stackrel{}{\mu }_s(t),$$
(8)
where $`\tau `$ is the electronic relaxation scattering time, and $`f`$ is the joint probability distribution function (see below).

To obtain an explicit expression for the electric-field dependent magnetization, we consider a site positional disorder that allows for small random radial displacements. Namely, the sites in a 3D cubic lattice are assumed to move from their equilibrium positions according to the normalized (separable) distribution function

$$f(\stackrel{}{r}_{i,j}\stackrel{}{R}_{i,j})f_r(\stackrel{}{r}_{i,j})f_R(\stackrel{}{R}_{i,j})$$
(9)
It can be shown that the main qualitative results presented here do not depend on the particular choice of the probability distribution function. For simplicity here we assume an exponential distribution law for the distance between grains, $`f_r(\stackrel{}{r})=f(x_1)f(x_2)f(x_3)`$ with $`f_r(x_j)=(1/d)e^{x_j/d}`$, and a short range distribution for the dependence of the center-of-mass probability $`f_R(\stackrel{}{R})`$ (around some constant value $`D`$). The specific form of the latter distribution will not affect the qualitative nature of the final result. (Notice that in fact the former distribution function $`f_r(\stackrel{}{r})`$ reflects a short-range character of the Josephson coupling in granular superconductor. Indeed, according to the conventional picture  the Josephson coupling $`J(\stackrel{}{r}_{ij})`$ can be assumed to vary exponentially with the distance $`\stackrel{}{r}_{ij}`$ between neighboring grains, i.e., $`J(\stackrel{}{r}_{ij})=Je^{\stackrel{}{\kappa }\stackrel{}{r}_{ij}}`$. For isotropic arrangement of identical grains, with spacing $`d`$ between the centers of adjacent grains, we have $`\stackrel{}{\kappa }=(\frac{1}{d},\frac{1}{d},\frac{1}{d})`$ and thus $`d`$ is of the order of an average grain size.) Taking the applied electric field along the $`x`$-axis, $`\stackrel{}{E}=(E_x,0,0)`$, we get finally

$$M_z(E_x)=\frac{B_z(E_x)}{\mu _0}H_z(E_x),$$
(10)
for the induced transverse magnetization (along the $`x_3=z`$-axis), where

$$B_z(E_x)=\mu _0M_0\frac{E_x/E_0}{1+(E_x/E_0)^2},$$
(11)
and

$$H_z(E_x)=M_0\left(\frac{E_0}{E_x}\right)\mathrm{log}\sqrt{1+\left(\frac{E_x}{E_0}\right)^2},$$
(12)
stand for the electric-field induced magnetic induction $`B_z(E_x)`$ and magnetic field $`H_z(E_x)`$, respectively. The induced Josephson current $`I(E_x)`$ is simply given by Ampere’s law $`I(E_x)=H_z(E_x)d`$. In these equations, $`M_0=2\pi eJNdD/\mathrm{}`$, with $`N`$ the total number of grains and $`E_0=\mathrm{}/2de\tau `$. Eq.(10) is the main result of this Section, which we proceed to analyze below. As is seen from Eq.(10), the behavior of the magnetization in the applied electric field is determined by the competition between the two contributions, the magnetic induction $`B_z(E_x)`$ and the current induced magnetic field $`H_z(E_x)`$ (or the corresponding Josephson current $`I(E_x)`$). Namely, below a critical (threshold) field $`E_c1.94E_0`$ where $`B_z(E_x)>\mu _0H_z(E_x)`$, a paramagnetic phase of the MEE is due to the modification of the magnetic induction in the applied electric field. On the other hand above the threshold field (when $`E_x`$ becomes larger than $`E_c`$) the Josephson current $`I(E_x)`$ induced contribution starts to prevail, leading to the appearance of the diamagnetic signal (seen as a small negative part of the induced magnetization in Fig.1).

Such an electric field induced paramagnetic-to-diamagnetic transition has been actually observed \[17-19\] for behavior of the critical current in ceramic HTCS and it was attributed  to a proximity-mediated enhancement of the superconductivity in a granular material in a strong enough electric field. Explicitly, the critical current in a $`YBCO`$ sample was found to reach a maximum at $`E=4\times 10^7V/m`$. To relate this experimental value with the model parameters, first of all, we need to estimate an order of magnitude of the relaxation time $`\tau `$ in a zero applied electric field. This will provide an upper limit for the relevant $`\tau `$-distribution in our system. It is reasonable to connect zero-field $`\tau \tau (0)`$ with the Josephson tunneling time  $`\tau _J=(R_0/R_n)(\mathrm{}/J)`$ (where $`R_0=h/4e^2`$, and $`R_n`$ is the normal state resistance between grains). Typically, for HTCS ceramics $`J/k_B90K`$ and $`R_n/R_010^3`$, so that $`\tau _J10^{10}s`$. At the same time, at high enough electric fields where the MEE becomes strongly nonlinear, we can expect quite a tangible decrease of the relaxation time. Indeed, for an average grain size $`d1\mu m`$, the characteristic field $`E_0=4\times 10^7V/m`$ (which corresponds to the region where a prominent enhancement of the critical current was observed \[17-19\]) introduces a substantially shorter relaxation time $`\tau (E_0)=\mathrm{}/2deE_010^{16}s`$, in agreement with observations.

To estimate the relative magnitude of the superconducting analog of the MEE predicted here, we can compare it with the normal (Ohmic) contribution to the magnetization

$$\stackrel{}{M}_n(\stackrel{}{E})\overline{\stackrel{}{\mu }_n(t)}=\frac{1}{\tau }\underset{0}{\overset{\tau }{}}𝑑t\underset{0}{\overset{\mathrm{}}{}}𝑑\stackrel{}{r}_{i,j}𝑑\stackrel{}{R}_{i,j}f(\stackrel{}{r}_{i,j},\stackrel{}{R}_{i,j})\stackrel{}{\mu }_n,$$
(13)
where

$$\stackrel{}{\mu }_n=\underset{i,j}{}I_{i,j}^n(\stackrel{}{r}_{i,j}\stackrel{}{R}_{i,j})$$
(14)
Here $`I_{i,j}^n=V_{i,j}/R_n`$ is the normal current component due to the applied electric field $`\stackrel{}{E}`$, with $`V_{i,j}=\stackrel{}{E}\stackrel{}{r}_{i,j}`$ being the induced voltage, and $`R_n`$ the normal state resistance between grains. As a result, the normal state contribution (for $`\stackrel{}{E}`$ along the $`x`$-axis) reads $`M_n=\alpha _nE_x`$, with $`\alpha _n=\pi d^2DN/R_n`$. Similarly, according to Eq.(10), the low field contribution to the superconducting MEE gives $`M_s\alpha _sE_x`$ with $`\alpha _s=2\pi e^2JN\tau (0)d^2D/\mathrm{}^2`$. Thus, at low enough applied fields (when $`E_xE_0`$)

$$\frac{\alpha _s}{\alpha _n}\frac{\tau (0)}{\tau _J},$$
(15)
where $`\tau _J=(R_0/R_n)(\mathrm{}/J)`$ with $`R_0=h/4e^2`$. According to our previous discussion on the relevant relaxation-time distribution spectrum in our model system, we may conclude that $`\tau (0)\tau _J`$. So, we arrive at the following ratio between the coefficients of the superconducting to normal MEEs, namely $`\alpha _s/\alpha _n1`$. The above estimate of the weak-links induced MEE (along with its rather specific field dependence, see Fig.1) suggests quite an optimistic possibility to observe the predicted effect experimentally in HTCS ceramics or in a specially prepared system of arrays of superconducting grains.

We note that in the present analysis we have not explicitly considered the polarization effects (due to the interaction between the applied electric field and the grain’s charges) which may become important at high enough fields (or for small enough grains), leading to more subtle phenomena (like Coulomb blockade and reentrant-like behavior) that will demand the inclusion of charging energy effects in the analysis.

3. MAGNETIC FIELD INDUCED CHARGING EFFECTS

This Section addresses a related phenomenon which is actually dual to the above-discussed analog of magnetoelectric effect. Specifically, we analyze a possible appearance of a non-zero electric polarization and the related change of the charge balance in the system of weakly-coupled superconducting junctions (modelled by the random $`3D`$ JJAs) under the influence of an applied magnetic field. Since the field-induced effects considered in this Section are expected to manifest themselves in high enough applied magnetic fields (with a nearly homogeneous distribution of magnetic flux along the junctions) and as long as the Josephson penetration length $`\lambda _J`$ exceeds the characteristic size of the Josephson network $`d`$ (which is related to the projected junction area $`S`$, where the field penetration actually occurs, as follows $`S=\pi d^2`$), the Josephson current-induced ”self-field” effects (which are important for a large-size junctions and/or small applied magnetic fields) may be safely neglected . Besides, it is known  that in discrete JJAs pinning (by a single junction) actually concurs with the ”self-field” effects. Specifically, it was found  that the ratio $`d/\lambda _J`$ is related to the dimensionless pinning strength parameter $`\beta `$ as $`d/\lambda _J=\sqrt{\beta }`$ suggesting that a weak pinning regime (with $`\beta 1`$) simultaneously implies a smallness of the ”self-field” effect and vice versa. And since artificially prepared Josephson networks allow for more flexibility in varying the experimentally-controlled parameters, it is always possible to keep the both above effects down by appropriately tuning the ratio $`d/\lambda _J`$. Typically, in this kind of experiments $`S=0.010.1\mu m^2`$ and $`d\lambda _J`$.

In what follows, we are interested in the magnetic field induced behavior of the electric polarization in a $`3D`$ JJA at zero temperature. Recall that a conventional (zero-field) pair polarization operator within the model under discussion reads

$$\stackrel{}{p}=\underset{i=1}{\overset{N}{}}q_i\stackrel{}{r}_i,$$
(16)
where $`q_i=2en_i`$ with $`n_i`$ the pair number operator, and $`r_i`$ is the coordinate of the center of the grain.

In view of Eqs.(1)-(6), and taking into account a usual ”phase-number” commutation relation, $`[\varphi _i,n_j]=i\delta _{i,j}`$, the evolution of the pair polarization operator is determined via the equation of motion

$$\frac{d\stackrel{}{p}}{dt}=\frac{1}{i\mathrm{}}[\stackrel{}{p},]=\frac{2e}{\mathrm{}}\underset{ij}{\overset{N}{}}J\mathrm{sin}\varphi _{i,j}(\stackrel{}{H})\stackrel{}{r}_{i,j}$$
(17)
Resolving the above equation, we arrive at the following net value of the magnetic-field induced polarization (per grain)

$$\stackrel{}{P}(\stackrel{}{H})\frac{1}{N}\overline{<\stackrel{}{p}(t)>}=\frac{2eJ}{\mathrm{}\tau N}\underset{0}{\overset{\tau }{}}𝑑t\underset{0}{\overset{t}{}}𝑑t^{}\underset{ij}{\overset{N}{}}<\mathrm{sin}\varphi _{i,j}(\stackrel{}{H})\stackrel{}{r}_{i,j}>,$$
(18)
where $`<\mathrm{}>`$ denotes a configurational averaging over the grain positions, while the bar means a temporal averaging (with a characteristic time $`\tau `$). To consider a field-induced polarization only, we assume that in a zero magnetic field, $`\stackrel{}{P}0`$, implying $`\varphi _{i,j}(0)0`$.

Since, as usual, the main qualitative results presented here do not depend on the particular choice of the probability distribution function, for a change, the following law $`f(\stackrel{}{r},\stackrel{}{R})=f_r(\stackrel{}{r})f_R(\stackrel{}{R})`$ with $`f_r(\stackrel{}{r})=\frac{1}{(2\pi d^2)^{3/2}}e^{(x^2+y^2+z^2)/2d^2}`$ and $`f_R(\stackrel{}{R})=\delta (Xd)\delta (Yd)\delta (Zd)`$ will be assumed. As a result, we observe that the magnetic field $`\stackrel{}{H}=(0,0,H_z)`$ (applied along the $`z`$-axis) will induce a non-vanishing longitudinal (along $`x`$-axis) electric polarization

$$P_x(H_z)=P_0G(H_z/H_0),$$
(19)
with

$$G(z)=ze^{z^2}$$
(20)
Here $`P_0=ed\tau J/\mathrm{}`$, $`H_0=\mathrm{\Phi }_0/S`$ with $`S=\pi d^2`$ being an average projected area of a single junction, and $`z=H_z/H_0`$.

For small applied fields ($`z1`$), the induced polarization $`P_x(H_z)\alpha _0H_z`$, where $`\alpha _0=e\tau dJ/\mathrm{}H_0`$ is the so-called linear magnetoelectric coefficient . However, as we mentioned in the very beginning, to correctly describe any induced effects in very small external fields, both Josephson junction pinning and ”self-field” effects have to be taken into account. At the same time, in view of Eq.(16) the induced polarization is related to the corresponding change of the effective charge $`\delta Q(1/d)_i<q_ix_i>`$ in applied magnetic field as follows

$$\delta Q(H_z)=\frac{P_x(H_z)}{d}=Q_0G(H_z/H_0),$$
(21)
where $`Q_0=e\tau J/\mathrm{}`$.

It is of interest also to consider the related field behavior of the effective flux capacitance $`\delta C\tau d\delta Q(H_z)/d\mathrm{\Phi }`$ which in view of Eq.(21) reads

$$\delta C(H_z)=C_0\left(12\frac{H_z^2}{H_0^2}\right)e^{H_z^2/H_0^2},$$
(22)
where $`\mathrm{\Phi }=SH_z`$, and $`C_0=\tau Q_0/\mathrm{\Phi }_0`$. Figure 2 shows the behavior of the induced effective charge $`\delta Q/Q_0`$ in applied magnetic field $`H_z/H_0`$. As is seen, at $`H_z/H_00.8`$ the effective charge reaches its maximum (while the capacitance changes its sign at this field), suggesting a significant redistribution of the junction charge balance in a model system under the influence of an applied magnetic field, near the Josephson critical field $`H_0`$. Note that a somewhat similar behavior of the magnetic field induced charge (and related capacitance) has been observed in 2D electron systems .

Taking $`\tau =10^{10}s`$ for the Josephson relaxation time (which is related to the Josephson plasma frequency, $`\omega _p\tau ^1`$, known to be the characteristic frequency of the system for $`E_CJ`$ regime, with $`E_C`$ being the Coulomb grain’s charge energy; typically $`\omega _p=10^910^{11}Hz`$) and $`J/k_B=90K`$ for a zero-temperature Josephson energy in $`YBCO`$ ceramics, we arrive at the following estimates of the effective charge $`Q_010^{16}C`$, flux capacitance $`C_010^{11}F`$, the equivalent current $`I_0Q_0/\tau 10^6A`$, and voltage $`V_0Q_0/C_010^5V`$. We note that the above set of estimates fall into the range of parameters used in typical experiments to study the charging effects both in single JJs (with the working frequency from RF range of $`\omega 10GHz`$ used to stimulate the system ) and JJAs  suggesting thus quite an optimistic possibility to observe the above-predicted field induced effects experimentally, using a specially prepared system of arrays of superconducting grains.

4. STRESS INDUCED PARAMAGNETIC MEISSNER EFFECT

The possibility to observe tangible piezoeffects in mechanically loaded grain boundary Josephson junctions (GBJJs) is based on the fact that under plastic deformation, grain boundaries (GBs) (which are the natural sources of weak links in HTCS), move rather rapidly via the movement of the grain boundary dislocations (GBDs) comprising these GBs \[38-40\]. Using the above evidence, in Ref.8 a piezophase response of a single GBJJ (created by GBDs strain field $`ϵ_d`$ acting as an insulating barrier of thickness $`l`$ and height $`U`$ in a $`SIS`$-type junction with the Josephson energy $`Je^{l\sqrt{U}}`$) to an applied stress was considered. To understand how piezoeffects manifest themselves through GBJJs, let us invoke an analogy with the so-called thermophase effect  (a quantum mechanical alternative for the conventional thermoelectric effect) in JJs. In essence, the thermophase effect assumes a direct coupling between an applied temperature drop $`\mathrm{\Delta }T`$ and the resulting phase difference $`\mathrm{\Delta }\varphi `$ through a JJ. When a rather small temperature gradient is applied to a JJ, an entropy-carrying normal current $`I_n=L_n\mathrm{\Delta }T`$ (where $`L_n`$ is the thermoelectric coefficient) is generated through such a junction. To satisfy the constraint dictated by the Meissner effect, the resulting supercurrent $`I_s=I_c\mathrm{sin}[\mathrm{\Delta }\varphi ]`$ (with $`I_c=2eJ/h`$ being the Josephson critical current) develops a phase difference through a weak link. The normal current is locally canceled by a counterflow of supercurrent, so that the total current through the junction $`I=I_n+I_s=0`$. As a result, supercurrent $`I_c\mathrm{sin}[\mathrm{\Delta }\varphi ]=I_n=L_n\mathrm{\Delta }T`$ generates a nonzero phase difference leading to the linear thermophase effect  $`\mathrm{\Delta }\varphi L_{tp}\mathrm{\Delta }T`$ with $`L_{tp}=L_n/I_c(T)`$.

By analogy, we can introduce a piezophase effect (as a quantum alternative for the conventional piezoelectric effect) through a JJ . Indeed, a linear conventional piezoelectric effect relates induced polarization $`P_n`$ to an applied strain $`ϵ`$ as  $`P_n=d_nϵ`$, where $`d_n`$ is the piezoelectric coefficient. The corresponding normal piezocurrent density is $`j_n=dP_n/dt=d_n\dot{ϵ}`$ where $`\dot{ϵ}(\sigma )`$ is a rate of plastic deformation which depends on the number of GBDs of density $`\rho `$ and a mean dislocation rate $`v_d`$ as follows  $`\dot{ϵ}(\sigma )=b\rho v_d(\sigma )`$ (where $`b`$ is the absolute value of the appropriate Burgers vector). In turn, $`v_d(\sigma )v_0(\sigma /\sigma _m)`$. To meet the requirements imposed by the Meissner effect, in response to the induced normal piezocurrent, the corresponding Josephson supercurrent of density $`j_s=dP_s/dt=j_c\mathrm{sin}[\mathrm{\Delta }\varphi ]`$ should emerge within the contact. Here $`P_s=2enb`$ is the Cooper pair’s induced polarization with $`n`$ the pair number density, and $`j_c=2ebJ/\mathrm{}V`$ is the critical current density. The neutrality conditions ($`j_n+j_s=0`$ and $`P_n+P_s=const`$) will lead then to the linear piezophase effect $`\mathrm{\Delta }\varphi d_{pp}\dot{ϵ}(\sigma )`$ (with $`d_{pp}=d_n/j_c`$ being the piezophase coefficient) and the concomitant change of the pair number density under an applied strain, viz., $`\mathrm{\Delta }n(ϵ)=d_{pn}ϵ`$ with $`d_{pn}=d_n/2eb`$.

To adequately describe magnetic properties of a granular superconductor, we again employ a model of random three-dimensional (3D) overdamped Josephson junction array which is based on the familiar tunneling Hamiltonian given by Eq.(1) (see Section 2). According to the above-discussed scenario, under mechanical loading the superconducting phase difference will acquire an additional contribution $`\delta \varphi _{i,j}(\sigma )=B\stackrel{}{\sigma }\stackrel{}{r}_{i,j}`$, where $`B=d_n\dot{ϵ}_0/\sigma _mj_cb`$ with $`\dot{ϵ}_0=b\rho v_0`$ being the maximum deformation rate and the other parameters defined earlier. If, in addition to the external loading, the network of superconducting grains is under the influence of an applied frustrating magnetic field $`\stackrel{}{H}`$, the total phase difference through the contact reads

$$\varphi _{i,j}(\stackrel{}{H},\stackrel{}{\sigma })=\varphi _{i,j}^0+\frac{\pi }{\mathrm{\Phi }_0}(\stackrel{}{r}_{i,j}\stackrel{}{R}_{i,j})\stackrel{}{H}B\stackrel{}{\sigma }\stackrel{}{r}_{i,j}$$
(23)
Once again, to neglect the influence of the self-field effects in a real material, the corresponding Josephson penetration length $`\lambda _J`$ must be much larger than the junction (or grain) size. Likewise, to ensure the uniformity of the applied stress $`\sigma `$, we also assume that $`d\lambda _\sigma `$, where $`\lambda _\sigma `$ is a characteristic length over which $`\sigma `$ is kept homogeneous.

When the Josephson supercurrent $`I_{i,j}^s=I_c\mathrm{sin}\varphi _{i,j}`$ circulates around a set of grains, it induces a random magnetic moment $`\stackrel{}{\mu }_s`$ of the Josephson network which results in the stress induced net magnetization (Cf. Eq.(8))

$$\stackrel{}{M}_s(\stackrel{}{H},\stackrel{}{\sigma })<\stackrel{}{\mu }_s>=\underset{0}{\overset{\mathrm{}}{}}𝑑\stackrel{}{r}_{i,j}𝑑\stackrel{}{R}_{i,j}f(\stackrel{}{r}_{i,j},\stackrel{}{R}_{i,j})\stackrel{}{\mu }_s$$
(24)
To capture the very essence of the superconducting piezomagnetic effect, in what follows we assume for simplicity that an unloaded sample does not possess any spontaneous magnetization at zero magnetic field (that is $`M_s(0,0)=0`$) and that its Meissner response to a small applied field $`H`$ is purely diamagnetic (that is $`M_s(H,0)H`$). According to Eq.(23), this condition implies $`\varphi _{i,j}^0=2\pi m`$ for the initial phase difference with $`m=0,\pm 1,\pm 2,..`$.

Taking the applied stress along the $`x`$-axis, $`\stackrel{}{\sigma }=(\sigma ,0,0)`$, normally to the applied magnetic field $`\stackrel{}{H}=(0,0,H)`$, and assuming an exponential distribution law for the distance between grains, $`f_r(\stackrel{}{r})=f(x)f(y)f(z)`$ with $`f(x_j)=(1/d)e^{x_j/d}`$ (see Section 2 for detailes), we get finally

$$M_s(H,\sigma )=M_0(\sigma )\frac{H_{tot}(H,\sigma )/H_0}{[1+H_{tot}^2(H,\sigma )/H_0^2]^2},$$
(25)
for the induced transverse magnetization, where $`H_{tot}(H,\sigma )=HH^{}(\sigma )`$ is the total magnetic field with $`H^{}(\sigma )=[\sigma /\sigma _0(\sigma )]H_0`$ being a stress-induced contribution. Here, $`M_0(\sigma )=I_c(\sigma )SN/V`$ with $`S=\pi dD`$ being a projected area around the JJ, $`H_0=\mathrm{\Phi }_0/S`$, and $`\sigma _0(\sigma )=\sigma _m[j_c(\sigma )/j_d](b/d)`$ with $`j_d=d_n\dot{ϵ}_0`$ and $`\dot{ϵ}_0=b\rho v_0`$ being the maximum values of the dislocation current density and the plastic deformation rate, respectively.

According to the recent experiments , the tunneling dominated critical current $`I_c`$ (and its density $`j_c`$) in HTCS ceramics was found to exponentially increase under compressive stress, viz. $`I_c(\sigma )=I_c(0)e^{\beta \sigma }`$ with $`\beta 1/\sigma _m`$. Specifically, the critical current at $`\sigma =9kbar`$ was found to be three times higher its value at $`\sigma =1.5kbar`$, clearly indicating a weak-links-mediated origin of the phenomenon (in the best defect-free thin films this ratio is controlled by the stress induced modifications of the carrier number density and practically never exceeds a few percents ). Strictly speaking, the critical current will also change (decrease) with applied magnetic field. However, for the fields under discussion (see below) in the first approximation this effect can be neglected. In view of Eq.(25), dependence of $`I_c`$ on $`\sigma `$ will lead to a rather strong piezomagnetic effects. Indeed, Fig.3 shows changes of the initial stress-free diamagnetic magnetization $`M_s/M_c`$ (solid line) under an applied stress $`\sigma /\sigma _c`$. Here $`M_cM_0(0)`$ and $`\sigma _c\sigma _0(0)`$ (see below for estimates). As we see, already relatively small values of an applied stress render a low field Meissner phase strongly paramagnetic (dotted and dashed lines) simultaneously increasing the maximum of the magnetization and shifting it towards higher magnetic fields. According to Eq.(25), the initially diamagnetic Meissner effect turns paramagnetic as soon as the piezomagnetic contribution $`H^{}(\sigma )`$ exceeds an applied magnetic field $`H`$. To see whether this can actually happen in a real material, let us estimate a magnitude of the piezomagnetic field $`H^{}`$. Typically, for HTCS ceramics $`S10\mu m^2`$, leading to $`H_01G`$. To estimate the needed value of the dislocation current density $`j_d`$, we turn to the available experimental data. According to Ref.39, a rather strong polarization under compressive pressure $`\sigma /\sigma _m0.1`$ was observed in $`YBCO`$ ceramic samples at $`T=77K`$ yielding $`d_n=10^2C/m^2`$ for the piezoelectric coefficient. Usually, for GBJJs $`\dot{ϵ}_010^2s^1`$, and $`b10nm`$ leading to $`j_d=d_n\dot{ϵ}_01A/m^2`$ for the maximum dislocation current density. Using the typical values of the critical current density $`j_c(\sigma )=10^4A/m^2`$ (found  for $`\sigma /\sigma _m0.1`$) and grain size $`d1\mu m`$, we arrive at the following estimate of the piezomagnetic field $`H^{}10^2H_0`$. Thus, the predicted stress induced paramagnetic Meissner effect should be observable for applied magnetic fields $`H10^2H_00.01G`$ which correspond to the region where the original PME was first registered \[1-3\]. In turn, the piezoelectric coefficient $`d_n`$ is related to a charge $`Q`$ in the GBJJ as  $`d_n=(Q/S)(d/b)^2`$. Given the above-obtained estimates, we get $`Q10^{13}C`$ for an effective charge accumulated by the GBs. Notice that the above values of the aplied stress $`\sigma `$ and the resulting effective charge $`Q`$ correspond (via the so-called electroplastic effect ) to an equivalent applied electric field $`E=b^2\sigma /Q10^7V/m`$ at which rather pronounced electric-field induced effects in HTCS have been recently observed \[17-21\].

Besides, according to Ref.43 the Josephson projected area $`S`$ slightly decreases under pressure thus leading to some increase of the characteristic field $`H_0=\mathrm{\Phi }_0/S`$. In view of Eq.(25), it means that a smaller compression stress is needed to actually reverse the sign of the induced magnetization $`M_s`$. Furthermore, if an unloaded granular superconductor already exhibits the PME, due to the orbital currents induced spontaneous magnetization resulting from an initial phase difference $`\varphi _{i,j}^0=2\pi r`$ in Eq.(23) with fractional $`r`$ (in particular, $`r=1/2`$ corresponds to the so-called  $`\pi `$-type state), then according to our scenario this effect will either be further enhanced by applying a compression (with $`\sigma >0`$) or will disappear under a strong enough extension (with $`\sigma <0`$). Given the markedly different mechanisms and scales of stress induced changes in defect-free thin films  and weak-links-ridden ceramics , it should be possible to experimentally register the suggested here piezophase effects.

ACKNOWLEDGMENTS

This work was conceived and partially done during my stay at the Universidade Federal de São Carlos (Brazil) where it was funded by the Brazilian Agency FAPESP (Projeto 2000/04187-8). I thank Wilson Ortiz and Fernando Araujo-Moreira for hospitality and stimulating discussions on the subject. I am also indebted to Anant Narlikar for his invitation to make this contribution for the current volume of the Studies.


format: text

text:
# Non-pulsing emission from X-ray pulsars

## 1. Her X-1

The “off” state of the X-ray pulsar Her X-1 was observed with the telescope ART-P on board GRANAT three times in 1990-1991. The source persistent flux in this state was detected at the level of few mCrabs (several percents of its high-state flux) and was characterized by the absence of X-ray pulsations — the $`3\sigma `$ upper limit on the pulse fraction in the 3-20 keV band varied between 2.9% and 18% in different observations. The source’s “off” state spectra are well described by a simple power law without evidence of a cutoff in the hard energy band or of emission lines. The spectral slopes differ by a factor of $`2`$, with the photon flux in the energy band 3-20 keV being virtually the same. Her X-1 is observed in its “off” state if the compact source is hidden behind the outer edge of a warped (or tilted to the orbital plane) accretion disk and, as it is believed, the source emission can greatly scattered in a hot corona (it’s height is $`H10^{11}`$ cm) above the outer parts of the disk. Such scattering of emission from the compact source by coronal electrons can account for the existence of an X-ray flux in the “off” state and it allows to explain the absence of pulsations in this state, as any information about variability on a time scale shorter than $`H/c3`$ s, where $`c`$ is the speed of light, must be lost.

## 2. Cen X-3

A significant declination (by a factor of 6) in the 4-20 keV source flux was detected during the observation on Aug 19, 1990. This declination occurred in two steps, latter of which corresponded to an orbital phase $`\varphi 0.89`$ at which the X-ray eclipse ingress was observed in other experiments. It is important to note that the intensity decline during the eclipse ingress mainly occurred in the hard energy band and the corresponding hardness (ratio of the 10-20 keV and 4-6 keV fluxes) decrease was equal to $`40`$%. After the eclipse ingress a non-zero and non-pulsing X-ray emission was detected by ART-P from the source at the signal-to-noise level of $`S/N13\sigma `$. The pulsar’s spectrum during the eclipse was measured with large statistical errors; nevertheless, a line at energy of $`7`$ keV was very reliable. By assuming that this is the 6.7-keV line of helium-like iron, we estimated its intensity, $`I_{6.7}(2.7\pm 1.4)\times 10^3`$phot$``$cm-2s-1 , and equivalent width, $`EW_{6.7}(1.2\pm 0.6)`$ keV. Although the source itself was already occulted by the disk of the optical star at this time, part of the surface of the scattering cloud remained in the visibility zone, and we observed precisely this emission scattered in the cloud. Note also that the measured spectrum at energies below 5 keV shows an excess of soft X-ray emission, which was previously revealed from the light curve analysis.

## 3. Vela X-1

The observation of this X-ray pulsar on June 15, 1992 (orbital phases 0.918–0.953) was carried out mainly during an eclipse. We divided it into two parts, in accordance with the beginning of the source’s eclipse egress. A non-zero flux of approximately 1/10 of the bright-state flux at the confidence levels 4.3 and 4.6$`\sigma `$ was observed from the source during the eclipse and eclipse egress, respectively. During the second part of the observation X-ray pulsations with the period $`P=283.33\pm 1.65`$ s were detected from the source at a confidence level of $`10.2\sigma `$. The pulse fraction was equal to $`33.9\pm 6.4`$% and it was a factor of 2 lower than one measured in the source’s bright state several days earlier. During the eclipse we failed to detect any pulsations — the 3$`\sigma `$ limit on the pulse fraction was 9.6%. The persistent flux, detected during an X-ray eclipse, may be the emission scattered in a fairly dense stellar wind or in the extended atmosphere of the optical star. Given that the binary characteristic size, $`la/c120`$ light seconds, is half the pulsation period, photons must undergo several scatterings for the pulse profile to be smeared enough. For the scattering optical depth of the wind $`\tau _\mathrm{T}1`$, the scattered-to-incident flux ratio allows us to estimate the solid angle $`\mathrm{\Omega }0.4\pi a^2`$ at which the scattering-envelope ring is seen from the compact source. Accordingly, the envelope height above the optical-star surface is $`H14R_{}.`$

More comprehensive discussion of the results of the observations non-pulsing emission from X-ray pulsars with the ART-P telescope are presented elsewhere (see Lutovinov et al. 1999, 2000a,b).

## References

Lutovinov, A.A., Grebenev, S.A., & Sunyaev, R.A. 1999, Astron. Lett., 25, 59.

Lutovinov, A.A., Grebenev, S.A., Pavlinsky, M.N., & Sunyaev, R.A. 2000a, Astron. Lett., 26, 691.

Lutovinov, A.A., Grebenev, S.A., Pavlinsky, M.N., & Sunyaev, R.A. 2000b, Astron. Lett., 26, 765.


text:
# Singularity Formation in 2+1 Wave Maps

## I Introduction

While it has been shown that wave maps on a $`1+1`$ dimensional Minkowski spacetime base evolved from smooth initial data exist for all time , and that those on an $`m+1`$ $`(m3)`$ Minkowski spacetime base can blow up in finite time , global existence for the $`2+1`$ case remains as yet unresolved. Scaling considerations identify $`2+1`$ as the critical dimension for wave maps, and so there is considerable interest in determining if indeed $`2+1`$ wave maps developed from smooth initial data can become singular in finite time or not. Here, we describe numerical work which strongly supports the contention that, at least for some sets of smooth initial data, they can.

There are special classes of $`2+1`$ wave maps for which global existence has been shown to hold: a) spherically equivariant wave maps with convex , or slightly more general targets , b) spherically symmetric wave maps with compact targets (plus a further technical condition on the target) , c) general wave maps (general target) with sufficiently small energy.

Not included in any of these three classes are spherically equivariant wave maps from $`2+1`$ Minkowski spacetime into the round two-sphere with initial data of arbitrary energy. Shatah and Struwe  have conjectured that singular behavior should be found in this class. Our numerical results reported here strongly support the validity of this conjecture.

We examine one-parameter families of data, with small values of the parameter corresponding to small energy data and therefore global existence, and with large values of the parameter corresponding to data possibly leading to singularity formation. One might hope to find especially interesting wave map development for data at or near the transition between small and large values. While this sort of “critical” behavior has been seen and studied in $`3+1`$ wave maps , we have not found nearly as clear an indication of universal critical behavior for the present $`2+1`$ case. This criticality issue needs further study, and is not treated in this paper. Here, our focus is on numerical evidence for singular wave map evolution from regular initial data.

We note that our studies of singularity formation in $`2+1`$ wave maps have been carried out independently of the work of Bizoń, Chmaj, and Tabor  using numerical algorithms which differ from theirs. However their results and ours agree substantially.

## II The Equations

Generally a wave map is defined to be a map $`\varphi ^A`$ from a spacetime (the “base”) into a Riemannian geometry (the “target”), with $`\varphi ^A`$ a critical point for the action

$$S[\varphi ]=_{M^{m+1}}\eta ^{\mu \nu }g_{AB}\left(\varphi \right)\left(_\mu \varphi ^A_\nu \varphi ^B\right)$$
(1)
where $`g_{AB}`$ is the Riemannian metric on the target manifold $`N^n`$, and $`\eta ^{\mu \nu }`$ is the (inverse) Lorentz-signature metric on the spacetime $`M^{m+1}`$. The Euler-Lagrange equations for this action take the form

$$^\mu _\mu \varphi ^A+\mathrm{\Gamma }_{BC}^A_\mu \varphi ^B^\mu \varphi ^C=0$$
(2)
where $`\mathrm{\Gamma }_{BC}^A`$ represents the Christoffel symbols corresponding to the target metric $`g_{AB}`$. This is a semilinear hyperbolic PDE system for $`\varphi ^A`$. We note that for certain targets, wave maps are known to physicists as “nonlinear sigma models.”

As noted above, the case of primary interest here is $`2+1`$ Minkowski spacetime for the base and the round two sphere for the target. In this case, the wave map PDE system (2) may be rewritten in the following form

$$\mathrm{}\varphi ^a+\left(_\mu \varphi ^b^\mu \varphi ^c\right)\delta _{bc}\varphi ^a=0.$$
(3)
where the indices $`a,b,c`$ take the values $`\{1,2,3\}`$ (indexing the ambient Euclidean 3-space for the target two sphere), and $`\delta _{bc}`$ is the metric for this ambient space. If we now impose the condition that the maps $`\varphi ^a`$ be spherically equivariant with angular wrapping number k, and write $`\varphi ^a(r,\theta ,t)`$ in the “hedgehog” form

$$\varphi ^a=\left(\begin{array}{c}\mathrm{sin}\chi (r,t)\mathrm{sin}k\theta \\ \mathrm{sin}\chi (r,t)\mathrm{cos}k\theta \\ \mathrm{cos}\chi (r,t)\end{array}\right)$$
(4)
where $`r`$ is the radial distance from the origin and $`\theta `$ is the azimuthal angle, then the wave map PDE system (2) reduces to the single equation

$$\ddot{\chi }=\frac{1}{r}\left(r\chi ^{}\right)^{}\frac{k^2\mathrm{sin}2\chi }{2r^2}$$
(5)
where a prime and a dot denote partial derivatives with respect to $`r`$ and $`t`$ respectively. Thus the study of the Cauchy problem for $`2+1`$ spherically equivariant ($`k`$-wrapped) wave maps into the round two sphere focuses on finding solutions $`\chi (r,t)`$ to Eq. (5) with regular initial data $`\chi (r,0)`$, $`\dot{\chi }(r,0)`$. Note that regularity at $`r=0`$ requires that we set $`\chi (0,t)=0`$ for all $`t`$.

While it may be interesting to examine if there is any variation of the behavior of solutions for wrapping numbers $`k`$ greater than one, we restrict our attention here to the single angular wrapping case $`k=1`$.

As for any field theory on Minkowski space, there is a divergence-free stress-energy tensor $`T_{\mu \nu }`$ associated with wave maps. From $`T_{\mu \nu }`$, we obtain the energy density function for spherically equivariant wave maps

$$\rho (r,t)=\frac{1}{2}\left[\dot{\chi }^2+\left(\chi ^{}\right)^2\right]+\frac{\mathrm{sin}^2\chi }{2r^2}$$
(6)
whose integral

$$E(t)=_r\rho (r,t)r𝑑r$$
(7)
is conserved (ie, $`E(t)=E(0)`$ for all $`t`$). The energy is a useful monitor of numerical accuracy, as discussed below.

## III Numerical Studies of Singularity Formation

Our numerical experiments consist of specifying parametrized families of initial data $`\{\chi _\lambda (r,0),\dot{\chi }_\lambda (r,0)\}`$ and numerically evolving a number of sets of such data in each family. A typical family–one of the simplest–is the approximately ingoing Gaussian pulse

$`\chi (r,0)`$ $`=`$ $`Ae^{\left(rR_0\right)^2/\delta ^2}`$ (8)
$`\dot{\chi }(r,0)`$ $`=`$ $`\chi ^{}(r,0).`$ (9)
This family has three parameters $`A`$, $`R_0`$, and $`\delta `$, with the most important one for our discussion being the scale parameter $`A`$. Note that the ingoing character of these solutions, which results from the choice of $`\dot{\chi }(r,0)`$, minimizes outer boundary effects. Note also that while, analytically, $`\chi (0,0)`$ is not zero, for the choices of $`R_0`$ and $`\delta `$ which we make, we can force $`\chi (0,0)`$ to be zero and retain smoothness to within numerical accuracy.

We evolve using a second order finite difference approximation to Eq. (5). We use an iterative Crank-Nicholson scheme implemented with RNPL , and also make use of the adaptive mesh framework developed by Choptuik . We have verified that the code generates solutions which converge quadratically in the grid spacing and conserve energy. In arguing that we are indeed generating singularities, we will discuss the convergence and energy conservation tests in more detail below.

For a general set of ingoing Gaussian pulse data, regardless of amplitude, the wave map evolution has the pulse maximum and energy density maximum initially moving inward (decreasing $`r`$). For small (subcritical) values of $`A`$, this inward motion of the maximum proceeds for a finite time, after which the maximum “bounces” away from the origin and begins to move outward (see Fig. 1). There is a general dispersal of the energy density; and for large $`t`$, there is very little energy density remaining near the origin.

For large (supercritical) values of $`A`$, the behavior of the evolving wave map is qualitatively the same initially. However, rather than bouncing away from the origin, the maxima for supercritical data continue to approach the origin (Fig. 1), with the concentration of energy around the origin appearing to grow without bound. As the energy density and the gradient of the function $`\chi `$ grow very large at the origin, the numerical evolution inevitably becomes unable to resolve the gradient, and the solution becomes sufficiently non-smooth to cause the numerical evolution to stop. If this accumulation is indeed a singularity forming, there is no hope for the numerical evolution to resolve it, being itself of finite resolution. The task then is to examine the behavior of the numerical solution up to this point.

Before doing so, we first discuss a couple of standard tests of a numerical solution. We let $`\chi (r,t)`$ be some solution to the (continuum) partial differential equation (5) and let $`\stackrel{~}{\chi }_h(r,t)`$ be the solution to a discrete form of that equation, for corresponding initial data, on a grid spacing $`h\mathrm{\Delta }r`$. The hope is that, as the grid spacing $`\mathrm{\Delta }r`$ gets smaller, the solutions to the discrete equation generated by the evolution code converge to the solutions of the PDE, $`\stackrel{~}{\chi }(r,t)\chi (r,t)`$. Because in general the explicit solutions to the PDE are unknown, we instead consider a series of numerical solutions on grids of increasing resolution, say $`\stackrel{~}{\chi }_{4h}`$, $`\stackrel{~}{\chi }_{2h}`$, $`\stackrel{~}{\chi }_h`$. If these are to converge to the PDE solution, then they must converge themselves. To examine this convergence, we define a convergence factor ($`Q`$) as follows

$$\text{Q}\frac{|\stackrel{~}{\chi }_{4h}\stackrel{~}{\chi }_{2h}|_2}{|\stackrel{~}{\chi }_{2h}\stackrel{~}{\chi }_h|_2},$$
(10)
where the norms are the $`l_2`$ norm. For these solutions to converge, the difference between solutions for increasing resolution must decrease and hence $`Q`$ must be greater than one. For second-order schemes, $`Q`$ is expected to be $`4`$.

Another common test of numerical accuracy focuses on the degree to which energy is conserved by the numerical evolution. The evolution governed by the PDE (5) does conserve energy; the question is whether this remains true for the numerical evolution. Letting $`E_{num}(t)`$ denote the energy calculated from the numerical solution at time $`t`$ (on the finite grid), and setting $`\mathrm{\Delta }(t)\mathrm{ln}\left|\frac{E_{num}(t)E_{num}(0)}{E_{num}(0)}\right|`$, we monitor $`\mathrm{\Delta }(t)`$ for different choices of grid spacing. The expectation is that $`\mathrm{\Delta }(t)`$ should decrease with increasing resolution; if we observe this, our confidence in the accuracy of our numerical solution is enhanced.

In Fig. 2, we show the evolution in time of three quantities–$`\mathrm{ln}\chi ^{}(0,t)`$, $`Q(t)`$, and $`\mathrm{\Delta }(t)`$–for numerical runs of supercritical ingoing Gaussian pulse data, done with nine different grid spacings. In the top frame, we show the behavior of the derivative of $`\chi `$ at the origin as a function of time. The figure shows that as the pulse travels inward, the derivative increases. Until just before $`t=8`$, all the resolutions show the same behavior as would be expected for a convergent evolution. However, near the blow-up time, the solutions diverge with higher resolutions providing a larger derivative. The convergence factor $`Q(t)`$ is shown in the middle frame; it likewise shows second-order convergence up to times close to the blow-up time. In the bottom frame, the change in energy $`\mathrm{\Delta }(t)`$ is shown. We see that as the resolution is increased, energy conservation improves.

What does this tell us about singularity formulation in wave maps evolved from (supercritical) ingoing Gaussian pulse data? We first argue that these results are consistent with what would be expected for such formation. As the singularity forms, higher and higher frequency components become important, and they are represented numerically only if one uses higher and higher grid resolutions. Hence, the behavior of the derivative of $`\chi `$ as the resolution improves would be expected to show larger and larger gradients, as seen in Figure 2a. Next, we note that the formation of a singularity should not hinder convergence except quite near the formation time, as is seen in Figure 2b.. Finally, energy conservation should be fine until the high frequency components play their role, as we see in Figure 2c. Hence, the results observed appear to be consistent with a singularity forming near $`t=8`$.

This does not guarantee that a singularity forms in these wave maps. There are other effects that might produce the apparently unbounded growth of the derivative of $`\chi `$ and of the energy density near the origin in these numerical simulations. For example, perhaps some unphysical, unstable mode grows because of the particulars of our chosen evolution scheme. We believe that this is not the case, for a number of reasons. First, the presence of such a mode would likely cause much larger growth in $`\mathrm{\Delta }(t)`$ than we see. Second, such modes would have to be excited only after some time (roughly independent of resolution) and only for families of sufficiently large energy. This is not consistent with our observations. Third, the excitation of this sort of instability would almost certainly depend critically on the precise finite difference scheme. Because Bizoń and his collaborators  observe similar behavior, using a different numerical evolution scheme, this does not appear to be the case. Thus we believe it very unlikely that the effects we are seeing are the result of a numerically unstable nonsingular mode.

Another situation in which one might numerically observe the formation of singularities that do not in fact evolve analytically from the corresponding data is if the continuum PDE solution is regularized by high frequency components which cannot be seen by the finite grid resolutions we use. The rather strong convergence behavior we see in our numerical solutions leads us to believe that this is not happening. We note in particular that such unresolved components would have to be separated in frequency space from the nontrivial low frequency components by a substantial margin, with a large range of dynamically irrelevant frequencies separating the two regimes. This seems to be very unlikely.

## IV Conclusion

The numerical studies we present here very strongly support the contention, previously conjectured by Shatah and Struwe , that smooth initial data for wave maps from $`2+1`$ Minkowski spacetime into the round two sphere can develop singularities (with unbounded derivatives) in finite time. As we note, there are many ways in which the numerical exploration of possible singularity formation might produce misleading indications. However, we believe that as a consequence of the numerical tests we have carried out, together with those done independently by Bizoń and his collaborators , the formation of singularities is the most likely conclusion.

There is much more one would like to know about these spatially equivariant wave maps, as well as about those without such symmetry. One would like to know, for example, if the solutions assume any universal form as one approaches the singularity. Our work (see Fig. 3) supports the results of Bizoń et al  which indicate that indeed the family of static spherically equivariant wave maps

$$\chi (r)=\pm 2\mathrm{arctan}\left(\lambda r\right)$$
(11)
does serve as a sort of universal model for singularity formation. This needs to be studied further. One would also very much like to understand the behavior of the wave maps which evolve from initial data near the transition from subcritical to supercritical data. The recent numerical work of Bizoń et al  suggests that the static solutions (11) play a central role in the evolution of the transitional wave maps as well in that of supercritical ones; however, this issue needs further investigation. (Note the absence of any self-similar solutions to the 2+1 wave map equations; for 3+1 wave maps, such solutions play a key role in behavior of solutions evolving from critical or near critical data).

## Acknowledgments

We thank Piotr Bizoń for helpful discussions. Partial support for this work has come from NSF Grant PHY-9800732 at the University of Oregon. SLL is appreciative of the support of NSF PHY-9900644 and of the financial support of Southampton College.


text:
# Spin susceptibility of neutron matter at zero temperature

## Acknowledgements

We are indebted to A. Fabrocini for providing us the results of the second order perturbative corrections to the energy per particle of neutron matter for the AU6’ interaction and to C.J. Pethick for insightful discussions. Portions of this work were supported by MURST-National Research Projects and CINECA computing center.


text:
# Multiplicity distribution and spectra of negatively charged hadrons in Au+Au collisions at √𝑠__NN = 130 GeV

## Abstract

The minimum bias multiplicity distribution and the transverse momentum and pseudorapidity distributions for central collisions have been measured for negative hadrons ($`h^{}`$) in Au+Au interactions at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 130 GeV. The multiplicity density at midrapidity for the 5% most central interactions is $`dN_h^{}/d\eta |_{\eta =0}=280\pm 1(\text{stat})\pm 20(\text{syst})`$, an increase per participant of 38% relative to $`p\overline{p}`$ collisions at the same energy. The mean transverse momentum is $`0.508\pm 0.012`$ GeV/c and is larger than in central Pb+Pb collisions at lower energies. The scaling of the $`h^{}`$ yield per participant is a strong function of $`p_{}`$. The pseudorapidity distribution is almost constant within $`|\eta |<1`$.

The collision of high energy heavy ions is a promising laboratory for the study of nuclear matter at high energy density and the possible creation and decay of the Quark Gluon Plasma . A new era in the study of high energy nuclear collisions began in the year 2000 with the first operation of the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory.

The multiplicity and inclusive single particle transverse momentum ($`p_{}`$) distributions of hadrons are important tools for understanding the evolutionary path of the system created in the collision and help to determine the characteristics of the early, hot and dense phase. In this Letter we present the minimum-bias multiplicity distribution and the pseudorapidity ($`\eta `$) and $`p_{}`$ distributions for the 5% most central collisions for negative hadrons ($`h^{}`$) in Au+Au interactions at a center-of-mass energy of $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 130 GeV per nucleon pair, measured by the STAR detector at RHIC. The results are compared with reference data from $`p\overline{p}`$ collisions at a similar energy and collisions of heavy nuclei at a lower energy.

The main tracking detector for STAR is a large Time Projection Chamber (TPC), which measures charged particles in the pseudorapidity range $`|\eta |<1.8`$ with complete azimuthal acceptance. It is placed inside a uniform solenoidal magnetic field of strength 0.25 T. The trigger detectors are an array of scintillator slats (CTB) arranged in a barrel surrounding the TPC, and two hadronic calorimeters (ZDCs) at $`\pm 18`$ m from the detector center and at zero degrees relative to the beam axis. The ZDCs intercept spectator neutrons from the collision and provide a measure of the collision centrality. Further details on the apparatus can be found in .

During the Summer 2000 run, RHIC delivered collisions between Au nuclei at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 130 GeV. The data presented here are from a minimum-bias sample, triggered by a coincidence of signals above threshold in both ZDCs with the RHIC beam crossing. The ZDC threshold was set to ensure efficient detection of single spectator neutrons. The efficiency of the ZDC coincidence trigger for central events was measured using a high-threshold CTB trigger. The trigger efficiency was found to be above 99% for the entire range of multiplicities reported in this Letter.

The offline reconstruction found a primary vertex for each event by propagating the measured tracks through the field towards the beamline. The vertex resolution for high multiplicity events is approximately 150 $`\mu `$m, both perpendicular and parallel to the beam axis. Events used in the analysis have a vertex within $`\pm 95`$ cm of the center of the TPC along the beam axis. The vertex finding efficiency is 100% for events with more than 50 primary tracks in the TPC acceptance, decreasing to 60% for those with fewer than 5 primary tracks. 60,000 minimum-bias Au+Au events were used for this analysis.

Particle production was studied through the yield of primary negative hadrons, comprising mostly $`\pi ^{}`$ with an admixture of $`K^{}`$ and $`\overline{p}`$. The $`h^{}`$ distribution includes the products of strong and electromagnetic decays. Negatively charged hadrons were studied in order to exclude effects due to participant nucleons. Charged particle tracks reconstructed in the TPC were accepted for this analysis if they fulfilled requirements on number of points on the track and on pointing accuracy to the event vertex. The measured raw distributions were corrected for acceptance, track reconstruction efficiency, and contamination due to interactions in material, misidentified non-hadrons, and the products of weak decays. The reconstruction efficiency was determined by embedding simulated tracks into real events at the raw data level, reconstructing the full events, and comparing the simulated input to the reconstructed output. This technique requires a precise simulation of isolated single tracks, achieved by a detailed simulation of the STAR apparatus based on GEANT  and a microscopic simulation of the TPC response. The multiplicity of the embedded tracks was limited to 5% of the multiplicity of the real event in the same phase space as the simulated data, thereby perturbing the real event at a level below the statistical fluctuations within the event sample.

The acceptance is on average 95% for tracks within the fiducial volume having $`p_{}>300`$ $`\mathrm{MeV}/c`$. The tracking efficiency ranges between 70%–95%, depending on $`p_{}`$ and the multiplicity of the event. For tracks with $`p_{}>200`$ $`\mathrm{MeV}/c`$ the efficiency is above 85%. Accepted tracks for this analysis have $`0.1<p_{}<2`$ $`\mathrm{GeV}/c`$ and $`|\eta |<1.0`$.

Instrumental backgrounds due to photon conversions and secondary interactions with detector material were estimated using the detector response simulations mentioned above, together with events generated by the HIJING model . The simulations were calibrated using data in regions where background processes could be directly identified. The measured yield also contains contributions from the products of weak decays, primarily $`K_S^0`$, that were incorrectly reconstructed as primary tracks. The background fraction of the raw signal is approximately 20% at $`p_{}`$ = 100 $`\mathrm{MeV}/c`$, decreasing with increasing $`p_{}`$. The average fraction of background tracks in the uncorrected sample is 7%. All corrections were calculated as a function of the uncorrected event multiplicity. The systematic uncertainty due to the corrections was estimated by studying the variation in the final spectra due to both a large variation in the track quality cuts with corresponding recalculation of the correction factors, and a small variation in the track quality cuts with no adjustment of the correction factors.

Figure 1 shows the corrected, normalized multiplicity distribution within $`|\eta |<0.5`$ and $`p_{}>100\mathrm{MeV}/c`$ for minimum bias Au+Au collisions. The data were normalized assuming a total hadronic inelastic cross section of 7.2 barn for Au+Au collisions at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 130 GeV, derived from Glauber model calculations . The multiplicity bin below $`N_h^{}`$ = 5 is not shown, due to large systematic uncertainties in the vertex reconstruction efficiency and large background contamination. Its relative contribution to the total cross section was estimated to be 21% by normalizing the HIJING multiplicity distribution to the measured data in the region $`5<N_h^{}<25`$. This procedure relies on the assumption that very peripheral interactions are well described by the superposition of a few nucleon-nucleon collisions in the geometry of a nuclear collision, and can therefore be accurately modelled by HIJING. The systematic error on the vertical scale is estimated to be 10% and is dominated by uncertainties in the total hadronic cross section and the relative contribution of the first bin. The systematic error on the horizontal scale is 6% for the entire range of multiplicity and is depicted by horizontal error bars on a few data points only.

The shape of the $`h^{}`$ multiplicity distribution is dominated over much of the $`N_h^{}`$ range by the nucleus-nucleus collision geometry, consistent with findings at lower energies. However, the shape of the tail region at large $`N_h^{}`$ is determined by fluctuations and acceptance. These overall features are also observed in the HIJING calculation, shown as histogram in Fig. 1. The distribution for the 5% most central collisions (360 mb), defined using ZDC coincidence, is shown as the shaded area in Fig. 1.

Figure 2, upper panel, shows the transverse momentum distribution of negatively charged hadrons for central Au+Au collisions at midrapidity ($`|\eta |<0.1`$) within $`0.1<p_{}<2\mathrm{GeV}/c`$. Statistical errors are smaller than the symbols. The correlated systematic error is estimated to be below 6%. The data are fit in the range $`0.2<p_{}<2`$ $`\mathrm{GeV}/c`$ by a QCD inspired power-law function of the form $`d^2N_h^{}/dp_{}^2d\eta =A(1+p_{}/p_0)^n`$ where $`A`$, $`n`$, and $`p_0`$ are free parameters. The upper panel of Fig. 2 also shows the $`p_{}`$-distributions of negatively charged hadrons for central Pb+Pb collisions at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 17 GeV from NA49  and for minimum-bias $`p\overline{p}`$ collisions at $`\sqrt{s}`$ = 200 GeV from UA1 , fitted with the same function. The NA49 distribution, which was reported in units of pion rapidity, was transformed to units of pseudorapidity. The UA1 invariant cross section $`Ed^3\sigma /d^3p`$ reported in Ref.  was scaled by $`2\pi /\sigma _{\mathrm{inel}}`$, where $`\sigma _{\mathrm{inel}}=42`$ mb . The power law fits all three datasets well. The mean $`p_{}`$ can be derived from the fit parameters as $`p_{}=2p_0/(n3)`$. The fit to the STAR data gives $`p_0=3.0\pm 0.3\mathrm{GeV}/c`$, $`n=14.8\pm 1.2`$, and $`p_{}=0.508\pm 0.012\mathrm{GeV}/c`$. The strong correlation of fit parameters $`p_0`$ and $`n`$ must be taken into account when calculating the error on $`p_{}`$. The $`p_{}`$ from STAR is larger than that from both central collisions of heavy nuclei at much lower energy ($`p_{}_{\mathrm{NA49}}0.429\mathrm{GeV}/c`$ ) and nucleon-nucleon collisions at a comparable energy ($`p_{}_{\mathrm{UA1}}=0.392\pm 0.003\mathrm{GeV}/c`$).

The PHENIX collaboration has reported that $`E_{}`$/$`N_{ch}`$ for central collisions of heavy nuclei is constant between SPS and RHIC energies, with a value of $`0.8`$ GeV . Under the assumption that the particle composition does not change significantly between the SPS and RHIC, this finding is in apparent disagreement with the above observation that $`p_{}`$ increases by 18%. It should be noted, however, that the systematic error on the SPS measurement of $`E_{}`$/$`N_{ch}`$ is 20% , so that these results are consistent within errors.

Figure 2, lower panel, shows the ratio of the STAR and UA1 $`p_{}`$-distributions. Since the UA1 distribution is measured at $`\sqrt{s}`$ = 200 GeV, $`d\sigma /dp_{}`$ is scaled by two factors for quantitative comparison to the STAR data at 130 GeV: (i) $`R(130/200)`$, the $`p_{}`$-dependent ratio of invariant cross sections for charged particle production in $`p\overline{p}`$ collisions at $`\sqrt{s}`$ = 130 and 200 GeV, and (ii)  $`\mathrm{T}_{\mathrm{AA}}`$= $`26\pm 2`$ mb<sup>-1</sup>, the nuclear overlap integral  for the 5% most central Au+Au collisions. $`R`$ varies from 0.92 at $`p_{}=0.2\mathrm{GeV}/c`$ to 0.70 at $`p_{}=2.0\mathrm{GeV}/c`$, and was derived using scaling laws for $`p_{}`$ and $`dN_{\mathrm{ch}}/d\eta `$ as a function of $`\sqrt{s}`$  together with the extrapolation to 130 GeV of power-law parameterizations at $`\sqrt{s}`$ = 200–900 GeV . The shaded boxes show the total error of the ratio, which is the linear sum of the errors of the measured data, depicted by the error bars, and the systematic error due to uncertainties in the scaling with $`\mathrm{T}_{\mathrm{AA}}`$ and $`R`$.

There are two simple predictions for the scaled ratio. In lower energy hadronic and nuclear collisions, the total pion yield due to soft (low $`p_{}`$) processes scales as the number of participants (“wounded” nucleons) in the collision (e.g. ). The scaled ratio in this case is 0.164, assuming 172 participant pairs  and a mean number of binary collisions of 1050 ($`=\sigma _{\mathrm{inel}}\mathrm{T}_{\mathrm{AA}}`$, where for $`\sqrt{s}=130\mathrm{GeV}\sigma _{\mathrm{inel}}40.5`$ mb ) for the 5% most central Au+Au events. In contrast, if hadron production is due to hard (high $`p_{}`$) processes and there are no nuclear-specific effects (see below), the hadron yield will scale as the number of binary nucleon-nucleon interactions in the nuclear collision and the value of the ratio is unity. There are important nuclear effects which may alter the scaling as a function of $`p_{}`$ from these simple predictions, including initial state multiple scattering , shadowing , jet quenching , and radial flow . The scaled ratio exhibits a strong $`p_{}`$ dependence, rising monotonically with increasing $`p_{}`$ from Wounded Nucleon scaling at low $`p_{}`$ but not reaching Binary Collision scaling at the highest $`p_{}`$ reported. This is consistent with the presence of radial flow, as well as the onset of hard scattering contributions and initial state multiple scattering with rising $`p_{}`$.

Figure 3 shows the normalized pseudorapidity distribution of $`h^{}`$ for the 5% most central collisions within $`|\eta |<1.0`$, both for $`p_{}>100\mathrm{MeV}/c`$ and for all $`p_{}`$. The latter was obtained by fitting a power-law function in the range $`0.1<p_{}<2\mathrm{GeV}/c`$ and extrapolating to $`p_{}=0`$ in order to estimate the content of the first bin. The error bars indicate the uncorrelated systematic errors. The statistical errors are negligible. The correlated systematic error applied to the overall normalization is estimated to be below 6% for $`p_{}>100`$ $`\mathrm{MeV}/c`$ and 7% for all $`p_{}`$.

The $`\eta `$ distribution is almost constant within $`|\eta |<1`$, exhibiting a small rise at larger $`\eta `$. This shape is expected from a boost invariant source (i.e., constant in rapidity), taking into account the transformation from $`y`$ to $`\eta `$. It should be noted that in Pb+Pb collisions at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 17 GeV the pseudo-rapidity distribution of charged hadrons  and the rapidity distribution of negative hadrons (assuming the pion mass)  were found to peak at midrapidity, suggesting a significant change in the longitudinal phase space distribution between the SPS and RHIC. Measurement of the rapidity distribution of identified particles is needed to establish the boost invariance at RHIC.

The $`h^{}`$ density at midrapidity for $`p_{}>100\mathrm{MeV}/c`$ is $`dN/d\eta |_{\eta =0}=261\pm 1(\text{stat})\pm 17(\text{syst})`$. Extrapolation to $`p_{}`$ = 0 yields $`dN/d\eta |_{\eta =0}=280\pm 1(\text{stat})\pm 20(\text{syst})`$. Assuming an average of 172 participant pairs per central Au+Au collision, this corresponds to $`1.63\pm 0.12`$$`h^{}`$ per participant nucleon pair per unit pseudorapidity, a 38% increase over the yield in $`p\overline{p}`$ collisions extrapolated to the same energy  (we neglect isospin correction factors of order 1–3%) and a 52% increase over Pb+Pb collisions at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 17 GeV .

For the total charged multiplicity density in Au+Au interactions at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 130 GeV, the PHOBOS collaboration has reported $`dN_{ch}/d\eta |_{|\eta |<1}=555\pm 12(\mathrm{stat})\pm 35(\mathrm{syst})`$ (6% most central collisions) , while the PHENIX collaboration has reported $`dN_{ch}/d\eta |_{\eta =0}=622\pm 1(\mathrm{stat})\pm 41(\mathrm{syst})`$ (5% most central) . To compare to these results, positive charged particles were analysed in the framework described above. For the 5% most central collisions, STAR measures a total charged multiplicity density of $`dN_{ch}/d\eta |_{\eta =0}=567\pm 1(\mathrm{stat})\pm 38(\mathrm{syst})`$.

In conclusion, we find that particle production per participant in central Au+Au collisions at $`\sqrt{s_{_{\mathrm{NN}}}}`$ =130 GeV increases by 38% relative to $`p\overline{p}`$ and 52% compared to nuclear collisions at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 17 GeV. The $`p_{}`$ distribution is harder than that of the $`p\overline{p}`$ reference system for the $`p_{}`$ region up to 2 $`\mathrm{GeV}/c`$. Scaling of produced particle yield with number of participants shows a strong dependence on $`p_{}`$, with Wounded Nucleon scaling achieved only at the lowest measured $`p_{}`$. The $`h^{}`$ pseudorapidity distribution is almost constant within $`|\eta |<1`$.

Acknowledgments: We wish to thank the RHIC Operations Group at Brookhaven National Laboratory for their tremendous support. This work was supported by the Division of Nuclear Physics and the Division of High Energy Physics of the Office of Science of the U.S.Department of Energy, the United States National Science Foundation, the Bundesministerium fuer Bildung und Forschung of Germany, the Institut National de la Physique Nucleaire et de la Physique des Particules of France, the United Kingdom Engineering and Physical Sciences Research Council, Fundacao de Amparo a Pesquisa do Estado de Sao Paulo, Brazil, and the Russian Ministry of Science and Technology.


text:
# Is the particle current a relevant feature in driven lattice gases?

## Abstract

By performing extensive MonteCarlo simulations we show that the infinitely fast driven lattice gas (IDLG) shares its critical properties with the randomly driven lattice gas (RDLG). All the measured exponents, scaling functions and amplitudes are the same in both cases. This strongly supports the idea that the main relevant non-equilibrium effect in driven lattice gases is the anisotropy (present in both IDLG and RDLG) and not the particle current (present only in the IDLG). This result, at odds with the predictions from the standard theory for the IDLG, supports a recently proposed alternative theory. The case of finite driving fields is also briefly discussed.

PACS numbers: 64.60.-i, 05.70.Fh

The Ising model exhibits a prototypical equilibrium phase transition, and the associated $`\varphi ^4`$ Ginzburg-Landau theory is a paradigm of continuous theory for equilibrium critical phenomena . However, thermodynamic equilibrium is exceptional in nature where stationary states are typically away from equilibrium . With the purpose of defining simple lattice models describing generic non-equilibrium phase transitions, different tentatives have been made in the last two decades. Among them, perhaps the most intriguing example is the driven lattice gas (DLG) . (Other interesting examples are the directed percolation model  and the Kardar-Parisi-Zhang equation). The DLG, being a straightforward extension of the Ising model, has in fact become a workbench for emergent non-equilibrium theories and field theoretical approaches.

The DLG is a d-dimensional kinetic Ising model with conserved dynamics, in which transitions in the direction (against the direction) of an externally applied field, $`\stackrel{}{E}`$, are favored (unfavored) , while transitions perpendicular to the field are unaffected by it. The field induces two main non-equilibrium effects: (i) the presence of a net current of particles along its direction, and (ii) strongly anisotropic configurations . At high temperatures, the system is in a disordered phase while, for half-filled lattices (the only case we refer to in what follows) there is a second-order critical point, below which the DLG segregates into (two) high and low density aligned-with-the-field stripes. Establishing unambiguously the DLG universality class is an important issue in the way to rationalize the behavior of non-equilibrium systems.

Continuous approaches such as Langevin and associated field theories  have been most useful in studying universality issues in equilibrium critical phenomena. In particular, coarse-grained approaches combined with renormalization group (RG) techniques provide a method for the classification of the different possible terms (operators) as relevant, irrelevant or marginal. In fact, Langevin equations are more illuminating than other (even more rigorous) approaches, as they permit to understand systematically how possible perturbations or model variations would affect critical properties. Consequently, many studies have focused on the DLG and its universality by using both non-equilibrium continuous approaches and computer simulations (unfortunately, general exact solutions are not available). Within this perspective, it is somewhat deceptive that after many computer and analytical studies, the universality class of the DLG remains a debated issue .

A phenomenological Langevin equation intended to capture the relevant physics of the DLG at criticality was proposed and renormalized more than a decade ago . This equation, referred to as driven diffusive system (DDS), is a natural extension of the conserved $`\varphi ^4`$ theory for the Ising equilibrium transition (model B ) and seems to capture the main symmetries and conservation laws of the discrete DLG. It includes a particle current term (which from naive power counting turns out to be the most relevant nonlinearity) as well as anisotropic coefficients. It certainly is a suitable and very reasonable candidate to be the canonical coarse-grained model representative of the DLG universality class. The DDS Langevin equation reads:

$`_t\varphi (𝐫,t)=`$ $`\tau _{}_{}^2\varphi _{}^4\varphi +{\displaystyle \frac{\lambda }{6}}_{}^2\varphi ^3`$ (2)
$`+\tau _{}_{}^2\varphi \alpha _{}\varphi ^2+\zeta (𝐫,t),`$
where $`\varphi `$ is the coarse grained field, $`\zeta `$ is a conserved Gaussian noise and the cubic term, being a dangerously irrelevant variable , is kept in order to ensure stability . $`\tau _{}`$, $`\tau _{}`$, $`\lambda `$ and $`\alpha `$ are model parameters. The most emblematic (exact) prediction derived from the DDS RG-analysis, namely, the mean field behavior of the order parameter critical exponent, $`\beta =1/2`$,  has eluded a large number of Monte Carlo (MC) analysis aimed at probing it , however. In particular, systematic deviations from the predicted scaling are observed both in $`d=2`$  and in $`d=3`$ . Indeed, different MC analysis (performed using a variety of aspect-ratios and order parameters) lead systematically to a value of $`\beta `$ close to $`0.3`$ (in $`d=2`$), with error bars excluding apparently the value $`\beta =1/2`$ (see  for a critical review of simulation analysis). This is a main indication that, strikingly enough the DDS equation does not describe properly the infinitely fast driven DLG (IDLG) critical properties. Moreover, there are some other hints suggesting strongly that the differences between the predictions of the standard Langevin approach and MC results are more fundamental than a simple discrepancy in the value of $`\beta `$. In particular, the intuition developed from MC simulations of the DLG and variants of it  suggests that, contrarily to what the DDS equation establishes, it is the anisotropy and not the presence of a current the basic ingredient controlling the critical behavior . For instance, in a modified DLG in which anisotropy is included by means others than a current , the scaling behavior at criticality remains unaltered upon switching on an (infinite) driving (see ). Other compelling evidences supporting this hypothesis can be found in .

In an attempt to clarify this puzzling situation, and reconcile continuous approaches with numerics, different scenarios have been explored. In particular, an alternative route to build up Langevin equations starting from generic microscopic master equations was recently proposed . By applying this approach to the DLG, one observes that, owing to a transition-rates saturation effect, the coefficient $`\alpha `$ of the non-linear current term, $`_{}\varphi ^2`$, vanishes in the limit of infinite driving fields and, therefore, it does not appear in the final Langevin equation nor it is generated perturbatively . The resulting theory (alternative to Eq.(2)) is:

$$_t\varphi (𝐫,t)=\tau _{}_{}^2\varphi _{}^4\varphi +\frac{\lambda }{6}_{}^2\varphi ^3+\tau _{}_{}^2\varphi +\zeta $$
(3)
plus higher order irrelevant contributions (note that a linear current term has been eliminated by employing a Galilean transformation ). This equation, named below anisotropic diffusive system (ADS), is a well known one: it coincides with the Langevin equation representing the random DLG (RDLG)  (for which the driving field takes values $`\mathrm{}`$ and $`\mathrm{}`$ in a random unbiased fashion, generating anisotropy but not an overall current). This theory has been extensively studied in ; its critical dimension is $`d_c=3`$ (instead $`d_c=5`$ for the DDS) and the critical exponents and finite size scaling (FSS) properties are now well known. Other systems in this universality class are the two-temperature model  and the ALGA model . This theory for the IDLG includes anisotropy as its basic non-equilibrium ingredient. Instead - for non-saturating, finite, driving fields - the cancellation of the nonlinear current term does not occur, and our method recovers the standard DDS equation.

Aiming at further clarifying these issues, we report here on extensive MC simulations of the IDLG and the RDLG in $`d=2`$. The main objectives are: (i) trying to conclude whether the IDLG and the RDLG share the same critical behavior or not; and (ii) measuring the critical exponents by performing systematic anisotropic finite-size scaling (FSS). In fact, we perform FSS analysis for both the IDLG and the RDLG by following the anisotropic FSS scheme proposed in  consistent with the ADS theory; this allows us to analyze systematically possible scaling differences between both models. We also report on the case of finite-driving DLG.

We consider rectangular lattices of size $`L_{}\times L_{}`$ with periodic boundary conditions and random sequential updating ; the external field $`\stackrel{}{E}`$ acts in the $`x`$ (parallel) direction. Particles jump to a randomly chosen nearest neighbor site (provided that it is empty) with probability: $`\mathrm{min}(1,\mathrm{exp}[\beta (\mathrm{\Delta }H+E\mathrm{\Delta }j)])`$, where $`\mathrm{\Delta }H`$ is the energy (Ising Hamiltonian) variation, and $`\mathrm{\Delta }j=(1,0,1)`$ for jumps along, against, and orthogonal to the direction of the field, respectively. Following  the order parameter is chosen as the structure factor $`S(0,2\pi /L_{})`$. In order to perform a systematic anisotropic FSS we considered system sizes $`20\times 20`$, $`45\times 30`$, $`80\times 40`$, and $`125\times 50`$. These aspect ratios satisfy $`L_{}^{\nu _{}/\nu _{}}=0.2236\times L_{}`$, where $`\nu _{}/\nu _{}1/2`$ consistent with ADS anisotropic spatial scaling . The number of MC steps considered varied between $`1.8\times 10^8`$ and $`2.4\times 10^8`$, much larger than in any previously reported MC simulations. The total CPU time employed is about eight months in a pentiumIII 400Mhz machine. The critical temperature is determined by using the fourth (Binder) cumulant method . For the IDLG, the critical temperature is found to be $`T_c^I=1.396(4)T_o`$ ($`T_o`$ is the Onsager temperature), slightly below previously reported values , while we find $`T_c^R=1.390(4)T_o`$ for the RDLG (see insets of Fig.3). These critical values were employed for the FSS analysis. In Fig.1 we plot the order parameter, rescaled by a factor $`L_{}^{\beta /\nu _{}}`$, versus $`ϵL_{}^{1/\nu _{}}`$, where $`ϵ`$ is the distance to the critical point, for different system sizes $`L_{}`$. A nearly perfect data-collapse is obtained by fixing $`\nu _{}=1.25`$ and $`\beta =0.33`$. The collapse gets worse upon slightly changing these values; more precise estimates of the associated error-bars is a difficult and not essential issue in this context. Notice that we are plotting in the same graph data for the IDLG and for the RDLG, implying that the FSS scaling function is precisely the same for both models. Furthermore, the slopes of the asymptotic branches are approximately $`1/3`$ and $`0.61`$, consistent both with the order parameter exponent being $`\beta 0.33`$, and $`\gamma 1.22`$ (see below) . In general, even when the scaling functions are universal their corresponding amplitudes are not expected to be so. For this reason, usually one has to introduce the, so called, metric factors (varying amplitudes)  in order to obtain superposition of scaling functions within the same universality class. Contrarily to this expectation, the magnetization scaling functions for IDLG and RDLG overlap perfectly. Therefore, it comes as a surprise that not only the scaling functions and the $`\beta `$ exponent are universal in both models, but even the amplitudes coincide. A similar situation has been recently reported for a different type of anisotropic FSS .

We have also computed the system susceptibilities, defined as the relative fluctuations of the order parameter: $`\chi =\frac{L_{}}{sin(\pi /L_{})}[<m^2><m>^2]`$. In Fig.2, we plot the susceptibility times $`L_{}^{\gamma /\nu _{}}`$ as a function of the rescaled distance to the critical point, $`ϵL_{}^{1/\nu _{}}`$. The best data collapse is obtained by employing the values $`\gamma =1.22`$ and $`\nu _{}=1.25`$ for both models with, again, coinciding amplitudes. It should be stressed that this is the first time a really good collapse is observed below the critical point for anisotropic scaling of the IDLG. Plotting the dimensionless Binder cumulant as a function of the rescaled distance to the critical point with $`\nu _{}=1.25`$, again, nearly perfect data collapse is obtained for both models and all system sizes (Fig. 3). We also performed simulations in square lattices ($`128\times 128`$) as in some previous studies . Monitoring $`m^{1/\beta }`$ as a function of $`T/T_o`$ we see no appreciable systematic difference between the curves for IDLG and RDLG, that have the same slope within numerical accuracy. The best linear fit correlation is obtained for $`\beta 0.33`$ in both cases, providing an extra consistency check for our results. Moreover, eye inspection of IDLG and RDLG configurations, for any geometry, at a fixed relative temperature, does not permit to distinguish one from the other. In particular, the interfacial properties look identical. Let us also stress that all the obtained exponent values are compatible with previous measures for the RDLG, as well as with the exponents obtained within an $`ϵ`$expansion of the ADS theory .

In conclusion, MC results support strongly that both the IDLG and the RDLG belong in the same universality class, and share not only critical exponents and scaling functions, but also the scaling amplitudes. This universality class is described by the ADS equation, Eq.(3). There is absolutely no hint of any difference in the asymptotic behavior between the model with a current (IDLG) with respect to the current-less one (RDLG). All the numerical evidence confirms that it is the anisotropy and not the net current (if any) the most relevant non-equilibrium ingredient of driven systems. As discussed in the introduction, this is striking from a field theoretical perspective given that the nonlinear current term, $`_{}\varphi ^2`$, is naively a relevant perturbation at the ADS fixed point. In an alternative approach, the coefficient of such a term vanishes. In this picture, the fast drive limit corresponds to a sort of multicritical point in which an a priori relevant operator is absent due to a cancellation of its coefficient and, consequently the usual ”up-down” Ising symmetry (i.e. the three-point correlation functions vanish) is restored at criticality. In any case, it should be stressed that, from a more general perspective, field theoretical descriptions of non-equilibrium systems are much more delicate and subtle than their equilibrium counterparts, and an extremely careful inspection of the system symmetries, conservation laws, and dynamical features is required before venturing to make predictions on nonequilibrium universality issues. For example, fermionic and bosonic non-equilibrium systems with the same dynamics, symmetries and conservation laws have recently reported to belong to different universality classes . One could wonder whether the DLG hard-core interaction should be taken into account into its Langevin representation.

Elucidating the critical behavior for finite $`E`$ remains a challenging and interesting objective. Both, our alternative Langevin-building approach and the standard one, Eq.(2), include a relevant current term in this case and, consequently, predict $`\beta =1/2`$. Obtaining clear-cut results in this case is a computationally expensive task since: (i) As the external field appears in the argument of an exponential, even relatively small values of $`E`$ generate situations close to saturation, and strong crossover effects could hide the true asymptotic regime. (ii) If the field value is taken too small, crossovers from the equilibrium regime may also burden observations. A possible scenario that could follow from MC analysis is that finite fields show mean field behavior; that would be a strong backing for our theory  that predicts the finite and the infinite driving cases to be qualitatively different. If, instead, scaling happens to be that of the ADS (as our preliminary MC results for $`E=3`$ and $`E=0.5`$ seem to indicate; for $`E=0.25`$ results do not quite fit this indication) it would prove that it is for any arbitrary value of the driving field that anisotropy is the most relevant ingredient of driven systems. This scenario would uncover a new puzzling situation and would certainly call for deeper theoretical understanding. Huge and careful simulations would be required to extract neat conclusions overcoming difficulties (i) and (ii) above.

Summing up, we have performed extensive MC simulations of the IDLG and the RDLG. By using anisotropic finite size scaling techniques we have shown that both models belong to the same universality class: their critical exponents, scaling functions and amplitudes are undistinguishable and coincide with those of the ADS equation. This result supports the conclusion that it is the presence of anisotropic coefficients, and not the particle current the most relevant ingredient in these non-equilibrium driven problems (at least in the fast drive limit). Further theoretical efforts are certainly required in order to (i) to sort out if our alternative Langevin approach is correct and what are its possible limitations, and (ii) to further clarify the universality issues of this quintessential non-equilibrium problem. Finally, it would also be very interesting to combine the powerful finite size methods recently introduced in this context by Caracciolo et al. in a nice recent work  with our alternative theory to verify if they lead to better data collapse than when used to test the standard DDS equation (hopefully without having to introduce strong corrections to scaling and providing good order-parameter scaling).

ACKNOWLEDGMENTS- We are grateful to S. Caracciolo for his very valuable comments. We also thank F. de los Santos and Prof. A. Sekkaki. This work has been supported by the European Network Contract ERBFMRXCT980183 and by the MCYT under project PB97-0842.