text:
# Abstract

## Abstract

Recent numerical simulations of MHD jets from accretion disks are briefly reviewed with emphasis on the scaling law for jet speed and the role of magnetic reconnection in relation to time variability in accretion disks, jets, and flares. On the basis of these studies, possible interpretation is given on why statistical properties of peak intensity, peak interval, and peak duration of gamma ray bursts (log-normal distribution) are different from those in solar flares and black hole accretion disks (power-law distribution). From these considerations, a new model, “magnetized plasmoid model”, is proposed for a central engine of gamma ray bursts.

## 1 Introduction

Recent development of astronomical observations has revealed that our universe is full of enigmatic explosive phenomena, such as jets, bursts, and flares. Jets ejected from active galactic nuclei (AGN) are probably among the biggest in size and the most energetic in total energy. Similar jets on much smaller scale have been found in young stellar objects (YSO) as well as in close binary system. These active objects show vigorous time variability, often called bursts or flares, in optical, radio waves, X-rays, etc., in almost all electromagnetic spectrum. Though the central objects and jets in AGNs, YSOs, and binary system are quite different in mass and size, there are many similar properties in them, such as time variability, morphology of jets, and existence of accretion disks.

Gamma ray bursts (GRBs) were discovered nearly 30 years ago. Since then, they have remained the most enigmatic object in our universe. Recent rapid development of observations, however, is uncovering a part of enigma of GRBs (e.g., Fishman and Meegan 1995, Piran 1999, Meszaros 2002); (1) they occur in cosmological distance (i.e., most luminous in our universe, $`10^{5153}`$ erg/s, but different from AGNs), (2) there is evidence that they are emitted from relativistic jets via synchrotron emission, and (3) some GRBs (long bursts) seem to be associated with supernovae. Some properties are similar to those of AGN jets (blazars), suggesting common physics in blazars and GRBs.

On the other hand, recent space observations of our Sun, the nearest star, have revealed that the solar corona is full of jets and flares. Though the total energy of the solar jets and flares is much smaller than those in comic jets and flares, the spectrum and time variability of electromagnetic waves emitted from solar flares are quite similar to those of cosmic flares, suggesting common physical origin. In the case of the solar flares, it has been established that magnetic field is the source of energy, so that the knowledge of solar jets and flares will be useful for understanding the role of magnetic field in cosmic jets and flares in distant stars and galaxies.

In this article, we first briefly review recent understanding of magnetically driven jets from accretion disks (as a model of AGN jets). We then discuss flares and magnetic reconnection associated with production of jets, especially in the case of protostars. Finally, we discuss similarity and differences in time variability of GRBs, black hole accretion disks, and solar flares, and on the basis of these observations, we propose a new model for central engine of GRBs.

## 2 MHD Jets

Magnetically driven jets from accretion disks (Fig. 1) were first proposed to explain jets from active galactic nuclei (Blandford 1976, Lovelace 1976, Blandford and Payne 1982). After the discovery of CO molecular bipolar flows in star forming regions (Snell et al. 1980), the MHD jet model started to be applied to bipolar flows and jets from young stars (Uchida and Shibata 1985, Pudritz and Norman 1986). The first time-dependent numerical simulation of MHD jets from accretion disks were carried out by Uchida and Shibata (1985) and Shibata and Uchida (1986).

The MHD jet model (e.g., Tajima and Shibata 1997, Ferrari 1998, for a review) has the following merits: (1) the magnetic force not only accelerate plasmas from disk surface to form bipolar jets but also extract angular momentum from accretion disks, enabling efficient accretion of plasma onto central objects (stars or black holes), (2) the magnetic force due to toroidal fields collimate jets by pinching effect.

Shibata and Uchida (1986) and subsequent many 2.5D MHD numerical simulations (Stone and Norman 1994, Matsumoto et al. 1996, Hirose et al. 1997, Ouyed and Pudritz 1997, Kudoh et al. 1998, 2002, Kuwabara et al. 1999, Ustyugova et al. 1999, Kato et al. 2002) and 3D MHD simulations (Matsumoto and Shibata 1997, Ouyed et al. 2003) have revealed that the velocity of MHD jets is of order of Keplerian velocity at the footpoint of MHD jets. Kudoh et al. (1998) and Kato et al. (2002) have shown that the jet velocity has a weak dependence on $`B`$ (Fig. 2);

$$V_{jet}E_{mg}^{}{}_{}{}^{1/6}B_{p}^{}{}_{}{}^{1/3}.$$
$`(1)`$
Here,$`E_{mg}=(V_A/V_k)^2B_{p}^{}{}_{}{}^{2}`$, and $`B_p`$ is the initial poloidal magnetic field strength.

This is consistent with the results of one-dimensional steady jet theory (Kudoh and Shibata 1995, 1997), and also can be derived semi-analytically using the Michel (1969)’s relation;

$$V_{\mathrm{}}\left(\frac{\mathrm{\Omega }^2B_p^2r^4}{\dot{M}}\right)^{1/3}.$$
$`(2)`$
Here, $`\mathrm{\Omega }`$ is the angular speed of the disk at the footpoint of the jet, $`r`$ is the radial distance from the central object, and $`\dot{M}`$ is the mass flux of the jet.

It is very important to note that the terminal velocity strongly depends on the mass flux $`\dot{M}`$. In our problem, the mass flux is given by

$$\dot{M}4\pi \rho _{slow}V_{slow}r^24\pi \rho _{slow}C_s\frac{B_p}{B}r^2,$$
$`(3)`$
where $`\rho _{slow}`$ is the mass density at the slow magnetosonic point, $`V_{slow}`$ is the slow magnetosonic speed, $`C_s`$ is the sound speed, and $`B=(B_p^2+B_\phi ^2)^{1/2}`$. In a cold disk, the slow magnetosonic point corresponds to the local maximum of the effective gravitational (Blandford-Payne) potential, and is located near the disk plane ($`\rho _{slow}/\rho _00.1`$ in our case). For stronger fields ($`E_{mg}>10^2`$), the magnetic field lines become straight $`BB_pB_\phi `$ so that the mass flux does not depend on $`B_p`$, but for weaker fields ($`E_{mg}<10^2`$), field lines are highly twisted and the azimuthal component becomes dominant $`BB_\phi B_p`$ near the disk so that the mass flux is in proportion to $`B_p`$ (Kudoh and Shibata 1995, 1997). Considering these effects, the terminal speed of the MHD jet becomes

$$\frac{V_{\mathrm{}}}{V_k}\left(\frac{\rho _{slow}}{\rho _0}\right)^{1/6}E_{th}^{1/6}E_{mg}^{1/6},$$
$`(4)`$
where $`E_{th}=(C_s/V_k)^2`$ = (thermal energy)/ (gravitational energy), $`E_{mg}=(V_A/V_k)^2`$ = (magnetic energy)/ (gravitational energy). Since $`\rho _{slow}/\rho _00.1`$ and $`E_{th}/E_{mg}\beta =`$ gas pressure / magnetic pressure $`110`$ in the disk, we find that the terminal speed of the jet is comparable to the Keplerian velocity for wide range of poloidal magnetic field strength, and $`V_{\mathrm{}}B_p^{1/3}`$. This explains the results (eq. 1) of the previous 2.5D MHD numerical simulations very well.

We should emphasize again that even if the magnetic field strenth is very weak in accretion disks, the jet velocity is roughly comparable to Keplerian speed ($`V_{jet}0.11.0V_k`$ for $`E_{mg}10^810^2`$). The physical reason is that magnetic field lines are highly twisted by the differential rotation of the disk until the local magnetic energy density ($`B_\phi ^2/8\pi `$) becomes comparable to the rotational energy ($`\rho V_k^2/2\rho GM/r`$ gravitational energy) at the surface of the disk. Since the kinetic energy of the jet ($`\rho V_{jet}^2/2`$) comes from the magnetic energy, it eventually becomes comparable to the gravitational energy ($`\rho V_k^2/2`$) at the disk surface (i.e., at the slow magnetosonic point). This process is similar to the magneto-rotational instability (Balbus and Hawley 1991) in the sense that the magnetic effect becomes eventually important even if the initial magnetic field is very weak.

Recently, Koide et al. (1998, 2000, 2002) have succeeded to extend these newtonian MHD simulations of jets to general relativistic MHD versions. They have shown that the accretion (and the jet ejection) become more violent near black hole than in the newtonian case. They confirmed the extraction of rotational energy from a Kerr hole by the effect of magnetic field. The maximum Lorentz factor of jets in these simulations is still of order of 2, much smaller than those observed for AGN jets (Lorentz factor 10-100) and GRBs (Lorentz factor 100-1000). We do not know yet whether this result (small maximum Lorentz factor) is simply a result of numerical limitations, or a result of physics.

It is important to note that the jet ejection has never reached steady state, even if the velocity and mass flux of the jet is well explained by steady theory (Kudoh et al. 1998). Recently, Sato et al. (2003) have succeeded to run the MHD simulation of jets including accretion disk self-consistently for many orbital periods (up to 15-20 orbits), and revealed that the jet ejection is intermittent and often associated with transient accretion events. They found that accretion disk is fully turbulent (due to magneto-rotational instability), and full of reconnection event. Often the jets are ejected in association with such reconnection events. The time variability of mass accretion rate in a simulated accretion disk show power law in the power density spectrum as first noted by Kawaguchi et al. (2000).

## 3 Flares : Magnetic Reconnection

Recent space solar observations such as Yohkoh, SOHO, TRACE have revealed that solar corona is much more dynamic than had been thought, and is full of flares, microflares, nanoflares, jets, and various mass ejections. Among them, the largest mass ejectins are called coronal mass ejections (CMEs). It has also been revealed that that the reconnection plays essential role not only in large scale flares and CMEs, but also small scale flares and jets, leading to unified model (e.g., Shibata 1999).

Hayashi, Shibata, Matsumoto (1996) presented a magnetic reconnection model of protostellar flares, by performing 2.5D time dependent MHD nummerical simulations of interaction between an accretion disk and stellar dipole magnetic field. Figure 3 shows one of their simulation results. They assumed that an accretion disk is penetrated by stellar dipole field at t = 0, and examined the subsequent evolution of the interaction between a rotating disk and a stellar dipole field. The initial process occurring near the disk is basically the same as those in the nonsteady MHD jet model (e.g., Shibata and Uchida 1986). The magnetic field is twisted by the rotating disk, and the J x B force associated with the twist accelerates the plasma in the surface layer of the disk to form an MHD jet in bipolar directions. In this case, the magnetic twist is accumulated in a closed loop, increasing the magnetic pressure of the loop, which eventually leads to the ejection of the magnetic loop after about one orbit. After the ejection of the loop, a current sheet is created inside the loop, leading to fast reconnection there. This process is similar to that occurring in solar coronal mass ejections, and basic reconnection mechanism is the same as in solar flares (e.g., Tsuneta et al. 1992, Shibata 1999, Yokoyama and Shibata 2001). The reconnection releases huge amount of magnetic energy of order of $`10^{36}`$ erg (about $`10^4`$ times more energetic than solar flares) stored in a sheared loop with a size of $`L10^{11}`$ cm.

The temperature of super hot plasma created by reconnection amounts to

$$T10^8\left(\frac{B}{100\mathrm{G}}\right)^{6/7}\left(\frac{n_0}{10^9\mathrm{cm}^3}\right)^{1/7}\left(\frac{L}{10^{11.5}\mathrm{cm}}\right)^{2/7}\mathrm{K},$$
$`(5)`$
which is based on the balance between reconnection heating and conduction cooling (Yokoyama and Shibata 2001, Shibata and Yokoyama 2002). Here $`n_0`$ is the pre-flare coronal density. These results explain characteristics of protostellar flares observed by ASCA and ROSAT (e.g., Koyama et al. 1996, Shibata and Yokoyama 2002).

## 4 Gamma Ray Bursts

It is well known that the gamma ray burst light curves are very similar to those of solar flare gamma rays and hard X-rays. Hence, it was often considered that there may be some physical similarity between gamma ray bursts and solar flares. Recently, it has in fact been revealed that there is evidence that gamma ray bursts are emitted from collimated jets, which might be similar to relativistic jets from active galactic nuclei. If this is so, MHD jet models developed for AGN jets can be applied to gamma ray bursts. There is also evidence that some neutron stars have very strong magnetic fields up to $`10^{15}`$ G, called magnetars (Duncan and Thompson 1992, Kouveliotou et al. 1999). Hence, the GRB model including strong magnetic field has been proposed (e.g., Kluzniak and Ruderman 1998). On the other hand, recent theory of accretion disks showed that magnetic field is essential for generating viscosity through the occurrence of magneto-rotational instability (Balbus and Hawley 1991). Altogether, the role of magnetic field in gamma ray bursts is considered to be more and more important than had been thought.

Here we shall discuss analogy between solar flares/coronal mass ejections (CME) and gamma ray bursts, emphasizing the basic MHD physics of solar flares/CMEs.

As we discussed above, the time variability of gamma ray burst light curve is similar to those of solar flares. The power spectrum analysis of time variability of gamma ray bursts (Beloborodov et al. 2000) and solar X-ray emission (Ueno et al. 1997) show that both show power-law disribution with index $`\alpha 1.51.8`$. Interestingly, the X-ray light curve of Cygnus X-1 (black hole candidate) also show similar time variability with power law spectrum (Miyamoto et al. 1992, Negoro 1992).

However, there is a fundamental difference in statistical properties between gamma ray bursts and solar flares. That is, many physical quantities in solar flares, such as duration, peak intensity, peak interval, show power-law frequency distributions, while these quantities in gamma ray bursts do not show power-law distribution but show log-normal distribution (Li and Fenimore 1996; see Fig. 4). In other words, there is no characteristic time and intensity in solar flares, whereas there is characteristic time and intensity in gamma ray bursts. What does this mean ? In spite of apparent similarity between solar flares and gamma ray bursts, does this suggest that basic physics of gamma ray bursts is different from that of solar flares ?

In relation to this, Negoro and Mineshige (2002) recently found interesting fact on statistical properties of X-ray shots (sub-bursts) emitted from black hole accretion disks: If only large shots are picked up from the time variability of X-ray emission of accretion disk, the distribution of shot peak become log-normal. They argued that the X-ray shots are occurring at the surface of accresion disks, and hence may be similar to solar flares and do not have characteristic time scale (Ueno et al. 1997, Kawaguchi et al. 2000). On the other hand, according to internal shock model of gamma ray bursts (e.g., Meszaros 2002), the burst emissions are from internal shocks in the jet ejected from central engines. If the jets are ejected in association with X-ray shots, only large shots can eject enough mass to generate internal shock. This explains log-normal distribution in the duration and peak intensity of GRB.

Similar discussion may be applied to solar flares. In the case of solar flares, there is a tendency that larger flares produce larger mass ejections, called coronal mass ejections (CMEs). It is difficult for small flares to eject much mass from solar surface, because they do not have enough energy to escape from magnetically confined solar active regions. Aoki et al. (2003) examined this property by using actual data. Figure 5 shows number of CMEs associated with solar flares versus peak X-ray flux of those flares, indicating that log-normal distribution roughly holds. The number of CMEs versus CME speeds also show similar log-normal distribution.

Hence we propose a new model for gamma ray bursts (see Figure 6), which is basically in the same line of thought as that of Negoro and Mineshige (2002): magnetic reconnection (flares) occurs everywhere in the surface of accretion disks, whose spatial distribution is fractal, and time variability is power-law, both of which are quite similar to those of solar flares. Only ejecta from energetically and spatially large reconnection events (flares or shots) can escape from magnetosphere of accretion disks to form jets. In fact, nonsteady MHD simulations of magnetically driven jets from accretion diks (Kudoh et al. 2002, Kuwabara et al. 2002, Sato et al. 2003) show that the ejection of jets is highly time variable, intermittent, and base of jets is full of reconnection events; when large reconnection events occur, large plasmoid ejection occurs, like the protostellar flare model (Hayashi et al. 1996, Miller and Stone 1997, Goodson et al. 1999). Jets have non-uniform density distribution, consisting of intermittently ejected plasmoid (confined in magnetic island or helical field in 3D space) like solar coronal mass ejections. It is natural that these intermittently ejected magnetized plasmoid produce lots of internal shocks, which are site of gamma ray emission as modeled in internal shock model of gamma ray bursts.

There is another merit in this “magnetized plasmoid model” (or “flare/CME model”). Ejected plasmoid has a structure similar to spheromak. It is well known that spheromak is unstable to tilting instability (e.g., Hayashi and Sato 1984), so that ejected spheromak-like plasmoids would have various angles to the direction of the ejection as they propagate. Hence they collide with each other, and there is a high possibility that tangent magnetic field has an opposite component, thus leading to magnetic reconnection. This means that even if the initial energy conversion from magnetic energy to plasma energy (kinetic and internal) is not efficient, it is possible that all magnetic energy (Poynting flux) contained in magnetized plasmoids would eventually be converted to plasma energy through magnetic reconnection. It should also be noted that magnetic reconnection generate lots of MHD shocks: Petschek slow shocks are formed just sides of reconnection jets, and fast shocks are created when reconnection jets collide with ambient medium (e.g., Yokoyama and Shibata 2001). These situations, containing lots of shocks as well as X-type and O-type neutral points, are similar to fractal MHD turbulence, and very suitable for high energy particle acceleration. Such fractal structure or fractal reconnection (Shibata and Tanuma 2001, Tanuma et al. 2001) may be the origin of power-law time variability spectrum of solar flares, black hole accretion disks, and GRBs.

We are grateful to T. Kudoh, S. X. Kato, K. Sato, S. Mineshige, H. Negoro, T. Murakami, and T. Ishii for various help and useful discussions.

Discussion

J. Rhoads: How do you generalize the relation $`V_{jet}=v_{escape}=v_{Kepler}`$ to the highly relativistic case ? Here $`V_{jet}=v_{escape}=c`$ is obvious, but it is not clear if $`\mathrm{\Gamma }_{jet}=\mathrm{\Gamma }_{escape}`$. For the Newtonian case, I assume that $`V_{escape}`$ is measured in the region where the jet is launched. Trying to apply the same idea to the $`\mathrm{\Gamma }=100`$ flow in a GRB would imply the jet is launched barely outside the event horizon. This seems improbable. How do we explain this ?

K. Shibata: You are right. In the Newtonian case, our theory explains observed speeds of various jets very well, such as protostellar jets, SS433 jets, jets from cataclysmic variables, and so on. However, in the relativistic case, we have not yet succeeded to explain observed large Lorentz factor of AGN jets and gamma ray bursts.

M. Lyutikov: In case of a rotating BH the flow is generated on field lines that do not cross BH horizon, but thread the ergosphere.

K. Shibata: Koide et al. (2002) successfully simulated such case.

A. Beloborodov: Keplerian timescale is very short ( ms) to explain the main variability in GRBs (0.1-10 s). Should the accretion timescale appear as a characteristic time in the model variability spectra ?

K. Shibata: MHD numerical simulations of accretion disks and jets revealed that the time variability of accretion rate show temporal $`1/f^\alpha `$ fluctuation (Kawaguchi et al. 2000, Sato et al. 2003), i.e., there are continuous distribution on time scales, which includes the one much longer than Keplerian time scale. The longest time scale of variability is probably related to duration of main accretion and/or dynamo in our magnetized accretion disks.

S. Woosley: Getting 10% of the accreted mass into the jet is too much in the context of GRBs. $`\mathrm{\Gamma }`$ of 2 may be OK initially but very important is the energy loading of that matter. Do you follow energy generation or are your jets all ”cold” ? Do you see the energy loading might be increased ?

K. Shibata: The jet with 10% of the accreted mass is a cold, dense jet directly ejected from cold part of the disk. We found another component of a jet ejected from hot corona with less material. Since our mechanism can produce enough poynting flux to explain total energy of GRBs, we suggest magnetic reconnection in the corona and jet might explain the energy loading. Such simulations of jets including reconnection in general relativistic regime should be done in future.

A. Hujeirat: The Blandford & Payne 1982 and my calculations predict the existence of a super-Keplerian transition layer between the disk and the corona. Why do your calculation don’t predict such a layer ?

K. Shibata: Our calculations also show such super-Keplerian layer just above the disk (see e.g., Shibata and Uchida 1986, Kudoh and Shibata 1997).

A. Brandenburg: Why does the magnetic field not penetrate the black hole horizon? This seems to be quite a general phenomena in relativistic electrodynamics. What is the force preventing this ?

K. Shibata: In our simulations, we adopt the coordinate system such that the time proceeds slowly near the event horizon, so that accreting plasmas and magnetic field lines have not yet penetrated into the event horizon in our results.

C. Fendt: Your statement about the highly intermittant jet formation was derived from ideal MHD simulations of the disk. What happens if you include turbulence or diffusion ? Will it stabilize the jet formation?

K. Shibata: Kuwabara et al. (2002) studied the effect of resistivity on jet formation. According to them, if the resistivity is large ($`R_m=rV_k/\eta <100`$), the jet formation is suppressed. In the interemediate regime ($`100<R_m<1000`$), the quasi-steady jet is formed. In the case of weak resistivity ($`R_m>1000`$), the results show formation of jets with high time variability and intermittency.


format: text

text:
# References

EVENT GENERATOR BENCHMARKING FOR
PROTON RADIOGRAPHY APPLICATIONS

Stepan G. Mashnik, Richard E. Prael, Arnold J. Sierk
Los Alamos National Laboratory, Los Alamos, NM 87545, USA

Konstantin K. Gudima
Institute of Applied Physics

Academy of Science of Moldova, Kishinev, MD-2028, Moldova

Nikolai V. Mokhov
Fermi National Accelerator Laboratory, Batavia, IL 60510, USA

Abstract

We have benchmarked the QGSM code and event generators of the MARS and LAHET3 codes as potential candidates for high-energy programs to be used in simulations for the Proton Radiography (PRad) Project. We have compiled from the literature experimental data on spectra of particles emitted from proton-induced reactions at incident energies from 30 GeV to 70 GeV on different nuclei and have performed calculations for all reactions for which we found data with these three codes without any modifications and using only default parameters and standard inputs. Our results (514 plots) show that all three codes describe reasonably most of the studied reactions, though all of them should be further improved before becoming reliable tools for PRad. We present here our conclusions concerning the relative roles of different reaction mechanisms in the production of specific secondary particles. We comment on the strengths and weaknesses of QGSM, MARS, and LAHET3 and suggest further improvements to these codes and to other models.

Introduction

The process of determining the feasibility of Proton Radiography (PRad) \[1-3\] as the radiographic probe for the Advanced Hydrotest Facility as well as its design and operation require information about spectra of secondary particles produced by high energy protons interacting in the target and structural materials. Reliable models and codes are needed to provide such data. We studied the literature and chose three potential candidates for high-energy codes that may be used in simulations for PRad, namely the Quark-Gluon String Model (QGSM) as developed by Amelin, Gudima, and Toneev , the MARS code by Mokhov et al. , and a version of the Los Alamos National Laboratory (LANL) transport code LAHET , known as LAHET3 .

The energy of the proton beam at PRad is supposed to be about 50 GeV. Unfortunately, there are very few measurements of particle spectra for proton-induced reactions exactly at 50 GeV or very close energies. In fact, we found only one published work at 50 GeV, namely spectra of $`\pi ^{}`$ and $`\pi ^+`$ measured at 159 from p(50 GeV) + W published in Russian together with pion spectra for other energies and targets, in a Joint Institute for Nuclear Research (Dubna) Communication by Belyaev et al. .

With only a few data available at 50 GeV, we benchmarked QGSM, MARS, and LAHET3 against measured spectra of particles emitted from interaction of protons with energies $`50\pm 20`$ GeV, i.e., from 30 to 70 GeV, with all targets for which we found experimental data. Independently of how many spectra were measured in an experiment, we performed calculations with the standard versions of QGSM, MARS, and LAHET3 without any modifications or adjustments, using only default parameters in the input of codes, and calculated double differential cross sections at 0, 4.75, 9, 13, 20, 45, 60, 90, and 159 degrees, angle-integrated energy spectra, and mean multiplicities for emission of n, p, d, t, 3He, 4He, $`\pi ^+`$, $`\pi ^{}`$, $`K^+`$, $`K^{}`$, and $`\overline{p}`$ for all cases listed below in Table 1. The next Section presents a brief description of the benchmarked codes, followed by results, discussion, and conclusions in the last two Sections.

Benchmarked Codes

QGSM: The core of the QGSM is built on a time-dependent version of the intranuclear cascade model developed at Dubna to describe both particle- and nuclei-induced reactions, often referred in the literature simply as the Dubna intranuclear Cascade Model (DCM) (see  and references therein). The DCM models interactions of fast cascade particles (“participants”) with nucleon spectators of both the target and projectile nuclei and includes interactions of two participants (cascade particles) as well. It uses experimental cross sections (or those calculated by the Quark-Gluon String Model for energies above 4.5 GeV/nucleon) for these elementary interactions to simulate angular and energy distributions of cascade particles, also considering the Pauli exclusion principle. When the cascade stage of a reaction is completed, QGSM uses the coalescence model described in  to “create” high-energy d, t, 3He, and 4He by final state interactions among emitted cascade nucleons, already outside of the colliding nuclei. After calculating the coalescence stage of a reaction, the QGSM moves to the description of the last slow stages of the interaction, namely to preequilibrium decay and evaporation, with a possible competition of fission using the standard version of the Cascade Exciton Model (CEM) . But if the residual nuclei have atomic numbers with $`A13`$, QGSM uses the Fermi break-up model to calculate their further disintegration instead of using the preequilibrium and evaporation models.

MARS: The MARS Monte-Carlo code system, being developed over 29 years, allows fast and reliable inclusive and exclusive simulation of three-dimensional hadronic and electromagnetic cascades in shielding, accelerator and detector components in the energy range from a fraction of an electron-volt up to about 100 TeV . It is under continuous development. The reliable performance of the code has been demonstrated in numerous applications at Fermilab, CERN, KEK and other centers as well as in special benchmarking studies. Description of elastic and inelastic $`hN`$, $`hA`$, $`\gamma A`$ and $`\nu A`$ cross sections is based on the newest compilations and parameterizations . At high energies (5 GeV$`<`$E$`<`$100 TeV), $`\sigma _{tot}`$, $`\sigma _{in}`$, $`\sigma _{prod}`$ and $`\sigma _{el}`$ are calculated in the framework of the Glauber multiple scattering theory with the $`\sigma _{hN}`$ as an input. The nucleon density distribution in nuclei is represented as the symmetrized Fermi function with the parameters of  for medium and heavy nuclei ($`Z>10`$) and the ones of  for $`Z<10`$. Modern evaluated nuclear data as well as fitting formulae are used to simulate hadron-nucleus elastic scattering. For protons, nuclear, Coulomb elastic scattering, and their interference is taken into account. At E$`>`$5 GeV, a simple analytical description used in the code for both coherent and incoherent components of $`d\sigma /dt`$ is quite consistent with experiment. A version of the Cascade-Exciton Model of nuclear reactions  as realized in the code CEM95  and containing also several recent refinements  is now implemented in the 1998 version of MARS  as default for 1-10 MeV $`<`$ E $`<`$ 3-5 GeV. A set of phenomenological models, as described in Ref. , is used for inclusive production of secondary particles in $`hA`$, $`dA`$, $`\gamma A`$ and $`\nu A`$ interactions at projectile energies from 5 GeV to 100 TeV. The 2001 version  of the MARS code was employed in the present benchmark.

LAHET3: LAHET is a Monte-Carlo code for the transport and interaction of nucleons, pions, muons, light ions, and antinucleons in complex geometry ; it may also be used without particle transport to generate particle production cross sections. LAHET allows one to choose one of several options for the Intra-Nuclear Cascade (INC) and fission models to be employed in calculations; it is widely used and well known in the applied nuclear physics community; therefore, we do not describe it here (a comprehensive description of LAHET may be found in  and references therein). The version of LAHET realized in the code LAHET3  uses a version of the code FLUKA, known in the literature as FLUKA96  to describe the first, INC stage of reactions, and its own Multistage Preequilibrium Model (MPM)  to describe the following intermediate preequilibrium stage, followed by evaporation/fission slow processes (or by the Fermi break-up model after the cascade instead of preequilibrium and evaporation/fission, if the residual nuclei have atomic numbers with $`A13`$ and for $`14A20`$ with excitation energy below 44 MeV), as described in . We mention again that only the high-energy event generator from FLUKA96 is employed here, as implemented in LAHET3; the default preequilibrium, evaporation and Fermi break-up models of LAHET3 are used for low energy nucleon and complex particle emission. More details and further references on LAHET3 with FLUKA96 can be found in .

Results and Discussion

Table 1 lists the cases we calculated with QGSM, MARS, and LAHET3, and provides references to experimental works where at least one spectrum of a secondary particle (from the ones listed in Introduction) was measured. A detailed report of the study containing 514 plots with spectra and multiplicities of secondary particles from reactions listed in Tab. 1 is now in preparation. Here, we present only our main conclusions and several

Table 1. Proton energy and target list covered by the present benchmark

| $`T_p`$ (GeV) | Nuclei | Measurements |
| --- | --- | --- |
| 30 | 9Be, 27Al | \[8, 22–24\] |
| 47 | 12C |  |
| 50 | 184W |  |
| 51 | 9Be, 48Ti |  |
| 53 | 27Al |  |
| 54 | 96Mo |  |
| 70 | 12C, 27Al, 64Cu, 118Sn, 208Pb | \[26–30\] |

examples of results from the study.

Our analyses have shown that all three codes tested here describe reasonably most of the secondary particle spectra. As a rule the higher the incident proton energy, the better the calculated spectra agree with experimental data. Several reaction mechanisms participate in the production of secondary nucleons and complex particles. These mechanisms are: 1) Fast INC processes; 2) preequilibrium emission from residual nuclei after the cascade stage of reactions; 3) evaporation of particles from compound nuclei after the preequilibrium stage, or from fission fragments, if the compound nucleus was heavy enough to undergo fission; 4) Fermi break-up of light excited nuclei formed after the cascade stage of reactions; 5) coalescence of complex particles by final state interactions among emitted cascade nucleons; 6) fast complex particle emission via knock out and pick up processes; 7) Multifragmentation of highly-excited residual nuclei after the INC. Their relative roles change significantly with the changing atomic mass number of the targets, and are different for different energies and angles of emission of secondary particles. Different codes describe these spectra better, worse, or do not describe them at all, depending of how these reaction mechanisms are (or are not) implemented into a specific code.

As an example, Fig. 1 shows spectra of p, d, t, and $`\pi ^{}`$ emitted at 9.17 deg from the reaction p(70 GeV) + 208Pb. Results for other reactions at 70 GeV are similar. One can see that all three codes describe the proton spectra well. The agreement for the pion spectra is not so good but is still reasonable, with some underestimation of the high-energy tails of spectra by QGSM and some overestimation by MARS. Note that as the angle of pion emission changes the situation is reverses: we observe that most of the high-energy tails of pion spectra at 159 deg, and to a lesser extent at 90 deg, are over-predicted by LAHET3 and underestimated by MARS.

The situation with the deuteron and tritium spectra is quite different. We see that deuterons with momentum of up to about 15 GeV/c and tritium with momenta up to 19 GeV/c are emitted and measured in this particular reaction. Utilizing the coalescence mechanism for complex particle emission, QGSM is able to describe high-energy deuteron production, and agrees well with the measurement. LAHET3 does not consider the coalescence of complex particles and therefore describes emission of only evaporative and preequilibrium deuterons with momenta not higher than 1 GeV/c. MARS does not consider emission of complex particles at such high incident proton energies, therefore no d and t spectra by MARS are shown in Fig. 1.

For tritium, the situation is worse since LAHET3, as is the case of deuterons, describes only preequilibrium emission and evaporation of tritons with momenta not higher than 1 GeV/c and QGSM, even taking into account coalescence of tritium, describes emission of t from this reaction up to only 2.5 GeV/c while the experimental spectrum of t extends to 19 GeV/c. This deficiency can be understood by considering the coalescence mechanism: It is much more probable to emit two cascade nucleons with very similar momenta that can coalesce into a deuteron than to get three INC nucleons with very similar momenta that can coalesce into a triton. The experimental values of high-energy triton spectra are several orders of magnitude below the corresponding values of the deuteron spectra, and the statistics of our present QGSM simulation could be simply too small to get such high-energy tritium via coalescence. There is also a possibility that knock out processes of preformed clusters (or fluctuations of nuclear density, leading to “fluctons”) by bombarding protons are seen in these experimental d and t spectra, but are not taken into account by any of the tested codes, providing the observed difference in the t spectrum and less pronounced, in the d spectrum. A third possible mechanism of complex particle emission with greater than 1 GeV/c momenta would be multifragmentation of highly-excited residual nuclei after the INC. This mechanism is not taken into account by any of the tested codes and we cannot estimate its contributation.

Fig. 2 shows examples of $`\pi ^+`$ spectra at 159 deg from 51 GeV proton collisions with 9Be and 48Ti. As already mentioned above for $`\pi ^{}`$, we see that LAHET3 overestimates the high-energy tails of pion spectra and MARS underestimates them a little. Similar results were obtained for other targets and incident proton energies.

Fig. 3 shows an example of how calculated proton spectra depend on the angle of emission, for the reaction p(30 GeV) + 9Be. We see that at 30 GeV, the agreement of calculated proton spectra with the data is not so good as we have in Fig. 1 for 70 GeV. The shapes and absolute values of proton spectra predicted by different codes depend significantly on the angle of detection, as does the agreement with the data. Similar results were obtained for other secondary particles and for other targets and incident energies.

Fig. 4 shows an example of mean multiplicity of secondary n, p, d, and $`\pi ^+`$ predicted by the tested codes for interaction of protons of about 50 GeV with different nuclei as functions of the mass number of targets. We see that predicted particle multiplicities differ significantly from each other, and the differences increase with increasing mass number of the target. The observed differences point to a quite significant difference in the treatment by the codes of the cascade stage of reactions (pions are emitted only at the cascade stage of reactions) and of the subsequent preequilibrium, evaporation, and Fermi break-up stages as well (we recall that at these incident energies MARS uses its own approximations for the total particle spectra without considering separately contributions from different mechanisms of nuclear reactions). These differences indicate that further experimental data are necessary at these incident proton energies and further development and improvement of the codes is required.

Further Work

Our study shows that all three codes describe reasonably well many of the secondary particle spectra analyzed here, though all of them should be further improved before becoming reliable tools for PRad.

For instance, we find that QGSM has some problems in a correct description of several pion spectra and does not describe sufficiently the high-energy tails of measured t and 3He spectra. Nevertheless, QGSM is the only code tested here that accounts for coalescence of complex particles from cascade nucleons and provides production of high-energy complex particles.

MARS overestimates the high-energy tails of some pion, kaon, and proton spectra at small angles (4.75, 9, and 13) and underestimates them a little at large angles (90 and 159). At these incident energies, MARS does not calculate complex particle production. However, MARS has one significant advantage in comparison with the two other codes: It is several orders of magnitude faster and requires almost no computing time, providing meanwhile reliable results for many applications.

LAHET3 overestimates the high-energy tails of practically all pion spectra at 159 and some nucleon and complex particle spectra in the preequilibrium energy region. It does not consider coalescence of complex particles and does not describe production of high-energy complex particles.

We observe also big differences between predicted high-energy tails of both neutron and proton spectra at 0 and for the mean multiplicities of almost all secondary particles, though no experimental data for these quantities are available at present for the reactions studied here.

We note that many of the problems we observe in our study for particular codes have already been solved, since all benchmarked event generators are under continuous development and improvement and all of them have been further improved in comparison with the versions we use in this study.

On the basis of QGSM, we have developed the Los Alamos version of the Quark-Gluon String Model code, LAQGSM . LAQGSM differs from QGSM by replacing the preequilibrium and evaporation parts of QGSM described according to the standard CEM  with the new physics from CEM2k  and has a number of improvements and refinements in the cascade and Fermi break-up models. Originally, both QGSM and LAQGSM were not able to describe fission reactions and production of light fragments heavier than 4He, as they had neither a high-energy-fission nor a fragmentation model. Recently, we addressed these problems  by further improving CEM2k and LAQGSM and by merging them with the Generalized Evaporation Model code GEM2 developed by Furihata . The improved LAQGSM+GEM2 code describes both spectra of secondary particles and yields of produced nuclides much better than QGSM does; exemplary results by LAQGSM and further references may be found in .

The MARS code system is being continuously developed and improved. For instance, a new version of the code, MARS14(2002) was completed after we started the present work. It contains a large number of improvements and refinements and provides better results in comparison with the version used here. Recently, the authors of MARS started to develop new and better approximations for the double differential cross sections of inelastic $`hN`$ and $`hA`$ interactions. The new systematics allow to solve the mentioned above problems with the pion, kaon, and proton spectra at forward and large angles and describe the experimental data much better.

The FLUKA code has also been updated very significantly (see e.g.,  and references therein) since the version FLUKA96 was incorporated into LAHET3 as used here; no updated version is yet incorporated into LAHET.

Our study points to the importance of taking into account coalescence in high-energy complex-particles production. We find it appropriate and easy to implement these processes into MARS and LAHET, as well as into other codes that do not now consider coalescence.

We think that at such high incident energies, multifragmentation of highly-excited heavy nuclei may also be of significance and should be taken into account in these event generators and in other codes.

Acknowledgment

The study was supported by the U. S. Department of Energy and by the Moldovan-U. S. Bilateral Grants Program, CRDF Project MP2-3025.


text:
# Smearing of the phase transition in Ising systems with planar defects

## 1 Introduction

Phase transitions in random systems remain one of the important and only partially solved problems in statistical physics. Originally, it was thought that even weak disorder destroys any critical point because the disordered system divides itself up into spatial regions which undergo the transition at different temperatures, leading to a smeared transition in which there would be no sharp singularities in thermodynamic quantities (see  and references therein). However, it was soon found that classical continuous phase transitions are generically sharp in the presence of weak, short-range correlated disorder. If a clean critical fixed point (FP) fulfills the Harris criterion  $`\nu 2/d`$, where $`\nu `$ is the correlation length exponent and $`d`$ is the spatial dimensionality, the disorder decreases under coarse graining. The system becomes asymptotically homogeneous at large length scales. Thus, the clean critical point is perturbatively stable against disorder, and the critical behavior of the random system is identical to that of the corresponding clean system. Macroscopic observables are self-averaging at the critical point, i.e., the relative width of their probability distributions goes to zero in the thermodynamic limit .

Even if the Harris criterion is violated, the transition will generically remain sharp, but the critical behavior will be different from that of the clean system. Depending on the fate of the disorder under coarse graining it can either be of finite or of infinite disorder type. In the first case the system remains inhomogeneous at all length scales with the relative strength of the inhomogeneities approaching a finite value for large length scales. These transitions are controlled by renormalization group (RG) fixed points with finite disorder. Macroscopic observables are not self-averaging, the relative width of their probability distributions approaches a size-independent constant . Examples of critical points in this class include the dilute three-dimensional Ising model or classical spin glasses. The second case occurs when the relative magnitude of the inhomogeneities increases without limit under coarse graining. The corresponding fixed points are called infinite-disorder fixed points. Here, the probability distributions of macroscopic variables become very broad (on a logarithmic scale) with the width increasing with system size. This can lead to activated rather than the usual power-law scaling at the critical point. A famous classical example of an infinite-disorder critical point occurs in the McCoy-Wu model .

In recent years, another aspect of phase transitions in disordered systems, the Griffiths phenomena, has regained a lot of attention. Griffiths phenomena are non-perturbative effects of rare disorder fluctuations in the vicinity of a phase transition. They can be understood as follows: Disorder in general decreases the critical temperature $`T_c`$ from its clean value $`T_c^0`$. In the temperature region $`T_c<T<T_c^0`$ the system does not display global order. However, in an infinite system one will find arbitrarily large regions that are devoid of impurities with a small but non-zero probability. These rare regions or ’Griffiths islands’ are locally in the ordered phase while the bulk system is still in the disordered phase. Their fluctuations are very slow because flipping them requires changing the order parameter in a large volume. Griffiths  was the first to show that this leads to a non-analytic free energy everywhere in the region $`T_c<T<T_c^0`$, which is now known as the Griffiths phase  or the Griffiths region. In generic classical systems, the contribution of the Griffiths singularities to thermodynamic (equilibrium) observables is very weak since the singularity in the free energy is only an essential one . In contrast, the consequences for the dynamics are much more severe with the rare regions dominating the behavior for long times. In the Griffiths region, the spin autocorrelation function $`C(t)`$ decays very slowly with time $`t`$, like $`\mathrm{ln}C(t)(\mathrm{ln}t)^{d/(d1)}`$ for Ising systems , and like $`\mathrm{ln}C(t)t^{1/2}`$ for Heisenberg systems . More recently, these results have been confirmed by rigorous calculations for the equilibrium  and dynamic  properties of disordered Ising models.

In many real systems, the disorder is not generated by point-like defects but by dislocations or disordered layers or grain boundaries. Often, extended impurities in a $`d`$-dimensional system can be modeled by disorder perfectly correlated in $`d_C=`$ 1 or 2 dimensions but uncorrelated in the remaining $`d_{}=dd_C`$ dimensions. While it is generally accepted that long-range disorder correlations increase the effects of the impurities, their influence on the phase transition has been controversial. Early RG work  based on a single expansion in $`ϵ=4d`$ did not produce a critical fixed point. This was interpreted as a smeared transition or a first order one . However, an additional expansion in the number $`d_C`$ of correlated dimensions  cured this problem and gave rise to critical fixed points with conventional power-law scaling. Further support for the sharp transition scenario came from Monte-Carlo simulations of a 3D Ising model with planar defects . All the perturbative RG studies miss, however, the effects of rare regions and Griffiths phenomena. These effects have been studied in detail in the already mentioned McCoy-Wu model , a disordered 2D Ising model in which the disorder is perfectly correlated in one dimension and uncorrelated in the other. The transition in this model was originally thought to be smeared. It was later found that it is sharp but of infinite-disorder type . Similar behavior has been found in several random quantum systems which are related to classical systems with linear defects . Based on these results it was believed for a long time that phase transitions generically remain sharp even in the presence of extended disorder.

In this paper we show, however, that this belief is not generally true. Specifically, we show that for Ising order parameter symmetry, planar defects destroy the sharp continuous phase transition via effects which are similar to but stronger than the usual Griffiths phenomena: In systems with planar defects true static long-range order can develop on isolated rare regions. As a result, the order parameter is very inhomogeneous in space close to the smeared transition, and the thermodynamic order parameter depends exponentially on temperature. We also show that this disorder-induced smearing generically occurs in an entire class of transitions in systems with extended disorder. The paper is organized as follows. In Sec. 2, the model is introduced, and the mechanism for the smearing of the transition is explained. In Sec. 3 we use Lifshitz tail arguments to determine the behavior in the ’tail’ of the rounded transition, i.e., for a very low concentration of ordered islands. Section 4 is devoted to a numerical study of a three-dimensional infinite-range Ising model which illustrates the smearing of the phase transition. Finally, in Sec. 5 we discuss the generality of the disorder-induced smearing, the relation to Griffiths phenomena, favorable conditions for observing the smearing, and the influence of the order parameter symmetry.

## 2 Rare regions, inhomogeneous order, and smearing

For definiteness we consider a $`d`$-dimensional $`\varphi ^4`$ theory with random-$`T_c`$ type disorder completely correlated in the $`d_C=2`$ directions $`x_1`$ and $`x_2`$ but uncorrelated in the remaining $`d_{}=dd_C`$ directions $`x_3,\mathrm{},x_d`$. The action is given by

$$S=d^dr\varphi (𝐫)\left[t+\delta t(𝐫_{})_𝐫^2\right]\varphi (𝐫)+ud^dr\varphi ^4(𝐫).$$
(1)
Here $`\varphi (𝐫)`$ is a scalar field, $`t`$ is the dimensionless distance from the clean critical point and $`\delta t(𝐫_{})`$ introduces the quenched disorder. $`𝐫_{}`$ is the projection of $`𝐫`$ on the uncorrelated directions $`x_3,\mathrm{},x_d`$. We consider two different types of disorder. In the first type, Poisson (or dilution) disorder, $`\delta t(𝐫_{})`$ is given by

$$\delta t(𝐫_{})=\underset{i}{}V[𝐫_{}𝐫_{}(i)]$$
(2)
where $`𝐫_{}(i)`$ are the random positions of planar impurities of spatial density $`c`$, and $`V(𝐫_{})`$ is a non-negative short-ranged impurity potential. The second type of disorder is a Gaussian distribution with zero mean and a correlation function

$$\delta t(𝐫_{})\delta t(𝐫_{}^{}{}_{}{}^{})=\mathrm{\Delta }^2\delta (𝐫_{}𝐫_{}^{}{}_{}{}^{}).$$
(3)
Now consider the effects of rare disorder fluctuations in this system: Analogous to the Griffiths phenomena, there is a small but finite probability for finding large spatial regions in $`𝐫_{}`$ direction which are more strongly coupled than the bulk system. For Poisson disorder, these are regions devoid of impurities while for Gaussian disorder it would be regions with a negative average $`\delta t(𝐫_{})`$. These rare regions can be locally in the ordered phase, even if the bulk system is still in the disordered phase. For Poisson disorder which fulfills $`\delta t(𝐫_{})>0`$ this starts to happen right below the clean transition, i.e., for $`t<0`$  while for the unbounded Gaussian disorder, locally ordered rare regions exist for all $`t`$. Since the defects in our system are planar, the rare regions are infinite in the $`x_1`$ and $`x_2`$ directions but finite in the remaining directions. This is a crucial difference to systems with uncorrelated disorder, where the rare regions are finite objects (see Sec. 5). In our system, each rare region is therefore equivalent to a two-dimensional Ising system and can undergo a real phase transition independently of the rest of the system. Thus, for planar defects, those rare regions which are in the ordered phase will develop a non-zero static order parameter which can be aligned by an infinitesimally small interaction or an infinitesimally small field.

The resulting phase transition is very different from conventional continuous phase transitions, be it transitions controlled by clean or finite-disorder or even infinite-randomness fixed points. In all of those cases, which represent “sharp” phase transitions, a non-zero order parameter develops as a collective effect of the entire system, accompanied by a diverging correlation length in all directions. In contrast, in a system with planar defects, the order parameter develops very inhomogeneously in space with different parts of the system (i.e., different $`𝐫_{}`$ regions) ordering independently at different $`t`$. Correspondingly, the correlation length in $`𝐫_{}`$ direction remains finite for all temperatures. This defines a smeared transition. We thus conclude that planar defects lead to a smearing of the phase transition.

Above, we have defined a smeared continuous phase transition via the behavior of correlation length. We now characterize it via the behavior of the free energy density or other global thermodynamic variables.111Note that global thermodynamic variables (which average over the entire sample) do not resolve the strong spatial inhomogeneities which are the hallmark of a smeared transition, and thus do not provide complete information about the transition Two qualitatively different cases of smeared transitions must be distinguished. If the disorder distribution is unbounded in the sense that it permits regions with a local $`T_c=\mathrm{}`$ (system (1) with Gaussian disorder falls into this class), the total (average) order parameter is non-zero for all temperatures, and the free energy density is analytic. In contrast, if the disorder distribution is bounded, there is a true paramagnetic phase with zero order parameter at high temperatures. At some temperature $`T_0`$, a non-zero order parameter starts to develop on rare spatial regions, accompanied by an essential singularity in the free energy density (which stems from the probability for finding a rare region). This singularity is weaker than the power-law singularities at the sharp transitions mentioned above. The model (1) with Poisson disorder falls into this second class, with $`T_0`$ being identical to the critical temperature of the pure system.

## 3 Lifshitz tail arguments

In this section we use extremal statistics to derive the leading thermodynamic behavior in the ’tail’ of the smeared transition, i.e., in the parameter region where a few islands have developed static order but their density is so small that they can be treated as independent. The approach is similar to that of Lifshitz  and others for the description of the tails in the electronic density of states. We first consider the Poisson type of disorder and later describe the differences for Gaussian disorder.

### 3.1 Poisson disorder

The probability $`w`$ of finding a large region of linear size $`L_R`$ (in $`𝐫_{}`$-space) devoid of any impurities is, up to pre-exponential factors, given by

$$w\mathrm{exp}(cL_R^d_{}).$$
(4)
As discussed in Sec. 2, such a rare region develops static long-range order at some $`t_c(L_R)`$ below the clean critical point $`t_c^0=0`$.222We use the bare value, $`t_c^0=0`$, for the clean critical reduced temperature. For our purpose, the distinction between the bare and the renormalized $`t`$ is of no importance. The value of $`t_c(L_R)`$ varies with the size of the region: The largest islands develop long-range order closest to the clean critical point. Finite size scaling (for the clean system, because the island is devoid of impurities) yields

$$t_c^0t_c(L_R)=|t_c(L_R)|=AL_R^\varphi $$
(5)
where $`\varphi `$ is the finite-size scaling shift exponent and $`A`$ is the amplitude for the crossover from $`d`$ dimensions to a slab geometry infinite in two dimensions but finite in $`d_{}=d2`$ dimensions. Combining (4) and (5) we obtain the probability for finding a rare region which becomes critical at $`t_c`$ as

$$w(t_c)\mathrm{exp}(B|t_c|^{d_{}/\varphi })(\text{for }t0)$$
(6)
where the constant $`B`$ is given by $`B=cA^{d_{}/\varphi }`$. The total (or average) order parameter $`m`$ is obtained by integrating over all rare regions which are ordered at $`t`$, i.e., all rare regions having $`t_c>t`$. Since the functional dependence on $`t`$ of the order parameter on a given island is of power-law type it does not enter the leading exponentials but only the pre-exponential factors. Therefore we obtain to exponential accuracy

$`m\mathrm{exp}(B|t_c|^{d_{}/\varphi })(\text{for }t0).`$ (7)
The homogeneous magnetic susceptibility consists of two different contributions, one from the islands on the verge of ordering and one from the bulk system still deep in the disordered phase. The bulk system provides a finite, noncritical background susceptibility throughout the tail region of the smeared transition. In order to estimate the contribution of the islands consider the onset of local ordering at the clean critical temperature $`t_c^0=0`$. Using (6) for the density of the islands we can estimate

$$\chi _0^\mathrm{\Lambda }𝑑tt^\gamma \mathrm{exp}(Bt^{d_{}/\varphi })$$
(8)
which is finite because the exponentially decreasing island density overcomes the power-law divergence of the susceptibility of an individual island. Here $`\gamma `$ is the clean susceptibility exponent and $`\mathrm{\Lambda }`$ is related to a lower cutoff for the island size. Once ordered islands exist they produce an effective background magnetic field everywhere in space which cuts off any possible divergence. Therefore we conclude that the homogeneous magnetic susceptibility does not diverge anywhere in the tail region of the smeared transition. However, there is an essential singularity at the clean critical point produced by the vanishing density of ordered islands.

The spatial magnetization distribution in the tail of the smeared transition is very inhomogeneous. On the already ordered islands, the local order parameter $`m(𝐫)`$ is of the same order of magnitude as in the clean system. Away from these islands it decays exponentially with the distance from the nearest island. The probability distribution $`P[\mathrm{log}m(𝐫)]`$ will therefore be very broad, ranging from $`\mathrm{log}m(𝐫)=O(1)`$ on the largest islands to $`\mathrm{log}m(𝐫)\mathrm{}`$ on sites very far away from an ordered island. The typical local order parameter $`m_{typ}`$ can be estimated from the typical distance of any point to the nearest ordered island. From (6) we obtain

$$r_{typ}\mathrm{exp}(B|t|^{d_{}/\varphi }/d_{}).$$
(9)
At this distance from an ordered island, the local order parameter has decayed to

$$m_{typ}e^{r_{typ}/\xi _0}\mathrm{exp}\left[C\mathrm{exp}(B|t|^{d_{}/\varphi }/d_{})\right]$$
(10)
where $`\xi _0`$ is the bulk correlation length (which is finite and changes slowly throughout the tail region of the smeared transition) and $`C`$ is constant. A comparison with (7) gives the relation between $`m_{typ}`$ and the thermodynamic order parameter $`m`$,

$$|\mathrm{log}m_{typ}|m^{1/d_{}}.$$
(11)
Thus, $`m_{typ}`$ decays exponentially with $`m`$ indicating an extremely broad local order parameter distribution. In order to determine the functional form of the local order parameter distribution, first consider a situation with just a single ordered island at the origin of the coordinate system. For large distances $`𝐫`$ the magnetization falls off exponentially like $`m(𝐫)=m_0e^{r_{}/\xi _0}`$. The probability distribution of $`x=\mathrm{log}[m(𝐫)]=\mathrm{log}m_0r_{}/\xi _0`$ can be calculated from

$$P(|x|)=\left|\frac{dN}{dx}\right|=\frac{dN}{dr_{}}\left|\frac{dr_{}}{dx}\right|=\xi _0\frac{dN}{dr_{}}\xi _0r_{}^{d1}$$
(12)
where $`dN`$ is the number of sites at a distance from the origin between $`r_{}`$ and $`r_{}+dr_{}`$ or, equivalently, having a logarithm of the local magnetization between $`x`$ and $`x+dx`$. For large distances we have $`|x|r_{}`$. Therefore, the probability distribution of $`\mathrm{log}m`$ generated by a single ordered island takes the form

$$P[\mathrm{log}(m)]|\mathrm{log}(m)|^{d_{}1}(\text{for }m1).$$
(13)
In the tail region of the smeared transition the system consists of a few ordered islands whose distance is large compared to $`\xi _0`$. The probability distribution of $`\mathrm{log}[m(𝐫)]`$ thus takes the form (13) with a lower cutoff corresponding to the typical island-island distance and an upper cutoff corresponding to a distance $`\xi _0`$ from an ordered island.

### 3.2 Finite-size effects

It is important to distinguish effects of a finite size $`L_C`$ in the correlated directions and a finite size $`L_{}`$ in the uncorrelated directions. If $`L_{}`$ is finite but $`L_C`$ is infinite static order on the rare regions can still develop. In this case, the sample contains only a finite number of islands of a certain size. As long as the number of relevant islands is large, finite size-effects are governed by the central limit theorem. However, for $`t0`$ very large and rare islands are responsible for the order parameter. The number $`N`$ of islands which order at $`t_c`$ behaves like $`NL_{}^d_{}w(t_c)`$. When $`N`$ becomes of order one, strong sample-to-sample fluctuations arise. Using (6) for $`w(t_c)`$ we find that strong sample to sample fluctuations start at

$`|t_L_{}|\left({\displaystyle \frac{d_{}}{B}}\mathrm{log}(L_{})\right)^{\varphi /d_{}}.`$ (14)
Thus, finite size effects are suppressed only logarithmically.

Analogously, one can study the onset of static order in a sample of finite size $`L_{}`$ (i.e., the ordering temperature of the largest rare region in this sample). For small sample size $`L_{}`$, the probability distribution $`P(t_c)`$ of the sample ordering temperatures $`t_c`$ will be broad because some samples do not contain any large islands. With increasing sample size the distribution becomes narrower and moves towards the clean $`t_c^0`$ because more samples contain large islands. The maximum $`t_c`$ coincides with $`t_c^0`$ corresponding to a sample without impurities. The lower cutoff corresponds to an island size so small that essentially every sample contains at least one of them. Consequently, the width of the distribution of critical temperatures in finite-size samples is governed by the same relation as the onset of the fluctuations,

$`\mathrm{\Delta }t_c\left({\displaystyle \frac{d_{}}{B}}\mathrm{log}(L_{})\right)^{\varphi /d_{}}.`$ (15)
We now turn to finite-size effects produced by a finite extension $`L_C`$ of the system in the correlated directions. In this case, the true static order on the rare regions is destroyed. Whether or not the system develops long-range order depends on the number $`d_{}`$ of uncorrelated dimensions. If $`d_{}=1`$ no static long-range order can develop, i.e. the transition is rounded by conventional finite size effects in addition to the disorder-induced smearing. In contrast, for $`d_{}>1`$ a true phase transition is possible but requires a finite interaction between the islands of the order of the temperature. This restores a sharp phase transition at a reduced temperature $`t_c(L_C)`$. To estimate the relation between $`L_C`$ and $`t_c(L_C)`$ we note that the interaction between two planar rare regions of linear size $`L_C`$ is proportional to $`L_C^2`$ and decays exponentially with their spatial distance $`r`$, $`E_{int}L_C^2\mathrm{exp}(r/\xi _0)`$, where $`\xi _0`$ is the bulk correlation length. With the typical distance given by (9) we obtain a double exponential dependence between $`L_C`$ and the critical reduced temperature $`t_c(L_C)`$:

$`\mathrm{log}(\mathrm{log}L_C)`$ $``$ $`|t_c|^{d_{}/\varphi }.`$ (16)
### 3.3 Gaussian disorder

In this subsection we apply Lifshitz tail arguments to the case of $`\delta `$-correlated Gaussian disorder . In contrast to the positive definite Poisson disorder, the Gaussian disorder is unbounded for $`\delta t(𝐫)\mathrm{}`$. Therefore, locally ordered islands can exist for all $`t`$, and, in principle, the tail region of the smeared transition stretches to $`t\mathrm{}`$. In order to determine the probability for finding a rare region which orders at a certain $`t_c`$, i.e., the Gaussian equivalent to (6), consider an island of linear size $`L_R`$. The average disorder value $`\delta t`$ in this region has the distribution

$$P[\delta t]\mathrm{exp}\left[L_R^d_{}\delta t^2/(2\mathrm{\Delta })\right].$$
(17)
According to finite size scaling, an island of size $`L_R`$ with an average disorder value of $`\delta t`$ will develop static order at $`t_c=\delta tAL_R^\varphi `$. For fixed island size, the probability distribution of $`t_c`$ reads

$$P(t_c)\mathrm{exp}\left[L_R^d_{}(t_c+AL_R^\varphi )^2/(2\mathrm{\Delta })\right].$$
(18)
We now calculate the stationary point of the exponent with respect to $`L_R`$ (saddle point approximation) to determine which island size gives the dominating contribution. For $`(2d_{}/\varphi )>0`$, this leads to

$$(2d_{}/\varphi )AL_r^\varphi =t_cd/\varphi .$$
(19)
This implies that small but very deep fluctuations are responsible for the tail of the smeared transition (in contrast to Poisson disorder where the tail is produced by the largest islands). Inserting into (18) gives the desired probability as a function of $`t_c`$,

$$w(t_c)\mathrm{exp}(\overline{B}t_c^{2d_{}/\varphi })$$
(20)
where $`\overline{B}`$ is a constant. In contrast, for $`(2d_{}/\varphi )<0`$, the main contribution comes from arbitrarily small islands. In this case the $`\delta `$-correlated disorder is unphysical and should be replaced by a distribution with a finite correlation length. The resulting $`w(t_c)`$ is then purely Gaussian. In the following we will not consider this case.

Starting from (20), one can derive the leading behavior in the tail of the smeared transition ($`t\mathrm{}`$). Integrating over all ordered islands gives the total magnetization

$$m\mathrm{exp}(\overline{B}t_c^{2d_{}/\varphi }).$$
(21)
The typical local magnetization $`m_{typ}`$ decays exponentially with the average magnetization $`m`$, i.e. (11) also holds in the Gaussian case as can be easily seen by starting the derivation from (20) rather than (6). For samples of finite size in the uncorrelated directions the onset of sample-to-sample fluctuations and the width of the sample ordering temperatures follow

$$t_L_{}\mathrm{\Delta }t_c(\mathrm{log}L_{})^{1/(2d_{}/\varphi )}.$$
(22)
If the samples are of finite size in the correlated directions but infinite in the uncorrelated ones (and if $`d_{}>1`$) a sharp phase transition is restored. The relation between $`L_C`$ and $`t_c`$ is double-exponential:

$$\mathrm{log}(\mathrm{log}L_C)t_c^{2d_{}/\varphi }.$$
(23)
A comparison of the results for Gaussian and Poisson disorder shows that the functional form of the thermodynamic relations at the smeared transition depends on the form of the disorder and is thus not universal.

## 4 Numerical results

### 4.1 Model and method

To illustrate the smearing of the phase transition we now show results from a numerical simulation of a 3D Ising model with planar bond defects, i.e., $`d_C=2`$ and $`d_{}=1`$. As discussed above, the smearing of the phase transition is the result of exponentially rare events. Therefore large system sizes are required to observe it. In order to reach sufficiently large sizes with a reasonable effort and to retain the possibility of static order on an isolated rare region, we consider a model with infinite-range interactions in the correlated directions (parallel to the defects) but short-range interactions in the uncorrelated direction (perpendicular to the defect planes). While the infinite-range model will not be quantitatively comparable to a short-range model, it provides a simple example for the rounding mechanism introduced in Sec. 2. Moreover, since the rounding mechanism only depends on the existence of static order on the rare regions but not on any details, we expect that the results will qualitatively valid for a short-range model, too (with the appropriate changes to the exponents in the relations derived in Sec. 3). Indeed, preliminary results  of Monte-Carlo simulations of a 3D short-range Ising model with planar defects fully agree with this expectation.

The Hamiltonian considered in this section reads

$$H=\frac{1}{L_C^2}\underset{x,y,z,y^{},z^{}}{}S_{x,y,z}S_{x+1,y^{},z^{}}\frac{1}{2L_C^2}\underset{x,y,z,y^{},z^{}}{}J_xS_{x,y,x}S_{x,y^{}z^{}}h\underset{x,y,z}{}S_{x,y,x}.$$
(24)
Here $`x,y,z`$ are the integer-valued coordinates of the Ising spins. The magnetic field $`h`$ acts as a symmetry breaker, it will be set to a small but non-zero value, typically $`h=10^{10}`$. $`L_C`$ is the system size in the $`y`$ and $`z`$ directions. The planar defects are parallel to the $`(y,z)`$-plane. They are introduced via disorder in the in-plane coupling $`J_x`$ which is a quenched binary random variable with the distribution $`P(J)=(1c)\delta (J1)+c\delta (JJ_d)`$ with $`J_d<1`$. Here, $`c`$ is the defect concentration and the deviation of $`J_d`$ from 1 measures the strength of the defects. For our calculations we use $`J_d=0`$, i.e., very strong defects. The fact that one can independently vary concentration and strength of the defects in an easy way is the main advantage of this binary disorder distribution. However, it also has unwanted consequences, viz. log-periodic oscillations of many observables as functions of the distance from the critical point . These oscillations are special to the binary distribution and unrelated to the smearing considered here; we will not discuss them further.

Because the interaction is infinite-ranged in the correlated ($`y`$ and $`z`$) directions, these dimensions can be treated exactly: We introduce the plane magnetizations

$$m_x=\frac{1}{L_C^2}\underset{y,z}{}S_{x,y,z}$$
(25)
into the partition function and integrate out the original Ising variables $`S_{x,y,z}`$. In the thermodynamic limit, $`L_C\mathrm{}`$, the resulting integral can be solved using saddle point approximation. This leads to a set of coupled local mean-field equations

$$m_x=\mathrm{tanh}\left[(m_{x1}+J_xm_x+m_{x+1}+h)/T\right],$$
(26)
where $`T`$ is the temperature (we work in units where the Boltzmann constant is $`k_B=1`$). In the clean case, $`c=0`$, all $`J_x`$ are equal to 1. Consequently, all local mean-field equations are identical, reading $`m=\mathrm{tanh}[(3m+h)/T]`$. The clean critical temperature of the model (24) is therefore $`T_c^0=3`$.

In the presence of disorder, the local mean-field equations (26) have to be solved numerically. We employ a simple self-consistency cycle: All $`m_x`$ are initially set to $`m_x=0.1`$. We then insert the current values of $`m_x`$ into the r.h.s. of (26) and calculate new values $`m_x^{\mathrm{new}}`$ which are then put into the r.h.s. and so on. This is repeated until the change $`|m_x^{\mathrm{new}}m_x|`$ falls below a threshold of $`10^{12}`$ for all $`x`$. This method is very robust, it never failed to converge, but it can be very slow in the tail region of the transition where ordered islands coexist with large disordered regions. In this region sometimes several 1000 iterations were necessary. However, in contrast to more sophisticated methods like the Newton-Raphson method no matrix inversions are required. Therefore we were able to simulate very large systems up to $`L_{}=10^6`$.

We also note that the model (24) is similar to the mean-field McCoy-Wu model considered in Ref. . The raw data obtained in this paper ($`L_{}`$ up to 1000) are very similar to ours, but they were interpreted in terms of conventional power-law critical behavior with an unusually large order parameter exponent $`\beta 3.6`$. This may be partially due to an unfortunate choice of parameters (in particular, a higher impurity concentration of $`c=0.5`$) which makes it harder to observe the exponential tail of the smeared transition. We will come back to this question in Sec. 5.

### 4.2 Thermodynamic magnetization and susceptibility

In this subsection we present numerical results for the total magnetization $`m=L_{}^1_xm_x`$ and the homogeneous susceptibility $`\chi =m/h`$. The left panel of figure 1 gives an overview over the behavior of $`m`$ and $`\chi `$ as functions of temperature $`T`$ for size $`L_{}=1000`$ and impurity concentration $`c=0.2`$. The data are averages over 1000 disorder realizations. (Thermodynamic quantities involve averaging over the whole system. Thus ensemble averages rather than typical values give the correct infinite system approximation.)

At a first glance these data suggest a sharp phase transition close to $`T=2.96`$. However, a closer inspection of the transition region in the right panel of figure 1 shows that the singularities are rounded. If this smearing was a conventional finite-size effect the magnetization curve should become sharper with increasing $`L_{}`$ and the susceptibility peak should diverge. This is not the case here, both the magnetization and the susceptibility data are essentially size-independent. We conclude that the smearing is an intrinsic effect of the infinite system.

In order to compare the total magnetization data with the analytical Lifshitz tail prediction (7) we need to determine the value of the clean finite-size scaling shift exponent for the model (24). Finite-size scaling is governed by the first appearance of a non-vanishing solution of the linearized clean mean-field equation in a slab geometry of varying thickness. This equation is equivalent to a one-dimensional Schrödinger equation in a potential well, leading to a clean shift exponent of $`\varphi =2`$ and $`d_{}/\varphi =1/2`$. In the left panel of figure 2, we therefore plot the logarithm of the total magnetization, averaged over 300 samples, as a function of $`(T_c^0T)^{1/2}`$ for system size $`L_{}=10000`$ and several impurity concentrations $`c`$.

For all $`c`$, the data follow eq. (7) over several orders of magnitude in $`m`$. Fits of the data to (7) are used to determine the decay constants $`B`$. The right panel of figure 2 shows that these decay constants depend linearly on $`\mathrm{log}(1c)`$. This is exactly the expected behavior in a lattice model with binary disorder since the probability for finding an island of size $`L_R`$ devoid of impurities behaves like $`(1c)^{L_R}e^{L_R\mathrm{log}(1c)}`$.

### 4.3 Local magnetization

In the tail region of the smeared transition the system consists of a few rare ordered regions (on sufficiently large islands devoid of impurities). These rare regions are far apart, and between them the local magnetization is exponentially small. This behavior is illustrated in figure 3.

The left panel shows the local magnetization of a segment of a large sample in the tail region of the smeared transition. The data show that the local magnetization on the islands is large ($`0.1`$) but it drops very rapidly with increasing distance from the islands. This drop-off can be used to estimate the bulk correlation length which is very small in this parameter region, $`\xi _03\mathrm{}4`$. Note that the seeming saturation of the magnetization decrease at $`m10^9`$ is a result of the finite external field $`h=10^{10}`$. Without field, the magnetization would drop much further. A comparison of the local magnetization $`m_x`$ and the local coupling constant $`J_x`$ in the right panel of figure 3 shows that magnetic order only exists on a sufficiently large island devoid of any impurities.

In order to quantify these observations, we have also calculated the probability distribution $`P(\mathrm{log}m_x)`$ of the local magnetization values, or, more precise, their logarithms. Figure 4 shows this distribution for a single large sample of $`L_{}=200000`$, impurity concentration $`c=0.2`$, and temperatures ranging from $`T=2.8`$ (deep in the ordered phase) to $`T=2.99`$ (in the tail of the smeared transition).

Clearly, with increasing temperatures the distribution becomes very broad, even on this logarithmic scale. Two qualitatively different regimes can be distinguished. For temperatures not too close to the clean critical point ($`T=2.80`$ to 2.96), the distribution is dominated by an exponential decay towards small magnetization values (i.e., towards negative $`\mathrm{log}m`$) with the decay constant decreasing with increasing temperature. However, once the system enters the rare-region dominated tail region of the smeared transition (which, as can be seen from figures 1 and 2, happens roughly at $`T=2.96\mathrm{}2.97`$ for $`c=0.2`$) the behavior changes: In agreement with (13) the probability distribution $`P(\mathrm{log}m_x)`$ becomes constant except for a peak at large magnetization values which we attribute to the behavior of $`m_x`$ on the ordered islands (eq. (13) was derived from the behavior of the sites not belonging to one of the rare ordered islands).

### 4.4 Finite-size effects and sample-to-sample fluctuations

In this subsection we present numerical data for the effects of a finite sample size $`L_{}`$ in the uncorrelated direction, i.e., perpendicular to the impurities. Effects of a finite size $`L_C`$ in the correlated directions cannot easily be studied within our approach because the limit $`L_C\mathrm{}`$ has been taken in deriving the local mean-field equations (26).

To investigate sample-to-sample fluctuations in a finite-size system we compare the average sample magnetization $`m_{av}=m`$ and the typical sample magnetization $`m_{typ}=\mathrm{exp}\mathrm{log}m`$ where $`m`$ is the total magnetization of a sample and $``$ denotes the average over the disorder. A significant difference between average and typical magnetizations indicates strong sample-to-sample fluctuations. The left panel of figure 5 shows the typical magnetizations for four different sizes $`L_{}`$ together with the average magnetization (which is size-independent within the accuracy of this plot) for 300 to 1000 samples, depending on $`L_{}`$.

In the tail region of the smeared transition the typical and average magnetizations deviate from each other by several orders of magnitude. From these data we define a temperature $`T_L_{}`$ which describes the onset of the sample-to-sample fluctuations via $`m_{av}(T_L_{})/m_{typ}(T_L_{})=\mathrm{e}`$. This temperature shifts to larger values with increasing system size $`L_{}`$. The inset shows $`t_L_{}=T_c^0T_L_{}`$ as a function of $`(\mathrm{log}L_{})^2`$. The plot shows that $`t_L_{}`$ follows eq. (14) in good approximation.

The right panel of figure 5 shows the probability distribution of the sample critical temperature which is defined as the ordering temperature of the largest rare region in a particular sample (600 to 2000 samples used). As predicted in subsection 3.2, the distribution becomes narrower with increasing size $`L_{}`$, and it shifts towards the clean critical temperature $`T_c^0=3`$. The width $`\mathrm{\Delta }T_c`$ of this distribution has been determined from the values where $`P`$ has decreased by a factor of 1/e. The inset shows that $`\mathrm{\Delta }T_c`$ is proportional to $`(\mathrm{log}L_{})^2`$ as predicted in eq. (15).

## 5 Discussion

To summarize, we have shown that true static order can develop on an isolated rare region in an Ising model with planar defects. As a result, different parts of the system undergo the phase transition at different temperatures, i.e., the sharp transition is smeared by the defects. In this final section we discuss the generality of our findings and their relation to Griffiths phenomena and the Harris criterion. We also discuss favorable conditions for observing the smearing in experiments or simulations.

The origins of the smearing of the phase transition introduced in this paper and of Griffiths phenomena are very similar, both are caused by rare large spatial regions which are locally in the ordered phase. The difference between the two effects is a result of disorder correlations. For uncorrelated or short-range correlated disorder, the rare regions are of finite size. Thus, they cannot develop true static long-range order. Instead, the order parameter fluctuates slowly leading to the Griffiths singularities  discussed in Sec. 1. In contrast, if the rare regions are infinite in at least two dimensions, a stronger effect occurs: the rare regions which are locally in the ordered phase actually develop true static order, and this leads to a smeared transition. In other words, the same rare regions which would usually produce Griffiths singularities are responsible for the smearing of the transition in a system with planar defects. The tail of the smeared transition stretches to where the Griffiths temperature would normally be, i.e., to $`T=\mathrm{}`$ for Gaussian disorder or $`T=T_c^0`$ for Poisson disorder. So, the smeared transition replaces not only the conventional critical point but also the Griffiths phase.

We now discuss under which conditions a phase transition will be smeared by disorder. The arguments given in Sec. 2 in favor of a smeared transition only relied on the dimensionality of the rare regions and the order parameter symmetry. It is therefore expected that all phase transitions with discrete order parameter symmetry are smeared by quenched disorder if the defects are perfectly correlated in at least two dimensions. In systems with linear defects the transition will remain sharp if the interactions are short-ranged. However, if the interactions in the correlated direction fall of as $`1/r^2`$ or slower, even linear defects can lead to smearing because a 1D Ising model with $`1/r^2`$ interaction has an ordered phase . This is particularly important  for quantum phase transitions in itinerant electronic systems which can be mapped onto classical systems with $`1/\tau ^2`$ interaction in imaginary time direction.

Phase transitions with continuous order parameter symmetry are more stable against smearing. In systems with short-range interactions the dimensionality of the defects has to be at least three, and in systems with linear or planar defects long-range interactions are required for smearing. It is known  that classical XY and Heisenberg systems in dimensions $`d=1,2`$ develop long-range order at finite $`T`$ only if the interaction falls off more slowly than $`1/r^{2d}`$. Consequently, in a system with linear defects the phase transition will only be smeared if the interactions in defect direction fall off more slowly than $`1/r^2`$ (or more slowly than $`1/r^4`$ for planar defects).

A third important question concerns the universality of the thermodynamic behavior in the vicinity of the smeared transition. As can be seen from comparing the results of Poisson and Gaussian disorder, the functional dependence of the magnetization and other observables on the temperature is not universal, it depends on details of the disorder distribution. Therefore, only the presence or absence of smearing is universal in the sense of critical phenomena (i.e., depending on dimensionality and symmetry only) while the thermodynamic relations are non-universal.

Our fourth remark deals with the relation of the disorder-induced smearing and the Harris criterion which in the case of extended defects reads  $`\nu >2/d_{}`$. We emphasize that the phase transition can be smeared by planar defects even if the corresponding clean critical point fulfills the Harris criterion and appears to be stable. The reason is that the Harris criterion assumes a homogeneous transition and studies the behavior of the coarse-grained (root-mean-square) disorder at large length scales. However, the formation of static order on an isolated finite-size rare region is a non-perturbative finite-length scale effect in the tail of the disorder distribution. This type of effects is not covered by the Harris criterion.

We now turn to the question which conditions are favorable for the observation of the smearing in experiments or simulations. A significant number of large rare regions will only exist if the defect concentration is small, as can be seen from (6). On the other hand, the impurities have to be sufficiently strong so that the bulk system is still far away from criticality when the first rare regions start to order. Thus, the most favorable type of disorder for the observation of the smearing is a small density of strong impurities. In contrast, for weak impurities the rounding is restricted to a very narrow temperature range below the clean critical temperature and maybe masked by the remainder of the bulk critical fluctuations. If the impurity concentration is too high, the exponential drop-off is very fast making it very hard to observe the smearing. This may also explain why no smearing was observed in earlier simulations. Specifically, in Ref.  weak impurities ($`\mathrm{\Delta }J/J=0.1`$) of a high concentration $`c=0.5`$ were used. These are unfavorable conditions, and the maximum system size of $`L=27`$ used in this simulation was probably too small to observe the smearing.

Finally, we add one remark about the extreme tail of the smeared transition. In this parameter region, the ordered islands are very far apart and the aligning effect of the direct island-island interactions is very weak. It is therefore possible that in a real system other interactions play an important role in determining the relative orientation of the order parameter on the islands. This could lead, e.g., to spin glass order in the extreme tail of the smeared transition.

In conclusion, perfect disorder correlations in one or more directions dramatically increase the effects of the disorder on a phase transition with discrete order parameter symmetry. Zero-dimensional defects typically lead to a conventional critical point with power-law scaling. For linear defects the generic behavior seems to be an infinite-randomness critical point  while we have shown here that planar defects generically destroy the sharp transition by smearing.

We acknowledge support from the University of Missouri Research Board. Part of this work has been performed at the Aspen Center for Physics.

## References