text:
# 1 Introduction

## 1 Introduction

The cosmic ray spectrum is usually described by a superposition of power laws with index $`2.7`$ below and $`3.2`$ above the so-called “knee” at a few $`10^{15}`$eV. In the region of $`10^{18}10^{19}`$ eV the spectrum flattens, and has an index again close to $`2.7`$. There are also suggestions that the shape of the spectrum in the knee region has a fine structure which may reflect acceleration of specific groups of nuclei, e.g helium, oxygen, and iron, prominent in the Galactic abundances (see  and references therein). Recent experimental results show that the cosmic ray composition changes through the knee region of the spectrum, becoming heavier with increasing energy . For example, the CASA-MIA experiment measures the mass composition in the energy range $`10^{14}10^{16}`$ eV , and shows that above $`10^{15}`$ eV the iron group nuclei should dominate. The analysis of hybrid data from the High Resolution Fly’s Eye prototype and the MIA muon array indicates that the mass composition becomes lighter again in the region $`10^{17}10^{18}`$ eV , and is proton dominated at around $`10^{19}`$ eV .

It is usually argued that cosmic rays with energies below $`10^{14}`$ eV are accelerated at shocks of expanding supernovae (see e.g. ref. ). It seems that the acceleration mechanism of the higher energy cosmic rays may be different, although supernova shocks of various types have also been proposed  and also lead to fine structure in the spectrum and composition in the knee region . This idea is consistent with the recently observed variation in the composition of cosmic rays in the region $`10^{15}10^{16}`$eV mentioned above. However, it is uncertain what fraction of the cosmic rays above $`10^{14}`$ eV, are actually accelerated by supernova shocks, and so other mechanisms able to accelerate cosmic rays, particularly heavy nuclei, to energies above $`10^{15}`$eV may be required. It has also been suggested that protons and heavy nuclei could be accelerated to such energies in pulsar winds .

Nuclei can be also accelerated in pulsar magnetospheres as discussed by Cheng, Ho & Ruderman , Cheng & Chi  and Bednarek & Protheroe . In our previous work , we considered acceleration of iron nuclei in the magnetospheres of very young pulsars, and Crab-type pulsars, in the context of $`\gamma `$-ray and neutrino emission by supernova remnants. We found that, although the photo-disintegration process of iron nuclei during their acceleration and propagation through the neutron star magnetosphere is important, heavy nuclei would be injected into the medium of expanding supernova remnant. In the present paper we compute the possible contribution of nuclei injected by the galactic population of gamma-ray pulsars to the observed cosmic ray spectrum. We show that gamma-ray pulsars can contribute to the cosmic ray spectrum with energies above $`10^{15}`$ eV.

## 2 The pulsar spin evolution

Neutron stars are likely to be born with initial periods $`P_010`$ ms if we assume that the angular momentum of the pre-supernova star is conserved during the collapse. However, the initial angular momentum may be lost in the case of hot rapidly rotating neutron stars by emission of gravitational radiation due to the r-mode instability . As a result, the star would be left with a period of about 10 ms at about 1 year after its formation. After that, the period of the neutron star would change, mainly as a result of emission of magnetic dipole radiation with power given by

$`L_{\mathrm{em}}=B^2R_{\mathrm{NS}}^6\mathrm{\Omega }^4\mathrm{sin}^2i/6c^3,`$ (1)
where $`B`$ is the surface magnetic field at the magnetic poles, $`R_{\mathrm{NS}}`$ is the radius of the neutron star, $`i`$ is the angle between rotational and magnetic axes, and $`\mathrm{\Omega }=2\pi /P`$, where $`P`$ is the period of the star. Also, gravitational radiation may be important for the case of stars having significant ellipticity $`k`$,

$`L_{\mathrm{gr}}=32GI^2k^2\mathrm{\Omega }^6/5c^5,`$ (2)
where $`G`$ is the gravitational constant, $`I`$ is the moment of inertia .

We have investigated the population of young pulsars observed with age $`<10^6`$ years (about 100 pulsars) from the catalogue of Taylor, Manchester & Lyne  in order to estimate the neutron star ellipticity. By starting with an assumed initial period $`P_0=10`$ ms, perpendicular component of magnetic field $`B\mathrm{sin}i`$, and ellipticity $`k`$ we may use Eqs. 1 and 2 giving the rate of loss of rotational kinetic by magnetic dipole radiation and gravitational radiation, respectively, to calculate the period $`P`$ and period derivative $`\dot{P}`$ as a function of time since birth. For each of the young pulsars considered, and assuming standard values of $`I`$ and $`R_{\mathrm{NS}}`$, and the inferred value of $`B\mathrm{sin}i`$, we obtain the value of $`k`$ such that $`P`$ and $`\dot{P}`$ match those observed . We note that the values of $`k`$ obtained are rather insensitive to the assumed initial period for pulsars older than $`10^3`$ years. The results are plotted in Fig. 1, and show that a pulsar’s ellipticity must be correlated with the perpendicular component of its surface magnetic field multiplied by its radius and divided by its moment of inertia. Since the theoretical radii and masses of neutron stars lie in a rather narrow range, we conclude that there exists a correlation between the ellipticity and the surface magnetic field for neutron stars. The correlation coefficient obtained is equal to $`0.65`$, which for 100 pulsars means that the correlation is highly significant. The inferred dependence of pulsar ellipticity on surface magnetic field can be well approximated by

$`\mathrm{log}k=22.3+1.57\mathrm{log}B\mathrm{sin}i,`$ (3)
where $`B`$ is in Gauss.

The above result suggests that the ellipticity of a pulsar may be due to the deformation of the shape of the neutron star by its strong interior magnetic field (see e.g. ). The prescription for pulsar spin evolution we apply in our present work is based on Eqs. 13. Note, however, that in the case of some specific pulsars, the difference between the result from Eq. 3 and estimates based on the known age of the pulsar is significant, e.g., in the case of the Crab pulsar for which the estimated ellipticity is $`3\times 10^4`$ . In the case of the recently-discovered pulsar in the historical (AD 386) SNR G11.2-03, the standard formula  gives for the pulsar $`B=1.7\times 10^{12}`$ G, and using the historical age $`t=1612`$ yr a pulsar initial period of $`62`$ ms is found . If gravitational radiation were included, the initial period could be much shorter. We find that for the same magnetic field, an ellipticity $`k=1.49\times 10^3`$ would give an initial period of 10 ms for AD 386. This value of $`k`$, together with $`B=1.7\times 10^{12}`$ G, is consistent with the correlation described by Eq. 3 and seen in Fig. 1.

## 3 Acceleration of nuclei by young pulsars

Periodic emission of photons from radio to $`\gamma `$-ray energies gives direct evidence of the existence of electric potentials in neutron star magnetospheres which are able to accelerate particles to energies of at least $`100`$ GeV. To explain these observations, a number of models have been proposed including the polar gap model, the slot gap model, and the outer gap model (see e.g. Arons ). In the present paper, we concentrate on the outer gap model which is more attractive from the point of view of acceleration of nuclei by pulsars because it postulates induction of larger electric potentials in the pulsar’s magnetosphere. The inner gap (polar cap) models (e.g. Ruderman & Sutherland , Harding & Muslinov ) predict that the electric potentials induced in the inner gaps are significantly lower (of the order of $`10^{13}`$ V), than expected in the outer gaps for pulsars with the same parameters. Nuclei accelerated through such potentials would not suffer efficient photo-disintegration in collision with thermal photons with a characteristic temperature $`10^6`$ K emitted by the polar cap or from the whole surface of a classical radio pulsar (only very young pulsars with ages of days to a month are sufficiently hot to photo-disintegrate Fe nuclei ). Hence, we neglect this stage of pre-acceleration of Fe nuclei during their passage through the inner gap. Moreover, the inner and outer gaps do not need to extend along these same magnetic field lines, and so they do not need to be physically connected.

Nuclei, most probably iron nuclei, can be stripped off the hot surface of a neutron star due to the strong electric fields and the bombardment by particles produced in the pulsar magnetosphere . This can occur where the magnetic and rotation axes are counter-aligned, and therefore where the co-rotation charge is positive above the poles, and can pull Fe nuclei from the surface. If one assumes that pulsars have random orientations of these axes, then Fe nuclei can be extracted from a significant fraction of pulsars, but if this is not the case then our results could be significantly affected.

These nuclei are then accelerated by the pulsar outer gap electric field to high energies. We assume that the process of extraction of iron nuclei from the pulsar surface occurs during the $`\gamma `$-ray emitting stage of the pulsar as the production of $`\gamma `$-rays in the outer gap requires efficient production of relativistic $`e^\pm `$ plasma, and this also heats the polar cap region. Chen & Ruderman  distinguish two types of gamma-ray pulsars, so called Crab-like (Crab-type) and Vela-like (Vela-type). Beyond the Vela-like death line the gap cannot close, and this prevents copious production of secondary $`e^\pm `$ pairs, and causes efficient production of gamma-ray photons to cease. However, electric potentials inside the outer gap remain, and the acceleration of primary particles is still possible. In more recent papers (Romani ; Zhang & Cheng ; Cheng & Zhang ) a third type of pulsar, so called Geminga-like pulsars, is distinguished. In these pulsars the production of secondary $`e^\pm `$ pairs inside the outer gap is sustained by the soft X-ray emission from the polar cap region. The death line for Geminga-like pulsars occurs when there is insufficient flux of thermal photons from the polar cap region of the pulsar. However, primary particles can be still accelerated inside the outer gap even after a pulsar has crossed the Geminga-like death line (only pair production stops). We expect that the injection of the Fe nuclei from the region of the polar cap is related to the thermal radiation field produced by this region of the neutron star surface. The thermal radiation field drops significantly at a neutron star age $`10^6`$ yr (see Fig. 2 in Romani ). We estimate that at this age the pulsar stops injecting Fe nuclei into the magnetosphere.

### 3.1 Model for the electric field in the outer gap

During the $`\gamma `$-ray production stage of pulsars, three different types of outer gap can be distinguished, depending on the relative importance of different processes of $`\gamma `$-ray production: the Crab-type and Vela-type gaps (Cheng, Ho, Ruderman ), and the Geminga-type gap (Romani , Zhang & Cheng ). The separation between the Crab-type and the Vela-type outer gaps depends on the pulsar period and surface magnetic field, and is given by Eq. 25 of Chen & Ruderman

$`P_{\mathrm{CV}}0.04B_{12}^{2/5},`$ (4)
where $`B_{12}`$ is the pulsar surface magnetic field in $`10^{12}`$ G, and the period $`P_{\mathrm{CV}}`$ is in seconds. The maximum potential difference across different types of outer gap can be approximated by the simple formula (Cheng , private communication),

$`\mathrm{\Phi }(B,P)6.6\times 10^{12}f^2B_{12}P^2\mathrm{V}.`$ (5)
This formula gives the potential drop in outer gaps a factor of two lower than derived in original paper by Cheng, Ho & Ruderman , and is more consistent with the recent works on the outer gap model . However, we note that this simple formula may not match all observational data well because of uncertainties of viewing angle and inclination angle which are very important in determining the size of the outer gap. The parameter $`f`$ is used to parametrize the gap voltage and power, and is defined by Cheng, Ho & Ruderman

$`f82B_{12}^{13/20}P^{33/20},`$ (6)
Formulae (5) and (6) give the potential difference $`3\times 10^{14}`$ eV for the Crab pulsar. For Geminga-type gaps (Zhang & Cheng )

$`f5.5B_{12}^{4/7}P^{26/21}.`$ (7)
For Vela-type outer gaps a simple definition of the parameter $`f`$ is not available. We assume that in this case the potential drop in the outer gap is equal to a constant fraction of the total potential drop in the outer gap given by Eq. (5) with $`f^2=0.15`$ for all periods (see ). The dependence of the potential difference across the outer gap as a function of the pulsar period is plotted in Fig. 2 for three different values of the surface magnetic field of the pulsar.

According to Cheng, Ho & Ruderman  the potential drop across the gap depends also on the position of injection into the outer gap field, $`z`$ ($`z`$ ranges between 0 and $`a`$), in relation to the maximum gap width $`a`$. Therefore, we assume that the potential drop for position $`z`$ is given by (see definitions below eq. 3.10 in ref. )

$`\mathrm{\Phi }(z)=4\mathrm{\Phi }(B,P){\displaystyle \frac{z}{a}}[1{\displaystyle \frac{z}{a}}].`$ (8)
### 3.2 Radiation field inside the outer gap

Iron nuclei accelerated in the outer gap potential can suffer significant photo-disintegration in collisions with photons produced inside the gap. We estimate the density of soft photons inside outer gaps of pulsars with a given period and magnetic field by scaling from the density of soft photons inside the outer gap of the Crab pulsar assuming the outer gap model. This density can be estimated from the observed pulsed flux of soft photons, assuming that the dimension of the outer gap is of the order of the light cylinder radius. Ho  has shown, by analyzing a sample of Crab-type pulsars with periods $`15`$ ms, that the luminosity of the photons produced in the gap does not change drastically for Crab-type outer gaps. Based on the results of Ho , we interpolate/extrapolate the density of photons in the present Crab outer gap for other periods. The $`\gamma `$-ray luminosity of pulsars scales as $`f^3`$ times the total spin down pulsar luminosity , and so the $`\gamma `$-ray luminosity of the Crab type pulsars depends on the magnetic field as $`B^{1/20}`$ (see Eq. 6), i.e. it is only very weakly dependent on $`B`$. Since energies of synchrotron photons are proportional to the magnetic field in the gap, then the photon density in Crab-type outer gaps should be inversely proportional to the pulsars’ magnetic fields.

We note that the observations seem to suggest that the $`\gamma `$-ray luminosities of different types of pulsar are proportional to $`BP^2`$ (e.g. Fig. 7 in Arons ), although different models predict different dependencies of $`\gamma `$-ray luminosities on $`B`$ and $`P`$ . However the above dependence should be considered with some caution because it interpolates between the Geminga pulsar (and PSR 1055-52) and the Vela pulsar (and its twin PSR 1706-44). The pulsar PSR 1951+32 deviates by an order of magnitude and the Crab pulsar by a factor of two, and if the distance to the Vela pulsar is $`250`$ pc but not $`500`$ pc , then Vela pulsar deviates by a factor of four from this dependence. Also, the distance to Geminga pulsar is not well known.

The density of soft photons in Vela-type outer gaps is obtained based on the observed change of the pulsed luminosity of the Vela type pulsars (e.g. ), and our estimate of $`f`$ in the case of Vela-type pulsars, with the normalization to photon density in the Crab-type outer gap at the transition period given by Eq. (4). A good fit to the luminosity of the Vela-type pulsars is obtained if the density of photons in the outer gap drops with period as $`P^3`$. Therefore, we apply the following model for the soft photon density in Crab and Vela-type pulsars for different periods and surface magnetic fields

$`n_{\mathrm{ph}}n_{\mathrm{ph}}^{\mathrm{Crab}}\{\begin{array}{cc}\frac{P_{\mathrm{Crab}}B_{\mathrm{Crab}}}{PB},\hfill & \text{the Crab-type;}\hfill \\ \frac{P_{\mathrm{Crab}}B_{\mathrm{Crab}}}{PB}\left(\frac{P_{\mathrm{CV}}}{P}\right)^3,\hfill & \text{the Vela-type,}\hfill \end{array}`$ (9)
where $`P_{\mathrm{Crab}}=0.033`$ seconds, $`B_{\mathrm{Crab}}=4\times 10^{12}`$ G and $`n_{\mathrm{ph}}^{Crab}`$ is the density of photons in the outer gap of the Crab pulsar at present (Cheng, Ho & Ruderman 1986).

The Geminga-type outer gap is determined by the thermal photons produced on the neutron star surface. This thermal radiation field for pulsars with ages corresponding to Geminga-type pulsars is too low to create a significant target for photo-disintegration of Fe nuclei. The epoch of gamma-ray emission of Geminga-type pulsars is determined by the temperature of the neutron star surface which, according to models of neutron star evolution, drops rapidly after $`10^6`$ years (see e.g. Fig. 2 in Romani ). Therefore, Geminga-type outer gaps should not be able to operate after $`10^6`$ years from the explosion of supernova. Here we assume that the process of extraction of Fe nuclei from the neutron star surface does not occur after this age because, as we noted above, it is probably caused by heating of the polar cap region by $`e^\pm `$ pairs created in the pulsar magnetosphere.

### 3.3 Number of nuclei accelerated in the outer gap

As in our previous paper , we assume that the number of Fe nuclei injected per unit time, $`\dot{N}_{\mathrm{Fe}}`$, is simply related to the total power output of the pulsar $`L_{\mathrm{em}}(B,P)`$ ,

$`\dot{N}_{\mathrm{Fe}}=\xi L_{\mathrm{em}}(B,P)/Ze\mathrm{\Phi }(B,P),`$ (10)
where $`\xi `$ is a free parameter describing the fraction of the total power taken by relativistic nuclei accelerated in the outer gap, $`Z=26`$ is the atomic number of Fe, and $`\mathrm{\Phi }(B,P)`$ is defined by Eq. 5. The parameter $`\xi `$ has to be always less than one since the power of accelerated Fe nuclei can never be greater than the spin-down power of the pulsar. The charge density estimated from Eq. 10 exceeds the ‘classical’ Goldreich-Julian charge density close to the polar caps $`\rho =(\mathrm{\Omega }B/2\pi c)[1(\mathrm{\Omega }r/c)^2\mathrm{sin}^2\theta ]^1\mathrm{\Omega }B/2\pi c,`$ where the approximate formula is valid close to the surface where $`\mathrm{\Omega }r\mathrm{sin}\theta c`$. However, as noted by Goldreich & Julian  below their Eq. (8), the above formula is only valid in the co-rotation portion of the magnetosphere. If there is additional negative current onto the surface accompanying the charge injected from the surface in the form of Fe nuclei, then the above formula is not valid. Positive charge may flow onto the surface along the same magnetic lines as Fe nuclei, and could be produced somewhere in the pulsar magnetosphere, not necessary close to the pulsar’s surface (e.g. possibly in the outer gap). Pair production just above the polar cap does not need to occur in this case.

### 3.4 Number of neutron stars with different parameters in the Galaxy

As discussed above, we assume that all pulsars are born with periods of about 10 ms, or that they reach this period soon after formation by emission of gravitational waves due to the r-mode instability. As we argued in Section 2, pulsars need to have a significant ellipticity $`k`$ in order to reach their observed periods from a period at birth of 10 ms.

The distribution of the pulsar magnetic fields can be evaluated from observations, assuming that the surface magnetic field does not decay significantly with pulsar age. For example, Narayan  approximates the distribution of surface magnetic fields of pulsars at birth by

$`dN/d(\mathrm{log}B)0.065(10^{12}G/B)\mathrm{yr}^1,`$ (11)
for $`B>2\times 10^{12}`$ G. We estimate the distribution of pulsars with surface magnetic fields below $`2\times 10^{12}`$ G based on the Fig. (13) in Narayan . It is likely that this approximation underestimates the number of pulsars with relatively low magnetic fields because of observational selection effects. The above formula is normalized to a formation rate of pulsars in the Galaxy equal to 1 per 70 years.

## 4 Acceleration and propagation of nuclei from the neutron star surface

Nuclei extracted from the neutron star surface are accelerated and photo-disintegrated in the outer gap. Their disintegration products, and those which escape without any disintegration, are injected into the nebula and lose energy in collisions with matter ejected during the supernova explosion, and as a result of the expansion of nebula. All these processes are discussed below.

### 4.1 Propagation of nuclei through pulsar magnetospheres

The importance of photo-disintegration of nuclei with different mass numbers during propagation of nuclei in the radiation field of the outer gap is determined by the reciprocal mean free path, which we calculate using cross sections of Karakuła & Tkaczyk  (with improvements given in ref. ). The results of such a calculation for an outer gap with parameters applicable to the Crab pulsar are shown in Fig. 1 of Bednarek & Protheroe . It is evident that multiple photo-disintegrations of primary Fe nuclei will occur for this pulsar. In our present calculations, the photo-disintegration of nuclei in the thermal radiation field emitted by the pulsar surface is neglected since it may be only important during a few days after creation of the pulsar when its surface temperature is of the order of $`10^7`$ K , and consider only photodisintegration in the outer gap’s non-thermal radiation field.

We simulate the photo-disintegration process of nuclei accelerated by pulsars using the evolution model, initial parameters of pulsars, the outer gap model for acceleration, and the radiation field discussed in Sections 2 and 3. Multiple photo-disintegrations occur only for pulsars with Crab and Vela-type outer gaps. We simulate the acceleration and propagation of iron nuclei in pulsar outer gaps in order to obtain the energy spectrum of neutrons and protons (extracted from Fe nuclei) and the spectra of escaping secondary nuclei. Note, however, that the fate of neutral particles (neutrons) and charged particles (protons and secondary nuclei) is different. The neutrons move ballistically through the surroundings of the pulsar, and decay at distances determined by their Lorentz factors, and this can be inside or outside of the expanding envelope of the supernova. We compute the spectra of neutrons (and so protons from their decay) which decay inside and outside the Nebula as in Bednarek & Protheroe . In the case of protons from neutrons decaying outside the nebula, we assume that all of them escape into the Galactic medium without being overtaken by the expanding nebula (we drop the second term in brackets in Eq. (5) in Bednarek & Protheroe  which represents the effects of protons from neutrons decaying outside the nebula being overtaken and trapped by the expanding nebula at some later time). This term can be neglected since the minimum diffusion distance of these protons, with Lorentz factor $`\gamma _\mathrm{p}=10^5`$, in the Galactic magnetic field $`5\times 10^6`$ G is comparable to the radius of nebula during its time of free expansion .

Protons and secondary nuclei injected into the Nebula are trapped by the magnetic field of the nebula, and suffer adiabatic and collisional energy losses during the time of expansion of the nebula. These effects are considered in the next subsection.

### 4.2 Interaction of nuclei with matter inside nebulae

Soon after the supernova explosion, when the nebula was relatively small, nearly all the energetic neutrons would be expected to decay outside the nebula. However, at early times we must take account of collisions with dense matter of the expanding nebula. The optical depth for neutrons moving ballistically through the envelope may be estimated from $`\tau _{nH}\sigma _{np}n_Hr8.6\times 10^{14}M_1v_8^2t^2`$, where $`M=M_1M_{}`$ is the mass ejected during the supernova explosion in units of solar masses, $`r=vt`$, $`n_H=M/(4/3\pi r^3m_p)`$ is the number density of target nuclei, and $`v=10^8v_8`$ cm s<sup>-1</sup> is the expansion velocity of the nebula. We note that $`\tau _{nH}=1`$ at $`t0.93M_1^{1/2}v_8^1`$ yr.

Charged nuclei (protons, protons from decay of neutrons, and secondary nuclei) injected by the pulsar are captured by the nebula’s magnetic field and interact with the dense material of the expanding supernova remnant. The optical depth for these nuclei depends on the time, $`t`$, of the injection of nuclei, and can be estimated from

$`\tau _{\mathrm{Ap}}(t)1.3\times 10^{17}A^{2/3}M_1v_8^3t^2.`$ (12)
Interactions of nuclei with matter inside the nebula become negligible when the optical depth becomes less than one. For the parameters derived for the Crab, $`v_8=2`$ and $`M_1=3`$, and the average value of the mass number of nuclei injected from the pulsar magnetospheres into the nebula $`<A>48`$, the nebula becomes transparent for nuclei at times $`t>3.6\times 10^8A^{1/3}M_1^{1/2}v_8^{3/2}22`$ years. Nuclei injected at earlier times should suffer significant fragmentation in collisions with the matter inside the nebula. We have simulated this process of fragmentation of nuclei injected by the pulsar applying the semi-empirical formulae for the cross sections for the fragmentation of nuclei of Silberberg & Tsao  with parameters given in their tables 1A,B,C,D. In our simulation, we assumed that in a single collision a nucleus fragments into two nuclei with different mass numbers.

Protons created during complete fragmentation of nuclei are considered as being captured by the magnetic field of the nebula. In contrast, neutrons which move ballistically through the nebula have some chance of decaying inside or outside the nebula. These effects have been included in our simulation code. We consider the interactions of these secondary protons and neutrons with matter inside the nebula in the same way as for protons and neutrons from photo-disintegration of iron nuclei in the pulsar’s magnetosphere.

### 4.3 Escape of nuclei from nebulae

We assume that nebulae expand during the free expansion phase with the velocity observed in the case of the Crab Nebula, i.e. 2000 km s<sup>-1</sup>. When the amount of matter ejected during the supernova explosion becomes comparable to the amount of interstellar matter swept up by the nebula, then the nebula reaches the Sedov phase. The Sedov phase starts when the nebula reaches radius

$`R_\mathrm{S}=6.8\times 10^{18}(M_1/n)^{1/3}\mathrm{cm},`$ (13)
and this occurs at time

$`t_\mathrm{S}=R_\mathrm{S}/v,`$ (14)
after the explosion. Here $`n`$ is the density of the surrounding medium (cm<sup>-3</sup>), and $`v`$ is the expansion velocity of supernova during its phase of free expansion. For a supernova with an ejecta mass of $`M_{\mathrm{SN}}=3M_{}`$, a velocity of free expansion $`v=2000`$ km s<sup>-1</sup> as observed for the Crab Nebula, and a density of the interstellar medium $`n=0.3`$ cm<sup>-3</sup>, the Sedov phase should begin when the nebula reaches a radius of $`R_\mathrm{S}4.8`$ pc at about $`t_\mathrm{S}2000`$ years. This is only about twice as large as the present radius (and age) of the Crab Nebula. We assume the following approximation for the radius of a nebula at time $`t`$

$`R_{\mathrm{Neb}}(t)=\{\begin{array}{cc}vt,\hfill & \text{if }tt_\mathrm{S}\text{;}\hfill \\ vt_\mathrm{S}^{3/5}t^{2/5},\hfill & \text{if }t>t_\mathrm{S}\text{.}\hfill \end{array}`$ (15)
The magnetic field inside the nebula changes with time and depends on the distance from the center of the nebula, i.e. the pulsar. We estimate it applying the results of the Kennel & Coroniti  model, which gives the distribution of the magnetic field inside the nebula as a function of the parameter $`\sigma `$ (the ratio of the magnetic energy flux to the particle energy flux at the location of the pulsar wind shock at radius $`R_{\mathrm{sh}}`$). We estimate the magnetic field inside the nebula (in units of the magnetic field strength $`B_{\mathrm{sh}}`$ at the shock location) as a function of its radius based on the figs. 3 and 7 in Kennel & Coroniti . For nebulae with $`\sigma 0.003`$, the magnetic field drops significantly with the distance from the center of the nebula (fig. 3 in ), and can reach a value typical for the interstellar medium ($`5\times 10^6`$ G) at distance $`R_\mathrm{B}`$ which is less than the radius of the expanding nebula $`R_{\mathrm{Neb}}`$, and so we assume that beyond radius $`R_\mathrm{B}`$, the magnetic field is equal to $`5\times 10^6`$ G.

The magnetic field at the shock, $`B_{\mathrm{sh}}`$, can be found from the dependence of the pulsar’s magnetic field in the region of the inner magnetosphere ($`BR^3`$), and the pulsar wind zone ($`BR^1`$). It is given by

$`B_{\mathrm{sh}}\sqrt{\sigma }B_\mathrm{s}\left({\displaystyle \frac{R_\mathrm{s}}{R_\mathrm{L}}}\right)^3\left({\displaystyle \frac{R_\mathrm{L}}{R_{\mathrm{sh}}}}\right),`$ (16)
where $`R_\mathrm{s}`$ and $`B_\mathrm{s}`$ are the radius and the surface magnetic field of the pulsar, and $`R_\mathrm{L}=cP/2\pi `$ is the light cylinder radius of the pulsar. The value of $`\sigma `$ is estimated for the Crab pulsar to be $`0.003`$ (e.g. Kennel & Coroniti ). However, it has recently been found that for the Vela pulsar the value of $`\sigma `$ is rather close to 1 (Helfand et al. ). To our knowledge, there is no available information on $`\sigma `$ for other plerions. Therefore, in order to fulfill observational constraints of these two pulsars, we assume that pulsars with periods lower than the period of the Crab pulsar have $`\sigma 0.003`$ and pulsars older than the Vela pulsar have $`\sigma 1`$. For pulsars having a rate of the rotational energy loss between that of the Crab and Vela pulsars, $`\sigma `$ is obtained by linear interpolation between these two values.

The location of the wind shock, $`R_{\mathrm{sh}}`$, can be estimated from the comparison of the wind energy flux at the shock with the pressure of the outer nebula, which is determined by the supply of the magnetic energy to the nebula by the pulsar over its all lifetime (Rees & Gunn ). The location of the shock at an arbitrary time $`t_\mathrm{p}`$ after explosion of the supernova can be estimated from

$`{\displaystyle \frac{L_{\mathrm{em}}(t_\mathrm{p})}{4\pi R_{\mathrm{sh}}^2c}}{\displaystyle \frac{_0^{t_\mathrm{p}}\sigma L_{\mathrm{em}}(t)𝑑t}{\frac{4}{3}\pi R_{\mathrm{Neb}}^3}}.`$ (17)
As an example, we show in Fig. 3 the characteristic radii in the nebula, $`R_{\mathrm{Neb}}`$, $`R_{\mathrm{sh}}`$, and $`R_\mathrm{B}`$, for the following parameters: $`B_\mathrm{s}=3\times 10^{12}`$ G, $`k=3\times 10^4`$, $`v=2000`$ km s<sup>-1</sup>, $`M_{\mathrm{SN}}=3M_{}`$, and $`n=0.3`$ cm<sup>-3</sup>. The location of the shock, estimated from the above formula, is consistent with the observations of the shock in the Crab  and Vela  nebulae. For older nebulae $`R_\mathrm{B}`$ becomes less than $`R_{\mathrm{Neb}}`$, i.e. the magnetic field inside the outer regions of expanding nebula drops to $`5\times 10^6`$ G.

Charged particles, nuclei and protons, injected into the nebula will diffuse, and finally escape from it into interstellar space. We assume that particles escape from a nebula at time $`t_{\mathrm{max}}`$ when the diffusion distance of particles injected at the time $`t`$ becomes comparable to the radius of the nebula $`R_{\mathrm{Neb}}`$ at the time $`t_{\mathrm{max}}`$. In order to estimate this time at which they escape from the nebula, we compare the diffusion distance for nuclei with the mass number $`A`$ and the Lorentz factor $`\gamma _\mathrm{A}`$ in the time dependent magnetic field inside an expanding nebula with radius given by Eq. 15. The distance diffused depends on the diffusion coefficient which we take to be

$`D=D_0\rho ^{1/3},`$ (18)
as appropriate to a Kolmogorov spectrum of turbulence. Here, $`\rho E/eZ`$ is the particle rigidity, and $`D_0`$ is normalized in such a way that $`D=r_\mathrm{L}c/3`$ when the Larmor radius, $`r_\mathrm{L}`$, is equal to the radius of the remnant (largest scale) at time t (given by Eq. 15). This is a conservative assumption as it will give rise to the greatest trapping of nuclei in the nebula (and highest adiabatic losses) while being consistent with a Kolmogorov spectrum of turbulence. However, the precise form adopted for $`D`$ has little effect on our results as most particles are released at $`t_{\mathrm{max}}`$. The diffusion distance is obtained by integration over time from the injection of nuclei at $`t_{\mathrm{inj}}`$ up to escape at $`t_{\mathrm{esc}}`$

$`r_{\mathrm{dif}}={\displaystyle _{t_{\mathrm{inj}}}^{t_{\mathrm{esc}}}}\sqrt{{\displaystyle \frac{3D}{2t^{}}}}𝑑t^{}.`$ (19)
Note that the diffusion coefficient depends on the age of the nebula, and this is included when integrating Eq. 19. For nuclei with sufficiently low energies, $`t_{\mathrm{max}}`$ may be shorter than $`t_{\mathrm{esc}}`$, and we then take $`t_{\mathrm{esc}}=t_{\mathrm{max}}`$. The adiabatic losses of particles inside nebula are considered up to the moment when the expansion velocity of the nebula during the Sedov phase, $`v_{\mathrm{exp}}=v(t_\mathrm{S}/t)^{3/5}`$, drops to a few tens km $`s^1`$ (i.e. to the characteristic velocities of the clouds in the interstellar medium). This happens at the age of $`10^6`$ yr. Note that at this age the magnetic field inside the nebula is larger than $`5\times 10^6`$ G only at small region above the radius of the shock ($`R_\mathrm{B}2R_{\mathrm{sh}}`$, see Fig. 3).

Adiabatic and interaction energy losses of nuclei during the expansion of the nebula are included by estimating the Lorentz factor of nuclei at arbitrary time, $`t`$,

$`\gamma (t)=\gamma (t_{\mathrm{inj}}){\displaystyle \frac{t_{\mathrm{inj}}+t}{2t_{\mathrm{inj}}K^{\tau _{\mathrm{Ap}}(t)}}},`$ (20)
where $`K`$ is the inelasticity coefficient, assumed equal to 0.5 for protons and neutrons. For nuclei with mass number $`A`$, we take the inelasticity coefficient to be $`K0.5/A^{1/3}`$ since the probability of a collision of a single nucleon within the energetic nucleus, with a proton of the remnant nebula is, on average, proportional to $`(\sigma _{Ap}/A)(A^{2/3}/A)A^{1/3}`$.

## 5 Injection rate of nuclei into the Galaxy

Taking into account all of the effects discussed above, we compute the injection rate of nuclei of different mass numbers by the galactic population of pulsars into the interstellar medium, and this is given in Fig. 4(a). The spectrum shows two distinct maxima corresponding to heavy nuclei (dot-dot-dot-dash curve peaking at a few $`10^{15}`$eV) and to protons (solid curve peaking at $`10^{14}`$eV). The contribution of medium mass nuclei is visible at energies between the proton maximum and the heavy nuclei maximum. The main contribution to the spectrum of nuclei below the high-energy peak is from Crab-type pulsars. This is clearly seen in Fig. 4(b) in which we present the contribution of pulsars having different types of outer gap to the total spectrum of nuclei injected into the Galaxy the from pulsars. The main contribution to the maximum at $`3\times 10^6`$ GeV comes from Vela-type and Geminga-type pulsars. However the nuclei with the highest energies are injected by the Vela and Crab-type pulsars as expected from inspection of the electric potentials shown in Fig. 2.

## 6 Cosmic ray propagation

We assume that the radial distribution of the cosmic ray emission from gamma-ray pulsars follows the distribution of old neutron stars in the disk of the galaxy as given by Paczynski

$$q(r)2\pi rdr=a_R\mathrm{exp}(r/R_{\mathrm{exp}})R_{\mathrm{exp}}^22\pi rdr$$
(21)
with $`a_R=1.0683/2\pi `$, $`R_{\mathrm{exp}}=4.5`$ kpc, and $`r`$ being the the distance from the centre of the Galaxy of radius $`r_{\mathrm{max}}=20`$ kpc, and that they are uniformly distributed in a disk of thickness $`2h_g`$. For $`h_g`$ we take an estimate of the scale height of pulsar birth locations, $`h_g=130`$ pc , since random pulsar velocities $`200`$ km s<sup>-1</sup> would not significantly increase the scale height for pulsars much younger than $`10^6`$ y. Combining $`q(r)`$ with the calculated spectra of cosmic rays of atomic number $`Z`$, $`Q_Z(E)=dN/dE/dt`$ presented in Fig. 4(a), we obtain the spatial distribution of cosmic ray injection

$$\zeta (r,E)=q(r)Q_Z(E)/(2h_g).$$
(22)
We assume that at energies above $`10^3`$ GeV we can reasonably neglect interactions and energy losses of cosmic ray nuclei, and that diffusion is the dominant process in shaping the observed cosmic ray spectrum. We consider two different models (Ptuskin et al ) for propagation of cosmic rays in the galaxy which are consistent with the observed secondary to primary ratios: (i) the minimal reacceleration model, and (ii) the standard diffusion model. We use the best-fitting parameters found by Ptuskin et al  for a halo half-thickness $`h_h=5`$ kpc (we note here that the cosmic ray intensity at Earth obtained is rather insensitive to the value of $`h_g`$ chosen provided that $`h_gh_h`$). For the minimal reacceleration model, the diffusion coefficient used is

$$D=3.8\times 10^{29}(\rho /1\mathrm{GV})^{1/3}\mathrm{cm}^2\mathrm{s}^1$$
(23)
where $`\rho E/Z`$ is the particle rigidity. For the standard diffusion model, the diffusion coefficient used is

$$D=3.0\times 10^{28}(\rho /1\mathrm{GV})^{0.54}\mathrm{cm}^2\mathrm{s}^1.$$
(24)
Assuming cylindrical geometry, we then use Eq. 3.22 of Berezinsky et al. , to calculate the cosmic ray density at Earth assumed to be in the Galactic plane 8 kpc from the centre of the Galaxy.

The spectra after propagation are plotted in Fig. 5 for the various species where they are compared with the world data set of observed cosmic ray intensity . Fig. 5(a) and (b) show results for the two propagation models. The results are rather insensitive to the propagation model used, and this can be easily understood after inspecting the injection spectra in Fig. 4(a) and Eqs. 2324. The peak at $`10^5`$ GeV is due to protons with rigidity $`10^5`$ GV, and the peak at $`3\times 10^6`$ GeV is due to iron nuclei with rigidity $`1.5\times 10^5`$ GV, approximately the same as the proton peak. Looking at Eqs. 23 and 24, we note that the diffusion coefficients for the two propagation models are equal at $`1.8\times 10^5`$ GV, and this is within less than a factor of 2 of the rigidities of both the proton peak and the iron peak, and so the two propagation models will give almost identical fluxes for the same pulsar injection spectrum. The propagated fluxes are interestingly close to those observed, with the region between $`3\times 10^6`$ GeV and $`2\times 10^7`$ GeV being about $`20`$% of the observed flux.

## 7 Discussion and Conclusion

As can be seen in Fig. 5, we can expect fine structure in the spectrum in the region of the knee, as well as subtle composition changes. Our calculations were made assuming 100% efficiency of conversion of spin-down power into accelerated iron nuclei, i.e. $`\xi =1`$ in Eq. 10. Other uncertainties entering our calculation relate the propagation of cosmic rays. The diffusion coefficients assumed were based on an analysis of lower energy cosmic ray data, and their extrapolation to higher rigidities may be questionable. This would increase our predicted flux if the actual diffusion coefficient were lower. Also, the pulsar formation rate is rather poorly known. With all these uncertainties, we feel justified is adopting $`\xi 1`$ in order see whether injection of cosmic rays accelerated in pulsar outer gaps could be potentially interesting.

We shall next examine the composition in the knee region of the spectrum. The mass composition below $`10^{15}`$ eV is characterized by an average value of $`<\mathrm{ln}A>1.2`$ . If we assume that this component of the spectrum extends to higher energies, and add to it the contribution of heavy nuclei (with average $`A45`$) as obtained for our pulsar acceleration model, this would change the average value of $`<\mathrm{ln}A>`$ to $`2.4`$ above the knee. Such a change in $`<\mathrm{ln}A>`$ has recently been reported for the CASA-MIA air shower array data .

Our calculations have assumed a steady-state production of pulsar-accelerated cosmic rays throughout the Galaxy with a smooth spatial distribution of sources following that inferred for pulsars. It neglects the effect of relatively nearby sources, such as any within 1 kpc releasing cosmic rays during the past $`10^4`$ years which would enhance the pulsar contribution, and so the effects of pulsar acceleration on the cosmic ray composition could be even more dramatic. Finally, we conclude that pulsar-accelerated cosmic rays could make a significant contribution of heavy nuclei to the cosmic rays observed in the knee region, and that this would manifest itself as fine structure in the spectrum, and a heavy composition at $`3\times 10^6`$$`2\times 10^7`$ GeV.

## Acknowledgements

We would like to thank the referee for useful comments and suggestions, K.S. Cheng for information on the details of the outer gap model, and T.A. Porter for discussions on cosmic ray propagation. The work of WB is supported by the Komitet Badań Naukowych through the grants 2P03D 001 14 and 2P03C 006 128. The research of RJP is supported by the Australian Research Council.


format: text

text:
# CLUSTERING IN DEEP (SUBMILLIMETRE) SURVEYS

## 1 A model for the angular clustering

Fig.1 shows fluctuations in the galaxy counts as measured in the APM Galaxy Survey . Symbols with errorbars show the square root of the variance $`\overline{w}_2(\theta )`$, e.g. $`\mathrm{\Delta }N/N\sqrt{\overline{w}_2}`$, measured in squared cells of area $`A\theta ^2`$. The data can be described by a power law:

$$\mathrm{\Delta }N/N\sqrt{\overline{w}_2}(A/A_0)^\beta $$

with $`A_07.6\times 10^5`$ sq.deg. and $`\beta 0.175`$, for $`A<2.7`$ sq. deg. For larger areas there is an exponential cut-off with a characteristic scale of $`A_c110`$ sq. deg. The complete model is represented by:

$$\frac{\mathrm{\Delta }N}{N}(A/A_0)^\beta \mathrm{exp}[A/A_c],$$
(1)
and is shown as a solid curve in Fig.1.

## 2 Shot-noise

The effects of shot-noise on the counts, e.g. due to the small number of sources in the sample, is easy to quantify. If the mean number density of sources at a given flux is $`𝒩`$, then the number of sources in a sample of area $`A`$ is: $`N=A𝒩`$. The total variance on such area is:

$$\overline{w}_2^{total}=\overline{w}_2+\frac{1}{N},$$
(2)
where $`\overline{w}_2`$ is the intrinsic variance (in the case of high density) . Thus, the number counts variations due to intrinsic clustering and shot-noise is:

$$\frac{\mathrm{\Delta }N}{N}\left[\overline{w}_2[\theta ;z;𝒟]+\frac{1}{A𝒩}\right]^{1/2}$$
(3)
where $`\overline{w}_2[\theta ;z;𝒟]`$ is the angular variance of depth $`𝒟`$ at redshift $`z`$. We next need to quantify the effects of projection and clustering evolution.

## 3 Projection effects

Let us first assume that the 2-pt function $`\xi _2(r_{12})`$ does not evolve in co-moving coordinates: $`\xi _2(r_{12})=\xi _2(x_{12})`$, where $`r_{12}=x_{12}/(1+z)`$. We then have that:

$$\overline{w}_2(\theta )=_V𝑑\text{x}_1𝑑\text{x}_2\psi (x_1)\psi (x_2)\xi _2(\text{x}_{12}),$$
(4)
where $`d\text{x}_i`$ is the co-moving volume, $`\psi (x)`$ is the normalized probability that a galaxy at a coordinate $`x`$ is included in the catalogue and $`V`$ is a cone (or pyramid) of radius (or side) $`\theta `$ and infinite depth (with solid angle equal to the sample area $`A`$). In the small angle approximation, e.g. when the transverse distances are much smaller than the radial ones, we can relate the angular clustering of samples projected at two different characteristic co-moving depths $`𝒟_2`$ and $`𝒟_1`$ by:

$$\overline{w}_2[\theta ;𝒟_1]=\frac{𝒟_2}{𝒟_1}\overline{w}_2[\theta 𝒟_2/𝒟_1;𝒟_2].$$
(5)
We take $`𝒟_2400h^1\mathrm{Mpc}`$ as the mean APM depth and $`𝒟_12500h^1\mathrm{Mpc}`$ as the mean co-moving depth of a typical submm sample (the actual number depends on both the mean redshift and the cosmological parameters, as given by the luminosity distance relation). Thus the curve that fits the APM measurements (solid curve) in Fig.1 must be moved left by a factor $`𝒟_2/𝒟_1`$ to account for the fact that, at larger radial distances, the same physical length subtends a smaller angle. The APM curve must also be moved down by a factor $`𝒟_2/𝒟_1`$ because more galaxies are seen in projection at greater radial distances. The resulting prediction in this case, scaled from Eq. to $`𝒟_12500h^1\mathrm{Mpc}`$ with Eq., is shown as the dotted line in Fig.1. The dotted-dashed line includes the shot-noise contribution in Eq. for a submm 850$`\mu `$m flux limit of 3 mJy, with surface density $`𝒩=2000`$/sq. deg., which are characteristic of the HDF counts. Note how in this case the variance of the counts are dominated by shot-noise at almost all scales. These predictions agree well with the source-count variations found in our submm simulations of the HDF (see Hughes & Gaztañaga in these proceedings), which also have clustering fixed in co-moving coordinates and the same mean depth.

## 4 The redshift evolution of clustering

We next model the redshift evolution of $`\xi _2`$ in proper coordinates $`r_{12}=x_{12}/(1+z)`$ as:

$$\xi (r_{12};z)(1+z)^{(3+ϵ)}\xi (r_{12};0)$$
(6)
For stable clustering (pattern fixed at proper separations) we have $`ϵ0`$ at small scales . For a power-law correlation with slope $`\gamma 1.7`$, linear theory gives $`ϵ=\gamma 10.7`$ and non-linear growth gives $`ϵ1`$. If the clustering is fixed in co-moving coordinates, then $`ϵ\gamma 31.3`$ which produces even less evolution than in the stable clustering regime. This fixed clustering model describes galaxies which are identified with high density-peaks: peaks move less than particles, which results in less evolution. Weak evolution is consistent with the strong clustering observed in Lyman-break galaxies at $`z3`$, which is comparable to the clustering of present-day galaxies . The above model of cluster evolution projects as:

$$\overline{w}_2[\theta ;z_1;𝒟]\left(\frac{1+z_2}{1+z_1}\right)^{3+ϵ\gamma }\overline{w}_2[\theta ;z_2;𝒟]$$
(7)
where $`\gamma `$ varies from $`1.6`$ over the smaller scales to $`\gamma 2`$ near the exponential break. Here we use $`z_20.15`$ for the APM depth and $`z_12`$ for the depth in sub-mm samples. The resulting predictions, scaled from Eq. to $`𝒟_12500h^1\mathrm{Mpc}`$, are shown as the dotted ($`ϵ=\gamma 3`$), short-dashed line ($`ϵ=0`$) and long-dashed lines ($`ϵ=\gamma 1`$) in Fig.1.

## 5 Discussion

In summary, our recipe for the angular variance of number count fluctuations $`(\frac{\mathrm{\Delta }N}{N})^2\overline{w}_2`$ in a galaxy sample of area $`A=\theta ^2`$, mean co-moving depth $`𝒟`$ and mean redshift $`z`$, is given by:

$$\overline{w}_2[\theta ;z;𝒟]\frac{400h^1\mathrm{Mpc}}{𝒟}\left(\frac{1.15}{1+z}\right)^{3+ϵ\gamma }(\theta /\theta _0)^{2\beta }e^{(\theta /\theta _c)^2}+\frac{1}{\theta ^2𝒩}$$
(8)
where $`\beta 0.175`$, $`\theta _c10.5`$ deg., $`\theta _08.7\times 10^3`$ deg., $`3+ϵ\gamma `$ is between 0 and 2, depending on clustering evolution, and $`𝒩`$ is the mean number density of sources at the given flux that produces shot-noise fluctuations (e.g. Eq.). Several cases for the above model are shown in Fig.1. Assuming the shot-noise is negligible, e.g. $`𝒩\mathrm{}`$, to reach a 1% level of fluctuation in $`N`$ (lower horizontal line in the figure) we require a sample of about $`A6`$ sq. .deg. for the model with co-moving evolution and about $`A1`$ sq. deg. for a model with strong clustering evolution. In the submm survey of the HDF , with $`A=0.001`$ sq. deg. (the left vertical line in the Fig. 1), we expect $`\mathrm{\Delta }N/N0.2`$ in the case of weak clustering evolution and $`0.05`$ for the strong evolution case. Thus, in principle, comparing number counts in a few more submm surveys, similar in area to the HDF, would provide a clear discrimination between clustering evolution models.

Nevertheless in small submm surveys the number counts for bright sources will be low, e.g. $`N=\theta ^2𝒩3`$ for $`S_{850\mu \mathrm{m}}3`$ mJy in the HDF survey , and hence the shot-noise correction dominates and masks any evolution in the clustering. This situation is shown as a dot-dashed line in Fig.1 for the submm HDF survey. In this case, we could still subtract the shot-noise contribution using Eq. although this will introduce large uncertainties. Thus, in the future, it is necessary to extend the surveys to larger areas and/or to lower flux densities in order to have a discriminating measurement of clustering. The Gran Telescopio Milimétrico/Large Millimeter Telescope  will be the optimal facility to provide this new generation of deep and wide mm surveys.


text:
# Magnetic susceptibility, specific heat and dielectric constant of hexagonal YMnO3, LuMnO3 and ScMnO3

## I Introduction

Rare-earth manganites of the perovskite-type structure REMnO3 have been discovered in the 1950s. Materials with small ionic radius (RE = Ho, Er, Tm, Yb, Lu, Y and Sc) crystallize in the hexagonal structure \[S.G.*P63cm*\] while the compounds with larger ionic radius (RE = La, Ce, Pr, Nd, Sm, Eu, Gd, Tb or Dy) are orthorombic \[S.G.*Pnma*. The orthorhombic compounds (S.G. *Pnma*) form the basis for the CMR materials and they have been studied widely in the recent years. Unlike their orthorhombic counterparts, investigations in hexagonal REMnO3 compounds are few despite their interesting properties. Systems belonging to the hexagonal class undergo a ferroelectric transition far above the ordering of Mn3+ spins. Such magnetic ordering occurs in both structures, however, ferroelectricity is possible only for the hexagonal phase which has the noncentrosymmetric phase group *P63cm*. Although hexagonal manganates have been studied for many years , very recently there has been a resurgence of interest in these materials . Most hexagonal REMnO3 compounds can exhibit two transitions, i.e., a very high temperature (TFE $``$ 900 K) ferroelectric distortion and a low temperature ($``$ 100 K) magnetic ordering. The development of two order parameters is a rarity in oxides  and opens the possibility for electric-magnetic interactions and ”tuning” with electric and magnetic fields.

Current efforts have focused on the magnetic transition and refining the magnetic structure . These investigations have been complicated by the formation of a ferromagnetic impurity phase, Mn3O4, with a Curie temperature TC $``$ 43 K. Amounts of order of 0.5 at.% of this phase are found in most samples reported to date , causing significant anomalies in the bulk properties. For example, the true REMnO3 antiferromagnetic transition is masked by the ferromagnetic signal in the susceptibility and the specific heat may exhibit a second strongly magnetic field dependent peak around 40 K with a tiny entropy content.

In order to examine the intrinsic magnetic transition, we have developed a preparation technique for fabrication of pure samples (impurity content less then 0.1 at.%) of YMnO3, LuMnO3 and ScMnO3. We have measured the magnetic susceptibility ($`\chi `$) up to 400 K in different external fields. Specific heat (Cp) and its field dependence were also determined. Using a two Debye - temperature model we can deduce the excess $`\mathrm{}`$Cp and the magnetic entropy. Finally, we have measured the dielectric constant ($`ϵ`$). Our results show that there is a single antiferromagnetic transition at TN in each of these compounds. However, there appears a continuous increase in $`\chi `$(T) and a broad low T maximum ($``$ 50 K) in $`\mathrm{}`$Cp/T with significant entropy as the temperature is reduced. We attribute these anomalous features at T $`<`$ TN to the frustration of Mn3+ spins on a triangular lattice. From our measurements of the dielectric response of our high purity samples we find a more enhanced ”S”-shaped curve than the one reported in a previous investigation . This signature of TN in the dielectric susceptibility at TN $``$ TFE we believe is caused by a small change in the ferroelectric domain wall mobility and not by a direct coupling between the magnetic and ferroelectric order parameter.

## II Sample preparation and Experimental techniques

The samples investigated were prepared by the solid-state reaction technique at ambient pressure. Cation oxides of Y2O3, Lu2O3, Sc2O3 (99.99%) and MnO2 (99.99%), obtained from Alpha Aesar , were mixed in a 1:2 molar ratio to achieve the stochiometry of hexagonal LnMnO3. The mixture was well ground and calcinated in O2-flow at 11000C for 24h. To ensure better homogeneity the mixtures were ground again, compacted into small pellets and reheated in air at 14000C for 48 h for YMnO3 and LuMnO3; and at 9000C in flowing O2 for ScMnO3. It is only when this high-temperature ”reheating” is done for YMnO3 and LuMnO3, and the exact procedures given in ref.  are followed for ScMnO3 that we obtain impurity-free samples. This optimum procedure unfortunately depends slightly on the starting materials but the magnetic and structural properties of these compounds seem not to be sensitive to the O2 content. We may have a small $`\delta `$ missing/excees in the oxygen stochiometry such as REMnO3±δ.

The powder x-ray diffraction data on all three compounds show that the polycrystalline samples have the correct hexagonal structure without any trace of impurity phases. The lattice parameters found for each compound are in agreement with previous reports . The radius of the pellets is 1.5 mm and their thickness varies from 0.5 mm to 0.8 mm.

Dc-susceptibility ($`\chi `$) was measured in a Quantum Design magnetometer (MPMS 5S) for each sample using two fields, 0.1 T and 2 T in the temperature range from 1.8 to 400 K. Heat-capacity measurements were performed using a Quantum Design PPMS system in the range from 1.8 to 250 K and in fields up to 9 T. A relaxation technique was used with a resolution of 2% and an accuracy better than 5%. All samples showed only *one* sharp anomaly in the specific heat at TN without any additional peaks that could be ascribed to an impurity phase. Comparing with the specific heat of the most common impurity (Mn3O4) this leads to an estimation of less than 0.05% of this impurity. We also estimated the amount of second phase in all samples by measuring magnetization in different fields, which leads to even lower amounts. Finally we measured the dielectric constant $`ϵ`$ as function of the temperature using a capacitance bridge operating at a frequency of 1 MHz.

## III Susceptibility

Magnetic susceptibility measurements on the hexagonal compounds previously investigated failed to show any clear anomaly at TN in the $`\chi `$ versus $`T`$ curve. In 1, 2 and 3 we show the results on our high purity samples. The anomaly at TN is evident from the graphs. Reciprocal susceptibilities are also shown as an inset of these figures.

As expected for triangular antiferromagnets the susceptibility at the lowest temperature (T$`<`$ TN) does not decrease to 2/3 of its value at TN as for the classical two-sublattice case. From the earliest neutron diffraction combined with group theory arguments for this hexagonal lattice, Berthaut et al.  concluded that the Mn magnetic moments are directed in the basal plane and are oriented with angles of 120o with respect to each other. If the field is perpendicular to the magnetic plane all spins will be canted out from the plane with a small angle $`\delta `$ in the direction of the applied field. For this field direction the zero-temperature susceptibility will be equal to $`\chi (T_N)`$. A simple calculation with classical spins arranged with 120o angles, shows that $`\chi (T=0)`$ is at least $`\chi (T_N)`$ for all directions of the field in the basal plane. A detailed calculation will be carried out for this intriguing behavior.

The insets of the figures 13 exhibit the reciprocal susceptibility versus temperature. By fitting with a Curie-Weiss law in the temperature range from 300 to 400 K the values of the effective Mn moments ($`\mu _{eff}`$) and the paramagnetic Curie temperatures ($`\theta _{CW}`$) could be obtained. I collects these parameters together with TN, taken as the temperature where the kink in $`\chi `$ versus T occurs.

Because of the large values of $`\theta _{CW}`$ the temperature range for the fits (300 K - 400 K) is still too low. This explains the slight discrepancies between the obtained effective moments and the expected value of the $`Mn^{3+}`$(4.9$`\mu _B`$). Also the absolute values of $`\theta _{CW}`$ are much larger then the ordering temperature $`T_N`$. Ramirez  defined the ratio between these parameters as the measure of the geometrical frustration of the antiferromagnetic system. For our in-plane triangular lattice these value appears to be of order 5 which puts our hexagonal system in the class of moderately frustrated systems.

Note that a ferromagnetic impurity (Mn3O4) content of less then 0.05% will be immediately visible from the susceptibility measurements in two fields (not shown), e.g. 0.1 and 2 T. Our samples do not show any sign of the presence of Mn3O4 meaning that its abundance is less then 100 ppm.

## IV Specific heat

The results of the specific heat measurements for all three compounds in 0 and 9 T magnetic fields are shown in the insets of the 4, 5 and 6 together with the estimates for their lattice contribution (dash lines). The lattice specific heat was estimated using two Debye functions for the ”heavy”(Mn+RE) and ”light”(O2) elements in our compounds. The fits were carried out considering the specific heat at temperatures above 1.5 times $`T_N`$. We found reasonable values for the Debye temperatures, i.e., $`\theta _D`$(”light”)$``$ 780 K - 800 K and $`\theta _D`$(”heavy”) $``$ 350 K - 420 K. Subtraction of this estimation of the lattice specific heat provides us the (excess) magnetic contribution of the system. These $`\mathrm{\Delta }`$Cp are shown as $`\mathrm{\Delta }`$Cp/T versus T curves in 45 and 6 (left y-axis). On the right ordinate we plot the entropy computed from the integral. All samples exhibit sharp anomalies at the corresponding transition temperature indicating the onset of antiferromagnetic order. Up to an applied field of 9 T the specific heat is independent of the field. Contrary to previous work  we could not detect a second lower temperature anomaly, which becomes visible when the sample is contaminated by more then 0.5% of impurity (Mn3O4). The resulting entropy is in excellent agreement with that expected for a Mn3+ ion (S=2).

$$\mathrm{\Delta }S=Rln(2S+1)=13.38J/molK$$

In the literature  discrepancies between the deduced and expected entropy up to a factor of 10 have been reported due to an inaccurate estimate of the lattice contribution. Our analysis gives values for $`\mathrm{\Delta }`$S of at least 90% of the expected one. Note that the magnitude of $`\mathrm{\Delta }`$Cp at $`T_N`$ is approximately the same for all three compounds (18 J/mol K).

Furthermore, there is an the interesting anomaly (”bump”) at low T ($`50K`$) in the excess specific heat ( $`\mathrm{\Delta }`$Cp/T). The amount of entropy removed directly at $`T_N`$ is less than half of the total and this relative amount decreases going from Y via Lu to Sc. The physical cause of this ”bump” in $`\mathrm{\Delta }`$Cp/T and the corresponding entropy involved is related, in our opinion, to the frustrated magnetism. Here it seems that there are still spin degrees of freedom available far below TN and the spin system is progressively searching for its true ground state. We believe the ”bump” is correlated with the constantly increasing susceptibility below TN (see figures 4 \- 6).

Our observations are in agreement with different magnetic structures found via neutron diffraction by Muñoz et al. , via Raman-scattering and infrared-active spectra studied by Iliev et al.  and via non-linear optical spectroscopy by Fröhlich et al. . Thermal expansion measurements on single crystals are in progress for these hexagonal frustrated systems in order to give us a better picture of this ”hidden” processes that occur at lower temperatures.

## V Dielectric susceptibility

In general, one expects a coupling between ferroelectricity and magnetism only when TFE is close to TN. However, a coupling between the ferroelectric and antiferromagnetic order parameters in YMnO3 has been recently claimed by Huang $`\text{e}tal.`$  based on the observation of a small ”S”-shaped anomaly ($``$2%) at the Neel temperature in the $`ϵ`$ versus $`T`$ curve. They have suggested that such a coupling could arise due to the magnetostrictive effect associated with a lattice change accompanying the antiferromagnetic ordering of Mn3+ spins. We also have measured the dielectric response of our high purity samples and our data reveal a more enhanced ”S”-shaped curve - see 7. Note that the dielectric susceptibility, $`ϵ`$(T) is probed at TN $``$ TFE and that the step at TN is rather small. These leads us to believe that we are observing a magnetic ordering effect on the mobility of the domain walls of the ferroelectric domains and not a direct coupling between the two order parameters.

## VI Discussions

For samples containing small amounts of impurity phase (Mn3O4 with TC = 43 K) strong anomalies occurs around 43 K. For example, sharp spikes with small entropy contents appear in the specific heat  and there are strong ”upturns” and field dependence in susceptibilities. Our specially prepared samples (see II), whose data are shown in the previous figures, are completely free of such anomalies. Therefore, we believe our samples are without the troublesome ferromagnetic impurity phase and the claims of unusual behavior around 40 K by other authors are the result of impure samples.

There is remarkable difference in the bulk properties of YMnO3 and LuMnO3 on one hand and those of ScMnO3 at the other. The magnetic susceptibility of the former two compounds shows the characteristic kink at TN, while the susceptibility of ScMnO3 increases much more rapidly as temperature is decreased through the Neel point- see 8. Also, the magnetic susceptibility of ScMnO3 is sensitive to magnetic history of the sample, i.e., there is a difference between the field cooled and zero-field cooled values even for magnetic fields as low as 50 mT. Such field dependent history is not found in low field $`\chi `$ measurement for the YMnO3 nor for LuMnO3.

Considering the specific heat at the lowest temperatures, one sees that the total specific heat of ScMnO3 is proportional to T2 up to a temperature of 30 K, while those of YMnO3 and LuMnO3 show a much stronger T dependence (see inset of 8\- logarithmic scale). We would expect that the low temperature specific heat of an antiferromagnetic material can be described by a T3 term due to phonons and antiferromagnetic spin-waves (possibly with an exponential term due to a gap at zero wave-vector in the spin-wave spectrum). Now, a T2 dependence at low temperatures (below 30 K) shows that the specific heat is dominated by magnetic effects rather then by phonon excitations. Such a magnetic contribution which deviates from the antiferromagnetic spin-wave type indicates a ferromagnetic origin where the specific heat due to spin waves should follow T3/2.

These subtle differences in the magnetic behavior are the subject of our neutron diffraction and muon spin rotation experiments on polycrystalline and single crystal materials of these fascinating Mn-oxides.

## VII Conclusions

We have synthesized for the first time polycrystalline samples of the hexagonal REMnO3 compounds with a sufficient purity such that the antiferromagnetic magnetic ordering is revealed in the magnetic susceptibility. We have measured the specific heat and we show that an adequate analysis leads to the full entropy of the magnetic Mn-ions. We find that a considerable part of this entropy is released only at temperatures of order of half of the Neel temperature. There appears to be a remarkable difference between the properties of YMnO3 and LuMnO3 on the one hand and those of ScMnO3 on the other. Finally, we have measured the dielectric susceptibilities, which shows a kink at the respective magnetic transition temperatures. We believe this to be due to an influence from the antiferromagnetic ordering on the mobility of the ferroelectric domain-walls.

The authors gratefully acknowledge R.W.A. Hendrikx, and T.J. Gortenmulder for their technical assistance in analyzing the samples, H.B. Brom and I.G. Romijn for the use of their dielectric susceptibility apparatus and R. Feyerherm for valuable discussions. This work was supported by the Dutch Foundation FOM.


text:
# Evidence for softening of first-order transition in 3𝐷 by quenched disorder

## Abstract

We study by extensive Monte Carlo simulations the effect of random bond dilution on the phase transition of the three-dimensional 4-state Potts model which is known to exhibit a strong first-order transition in the pure case. The phase diagram in the dilution-temperature plane is determined from the peaks of the susceptibility for sufficiently large system sizes. In the strongly disordered regime, numerical evidence for softening to a second-order transition induced by randomness is given. Here a large-scale finite-size scaling analysis, made difficult due to strong crossover effects presumably caused by the percolation fixed point, is performed.

The influence of random, confining geometries on first-order phase transitions has been the subject of exciting experimental studies in the past few years. The case of the isotropic to nematic transition of $`nCB`$ liquid crystals confined into the pores of aerogels consisting of multiply connected internal cavities has been particularly extensively studied and led to spectacular results: The first-order transition of the corresponding bulk liquid crystal is drastically softened in the porous glass and becomes continuous , an effect which was not attributed to finite-size effects but rather to the influence of random disorder.

The first attempt to reproduce such a softening scenario using Monte Carlo (MC) simulations was reported by Uzelac et al.  who studied a three-dimensional (3D) $`q`$-state Potts model with spin variables (taking $`q=3`$ and $`4`$ states per spin) located inside the randomly connected pores of an aerogel modeled by diffusion-limited cluster aggregation. Although in experimental studies  commonly random fields or random uniaxial anisotropies are suggested to explain the softening of the transition, the random disorder chosen in Ref.  is coupled to the energy density and thus more akin to bond-dilution.

The qualitative effect of random bond disorder on second-order phase transitions is well understood through the Harris relevance criterion , and a beautiful experimental confirmation was reported in a LEED investigation of a two-dimensional (2D) order-disorder transition . For systems with a first-order transition in the pure case, randomness generically softens the transition and, under certain circumstances, may even induce a second-order transition according to a picture first proposed by Imry and Wortis .

In 2D, the natural candidate for theoretical investigations is the $`q`$-state Potts model, since this model is exactly known to exhibit regimes with first- and second-order transitions , depending on the value of $`q`$. Many results were obtained in both regimes in the last ten years , including approximate analytic treatments, MC simulations, transfer-matrix calculations, and high-temperature series expansions. Among others also quite intricate problems such as self-averaging and multi-fractality have recently been studied in some detail .

In 3D, to date only the Ising model with site dilution has been studied extensively . In accordance with the Harris criterion the presence of random disorder was found to modify the critical exponents to values close to $`\nu =0.684(5)`$, $`\gamma /\nu =1.96(3)`$, and $`\beta /\nu =0.519(8)`$. Concerning the influence of random disorder on first-order transitions, even less is known in 3D. Apart from the exploratory work  of Uzelac et al. , only the site diluted 3-state Potts model, which in the pure 3D case has a very weak first-order transition , has recently been studied via large-scale MC simulations . This study led to the conclusion that the critical exponent $`\nu `$ governing the scaling behavior of the correlation length is compatible with that of the 3D site diluted Ising model, whereas the $`\eta `$ exponent is definitely different.

The purpose of this Letter is to present numerical evidence for softening of the transition when it is strongly of first order in the pure system, in order to be sensitive to disorder effects. The paradigm in 3D is the 4-state Potts model, since the correlation length at the transition temperature of the unperturbed system is small enough ($`\xi 3`$ in lattice spacing units ) to allow simulations of significantly large systems. For the pure 3D 5-state Potts model the first-order transition is already too strong . In the following we, therefore, consider the 4-state bond diluted Potts model on simple-cubic lattices of size $`V=L^3`$ with periodic boundary conditions. The Hamiltonian of the system with independent, quenched random interactions is written as $`\beta =_{(i,j)}K_{ij}\delta _{\sigma _i,\sigma _j}`$, where the spins take the values $`\sigma _i=1,\mathrm{},4`$ and the sum goes over all nearest-neighbor pairs $`(i,j)`$. The coupling strengths are allowed to take two different values $`K_{ij}=KJ/k_BT`$ and $`0`$ with probabilities $`p`$ and $`1p`$, respectively. The order parameter for a given realization of the $`K_{ij}`$ is defined by the majority orientation of the spins, $`m=\mu `$, where $`\mu =(q\rho _{\mathrm{max}}1)/(q1)`$ and $`\rho _{\mathrm{max}}`$ is the maximum value of the density of spins $`\rho _\sigma `$ in the $`q=4`$ possible spin states. The thermal average over the MC iterations is indicated by brackets $`\mathrm{}`$, and the physical quantities are then averaged over disorder realizations, e.g., $`\overline{m}=\overline{\mu }`$. For each disorder realization the susceptibility is defined as usual via the fluctuation-dissipation theorem, $`\chi =KV(\mu ^2\mu ^2)`$.

The present MC simulations consist of two parts. First, a scan of the dilution-temperature plane in order to determine the phase diagram, and second, a large-scale finite-size scaling (FSS) study at $`p=0.56`$ up to $`L=96`$ in the dilution regime exhibiting second-order transitions. The spin updates were performed with the cluster-flipping method  in the Swendsen-Wang formulation which turned out to be better behaved than the Wolff single-cluster version for high dilutions (small $`p`$), where small clusters connected by non-vanishing bonds are more likely to appear. We made sure that, by adjusting the length of the runs of cluster algorithms, at least 10 tunnelling events between the two coexisting phases of the pure system ($`p=1`$) were observed up to lattice size $`L=16`$. To improve the accuracy, for weak dilutions ($`p`$ between 1 and 0.68) where we obtained evidence for first-order transitions, we also used multicanonical algorithms . For the determination of the maxima of observables we applied standard histogram reweighting techniques in order to extrapolate the results over a temperature range around the simulation point.

At each probability $`p`$ and for each realization of the random couplings, between $`15\times 10^3`$ to $`30\times 10^3`$ MC sweeps per spin were performed, resulting in at least 250 (almost) independent measurements of the physical quantities for the largest lattice size considered. This turned out to be sufficient for reliable thermal averages. For the average over disorder realizations, between 2 000 and 5 000 samples were generated.

The phase diagram is determined from the locations of the maxima of the average susceptibility, $`K_{\mathrm{max}}`$, obtained for systems of increasing sizes up to $`L=16`$. Theoretically we expect that the transition remains of first order in the regime of low impurity concentrations. For increasing concentrations a regime of second-order transitions should appear from a (tri-) critical concentration until the percolation threshold is reached where the transition vanishes altogether. After this percolation threshold, no ordered phase can exist at any finite temperature. Two points are known in the $`pk_BT/J`$ plane: The transition temperature $`T_t`$ of the pure system  $`k_BT_t(p=1)/J1.59076`$, and $`k_BT_t(p=p_c)/J=0`$ at the percolation threshold $`p_c0.2488`$. As can be inspected in Fig. 1, the temperatures of the susceptibility maxima for different lattice sizes are very stable and already for $`10^3`$ spins an accurate transition line is obtained. Also shown is the result of a simple mean-field argument taking into account the average number of neighbors, $`k_BT_t(p)/\overline{z}J=\mathrm{const}`$, where the constant is chosen such that $`T_t`$ of the pure system is reproduced. This leads to a simple linear approximation of the transition line, $`k_BT_t(p)/J=1.59076\times p`$, which is surprisingly accurate over a significant range of $`p`$-values.

By monitoring the FSS behavior of various thermodynamic quantities as well as the (pseudo-) dynamics of the update algorithm, we estimate the tricritical point to be located between $`p=0.68`$ and $`0.84`$: The shape of the energy probability densities as well as the Binder cumulants suggest a first-order transition at $`p=0.84`$ and above while a clear second-order signal is observed at $`p=0.68`$ and below. For the investigation of the critical properties in the second-order regime we have chosen a dilution $`p=0.56`$ where the corrections to scaling are seemingly relatively small, since the effective transition temperatures corresponding to the susceptibility maxima remain almost constant in the range of sizes $`2L16`$ used for the determination of the phase diagram.

In order to convince ourselves that for $`p=0.56`$ the transition is indeed of second order, let us first consider the average probability densities of the energy, $`P(e)`$. In Fig. 2 their shapes close to $`K_{\mathrm{max}}`$ are depicted for various lattice sizes up to $`L=96`$. We see that the system exhibits for small sizes two distinct peaks which clearly collapse into a single peak when one approaches the thermodynamic limit. This is precisely what is expected at a second-order phase transition, while in the case of a first-order transition the double-peak structure should persist for all sizes and, in fact, should become even more pronounced when the system size increases. Physically, a two-peak structure would reveal the presence of two phases at the transition temperature, and this coexistence is the characteristic feature of a first-order transition. More quantitatively we have confirmed that the effective interface tension $`\sigma _{\mathrm{od}}`$ derived from the usual relation $`P_{\mathrm{min}}/P_{\mathrm{max}}\mathrm{exp}(2\sigma _{\mathrm{od}}(L)L^{d1})`$ vanishes in the infinite-volume limit; see the inset of Fig. 2.

Our final goal is a quantitative characterization of the critical behavior by providing estimates for the critical exponents of the transition. To this end we have performed a standard FSS analysis at $`p=0.56`$. As can be inspected in Fig. 3, the corrections to asymptotic FSS for $`\overline{\chi }_{\mathrm{max}}`$ seem to become quite small above $`L=30`$. The data are thus linearly fitted to $`a_\chi L^{\gamma /\nu }`$ for $`L`$ in the ranges $`L_{\mathrm{min}}>30`$ to $`L_{\mathrm{max}}=96`$, and the resulting exponents are collected in Table I. Selecting the fits with the smallest chi-squared per degree of freedom, $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}`$, we take as the final result the lines in bold face, e.g., $`\gamma /\nu =1.50\pm 0.02`$.

The quantity $`_K\mathrm{ln}\overline{m}`$ gives an estimation of the exponent $`\nu `$, $`(_K\mathrm{ln}\overline{m})_{K_{\mathrm{max}}}L^{1/\nu }`$. Here our analysis leads to an estimate of $`1/\nu =1.33\pm 0.03`$ or $`\nu =0.752\pm 0.014`$, in agreement with the stability condition of the random fixed point ($`\nu 2/D=0.666\mathrm{}`$) and significantly different from the estimate for the site diluted 3D Ising model ($`\nu =0.684(5)`$). The same procedure was applied to the magnetization evaluated at the temperature where the susceptibility is maximal. Judging the values of $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}`$ would lead to the result given in bold face in Table I, but the effective exponent for the magnetization is clearly not yet stable. We therefore also considered the FSS behavior of higher (thermal) moments of the magnetization, $`\overline{\mu ^n}`$ which should scale with a dimension $`n\beta /\nu `$. The results for the first moments exhibit, however, again much stronger corrections to scaling than we observed for $`\overline{\chi }`$ or $`_K\mathrm{ln}\overline{m}`$, leading to our final estimate of $`\beta /\nu =0.65\pm 0.05`$. .

From the log-log plots of the three quantities in Fig. 3 one can clearly observe a crossover from a percolation-type behavior at small sizes, characterized by the exponents  $`\gamma /\nu 2.05`$, $`1/\nu 1.124`$, and $`\beta /\nu 0.45`$, towards a new regime at large sizes, which presumably corresponds to the random fixed point, with exponents as given above. The numerical evidence for this interpretation is quite striking, but with the present system sizes we can of course not completely rule out the possibility of corrections to scaling, in particular since for the 3D disordered Ising model it is well known that such corrections at the random fixed point are strong (with a correction-to-scaling exponent around $`\omega =0.4`$). In order to investigate this question for the 3D 4-state Potts model, we tried to fit the physical quantities to the standard expression, e.g., $`a_\chi L^{\gamma /\nu }[1+b_\chi L^\omega +\mathrm{}]`$, including a sub-dominant correction-to-scaling term. Since 4-parameter non-linear fits are notoriously unstable, we performed linear fits where the exponents are kept fixed and only the amplitudes are free parameters. In Fig. 4, we show a 3D plot of the total $`\chi ^2`$ for the susceptibility fits as a function of $`\gamma /\nu `$ and $`\omega `$. We observe a clear, stretched valley which confirms that $`\gamma /\nu `$ is close to $`1.5`$, but obviously this does not allow any reliable estimation of the correction-to-scaling exponent $`\omega `$. The same procedure for $`1/\nu `$ gives qualitatively the same picture and confirms our previous estimate of $`1/\nu =1.33`$. For $`\beta /\nu `$, on the other hand, the $`\chi ^2`$-landscape turns out to be very flat and extremely sensitive to the fit range.

To conclude, from large-scale MC simulations of the 3D 4-state bond diluted Potts model we obtained a) the phase diagram in the $`pT`$ plane in very good agreement with (rescaled) mean-field theory, b) the approximate location of the tricritical point around $`p_{TCP}=0.76(8)`$, and c) at the dilution $`p=0.56`$ clear evidence for softening of the rather strong first-order phase transition in the pure case towards a continuous transition with estimates for the critical exponents of $`\nu =0.752(14)`$, $`\gamma /\nu =1.50(2)`$, $`\gamma =1.13(4)`$, $`\beta /\nu =0.65(5)`$, and $`\beta =0.49(5)`$. These are clearly different from the values for both the disordered 3D Ising \[$`\nu =0.684(5)`$, $`\gamma /\nu =1.963(5)`$\] and the $`3`$state Potts \[$`\nu =0.690(5)`$, $`\gamma /\nu =1.922(4)`$\] models.

###### Acknowledgements.

We gratefully acknowledge financial support by the DAAD and EGIDE through the PROCOPE exchange programme. C.C. thanks the DFG for financial support through the Graduiertenkolleg “Quantenfeldtheorie” in Leipzig. Work partially supported by the computer-time grants hlz061 of NIC, Jülich, C2000-06-20018 of the Centre Informatique National de l’Enseignement Supérieur (CINES), and C2000015 of the Centre de Ressources Informatiques de Haute Normandie (CRIHAN).


text:
# Universal quantum processors with arbitrary radix 𝑛≥2

(March 2001)

## Abstract

Here is discussed the Hamiltonian approach to construction of deterministic universal (in approximate sense) programmable quantum circuits with qubits or any other quantum systems with dimension of Hilbert space is $`n2`$.

Let us suggest, that universal quantum processor can be described as circuit with three “buses”: quantum, intermediate, and pseudo-classical. It is discussed with more details elsewhere (in  based on results about quantum control  and programmability ). Here are mentioned only topics related with theme of present paper about universality.

If we have finite set of universal (in approximate sense ) quantum gates, we may approximate any unitary operator with necessary precision . There exist sets of two-gates both for binary quantum circuits with qubits  and for non-binary ones with equal dimensions $`n>2`$ of Hilbert space of each elementary system (sometime it is called qu$`n`$it) .

Here is considering application of universal set of quantum gates in three-level design of quantum processor described above. Intermediate bus together with quantum data bus is example of quantum control , i.e., if state of intermediate bus is $`|p`$ then $`p`$-th gate $`U_p`$ from given set is applied to quantum data bus, $`|dU_p|d`$.

An approximation of given unitary operator $`U`$ is decomposition on product of universal gates $`U_k`$:

$$U=U_{p_l}\mathrm{}U_{p_2}U_{p_1}$$
(1)
And so it is necessary to apply quantum control $`l`$-times with value of intermediate bus is $`|p_l`$ on $`l`$-th step .

For generation of the sequence of numbers $`p_l`$ may be used  reversible111It can be shown, that using irreversible classical gates on intermediate and pseudo-classical buses together with inevitable problem with quantum description of composite system with three buses may also cause non-unitary evolution of quantum data bus classical (in terminology used by ) circuits acting on pseudo-classical bus together with intermediate bus. Due to this property of quantum processor it may be important sometime to have well defined algorithm of decomposition Eq. (1) for arbitrary unitary operator $`U`$. Otherwise we should keep all $`l`$ numbers as table and then size of pseudo-classical bus is proportional to maximally possible length $`l_{\mathrm{max}}`$ of simulation (“cyclic Q-ROM approach” ).

An approach with clear algorithm of decomposition uses structure of Lie algebra of unitary group  and from physical point of view related with Hamiltonian of quantum circuits. In such approach main object is not a unitary matrix $`U`$, but Hamiltonian of the evolution, i.e. Hermitian matrix $`H`$ with property

$$U(\tau )=e^{iH\tau },t,U=U(t)$$
(2)
for some real parameters $`t,\tau `$, time.

The exponent Eq. (2) makes possible for small $`\tau `$ to use approximation with sums instead of products. Say if Hermitian matrix $`H`$ can be decomposed on some basis $`H_I`$ in space of Hermitian matrices, then

$$U(\tau )=\mathrm{exp}(i\underset{I}{}n_IH_I\tau )\underset{I}{}\mathrm{exp}(in_IH_I\tau ),\tau 0.$$
(3)
Let $`U_I=\mathrm{exp}(iH_I\tau )`$ and $`n_I`$ may be approximated by natural number (integral part) if $`\tau `$ is small and $`n_I`$ are big. Then:

$$\mathrm{exp}(i\underset{I}{}n_IH_I\tau )\underset{I}{}U_I^{n_I},UU(t)\left(\underset{I}{}U_I^{n_I}\right)^{t/\tau }$$
(4)
where $`t/\tau `$ may be considered as natural number or approximated by it for small $`\tau `$.

But number of elements of basis of Hermitian matrix for quantum circuit is exponentially large. Say for system with $`k`$ qubits it is $`4^k`$ and it is $`n^{2k}`$ for other radix $`n`$. Due to such property, as Hamiltonians of universal quantum gates are chosen only small subset with property, that any other element of basis can be generated by using sequence of commutators

$$H_{[JK]}i(H_JH_KH_KH_J)$$
(5)
To produce $`U_{[JK]}=\mathrm{exp}(iH_{[JK]}\tau )`$ it is possible to use approximation:

$$U_{[JK]}U_J^{s_\tau }U_K^{s_\tau }U_J^{s_\tau }U_K^{s_\tau },s_\tau =1/\sqrt{\tau }$$
(6)
there for small $`\tau `$ $`s_\tau `$ again can be considered as a natural number (integral part). Here small indexes like $`k`$ are used for set of initial universal gates, $`U_1`$ and capital indexes like $`J`$ are compound, $`U_{[[12]3]}`$.

Here is a problem with negative power, $`s_\tau `$. It can be resolved for qubit by special design  where for any universal gate $`U_k(t_0)`$ is unit for same $`t_0`$ and so $`U_k^s=U_k^{t_0s}`$ with positive $`t_0s`$. For non-binary universal gates it is not necessary so, but it is possible simply to double amount of universal gates and have $`U_k^1`$ together with any gate $`U_k`$.

Let us describe now algorithm of approximation in general. It is chosen some small interval of time $`\tau `$, the less it the higher precision. For approximation of some gate $`U`$ with Hamiltonian $`H`$ as in Eq. (2) may be found linear decomposition of the $`H`$ by basis $`H_I`$ Eq. (3). Each component with compound indexes $`U_I`$ is represented by universal gates $`U_k`$ as it was shown in Eq. (6), and finally approximation of $`U`$ is result of nested cycles described by Eq. (3), Eq. (4).

The example shows, that approximation may demand many operations and number of steps quickly grows with refinement of precision. But if there is some method to implement structure of the approximating algorithm as reversible classical circuit, then size of program register (defined mainly by length of pseudo-classical bus) may be more appropriate. Due to it algebraic structures similar with introduced in  may be useful.

There is also important note about reversibility of algorithm of approximation. Formally it is possible to use some irreversible circuit and apply standard technique, i.e. instead of irreversible function $`f:af(a)`$ to use reversible one with property $`F:(a,0)(a,f(a))`$, but then each step of algorithm will produce new portion of “junk” and pseudo-classical bus again must have size proportional to $`l_{\mathrm{max}}`$, and such a case maybe even worst, than reversible cyclic Q-ROM register without any programming.

But even if such reversible algorithm is found, the property of universality may be rather formal. It is clear from consideration above, that number of steps is proportional to amount of non-vanishing terms $`H_I`$ in decomposition of Hamiltonian. So, if problem area is related with absolutely arbitrary unitary operators, then particular set of gates is not very essential and length of simulation is exponential on number of qubits.

Of course such situation would not be realistic if number of qubits is big enough and so set of basic universal gates should take into account particular set of possible problems. Say it may be task to simulate any possible quantum circuits composed by many different $`k`$-gates, with $`k`$ is not very large, i.e. “local” quantum circuits.