text:
# Thermal X-ray Emission and Cosmic Ray Production in Young Supernova Remnants

## 1 Introduction

The strong forward and reverse shocks in young supernova remnants (SNRs) heat the ambient interstellar medium (ISM) and supernova ejecta material to X-ray emitting temperatures. This X-ray emission contains information on the supernova (e.g., type, explosion energy, mass and composition of ejecta) and on the nature of the ambient medium (e.g., uniform, stellar wind).

In addition to heating the plasma, the forward and reverse shocks are thought to accelerate some fraction of the shocked material to cosmic ray energies. A large body of evidence suggests that this acceleration can be quite efficient with as much as 50% of the ram kinetic energy of the ejecta being placed in relativistic particles and removed from the thermal plasma (e.g., Blandford & Eichler ; Jones & Ellison ; Kang & Jones ; Berezhko, Ksenofontov, & Petukhov ), and some early work investigated its observational impact (Heavens , Boulares & Cox , Dorfi & Böhringer ). However, despite the expected efficiency of shock acceleration, the modeling and interpretation of X-ray observation from SNRs are generally done assuming that the shocks do not place a significant fraction of their energy in cosmic rays.

The recent discovery of nonthermal X-ray emission in shell-like SNRs like SN1006 (Koyama et al. ), Cas A (Allen et al. ), G347.3-0.5 (Slane et al. ) and RCW86 (Borkowski et al. ) has shown that SNR shocks do accelerate electrons to TeV energies and supports the case for efficient proton acceleration and the necessity of investigating its effects in the modeling of thermal X-ray emission.

Besides providing a more self-consistent SNR model, X-ray models of cosmic ray modified SNRs provide a test of the fundamental assumption that SNRs are the primary sources of galactic cosmic ray ions. Since shocks put more energy into accelerated ions than electrons, nonlinear effects seen in X-ray emission will be evidence for the efficient shock acceleration of ions as well.

Here, we investigate the effects of efficient particle acceleration on the hydrodynamics and thermal X-ray emission of young SNRs. We have developed a model which couples self-similar hydrodynamics (Chevalier ) to nonlinear diffusive shock acceleration (Berezhko & Ellison ; Berezhko et al. ; Ellison, Berezhko, & Baring ), and X-ray emission including non-equilibrium ionization effects (Decourchelle & Ballet ). We show that a modified self-similar approach is relevant in most cases and compare it to the standard unmodified test-particle case.

## 2 Nonlinear Shock Acceleration

There is convincing direct and indirect evidence from spacecraft observations and plasma simulations that collisionless shocks are efficient nonlinear accelerators. Consequently, they self-generate magnetic turbulence and develop a smooth shock precursor (e.g., Drury ; Blandford & Eichler ; Jones & Ellison ; Achterberg, Blandford, & Reynolds ; Terasawa et al. ; Ellison et al. ). When the acceleration is nonlinear, the shock compression ratio increases and the shocked temperature decreases compared to test-particle values. The entire proton spectrum must be accounted for self-consistently with a non-thermal tail connecting the quasi-thermal population to the energetic one; any change in the relativistic population must impact the thermal population. In contrast, the slope and normalization of test-particle power laws can be changed without changing the temperature and density of the shock heated gas. A description of nonlinear effects and the model used here is given in Berezhko & Ellison (), Ellison et al. ().

## 3  Modified hydrodynamic model

### 3.1 Model

The evolution of young SNRs can be described by self-similar solutions (Chevalier ; Nadezhin ) provided that the initial density profiles in the ejecta and in the ambient medium have power law density distributions, respectively $`\rho =g^nr^nt^{n3}`$ (with $`n>5`$) and $`\rho =qr^s`$ (with $`s<3`$), where $`r`$ is the radius, $`t`$ the SNR age and $`g`$ and $`q`$ are constant. $`g`$ can be expressed in terms of the explosion energy and ejected mass (Decourchelle & Ballet ). The effects of cosmic ray pressure on SNR dynamics were first considered by Chevalier () using a two-fluid, self-similar solution with a thermal gas (adiabatic index $`\gamma =5/3`$) and a relativistic gas ($`\gamma =4/3`$) having an arbitrary constant fraction of the total energy. Here, we use Chevalier’s () solutions for driven waves but with boundary conditions calculated from the shock acceleration model from Berezhko & Ellison (1999).

As a function of shock velocity, $`V_\mathrm{s}`$, ambient density, $`\rho _0`$, unshocked magnetic field, $`B_0`$, and ambient temperature $`T_0`$, as well as the maximum particle energy, $`E_{\mathrm{max}}`$ 111$`E_{\mathrm{max}}`$ is determined from the shock radius, age, and diffusive mean free path (e.g., Ellison, Berezhko, & Baring ). and injection efficiency, $`\eta _{\mathrm{inj},\mathrm{p}}`$ (i.e., the fraction of protons with superthermal energies), the nonlinear model yields downstream values for the density and pressures in the thermal and relativistic gases, at the forward and reverse shocks. Alfvén wave heating in the precursor is used, reducing the efficiency compared to adiabatic heating. Starting with test-particle solutions at a given age, we iterate until the modified solutions and nonlinear boundary conditions are satisfied.

Our solutions neglect cosmic ray diffusion (assuming they are spatially coupled to the gas – an excellent approximation for all but the highest energy particles), and assume that both the compression ratio and the relativistic to total pressure ratio are independent of time at either shock, which is required for self-similarity to be valid.

Unless explicitly stated, we take $`T_0=10^4`$ K and $`B_0=5`$ $`\mu `$G at both shocks. For the SNR parameters, we take: $`g=3.410^9`$ CGS, corresponding to $`E_{\mathrm{SN}}(M_{}/M_{\mathrm{ej}})^{2/3}=0.5310^{51}`$ erg, and $`q=\rho _0=0.42`$ amu/cm3.

### 3.2 Validity of the self-similar approach

We can a posteriori determine the validity of the self-similar assumption, by examining the variation versus time of the compression ratio and relativistic to total pressure ratio, at the reverse shock (RS) and at the forward shock (FS). As seen in Figure 1, these ratios are nearly constant for injection parameters $`\eta _{\mathrm{inj},\mathrm{p}}\begin{array}{c}>\hfill \\ \hfill \end{array}210^4`$ and negligible for $`\eta _{\mathrm{inj},\mathrm{p}}\begin{array}{c}<\hfill \\ \hfill \end{array}10^5`$, so that the self-similar solutions are well justified in these cases.

For values of $`\eta _{\mathrm{inj},\mathrm{p}}\begin{array}{c}<\hfill \\ \hfill \end{array}2\times 10^4`$, unmodified solutions can occur at very high sonic Mach numbers $`M_\mathrm{s}`$ (i.e., early ages) which show a rapid transition to nonlinear solutions as $`M_\mathrm{s}`$ drops (see Berezhko & Ellison ). The self-similar approach is valid if this rapid transition occurs early or late compared to the age of the SNR.

### 3.3 Results on the dynamics

Figure 2 shows the density and temperature as a function of radius for various values of $`\eta _{\mathrm{inj},\mathrm{p}}`$. As $`\eta _{\mathrm{inj},\mathrm{p}}`$ increases, the compression ratio becomes higher (up to 15 here instead of 4), while the shocked temperature drops (nearly a factor of 10 between the test-particle result and $`\eta _{\mathrm{inj},\mathrm{p}}=10^4`$). Such trends were expected from Chevalier (), but the actual nonlinear effects lead to larger modifications because of the influence of particle escape: the maximum compression ratio is no longer limited to 7, but can rise, in principle, to infinite values.

Notably, these modifications occur for $`\eta _{\mathrm{inj},\mathrm{p}}\begin{array}{c}>\hfill \\ \hfill \end{array}5\times 10^5`$, considerably lower than that deduced for the Earth’s bow shock and from plasma simulations (i.e., $`\eta _{\mathrm{inj},\mathrm{p}}10^3`$; Ellison et al. ). Table 1 gives the relative fractions of kinetic, $`E_{\mathrm{kin}}`$, thermal, $`E_{\mathrm{th}}`$, cosmic ray, $`E_{\mathrm{cr}}`$, and escaped energy, $`E_{\mathrm{esc}}`$, summed over the shocked ejecta and the shocked ambient medium. An important fraction of the energy can be carried off by escaping particles. Note that in young SNRs, the total energy in the shocked region increases with time as more and more ejecta is being shocked.

Cosmic ray modified dynamics lead to a thinner interaction region between the forward and reverse shocks and much larger density and pressure gradients. This suggests that the growth of hydrodynamic instabilities in the shocked material near the contact discontinuity will be faster than otherwise. With a forward shock closer to the contact discontinuity, this may help these turbulent effects extend to the forward shock and beyond.

Equally important, the higher densities and lower temperatures will lead to more rapid Coulomb equilibration than expected in a test-particle case.

### 3.4 Velocity versus shock temperature

Table 1. Temperature-Velocity Relation and energy fractions      $`\eta _{\mathrm{inj},\mathrm{p}}`$ $`x`$ $`E_{\mathrm{kin}}`$ $`E_{\mathrm{th}}`$ $`E_{\mathrm{cr}}`$ $`E_{\mathrm{esc}}`$    RS $`10^5`$ 0.17 0.85 0.13 0.01 0.00   $`10^4`$ 0.02 0.80 0.02 0.08 0.10   $`10^2`$ 0.01 0.78 0.01 0.08 0.12  FS $`10^5`$ 0.18 0.48 0.48 0.03 0.00   $`10^4`$ 0.02 0.42 0.04 0.26 0.28   $`10^2`$ 0.01 0.41 0.01 0.24 0.33

We emphasize that when particle acceleration is efficient, the relation between the shock velocity and the mean post-shock temperature of the shocked gas is no longer approximated by the expression for a strong test-particle shock, i.e., $`x=[1/(\mu m_\mathrm{H})](kT_\mathrm{s}/V_\mathrm{s}^2)3/16`$. Here $`\mu m_\mathrm{H}`$ is the mean particle mass and $`T_\mathrm{s}`$ is the post-shock temperature. Table 1 shows how $`x`$ drops well below $`3/16`$ for large values of $`\eta _{\mathrm{inj},\mathrm{p}}`$, and indicates that shock velocities derived from a measured $`T_s`$ may be underestimated by a large factor. Independent measurements of $`V_s`$ and $`T_s`$ provide strong constraints on the amount of energy going into relativistic particles. This was shown using Chandra observations of E0102.2-7219 SNR, by Hughes, Rakowski, & Decourchelle (2000), who concluded that the low observed post-shock electron temperature combined with the high observed shock velocity could only be understood if a non-negligible fraction of the shock ram energy went into cosmic rays.

## 4 Interpretation of X-ray thermal spectra

The increase in shocked density and decrease in shocked temperature from nonlinear effects, strongly impact the thermal X-ray emission. Compared to test-particle models, we predict that regions in young SNRs undergoing efficient acceleration will have a lower temperature, more intense thermal X-ray emission (due to the larger emission integral $`n_e^2𝑑V`$ and lower temperature), and a faster ionization timescale. In older remnants, the temperature should drop below X-ray emitting values earlier.

We have coupled our modified hydrodynamics model to the calculation of the non-equilibrium ionization and thermal X-ray emission. The ionization and recombination rates have been taken from Arnaud & Rothenflug () and Arnaud & Raymond (). The X-ray emission code is from Mewe, Gronenschild, & Van den Oord () and Mewe, Lemen, & Van den Oord (). The abundances in the interstellar medium are taken to be solar (Meyer ). We assume full equipartition between the electrons and ions ($`T_e=T_i`$) and neglect any contribution from the superthermal electrons to the ionization or line excitation. For efficient shock acceleration, this may be important and we are investigating these effects (Porquet et al. 2000).

As an indication of how the interpretation of X-ray data depends on acceleration efficiency, we model the X-ray observations of Kepler’s SNR obtained over a broad energy range with ASCA (Kinugasa & Tsunemi 1999) and RXTE (Decourchelle & Petre 1999; Decourchelle et al. 2000). For several values of $`\eta _{\mathrm{inj},\mathrm{p}}`$, we adjust the SNR parameters to fit the observed spectrum, which is dominated by the emission from the shocked ejecta as shown in Figure 3. All values of $`\eta _{\mathrm{inj},\mathrm{p}}`$ yield similar values for the reduced chi-squared ($`4`$), the ionization timescale, the reverse shock temperature $`3`$ keV, the hydrogen column density $`410^{21}`$ cm-2, and oversolar abundances of heavy elements in the shocked ejecta (O, Ne, Mg, Si, S, Ar, Ca, and Fe). The derived values of forward and reverse shock compression ratios ($`r_{\mathrm{FS}}`$ and $`r_{\mathrm{RS}}`$), FS and RS speeds ($`V_{\mathrm{FS}}`$ and $`V_{\mathrm{RS}}`$), ambient density, distance, and the energy-mass relation are, however, very different as shown in Table 2.

To get a temperature high enough to produce the iron K-line and the high energy continuum, all nonlinear fits require larger shock velocities and lower ambient densities. These imply more kinetic energy in the ejecta and a larger inferred distance, $`D`$, than the test-particle case. To obtain a distance consistent with observations (4.8 $`\pm `$ 1.4 kpc, Reynoso and Goss ), a low efficiency at the reverse shock is required. This can be done in three ways: with $`B_0\begin{array}{c}>\hfill \\ \hfill \end{array}25`$$`\mu `$G at the reverse shock large enough to dampen acceleration, with $`B_0\begin{array}{c}<\hfill \\ \hfill \end{array}0.1`$ $`\mu `$G and $`\eta _{\mathrm{inj},\mathrm{pRS}}\begin{array}{c}<\hfill \\ \hfill \end{array}10^4`$ to allow high Mach number, unmodified solutions, or with no reverse shock acceleration ($`\eta _{\mathrm{inj},\mathrm{pRS}}0`$). A low magnetic field is expected from the expansion of the ejecta, but magnetic field amplification could increase it. In all cases, a good fit is obtained for efficient acceleration at the forward shock and test-particle conditions at the reverse shock. Figure 3 and Table 2 (last row) show the case with $`\eta _{\mathrm{inj},\mathrm{pRS}}=0`$.

Spatially resolved spectroscopy is required to confirm both a low efficiency of the acceleration at the reverse shock and our prediction of efficient acceleration at the forward shock, which is expected for standard values of the upstream density and magnetic field. The arbitrary power law component required to fit the broad-band X-ray continuum (Decourchelle & Petre ) is likely related to the forward shock. In future work, this component will be replaced by including self-consistently the nonthermal continuum, and using additional constraints on the acceleration parameters from the radio data.

## 5 Conclusion

Nonlinear effects from efficient particle acceleration strongly influence the X-ray emission in SNRs, and may serve as a tracer of cosmic ray acceleration. Here, we have focused on the consequences of shock acceleration on the thermal X-ray emission and have presented the first nonlinear, non-equilibrium ionization modeling of the X-ray emission from a specific SNR. While our self-similar approach has limitations which must be carefully considered, we believe it contains the essential physics and allows the exploration of a large range of parameters. We find that typical source parameters of young SNRs should result in nonlinear acceleration at the forward shock, but yield near test-particle conditions at the reverse shock if the magnetic field is as low as expected from the ejecta expansion.

When cosmic ray production is efficient, shock heating is linked to particle acceleration and thus to the broad-band photon emission from energetic electrons and ions. The high energy electrons produce both radio and X-ray synchrotron and this will be included in future work.

In addition to a possible nonthermal synchrotron component, regions undergoing strong acceleration should exhibit stronger and lower temperature thermal emission, while regions with little acceleration should be fainter and at a higher temperature.

Also, other complications (e.g., the variation of the ambient magnetic field orientation (Reynolds ), shock interactions with clouds and dense knots, etc.) influence the full modeling of SNRs. The spatially resolved spectra from the Chandra and XMM-Newton satellites should greatly aid the investigation of these and other effects.

###### Acknowledgements.

We thank E. Berezhko, L. Drury, J. Hughes, J.P. Meyer and E. Parizot for fruitful discussions as well as the referee for helpful comments.


format: text

text:
# Dynamic Phase Transition, Universality, and Finite-size Scaling in the Two-dimensional Kinetic Ising Model in an Oscillating Field

## I Introduction

Metastability and hysteresis are widespread phenomena in nature. Ferromagnets are common systems that exhibit these behaviors , but there are also numerous other examples ranging from ferroelectrics  to electrochemical adsorbate layers  to liquid crystals . A simple model for many of these real systems is the kinetic Ising model (in either the spin or the lattice-gas representation). For example, it has been shown to be appropriate for describing magnetization dynamics in highly anisotropic single-domain nanoparticles and uniaxial thin films .

The system response to a single reversal of the “external field” is fairly well understood . In sufficiently large systems below the equilibrium critical temperature, $`T_c`$, the order parameter changes its value through the nucleation and growth of many droplets, inside which it has the equilibrium value consistent with the value of the applied field, as shown in Fig. 1. This is the multi-droplet regime of phase transformation . The well-known Avrami’s law  describes this process of homogeneous nucleation followed by growth quite accurately up to the time when the growing droplets coalesce and the stable phase becomes the majority phase . The intrinsic time scale of the system is given by the metastable lifetime, $`\tau `$, which is defined as the average first-passage time to zero magnetization. It is a measure of the time it takes for the system to escape from the metastable region of the free-energy landscape. In this paper we will use the magnetic language in which the order parameter is the magnetization, $`m`$, and its conjugate field is the external magnetic field, $`H`$. Analogous interpretations, e.g., using the terms polarization and electric field for ferroelectric systems , and coverage and chemical potential for adsorption problems , are straightforward.

It is natural next to ask, “what is the response to an oscillating external field?” The hysteretic behavior in ferromagnets has attracted significant experimental interest, mainly focused on the characteristic behavior of the hysteresis loop and its area. Its dependence on the field amplitude and frequency has been intensively studied and its scaling behavior (power law versus logarithmic) is still under investigation, both experimentally  and theoretically . For the kinetic Ising ferromagnet in two dimensions it has been recently shown  that the true behavior is in fact a crossover, approaching a logarithmic frequency dependence only for extremely low frequencies.

An important aspect of hysteresis in bistable systems, which is the focus of the present paper, is the dynamic competition between the two time scales in the system: the half-period of the external field $`t_{1/2}`$ (proportional to the inverse of the driving frequency) and the metastable lifetime $`\tau `$. For low frequencies, a complete decay of the metastable phase almost always occurs in each half-period, just like it does after a single field reversal. Consequently, the time-dependent magnetization reaches a limit cycle which is symmetric about zero \[Fig. 2(a)\]. For high frequencies, however, the system does not have enough time to switch during one half-period, and the symmetry of the hysteresis loop is broken. The magnetization then reaches an asymmetric limit cycle \[Fig. 2(b)\]. Avrami’s law  is a good approximation when the majority of the droplets do not overlap. Thus, it can be employed to estimate the time-dependent magnetization and the dynamic order parameter (period-averaged magnetization) in the low-frequency (see the Appendix) and in the high-frequency  limits. However, it cannot describe the “critical regime” where $`t_{1/2}`$ becomes comparable to $`\tau `$, and which is dominated by coalescing droplets.

This symmetry breaking between the symmetric and asymmetric limit cycles has been the subject of intensive research over the last decade. It was first observed during numerical integration of a mean-field equation of motion for the magnetization of a ferromagnet in an oscillating field . Since then, it has been observed and studied in numerous Monte Carlo (MC) simulations of kinetic Ising systems , as well as in further mean-field studies . It may also have been experimentally observed in ultrathin films of Co on Cu(001) . The results of these studies suggest that this symmetry breaking corresponds to a genuine continuous non-equilibrium phase transition. For recent reviews see Refs. . Associated with the transition is a divergent time scale (critical slowing down)  and, for spatially extended systems, a divergent correlation length . Estimates for the critical exponents and the universality class of the transition have recently become available  after the successful application of finite-size scaling techniques borrowed from equilibrium critical phenomena .

The purpose of the present paper is to extend preliminary results  and to provide more accurate estimates of the exponents for two-dimensional kinetic Ising systems in a square-wave oscillating field. The use of the square-wave field tests the universality of the dynamic phase transition (DPT) , and it also significantly increases computational speed, compared to the more commonly used sinusoidal field. We further explore the universal aspects of the transition by varying the temperature and field amplitude within the multi-droplet regime, and we study the cross-over to the strong-field regime where the transition disappears. In obtaining our results, we rely on dynamic MC simulations. Computational methods are always helpful, especially when theoretical ideas are largely missing. There are cases, however, when even the use of standard equilibrium techniques, such as finite-size scaling requires some insight and building analogies between equilibrium and non-equilibrium systems . This is the case for our present study. No effective “Hamiltonian” was known before the completion of this work for the dynamic order-parameter (in the coarse-grained sense), from which the long-distance behavior of the model could be derived. This is a typical difficulty when dealing with systems far from equilibrium . Recently, however, a coarse-grained Hamiltonian has been derived  for the dynamic order-parameter, supporting our results for the DPT. Similar to the previous work for sinusoidly oscillating fields , we perform large-scale simulations and finite-size scaling to investigate the universal properties of the DPT.

The remainder of the paper is organized as follows. In Sec. II we define the model and the observables of interest. Section III contains the Monte Carlo results and analyses. Conclusions and outlook are given in Sec. IV.

## II Model and Relevant Observables

To model spatially extended bistable systems in two dimensions, we study a nearest-neighbor kinetic Ising ferromagnet on a $`L`$$`\times `$$`L`$ square lattice with periodic boundary conditions. The model is defined by the Hamiltonian

$$=J\underset{ij}{}s_is_jH(t)\underset{i}{}s_i,$$
(1)
where $`s_i`$$`=`$$`\pm 1`$ is the state of the $`i`$th spin, $`J>0`$ is the ferromagnetic interaction, $`_{ij}`$ runs over all nearest-neighbor pairs, $`_i`$ runs over all $`L^2`$ lattice sites, and $`H(t)`$ is an oscillating, spatially uniform applied field. The magnetization per site,

$$m(t)=\frac{1}{L^2}\underset{i=1}{\overset{L^2}{}}s_i(t),$$
(2)
is the density conjugate to $`H(t)`$. The temperature $`T`$ is fixed below its zero-field critical value $`T_\mathrm{c}`$ ($`J/k_\mathrm{B}T_\mathrm{c}`$$`=`$$`\mathrm{ln}(1+\sqrt{2})/2`$ , where $`k_\mathrm{B}`$ is Boltzmann’s constant), so that the magnetization for $`H`$$`=`$$`0`$ has two degenerate spontaneous equilibrium values, $`\pm m_{\mathrm{sp}}(T)`$. For nonzero fields the equilibrium magnetization has the same sign as $`H`$, while for $`H`$ not too strong, the opposite magnetization direction is metastable and decays slowly towards equilibrium with time, as described in Sec. I.

The dynamic used in this study, as well as in Refs. , is the Glauber single-spin-flip MC algorithm with updates at randomly chosen sites. Note that the random sequential update scheme corresponds to independent Poisson arrivals for the update attempts (discrete events) at each site. Thus, the arrival pattern is strongly asynchronous. The time unit is one MC step per spin (MCSS). Each attempted spin flip from $`s_i`$ to $`s_i`$ is accepted with probability

$$W(s_is_i)=\frac{\mathrm{exp}(\beta \mathrm{\Delta }E_i)}{1+\mathrm{exp}(\beta \mathrm{\Delta }E_i)}.$$
(3)
Here $`\mathrm{\Delta }E_i`$ is the energy change resulting from acceptance, and $`\beta =1/k_\mathrm{B}T`$. For the largest system studied ($`L`$$`=`$$`512`$) we used a scalable massively parallel implementation of the algorithm for this asynchronous dynamics . The parallel discrete-event scheme ensures that the underlying dynamic is not changed (that is, the update attempts are identical, independent Poisson arrivals at each site), while a substantial amount of parallelism is exploited. The parallel implementation  was carried out on a Cray T3E, employing up to 256 processing elements.

The dynamic order parameter is the period-averaged magnetization ,

$$Q=\frac{1}{2t_{1/2}}m(t)𝑑t,$$
(4)
where $`t_{1/2}`$ is the half-period of the oscillating field, and the beginning of the period is chosen at a time when $`H(t)`$ changes sign. Although the phase of the field does not influence the results reported in this paper, the choice made here is convenient in studies of the hysteresis loop-area distributions and consistent with Refs. . Analogously we also define the local order parameter

$$Q_i=\frac{1}{2t_{1/2}}s_i(t)𝑑t,$$
(5)
which is the period-averaged spin at site $`i`$. For slowly varying fields the probability distribution of $`Q`$ is sharply peaked at zero . We shall refer to this as the dynamically disordered phase. It is illustrated by the evolution of the magnetization in Fig. 2(a) and by the $`Q`$$``$$`0`$ time series in Fig. 3. For rapidly oscillating fields the distribution of $`Q`$ becomes bimodal with two sharp peaks near $`\pm m_{\mathrm{sp}}(T)`$, corresponding to the broken symmetry of the hysteresis loops . We shall refer to this as the dynamically ordered phase. It is illustrated in Fig. 2(b) and by the $`Q`$$``$$`m_{\mathrm{sp}}`$$``$$`𝒪(1)`$ time series in Fig. 3. Near the DPT we use finite-size scaling analysis of MC data to estimate the critical exponents that characterize the transition. We also keep track of the normalized period-averaged internal energy (in units of $`J`$) ,

$$E=\frac{1}{2t_{1/2}}\frac{1}{L^2}\underset{ij}{}s_i(t)s_j(t)dt,$$
(6)
since it also exhibits important characteristics of the DPT.

Previous studies of the DPT have used an applied field which varies sinusoidally in time. While sinusoidal or linear saw-tooth fields are the most common in experiments and are necessary to obtain a vanishing loop area in the low-frequency limit , the wave form of the field should not affect universal aspects of the DPT. This should be so because the transition essentially depends on the competition between two time scales: the half-period $`t_{1/2}`$ of the applied field, and the average time it takes the system to leave the metastable region near one of its two degenerate zero-field equilibria when a field of magnitude $`H_0`$ and sign opposite to the magnetization is applied. This metastable lifetime, $`\tau (T,H_0)`$, is estimated as the average first-passage time to zero magnetization. In the present paper we use a square-wave field of amplitude $`H_0`$. This has significant computational advantages over the sinusoidal field variation since we can use two look-up tables to determine the acceptance probabilities: one for $`H=+H_0`$ and one for $`H=H_0`$.

In terms of the dimensionless half-period,

$$\mathrm{\Theta }=t_{1/2}/\tau (T,H_0),$$
(7)
the DPT should occur at a critical value $`\mathrm{\Theta }_\mathrm{c}`$ of order unity. Although $`\mathrm{\Theta }`$ can be changed by varying either $`t_{1/2}`$, $`H_0`$, or $`T`$, in a first approximation we expect $`\mathrm{\Theta }_\mathrm{c}`$ to depend only weakly on $`H_0`$ and $`T`$. This expectation will be confirmed in Sec. IV by simulations carried out at several values of $`H_0`$ and $`T`$ for different system sizes.

In many studies of the DPT the transition has been approached by changing $`H_0`$ or $`T`$ . While this is correct in principle, $`\tau (T,H_0)`$ depends strongly and nonlinearly on its arguments . We therefore prefer changing $`t_{1/2}`$ at constant $`H_0`$ and $`T`$ , as this gives more precise control over the distance from the transition.

We focus on systems which are not only larger than the critical droplet, but also significantly larger than the typical droplet separation . In this regime many supercritical droplets form and contribute to the decay of the metastable phase (the KJMA or Avrami theory for homogeneous nucleation ), as seen in Fig. 1. This is the only regime where the DPT is expected to exist. For small systems one observes subtle finite-size effects, not related to the DPT but rather to the stochastic single-droplet decay mode . In the single-droplet regime, subject to a periodic applied field, the system exhibits stochastic resonance .

## III Simulation results

### A Signs of the dynamic phase transition

We performed extensive simulations and finite-size scaling analysis of the data on square lattices with $`L`$ between $`64`$ and $`512`$ at $`T`$$`=`$$`0.8T_c`$ and $`H_0`$$`=`$$`0.3J`$. We also investigated the universality of the DPT within the multi-droplet regime for various fields and temperatures below the equilibrium critical temperature, using smaller systems with $`L`$ from $`64`$ to $`128`$. Typical runs near the DPT consist of $`2`$$`\times `$$`10^5`$ full periods. For example, at $`T`$$`=`$$`0.8T_c`$ and $`H_0`$$`=`$$`0.3J`$, where the critical half-period of the field is about $`70`$ MCSS, this corresponds to $`2.8`$$`\times `$$`10^7`$ MCSS. Away from the transition point, an order of magnitude shorter runs were sufficient to obtain high-quality statistics.

The system was initialized with all spins up and the square-wave external field started with the half-period in which $`H`$$`=`$$`H_0`$. After some relaxation the system magnetization would reach a limit cycle \[Fig. 2\] (except for thermal fluctuations). In other words, $`Q`$ \[Eq. (4)\] (together with other period-averaged quantities) becomes a stationary stochastic process \[Fig. 3\]. We discarded the first $`1000`$ periods of the time series to exclude transients from the stationary-state averages.

For large half-periods ($`\mathrm{\Theta }\mathrm{\Theta }_\mathrm{c}`$) the magnetization switches every half-period \[Fig. 2(a)\] and $`Q`$$``$$`0`$, while for small half-periods ($`\mathrm{\Theta }\mathrm{\Theta }_\mathrm{c}`$) the magnetization does not have time to switch during a single half-period \[Fig. 2(b)\], resulting in $`|Q|`$$``$$`m_{\mathrm{sp}}`$, as can be seen from the time series in Fig. 3. The transition between the high- and low-frequency regimes is characterized by large fluctuations in $`Q`$ near $`\mathrm{\Theta }_\mathrm{c}`$ \[Fig. 3\].

To illustrate the spatial aspects of the transition we also show configurations of the local order parameter $`\{Q_i\}`$ in Fig. 4. Below $`\mathrm{\Theta }_\mathrm{c}`$ \[Fig. 4(a)\] the majority of spins spend most of their time in the $`+1`$ state, i.e., in the metastable phase during the first half-period, and in the stable equilibrium phase during the second half-period (except for equilibrium fluctuations). Thus, most of the $`Q_i`$$``$$`+1`$. Droplets of $`s_i`$$`=`$$`1`$ that nucleate during the negative half-period and then decay back to $`+1`$ during the positive half-period show up as roughly circular gray spots in the figure. Since the spins near the center of such a droplet become negative first and revert to positive last, these spots appear darkest in the middle. Also, for not too large lattices one occasionally observes the full reversal of an ordered configuration $`\{Q_i\}`$$``$$`\{Q_i\}`$, typical of finite, spatially extended systems undergoing symmetry breaking. Above $`\mathrm{\Theta }_\mathrm{c}`$ \[Fig. 4(c)\] the system follows the field in every half-period (with some phase lag) and $`Q_i`$$``$$`0`$ at all sites $`i`$. Near $`\mathrm{\Theta }_\mathrm{c}`$ \[Fig. 4(b)\] there are large clusters of both $`Q_i`$$``$$`+1`$ and $`Q_i`$$``$$`1`$ separated by “interfaces” where $`Q_i`$$``$$`0`$. These large-scale structures remain reasonably stationary over several periods.

For finite systems in the dynamically ordered phase the probability density of $`Q`$ becomes bimodal. Thus, to capture symmetry breaking, one has to measure the average norm of Q as the order parameter, i.e., $`|Q|`$ . Figure 5(a) shows that this order parameter is of order unity for $`\mathrm{\Theta }<\mathrm{\Theta }_\mathrm{c}`$ and vanishes for $`\mathrm{\Theta }>\mathrm{\Theta }_\mathrm{c}`$, except for finite-size effects. To characterize and quantify this transition in terms of critical exponents we employ the well-known technique of finite-size scaling . The quantity analogous to the susceptibility is the scaled variance of the dynamic order parameter ,

$$X_L^Q=L^2\left(Q^2_L|Q|_L^2\right).$$
(8)
Note that for our system the field conjugate to $`Q`$ and a corresponding fluctuation-dissipation theorem are not known, hence we cannot measure the susceptibility directly. For finite systems $`X_L`$ has a characteristic peak near $`\mathrm{\Theta }_\mathrm{c}`$ \[see Fig. 5(b)\] which increases in height with increasing $`L`$, while no finite-size effects can be observed for $`\mathrm{\Theta }\mathrm{\Theta }_\mathrm{c}`$ and $`\mathrm{\Theta }\mathrm{\Theta }_\mathrm{c}`$. This implies the existence of a divergent length scale, possibly the correlation length which governs the long-distance behavior of the local order-parameter correlations $`Q_iQ_j`$. The location of the maximum in $`X_L^Q`$ also shifts with $`L`$, which gives further important information about the critical exponents.

The normalized stationary time-displaced autocorrelation function of the order parameter,

$$C_L^Q(n)=\frac{Q(j)Q(j+n)Q(j)^2}{Q^2(j)Q(j)^2},$$
(9)
provides further insights into the DPT as the system exhibits critical slowing down \[Fig. 6\]. This can be seen as increasing correlation times with increasing system sizes. In Sec. III.D we provide a quantitative analysis of the correlation times.

We also measured the period-averaged internal energy \[Eq. (6)\] and its fluctuations

$$X_L^E=L^2\left(E^2_LE_L^2\right),$$
(10)
as can be seen in Fig. 7. The peaks of these fluctuations exhibit a slow increase with the system size (compared to the order-parameter fluctuations), as one may anticipate by analogy with the equilibrium heat capacity.

### B Finite-size scaling

Scaling laws and finite-size scaling for equilibrium systems with an a-priori known Hamiltonian can be systematically derived using the concepts of the free energy and the renormalization group . The kinetic Ising model with the explicitly time-dependent Hamiltonian, Eq. (1), is driven far from equilibrium. Although the order-parameter distribution $`P(Q)`$ is stationary, the effective Hamiltonian controlling its fixed-point behavior has not been known until recently  (after the completion of this study). Motivated by the similarity of the finite-size effects shown in Figs. 5-7 to those characteristic of a typical continuous phase transition, we borrow the corresponding scaling assumptions from equilibrium finite-size scaling. For our model the quantity analogous to the reduced temperature in equilibrium systems (i.e., the distance from the infinite-system critical point) is

$$\theta =\frac{|\mathrm{\Theta }\mathrm{\Theta }_\mathrm{c}|}{\mathrm{\Theta }_\mathrm{c}}.$$
(11)
Finite-size scaling theory provides simple scaling relations for the observables for finite systems in the critical regime :

$`|Q|_L`$ $`=`$ $`L^{\beta /\nu }_\pm (\theta L^{1/\nu })`$ (12)
$`X_L^Q`$ $`=`$ $`L^{\gamma /\nu }𝒢_\pm (\theta L^{1/\nu })`$ (13)
$`X_L^E`$ $`=`$ $`c_1\mathrm{ln}\left(L𝒥_\pm (\theta L^{1/\nu })\right),`$ (14)
where $`_\pm `$, $`𝒢_\pm `$, and $`𝒥_\pm `$ are scaling functions and the $`+`$ ($``$) index refers to $`\mathrm{\Theta }>\mathrm{\Theta }_\mathrm{c}`$ ($`\mathrm{\Theta }<\mathrm{\Theta }_\mathrm{c}`$). The logarithmic scaling in $`X_L^E`$ is motivated by the very slow divergence of the scaled period-averaged energy variance \[Fig. 7(b)\]. The above formulation of scaling is explicitly based on the infinite-system critical point $`\mathrm{\Theta }_\mathrm{c}`$, which can be estimated with far greater accuracy than the location of the maximum of the order-parameter fluctuations for the individual finite system sizes. We use the fourth-order cumulant intersection method  to estimate the value of $`\mathrm{\Theta }_c`$ at which the transition occurs in an infinite system. In order to do this, we plot

$$U_L=1\frac{Q^4_L}{3Q^2_L^2}$$
(15)
as a function of $`\mathrm{\Theta }`$ for several system sizes as shown in Fig. 8. For the largest system ($`L`$=$`512`$) the statistical uncertainty in $`U_L`$ was too large to use it to obtain estimates for the crossing. Our estimate for the dimensionless critical half-period, based on the remaining five system sizes, is $`\mathrm{\Theta }_\mathrm{c}`$$`=`$$`0.918`$$`\pm `$$`0.005`$ with a fixed-point value $`U^{}`$$`=`$$`0.611`$$`\pm `$$`0.003`$ for the cumulant \[Fig 8(b)\].

Then at $`\mathrm{\Theta }_c`$ the scaling forms Eqs. (12-14) yield

$`|Q|_L`$ $``$ $`L^{\beta /\nu }`$ (16)
$`X_L^Q`$ $``$ $`L^{\gamma /\nu }`$ (17)
$`X_L^E`$ $``$ $`c_2+c_1\mathrm{ln}(L),`$ (18)
which enable us to estimate the exponent ratios $`\beta /\nu `$ and $`\gamma /\nu `$, and to directly check the postulated logarithmic divergence in the period-averaged energy fluctuations. Plotting $`|Q|_L`$ and $`X_L^Q`$ at $`\mathrm{\Theta }_c`$ and utilizing a weighted linear least-squares fit to the logarithmic data yields $`\beta /\nu `$$`=`$$`0.126\pm 0.005`$ \[Fig. 9(a)\] and $`\gamma /\nu `$$`=`$$`1.74\pm 0.05`$ \[Fig. 9(b)\]. Note that these values are extremely close (within statistical errors) to the corresponding ratios for the equilibrium two-dimensional Ising universality class, $`\beta /\nu `$$`=`$$`1/8`$$`=`$$`0.125`$ and $`\gamma /\nu `$$`=`$$`7/4`$$`=`$$`1.75`$. Further, the straight line in Fig. 9(c) indicates the slow logarithmic divergence of $`X_L^E`$ at the critical point. In addition to the scaling at $`\mathrm{\Theta }_c`$, we also checked the divergences of the peaks of the fluctuations, $`(X_L^Q)_{\mathrm{peak}}`$ and $`(X_L^E)_{\mathrm{peak}}`$, since they asymptotically should follow the same scaling laws, Eqs. (17) and (18), respectively. The measured exponent $`\gamma /\nu `$$`=`$$`1.78\pm 0.05`$ for $`(X_L^Q)_{\mathrm{peak}}`$ and the logarithmic divergence for $`(X_L^E)_{\mathrm{peak}}`$ agree to within the statistical errors with the results obtained at $`\mathrm{\Theta }_c`$, as can be seen in Fig. 9(b) and (c), respectively.

From the finite-system shifting of the transition one can estimate the correlation-length exponent $`\nu `$ by tracking the shift in the location of the maximum in $`X_L^Q`$:

$$|\mathrm{\Theta }_c(L)\mathrm{\Theta }_c|L^{1/\nu },$$
(19)
where $`\mathrm{\Theta }_c(L)`$ is the location of the peak for finite systems. However, the precision of this method for our data is very poor, due to limited resolution in finding the locations of the maxima and consequently the large relative errors in $`|\mathrm{\Theta }_c(L)\mathrm{\Theta }_c|`$. Excluding the smallest (due to strong corrections to scaling) and the largest systems (due to very poor resolution and extremely large statistical error), we obtain $`\nu `$$`=`$$`0.87\pm 0.4`$, but the large error estimate obviously implies rather poor accuracy \[Fig. 10\].

To obtain a more complete picture of how well the scaling relations in Eqs. (12) and (13) hold, we plot $`|Q|_LL^{\beta /\nu }`$ \[Fig. 11(a)\] and $`X_L^QL^{\gamma /\nu }`$ \[Fig. 11(b)\] vs $`\theta L^{1/\nu }`$ . For the exponent ratios we used $`\beta /\nu `$$`=`$$`1/8`$$`=`$$`0.125`$, and $`\gamma /\nu `$$`=`$$`7/4`$$`=`$$`1.75`$, since our estimate for those (within small statistical errors) implied that they take on the equilibrium two-dimensional Ising universal values. Most importantly, we used various values of $`\nu `$ between $`0.5`$ and $`1.2`$ to find the best data collapse as observed visually, since our estimate for this exponent was far from reliable. The “optimal” value obtained this way (by showing scaling plots to group members who did not know the particular values of $`\nu `$ used), and used in Fig. 11(a) and (b), is $`\nu `$$`=`$$`0.95\pm 0.15`$. Full scaling plots using the exact Ising exponents are also shown in Fig. 11(c) and (d), and they result in similarly good data collapse.

### C Order-parameter histograms at criticality

We devote this subsection to analyzing the universal characteristics of the full order-parameter distribution, $`P(Q)`$, at the critical point. This distribution is bimodal for finite systems if observed for sufficiently long times (Fig. 12) . It is more convenient to focus on the distribution of $`|Q|`$, avoiding the effect of the insufficient number of switching events between the two symmetry-broken phases for large systems, which causes the skewness in Fig. 12(b). Figure 13(a) shows the order-parameter distributions $`P_L(|Q|)`$ at the critical point for various system sizes. Finite-size scaling arguments  suggest that at $`\mathrm{\Theta }_\mathrm{c}`$

$$P_L(|Q|)=L^{\beta /\nu }𝒫(L^{\beta /\nu }|Q|).$$
(20)
Thus, the scaled distributions, $`L^{\beta /\nu }P_L(|Q|)`$ vs $`x`$$`=`$$`|Q|L^{\beta /\nu }`$, should fall on the same curve $`𝒫(x)`$ for different system sizes. Again, we used $`\beta /\nu =1/8`$. The quality of the data collapse is quite impressive \[Fig. 13(b)\], with deviations mainly observed for the smallest $`L`$ and the largest values of $`|Q|`$ \[Fig. 13(c)\], possibly as a result of corrections to scaling.

What we find somewhat surprising, is that the distribution appears to be identical (except for stronger corrections to scaling at the DPT) to that of the equilibrium two-dimensional Ising model on a square lattice with periodic boundary conditions at criticality, without a need for any additional scaling parameters. We checked this by performing standard equilibrium two-dimensional Ising simulations with Glauber dynamics and system sizes ranging from $`L`$$`=`$$`64`$ to $`L`$$`=`$$`128`$, and also by comparing our scaled DPT order-parameter histograms to the high-precision two-dimensional equilibrium Ising MC data of Ref.  (Fig. 13). We had expected the shapes of the distributions to be identical for the DPT and equilibrium Ising model, as a consequence of the identical values for the cumulant fixed-point value $`U^{}`$. However, it is not obvious to us why the microscopic length scales in the DPT and the equilibrium Ising model also appear to be identical, as evidenced by the absence of the need for an additional scaling, $`L`$$``$$`L/a`$ ($`a`$ being the microscopic length scale in the DPT).

### D Critical slowing down

Computing the stationary autocorrelation function given by Eq. (9), we already pointed out that at $`\mathrm{\Theta }_\mathrm{c}`$ the correlation time increases fast with system size \[Fig. 6\]. Correlation times are typically extracted from an exponential decay as

$$C_L^Q(n)e^{n/\tau _L^Q},$$
(21)
and they are expected to be finite for finite systems. The correlation time $`\tau _L^Q`$ is also well defined in the $`L`$$``$$`\mathrm{}`$ limit away from the transition. However, it diverges with $`L`$ at the transition point as

$$\tau _L^QL^z,$$
(22)
where $`z`$ is the dynamical critical exponent. For not too late times we had reasonable statistics including the larger systems (up to $`L`$$`=`$$`256`$) to fit the usual exponential decay \[Fig. 14(a)\]. Then plotting the correlation times $`\tau _L^Q`$ vs $`L`$ yields the dynamic exponent $`z=1.91\pm 0.15`$, as shown in Fig. 14(b). This value is within two standard deviations of most estimates for the dynamic exponent of the two-dimensional equilibrium Ising model with local dynamics .

### E Universality for various temperatures and fields and cross-over to the strong-field regime

The underlying ingredient for the spatially extended bistable systems exhibiting a DPT is the local metastability (and the corresponding characteristic time spent in the metastable “free-energy well”) in the presence of an external field. This, in turn, provides a competition between time scales if the system is driven by a periodic field. Based on this, we expect that sufficiently large systems (in which many droplets contribute to the decay of the metastable phase) exhibit the DPT at a half-period $`t_{1/2}`$ comparable to the metastable lifetime $`\tau (T,H_0)`$. In other words, we expect the critical dimensionless half-period to be of order one, $`\mathrm{\Theta }_\mathrm{c}𝒪(1)`$.

To test this expectation, we performed simulations at $`T`$$`=`$$`0.8T_c`$ for field amplitudes ranging from $`0.3J`$ to infinity with system sizes $`L`$$`=`$$`64`$, $`90`$, and $`128`$. \[$`H_0`$$`=`$$`\mathrm{}`$ corresponds to the Glauber spin-flip probabilities Eq. (3) being equal to $`0`$ (1) depending on whether the spin is parallel (anti parallel) to the external field, with no influence from the configuration of the neighboring spins.\] We further performed runs at $`T`$$`=`$$`0.9T_c`$, $`T`$$`=`$$`0.6T_c`$, and $`T`$$`=`$$`0.5T_c`$ for various field amplitudes and system sizes $`L=64`$ and $`90`$. The typical run length was $`2`$$`\times `$$`10^4`$ periods. The purpose of these runs was to explore the universal nature of the DPT in the multi-droplet regime, and the crossover to the strong-field regime where the DPT should disappear.

In the strong-field regime the droplet picture breaks down since the individual spins are decoupled. Thus, the metastable phase no longer exists, and the decay of the phase having opposite sign to the external field approaches a simple exponential form (which becomes exact in the $`H_0`$$``$$`\mathrm{}`$ limit). In the Appendix we show that under these conditions the system magnetization always relaxes to a symmetric limit cycle with $`Q`$$`=`$$`0`$ for all frequencies, thus, no DPT can exist.

Figures 15(a), (b), and (c) show the order parameter vs the dimensionless half-period for $`L=64`$ and a range of field amplitudes at $`T`$$`=`$$`0.8T_c`$, $`T`$$`=`$$`0.6T_c`$, and $`T`$$`=`$$`0.5T_c`$, respectively. The typical order-parameter profile where the system exhibits the DPT prevails up to some temperature dependent cross-over field amplitude $`H_\times (T)`$ \[filled symbols in Fig. 15\]. For $`H_0>H_\times (T)`$ the underlying decay mechanism belongs to the strong-field regime, and correspondingly the DPT disappears as expected.

One can trace this crossover from the multi-droplet to the strong-field regime by plotting $`\mathrm{\Theta }_\mathrm{c}`$ vs $`H_0/J`$ \[Fig. 16(a)\] and vs $`\tau `$ \[Fig. 16(b)\] for fixed temperatures. Here again we employed the fourth-order cumulant intersection method using $`L`$$`=`$$`64`$ and $`L`$$`=`$$`90`$ to identify the infinite-system transition point, $`\mathrm{\Theta }_\mathrm{c}`$. The crossover to the strong-field regime is indicated by the drop in $`\mathrm{\Theta }_\mathrm{c}`$ for large fields (small lifetimes).

The DPT in the multidroplet regime ($`H_0`$$`<`$$`H_\times (T)`$) appears to be universal, as can be seen from the scaled order-parameter plot in Fig. 17. Note that both graphs contain $`28`$ DPT data sets: three different temperatures, with four different field amplitudes for each, and at least two different system sizes ($`L=64,90`$, and $`128`$ at $`T`$$`=`$$`0.8T_c`$, and $`L=64`$ and $`90`$ at $`T`$$`=`$$`0.6T_c`$ and $`T`$$`=`$$`0.5T_c`$) for all these parameters! While the slopes in the asymptotic scaling regime appear to be the same, the small parallel shift may be the result of the non-universal critical amplitudes at different temperatures and fields, or simply our inaccuracy in determining $`\mathrm{\Theta }_\mathrm{c}(T,H_0)`$, due to the relatively short runs.

On the other hand, as expected, the system shows no singularity and $`|Q|`$$`\stackrel{L\mathrm{}}{}`$$`0`$ for any non-zero frequency in the strong-field regime, $`H_0`$$`>`$$`H_\times (T)`$ (see the Appendix). Here, the finite-size effects simply reflect the central-limit theorem, i.e., $`Q^2`$$``$$`𝒪(1/L^2)`$, or $`|Q|`$$``$$`𝒪(1/L)`$. Figures 18(a) and (b) illustrate this at $`T`$$`=`$$`0.8T_c`$, $`H_0`$$`=`$$`1.5J`$ and at $`T`$$`=`$$`0.8T_c`$, $`H_0`$$`=`$$`\mathrm{}`$, respectively. One also expects that a similar behavior prevails at $`T`$$`>`$$`T_\mathrm{c}`$ for any field amplitude, i.e., $`H_\times (T)`$ vanishes for $`T`$$`>`$$`T_\mathrm{c}`$. We checked this for $`T`$$`=`$$`1.1T_c`$ and $`H_0`$$`=`$$`0.05J`$, and the results confirm the $`|Q|`$$``$$`𝒪(1/L)`$ scaling of the order parameter for all frequencies, showing none of the characteristic finite-size effects of a DPT \[Fig. 18(c)\]. This result is in distinct disagreement with older work such as Refs. , which report observations of the DPT at temperatures considerably above $`T_\mathrm{c}`$. For smaller system sizes, however, one may observe nontrivial resonance .

Due to the expected complications of a divergent equilibrium correlation length, we did not perform simulations at $`T_\mathrm{c}`$. However, we conjecture that $`H_\times (T)`$ vanishes at $`T_\mathrm{c}`$. We expect this conjecture to be extremely difficult to prove or disprove numerically, due to very large and possibly complicated finite-size effects.

## IV Conclusions and Outlook

In this paper we have studied the hysteretic response of a spatially extended bistable system exhibiting a DPT. Our model system is the two-dimensional kinetic Ising ferromagnet below its equilibrium critical temperature subject to a periodic square-wave applied field. The results indicate that for field amplitudes and temperatures such that the metastable phase decays via the multi-droplet mechanism, the system undergoes a continuous dynamic phase transition when the half-period of the field, $`t_{1/2}`$, is comparable to the metastable lifetime, $`\tau (T,H_0)`$. Thus, the critical value $`\mathrm{\Theta }_\mathrm{c}`$ of the dimensionless half-period defined in Eq. (7) is of order unity. As $`\mathrm{\Theta }`$ is increased beyond $`\mathrm{\Theta }_\mathrm{c}`$, the order parameter $`|Q|`$ (the expectation value of the norm of the period-averaged magnetization) vanishes \[Fig. 5(a)\], displaying singular behavior at the critical point, as shown in Fig. 5(b) and (c).

The characteristic finite-size effects in the order parameter and its fluctuations indicate that there is a divergent correlation length associated with the transition. We used standard finite-size scaling techniques adopted from the theory of equilibrium phase transitions. We estimated $`\mathrm{\Theta }_\mathrm{c}`$ and the critical exponents $`\beta `$, $`\gamma `$, and $`\nu `$ from relatively high precision data for system sizes between $`L=64`$ and $`512`$ at $`T`$$`=`$$`0.8T_\mathrm{c}`$ and $`H_0`$$`=`$$`0.3J`$. Our best estimates are $`\beta /\nu `$$`=`$$`0.126\pm 0.005`$, $`\gamma /\nu `$$`=`$$`1.74\pm 0.05`$, and $`\nu `$$`=`$$`0.95\pm 0.15`$. These values agree within statistical errors with those previously obtained with a sinusoidally oscillating field , providing strong evidence that the shape of the field oscillation does not affect the universal aspects of the DPT. Observing the stationary autocorrelation function we also saw that at the transition point the system exhibits critical slowing down governed by the dynamic exponent $`z`$$`=`$$`1.91\pm 0.15`$. This is also very close to the corresponding exponent $`z`$$`=`$$`2.12(5)`$, measured in standard two-dimensional Ising simulations with local dynamics . Our best values for the exponent ratios $`\beta /\nu `$ and $`\gamma /\nu `$ are given with relatively high confidence, while for $`\nu `$ it is rather poor. In this sense tracking down the exponent $`\nu `$ and obtaining an accurate estimate for it remains elusive. Note, however, that we could only rely on the standard (single spin-flip) MC algorithm, since we had to preserve the underlying dynamics. Using more sophisticated algorithms to avoid critical slowing down as seen in Figs. 6 and 14, would require an underlying “Hamiltonian” for the corresponding local order parameter $`\{Q_i\}`$, which is not yet known. While in a coarse-grained/universal sense a $`\varphi ^4`$ Hamiltonian is supported by our data, it does not point to any one particular microscopic Hamiltonian for the microscopic order parameter $`\{Q_i\}`$. However, a $`\varphi ^4`$ coarse-grained Hamiltonian for $`\{Q_i\}`$ has been recently derived starting from the time-dependent Ginzburg-Landau equation for the magnetization.

Of the known universality classes, our exponent estimates for the DPT are closest (and within the statistical errors) to those of the the two-dimensional equilibrium Ising model: $`\beta /\nu `$$`=`$$`1/8`$$`=`$$`0.125`$, $`\gamma /\nu `$$`=`$$`7/4`$$`=`$$`1.75`$, and $`\nu `$$`=`$$`1`$. Consequently, our measured exponent ratios satisfy the hyperscaling relation

$$2(\beta /\nu )+\gamma /\nu =1.99\pm 0.05d,$$
(23)
where $`d`$$`=`$$`2`$ is the spatial dimension . Further, the fixed-point value of the fourth-order cumulant, $`U^{}`$$`=`$$`0.611`$$`\pm `$$`0.003`$, is also extremely close to that of the Ising model, $`U^{}`$$`=`$$`0.6106901(5)`$ . These findings provide conclusive evidence that the DPT indeed corresponds to a non-trivial fixed point. We tested the full data collapse for the scaled order parameter and its variance, as shown in Fig. 11, and it confirmed the existence of the universal scaling functions given by Eqs. (12) and (13). Also, at the critical frequency, the order-parameter distributions follow finite-size scaling predictions, Eq. (20), as shown in Fig. 13. More surprisingly, the critical DPT order-parameter distributions fall on that of the two-dimensional equilibrium Ising model at the critical temperature (except for stronger corrections to scaling), without any additional fitting of the underlying microscopic length scale.

While our finite-size scaling data clearly indicate the existence of a divergent length scale, we did not measure the correlation length for the local order parameter directly. Future studies may include extracting the correlation length, $`\xi _Q`$, from the $`Q_iQ_j`$ correlations (or from the corresponding structure factor). This approach would also provide another way to measure the exponent $`\nu `$ by plotting $`\xi _Q`$ vs $`\theta `$ in the critical regime for large systems, and assuming $`\xi _Q`$$``$$`\theta ^\nu `$.

We also studied the universal aspects of the DPT at other temperatures and fields. Shorter runs at $`T`$$`=`$$`0.8T_\mathrm{c}`$, $`T`$$`=`$$`0.6T_\mathrm{c}`$, and $`T`$$`=`$$`0.5T_\mathrm{c}`$ for field amplitudes $`H_0`$$`<`$$`H_\times (T)`$ also confirmed scaling and the universality of the DPT \[Fig. 17\]. The condition for the field amplitude implies that the system only exhibits a DPT in the multi-droplet regime. For $`H_0`$$`>`$$`H_\times (T)`$ strong-field behavior governs the decay of the magnetization, and the DPT disappears, as indicated by Fig. 16 and Fig. 18(a),(b). We also found that the high-temperature phase is qualitatively similar to the strong-field regime in that there is no sign of a DPT for $`T`$$`>`$$`T_\mathrm{c}`$ \[Fig. 18(c)\].

One may ask how general the phenomenon of a DPT is in spatially extended bistable systems, subject to a periodic applied “field” which drives the system between its metastable and stable “wells.” It is possible that having up-down symmetry of the period-averaged “magnetization” is sufficient for possessing a Hamiltonian at the coarse-grained level, even if the system is driven far away from equilibrium and is microscopically irreversible . Future research can address this question by studying other systems (not necessarily ferromagnets) that exhibit hysteresis.

###### Acknowledgements.

We thank S. W. Sides, S. J. Mitchell, and G. Brown for stimulating discussions. We would like to thank W. Janke for providing us with data from Ref.  for comparison of the critical order-parameter distributions. We acknowledge support by the US Department of Energy through the former Supercomputer Computations Research Institute, by the Center for Materials Research and Technology at Florida State University, and by the US National Science Foundation through Grants No. DMR-9634873, DMR-9871455, and DMR-9981815. This research also used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098.

## Low-frequency and Strong-field Mean-field Approximation

For simplicity, in the following we assume that the magnetization decays from $`m`$$`=`$$`+1`$ to $`m`$$`=`$$`1`$ after a single field reversal ($`H_0`$$``$$`H_0`$). This is a good approximation below the equilibrium critical temperature ($`m_{\mathrm{sp}}`$$``$$`1`$), and for any temperature when $`H_0`$$``$$`\mathrm{}`$. Further, we assume that the volume fraction of meta- or unstable spins follows a simple monotonic decay, $`\widehat{\phi }(t)`$.

In terms of the volume fraction of positive spins, $`\varphi (t)`$, the magnetization can be written as

$$m(t)=2\varphi (t)1.$$
(24)
Subject to a square-wave field,

$$H(t)=\{\begin{array}{cc}H_0\hfill & 0t<t_{1/2}\hfill \\ +H_0\hfill & t_{1/2}t<2t_{1/2}\hfill \end{array},$$
(25)
in the first (second) half-period the volume fraction of the positive (negative) spins decays according to $`\widehat{\phi }(t)`$. Thus, in each period (measuring time $`t`$ from the beginning of the period)

$$\varphi (t)\{\begin{array}{cc}\varphi (0)\widehat{\phi }(t)\hfill & 0tt_{1/2}\hfill \\ 1\left[1\varphi (t_{1/2})\right]\widehat{\phi }(tt_{1/2})\hfill & t_{1/2}t2t_{1/2}\hfill \end{array}.$$
(26)
Using this approximation, one directly obtains a linear mapping

$$\varphi _{n+1}=1[1\varphi _n\widehat{\phi }(t_{1/2})]\widehat{\phi }(t_{1/2}),$$
(27)
where $`\varphi _n`$$``$$`\varphi (2nt_{1/2})`$ is the volume fraction of the positive spins at the beginning of the $`n`$th period, $`n=0,1,2,\mathrm{}`$. The stationary value of this quantity is

$$\varphi ^{}=\underset{n\mathrm{}}{lim}\varphi _n=\frac{1}{1+\widehat{\phi }(t_{1/2})}.$$
(28)
Consequently, the magnetization reaches a stationary limit cycle

$$m(t)\{\begin{array}{cc}\frac{2}{1+\widehat{\phi }(t_{1/2})}\widehat{\phi }(t)1\hfill & 0tt_{1/2}\hfill \\ 1\frac{2}{1+\widehat{\phi }(t_{1/2})}\widehat{\phi }(tt_{1/2})\hfill & t_{1/2}t2t_{1/2}\hfill \end{array}.$$
(29)
In this limit cycle the magnetization oscillates about zero,

$$m(0)=m(t_{1/2})=\frac{1\widehat{\phi }(t_{1/2})}{1+\widehat{\phi }(t_{1/2})}$$
(30)
and the symmetry of the magnetization, $`m(t\pm t_{1/2})`$$`=`$$`m(t)`$ implies

$$Q=\frac{1}{2t_{1/2}}m(t)𝑑t=0.$$
(31)
This corresponds to the symmetric (dynamically disordered) phase. Note that this symmetric phase is always reached when $`\widehat{\phi }`$ decreases monotonically from unity at $`t`$$`=`$$`0`$ to zero as $`t`$$``$$`\mathrm{}`$

In the multi-droplet regime, the volume fraction of the metastable phase decays according to Avrami’s law . In the low-frequency limit, $`t_{1/2}`$$``$$`\tau `$ ($`\mathrm{\Theta }`$$``$$`1`$), each half-period almost always contains a complete metastable decay \[Fig. 2(a)\]. Avrami’s law for the metastable volume fraction in each half-period can then be directly applied using $`\widehat{\phi }(t)`$$`=`$$`\phi _{\mathrm{ms}}(t)`$$``$$`e^{(\mathrm{ln}2)t^3/\tau ^3}`$ . This functional form for $`\phi _{\mathrm{ms}}(t)`$ breaks down when $`t_{1/2}`$ becomes comparable to $`\tau `$ ($`\mathrm{\Theta }`$$``$$`1`$), thus this simple mean-field approximation cannot predict any instability related to the DPT.

In the $`H_0`$$``$$`\mathrm{}`$ limit the individual spins become decoupled. Then one can obtain $`\widehat{\phi }(t)`$$`=`$$`e^{(\mathrm{ln}2)t/\tau }`$ which is exact for all frequencies. Further, this exponential decay is a good approximation everywhere in the strong-field regime, thus, no DPT can exist there.


text:
# Nuclear Aspects of Nucleosynthesis in Massive Stars

## 1 Introduction

Nuclear reactions play a major role not only in the nucleosynthetic processes determining the elemental abundances in the solar system and the Galaxy but also for determining structure and final fate of a star. Massive stars ($`>8M_{}`$) experience a number of burning phases before they explode as type II supernovae after the collapse of the Fe core. Important nuclear reactions in the late burning stages and in the explosion proceed on isotopes experimentally not sufficiently well investigated or on unstable nuclei which cannot be studied in the laboratory. Thus, astrophysics requires a sound theoretical understanding of nuclear reactions and tests our knowledge at the extremes.

Numerous studies have been devoted to the evolution of such stars and their nucleosynthetic yields. However, our knowledge of both the input data and the physical processes affecting the development of these objects has improved dramatically in recent years. Thus, it became worthwhile to attempt to improve on and considerably extend the previous investigations on pre– and post–collapse evolution and nucleosynthesis. Here we present first results for a 15 $`M_{}`$ stellar model with improved stellar and nuclear physics. In this report we mainly concentrate on a few of the nuclear physics issues involved, a more extended report including all details of the simulation will be published elsewhere. Below, we discuss the prediction of nuclear rates in the statistical model and, specifically, the treatment of $`\alpha `$-particle capture on isospin conjugated targets. The importance of obtaining more information on $`\alpha `$+nucleus potentials is emphasized separately.

## 2 Nuclear Reactions

The nuclear reaction network during the explosive phase (as the most extreme case in our calculation) contains about 2350 isotopes and is shown in Fig. 1. It is evident that there are many isotopes far off stability included for which experimental information is scarce. This is even more true in more exotic scenarios such as the r– and rp–processes which involve isotopes close to the neutron– and proton–dripline, respectively. (In our calculations we do not follow the proposed r–process in the $`\nu `$–wind emanating from the proto–neutron star shortly after the collapse of the Fe core.)

Important are weak reactions and nuclear reactions with nucleons and $`\alpha `$ particles. Decay data is also available further off stability whereas nuclear reaction cross sections involving nucleons and light ions are practically known only for stable nuclei. The majority of the latter reactions can be described in the framework of the statistical model (Hauser-Feshbach theory) which describes the reaction proceeding via the formation of a compound nucleus and averages over resonances. Many nuclear properties enter the computation of the HF cross sections: mass differences (separation energies), optical potentials, Giant Dipole Resonance widths, level densities. The resulting transmission coefficients can be modified due to pre-equilibrium effects which are included in width fluctuation corrections (see also a previous paper and references therein) and by isospin effects. It is in the description of the nuclear properties where the various HF models differ. In astrophysical applications usually different aspects are emphasized than in pure nuclear physics investigations. Many of the latter in this long and well established field were focused on specific reactions, where all or most ”ingredients” were deduced from experiments. As long as the reaction mechanism is identified properly, this will produce highly accurate cross sections. For the majority of nuclei in astrophysical applications such information is not available. The real challenge is thus not the application of well-established models, but rather to provide all the necessary ingredients in as reliable a way as possible, also for nuclei where no such information is available.

### 2.1 Statistical Model Rates

As the basis for the creation of our reaction rate set we used statistical model calculations obtained with the NON-SMOKER code. A library of theoretical reaction rates calculated with this code and fitted to an analytical function — ready to be incorporated into stellar model codes — was published recently. It includes rates for all possible targets from proton– to neutron–dripline and between Ne and Bi, thus being the most extensive published library of theoretical reaction rates to date. It also offers rate sets for a number of mass models which are suited for different purposes. For the network described here we utilized the rates based on the FRDM set.

### 2.2 $`\alpha `$ Particles: Isospin Effects

The consideration of isospin effects has two major effects on statistical cross sections in astrophysics: the suppression of $`\gamma `$ widths for reactions involving self-conjugate nuclei and the suppression of the neutron emission in proton-induced reactions. Here, we only discuss the former. In the case of ($`\alpha `$,$`\gamma `$) reactions on targets with $`N=Z`$, the cross sections will be heavily suppressed because $`T=1`$ states cannot be populated due to isospin conservation. A suppression will also be found for capture reactions leading into self-conjugate nuclei, although somewhat less pronounced because $`T=1`$ states can be populated according to the isospin coupling coefficients. In previous reaction rate calculations the suppression of the $`\gamma `$–widths was treated completely phenomenologically by employing arbitrary and mass-independent suppression factors. In the NON-SMOKER code, the appropriate $`\gamma `$ widths are automatically obtained, by explicitly accounting for $`T^<`$ and $`T^>`$ states. The astrophysical importance of $`\alpha `$ capture on target nuclei with $`N=Z`$ is manifold. In the Ne- and O-burning phase of massive stars, alpha capture reaction sequences are initiated at 24Mg and 28Si, respectively, and determine the abundance distribution prior to the Si-burning phase. Nucleosynthesis in explosive Ne and explosive O burning in type II supernovae depend on reaction rates for $`\alpha `$ capture on 20Ne to 36Ar. An $`\alpha `$ capture chain on such self-conjugate nuclei actually determines the production of 44Ti , which contributes to the light curve by its $`\beta `$ decay to 44Ca via 44Sc.

### 2.3 $`\alpha `$ Particles: Optical $`\alpha `$+Nucleus Potentials

A further complication in the treatment of reactions on intermediate and light nuclei involving $`\alpha `$ particles is the limited success in defining an appropriate optical potential, especially for the low energies typical for astrophysical environments. Early astrophysical studies (e.g.) made use of simplified equivalent square well potentials and the black nucleus approximation. It is equivalent to a fully absorptive potential, once a particle has entered the potential well and therefore does not permit resonance effects. This leads to deviations from experimental data at low energies, especially in mass regions where broad resonances in the continuum can be populated. An additional effect, which is only pronounced for $`\alpha `$ particles, is that absorption in the Coulomb barrier is neglected in this approach. Improved calculations have to employ appropriate global optical potentials which make use of imaginary parts describing the absorption. In the case of $`\alpha `$-nucleus potentials, there were only very few global parametrizations attempted at astrophysical energies, also due to the scarcity of experimental data in the energy region of interest. The high Coloumb barrier makes a direct experimental approach very difficult at low energies. Current astrophysical calculations mostly employ a phenomenological Saxon–Woods potential based on extensive data. Future improved $`\alpha `$ potentials have to take into account the mass- and energy-dependence of the potential. Extended investigations of $`\alpha `$ scattering data have shown that the data can best be described with folding potentials. Few attempts have been made to construct such an improved global potential. Nevertheless, the postulation of such an optical potential close to or below the Coulomb barrier remains one of the major challenges. More experimental data is clearly needed.

## 3 Further Inputs

Further nuclear input were updated rates from experimental cross sections for light as well as heavy nuclei and updated beta-decay rates. New predictions of weak rates were also included. In respect to earlier simulations we also used updated neutrino loss rates and opacity tables (OPAL95), and consider mass loss due to stellar wind. For further details, refer to our other papers.

## 4 Results and Summary

For the first time, we studied consistently the production of all isotopes up to Bi during the pre-supernova evolution and the type II supernova explosion of a massive star. Exemplary for our results, the production factors of a 15 $`M_{}`$ star are shown in Fig. 2. The revised weak rates introduce an important change mainly during core silicon burning and thereafter leading to an increase of the central $`Y_e`$ and smaller Fe core masses at the onset of core collapse. In addition to the well-known strong dependence of the stellar structure on the 12C($`\alpha `$,$`\gamma `$)$`{}_{}{}^{16}O`$ rate, we also found the ($`\alpha `$,n)–($`\alpha `$,$`\gamma `$) branching on 22Ne to be an important candidate for further laboratory study. It sensitively determines the strength of the s–process in the SN models.

Summarizing, the progress in predictions of nuclear reactions has made it possible to consistently study the nucleosynthesis in a type II supernova model over a wide range of nuclear masses. The new investigations also underline the importance of new experimental and theoretical studies of specific nuclear properties.

## Acknowledgments

This research was supported, in part, by DOE (W-7405-ENG-48), NSF (AST 97-31569, INT-9726315), and the Alexander von Humboldt Foundation (FLF-1065004). T. R. acknowledges support by a PROFIL professorship from the Swiss NSF (grant 2124-055832.98).


text:
# Adsorbed 3d transition metal atoms and clusters on Au(111): Signatures derived from one-electron calculations.

## Abstract

The spectroscopic characteristics of systems with adsorbed $`d`$ impurities on noble metal surfaces should depend on the number and geometric arrangement of the adsorbed atoms and also on their $`d`$ band filling. Recent experiments using scanning tunneling microscopy have probed the electronic structure of all $`3d`$ transition metal impurities and also of Co dimers adsorbed on Au(111), providing a rich variety of results. In this contribution we correlate those experimental results with ab-initio calculations and try to establish necessary conditions for observing a Kondo resonance when using the single impurity Anderson model. We find that the relevant orbitals at the STM tip position, when it is on top of an impurity, are the $`d`$ orbitals with $`m=0`$ and that the energy of these levels with respect to the Fermi energy determines the possibility of observing a spectroscopic feature due to the impurity.

The electronic structure of $`3d`$ transition-metal impurities on the (111) surface of copper and gold has been studied recently with scanning tunneling microscopy . The first measurements, performed only for Co impurities, gave rise to a direct spectroscopic observation of the Kondo resonance for an isolated impurity . Further experiments with Co dimers  and with other $`3d`$ transition metal atoms  show that there is a different behavior of $`dI/dV`$ for different impurity $`d`$ filling, thus varying significantly across the $`3d`$ series in the periodic table. These results provide a rich phenomenology that we believe can be correlated with results of LSDA calculations.

In a previous paper  we calculated the electronic densities of states of a periodic system that simulates the experimental situation described above. It consisted of repeated slabs of five layers of Au(111) with Co atoms deposited on both sides of them forming a dilute adsorbed layer. For this purpose we used the FP-LAPW method and the LSDA approximation with the WIEN97 code . The main result of that paper was that the minority spin contribution to the density of states due to the $`3d`$ orbitals of Co with $`m=0`$ is a very narrow band located precisely at $`E_F`$, which could be related to the observed experimental feature.

In the present work we have extended the calculation to other transition metals and also to dimers. The unit cell is the same for all the impurity cases, it has five layers of Au with three atoms in each layer. One impurity atom is located on each side of the slab at one of the hollow sites of the Au(111) slab structure forming a $`\sqrt{3}\times \sqrt{3}`$ adsorbed layer. The slabs are separated by enough empty space so as to simulate non interacting surfaces. The distance between impurity atoms on the same plane is 5 Å, and their distance to the surfaces is relaxed so as to minimize the energy.

Magnetism can be quantified in these calculations by giving the magnetic moment inside the muffin tin spheres or the $`d`$ band splitting for the different impurities. We show these results in Table I, they agree with previous calculations for similar adsorbed systems . The given splittings are those corresponding to the separation between majority and minority centers of the $`m=0`$ narrow band, that we assume is the more relevant one. Both magnitudes, $`\mu `$ and $`E_{up}E_{dn}`$, are seen to be roughly proportional. In order to simplify the analysis and to compare with results from model Hamiltonian calculations we have made the working assumption that it is enough to consider just the single $`3d`$ level with $`m=0`$ in order to correlate the calculated results with the STM experimental ones. These experiments were performed in a small energy range around the Fermy energy of Au, using a bias voltage of less than 0.1 eV. Our calculations indicate that only the majority band of V and the minority band of Co fall within that range. One would then expect similar experimental results for both impurities, but this is not the case .

Local densities of states at the impurity sites, of $`d`$ symmetry and for $`m=0`$ in particular, are shown in Figs. 1 to 6. Also shown in that figure are the charge densities in the vacuum region, plotted along the plane perpendicular to the slab at the impurity site, corresponding to two small energy windows close to $`E_F`$. According to the usual theory of STM spectroscopy  the tunneling current is directly related to the charge density at the tip position. We see from Figs. 1 to 6 that the magnitude of this charge density along the normal to the surface is in all cases significantly larger for the energy window that contains the $`3d`$ orbital with $`m=0`$ symmetry, namely $`w1`$. The largest hybridization with the conduction electron sea therefore takes place at the energies where the $`m=0`$ symmetry provides the largest contribution. Thus, tunneling conductance should be essentially affected by the presence of these symmetries on the impurity atoms if they happen to fall in the narrow energy window where the experiments are made. This effect should be observed both for minority or majority spins, if the energy conditions were given.

When analyzing each case in more detail we see that for Ti the majority spin electron density in window $`w1`$ is about one order of magnitude larger than the minority contribution at the tip site (4 to 5 Å, from the surface) and also larger than the contributions of both majority and minority character in the other window ($`w2`$). This agrees very well with the observed increase in the STM signal when the bias voltage sweeps energies from below to above the Fermi level, going across this last one. We believe this reflects the increased density of states just above the Fermi level, when electrons start to flow from the tip into the empty states added by the impurity. Because of the experimental width of this feature it is possibly not a manifestation of the Kondo effect. The second largest value for the ratio of charge density in window $`w1`$ versus the one in window $`w2`$ is obtained for Ni, and may be interpreted in a similar way as that of Ti, although the magnitudes are much smaller. This system could not show a Kondo effect as it is non-magnetic. For Cr and Fe the relevant states are not close to $`E_f`$ so that their influence due to hybridization should not be visible. The cases of Co and V, where we find a narrow band peaked at $`E_f`$, present very similar values of the ratio of charge densities in both windows,being this ratio of the order of 2. Besides this, the absolute charge density values are also very small. We shall analyze below in which cases one could expect a signature of the Kondo peak.

An important one-electron feature, not considered in our calculations, is the interaction of the adsorbed atoms with the (111) surface state at the $`\mathrm{\Gamma }`$ point of noble metals. This state is responsible for the appearance of standing wave patterns in quantum corrals, and also possibly of recently observed mirages . Our representation of the surface by slabs that are only five layers thick would give a poor approximation of the surface state but the periodicity of the adsorbed atoms, that are not too far away from each other, probably destroys the surface state completely. As a first approximation we assume that this interaction will not depend on $`d`$ band filling.

To study the importance of many body effects not considered in the LSDA calculations the single impurity Anderson Hamiltonian (SIAM) is one of the usual approaches . It is known to be equivalent, for most sets of parameters, to the Kondo model Hamiltonian, through the Schrieffer-Wolf transformation . It consists of a conduction band, that provides the Fermi energy or zero energy, and a single magnetic impurity level at $`E_0`$. This level interacts with the band, thus acquiring a width $`\delta `$, and electrons with opposite spins occupying that level feel a Coulomb repulsion $`U`$. The mean field solution of this model Hamiltonian results in two narrow bands, centered at $`E_0`$ and $`E_0+U`$, which could be parametrized for the different adsorbed systems using the results of our LSDA calculations, for example setting for $`E_0`$ and $`E_0+U`$ $`E_{up}`$ and $`E_{dn}`$ of Table I, respectively. The many body solution of this model Hamiltonian gives rise to a narrow peak at the Fermi level of the spectral function. This Kondo peak is always present at T=0, except when $`E_0`$ (or $`E_0+U`$) is closer than $`\delta `$ to the Fermi energy. In this situation, the mixed valence case, the many body feature is masked inside the broader one-electron peak . The possibility of actually seing the Kondo peak experimentally is limited by the exponential dependence of the Kondo temperature on the separation of $`E_0`$ (or $`E_0+U`$) from the Fermi energy. This may be the reason why it is not seen for most of the $`3d`$ transition metal impurities. For Co and V the present LSDA calculations give the $`m=0`$ level very close to $`E_F`$ and would therefore predict the possibility of observing a Kondo peak in both cases. However, in the case of Co the density of states for $`m=0`$ shows narrower features that for V, this last one spreading more uniformly within an energy range of more than 0.2 eV. This could explain why no special STM feature is observed in V. It is well known that the LDA approximation enhances the width of hybridization peaks, so that we are not able to decide if the observed feature in the case of Co is essentially many body (Kondo) or not. The temperature dependence would of course help to answer this question. In any case, only when one of the two levels, at $`E_0`$ or at $`E_0+U`$, is close to $`E_F`$ there is a possibility of observing this narrow feature.

When two Co impurities are brought close together  the experimental feature of the single Co impurity dissappears. In a simple tight-binding model one would expect a splitting of both the majority and the minority spin levels due to the interaction between the two impurities, giving rise to pairs of bonding and antibonding states. To perform an LSDA calculation similar to the previous one requires a much larger unit cell in order to prevent each dimer from interacting with the others. For this last reason we performed a simplified calculation, using only one monolayer of Au(111) and depositing dimer rows on one side of it. The two atoms in the dimer were placed at a separation of $`2.87\AA `$, in two next nearest neighbour hollow sites of the surface, and atoms in different dimers were situated at least $`5\AA `$ apart. For Co dimers in the ferromagnetic configuration we found a dimer magnetic moment of $`4\mu _B`$ and the $`m=0`$ majority levels appear at -3.10 and -1.27 eV, while the minority levels are at -1.54 and +0.43 eV with respect to the Fermi level of gold. None of the levels is within the experimental energy range and for this reason the Kondo temperature would be an order of magnitude or more lower than for an isolated Co impurity. The antiferromagnetic configuration of the dimer has a very small energy difference with the ferromagnetic one, and the $`3d`$ levels are also far from $`E_F`$. Therefore, the reason for not observing the narrow feature at $`E_f`$ in the dimer case may be either that the Kondo temperature is too low or that the system is antiferromagnetic, in which case there would be no Kondo peak.

In summary, we have related results of LSDA calculations with features appearing near the Fermi level when transition metal atoms or clusters are adsorbed on Au(111). The proximity of the impurity $`d`$ level with $`m=0`$ to the Fermi energy of the conduction band seems to be a necessary but not sufficient condition for observing an abrupt change in the tunneling current.

###### Acknowledgements.

We thank Dr.Javier Guevara for helpful discussions, and acknowledge CONICET and SECYT (grant PICT 03-00105-02043) of Argentina for partial support.


text:
# IIT-HEP-00/1 Energy Absorber R&D

## 1 Introduction

Ionization cooling of a muon beam can be accomplished by passing the beam through energy-absorbing material and accelerating structures, both embedded within a focusing magnetic lattice. The rate of change of the normalized transverse emittance with path length is then approximately described by

$$\frac{dϵ_n}{ds}=\frac{1}{\beta ^2}\frac{dE_\mu }{ds}\frac{ϵ_n}{E_\mu }+\frac{1}{\beta ^3}\frac{\beta _{}(0.014\mathrm{GeV})^2}{2E_\mu m_\mu L_R},$$
(1)
where $`s`$ is the path length, $`E_\mu `$ the muon energy, $`L_R`$ the radiation length of the absorber medium, $`\beta =v/c`$, and $`\beta _{}`$ is the betatron function of the beam (where the size of the beam is given by $`\sigma _x=\sigma _y=\sqrt{ϵ_n\beta _{}/\beta \gamma }`$).

In Eq. 1 we see, in addition to the $`dE/ds`$ transverse cooling term, a transverse heating term due to multiple Coulomb scattering of the muons in the absorbers. Since cooling ceases once the heating and cooling terms are equal, Eq. 1 implies an equilibrium emittance, which in principle (neglecting other limiting effects) would be reached asymptotically were the cooling channel continued indefinitely. Since the heating term is proportional to $`\beta _{}`$ and inversely proportional to the radiation length of the absorber medium, the goal of achieving as small an equilibrium emittance as possible requires us to locate the absorber only in low-$`\beta _{}`$ regions and to use a medium with the longest possible radiation length, namely hydrogen. To achieve low $`\beta _{}`$, we want the strongest possible focusing elements. We are thus led to superconducting solenoids filled with liquid hydrogen as possibly the optimal solution.111However, lithium lenses might give an even lower equilibrium emittance than solenoids with liquid hydrogen, since stronger focusing fields may be feasible with liquid-lithium lenses than with magnets, and this may overcome the radiation-length advantage of hydrogen.

## 2 Cooling channel designs

A variety of solenoidal focusing lattices have been considered for muon-beam cooling, including the so-called Alternating-Solenoid, FOFO, DFOFO, SFOFO, and Single- and Double-Flip designs . In all of these, liquid-hydrogen (LH2) absorbers are used to minimize multiple scattering. The specifications of the absorbers for some representative cases are given in Table 1.

## 3 Absorber power handling

A key problem is handling the heat deposited in the hydrogen by the muon beam. As shown in Table 1, this can exceed 100 watts per absorber. LH2 targets have been successfully operated in such a heat-deposition regime, notably the target for the SAMPLE experiment at Bates . These high-power target designs have followed the approach developed at SLAC , in which the liquid is pumped around an external cooling loop (see Fig. 1). The loop includes a heat exchanger, within which the deposited heat is transferred to a supply of cold helium gas, as well as a heater for fast regulation of temperature. The difficulty in such designs is assuring flow of the LH2 within the absorber cell that is predominantly transverse to the beam direction, since any little parcel of LH2 that flowed along the beam would be quickly heated above the boiling point. In the SAMPLE target (Fig. 1b) transverse flow was accomplished using a specially-designed perforated baffle, and for the proposed SLAC Experiment 158 target (Fig. 1a) the use of asymmetrically-perforated screens is under investigation . We are developing a design in which transverse flow within the absorber cell is produced by variously-angled nozzles placed around the periphery.

In addition to the external-loop approach, we are considering as an alternative heat exchange internal to the absorber vessel, using driven convection to provide mixing and transverse flow (Fig. 2). This has the obvious virtue of fewer moving parts. Moreover, the transverse flow arises naturally, rather than being imposed by a clever but complicated design. The key question is whether there will be sufficient convection at the anticipated power level to give the necessary rate of heat exchange via the cooling tubes at the periphery of the cell. We are addressing this by a series of computational-fluid-dynamics (CFD) calculations that numerically solve the Navier-Stokes equation on a suitable two-dimensional222Given the deflection of the absorber windows (discussed below) a three-dimensional calculation would be preferable but is impractical for any reasonable amount of computing resources. grid to evaluate the Nusselt number (the dimensionless parameter that characterizes the rate of convective heat transport) vs. the Rayleigh number (the dimensionless parameter, proportional to the power dissipation, that characterizes the degree of turbulence). Preliminary indications are that there will be sufficient convection; this will be investigated further numerically and tested experimentally on the bench.

## 4 Absorber windows

The thickness of the absorber windows is a critical parameter. They must be thick enough to sustain the pressure from the LH2, yet as thin as possible so as to minimize multiple scattering. The window thicknesses in Table 1 have been chosen based on the ASME standard for pressure vessels . This choice also satisfies the Fermilab safety code for liquid-hydrogen targets . As illustrated in Fig. 3, ASME  specifies three standard window profiles: hemispherical, ellipsoidal, and torispherical. The minimum thickness required in each case is

$`(\mathrm{hemispherical})t`$ $`=`$ $`{\displaystyle \frac{0.5PR}{SE0.1P}},s=R=0.5D`$ (2)
$`(\mathrm{ellipsoidal})t`$ $`=`$ $`{\displaystyle \frac{0.5PD}{SE0.1P}},s=0.25D`$ (3)
$`(\mathrm{torispherical})t`$ $`=`$ $`{\displaystyle \frac{0.885PD}{SE0.1P}},s=0.169D,`$ (4)
where $`P`$ is the pressure differential, $`R`$ the radius of curvature (for hemispherical windows), $`D`$ the length of the major axis (for ellipsoidal or torispherical windows), i.e. the absorber diameter, $`S`$ the maximum allowable stress, and $`E`$ the weld efficiency. In the above equations we give also $`s`$, the sagitta of the window at its center, which determines which window shapes can be used for absorbers of given dimensions: since each absorber has two windows, the absorber length $`L`$ must satisfy $`L>2s`$. (Note that for ellipsoidal windows the ASME code considers only the case where the length of the major axis is twice that of the minor axis, and for the torispherical case the radius of curvature of the “knuckle” is 6% that of the main cap; see Fig. 3.) The maximum allowable stress is the smaller of (ultimate stress)/4 or (yield stress)/1.5 . In practice we find that it is the ultimate stress that is the limit.

The hemispherical window shape minimizes the needed thickness. However, a hemispherical window is practical only for absorbers whose length exceeds their diameter. Since our absorbers are typically wider in diameter than they are long, we choose either the ellipsoidal or torispherical window profile. As an example, we show in Fig. 4 the mechanical design for the “FOFO1” absorber of Table 1. Given the 15 cm radius of the absorber, both hemispherical and ellipsoidal windows are ruled out. Torispherical windows with 5.1 cm sagitta leave just 2.4 cm for the flange and manifold assembly that joins the two windows and connects to the hydrogen inlets and outlets. Our solution is to machine each window with an integral flange out of a single block of material. The two flanges bolt together in “clamshell” fashion to form the manifold.

In the Neutrino Factory Feasibility Study  we consider FOFO and Single-Flip options for the cooling channel. In both designs the cooling performance is significantly limited by scattering in the absorber windows. To minimize this effect we propose to operate the absorbers at 1 atm pressure.333Operation of high-power LH2 absorbers at 1 atm pressure is not an established technique, 2 atm having been used in previous high-power target designs . There may be safety concerns that will prevent such operation ; this is a topic for R&D. In the Single-Flip design the absorbers are long enough to permit the use of 300 $`\mu `$m ellipsoidal windows.44footnotemark: 4 In the FOFO design two cooling sections are used having two different absorber sizes. In the first (“FOFO1”) section, 400 $`\mu `$m torispherical windows are used as discussed above.44footnotemark: 4 By the second section, the beam has become small enough to permit a reduction in absorber diameter, allowing use of 200 $`\mu `$m ellipsoidal windows.444We assumed a weld efficiency $`E=0.9`$ in specifying these window thicknesses, however in the integral-flange approach discussed above, $`E=1`$, allowing 10% thinner windows, or alternatively, operation at 1.1 atm with the thicknesses given. The reduction in window thickness results in a lowering of the equilibrium emittance from 2.6$`\pi `$ to 2.2$`\pi `$ mm$``$rad and a corresponding increase in the cooling rate.

We have also begun to explore the option of customizing the thickness profile of the window in order to minimize the thickness at the center while maximizing strength. An ANSYS finite-element calculation has been carried out that shows that the stresses in a torispherical window are greatest near the edge, in the region in which the window curvature under pressure exhibits a point of inflection (Fig. 5). By thickening the material near the edge one can reduce the maximum stresses substantially, allowing the material near the center to be thinner by perhaps as much as a factor of 5. Of course manufacturability will also impose a limit on how thin the center of a machined window can be; we will explore this soon by building and testing prototypes.

## 5 R&D issues

Beryllium or AlBeMet (a composite of 62% beryllium/38% aluminum) could reduce the impact of the windows on the cooling performance. However, based on the CEA bubble-chamber accident, beryllium is believed to be incompatible with liquid hydrogen, and an R&D program will be required to establish safe design parameters for these materials. With 40% greater strength than aluminum and 2.1 times the radiation length, AlBeMet has the potential to lower the total radiation-length fraction per absorber from 2.4% to 1.8% or less, depending on the detailed optimization of absorber dimensions. (While beryllium windows may also be feasible, there may be little additional gain in going beyond AlBeMet.) Other cooling scenarios (e.g. SFOFO) use absorbers that are thicker compared to their diameter. Here effects of windows on cooling performance are reduced, and aluminum windows may be adequate. Whether R&D on exotic window materials is worthwhile may thus depend on which cooling approach prevails.

In all scenarios the specific power dissipation in the absorbers is large and represents a substantial portion of the cryogenic load of the cooling channel. Handling this heat load is a significant design challenge. An R&D program is already in place at IIT to understand the thermal and fluid-flow aspects of maintaining a constant temperature within the absorber volume despite the large spatial and temporal variations in power density. This program is beginning with CFD studies and is planned to proceed to bench tests and high-power beam tests of absorber prototypes over the next year.

In some scenarios (especially those with emittance exchange), lithium hydride (LiH) absorbers may be called for. Since it is a solid, LiH in principle can be fabricated in arbitrary shapes. In emittance-exchange channels, dispersion in the lattice spatially separates muons according to their energies, whereupon specially shaped absorbers can be used to absorb more energy from muons of higher energy and less from those of lower energy. However, solid LiH shapes are not commercially available, and procedures for their fabrication would need to be developed. Such an effort is challenging since LiH reacts with water, releasing hydrogen gas and creating an explosion hazard.

## 6 Acknowledgements

We thank M. Boghosian of IIT who carried out the CFD calculations and Z. Tang of Fermilab who carried out the ANSYS window calculations. This work was supported in part by the U.S. Dept. of Energy, the National Science Foundation, and the IIT Research Institute.