text:
# The Galactic population and birth rate of radio pulsars

## 1. Selection effects in pulsar surveys

Fig. 1 shows the current sample of 1300 pulsars in the ATNF on-line catalogue (Hobbs et al. these proceedings) projected onto the Galactic plane.

Rather than being distributed about the Galactic centre, the majority of pulsars are local objects. Far from being representative of the true population, this sample is heavily biased by a number of selection effects which we now outline below.

The inverse square law. Like all astronomical sources, observed pulsars of a given luminosity $`L`$ are strongly selected by their apparent flux density, $`S`$. For pulsars, which beam to a certain fraction $`f`$ of $`4\pi `$ sr111In the absence of prior knowledge about beaming, geometrical factors are usually ignored and the resulting “pseudoluminosity” is quoted at some standard observing frequency; e.g., at 1400 MHz, given a measured flux density and distance, the pseudoluminosity $`L_{1400}=S_{1400}d^2`$. , $`S=L/(f4\pi d^2)`$, where $`d`$ is the distance to the pulsar. Since all pulsar surveys have some limiting flux density only those objects bright or close enough will be detectable.

The radio sky background. One limit to pulsar search sensitivities is the thermal noise in the receiver, i.e. the “system temperature”, $`T_{\mathrm{sys}}`$. While every effort is made to minimize $`T_{\mathrm{sys}}`$ at the telescope, synchroton radiating electrons in the Galactic magnetic field contribute significantly with a “sky background” component, $`T_{\mathrm{sky}}`$. At observing frequencies $`\nu 0.4`$ GHz, $`T_{\mathrm{sky}}`$ dominates $`T_{\mathrm{sys}}`$ along the Galactic plane. Fortunately, $`T_{\mathrm{sky}}\nu ^{2.8}`$ so this is effect is significantly reduced at higher frequencies; e.g. for the PM system, $`\nu =1.4`$ GHz.

Propogation effects in the intersellar medium (ISM). The ISM is a mixed blessing for pulsar astronomers. On one hand, the dispersion of pulses caused by the ionized component of the ISM permits an estimate of $`d`$ through the dispersion measure. Conversely, dispersion and scatter-broadening of the pulses conspire against detection of short period and/or distant objects. The effects of scattering are shown in Fig. 2. Fortunately, like $`T_{\mathrm{sky}}`$, the scatter-broadening time $`\tau _{\mathrm{scatt}}`$ has a strong frequency dependence, scaling roughly as $`\nu ^4`$.

As shown in Fig. 2, there is a transition frequency of about 1 GHz, below which scattering can hide a large fraction of the population. Another factor is scintillation, the modulation of apparent flux densities by refractive or diffractive “screens” of material along the line of sight (Rickett 1970). This is particularly important for nearby pulsars, where the apparent flux densities can vary significantly. For example, two northern sky surveys carried out 20 years apart with comparable sensitivity (Damashek, Taylor & Hulse 1978; Sayer, Nice & Taylor 1997) detected a number of pulsars above and below the nominal search thresholds of one experiment but not the other. Ideally, surveying the sky multiple times minimizes the effects of scintillation against pulsars nominally above the threshold, and maximizes the detections of faint pulsars through favourable scintillation.

Finite size of the emission beam. As mentioned above, the fact that pulsars do not beam to $`4\pi `$ sr means that we see only a fraction $`f`$ of the total active population. For a circular pencil beam, Gunn & Ostriker (1970) estimated $`f1/6`$. A consensus on the precise shape and evolution of the emission beam, however, has yet to be reached. Narayan & Vivekanand (1983) argued that the beam shape is elongated in the meridonial direction. Lyne & Manchester (1988), on the other hand, favour a circular beam. Using the same database, Biggs (1990) presented evidence in favour of meridonal compression! All of these studies do agree that the beam size is period dependent, with shorter period pulsars having larger beaming fractions. For example, Tauris & Manchester (1998) found that $`f=0.09\left[\mathrm{log}(P/\mathrm{sec})1\right]^2+0.03`$, where $`P`$ is the period. Undoubtedly, a complete model for $`f`$ needs to account for other factors, such as evolution of the inclination angle between the spin and magnetic axes. Given the uncertainties, most authors quote results with and without a beaming correction.

Pulse nulling. The abrupt cessation of the pulsed emission for many pulse periods, was first identified by Backer (1970). Ritchings (1976) presented evidence that the incidence of nulling became more frequent in older long-period pulsars, suggesting that it signified the onset of the final stages of the neutron star’s life as an active radio pulsar. Since most pulsar surveys have short ($`<`$ few min) integration times, there is an obvious selection effect against nulling objects. Means of combatting this effect are to look for individual pulses in search data (see e.g. Nice 1999), survey the sky many times, or use longer integrations. The PM survey, which employs 35-min pointings, is proving particularly effective at detecting nulling pulsars and should soon be able to better quantify this population and provide a more satisfactory understanding of nulling pulsars.

Orbital motion. In standard pulsar searches, where time series are Fourier transformed, the signal from a binary pulsar can be Doppler shifted over several bins in the Fourier domain. In extreme cases, where the survey integration time is a significant fraction of the orbital period, this results in a loss of sensitivity. For example the loss of signal-to-noise of the original binary pulsar B1913+16 during a 35-min observation of the PM survey can be as much as 90%. Correcting for this effect using “acceleration searches”, is now becoming more routine, thanks to the ever-increasing availability of high-speed computer resources. Since it is quite possible that young pulsars in tight binary orbits exist, the very deep ($``$ 2–10 hr) searches for young pulsars described by Camilo in this volume are now being re-analysed with full acceleration searches.

To get an idea of how biased the sample is due to the above effects, Fig. 1 shows the cummulative distribution of pulsars as a function of distance from the Sun projected onto the Galactic plane. Also shown is the expected distribution for a simulated population in which there are no selection effects. As can be seen, the two samples are closely matched only out to a kpc or so before the selection effects become significant. From these curves, we deduce that less than 10% of the potentially observable population in the Galaxy are currently detectable. Rigorous conclusions about the true pulsar population can only be made after properly accounting for these selection effects.

## 2. Landmark papers in pulsar statistics

Shortly after the discovery of pulsars, a number of authors began to consider their implications for the Galactic population of neutron stars. A widely cited paper from that era is the work of Gunn & Ostriker (1970; hereafter GO). To tackle the problem analytically, GO made two simplifying assumptions based on the sample of 52 pulsars known at that time: (a) the relationship $`LB^2`$, where $`L`$ is the radio luminosity, and $`B`$ is the dipole magnetic field strength; (b) the evolution of $`B`$ being an exponential decay with a time constant $`t_d`$. With these in hand, and using the dipolar spin-down expression $`B^2P\dot{P}`$, GO derived expected distributions for the observed population. In particular, they showed that the observed number of pulsars $`N_{\mathrm{obs}}=\pi D^2F\mathrm{\Sigma }t_d\mathrm{exp}(2\sigma _P^2)/2`$, where $`D`$ is the mean distance of the observed sample, $`F`$ is a parameter relating to the completeness of the surveys, $`\mathrm{\Sigma }`$ is the local birth rate of potentially observable pulsars and $`\sigma _P`$ is derived from their model fit to the observed period distribution. Assuming $`f=1/6`$, and extrapolating the local birth rate over the Galaxy, they arrived at a Galactic birth rate $`=1/(30\mathrm{yr})`$, in good agreement with the best estimates of the rate of type-II supernovae at that time (e.g. Blaauw 1961).

Although mathematically appealing, GO’s analytical approach required a number of simplifying assumptions about the pulsar population and the Galaxy itself. One of these assumptions, the spontaneous decay of the magnetic field, continues to be extremely controversial. Since GO’s original study, many papers have been written presenting arguments for and against field decay (see e.g. Lyne, Manchester & Taylor 1985; Bailes 1989; Narayan & Ostriker 1990; Bhattacharya et al. 1992; see also the contribution by van Leeuwen et al. in this volume). A less model-dependent approach to the problem, first developed by Large (1971), can be summarized by the following expression:

$$dn(P,z,R,L)=V(P,z,R,L)\rho (P,z,R,L)dPdzdRdL.$$

Here $`N`$ is the observed population of pulsars as a function of period, $`P`$, distance from the Galactic plane, $`z`$, Galactocentric radius, $`R`$ and luminosity, $`L`$. The quantity $`V`$ represents the volume of the Galaxy effectively searched and $`\rho `$ is the underlying (true) distribution of the population. Since we know $`n`$ and can estimate $`V`$ on the basis of (hopefully) well-understood survey sensitivities, we can invert the above expression to solve for $`\rho `$. The only simplification required to do this is to assume222Following a question about this assumption after my talk, I examined the current pulsar sample. Apart from a weak coupling between $`P`$ and $`z`$, there are no significant relationships between any of these parameters. The assumption of independent distributions seems to be well founded. that $`P`$, $`z`$, $`R`$ and $`L`$ are independent quantities. The problem then reduces to four equations which can be solved for the underlying distributions of interest: $`\rho _P(P)`$, $`\rho _z(z)`$, $`\rho _R(R)`$ and $`\rho _L(L)`$.

Large’s method was somewhat ahead of its time: back in 1971, pulsar surveys of the Galaxy were still in their early stages so that $`V`$ was not well determined. By the late 1970s, however, a number of large-scale searches had been carried out and Taylor & Manchester (1977; hereafter TM) applied the above technique to the sample of $`150`$ pulsars then known (see also Davies, Lyne & Seiradakis 1977). Integrating the derived distribution functions over the Galaxy, and assuming $`f=0.2`$, TM estimated the total number of active pulsars in the Galaxy to range between 60,000 and 850,000. Here the principle uncertainty is the assumed distance model. TM considered a uniform electron content with a mean electron density in the range $`0.02<n_e<0.03`$ cm-3. To calculate the Galactic birth rate, $``$, TM required an estimate of the mean pulsar lifetime, $`T`$, which they obtained from an analysis of the $`z`$ distribution as a function of characteristic age $`\tau =P/(2\dot{P})`$. In the $`z\tau `$ diagram (see their Fig. 7), TM argued that the discrepancy between the expected and observed $`z`$ for characteristic ages larger than a few Myr set a limit to $`T=4`$ Myr. This leads to an implied birthrate $`=1/(6\mathrm{yr})`$! Davies et al. (1977) reached similar conclusions.

The rather high birth rates from these and other analyses prompted Phinney & Blandford (1981) and in particular Vivekanand & Narayan (1981; hereafter VN) to consider a less model-dependent approach to estimating $``$. The method involves binning the observed sample as a function of spin period, $`P`$, and in each bin computing the flow or “current” of pulsars, $`J`$, through the bin:

$$J(P)=\frac{1}{\mathrm{\Delta }P}\underset{i=1}{\overset{n_{\mathrm{bin}}}{}}\frac{\dot{P}_i\xi _i}{f_i}.$$

Here, $`n_{\mathrm{bin}}`$ is the number of pulsars in a period bin of width $`\mathrm{\Delta }P`$, $`\xi _i`$ and $`f_i`$ are the “scale factor” and beaming fraction of the $`i^{\mathrm{th}}`$ pulsar respectively. As discussed earlier, $`f_i`$ is based on some beaming model. For a given pulsar, its scale factor $`\xi _i`$ represents the number of pulsars with similar parameters in the Galaxy. In principle, this is very similar to the $`V/V_{\mathrm{max}}`$ first used to correct observationally biased samples of quasars (Schmidt 1968). In practice $`\xi `$ is computed using a Monte Carlo simulation of $`N`$ pulsars with identical periods and luminosities. Using accurate models for the various pulsar surveys, it is relatively straightforward to calculate the number of pulsars $`n`$ that are detectable from that population. As a result, $`\xi =N/n`$. Detailed simulations to test this approach show that the scale factors give reliable results about the population of pulsars with luminosities above $`L_{\mathrm{min}}`$, the minimum luminosity in the sample (see Lorimer et al. 1993). We note in passing that it is possible to form the luminosity function for $`L>L_{\mathrm{min}}`$ by weighting the observed luminosities by the appropriate scale factors. Integrating this distribution then yields $`N(>L_{\mathrm{min}})`$.

The beauty of the pulsar-current analysis is that it makes just two fundamental assumptions about pulsars: (i) they are a steady-state population; (ii) they are spinning down steadily from short to long periods. The first of these is justified since the ages of pulsars, while not well known ($`10^{78}`$ yr), are certainly less than the age of the Galaxy, $`10^{10}`$ yr. The second is, of course, well in accord with timing observations. The birth-rate can be computed from this analysis by simply plotting $`J`$ as a function of $`P`$. In the standard model where pulsars are born spinning rapidly, there should be a peak in the current at short periods followed by a decline in the current as pulsars end their life with longer periods. The birth rate is then just the height of this peak. VN derived $`=1/(1627\mathrm{yr})`$, in much better agreement with the supernova rate.

A by-product of VN’s analysis was their conclusion that a significant fraction of pulsars are “injected” into the population with relatively long initial spin periods ($`P_00.5`$ s), rather than the standard picture of birth with $`P_020`$ ms. Rather like field decay, arguments for and against injection have been presented ever since. Following criticism by Lyne, Manchester & Taylor (1985; hereafter LMT) that VN’s analysis had not properly taken account of selection effects such as scattering, Narayan (1987) performed a more detailed analysis which provided further support for injection. Lorimer et al. (1993) used Monte Carlo simulations to investigate the validity of the pulsar current approach and found that it was prone to systematic errors induced by faint, nearby pulsars. Excluding these objects from the analysis significantly reduces the case for injection. Perhaps the final word on this issue will be a detailed pulsar current analysis of the PM survey sample (Vranesevic et al. this volume).

A major effort to quantify the Galactic distribution was the analysis of LMT who combined the theoretical framework of GO with the numerical approach of Large to the sample of 316 pulsars detected in four major 400-MHz surveys. Two key results from their paper are shown in Fig. 3.

As a result of their numerical analysis, LMT were able to place better constraints on the radial distribution than TM. However, due to the extreme selection effects on 400-MHz surveys towards the inner regions of the Galaxy, the corrected radial distribution becomes very uncertain below 4 kpc (see Fig. 3).

To derive the birth rate, $``$, rather than applying a pulsar current anlalysis, LMT followed GO’s approach to model the observed distributions and concluded that the best-fit timescale for exponential field decay to be 9.1 Myr. Since their analysis assumed GO’s $`LB^2`$ luminosity law, magnetic field decay provided a mechanism for luminosity decay. By considering their corrected luminosity function as a steady-state population where pulsars flowed into successively fainter luminosity bins, LMT were able to derive $``$ as a function of $`L`$. The result of this analysis is shown in Fig. 3, from which LMT concluded $`=1/(30250\mathrm{yr})`$.

While an elegant approach to the problem, LMTs method relies on the field decay hypothesis being correct. As mentioned previously, the issue of magnetic field decay remains unresolved. Since magnetic field evolution underpins much of pulsar statistics, perhaps the single most important breakthrough in this area would be a comprehensive reassesment of this issue.

## 3. Recent progress in the Galactic distribution of pulsars

Although much effort has gone into improving the Monte Carlo simulations of pulsar population modeling since LMT, relatively little progress has been made in improving our knowledge of the Galactic distribution. Of particular interest is the underlying density $`\rho _R`$ of pulsars as a function of Galactocentric radius, $`R`$. LMT and others were only able to place poor constraints on this function (Fig. 3) and most subsequent work has assumed a Gaussian distribution for $`\rho _R`$ (e.g. Narayan 1987). In fact, as pointed out by Bailes & Kniffen (1991), there is no reason to prefer this function over one where $`\rho _R0`$ as $`R0`$.

Improving our understanding of $`\rho _R`$ at small $`R`$ requires better statistics of the inner-Galaxy pulsars. Due to the propogation and sky-background selection effects mentioned earlier, low-frequency (0.4 GHz) surveys are very poor probes of this population. One of the main motivations for high-frequency (1.4 GHz) surveys of the Galactic plane is that they are less prone to these effects. Johnston (1994) analysed two such surveys carried out in the late 1980s (Clifton et al. 1992; Johnston et al. 1992) and found that models with a deficit of pulsars in the inner Galaxy were strongly preferred over a simple Gaussian profile for $`\rho _R`$.

The sample of $`150`$ pulsars detected by these two early surveys has now been completely surpassed by the phenomenal success of the PM survey: over 700 new pulsars have been found in the Galactic plane search so far. Together with re-detections of known pulsars, the PM sample amounts to 914 pulsars. I have applied Large’s numerical approach to this new sample

to derive new distribution functions in $`R`$, $`L`$, $`z`$ and $`P`$. Preliminary results are shown in Fig. 4. The corrected radial density function clearly supports Johnston’s conclusion for a deficit of pulsars in the inner Galaxy. More work is required in quantifying the significance of this result, particularly in the $`R=0.5`$ kpc bin where only one pulsar is presently known! While a population of inner-Galaxy pulsars could be masked from the PM survey by severe scattering not currently taken into account, it is interesting that the now completed Effelsberg 5-GHz Galactic centre survey (Klein et al. these proceedings) has not found a single pulsar.

Although the PM survey has been the most prolific probe of the Galactic population to date, we are still only scratching the surface. One of many uncertain areas about the pulsar population is the shape of the luminosity function. From Fig. 4 we see a clear departure from the $`d\mathrm{log}N/d\mathrm{log}L=1`$ relationship at low $`L`$. Integrating Fig. 4 results in a Galactic population of $`25,000\pm 2000`$ potentially observable pulsars with $`L_{1400}>0.3`$ mJy kpc2. Below this limit, the population is essentially unknown. Future surveys with the SKA, which should easily be able to detect $`15,000`$ pulsars with $`L_{1400}>0.3`$ mJy kpc2 should also detect many fainter objects and truly constrain the pulsar luminosity function.

### Acknowledgments.

I thank the Particle Physics and Astronomy Research Council and the Royal Society for supporting my attendance at this meeting.

## References

Backer D. C., 1970, Nature, 228, 42

Bailes M. & Kniffen D. A., 1992, ApJ, 391, 659

Bailes M., 1989, ApJ, 342, 917

Bhattacharya D., et al. 1992, A&A, 254, 198

Biggs J. D., 1990, MNRAS, 245, 514

Blaauw A., 1961, Bull. Astr. Inst. Netherlands, 15, 265

Clifton T. R., et al. 1992, MNRAS, 254, 177

Damashek M., Taylor J. H., & Hulse R. A., 1978, ApJ, 225, L31

Davies J. G., Lyne A. G., & Seiradakis J. H., 1977, MNRAS, 179, 635

Gunn J. E. & Ostriker J. P., 1970, ApJ, 160, 979 (GO)

Johnston S., 1994, MNRAS, 268, 595

Johnston S., et al. 1992, MNRAS, 255, 401

Large M. I., 1971. IAU Symposium No. 46, Dordrecht, Reidel, p. 165

Lorimer D. R., et al. 1993, MNRAS, 263, 403

Lyne A. G. & Manchester R. N., 1988, MNRAS, 234, 477

Lyne A. G., Manchester R. N., & Taylor J. H., 1985, MNRAS, 213, 613 (LMT)

Narayan R. & Ostriker J. P., 1990, ApJ, 352, 222

Narayan R. & Vivekanand M., 1983, A&A, 122, 45

Narayan R., 1987, ApJ, 319, 162

Nice D. J., 1999, ApJ, 513, 927

Phinney E. S. & Blandford R. D., 1981, MNRAS, 194, 137

Rickett B. J., 1970, MNRAS, 150, 67

Ritchings R. T., 1976, MNRAS, 176, 249

Sayer R. W., Nice D. J., & Taylor J. H., 1997, ApJ, 474, 426

Schmidt M., 1968, ApJ, 151, 393

Tauris T. M. & Manchester R. N., 1998, MNRAS, 298, 625

Taylor J. H. & Manchester R. N., 1977, ApJ, 215, 885 (TM)

Vivekanand M. & Narayan R., 1981, JA&A, 2, 315 (VN)


format: text

text:
# Metropolis Algorithms in Generalized Ensemble

## I Introduction

In complex systems such as spin systems and protein systems, conventional simulations in the canonical ensemble will get trapped in states of energy local minima at low temperatures. We employ the generalized-ensemble algorithms in order to overcome this multiple-minima problem (for reviews, see Refs. RevHO RevSO ). In a generalized-ensemble simulation, each state is weighted by a non-Boltzmann probability weight factor so that a random walk in potential energy space may be realized. The random walk allows the simulation to escape from any energy barrier and to sample much wider configurational space than by conventional methods. Monitoring the energy in a single simulation run, one can obtain not only the global-minimum-energy state but also canonical ensemble averages as functions of temperature by the single-histogram FS1  and multiple-histogram FSWHAM  reweighting techniques.

One of the most well-known generalized-ensemble methods is perhaps multicanonical algorithm (MUCA) MUCA . (The method is also referred to as entropic sampling Lee , adaptive umbrella sampling MZ  of the potential energy BK , random walk algorithm Landau , and density of states Monte Carlo dePablo .) MUCA was first introduced to the molecular simulation field in Ref. HO . Since then MUCA has been extensively used in many applications in protein and related systems (for a review, see, e.g., Ref. RevHO ).

The replica-exchange method (REM) RE1 ; RE2  is another widely used generalized-ensemble algorithm. (Closely related methods were independently developed in Refs. RE3 JWK . REM is also referred to as multiple Markov chain method Whit  and parallel tempering STrev . For recent reviews with detailed references about the method, see, e.g., Refs. RevMSO ; IBArev .) REM has also been introduced to protein systems H97 Berne .

Both MUCA and REM are already very powerful, but we have also developed several new generalized-ensemble algorithms as further extensions of MUCA and/or REM SKO BNO03 .

In this article, we first describe the two familiar methods: MUCA and REM. We then present some of our new generalized-ensemble algorithms. The effectiveness of these methods is illustrated with a 2-dimensional Potts model, Lennard-Jones fluid system, and protein system.

## II Methods

In the regular canonical ensemble with a given inverse temperature $`\beta 1/k_BT`$ ($`k_B`$ is the Boltzmann constant), the probability distribution of potential energy $`E`$ is given by

$$P_B(E;T)n(E)W_B(E;T)n(E)e^{\beta E},$$
(1)
where $`n(E)`$ is the density of states. Since the density of states $`n(E)`$ is a rapidly increasing function of $`E`$ and the Boltzmann factor $`W_B(E;T)`$ decreases exponentially with $`E`$, the probability distribution $`P_B(E;T)`$ has a bell-like shape in general. A Monte Carlo (MC) simulation based on the Metropolis algorithm Metro  generates states in the canonical ensemble with the following transition probability from a state $`x`$ with energy $`E`$ to a state $`x^{}`$ with energy $`E^{}`$:

$$w(xx^{})=\mathrm{min}(1,\frac{W_B(E^{};T)}{W_B(E;T)})=\mathrm{min}(1,e^{\beta (E^{}E)}).$$
(2)
However, it is very difficult to obtain canonical distributions at low temperatures with this conventional Metropolis algorithm. This is because the thermal fluctuations at low temperatures are small and the simulation will certainly get trapped in states of energy local minima.

In the “multicanonical ensemble” MUCA , on the other hand, the probability distribution of potential energy is defined as follows so that a uniform flat distribution of $`E`$ may be obtained:

$$P_{mu}(E)n(E)W_{mu}(E)\mathrm{constant}.$$
(3)
Hence, the multicanonical weight factor $`W_{mu}(E)`$ is inversely proportional to the density of states, and the Metropolis criterion for the multicanonical MC simulations is based on the following transition probability:

$$w(xx^{})=\mathrm{min}(1,\frac{W_{mu}(E^{})}{W_{mu}(E)})=\mathrm{min}(1,\frac{n(E)}{n(E^{})}).$$
(4)
Because the MUCA weight factor $`W_{mu}(E)`$ is not a priori known, however, one has to determine it for each system by iterations of trial simulations.

After the optimal MUCA weight factor is obtained, one performs a long MUCA simulation once. By monitoring the potential energy throughout the simulation, one can find the global-minimum-energy state. Moreover, by using the obtained histogram $`N_{\mathrm{mu}}(E)`$ of the potential energy distribution $`P_{\mathrm{mu}}(E)`$, the expectation value of a physical quantity $`A`$ at any temperature $`T=1/k_\mathrm{B}\beta `$ can be calculated from

$$<A>_T=\frac{{\displaystyle \underset{E}{}}A(E)n(E)e^{\beta E}}{{\displaystyle \underset{E}{}}n(E)e^{\beta E}},$$
(5)
where the best estimate of the density of states is given by the single-histogram reweighting techniques (see Eq. (3)) FS1 :

$$n(E)=\frac{N_{\mathrm{mu}}(E)}{W_{\mathrm{mu}}(E)}.$$
(6)
The system for replica-exchange method (REM) RE1 ; RE2  consists of $`M`$ non-interacting copies, or replicas, of the original system in canonical ensemble at $`M`$ different temperatures $`T_m`$ ($`m=1,\mathrm{},M`$). We arrange the replicas so that there is always one replica at each temperature. Then there is a one-to-one correspondence between replicas and temperatures. Let $`X=\{\mathrm{},x_m^{[i]},\mathrm{}\}`$ stand for a state in this generalized ensemble. Here, the superscript $`i`$ and the subscript $`m`$ in $`x_m^{[i]}`$ label the replica and the temperature, respectively. A simulation of REM is then realized by alternately performing the following two steps. Step 1: Each replica in the canonical ensemble at a fixed temperature is simulated simultaneously and independently for a certain number of MC steps. Step 2: A pair of replicas, say $`i`$ and $`j`$, which are at neighboring temperatures, say $`T_m`$ and $`T_{m+1}`$, respectively, are exchanged: $`X=\{\mathrm{},x_m^{[i]},\mathrm{},x_{m+1}^{[j]},\mathrm{}\}X^{}=\{\mathrm{},x_m^{[j]},\mathrm{},x_{m+1}^{[i]},\mathrm{}\}`$. The transition probability of this replica exchange is given by the following Metropolis criterion:

$$w(XX^{})=\mathrm{min}(1,e^\mathrm{\Delta }),$$
(7)
where

$$\mathrm{\Delta }\left(\beta _{m+1}\beta _m\right)\left(E\left(q^{[i]}\right)E\left(q^{[j]}\right)\right).$$
(8)
From the results of a long REM production run, one can obtain the canonical ensemble average of a physical quantity $`A`$ as a function of temperature from Eq. (5), where the density of states is given by the multiple-histogram reweighting techniques FSWHAM  as follows. Let $`N_m(E)`$ and $`n_m`$ be respectively the potential-energy histogram and the total number of samples obtained at temperature $`T_m=1/k_\mathrm{B}\beta _m`$ ($`m=1,\mathrm{},M`$). The best estimate of the density of states is then given by FSWHAM

$$n(E)=\frac{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1N_m(E)}{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1n_me^{f_m\beta _mE}},$$
(9)
where

$$e^{f_m}=\underset{E}{}n(E)e^{\beta _mE}.$$
(10)
Here, $`g_m=1+2\tau _m`$, and $`\tau _m`$ is the integrated autocorrelation time at temperature $`T_m`$. Note that Eqs. (9) and (10) are solved self-consistently by iteration FSWHAM  to obtain the dimensionless Helmholtz free energy $`f_m`$ and the density of states $`n(E)`$.

We now introduce new generalized-ensemble algorithms that combine the merits of MUCA and REM. In the replica-exchange multicanonical algorithm (REMUCA) SO3 ; MSO03  we first perform a short REM simulation (with $`M`$ replicas) to determine the MUCA weight factor and then perform with this weight factor a regular MUCA simulation with high statistics. The first step is accomplished by the multiple-histogram reweighting techniques FSWHAM . Let $`N_m(E)`$ and $`n_m`$ be respectively the potential-energy histogram and the total number of samples obtained at temperature $`T_m=1/k_\mathrm{B}\beta _m`$ of the REM run. The density of states $`n(E)`$, or the inverse of the MUCA weight factor, is then given by solving Eqs. (9) and (10) self-consistently by iteration FSWHAM . The formulation of REMUCA is simple and straightforward, but the numerical improvement is great, because the weight factor determination for MUCA becomes very difficult by the usual iterative processes for complex systems.

While multicanonical simulations are usually based on local updates, a replica-exchange process can be considered to be a global update, and global updates enhance the sampling further. Here, we present a further modification of REMUCA and refer to the new method as multicanonical replica-exchange method (MUCAREM) SO3 ; MSO03 . In MUCAREM the final production run is not a regular multicanonical simulation but a replica-exchange simulation with a few replicas in the multicanonical ensemble. Because multicanonical simulations cover much wider energy ranges than regular canonical simulations, the number of required replicas for the production run of MUCAREM is much less than that for the regular REM, and we can keep the merits of REMUCA (and improve the sampling further). The details of REMUCA and MUCAREM can be found in Ref. RevMSO

Besides canonical ensemble, MC simulations in isobaric-isothermal ensemble mcd72  are also extensively used. This is because most experiments are carried out under the constant pressure and constant temperature conditions. The distribution $`\mathrm{P}_{NPT}(E,V)`$ for $`E`$ and $`V`$ is given by

$$\mathrm{P}_{NPT}(E,V)=n(E,V)\mathrm{e}^{\beta _0H}.$$
(11)
Here, the density of states $`n(E,V)`$ is given as a function of both $`E`$ and $`V`$, and $`H`$ is the “enthalpy”:

$$H=E+P_0V,$$
(12)
where $`P_0`$ is the pressure at which simulations are performed. This ensemble has bell-shaped distributions in both $`E`$ and $`V`$.

We now introduce the idea of the multicanonical technique into the isobaric-isothermal ensemble MC method and refer to this generalized-ensemble algorithm as the multibaric-multithermal algorithm OO03 . This MC simulation performs random walks in volume space as well as in potential energy space.

In the multibaric-multithermal ensemble, each state is sampled by a weight factor $`W_{\mathrm{mbt}}(E,V)\mathrm{exp}\{\beta _0H_{\mathrm{mbt}}(E,V)\}`$ ($`H_{\mathrm{mbt}}`$ is referred to as the multibaric-multithermal enthalpy) so that a uniform distribution in both potential energy and volume is obtained:

$$\mathrm{P}_{\mathrm{mbt}}(E,V)=n(E,V)W_{\mathrm{mbt}}(E,V)=\mathrm{constant}.$$
(13)
We call $`W_{\mathrm{mbt}}(E,V)`$ the multibaric-multithermal weight factor.

In order to perform the multibaric-multithermal MC simulation, we follow the conventional isobaric-isothermal MC techniques mcd72 . In this method, we perform Metropolis sampling on the scaled coordinates $`𝒔_i=L^1𝒓_i`$ ($`𝒓_i`$ are the real coordinates) and the volume $`V`$ (here, the particles are placed in a cubic box of a side of size $`L\sqrt[3]{V}`$). The trial moves of the scaled coordinates from $`𝒔_i`$ to $`𝒔_{}^{}{}_{i}{}^{}`$ and of the volume from $`V`$ to $`V^{}`$ are generated by uniform random numbers. The enthalpy is accordingly changed from $`H(E(𝒔^{(N)},V),V)`$ to $`H^{}(E(𝒔^{(N)},V^{}),V^{})`$ by these trial moves. The trial moves will be accepted with the probability

$$w(xx^{})=\mathrm{min}(1,\mathrm{exp}[\beta _0\{H^{}HNk_BT_0\mathrm{ln}(V^{}/V)\}]),$$
(14)
where $`N`$ is the total number of particles in the system.

Replacing $`H`$ by $`H_{\mathrm{mbt}}`$, we can perform the multibaric-multithermal MC simulation. The trial moves of $`𝒔_i`$ and $`V`$ are generated in the same way as in the isobaric-isothermal MC simulation. The multibaric-multithermal enthalpy is changed from $`H_{\mathrm{mbt}}(E(𝒔^{(N)},V),V)`$ to $`H_{\mathrm{mbt}}^{}(E(𝒔^{(N)},V^{}),V^{})`$ by these trial moves. The trial moves will now be accepted with the probability

$$w(xx^{})=\mathrm{min}(1,\mathrm{exp}[\beta _0\{H_{\mathrm{mbt}}^{}H_{\mathrm{mbt}}Nk_BT_0\mathrm{ln}(V^{}/V)\}]).$$
(15)
While MUCA yields a flat distribution in potential energy and performs a random walk in potential energy space, we can, in principle, choose any other variable and induce a random walk in that variable. One such example is the multi-overlap algorithm BNO03 . Here, we choose a protein system and define the overlap in the space of dihedral angles by, as it was already used in HMO97 ,

$$q=(nd)/n,$$
(16)
where $`n`$ is the number of dihedral angles and $`d`$ is the distance between configurations defined by

$$d=vv^1=\frac{1}{\pi }\underset{i=1}{\overset{n}{}}d_a(v_i,v_i^1).$$
(17)
Here, $`v_i`$ is our generic notation for the dihedral angle $`i`$, $`\pi <v_i\pi `$, and $`v^1`$ is the vector of dihedral angles of the reference configuration. The distance $`d_a(v_i,v_i^{})`$ between two angles is defined by

$$d_a(v_i,v_i^{})=\mathrm{min}(|v_iv_i^{}|,2\pi |v_iv_i^{}|).$$
(18)
We want to simulate the system with weight factors that lead to a flat distribution in the dihedral distance $`d`$, and hence to a random walk process in $`d`$:

$$d<d_{\mathrm{min}}d>d_{\mathrm{max}}\mathrm{and}\mathrm{back}.$$
(19)
Here, $`d_{\mathrm{min}}`$ is chosen sufficiently small so that one can claim that the reference configuration has been reached. The value of $`d_{\mathrm{max}}`$ has to be sufficiently large to introduce a considerable amount of disorder.

Moreover, we can define a weight factor that leads to a random-walk process between two configurations. This multi-overlap simulation allows a detailed study of the transition states between the two configurations, whereas a random walk in energy space of a regular MUCA simulation may miss the transition state (see Ref. BNO03  for details).

## III Results

We now present the results of our simulations based on the algorithms described in the previous section.

The first example is a spin system. We studied the 2-dimensional 10-state Potts model NSMO . The lattice size was $`34\times 34`$. This system exhibits a first-order phase transition Bax . In Fig. 1 we show the probability distributions of energy at three tempeartures (above the critical temperature $`T_C`$, at $`T_C`$, and below $`T_C`$) and average energy and specific heat as functions of inverse temperature. All these results imply that the system indeed undergoes a first-order phase transition.

Iterations of MUCAREM (and REMUCA) can be used to obtain an optimal MUCA weight factor MSO03 . In Fig. 2 we show the results of our MUCA weight factor determination by MUCAREM. We first made a REM simulation of 10,000 MC sweeps (for each replica) with 32 replicas (Fig. 2(a)). Using the obtained energy distributions, we determined the (preliminary) MUCA weight factor, or the density of states, by the multiple-histogram reweighting techniques of Eqs. (9) and (10). Because the trials of replica exchange are not accepted near the critical temperature for first-order phase transitions, the probability distributions in Fig. 2(a) for the energy range from $`1.5`$ to $`1.0`$ fails to have sufficient overlap, which is required for successful application of REM. This means that the MUCA weight factor in this energy range thus determined is of “poor quality.” With this MUCA weight factor, however, we made iterations of three MUCAREM simulations of 10,000 MC sweeps (for each replica) with 8 replicas (Fig. 2(b) for the first iteration and Fig. 2(c) for the third iteration). In Fig. 2(b) we see that the distributions are not completely flat, reflecting the poor quality in the phase-transition region. This problem is rapidly rectified as iterations continue, and the distributions are completely flat in Fig. 2(c), which gives an optimal MUCA weight factor in the entire energy range by the multiple-histogram reweighting techniques.

Besides MUCAREM the methods of Berg MUCAW  and Wang-Landau Landau  are also effective for the determination of the MUCA weight factor (or the density of states). In Fig. 3 we compare how fast these three methods converge to yield an optimal MUCA weight factor, or a flat distribution. Each figure shows three curves superimposed that correspond to the (hot-start) MUCA simulation with the weight factor that was “frozen” after 300,000 MC seeps, 500,000 MC sweeps, and 1,000,000 MC sweeps of iterations of the weight factor determination. While the results of MUCAREM and Berg’s method are similar, those of Wang-Landau method have quite different behavior. After 300,000 MC sweeps, MUCAREM and Berg’s method gives reasonably flat distribution for $`E>1.0`$ and $`E>1.3`$, respectively, whereas Wang-Landau method gives a distribution that is quite spiky (but already covers low-energy regions). After 500,000 MC sweeps, MUCAREM essentially gives a flat distribution in the entire energy range and Berg’s method for $`E>1.7`$, while Wang-Landau method also covers the entire range (though still very rugged). After 1,000,000 MC sweeps, the three methods all give reasonably flat distributions. Details will be published elsewhere NSMO .

We now present the results of our multibaric-multithermal simulation. We considered a Lennard-Jones 12-6 potential system. We used 500 particles ($`N=500`$) in a cubic unit cell with periodic boundary conditions. The length and the energy are scaled in units of the Lennard-Jones diameter $`\sigma `$ and the minimum value of the potential $`ϵ`$, respectively. We use an asterisk ($``$) for the reduced quantities such as the reduced length $`r^{}=r/\sigma `$, the reduced temperature $`T^{}=k_\mathrm{B}T/ϵ`$, the reduced pressure $`P^{}=P\sigma ^3/ϵ`$, and the reduced number density $`\rho ^{}=\rho \sigma ^3`$ ($`\rho N/V`$).

We started the iterations of the multibaric-multithermal weight factor determination from a regular isobaric-isothermal simulation at $`T_0^{}=2.0`$ and $`P_0^{}=3.0`$. In one MC sweep we made the trial moves of all particle coordinates and the volume ($`N+1`$ trial moves altogether). For each trial move the Metropolis evaluation of Eq. (15) was made. Each iteration of the weight factor determination consisted of 100,000 MC sweeps. In the present case, it was required to make 12 iterations to get an optimal weight factor $`W_{\mathrm{mbt}}(E,V)`$. We then performed a long multibaric-multithermal MC simulaton of 400,000 MC sweeps with this $`W_{\mathrm{mbt}}(E,V)`$.

Figure 4 shows the probability distributions of $`E^{}/N`$ and $`V^{}/N`$. Figure 4(a) is the probability distribution $`\mathrm{P}_{NPT}(E^{}/N,V^{}/N)`$ from the isobaric-isothermal simulation first carried out in the process (i.e., $`T_0^{}=2.0`$ and $`P_0^{}=3.0`$). It is a bell-shaped distribution. On the other hand, Fig. 4(b) is the probability distribution $`\mathrm{P}_{\mathrm{mbt}}(E^{}/N,V^{}/N)`$ from the multibaric-multithermal simulation finally performed. It shows a flat distribution, and the multibaric-multithermal MC simulation indeed sampled the configurational space in wider ranges of energy and volume than the conventional isobaric-isothermal MC simulation.

Figure 5 shows the time series of $`V^{}/N`$. In Fig. 5(a) we show the results of the conventional isobaric-isothermal simulations at $`(T^{},P^{})=(2.0,2.2)`$ and (2.0,3.8), while in Figure 5(b) we give those of the multibaric-multithermal simulation. The volume fluctuations in the conventional isobaric-isothermal MC simulations are only in the range of $`V^{}/N=1.31.4`$ and $`V^{}/N=1.51.6`$ at $`P^{}=3.8`$ and at $`P^{}=2.2`$, respectively. On the other hand, the multibaric-multithermal MC simulation performs a random walk that covers even a wider volume range.

We calculated the ensemble averages of potential energy per particle, $`<E^{}/N>_{NPT}`$, and density, $`<\rho ^{}>_{NPT}`$, at various temperature and pressure values by the reweighting techniques. They are shown in Fig. 6. The agreement between the multibaric-multithermal data and isobaric-isothermal data are excellent in both $`<E^{}/N>_{NPT}`$ and $`<\rho ^{}>_{NPT}`$.

The important point is that we can obtain any desired isobaric-isothermal distribution in wide temperature and pressure ranges ($`T^{}=1.62.6`$, $`P^{}=2.23.8`$) from a single simulation run by the multibaric-multithermal MC algorithm. This is an outstanding advantage over the conventional isobaric-isothermal MC algorithm, in which simulations have to be carried out separately at each temperature and pressure, because the reweighting techniques based on the isobaric-isothermal simulations can give correct results only for narrow ranges of temperature and pressure values.

The third example is a system of a biopolymer. A brain peptide Met-enkephalin has the amino-acid sequence Tyr-Gly-Gly-Phe-Met. We fix the peptide-bond dihedral angles $`\omega `$ to $`180^{}`$, which implies that the total number of variable dihedral angles is $`n=19`$. We neglect the solvent effects as in previous works. The low-energy configurations of Met-enkephalin in the gas phase have been classified into several groups of similar structures OKK92 ; MHO98 . Two reference configurations, called configuration 1 and configuration 2, are used in the following and depicted in Fig. 7. Configuration 1 has a $`\beta `$-turn structure with hydrogen bonds between Gly-2 and Met-5, and configuration 2 a $`\beta `$-turn with a hydrogen bond between Tyr-1 and Phe-4 MHO98 . Configuration 1 corresponds to the global-minimum-energy state and configuration 2 to the second lowest-energy state. The distance between the two configurations is $`d=6.62`$ and the values of the potential energy (ECEPP/2 ECEP3 ) for configuration 1 and configuration 2 are $`10.72`$ kcal/mol and $`8.42`$ kcal/mol, respectively.

We analyze the free-energy landscape HOO99  from the results of our multi-overlap simulation at 300 K that performs a random walk between configurations 1 and 2. We study the landscape with respect to some reaction coordinates (and hence it should be called the potential of mean force). In order to study the transition states between reference configurations 1 and 2, we first plotted the free-energy landscape with respect to the dihedral distances $`d_1`$ and $`d_2`$ of Eq. (17). However, we did not observe any transition saddle point. A satisfactory analysis of the saddle point becomes possible when the root-mean-square (rms) distance (instead of the dihedral distance) is used. Figure 8 shows contour lines of the free energy reweighted to $`T=250`$ K, which is close to the folding temperature (HMO97 ; BNO03 ). Here, the free energy $`F(r_1,r_2)`$ is defined by

$$F(r_1,r_2)=k_BT\mathrm{ln}P(r_1,r_2),$$
(20)
where $`r_1`$ and $`r_2`$ are the rms distances from the reference configuration 1 and the reference configuration 2, respectively, and $`P(r_1,r_2)`$ is the (reweighted) probability at $`T=250`$ K to find the peptide with values $`r_1,r_2`$. The probability was calculated from the two-dimensional histogram of bin size 0.06 Å$`\times `$ 0.06 Å. The contour lines were plotted every $`2k_BT`$ ($`=0.99`$ kcal/mol for $`T=250`$ K).

Note that the reference configurations 1 and 2, which are respectively located at $`(r_1,r_2)=(0,4.95)`$ and $`(4.95,0)`$, are not local minima in free energy at the finite temperature ($`T=250`$ K) because of the entropy contributions. The corresponding local-minimum states at A<sub>1</sub> and B<sub>1</sub> still have the characteristics of the reference configurations in that they have backbone hydrogen bonds between Gly-2 and Met-5 and between Tyr-1 and Phe-4, respectively.

The transition state C in Fig. 8(a) should have intermediate structure between configurations 1 and 2. In Fig. 8(b) we show a typical backbone structure of this transition state. We see the backbone hydrogen bond between Gly-2 and Phe-4. This is precisely the expected intermediate structure between configurations 1 and 2, because going from configuration 1 to configuration 2 we can follow the backbone hydrogen-bond rearrangements: The hydrogen bond between Gly-2 and Met-5 of configuration 1 is broken, Gly-2 forms a hydrogen bond with Phe-4 (the transition state), this new hydrogen bond is broken, and finally Phe-4 forms a hydrogen bond with Tyr-1 (configuration 2).

In Ref. MHO98  the low-energy conformations of Met-enkephalin were studied in detail and they were classified into several groups of similar structures based on the pattern of backbone hydorgen bonds. It was found there that below $`T=300`$ K there are two dominant groups, which correspond to configurations 1 and 2 in the present article. Although much less conspicuous, the third most populated structure is indeed the group that is identified to be the transition state in the present work.

In table 1 we list the numerical values of the free energy, internal energy, and entropy multiplied by temperature at the two local-minimum states (A<sub>1</sub> and B<sub>1</sub> in Fig. 8(a)) and the transition state (C in Fig. 8(a)). Here, the internal energy $`U`$ is defined by the (reweighted) average ECEPP/2 ECEP3  potential energy:

$$U(r_1,r_2)=<E(r_1,r_2)>_T.$$
(21)
The entropy $`S`$ was then calculated by

$$S(r_1,r_2)=\frac{1}{T}\left[U(r_1,r_2)F(r_1,r_2)\right].$$
(22)
The free energy was normalized so that the value at A<sub>1</sub> is zero.

The state A<sub>1</sub> can be considered to be “deformed” configuration 1, and B<sub>1</sub> deformed configuration 2 due to the entropy effects, whereas C is the transition state between A<sub>1</sub> and B<sub>1</sub>. Among these three points, the free energy $`F`$ and the internal energy $`U`$ are the lowest at A<sub>1</sub>, while the entropy contribution $`TS`$ is the lowest at C. The free energy difference $`\mathrm{\Delta }F`$, internal energy difference $`\mathrm{\Delta }U`$, and entropy contribution difference $`T\mathrm{\Delta }S`$ are 1.0 kcal/mol, 1.9 kcal/mol, and $`0.9`$ kcal/mol between B<sub>1</sub> and A<sub>1</sub>, 2.2 kcal/mol, 4.6 kcal/mol, and $`2.4`$ kcal/mol between C and A<sub>1</sub>, and 1.2 kcal/mol, 2.7 kcal/mol, and $`1.5`$ kcal/mol between C and B<sub>1</sub>. Hence, the internal energy contribution and the entropy contribution to free energy are opposite in sign and the magnitude of the former is roughly twice as that of the latter at this temperature.

## IV Conclusions

In this article we have described the formulations of the two well-known generalized-ensemble algorithms, namely, multicanonical algorithm (MUCA) and replica-exchange method (REM). We then introduced four new generalized-ensemble algorithms as further extensions of the above two methods, which we refer to as replica-exchange multicanonical algorithm (REMUCA), multicanonical replica-exchange method (MUCAREM), multibaric-multithermal algorithm, and multi-overlap algorithm.

With these new methods available, we believe that we now have working simulation algorithms for spin systems and biomolecular systems.

Acknowledgments
This work is supported, in part, by NAREGI Nanoscience Project, Ministry of Education, Culture, Sports, Science and Technology, Japan.


text:
# On the origin of the galaxy luminosity function

## 1 Introduction

The Cold Dark Matter theory of structure formation has enjoyed notable successes in recent years. Not only does it provide a unified interpretation of the cosmic microwave background (CMB) and galaxy clustering, but on smaller scales it has led to successful models of the Ly$`\alpha `$ forest and gravitational lensing by galaxies, both strong and weak. While controversy surrounds the compatibility of the CDM theory with the measured dark-matter densities deep within galaxies, it seems likely that problems in this area arise from our poor understanding of galaxy formation rather than pointing to genuine shortcomings in the CDM model (Binney 2003).

Since the CDM theory is a theory of the invisible, tests of the theory rely to a large extent on predictions of how dark matter affects the baryon content of the Universe. On large scales forces other than gravity can generally be neglected, so the physics is simple and reliable. On the smallest scales both strong and electromagnetic interactions between baryons are as important as gravitational interactions, so the physics is complex. We know that this physics somehow gives rise to stars and galaxies, but the details are obscure.

Given the importance of galaxy-formation theory for making observational predictions from the CDM model, strenuous efforts have been made to model galaxy formation in the presence of CDM. These efforts are strongly influenced the papers of Rees & Ostriker (1977) and White & Rees (1978), which took it as axiomatic that gas that falls into a gravitational potential well immediately shocks to the well’s virial temperature. Galaxies were assumed to form as this rapidly heated gas slowly cooled. I argue here that this picture is wrong. In reality galaxies are formed from the the significant fraction of gas that fails to heat as it falls into the potential well, and that the fraction that immediately heats to the virial temperature is unable to cool subsequently. As the clustering hierarchy of the CDM model develops, the fraction of gas that heats to the virial temperature increases, so the efficiency of galaxy formation falls. The argument draws on recent developments in the study of ‘cooling flows’, simulations of baryonic infall (Birnboim & Dekel 2003; Katz et al. 2003), and semi-analytic galaxy-formation theory (Benson et al. 2003).

## 2 Does gas virialize?

The influential papers of galaxy-formation theory by Rees & Ostriker (1977) and White & Rees (1978) considered it evident that when gas falls into a potential well, it would be shock heated to the well’s virial temperature. In my thesis work (Binney 1977) I studied the shock heating of infalling gas in circumstances that maximized the chances of the gas shock heating: inspired by Zel’dovich’s pancake theory of galaxy formation I considered the accumulation of gas in a sheet of uniform surface density following the collapse of an initially homogeneous spheroidal body of gas. Gas in the equatorial plane was stationary at all times but increased in density as gas fell onto it from the polar directions with steadily increasing speed. As the speed of the infalling gas increased, the post-shock temperature of the infallen gas rose. In analytic work I compared the cooling time of the post-shock gas with the time required for infalling gas to travel a distance equal to the distance of the shock from the equatorial plane. This comparison showed that if the cooling rate was dominated by bremsstrahlung, the shock would break away from the equatorial plane only if the radius of the collapsing ellipsoid exceeded $`30\mathrm{kpc}`$. I used one-dimensional hydrodynamical simulations to test the analytical work and found that it seriously understimated the importance of cooling because (i) the infall velocity only gradually rises to the characteristic virial value, and (ii) even in a pure hydrogen-helium plasma, cooling at low temperatures greatly exceeds the bremsstrahlung rate. In a galaxy-sized collapse bremsstrahlung never dominates the cooling rate, and even for a cluster-sized collapse bremsstrahlung is important only when the last mass fraction falls in. Consequently, the fraction of gas that heats to the virial temperature is negligible in a galaxy-sized collapse and much less than unity in a cluster-sized one.

I considered the impact on the heated gas fraction of inhomogeneities in the infalling gas, and concluded that they were likely to favour cooling by increasing the surface density at which gas impacts the equatorial plane. Hence I concluded that essentially no gas would heat to the virial temperature when a galaxy-sized halo collapses, in flat contradiction to the fundamental premise of the influential paper of Rees & Ostriker (1977).

When CDM replaced neutrinos as the favoured dark matter component in the 1980s, it became necessary to revisit my 1977 study because CDM gives rise to smaller and cuspier potential wells, favouring cooling, while it enables one to lower the baryon density, which favours heating. Demoralized by the scant attention paid to my 1977 study I did not do so, but recently Birnboim & Dekel (2003) have done something similar. They have performed one-dimensional hydrodynamical simulations of baryonic infall into a spherical potential well, and concluded that “most galactic haloes that have collapsed and virialized by $`z2`$ did not produce a virial shock. Haloes less massive that $`10^{11}\mathrm{M}_{}`$ never produce a shock even if the gas has zero metallicity. If the metallicity is non-negligible (e.g., $`Z0.05`$) this lower bound to shock formation rises to $`7\times 10^{11}\mathrm{M}_{}`$.” It seems that the bottom line has changed very little in 25 years.

Further support for my contention that infall is ineffective in heating gas to the virial temperature is provided by Katz et al. (2003). They used smooth-particle hydrodynamics (SPH) to follow the evolution of a mixture of CDM and baryons in cosmological clustering simulations. They measured the maximum temperature $`T_{\mathrm{max}}`$ reached by each gas particle within a given halo, and found that at infall redshifts $`z\genfrac{}{}{0pt}{}{_>}{^{}}1`$ the distribution of $`T_{\mathrm{max}}`$ is distinctly bimodal, with about equal masses associated with peaks centred on $`3\times 10^4\mathrm{K}`$ and the virial temperature. At smaller values of $`z`$, the gap in a histogram of $`T_{\mathrm{max}}`$ values fills in, but the distribution of $`T_{\mathrm{max}}`$ values remains extremely broad, so a great deal of matter continues to fall in without heating to near the virial temperature. The quantity of infalling gas diminishes strongly over time, so most gas falls in when the temperature distribution is strongly bimodal.

Future increases in the spatial resolution of the simulations is likely to increase the fraction of gas that falls in cold since increasing the resolution raises the density of the densest, fastest-cooling gas and thus makes it harder for gas to heat to the virial temperature.

## 3 Feedback and the galaxy luminosity function

The clustering of CDM has been exstensively simulated and the popular ‘semi-analytic model’ of galaxy formation has been used to associate with each dark-matter halo a galaxy of given luminosity and colours (Kauffmann, White & Guiderdoni, 1993; Lacey et al. 1993; Cole et al. 1994; Kauffmann et al. 1999; Somerville & Primack 1999; Cole et al. 2000; Benson et al. 2003). The algorithms employed in semi-analytic models start by assuming that each halo’s quota of baryons shock heats to the virial temperature and then cools on the timescale that follows from the density of the gas. As gas cools it forms a cold disk within which stars form at some rate, so the mass and age distribution of the stellar population of any halo is determined, and its luminosity and colours can be predicted. Cole et al. (2000) provide a lucid account of this very detailed work.

The mass function of dark-matter haloes $`\mathrm{\Psi }(M)=\mathrm{d}N/\mathrm{d}M`$ is predicted by simulations of dark-matter clustering (Jenkins et al. 2001). It is quite different from the galaxy luminosity function $`\mathrm{\Phi }(L)=\mathrm{d}N/\mathrm{d}L`$ in two respects: (i) while at the low-$`M`$ end $`\mathrm{\Psi }(M)`$ is rather close to a power law $`\mathrm{\Psi }(M)M^\beta `$ with $`\beta 2.17`$, the luminosity function of galaxies can be fitted by Schechter’s (1976) formula $`\mathrm{\Phi }(L)L^\alpha \mathrm{exp}(L/L^{})`$ with a flatter slope: $`\alpha 0.95`$ (Cole et al. 2001; Kochanek et al. 2001). (ii) While both $`\mathrm{\Psi }`$ and $`\mathrm{\Phi }`$ fall below these power laws at large values of their argument, the luminosity $`L^{}`$ of the break in $`\mathrm{\Phi }(L)`$ corresponds to galactic masses that are $`\genfrac{}{}{0pt}{}{_>}{^{}}100`$ times smaller than the mass of the break in $`\mathrm{\Psi }`$. Thus, if one predicted $`\mathrm{\Phi }(L)`$ by converting $`\mathrm{\Psi }(M)`$ using the mass-to-light ratio of a typical $`L^{}`$ galaxy, one would predict too many galaxies at both extremes of the luminosity scale.

Two mechanisms provide plausible explanations of the relative dearth of low-luminosity galaxies: upon reionization at $`z_{\mathrm{ion}}6`$, the temperature of the baryons would have risen to $`T2\times 10^4\mathrm{K}`$ even before any shock heating. Gas at this temperature cannot be effectively trapped in a potential well with a peak circular velocity $`v_{\mathrm{max}}`$ smaller than $`35\mathrm{km}\mathrm{s}^1`$ (Efstathiou 1992; Quinn Katz & Efstathiou 1996; Thoul & Weinberg 1996). Since only the more massive haloes have larger circular speeds, the only low-mass haloes with stars will be those that formed before $`z_{\mathrm{ion}}`$. This prediction chimes nicely with the observation that dwarf galaxies are overabundant in clusters of galaxies because on the average cluster galaxies formed before comparable field galaxies.

Another mechanism capable of explaining the dearth of low luminosity galaxies is feedback from supernovae: a burst of star formation is expected to initiate the formation of a galaxy. The massive stars in this burst will complete their lives and explode as core-collapse supernovae on a timescale that is short compared to the dynamical time of the embedding halo. The supernova blasts will heat surrounding gas to temperatures $`\genfrac{}{}{0pt}{}{_>}{^{}}10^6\mathrm{K}`$. If $`v_{\mathrm{max}}`$ is smaller than $`100\mathrm{km}\mathrm{s}^1`$, the hot gas will stream outwards and further star formation may be effectively suppressed (Dekel & Silk 1986).

Benson et al (2003) have recently re-examined the fit to the galaxy luminosity function that one obtains by applying the semi-analytic model to the known mass function of dark-matter haloes. They conclude that with the relatively large baryon density that is now mandated by studies of large-scale structure, especially the CMB, the semi-analytic model cannot simultaneously fit both the low- and the high-luminosity end of $`\mathrm{\Phi }(L)`$. The strength of feedback by supernovae is poorly constrained and may be treated as a free parameter. When strong feedback is chosen, the density of low-luminosity galaxies can be made sufficiently small. However, the baryons expelled from low-luminosity systems later fall into haloes of larger mass and make them more luminous. The consequence is that there are then too many high-luminosity galaxies.

## 4 Cooling flows

We observe systems that contain gas at the virial temperature: X-ray observations long ago showed that nearly all clusters of galaxies and most luminous elliptical galaxies have such gas. Thus we are in a position to check observationally an important ingredient of the semi-analytic model of galaxy formation.

The cooling time of the gas in galaxy clusters and elliptical galaxies decreases as one approaches the centre, and in most systems the central value of the cooling time is much shorter than the Hubble time. This observation led to the expectation that gas in these systems was steadily forming stars as it cooled (Cowie & Binney 1977; Fabian & Nulsen 1977). From the radial surface-brightness profile of the X-rays and the assumption that the ‘cooling flow’ had reached a steady state, one can derive the rate of star formation as a function of radius (Nulsen 1986). Massive-star formation is certainly not taking place at the rate predicted by a Salpeter-like initial mass function (IMF), and formation of predominantly low-mass stars could also be excluded (Prestwich et al. 1997). Nor was cold gas present in the quantity required if star formation was somehow suppressed (Donahue et al. 2000; Edge 2001; Edge et al. 2002). Data from the Chandra and XMM-Newton satellites has recently shown that very little if any gas is cooling in clusters of galaxies, notwithstanding their short cooling times (Tamura et al. 2001; Peterson et al. 2002). Appropriate data are not yet available for individual elliptical galaxies, but the presumption must be that their smaller ‘cooling flows’ behave in a similar way. It is evident that the gas is being heated.

The details of how the gas is heated are controversial and irrelevant to the present discussion. What matters is that (a) very little if any gas is cooling to temperatures lower than $`1/3`$ that of the main body of gas, (b) any cold gas is confined to the very centre of the system (Edge et al. 2002), and (c) this gas is dusty and is quite likely material that has fallen in cold and is unconnected with the hot gas (Sparks et al. 1989).

Another important implication of the study of extended X-ray sources is that the gas observed in these systems is not heated to the virial temperature by gravitational shocking, but by supernovae and/or an active galactic nucleus (AGN). This fact was first recognized by Kaiser (1991), who pointed out that in the contrary case not only would the number of bright X-ray sources at $`z1`$ be greater than observed, but the correlation between luminosity and temperature at the present epoch would be less steep than is observed (Kaiser 1986). Kaiser’s inference has been confirmed in many subsequent studies (e.g., Shimizu et al. 2003). Direct evidence for non-gravitational heating is provided by the large fraction ($`\genfrac{}{}{0pt}{}{_>}{^{}}50`$ percent) of the heavy elements in clusters of galaxies which are contained in the X-ray emitting gas but must have been synthesized within the galaxies.

## 5 Calling a halt to star formation

The study of ‘cooling flows’ teaches us that systems in which all gas is at the virial temperature do not form stars and do not accumulate cold gas from which stars could form at a later date. The evidence for cold infall suggests that star formation in such a system could be restarted if a quantity of cold gas fell in. An episode of cold infall will not necessarily lead to star formation, however, because a sufficiently dense hot atmosphere will ablate and reheat infalling gas before it can form stars. In fact, we can assume that cold gas enters the virialized system at a density that does not permit rapid star formation, since otherwise it would already have fragmented into stars. Hence to form stars it has to accumulate in the halo potential, probably on some approximately closed orbits. The extended HI disks that are frequently observed around late-type galaxies are regions where such accumulation of cold infalling material is taking place. A sufficiently dense, hot atmosphere will inhibit this accumulation by ablating and stripping the cold gas, both (i) as it falls towards the closed parking orbits, and (ii) within those orbits. It will be highly vulnerable to ablation in phase (i) because its surface density will be low, but stage (ii) can last longer. A credible ab initio calculation of the ablation rate is hard because it would require a quantitative understanding of the magnetized, turbulent boundary layer between the cold and hot gas, which differ by at least three orders of magnitude in density.

## 6 Putting it all together

The following scheme seems to account neatly for all the evidence summarized above. When a dark halo virializes, a fraction $`f`$ of the progenitor’s gas remains at $`T<10^4\mathrm{K}`$, while the remainder heats to the virial temperature. The fraction $`f`$ decreases as the halo mass increases. If the halo is massive enough, the fraction $`f`$ is retained as a disk in which star formation commences. The gravitationally heated fraction is further heated by dying stars and by any massive black hole that forms at the system’s nucleus. In a lower-mass halo this gas expands and flows out of the system. As we move up the sequence to larger masses, the temperature that must be reached to push the hot gas out increases, and at some mass $`M^{}`$, which we must identify with that of $`L^{}`$ galaxies, the gas cannot be pushed completely out of the halo, although it can be pushed so high into the halo that it takes several gigayears to fall back once the rate of heating by core-collapse supernovae has moderated (D’Ercole et al 1989). Once the hot gas cannot escape, its density steadily increases as dying stars inject fresh hot gas. At some point it becomes dense enough to ablate and absorb any infalling cold material before it can give rise to significant star formation.

The density of virial-temperature gas increases rapidly with halo mass near $`M^{}`$ because both the fraction $`1f`$ and the ability to retain supernova-heated gas increase with $`M`$.

Current semi-analytic models of galaxy formation form too many high-luminosity galaxies when appropriate feedback is used because they assume that gas heated to the virial temperature can cool and make stars. Once this error is corrected, the models will yield luminosity functions that cut off sharply above $`L^{}`$. Instead the revised models will predict that high-mass haloes will contain large masses of X-ray emitting virialized gas, in agreement with observation.

## 7 Conclusions

Observations of cooling flows lead to the conclusion that virial-temperature gas cannot cool and give rise to star formation. The inability of hot gas to cool when supernova heating is effective is obvious. The continued inability of this gas to cool once the supernova heating fades has surprised many, but can be readily explained and was predicted (Tabor & Binney 1993; Binney & Tabor 1995). The key point is that any very cold gas that forms, does so in the immediate viscinity of the central black hole, and induces accretion that promptly reheats the gas, largely shutting off further accretion. The details of this chain of events are inadequately understood, but the prospects for elucidating them are good (Reynolds et al, 2001; Omma et al 03).

Since gas heated to the virial temperature cannot cool and form stars, galaxy formation would be impossible if Rees & Ostriker were correct in assuming that gas is heats to the virial temperature when a system virializes. Three different simulations argue to the contrary that a significant fraction of the gas never reaches the virial temperature. Galaxies form from this under-heated fraction. In lower-mass systems the heated fraction is largely ejected. Some of this subsequently finds its way into more massive haloes. These haloes heat a larger fraction of their progenitor’s gas to the virial temperature and tend to retain both this gas and gas ejected by their stellar systems. Consequently, they develop dense atmospheres of virial-temperature gas. These atmospheres ablate any cold gas that falls into their haloes, so star formation ceases in these systems. Consequently, high-mass haloes are now associated with only modest luminosities from rather old stars and substantial masses of X-ray emitting virial-temperature gas.


text:
# The magnetic reversal in dot arrays recognized by the self-organized adaptive neural network

## 1 Introduction

Many interesting applications of the nanometer periodic magnetic dot arrays are expected in the magnetic recording  and magnetic sensors. The magnetic nanostructures represent model systems for study of interesting physical phenomena stemming from the well controlled dot shape anisotropy and dominance of magnetostatic interactions at the scales comparable or larger than dot size. The magnetostatics cause the coupling of domains over the inter-dot interfaces. The current perspective technologies require understanding and controlling of dynamic magnetization processes in fine particles on nanoscales, because by them the magnetic bits are written into recording media. Hence of particular interest is the dynamics of magnetization reversal.

The static and dynamic nonuniform magnetization states in magnetic nanostructures including the thin films and isolated small particles have been studied by micromagnetic  and Monte Carlo  simulations. On a qualitative level intra-dot nonuniformities can be characterized as a superpositions of vortices, flowers, domain walls, etc.. The problem of nonuniformities in exchange-dominated structures can be solved exactly .

Our recent activity is directed to algorithms allowing systematic and fast classification of simulated magnetic structures. The methodology presented in this paper is inspired by the progress in the theory of artificial neural networks that are nonlinear models suitable to reduce, classify or interpolate data structures representing the complex patterns . In the previous paper  we suggested the implementation of radial basis function networks to model the ordering of dots in arrays. In the next paper  the neural network training was applied to multiscale numerical computations. The lack of our previous two approaches was the absence of realistic dynamics and restriction to quasi-static limit. In the present paper our interest is focused to the networks which analyze spin configurations generated by the spin lattice dynamical simulations. It means that network solves here classical post-processing problem known in the micromagnetic context . Our approach has many in common with experiment treatment , where self-organized network was employed to identification of vortex modes of magnetic thin film media. In the presented paper the classification is applied to configurations generated by the magnetization reversal within the XY superlattice model of ultra-dense array. The motivation for study of ultra-dense spin structures was the expectation of extremal inter-dot effects.

The main reason for the choice of the XY spin model is relative computational simplicity of 2-D dynamics. But there are two regimes for which 2-D approximation is justified physically. (i) The first regime belongs to limit of very strong damping which suppresses precession and accelerates reversal . (ii) The second regime represents limit of strong in-plain  or shape anisotropy associated with demagnetization field of thin film sample . In this regime anisotropy causes deformation of spin trajectory and accelerates reversal . Nevertheless, one should to note that realistic damping is sufficiently small, and thus 2-D remagnetization is practically quasistatic. In the quasistatic limit the damping parameter has meaning of auxiliary parameter rescaling the time. Recent results obtained for discrete spin systems provide more optimistic view on the valuation of artificial dynamics of reversal. The example is the simulation of  demonstrating that switching time Monte Carlo dependences can coincide after some rescaling with predictions of the stochastic Landau-Lifshitz equation.

The plan of the paper is the following. The interaction picture and dynamics of array model are defined in section 2. The simulations of array reversal are presented and discussed in section 3. The details of adaptive neural network treatment of intra-dot configurations are given in subsection 3.3.

## 2 XY periodically modulated superlattice model of array

In this section we formulate the interaction picture and dynamics of ultra dense array model. The magnetic ordering is described by the systems of the classical planar spins located on the square two dimensional lattice. The spin on site $`\alpha `$ is affected by effective field

$$𝐡_\alpha ^{\mathrm{eff}}=\frac{D}{J}\underset{\genfrac{}{}{0pt}{}{\beta =0}{\beta \alpha }}{\overset{L^21}{}}\frac{3𝐧_{\alpha \beta }(𝐧_{\alpha \beta }𝐒_\beta )𝐒_\beta }{|𝐫_{\alpha \beta }|^3}+\underset{\genfrac{}{}{0pt}{}{\beta =0}{|𝐫_{\alpha \beta }|=1}}{\overset{L^21}{}}ϵ_{\alpha \beta }𝐒_\beta +\frac{𝐇}{J},$$
(1)
where $`𝐒_\alpha =S_{\alpha \mathrm{x}}𝐞_\mathrm{x}+S_{\alpha \mathrm{y}}𝐞_\mathrm{y}`$ is the planar spin, $`𝐇`$ is the external magnetic field in energy units, $`L`$ is the linear size of the square $`L\times L`$ spin lattice system. Here and in further $`𝐞_\mathrm{x}`$, $`𝐞_\mathrm{y}`$ denote the unit vectors along the main Cartesian axes. Rather convenient for a computer purposes is to use one dimensional enumeration scheme of spin sites $`\alpha \{0,1,2,\mathrm{}L^21\}`$. (All of the geometric details presented in this section are discussed in the Appendix A.1.) The first sum of Eq.(1) represents dipolar contribution proportional to constant $`D`$, where $`𝐫_{\alpha \beta }`$ is vector pointing from site $`\alpha `$ to $`\beta `$. The vector is measured in lattice constant units, $`𝐧_{\alpha \beta }=𝐫_{\alpha \beta }/|𝐫_{\alpha \beta }|`$ is the auxiliary unit vector. The second sum of Eq.(1) involves nearest neighbor exchange coupling term $`Jϵ_{\alpha \beta }`$ determined by single parameter $`J>0`$ multiplying the site-site connectivity matrix $`ϵ_{\alpha \beta }\{1,0\}`$ constrained by $`|𝐫_{\alpha \beta }|=1`$ condition; $`ϵ_{\alpha \beta }`$ is $`0`$ for any bond at inter-dot boundary and $`1`$ inside a dot. The site dependence of $`ϵ_{\alpha \beta }`$ takes into account square shape of dot of the linear size $`L_0`$. The spin system partitioned by $`ϵ_{\alpha \beta }`$ consists of $`L_\mathrm{d}^2=(L/L_0)^2`$ dots with minimum inter-dot distance equal to the unit lattice spacing and corresponding to extremal space filling.

In 2-D case Landau-Lifshitz equation reduces to damping term. The numerical integration scheme discussed in Appendix A.2 can be derived from recursive propagator form

$$𝐒_\alpha (\tau +\mathrm{\Delta }\tau )=\widehat{𝐔}\left(_\tau ^{\tau +\mathrm{\Delta }\tau }\omega _\alpha (\tau ^{})\text{d}\tau ^{}\right)𝐒_\alpha (\tau ),$$
(2)
where $`\omega _\alpha (\tau )=(𝐒_\alpha (\tau )𝐞_\mathrm{x})(𝐡_\alpha ^{\mathrm{eff}}(\tau )𝐞_\mathrm{y})(𝐒_\alpha (\tau )𝐞_\mathrm{y})(𝐡_\alpha ^{\mathrm{eff}}(\tau )𝐞_\mathrm{x})`$ is instant angular frequency of the spin rotation, $`\mathrm{\Delta }\tau `$ is the time integration step, $`\tau `$ is the reduced time related to real time $`t`$ as $`\tau =\mu _0\gamma M_\mathrm{s}t`$, where $`\mu _0`$ is permeability of free space, $`\gamma `$ is the gyromagnetic ratio and $`M_\mathrm{s}`$ is the saturation magnetization. The unitarity of spin rotation matrix

$$\widehat{𝐔}(\phi )=\left(\begin{array}{cc}\mathrm{cos}\phi \hfill & \hfill \mathrm{sin}\phi \\ \mathrm{sin}\phi \hfill & \hfill \mathrm{cos}\phi \end{array}\right)$$
(3)
implies that numerical scheme adapted from Eq.(2) conserves the spin vector size independently of $`_\tau ^{\tau +\mathrm{\Delta }\tau }\omega _\alpha (\tau ^{})\text{d}\tau ^{}`$ quadrature.

The time demanding task represents summation of dipolar fields. The substantial reduction of time allows the use of the hierarchical approach  with a detailed summation of near and rough summation of far field contributions. The specific geometry of array leads to the natural spin block definition

$$𝐒_j^\mathrm{d}=\underset{\genfrac{}{}{0pt}{}{\beta =0}{i(\beta )=j}}{\overset{L^21}{}}𝐒_\beta ,j\{0,1,2,\mathrm{},L_\mathrm{d}^21\},$$
(4)
where sum runs over the spins of index $`\beta `$ belonging to dot index $`i(\beta )=j`$ and dot position $`𝐫_{i(\beta )}^\mathrm{d}`$. In frame of two-level hierarchical approximation dipolar field is splitted to near $`𝐡_\alpha ^{\mathrm{eff},\mathrm{N}}`$ and far $`𝐡_\alpha ^{\mathrm{eff},\mathrm{F}}`$ field

$$𝐡_\alpha ^{\mathrm{eff},\mathrm{N}}=\frac{D}{J}\underset{\genfrac{}{}{0pt}{}{\beta =0,\beta \alpha }{|𝐫_{i(\alpha ),i(\beta )}^\mathrm{d}|r_{\mathrm{cut}}}}{\overset{L^21}{}}\frac{3𝐧_{\alpha \beta }(𝐧_{\alpha \beta }𝐒_\beta )𝐒_\beta }{|𝐫_{\alpha \beta }|^3}+\underset{\genfrac{}{}{0pt}{}{\beta =0}{|𝐫_{\alpha \beta }|=1}}{\overset{L^21}{}}ϵ_{\alpha \beta }𝐒_\beta +\frac{𝐇}{J},$$
(5)
$$𝐡_\alpha ^{\mathrm{eff},\mathrm{F}}=\frac{D}{J}\underset{\genfrac{}{}{0pt}{}{j=0}{|𝐫_{i(\alpha )j}^\mathrm{d}|>r_{\mathrm{cut}}}}{\overset{L_\mathrm{d}^21}{}}\frac{3𝐧_{i(\alpha )j}^\mathrm{d}(𝐧_{i(\alpha )j}^\mathrm{d}𝐒_j^\mathrm{d})𝐒_j^\mathrm{d}}{|𝐫_{i(\alpha )j}^\mathrm{d}|^3},$$
(6)
depending on the unit vectors $`𝐧_{i(\alpha )j}^\mathrm{d}=𝐫_{i(\alpha )j}^\mathrm{d}/|𝐫_{i(\alpha )j}^\mathrm{d}|`$, $`𝐫_{i(\alpha )j}^\mathrm{d}=𝐫_j^\mathrm{d}𝐫_{i(\alpha )}^\mathrm{d}`$. The summation rule is controlled by the scale-separation cut-off $`r_{\mathrm{cut}}`$.

## 3 The results of array reversal simulations

In this section we present the simulations of the square shaped finite array including $`40000`$ spins. To analyze the simulation outputs we applied: (i) the qualitative recognition of the configuration snapshots; (ii) the comparing of the complex array dynamics with the dynamics of much simpler systems; (iii) the clustering of intra-dot magnetic configurations using neural networks. The global information on the array ordering concerns the magnetization per site

$$𝐌(\tau )=\frac{1}{L^2}\underset{\alpha =0}{\overset{L^21}{}}𝐒_\alpha (\tau ).$$
(7)
The simulations were performed for uniform zero $`𝐒_\alpha |_{\tau =0}=𝐇/H=𝐞_\mathrm{x}`$ and non-zero bias $`𝐒_\alpha |_{\tau =0}=\mathrm{cos}(\phi _0)𝐞_\mathrm{x}+\mathrm{sin}(\phi _0)𝐞_\mathrm{y}`$, $`\phi _00`$ initial conditions. In both cases the characteristic time of array is the time when $`𝐇𝐌`$ projection gains at first the bound

$$𝐇𝐌(\tau _{\mathrm{ch}})=0.$$
(8)
### 3.1 The qualitative approach based on two-scale similarity

To obtain the preliminary qualitative information about the switching mechanism we run the simulation of reversal dynamics in the field $`H/J=1`$ and coupling $`D/J=0.1`$ for dot size $`L_0=20`$ and bias $`\phi _0=0`$. Under these conditions we obtained estimate $`\tau _{\mathrm{ch}}=8.516`$.

At the beginning of reversal the external magnetic field nucleates reversal domains at the corners of dots. The recognition of simulated systems shows that remagnetization is more homogeneous and easier when dots are localized at the corners of array, whereas nonuniformities occur due to mirror symmetry of effective fields typical for dots standing in positions closer to axial zone of array. This property is illustrated by snapshot taken at characteristic time as is shown in Fig. 1. The typical for nonuniform remagnetization is that reverse domains, which gain external field polarity, grow from the dot corners and pass toward the dot centers. At the intermediate stages domain walls - nonlinear waves  between reverse and initial polarity domains become narrow. In the context with nanoparticles and our small-scale simulations we prefer to use the term nonlinear wave instead of domain wall since nonlinear wave better expresses nonlocal character of arising spin structures. At the final stage of the inhomogeneous reversal the sum of external and dipolar field can, if sufficiently strong, homogenize intra-dot ordering by annihilating of wave/anti-wave pairs (WAP). Similar mechanism was studied in detail for linear spin chain . As we see from the scheme depicted in Fig. 2 the reversal mechanism in dots is position and time dependent for finite size arrays. The formation of WAP is an example of inhomogeneous remagnetization. The process is characteristic for later times and dots placed nearly the axial array zone as is depicted in Fig. 3. The parallelism between spin chain and array behavior inspired the monitoring of spin projections summarized in Fig. 4. The figure shows $`S_{\alpha y}`$ components taken from several superlattice site lines crossing array in external field direction. The time sequence of these projections reveals that each inter-dot interface acts as a phase-shifter of WAP formed in different dots. The consequence of inter-dot coupling is the large-scale WAP with shape resembling its small-scale intra-dot counterpart (see Fig. 5). It is rather surprising that for its build up magnetostatic forces suffice. The spatio-temporal two-scale similarity of magnetic ordering illustrated in Fig. 6 originates from spatial distribution of dipolar fields of square-shaped dots embedded into square-shaped array.

### 3.2 The quantitative analysis of switching, the effect of bias

#### 3.2.1 The dot-size dependences

We continue with the comparing of $`\tau _{\mathrm{ch}}(L_0)`$ dependences for patterned structure as well as for the isolated dots with zero and non-zero bias. The varying of $`L_0`$ is in some sense equivalent to varying of dipolar and exchange coupling contributions. The results for patterning $`L_0=5,10,20,25,40,50`$ for $`D/J=1`$ are depicted in Fig. 7, where the characteristic times of the interacting and isolated dots are also compared. From the simulations it follows that inhomogeneities grow as $`L_0`$ increases. The reverse phase nucleates preferentially at dot surfaces and reduces $`\tau _{\mathrm{ch}}`$ in larger dots. If inter-dot interactions are considered $`\tau _{\mathrm{ch}}(L_0)`$ dependences look smoother than dependences obtained for isolated dots. The interactions affect spatial dispersion of $`\omega _\alpha `$, which means that some dots are accelerated and some decelerated during reversal. The simulations of isolated dots show that differences between zero and non-zero bias are larger for smaller $`L_0`$. The decrease of $`\tau _{\mathrm{ch}}`$ due to uniform bias $`\phi _0>0`$ is a manifestation of the Stoner-Wohlfarth coherent regime.

#### 3.2.2  The interplay of $`D`$, $`J`$, $`H`$, $`\phi _0`$ parameters, minimal four-spin model

In this section minimal four-spin model and complete array model of magnetization reversal are discussed and compared. The minimal model consists of four spins arranged to $`2\times 2`$ square segment of unit size. In the segment four nearest pairs are coupled by the exchange interaction and all of the pairs by dipole-dipole interaction. Despite of substantial geometric truncation, all of the interactions considered for array take place in four-spin model. The simulation of the minimal model system provides $`\tau _{\mathrm{ch}}^\mathrm{f}`$ characteristics depicted in Fig. 8. They consist of two qualitatively different regimes separated by bifurcation point $`D=D_\mathrm{c}^\mathrm{f}`$ manifested by $`\tau _{\mathrm{ch}}^\mathrm{f}(D)`$ cusp. At $`D_\mathrm{c}^\mathrm{f}`$ several modes compete, the interaction and external fields are nearly balanced and frustrated spins need quite long time to select the winning mode. For $`D<D_\mathrm{c}^\mathrm{f}`$ the reversal is very susceptible to choice of initial bias. It is easy to check that even small $`\phi _0>0`$ accelerates reversal and stabilizes coherent or quasi-coherent rotation in coincidence with the Stoner-Wohlfarth scenario. Further simplification of four-spin dynamics which clarifies situation can be done assuming the total coherence of spins. The consequence is the elimination of three angular variables and analytic expression

$$\tau _{\mathrm{ch}}^\mathrm{f}=\frac{4\sqrt{2}J}{}_{\phi _0}^{\pi /2}\frac{\mathrm{d}\phi }{3D\mathrm{cos}(2\phi )+4\sqrt{2}H\mathrm{sin}\phi }.$$
(9)
where $`=4\pi \lambda /\gamma \mu _0M_\mathrm{s}`$ is undimensional parameter proportional to damping parameter $`\lambda `$. The expression obtained in the limit $`D0`$ is rather well known

$$\tau _{\mathrm{ch}}=\frac{J}{H}_{\phi _0}^{\pi /2}\frac{\mathrm{d}\phi }{\mathrm{sin}\phi }=\frac{J}{H}\mathrm{ln}(\text{cotg}(\phi _0/2)).$$
(10)
For $`D>0`$ the integral Eq.(9) leads to finite $`\tau _{\mathrm{ch}}^\mathrm{f}`$ even in the case $`\phi _0=0`$. It means that dipolar coupling removes logarithmic divergence $`\mathrm{ln}\phi _0`$ of $`\tau _{\mathrm{ch}}^\mathrm{f}`$ and substitutes in part the effect of $`\phi _0>0`$. In the regime $`D<D_\mathrm{c}^\mathrm{f}`$ the bias $`\phi _0=0`$ induces non-uniform flower mode belonging to characteristic time which is out of the scope of Eq.(9). The four-spin dynamics attained by simulation for $`D>D_\mathrm{c}^\mathrm{f}`$ is even more complex. In that case the choice of bias losses its dynamical relevance. The simulations show that quite short time is needed to form the transient long-living flower mode (long with respect to $`\tau _{\mathrm{ch}}^\mathrm{f}`$) which alters to vortex mode at final stages.

With the four-spin system behavior in mind the attention is turned now to the complex array dynamics. As we see from Fig. 9 $`\tau _{\mathrm{ch}}(D/J)`$ dependence does not contain sharp cusp. Instead of it, the slope variation of $`\tau _{\mathrm{ch}}(D/J)`$ arises (labeled as $`D_\mathrm{c}`$ in analogy to $`D_\mathrm{c}^\mathrm{f}`$). For $`L_0=20`$ and $`H/J=1`$ the position of $`D_\mathrm{c}`$ can be estimated as $`D_\mathrm{c}5J`$. The configuration snapshots show that near $`D_\mathrm{c}`$ the array consists of spin configuration loops of diameters approaching single dot size or size of several dots. For $`D>D_\mathrm{c}`$ the array reversal generates non-collinear antiferromagnetic phases. The common feature of oversimplified four-spin model and array simulation is that finite $`\tau _{\mathrm{ch}}^\mathrm{f}`$ and $`\tau _{\mathrm{ch}}`$ are attained even for $`\phi _0=0`$ due to non-collinearity of external and dipolar fields. From dependences depicted in Figs. 7 and 9 it follows that geometric variations in $`L_0`$ can produce less remarkable changes than $`D/J`$ parameters.

For completeness the external field dependences of characteristic time should be mentioned. From Eq.(9) one can obtain for $`D/H1`$ the asymptotic expansion of characteristic time $`\tau _{\mathrm{ch}}^\mathrm{f}=(J/(H))`$ $`_{n=0}^{\mathrm{}}(1)^n(D/H)^nc_n(\phi _0)`$, $`c_n(\phi _0)=3^n2^{5n/2}_{\phi _0}^{\pi /2}d\phi \mathrm{cos}^n(2\phi )(\mathrm{sin}(\phi ))^{n1}`$ suitable for proposition of reasonable fitting function covering the array data. The array simulations in Fig. 10 shows that $`\tau _{\mathrm{ch}}`$ coincides qualitatively with asymptotic predictions of four-spin system.

### 3.3  The intra-dot complex structures in representation of ART network

The artificial neural network models  inspired by the physiology are able to mimic action of neurons and synaptic connections of brain. The salient feature of neural systems is the associative recognition of the complex structures. This ability is attained by training of synaptic connections - neurons. In this paper we deal with networks based on the adaptive resonance theory (ART) developed by Carpenter and Grossberg  originally as a model explaining adaptive phenomena in visual systems. The family of ART algorithms belongs to unsupervised self-organized algorithms motivated by need to construct network sufficiently adaptive to novel inputs. The adaptivity means that network produces new neurons if diversity of inputs overcomes the predefined threshold represented by so called vigilance parameter. On the other hand redundancy of information leads to the annihilation of neurons. For summary of ART network applications see e.g.

#### 3.3.1 The encoding of intra-dot magnetic structure

The ideas of clustering and adaptivity of the ART networks led us to their application to magnetic intra-dot configurations generated during reversal. At first, an appropriate format for encoding of magnetic configurations is suggested. Similarly to formulation  the configurations are encoded by $`2N_\mathrm{c}`$ dimensional tuples

$$\stackrel{~}{m}_i[𝐦_{i1},𝐦_{i2},\mathrm{},𝐦_{iN_\mathrm{c}}],$$
(11)
where

$$𝐦_{in}=\frac{1}{N_\mathrm{S}}\underset{l\mathrm{}_{in}}{}𝐒_l,i=1,2,\mathrm{},L_\mathrm{d}^2,n=1,2,\mathrm{},N_\mathrm{c}$$
(12)
is the local average taken from $`N_\mathrm{S}`$ spins (see Fig. 11), where sum runs over site $`l`$ taken from $`n`$th segment of $`i`$th dot labeled as $`\mathrm{}_{in}`$. Each dot is subdivided to $`N_\mathrm{c}=L_0^2/N_\mathrm{S}`$ intra-dot square segments. The system of $`L_\mathrm{d}^2`$ intra-dot configurations is transformed into $`\stackrel{~}{m}_i`$ tuples incoming to ART, which compresses them to $`N_\mathrm{w}L_\mathrm{d}^2`$ neurons

$$\stackrel{~}{w}_j[𝐰_{j1},𝐰_{j2},\mathrm{},𝐰_{jN_\mathrm{c}}],j\{1,2,\mathrm{}N_\mathrm{w}\}$$
(13)
representing the basal types of intra-dot ordering. The self-organization mechanism of neurons discussed in the next is controlled by the coincidence between input $`\stackrel{~}{m}_i`$ and output $`\stackrel{~}{w}_j`$. As a proper measure of coincidence is chosen Euclidean distance

$$\stackrel{~}{w}_j\stackrel{~}{m}_i=\sqrt{\underset{n=1}{\overset{N_\mathrm{c}}{}}(𝐰_{jn}𝐦_{in})(𝐰_{jn}𝐦_{in})^\mathrm{T}},$$
(14)
where $`\mathrm{T}`$ labels transposition.

#### 3.3.2 The ART clustering algorithm

The training algorithm is described in the following points

1. Initialization setting $`N_\mathrm{w}=1`$, $`\stackrel{~}{w}_1=\stackrel{~}{m}_i`$ for random dot selection between $`i\{0,1,\mathrm{},L_\mathrm{d}^21\}`$.
2. Loop through the set of intra-dot magnetic configurations. For each randomly selected dot $`i`$ follows:

1. Presentation of $`\stackrel{~}{m}_i`$ to ART network. The training is realized for $`k=1,2,\mathrm{}`$ iteration steps.
    2. Computing of the actual index $`j^{}(k,i)`$ of the fired neuron $`\stackrel{~}{w}_{j^{}(k,i)}(k)`$ trained by $`\stackrel{~}{m}_i`$ pattern according to the competitive rule

    $$\begin{array}{cccc}\hfill j^{}(k,i)=& \text{arg}\hfill & \text{min}& \stackrel{~}{w}_j(k)\stackrel{~}{m}_i\hfill \\ & & ^{j=1,2,\mathrm{},N_\mathrm{w}}& \end{array}.$$
    (15)
    3. Comparing $`\stackrel{~}{w}_{j^{}(k,i)}(k)\stackrel{~}{m}_i`$ to the vigilance parameter $`\rho `$.

    1. Update of the neuron weights via Hebbian training rule

        $$\stackrel{~}{w}_{j^{}(k,i)}(k+1)=\stackrel{~}{w}_{j^{}(k,i)}(k)+\eta (k)\left(\stackrel{~}{m}_i\stackrel{~}{w}_{j^{}(k,i)}(k)\right)$$
        (16)
        is applied if $`\stackrel{~}{w}_{j^{}(k,i)}\stackrel{~}{m}_i\rho `$. The update shifts winning vector toward nearby input with rate proportional to plasticity parameter $`\eta (k)=\eta _0\mathrm{exp}\left(k/\tau _{\mathrm{train}}\right)`$ relaxing with training time $`\tau _{\mathrm{train}}`$.
        2. Creation of new neuron $`N_\mathrm{w}N_\mathrm{w}+1`$, $`\stackrel{~}{w}_{_{N_\mathrm{w}}}\stackrel{~}{m}_i`$ if $`\stackrel{~}{w}_{j^{}(k,i)}(k)\stackrel{~}{m}_i>\rho `$.
    4. Annihilation of neuron pair $`[\stackrel{~}{w}_{z_1^{}}(k),\stackrel{~}{w}_{z_2^{}}(k)]`$, ordered as $`z_1^{}<z_2^{}`$, with indices $`z_1^{}`$, $`z_2^{}`$ selected according to minimum distance

    $$\begin{array}{cccc}\hfill [z_1^{},z_2^{}]=& \text{arg}\hfill & \text{min}& \stackrel{~}{w}_{z_1}\stackrel{~}{w}_{z_2}\hfill \\ & & ^{z_1,z_2}& \end{array}$$
    (17)
    if $`\stackrel{~}{w}_{z_1^{}}\stackrel{~}{w}_{z_2^{}}<\rho `$. The product of annihilation $`[z_1^{},z_2^{}]z_1^{}`$, $`N_\mathrm{w}N_\mathrm{w}1`$ is the neuron $`\stackrel{~}{w}_{z_1^{}}`$ determined by the midpoint rule $`\stackrel{~}{w}_{z_1^{}}\frac{1}{2}(\stackrel{~}{w}_{z_1^{}}+\stackrel{~}{w}_{z_2^{}})`$. After this update of $`\stackrel{~}{w}_{z_1^{}}`$ the neurons having $`jz_2^{}`$ are rewritten: $`\stackrel{~}{w}_{j1}\stackrel{~}{w}_j`$.
3. Stop criterion is represented by the inequality

$$\frac{1}{L_\mathrm{d}^2}\underset{i=0}{\overset{L_\mathrm{d}^21}{}}\stackrel{~}{w}_{j^{}(k,i)}(k+1)\stackrel{~}{w}_{j^{}(k,i)}(k)<\epsilon ,$$
(18)
where $`\epsilon `$ is small parameter. If the above inequality is not fulfilled, the algorithm follows from the step 2 with $`k`$ incremented by $`1`$.

The stop criterion indicates when the network attains a fixed point of $`\stackrel{~}{w}_j,j=1,2,\mathrm{},N_\mathrm{w}`$ including the cluster of typical representative intra-dot configurations.

#### 3.3.3 The ART applied to different stages of reversal

The ART network has been applied independently to array configurations attained at selected stages of reversal. The results of self-organization are presented in Figs. 13-15. The ART self-organization was carried out for tuned parameters $`\eta _0=0.1`$, $`\tau _{\mathrm{train}}=20`$, $`\epsilon =10^6`$, $`\rho =1.15`$. The most important seems to be optimized choice of an appropriate vigilance connecting $`\rho `$ and $`N_\mathrm{w}`$. It stem from calibration in Fig. 12. During it $`\rho `$ was decreased slowly from the initial $`\rho =1.5`$ stabilizing single neuron. This decrease of $`\rho `$ was stopped for $`N_\mathrm{w}=L_\mathrm{d}^2=100`$ competing neurons.

The ART view point on reversal leads to intra-dot taxonomy depicted in Fig. 13. In this figure the most similar neurons \[in the sense of generalized distance Eq.(14)\] at different times are joined by the arrows. The ART model clearly distinguishes between branches of coherently rotating monodomain dots and dots with more complex order. The diversity of configurations occurring near the internal energy maxima is reflected by the cardinality of neuron population. The ART pruning at late time stages is associated with small intra-dot diversity within the array. How choice of $`\rho `$ affects the information content of network is clearly demonstrated in Fig. 14. The supplementary problem solved by ART represents classification of interfacial inter-dot configurations. In that case the configuration tuples $`\stackrel{~}{m}_i`$, $`\stackrel{~}{w}_j`$ have been combined from $`𝐦_{in}`$ segments \[see Eq.(12)\] of adjacent quadrants of four neighboring dots. The results depicted in Fig. 15 indicate preferential formation of inter-dot vortices and strong tendency to inter-dot pinning also for conditions when formation of intra-dot vortices is disadvantageous.

## 4 Conclusions

The magnetization reversal has been simulated for superlattice spin model of the ultra-dense magnetic dot array. Several aspects of reversal dynamics for different exchange and magnetostatic couplings, external fields, geometric parameters and initial conditions have been studied. It was found that the principal physical consequence of geometric modulation is the formation of two-scale WAP. The ART viewpoint seems to be fruitful in thinking how dot switches from one configuration to another within the large time intervals, and how to identify dominant channels of intra-dot evolution. The present experience creates believe that advantages of ART paradigm should be exploited namely in 3-D lattices, where visual classification meets natural bounds.

## 5 Acknowledgments

The work was partially supported by the grant No.1/6020/99 (Slovak Grant Agency VEGA). M. Gmitra expresses thanks to Dr. Z. Frait and Dr. V. Kamberský from Institute of Physics AS CR, Prague, for their hospitality and valuable discussions during stay supported by training network contract HPRN-CT-1999-00150.

## Appendix A Appendix

### A.1 One dimensional enumeration scheme of superlattice sites

In the Appendix site enumeration scheme mentioned in section 2 is described. In this scheme the superlattice relative position vector $`𝐫_{\alpha \beta }`$ pointing from site $`\alpha `$ to site $`\beta `$ is equal to

$`𝐫_{\alpha \beta }=\left[(\beta ,L)(\alpha ,L)\right]𝐞_\mathrm{x}+\left[(\beta ,L)(\alpha ,L)\right]𝐞_\mathrm{y},`$ (19)
where

$$(a,b)=a\text{mod}b,(a,b)=\left[a/b\right]_{\text{int}}.$$
(20)
Here $``$ is standard modulo function and $``$ is the integer part of its argument ($`L`$ is an integer multiple of $`L_0`$). The second part of Eq.(1) includes the exchange coupling $`Jϵ_{\alpha \beta }`$ generated by dependence

$$ϵ_{\alpha \beta }=\left(\mathrm{\hspace{0.17em}1}\delta _{\varphi (\beta )1,\varphi (\alpha )L_0}\right)\left(\mathrm{\hspace{0.17em}1}\delta _{\xi (\beta )1,\xi (\alpha )L_0}\right).$$
(21)
It gives rise to geometry of dense square patterning of the period $`L_0`$ in perpendicular $`𝐞_\mathrm{x}`$, $`𝐞_\mathrm{y}`$ directions. In Eq.(21) $`\delta `$ is the Kronecker symbol and

$$\varphi (\alpha )=((\alpha ,L),L_0),\xi (\alpha )=((\alpha ,L),L_0)$$
(22)
are dependences illustrated in Fig. 16. The position of some dot including site $`\alpha `$ can be also written in terms of $``$ and $``$ functions:

$$𝐫_{i(\alpha )}^\mathrm{d}=((\alpha ,L),L_0)𝐞_\mathrm{x}+((\alpha ,L),L_0)𝐞_\mathrm{y}.$$
(23)
Here

$$i(\alpha )=L_\mathrm{d}((\alpha ,L),L_0)+((\alpha ,L),L_0)$$
(24)
is one-dimensional map associating site index $`\alpha `$ with dot index $`i(\alpha )\{0,1,2,\mathrm{},L_\mathrm{d}^21\}`$.

### A.2 Integration scheme

To derive explicit integration scheme with variable time steps $`\kappa _1\mathrm{\Delta }\tau `$ and $`\kappa _2\mathrm{\Delta }\tau `$ it is supposed that integral of spin angular frequency $`\overline{\omega }_\alpha (\tau )`$ is approximated by

$$_\tau ^{\tau +\mathrm{\Delta }\tau }\omega _\alpha (\tau ^{})\text{d}\tau ^{}=\mathrm{\Delta }\tau \overline{\omega }_\alpha (\tau ),$$
(25)
where

$$\overline{\omega }_\alpha (\tau )=a_0\omega _\alpha (\tau )+a_1\omega _\alpha (\tau \kappa _1\mathrm{\Delta }\tau )+a_2\omega _\alpha (\tau \kappa _1\mathrm{\Delta }\tau \kappa _2\mathrm{\Delta }\tau )$$
(26)
is $`\tau `$local $`\mathrm{\Delta }\tau `$ averaged $`\omega _\alpha (\tau ^{})`$ expressed through coefficients

$`a_0`$ $`=`$ $`{\displaystyle \frac{2+6\kappa _1^2+3\kappa _2+6\kappa _1(1+\kappa _2)}{6\kappa _1(\kappa _1+\kappa _2)}},`$ (27)
$`a_1`$ $`=`$ $`{\displaystyle \frac{2+3\kappa _1+3\kappa _2}{6\kappa _1\kappa _2}},`$ (28)
$`a_2`$ $`=`$ $`{\displaystyle \frac{2+3\kappa _1}{6\kappa _2(\kappa _1+\kappa _2)}}`$ (29)
derived by equating to zero the series of difference $`[\mathrm{lhs}.\mathrm{Eq}.(\text{25})\mathrm{rhs}.\mathrm{Eq}.(\text{25})]`$ considered up to the order $`𝒪(\mathrm{\Delta }\tau ^4)`$. The reason for use of non-constant prefactors $`\kappa _{1,2}1`$ is the opportunity to change adaptively integration step via the instantaneous spin rotation extreme

$$\begin{array}{cc}\mathrm{max}& \overline{\omega }_\alpha (\tau \kappa _1\mathrm{\Delta }\tau ).\hfill \\ ^{\alpha =0,\mathrm{},L^21}& \end{array}$$
(30)
In the uniform case $`\kappa _1=\kappa _2=1`$, the quadrature is compatible with Adams-Bashworth 3th order integration formula. The initial two steps estimating $`\omega _\alpha (\tau \mathrm{\Delta }\tau )`$ and $`\omega (\tau 2\mathrm{\Delta }\tau )`$ can be performed using Euler scheme obtained for $`a_0=1`$, $`a_1=a_2=0`$.

## Appendix B Figure Captions