text:
# Neutrino Clustering in Cold Dark Matter Halos : Implications for Ultra High Energy Cosmic Rays

## I Introduction

The nature of cosmic rays above the Greisen-Zatsepin-Kuzmin (GZK) cutoff GZK  at $``$ 4 $`\times `$ 1019 eV is an unsolved problem in ultra high energy cosmic ray (UHECR) physics review . These events have been reported by the Akemo Giant Air Shower Array (AGASA) AG , Fly’s Eye FE , Havera Park HP , HiRes HR , and Yakutsk YK  collaborations. Interactions with the cosmic background photons $`\gamma _{cmb}`$ via photoproduction of pions ($`p\gamma _{cmb}p+N\pi ,n\gamma _{cmb}n+N\pi `$), photopair production ($`p\gamma _{cmb}pe^+e^{},\gamma \gamma _{cmb}e^+e^{}`$), and inverse Compton scattering ($`e^\pm \gamma _{cmb}e^\pm \gamma `$, $`p\gamma _{cmb}p\gamma `$) at high energies L94  constrain a $`10^{20}`$ eV cosmic ray to a few Mpc for the characteristic lengths of either charged cosmic rays or neutrons and photons. More specifically, the attenuation length of protons above the GZK cutoff is $``$ 50 Mpc. The lack of known processes to accelerate cosmic rays in small Galactic objects makes the Galactic origin of these ultra high energy particles unfeasible FS95 . Novel powerful acceleration mechanisms for light nuclei are required if these energetic particles are produced in nearby galaxies A02 . Exotic particles and dynamics have also been suggested, but these come with their own difficulties review .

One proposed explanation for the UHECRs is the Z-burst model, which tries to solve the puzzle without invoking new physics beyond the standard model of particle physics except for neutrino masses. Several recent experiments SK601 ; MAC01 ; SN01 ; TR01  have found evidence for non-zero neutrino mass. The Z-burst model hinges on the fact that ultra high energy neutrinos (and anti-neutrinos) produced at cosmological distances can reach the GZK zone unattenuated. Their resonant annihilation on the relic anti-neutrinos (and neutrinos) produces Z bosons, about 70% of which decay into hadrons within $`10^{25}`$ sec. The final state has fifteen pions and 1.35 baryon-antibaryon pairs on average FMS99 , where the fifteen pions decay into thirty high energy photons. The Z boson is highly boosted ($`10^{10}`$) W99 , resulting in a highly collimated beam with a half angle of $`10^{10}`$. This and the fact that the effect of magnetic fields at such high energies is negligible GMT01  ensure a high probability for the protons and photons to reach the observer if the Z-burst occurs in the direction of the Earth. The Z-burst model has been discussed in detail in many papers W82 ; R93 ; YSL98 ; GK ; CDF01 ; FKR01 ; FKR02 . The resulting cosmic ray flux has been shown to depend strongly on the density of the relic neutrinos YSL98 ; FMS99 ; W99 ; F02 ; FGSD01 , but the neutrino density in these calculations has been taken to be either the constant relic density from the big bang or some ad hoc value.

In this paper we perform a detailed calculation of the neutrino clustering in the local universe using realistic cosmological models and apply the results to the Z-burst model for UHECRs. Since the current constraints from cosmological observations and laboratory experiments indicate that the neutrino masses are small ($`<2`$ eV; see Sec II) and the CDM dominates the dark matter density ($`\mathrm{\Omega }_{\mathrm{cdm}}\mathrm{\Omega }_\nu `$), we do not expect the clustering of neutrinos to affect significantly that of the CDM. As a result, it is not essential to use full scale, time consuming $`N`$-body simulations. Instead, we solve the collisionless Boltzmann equation for the neutrino phase space distribution in a background potential given by the universal profile of CDM halos reported in recent high resolution simulations NFW96 . The Boltzmann equation is then linear in the neutrino density contrast and has tractable integral solutions. The advantage of this method over the conventional $`N`$-body simulations is that we can obtain the neutrino density profile much below the resolution scale ($`50`$ kpc) of large cosmological simulations by using as an input the CDM potential determined from much higher resolution simulations of individual halos. Moreover, the computation time required for our approach is negligible compared with numerical simulations, thereby allowing us to explore a large parameter space of neutrino masses and dark matter halo masses.

In Sec II the relevant Boltzmann equation and the integral solutions are derived. In Sec III results for the clustering of neutrinos for different neutrino masses and CDM halos are presented and compared with physical arguments based on neutrino free streaming and the Tremaine-Gunn constraint TG79 . The resulting neutrino density profiles are also compared with earlier $`N`$-body simulations KKPH96 , which show good agreement. In Sec IV the neutrino overdensity calculation is incorporated in the Z-burst model for UHECRs, where we make realistic predictions for the UHECR energy spectrum for different neutrino masses. We estimate the level of anisotropy in the UHECR flux by examining lines of sight towards five of the most massive clusters (Virgo, Perseus-Pisces, Hydra, Centaurus, and Coma) in the local universe (within 100 Mpc) where neutrinos are likely to be most clustered.

## II Boltzmann Equation for Neutrino Clustering in CDM Halos

In this section we develop an approach based on the collisionless Boltzmann equation to study how massive neutrinos cluster gravitationally in realistic cosmological models. We start by noting that the median velocity of unclustered background neutrinos of mass $`m_\nu `$ (in eV) at redshift $`z`$ is

$$\overline{v}=161(1+z)m_\nu ^1\mathrm{km}\mathrm{s}^1.$$
(1)
This implies that light neutrinos ($`m_\nu 2`$ eV) do not accrete significantly onto CDM protoclusters until $`z3`$ because the neutrino thermal velocities are greater than the velocity dispersion of a typical cluster or supercluster. We are then faced with the more tractable problem of how neutrinos cluster in the potential well of an existing CDM halo. Our approach is to use the collisionless Boltzmann equation for the neutrino phase space density $`f`$ and follow its evolution in a background CDM potential given by the approximate universal profile obtained in high resolution simulations of individual halos NFW96 . Note that the CDM potential is time-dependent in general. Earlier work has used the Boltzmann approach to study how neutrinos cluster around point masses in the context of cosmic string seeded galaxy formation BKT87 ; BW88 . This method allows us to calculate the neutrino density profiles in the inner part of the cluster ($`<10`$ kpc) that can not be resolved by large cosmological simulations. This will be seen to be important in the Z-burst model where a significant contribution to the cosmic ray flux comes from the inner regions of the halo.

In the Newtonian approximation and in physical coordinates, the collisionless Boltzmann equation takes the form

$$\frac{f}{t}+\dot{\text{r}}_rf+\dot{\text{p}}_pf=0.$$
(2)
Rewriting it in conformal time $`d\tau =dt/a`$ and in comoving position $`\text{x}=\text{r}/a`$ and momentum $`\text{q}=a\text{p}m_\nu \dot{a}\text{r}`$, we obtain

$`{\displaystyle \frac{1}{a}}{\displaystyle \frac{f}{\tau }}+{\displaystyle \frac{\text{q}}{m_\nu a^2}}_xfm_\nu \ddot{a}a\text{x}_qf`$
$`m_\nu _x\mathrm{\Phi }_qf=0,`$ (3)
where the Newtonian gravitational potential $`\mathrm{\Phi }`$ obeys $`\dot{\text{p}}=m_\nu \mathrm{\Phi }`$. At the time of decoupling the neutrino phase space density is given by the thermal Fermi-Dirac distribution $`f_0(q)(e^{q/T_{\nu ,0}}+1)^1`$ where $`T_{\nu ,0}=(4/11)^{1/3}T_{\gamma ,0}=1.676\times 10^4`$ eV is the temperature of the neutrino background today. Gravitational clustering distorts the spatially uniform $`f_0`$, so we write the full distribution as

$$f(\text{x},\text{q},\tau )=f_0(q)+f_1(\text{x},\text{q},\tau ).$$
(4)
The gravitational potential can also be written as

$$\mathrm{\Phi }(\text{x},\tau )=\mathrm{\Phi }_0(\text{x},\tau )+\mathrm{\Phi }_1(\text{x},\tau ),$$
(5)
where $`\mathrm{\Phi }_0`$ is related to the mean background comoving density $`\overline{\rho }=\overline{\rho }_{\mathrm{cdm}}+\overline{\rho }_\nu `$ by $`_x\mathrm{\Phi }_0=\frac{4\pi }{3}G\overline{\rho }a^2\text{x}`$, and $`\mathrm{\Phi }_1`$ is determined by the density contrast of both CDM and neutrinos in the halo:

$`_x^2\mathrm{\Phi }_1`$ $`=`$ $`4\pi Ga^2(\delta \rho _{\mathrm{cdm}}+\delta \rho _\nu )`$ (6)
Eq. (3) then becomes

$`{\displaystyle \frac{1}{a}}{\displaystyle \frac{f_1}{\tau }}+{\displaystyle \frac{\text{q}}{m_\nu a^2}}_xf_1m_\nu _x\mathrm{\Phi }_1_qf_0`$
$`m_\nu _x\mathrm{\Phi }_1_qf_1=0,`$ (7)
where we have used the Friedmann equation $`\ddot{a}=\frac{4\pi }{3}G\overline{\rho }a`$, which gives $`\ddot{a}a\text{x}+_x\mathrm{\Phi }_0=0`$. We note that eq. (7) is the full Boltzmann equation and no approximation has been made thus far. It is generally a nonlinear equation in $`f_1`$ where $`\mathrm{\Phi }_1`$ is related to $`f_1`$. For our problem, however, $`\mathrm{\Phi }_1`$ is mostly determined by the CDM whose potential has a well known, pre-determined form. Eq. (7) is then linear in $`f_1`$ and much easier to solve.

Furthermore, eq. (7) has a simple integral solution if the fourth term is neglected. For example, in earlier calculations that examined neutrino clustering onto point masses seeded by cosmic strings, this term was dropped in order to simplify the calculation BKT87 ; BW88 . We will also drop this term, but we justify this approach in two ways. First, we note that dropping this term requires $`_qf_1<_qf_0`$ and not $`f_1<f_0`$. The former is generally a less restrictive condition and can be satisfied even if $`f_1`$ is much larger than $`f_0`$ because on dimensional grounds, we have

$$\frac{_qf_1}{_qf_0}\frac{f_1/\sigma _v}{f_0/\overline{v}}\frac{\delta \rho _\nu /\sigma _v^4}{\overline{\rho }_\nu /\overline{v}^4}\delta _\nu \left(\frac{\overline{v}}{\sigma _v}\right)^4,$$
(8)
where $`\sigma _v`$ is the velocity dispersion of the gravitational potential $`\mathrm{\Phi }_1`$ and $`\overline{v}`$ is the median neutrino thermal velocity in eq. (1). Since only neutrinos with $`\overline{v}\sigma _v`$ are cold enough to fall into the gravitational wells, we expect the ratio $`_qf_1/_qf_0`$ to be much smaller than $`f_1/f_0`$, thereby making it easier to justify the dropping of the fourth term in eq. (7). (For example, $`_qf_1/_qf_0<0.2`$ for $`1`$ eV neutrinos in $`10^{14}M_{}`$ halos.) Eq. (8) further indicates that the neutrino overdensity is larger than $`f_1/f_0`$ by $`(\sigma _v/\overline{v})^3`$, a large factor in highly clustered regions. This explains qualitatively the large overdensity found in our numerical results to be presented in Sec III. In the next section we also provide further justification for ignoring the $`_qf_1`$ term by comparing our results with the neutrino density profiles of two halos obtained from earlier numerical simulations.

Following BKT87 , we convert eq. (7) into an ordinary differential equation by going into Fourier space and using a new time variable $`d\eta =d\tau /a`$:

$$\frac{\stackrel{~}{f_1}}{\eta }+\frac{i\text{k}\text{q}}{m_\nu }\stackrel{~}{f_1}+\frac{im_\nu }{k^2}4\pi Ga^4\stackrel{~}{\rho }\text{k}_qf_0=0,$$
(9)
where $`\stackrel{~}{f_1}`$ and $`\stackrel{~}{\rho }`$ are the Fourier transforms of $`f_1`$ and $`\rho `$. The solution is

$`\stackrel{~}{f}_1(\text{k},\text{q},\eta )`$ $`=`$ $`{\displaystyle \frac{im_\nu }{k^2}}4\pi G{\displaystyle _{\eta _0}^\eta }𝑑\eta ^{}e^{i\text{k}\text{q}(\eta \eta ^{})/m_\nu }`$ (10)
$`\times a^4(\eta ^{})\stackrel{~}{\rho }(\text{k},\eta ^{})\text{k}_qf_0,`$
where we have taken the initial neutrino phase space to be homogeneous, i.e., $`f_1(\eta _0)=0`$ and $`f(\eta _0)=f_0`$.

The comoving neutrino number density is given by

$$\stackrel{~}{n}_\nu (\text{k},\eta )\overline{n}_\nu =\frac{2}{h_p^3}d^3\text{q}\stackrel{~}{f_1}(\text{k},\text{q},\eta ),$$
(11)
which can be obtained from eq. (10) after integration by parts in $`q`$ and integration over angles:

$`\stackrel{~}{n}_\nu (\text{k},\eta )\overline{n}_\nu `$ $`=`$ $`{\displaystyle \frac{32\pi ^2Gm_\nu }{h_p^3k}}{\displaystyle _{\eta _0}^\eta }𝑑\eta ^{}a^4(\eta ^{})\stackrel{~}{\rho }(k,\eta ^{})`$ (12)
$`\times {\displaystyle _0^{\mathrm{}}}dqq{\displaystyle \frac{\mathrm{sin}[kq(\eta \eta ^{})/m_\nu ]}{e^{q/T_{\nu ,0}}+1}},`$
where $`\overline{n}_\nu 112`$ cm-3 is the cosmic mean comoving number density of one species of light neutrinos and antineutrinos, and $`h_p`$ is the Planck constant. Eq. (12) is the main equation that we will solve in this paper. It is a Volterra integral equation of the second kind that has the form $`f(t)=_a^t𝑑sK(t,s)f(s)+g(t)`$ and can be solved with the trapezoidal rule. It describes how neutrinos of a given mass $`m_\nu `$ cluster around a realistic CDM halo as a function of time. The density of the halo at a given time, $`\stackrel{~}{\rho }`$, should generally be the sum of the CDM and the neutrino components, but as we verify below, using the CDM potential alone is a very good approximation for the cosmological models of current interest.

## III Results for Neutrino Clustering in CDM Halos

In this section we present results for $`n_\nu `$ computed from eq. (12). We choose to integrate eq. (12) from $`z=3`$ to 0; the results differ by only about 5% if the initial redshift is pushed back to 5 because the neutrinos do not cluster appreciably at such early times as discussed above. We also need to specify the cosmological models and neutrino masses. For the cosmological parameters, we use the currently favored critically-flat model with matter density $`\mathrm{\Omega }_m=\mathrm{\Omega }_{\mathrm{cdm}}+\mathrm{\Omega }_\nu =0.35`$, cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$, and Hubble parameter $`h=0.7`$. Variations in these parameters at 10 to 20% level are not expected to alter our neutrino results significantly since the effect on the halo potential $`\mathrm{\Phi }_1`$ is small. For the neutrino masses, we consider four different models in which three models assume three degenerate massive species with masses 0.6, 0.3, and 0.15 eV respectively, and one model with one massive species with mass 1.8 eV. The corresponding density parameter in neutrinos is $`\mathrm{\Omega }_\nu 0.04,0.02,0.01`$ and 0.04, respectively, all much smaller than $`\mathrm{\Omega }_{\mathrm{cdm}}`$. This range of neutrino masses is chosen to span the current cosmological and laboratory limits. The most recent cosmological constraint comes from the galaxy clustering power spectrum of 2-Degree-Field Redshift Survey, which places an upper limit of 1.8 eV on the sum of the neutrino masses E02 . The Super-Kamiokande experiment SK601  provides strong evidence for oscillations between neutrino species with a mass difference of $`\delta m^2=(18)\times 10^3`$ eV, giving a minimum mass of $`0.07`$ eV if the neutrino masses are hierarchical. The choice of three degenerate neutrino masses is based on indications that if any of the mass eigenvalues is above 0.1 eV then all three masses are above 0.1 eV and almost degenerate BB02 .

Although the final neutrino density profile will depend strongly on the CDM profile, we do not expect the reverse to hold: CDM density will not be much affected by the clustering of neutrinos because all the models considered in this paper have small $`\mathrm{\Omega }_\nu `$/$`\mathrm{\Omega }_{\mathrm{cdm}}`$, so the CDM dominates the gravitational potential of a dark matter halo. We are therefore justified in using the universal CDM profile determined from the high resolution pure CDM simulations NFW96  as the input:

$$\rho _{\mathrm{cdm}}(r)=\frac{\overline{\rho }\overline{\delta }r_s^3}{r(r+r_s)^2},$$
(13)
where $`\overline{\delta }`$ and $`r_s`$ are given in terms of the concentration parameter Bet01 ; MF00

$`c`$ $`=`$ $`{\displaystyle \frac{9}{1+z}}\left({\displaystyle \frac{M}{1.5\times 10^{13}h^1M_{}}}\right)^{0.13},`$
$`\overline{\delta }`$ $`=`$ $`{\displaystyle \frac{200c^3}{3\left[\mathrm{ln}(1+c)c/(1+c)\right]}},`$ (14)
$`r_s`$ $`=`$ $`{\displaystyle \frac{1.63\times 10^5}{\mathrm{\Omega }_m^{1/3}c}}\left({\displaystyle \frac{M}{h^1M_{}}}\right)^{1/3}h^1\mathrm{Mpc}.`$
We evaluate the mass density $`\rho `$ on the right hand side of eq. (12) exactly by adding the CDM density given above and the neutrino density computed from previous time steps. We did find that approximating $`\rho `$ with the CDM profile alone (i.e. ignoring the neutrino contribution to the total potential) changes our results by no more than 10% for the light neutrino masses and cosmological parameters considered in this paper. We also tested the simplifying assumption made in BKT87 , which allowed them to reduce the integrals in eq. (12) analytically to a single integral by using $`(e^{q/T_{\nu ,0}}+1)^1=Ae^{q/T_{\nu ,0}}`$, i.e., by assuming a Maxwell-Boltzmann instead of Fermi-Dirac distribution. We found this simplification to cause a $``$ 20 % difference at small scales, so we do not use this approximation here.

Before presenting our results for the realistic models above, we first conduct a comparison study by checking our results for the neutrino density profile against those in the numerical simulations of Ref. KKPH96 , which investigated the clustering of CDM and neutrinos in two flat cosmological models: $`\mathrm{\Omega }_{\mathrm{cdm}}=0.8`$ and $`\mathrm{\Omega }_\nu =0.2`$ with two species of 2.3 eV neutrinos, and $`\mathrm{\Omega }_{\mathrm{cdm}}=0.7`$ and $`\mathrm{\Omega }_\nu =0.3`$ with one species of 7 eV neutrinos. Both models assumed $`h=0.5`$. These models are no longer consistent with current observations, but the simulation results provide a useful tool for us to test the validity of our Boltzmann approach. For a fair comparison, we use the CDM halo profile found in KKPH96  as an input:

$$\rho _{\mathrm{cdm}}(r)=\frac{C}{r(r+R)^\alpha },$$
(15)
where $`R`$ and $`\alpha `$ are 0.3 Mpc and 1.5 for the $`\mathrm{\Omega }_\nu =0.2`$ model and 0.11 Mpc and 1.1 for the $`\mathrm{\Omega }_\nu =0.3`$ model. We note that their outer profile is shallower than eq. (13) from the pure CDM simulations NFW96 . In Ref. KKPH96  this difference was attributed to the differing spectral index $`n`$ for the matter power spectrum of the neutrino models at the mass scale of the simulated halos ($`1.3\times 10^{15}`$ M): $`n1.36`$ for $`\mathrm{\Omega }_\nu `$=0.2 and $`n1.53`$ for $`\mathrm{\Omega }_\nu `$=0.3. This argument, however, appears inconsistent with the near universal nature of the halo density profile reported in Ref. NFW96 . A better understanding for the origin of halo profiles should help resolve this issue.

Fig. 1 compares the ratio of the neutrino and CDM density profiles from our approach vs. the two simulated halos in Ref. KKPH96 . For both halos we have used the same cosmological parameters, CDM profiles, and halo mass in our Boltzmann approach as in their simulations. We find a good agreement between the two methods for the outer parts of the cluster; whereas our results are lower by about 50% in the inner parts. This discrepancy may be due to our neglecting the fourth term in eq. (7), or due to statistical fluctuations in the substructure in their simulations. (Only two simulated halos are presented in Ref. KKPH96 .) The Boltzmann approach used here also allows us to explore how neutrinos respond to different CDM potentials. As an illustration of this, we show in Fig. 1 how the neutrino profile changes when the input CDM profile is changed from eq. (15) to the higher resolution profile of eq. (13). We conclude from Fig. 1 that we can obtain a reasonable estimate for neutrino clustering using eq. (12) instead of full scale $`N`$-body simulations.

We now turn to our results for the realistic cosmological models and neutrino masses given at the beginning of this section. The four panels in Fig. 2 show the neutrino overdensity $`n_\nu /\overline{n}_\nu `$ computed from eq. (12) for four models of neutrino masses. The more massive neutrinos clearly cluster more because of their lower thermal velocities. Within each panel, the four curves illustrate how $`n_\nu `$ increases with halo masses from $`10^{12}`$ to $`10^{15}M_{}`$ as a result of the deeper halo gravitational potentials. The growth of $`n_\nu `$ in the inner parts of the halo, where it is almost independent of $`r`$, is illustrated in Fig. 3 for 0.7 and 0.4 eV neutrinos in $`10^{15},10^{14}`$, and $`10^{13}M_{}`$ halos. Most of the clustering is seen to occur at low redshifts.

Unlike the CDM density that continues to rise towards the inner part of a halo as $`\rho 1/r`$, all curves in Fig. 2 show that $`n_\nu `$ flattens out at some radius. Similar features were also seen in Ref. BW88  for neutrinos clustered around cosmic strings. This relative suppression in the neutrino vs. CDM overdensity reflects neutrino free streaming, which dampens and retards perturbation growth on small length scales due to phase mixing. The neutrino damping scale, $`R_d`$, can be characterized by the length scale above which neutrinos behave like the CDM. The standard method to solve the Boltzmann equation for a fluid with pressure invloves the transformation of the Boltzmann equation into an infinite hierarchy of velocity moment equations MB95 , where the lowest three moments with $`l=0,1`$, and 2 correspond to the density, velocity, and shear of the fluid. The choice of the truncation of the hierarchy depends on the physical properties of the fluid and the length scales. For CDM, for example, all modes above $`l=1`$ are zero. The parameter $`R_d`$ gives the scale above which the Boltzmann hierarchy for neutrinos can be truncated at $`l=1`$ (as for the CDM), and below which more $`l`$-modes must be included to compute the neutrino damping effect accurately. This scale is given by DGS96

$$R_d(\tau )\frac{\tau }{\sqrt{1+[a(\tau )/a_{nr}]^2}},$$
(16)
where $`a_{nr}3T_{\nu ,0}/m_\nu `$ is the expansion factor at which the neutrinos become non-relativistic. For the cosmological models considered in this paper, the scale is $`R_d(z=0)5.8/m_\nu `$(eV) Mpc. From Fig. 1, we indeed find the ratio $`\rho _\nu /\rho _{cdm}`$ to be about the cosmic mean value $`\mathrm{\Omega }_\nu /\mathrm{\Omega }_{cdm}`$ at scales above $`R_d`$ and to decrease gradually at smaller radii (top panels), with a final flattening in the neutrino overdensity at $`r0.1R_d`$ (bottom panels). Fig. 2 shows that the radius at which $`\delta _\nu `$ flattens out depends weakly on the mass of the CDM halo. It occurs at smaller radii for less massive CDM halos primarily because the lower mass halos provide shallower potential wells. The damping scale $`R_d`$ of eq. (16) is to be contrasted with the neutrino free streaming distance, which is typically defined as the comoving distance traversed from the time of neutrino decoupling to $`a_{nr}`$: $`\lambda _{fs}_{\tau _i}^{\tau _{nr}}𝑑\tau ^{}/a(\tau ^{})600/m_\nu `$(eV) Mpc KT90 . The distance $`\lambda _{fs}`$ reflects the global streaming motion of neutrinos but not the local clustering properties of neutrinos around CDM after they become non-relativistic.

Another way to understand the results in Figures 1 and 2 is to compare the thermal velocities of neutrinos with the velocity dispersions of the CDM halos: neutrinos can cluster significantly only if their mean thermal velocity in eq. (1) is below the typical velocity of the host CDM halo. Fig. 4 compares these two characteristic velocities for a range of neutrino masses and halo masses. Since the NFW profile specifies only the spatial and not the velocity distribution of the CDM halo particles, two velocity ellipsoids are shown for comparison: isotropic, which is appropriate near the center of the halo, and the more radial distribution $`\beta =1v_t^2/v_r^2=0.5`$, which is appropriate for the outer regions. Fig. 4 illustrates that $`<0.15`$ eV neutrinos are too hot to be captured significantly by $`10^{14}M_{}`$ halos, while the more massive neutrinos can fall into progressively lower mass halos, a result consistent with that shown in Fig. 2.

Our results for neutrino clustering can be compared with the Tremaine-Gunn bound TG79 , which gives an upper limit on the neutrino density in the core of a halo based on the argument that the maximum coarse grained phase space density can not exceed the maximum initial phase space density due to phase mixing. Their results are not directly applicable to our problem because in their derivation, neutrinos are assumed as the sole constituent of dark matter, and the coarse grained neutrino distribution is assumed to be Maxwell-Boltzmann instead of Fermi-Dirac for computational convenience. More recent work KTB96  has extended the derivation to models including both CDM and neutrinos and obtained $`\rho _\nu |2\mathrm{\Phi }|^{3/2}m_\nu ^4/12\pi ^4`$, where $`\mathrm{\Phi }`$ is the gravitational potential of the system. For the NFW profile, we find

$`{\displaystyle \frac{n_\nu }{\overline{n}_\nu }}(r;m_\nu ,r_s,\overline{\delta })`$ $``$ $`40\left({\displaystyle \frac{m_\nu }{\mathrm{eV}}}\right)^3\left({\displaystyle \frac{r_s}{\mathrm{Mpc}}}\right)^3\overline{\delta }^{3/2}`$ (17)
$`\times \left[{\displaystyle \frac{r_s}{r}}\mathrm{ln}\left({\displaystyle \frac{r}{r_s}}+1\right)\right]^{3/2},`$
where $`r_s`$ and $`\overline{\delta }`$ are the CDM halo parameters given by eq. (14). For 1.8 eV neutrinos, for example, this formula gives $`n_\nu /\overline{n}_\nu <3.9\times 10^4,3.2\times 10^5,2.7\times 10^6,`$ and $`2.3\times 10^7`$ for $`10^{12},10^{13},10^{14}`$ and $`10^{15}M_{}`$ halos, respectively, at the scale radius $`r_s`$. One can see that this constraint is satisfied by at least three orders of magnitude for all neutrino overdensities in Figures 1 and 2.

## IV Implications for Ultra High Energy Cosmic Rays

In this section we apply the neutrino clustering results from Sec III to the Z-burst model for UHECRs. No previous work on the Z-burst model has included realistic calculations for $`n_\nu `$. Instead, the value of $`n_\nu `$ has been chosen based on certain observational constraints FKR02  or physical arguments W99 ; YSL98  and has differed greatly from $`n_\nu (1`$ \- $`10^5)\overline{n}_\nu `$. For instance, in FKR02 , it is inferred from the CDM distribution in our local universe, but the large smoothing scale $``$ 20 Mpc assumed in the calculation results in $`n_\nu \overline{n}_\nu `$. In contrast, our results from Sec III show that $`n_\nu `$ can be $`\overline{n}_\nu `$ in the inner $`1`$ Mpc of CDM halos. In Refs. W99 ; YSL98 , $`n_\nu `$ is approximated based on phase-space arguments similar to that of TG79 . While this approach gives an upper bound on the neutrino clustering, the actual overdensity can be significantly less, as we have discussed in the previous paragraph. In addition, the neutrino clustering scale of $``$ 5 Mpc assumed in YSL98  is much larger than what we find in our calculations. Our method gives specific predictions for $`n_\nu `$ as a function of halo radius, halo mass, and neutrino mass.

To estimate the cosmic ray flux we follow the standard assumption in the Z-burst model that the UHECRs above the GZK cutoff are produced by the resonant $`\nu \overline{\nu }`$ scattering, while the lower energy events are explained by protons originating from a uniform distribution of extragalactic sources. The latter appears consistent with the isotropic distribution of $`E<4\times 10^{19}`$ eV events detected in AGASA and HiRes BM99 . We compute the cosmic ray flux (in (eV m2 s sr)-1) from the Z-burst model with FKR02

$`F_Z(E)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}dE_p{\displaystyle _0^{R_{max}}}dr\left[{\displaystyle _0^{\mathrm{}}}dE_{\nu _i}F_{\nu _i}(E_{\nu _i},r)n_{\overline{\nu }_i}(r)+{\displaystyle _0^{\mathrm{}}}dE_{\overline{\nu }_i}F_{\overline{\nu }_i}(E_{\overline{\nu }_i},r)n_{\nu _i}(r)\right]`$ (18)
$`\times \sigma _{\nu \overline{\nu }}(s)\mathrm{Br}(Z\mathrm{hadrons}){\displaystyle \frac{\mathrm{d}N_{p+n}}{\mathrm{d}E_p}}\left|{\displaystyle \frac{P_p(r,E_p;E)}{E}}\right|.`$
Here $`F_{\nu _i}(E_{\nu _i},r)`$ is the flux of ultra high energy neutrinos with energy $`E_{\nu _i}`$ at distance $`r`$ and $`n_{\nu _i}(r)`$ is the physical number density of the relic neutrinos. (The repeated index $`i`$ is summed over different neutrino species.) The particle interactions are described by the cross section $`\sigma _{\nu \overline{\nu }}(s)`$ for the Z-boson production at the center-of-mass energy $`s=2m_\nu E_\nu `$, and by the branching ratio $`\mathrm{Br}(Z\mathrm{hadrons})=69.89\pm 0.07`$ % for the subsequent cascade of the Z-boson into hadrons PDGC00 . The factor $`\mathrm{d}N_{p+n}/\mathrm{d}E_p`$ gives the energy distribution of the produced protons and neutrons. The subsequent proton propagation is specified by $`P_p(r,E_p;E)`$, which gives the probability that a proton created at distance $`r`$ with energy $`E_p`$ arrives at Earth with an energy greater than $`E`$. It measures the amount of proton energy degradation due to the resonant photoproduction of pions and other processes discussed in Sec I. Specific values of $`P_p`$ has been calculated in FK01  for the range of parameters considered in this paper. We do not include in our UHECR flux estimate the contributions from the photons produced in the Z-decay because experimental data suggest that less than $`50\%`$ of the cosmic rays above $`4\times 10^{19}`$ eV are photons at the $`95\%`$ confidence level AHVWZ02 ; AG2 . Typical physical mechanisms used to explain the suppressed photon contributions are large universal radio background and sufficiently strong extragalactic magnetic fields ($`10^9`$ Gauss) FKR02 ; PB96 . The study of the effects of these parameters on the UHECR flux due to the Z-burst model can be a subject of future work.

A key feature of the Z-burst model is that the cross section $`\sigma _{\nu \overline{\nu }}(s)`$ for $`\nu \overline{\nu }Z^0`$ is enhanced by several orders of magnitude near the resonant energy in the rest frame of the relic neutrinos W82

$$E_{\nu _i}^{res}=\frac{M_Z^2}{2m_{\nu _i}}=4.2\times 10^{21}\mathrm{eV}\left(\frac{1\mathrm{eV}}{m_{\nu _i}}\right),$$
(19)
where $`M_Z`$ is the mass of the Z boson. The flux in eq. (18) to a good approximation therefore depends only on the neutrino resonant energy and not on the slope of the incident high energy neutrino spectrum. Eq. (18) can then be written as FKR02

$`F_Z(E)`$ $`=`$ $`\overline{\sigma }_{\nu \overline{\nu }}F_{\nu _i}(E_{\nu _i}^{res}){\displaystyle _0^{\mathrm{}}}𝑑E_p{\displaystyle _0^{R_{max}}}𝑑r`$ (20)
$`\times n_{\nu _i}(r)Q_p\left(y={\displaystyle \frac{4m_\nu E_p}{M_Z^2}}\right)`$
$`\times \left|{\displaystyle \frac{P_p(r,E_p;E)}{E}}\right|,`$
where $`n_{\nu _i}(r)`$ is the physical number density of neutrinos and antineutrinos at the Z-burst site at radial distance $`r`$, $`\overline{\sigma }_{\nu \overline{\nu }}=40.4`$ nb is the cross section for $`\nu \overline{\nu }Z^0`$ averaged over the width of the resonance, and $`F_{\nu _i}(E_{\nu _i}^{res})`$ is the incident flux of ultra high energy neutrinos at the resonant energy. The function $`Q_p`$ is the boosted momentum distribution from hadronic Z decays and can been calculated from experimental data FKR02 . It has a fairly broad peak at $`y10^2`$ and falls off approximately as $`y^7`$ for $`y0.5`$. The neutrino flux $`F_{\nu _i}(E_{\nu _i}^{res})`$ remains a free parameter in the Z-burst calculation since no successful astrophysical model yet exists to explain the production of $``$ 1021 eV neutrinos VMO97 ; WB ; GV02 . We do not attempt to model the effect of source evolution in our calculations since it is again an unknown quantity and is easy to incorporate once its nature is known.

We first present results for the cosmic ray spectrum $`F(E)`$ in the Z-burst model ignoring the spatial clustering in the neutrinos, i.e., we assume $`n_\nu =\overline{n}_\nu `$ in eq. (20). This assumption underestimates the flux in the Z-burst model, but we include the results here for comparison since this is a common assumption made in several Z-burst calculations FKR02 ; W99 ; GV02 . Fig. 5 shows the sensitivity of $`E^3F(E)`$ on the neutrino masses. The flux is higher at high energies for smaller $`m_\nu `$ because the momentum distribution $`Q_p`$ of the decay particles peaks at a higher energy for smaller $`m_\nu `$. For a given $`m_\nu `$, $`E^3F(E)`$ decreases rapidly at $`E>10^{21}`$ eV because $`Q_p`$ falls off as $`y^7`$ for $`y0.5`$. The integration is carried out to a maximum distance of $`R_{max}=2000`$ Mpc, but our results are insensitive to this choice as long as $`R_{max}`$ is sufficiently beyond the GZK zone of $`50`$ Mpc.

For ease of comparison, the curves in Fig. 5 are all normalized to the same incident neutrino flux of $`F_{\nu _i}(E_{\nu _i}^{res})=1.7\times 10^{35}`$ (eV m2 s sr)-1 for each of the three neutrino flavors. (For the one flavor $`0.07`$ eV model, the assumed flux is 3 times higher.) We do not attempt to determine this value by performing statistical fits to data because the UHECR spectrum from AGASA (square symbols) and HiRes (triangle symbols) disagree in both amplitude and shape. We do note that for models that have three degenerate neutrino masses of $`m_{\nu _i}<1`$ eV, this value for the neutrino flux is consistent with the existing upper bound from the Goldstone Lunar Ultra-high energy neutrino Experiment (GLUE) GLUE . The 0.07 eV model shown in Fig. 5, however, would need to be lowered by a factor of $``$ 4 in order to satisfy the GLUE upper limit. A better understanding of systematic effects in the GLUE experiment is needed before their results can be used to rule out models.

For comparison, the dotted curve in Fig. 5 shows the cosmic ray flux for protons originating from a uniform distribution of extragalactic sources with a constant comoving density. It is computed from

$`F_{EG}(E)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑E_p{\displaystyle _0^{R_{max}}}{\displaystyle \frac{dr}{R_{max}}}[1+z(r)]^3`$ (21)
$`\times F_p(E_p)\left|{\displaystyle \frac{P_p(r,E_p;E)}{E}}\right|,`$
where the unknown proton injection energy spectrum $`F_p(E_p)`$ is typically assumed to be a power law: $`F_p(E_p)=ϵ^1AE_p^\beta `$. The shape of $`F_{EG}`$ depends on the injection spectrum $`F_p`$, but for definiteness, we have assumed $`\beta =2.4`$ and $`A=5.98\times 10^{31}`$ (and an upper energy cutoff of $`E_p=10^{23}`$), which are found to be the best fit values FKR02  to the existing cosmic ray data that have a total experimental exposure of $`ϵ8\times 10^{16}`$ m2 s sr. The GZK cutoff is clearly seen at $``$ 4 $`\times `$ 1019 eV in the dotted curve. The flux rises beyond $`4\times 10^{20}`$ eV because the photoproduction of pions is a resonant process where the cross section peaks at $`E_p2.3\times 10^{20}`$ eV GZK  and decreases at higher energies, allowing a larger fraction of protons to reach us.

The predictions for the UHECR spectrum change significantly when we incorporate the neutrino overdensity computed in Sec III. To make realistic estimates for our local universe, we consider five lines of sight towards five of the most massive nearby clusters: Virgo, Centaurus, Hydra, Perseus-Pisces, and Coma, where the highest overdensity of neutrinos are expected. The distance, mass, rough angular extent, and equatorial coordinates of each of the clusters are listed in Table 1 VIZ ; NED . The cluster masses are taken from http://cfa-www.harvard.edu/$`\stackrel{~}{}`$huchra/clusters, where they are estimated from galaxy velocities and the virial theorem. We caution that these values have large error bars. The mass of the nearest cluster Virgo Bet94 ; TS ; SBB99 ; FSSB01 , for example, has been estimated to be $`1.56\times 10^{14}`$ M based on X-ray emission measurements by ROSAT Bet94 , to $`1.5\times 10^{15}`$ M based on the relativistic Tolman-Bondi method FSSB01 . (The Tolman-Bondi model is based on analytic solutions to the Einstein field equations for spherically symmetric pressure free overdensities in a homogeneous universe FSSB01 ; ELTPF99 .)

Fig. 6 shows our predictions for the cosmic ray flux (in (eV m2 s sr)-1) towards these five lines of sight for four different neutrino masses. Along each line of sight, high energy neutrinos from extragalactic sources are assumed to traverse a uniform sea of background neutrinos plus an overdensity of background neutrinos centered at the location of the given cluster, where $`n_\nu `$ is computed from eq. (12) for the mass of the cluster. We also include in our calculation a local $`n_\nu `$ for the Local Group of mass $`4\times 10^{12}M_{}`$ SS98 . Despite the smaller mass, the proximity of the Local Group to us leads to non-negligible contributions to the UHECR flux: about a factor $``$ 2 for the $`3\times 0.6`$ eV model, and up to a factor of $`10`$ at $`E10^{20}`$ eV for the 1.8 eV model. The difference is primarily due to the more efficient clustering of 1.8 eV neutrinos compared to 0.6 eV neutrinos in the Local Group.

Our main conclusion from Fig. 6 is that the flux of UHECRs in the Z-burst model should show significant anisotropy if $`m_\nu >0.3`$ eV, with the largest flux coming from the Virgo cluster. For $`m_\nu <0.1`$ eV, on the other hand, neutrinos are too hot to cluster appreciably even around the largest clusters in the universe, and the UHECR flux in the Z-burst model is nearly isotropic. We choose to plot in Fig. 6 the flux per steradian because the angular extent of the clusters cannot be precisely defined, but one can easily use the approximate angular extents of the clusters listed in Table 1 to estimate the expected anisotropy in the signal.

## V Conclusion and Discussion

We have introduced and tested a method based on the collisionless Boltzmann equation to calculate the gravitational clustering of massive neutrinos in CDM halos for realistic cosmological models. This method is valid for currently favored models with $`\mathrm{\Omega }_{\mathrm{cdm}}\mathrm{\Omega }_\nu `$ in which the clustering of neutrinos is mostly determined by the existing CDM halos while the clustering of the CDM is little affected by the neutrinos. One can then obtain the neutrino phase space distribution by solving the collisionless Boltzmann equation in a background potential given by the universal profile of CDM halos from high resolution simulations. The resulting Boltzmann equation is linear in the neutrino density contrast and has tractable intergral solutions that require negligible computational time in comparison with N-body simulations. This method has enabled us to obtain specific predictions for the neutrino overdensity as a function of halo radius, halo mass, and neutrino mass for a wide range of parameters.

Our calculation shows that neutrinos with masses $``$ 0.3 eV can cluster appreciably in CDM potential wells with masses $`10^{13}M_{}`$. The predicted neutrino overdensity increases with both the neutrino mass and the halo mass, ranging from $`10`$ for 0.3 eV neutrinos in $`10^{13}M_{}`$ halos to $`1500`$ for 1.8 eV neutrinos in $`10^{15}M_{}`$. Specific predictions are plotted in Figs. 2 and 3. Neutrino clustering has a strong impact on the Z-burst model that has been proposed as a possible explanation for the UHECR events. The predicted UHECR spectrum shown in Figs. 5 and 6 depends sensitively on the neutrino mass and overdensity, showing distinct spectral features towards nearby galaxy clusters if $`m_\nu 0.3`$ eV.

To illustrate the effects of neutrino mass and overdensity on the UHECR spectrum, we have chosen to normalize the flux in Figs. 5 and 6 with the same value (i.e. $`F_{\nu _i}(E_{\nu _i}^{res})=1.7\times 10^{35}`$ (eV m2 s sr)-1 for each flavor for the three degenerate mass models and three times higher for the one massive species model) instead of adjusting it by fitting individual spectrum to existing data. We have nonetheless included current data from the AGASA AG  and HiRes HR  experiments in Figures 5 and 6 for comparison. More events are needed to discriminate the different models and the directional dependence. The large increase in flux towards Virgo is an interesting signature of the Z-burst model for upcoming experiments such as Auger AUG  and OWL OWL99  that will provide an angular resolution of $``$ 1. Experimental limits on the anisotropy would in turn imply small neutrino inhomogeneities in the Z-burst model and can be used to place upper bounds on the neutrino mass.

A useful constraint on the Z-burst model is provided by the Energetic Gamma Ray Experiment Telescope (EGRET) measurement of the GeV $`\gamma `$-ray background flux, which must not be exceeded by the high energy photons produced in the Z-burst models once they cascade down to the GeV energy range. The result depends on the assumed redshift evolution of the sources that produce the incident high energy neutrinos, and on whether the sources themselves produce photons. The normalization of the neutrino flux cited in the previous paragraph rules out sources emitting a comparable flux in $`\gamma `$-rays because it leads to a conflict with the existing EGRET limits for the GeV $`\gamma `$-rays. For pure neutrino sources, calculations based on particle transport codes show that for neutrino masses of 0.1, 0.5, and 1 eV (ignoring neutrino clustering), the EGRET bound is met for $`\alpha <3`$, $`<0`$, and $`<3`$, respectively, where the source number density evolves as $`(1+z)^\alpha `$ FKR02 ; KKSS02 . When neutrino clustering is taken into account, our results from Fig. 2 show that the bound above for $`m_\nu <0.3`$ eV should be unaffected since they do not cluster appreciably in the Local Group. For larger neutrinos masses, however, we expect a less stringent bound on the source evolution due to local neutrino clustering. To derive quantitative constraints would require detailed transport calculations.

The implications of the neutrino clustering results presented in this paper extend beyond the problem of the UHECR spectrum. For UHECR, upcoming experimental results may converge on a spectrum that is consistent with the GZK cutoff and would therefore not require models such as the Z-burst. It is also likely that the Z-burst model is not the correct explanation for the UHECR events. However, the neutrino-anti-neutrino resonance scattering process remains one of few ways to detect the relic neutrinos, as first suggested in Ref. W82 . This paper has addressed neutrino clustering, a major uncertainty in all studies concerning relic neutrinos.

###### Acknowledgements.

We thank Ed Bertschinger, Tom Weiler, Jon Arons, Bhuvnesh Jain, and Nick Sarbu for useful discussions. We also thank Dmitry Semikoz, John Beacom, and the referee for comments that have helped to improve the manuscript. C.-P. M. acknowledges support of an Alfred P. Sloan Foundation Fellowship, a Cottrell Scholars Award from the Research Corporation, and NASA grant NAG 5-12173. A portion of the work was carried out at the Aspen Center for Physics. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration and the VizieR catalogue access tool, CDS, Strasbourg, France.


format: text

text:
# Plasticity of Ferromagnets near the Curie Point

## I Introduction

Various experimental data indicate a strong influence of the magnetic structure of ferromagnets on their plastic properties. Transition to a magnetically ordered state is accompanied by a strong change of plasticity (Zackay and Hazlett, 1953, Nabutovskaya, 1969,1971, Wolfenden, 1978, Bolling and Richman, 1969, Echigoya etal, 1973, Echigoya and Hayashi, 1979, Flor etal, 1980, Retat etal, 1985, Retat, 1987, Maksimova and Maiboroda, 1992, Gulayev and Svistunova, 1996). Measurements with small temperature intervals revealed an increase of plasticity in a close vicinity of the Curie point (Zackay and Hazlett, 1953, Nabutovskaya, 1969,1971, Wolfenden, 1978). This fact by itself is not so astonishing since many properties of ferromagnets change near the Curie point (Vonsovskii, 1974), however, the giant value of the observed effect is really surprising. According to Nabutovskaya (1969,1971) the microhardness of gadolinium changes near $`T_C`$ by a factor of two, whereas its elastic constants vary by a percent or so. For example, the Young modulus of gadolinium decreases in the Curie point only by 1.6% (Spichkin etal, 1999).

Recently we considered a mechanism of the influence of a magnetic field on plasticity of some nonmagnetic crystals (Molotskii and Fleurov, 1997, Molotskii, 2000, and references therein). The magnetic field induces transitions (singlet to triplet) between different spin states of the radical pairs formed by dangling bonds of a dislocation cores and paramagnetic obstacles in the course of the pair formation. The state of the pair may be either bonding or antibonding, depending on its spin configuration. The magnetic field can influence the relative occupations of these states and lead to an increase of probability of the dislocation detachment from paramagnetic obstacles. The crystal plasticity, as a result, increases.

This approach appeared to be very successful and allowed for an explanation of a large number of plasticity related phenomena in nonmagnetic crystals. We propose here to apply the same approach when discussing plasticity of ferromagnets. In this case, the role of the external magnetic field can be played by the local magnetic induction created by the spontaneously magnetised surrounding atoms. Recently, we demonstrated that this internal magnetic field may be important for the Invar hardening, leading to an increase of the critical resolved shear stress of Invar alloys with a lowering temperature (Molotskii and Fleurov, 2001). We expect to observe a temperature variation of various plasticity related characteristics of the crystal connected with the variation of the spontaneous magnetization with the temperature. When approaching the Curie point the role of the magnetization fluctuations may become more pronounced. The typical time of relaxation of the magnetization fluctuations is usually much larger than the characteristic times of formation or rupture of a dislocation - obstacle bond. Therefore we should consider this process under the influence of slowly varying magnetic inductions due to fluctuating magnetization. We will show below that this may result in a strongly increasing plasticity in the vicinity of the Curie point.

It is worth mentioning here that the theory, we develop in this paper, will not take into account an interaction of dislocations with domain walls, which may serve as efficient pinning centres for dislocation (Seeger etal, 1964). This assumption is justified as long as the typical domain size essentially exceeds the dislocation path length. According to the experiments of Zackay and Hazlett (1953), and of Nabutovskaya (1969, 1971) this path is limited by the average distance, $`L_f1/\sqrt{\rho _f}`$, between the forest dislocations. The density of the forest dislocations in those experiments was $`\rho _f10^9`$cm-2, hence the dislocation path length could not exceed $`L_f3\times 10^5`$cm, which was two to three orders of magnitude smaller than the typical domain sizes.

## II Influence of Internal Magnetic Induction on Plasticity of Magnets

Internal magnetic field of a ferromagnet due to spontaneous magnetization is rather large ($`1`$T) and is capable of influencing plastic properties of crystals. Such an important characteristic of crystal plasticity as the dislocation path length increases with the magnetic induction $`B`$ as (Molotskii, 2000)

$$L(B)=L_0(1+\frac{B^2}{B_0^2})$$
(1)
where $`B_0`$ is a constant characterizing the dislocation – obstacle bond. Usually its value lies in the range from 0.2 to 1T. In the case of nonmagnetic crystals $`B`$ is the value of the magnetic induction created by an external source. In the case of a ferromagnet we do not need an external magnetic field since a strong enough internal magnetic induction is always present. This internal magnetic induction decreases when temperature approaches the Curie point $`T_C`$ and starts strongly fluctuating. One may expect that in the critical region, at $`TT_C`$, these strong fluctuations of the internal magnetic induction will strongly influence plastic properties of the ferromagnet.

The magnetic field influences kinetics of the formation of the dislocation – obstacle bonds. The characteristic times of these processes are typically on the order of $`10^7`$s (Molotskii and Fleurov, 1997), which is shorter than the typical times of the large scale magnetization fluctuations determining the local internal magnetic induction $`B_{int}`$ at each bond (see, e.g., Ma, 1976). This allows one to consider the induction $`B_{int}`$ at each moment of time as created by a particular space distribution of the magnetization fluctuations in the ferromagnet.

For example, we estimate here the average dislocation free path length. According to equation (1) it is determined by the mean square of the internal magnetic induction,

$$B_{int}^2=\overline{B}_{int}^2+\mathrm{\Delta }B^2.$$
(2)
Here $`\overline{B}_{int}`$ is average value of the internal magnetic induction, whereas $`\mathrm{\Delta }B=B_{int}\overline{B}_{int}`$.

Fluctuations of the local internal magnetic induction acting on a particular dislocation – obstacle bond can be described using the theory of the second order phase transitions (see, e.g., Landau and Lifshitz, 1980). The local internal magnetic indiction of a homogeneously magnetised crystal can be estimated with the help of the Lorentz formula

$$𝐁=\frac{4\pi }{3}𝐌.$$
(3)
The probability that the magnetization $`M`$ deviates by $`\mathrm{\Delta }M`$ from its average value $`\overline{M}`$ in the volume $`V`$ around the bond is

$$w(\mathrm{\Delta }M,V)=\frac{1}{2}\frac{V_{}^{}{}_{}{}^{1/2}}{(2\pi k_BT\chi )^{3/2}}\{\begin{array}{cc}\mathrm{exp}\left\{\frac{\mathrm{\Delta }𝐌^2V}{2k_BT_C\chi }\right\},& V>V^{}\\ 0,& V<V^{}\end{array}$$
(4)
where $`\chi `$ is the magnetic susceptibility, $`k_B`$ is Boltzmann constant. The probability density function (4) is truncated at volumes smaller than certain characteristic volume $`V^{}`$ which will later serve as a fitting parameter. The general approach to the order parameter fluctuations, applied here, considers the ferromagnet as continuus medium. Hence, it holds at distances which are larger than typical interatomic distances. Therefore one may expect the typical scale of the excluded volume $`V^{}`$ not to exceed essentially several lattice spacings. Similar approach is implied by equation (3), which is obtained by considering a small empty volume surrounded by a magnetised medium.

The mean square fluctuations of the magnetization can be then calculated as

$$\mathrm{\Delta }𝐌^2=_V^{}^{\mathrm{}}d^3\mathrm{\Delta }M(\mathrm{\Delta }𝐌)^2w(\mathrm{\Delta }M,V)=\frac{k_BT_C\chi }{V^{}}$$
(5)
Now using equations (2) and (5) the mean square local magnetic field acting on the dislocation – obstacle bond can be obtained using the following equations

$$𝐁_{int}{}_{}{}^{2}=\overline{B}_{int}{}_{}{}^{2}+\mathrm{\Delta }B^2=\frac{16\pi ^2}{9}\overline{𝐌}{}_{}{}^{2}+\frac{16\pi ^2}{9}\frac{k_BT_C\chi }{V^{}}.$$
(6)
An inhomogeneous spatial distribution of the magnetization, caused by fluctuations, results in a demagnetization field, which should have been accounted for in equation (3). The demagnetization field is the strongest for the spherical fluctuations and is negligible for a fluctuation having the form of a plate. The exact account of the demagnetization is a very tedious task but finally leads to a numerical factor, smaller than one, in the second term in equation (6). This factor is absorbed in the volume $`V^{}`$, which anyhow serves as a fitting parameter.

The temperature dependence of the mean square local magnetic field in the vicinity of the Curie point can be readily calculated. For this we may use the temperature dependence of the magnetic susceptibility which has the form

$$\chi (T)=\{\begin{array}{cc}\frac{C}{2T_C}\left(\frac{T_C}{T_CT}\right)^\gamma ,& \text{at}T<T_C\\ & \\ \frac{C}{T_C}\left(\frac{T_C}{TT_C}\right)^\gamma & \text{at}T>T_C.\end{array}$$
(7)
Here $`\gamma `$ is a critical index,

$$C=\frac{np^2\mu _B^2}{3k_B}$$
(8)
is the Curie constant, in which $`n`$ is the particle density, $`p`$ is the effective number of the Bohr magnetons $`\mu _B`$ per atom.

The first term in equation (6) contains the average spontaneous magnetization $`\overline{M}(T)`$ whose temperature dependence is

$$\overline{M}(T)=\{\begin{array}{cc}M_0\left(1\frac{T}{T_C}\right)^\beta & \text{at}T<T_C,\\ & \\ 0& \text{at}T>T_C.\end{array}$$
(9)
Here $`M_0`$ is the limiting value of the magnetization at low temperature, $`\beta `$ is the corresponding critical index

## III Temperature dependence of the critical stress: nickel as a test case

An influence of the magnetization on plasticity was first discovered half a century ago by Zackay and Hazlett (1953), who studied the temperature dependence of the critical stress $`\sigma _c`$ for nickel. They have found that the dependence of $`\sigma _c(T)`$ passes a minimum near the Curie point $`T_C`$. Samples with 99.95% nickel content were used. There was also a small amount of iron (0.03%) and magnesium (0.02%) atoms. The paramagnetic iron atoms are efficient obstacles for dislocations in nickel, and they are most probably responsible for the observed sensitivity of the plasticity of nickel to the magnetic field.

The temperature dependence of the critical stress in a ferromagnet can be found, accounting for the part played by the internal magnetic fields. Dislocations can be bound to paramagnetic obstacles with the binding energy $`W_M`$. Then the critical stress is (Friedel, 1964)

$$\sigma _c(0)=\frac{|W_M|}{b^2l(B)}\left[1\left(\frac{T}{T_0}\right)^{\frac{2}{3}}\right]^{\frac{3}{2}}$$
(10)
where $`b`$ is the value of the dislocation Burgers vector, and $`l`$ is the average distance between the obstacles. The temperature dependent factor in equation (10) accounts for the thermal activation processes (Haasen, 1983). The parameter $`T_0`$ is proportional to the dislocation - obstacle binding energy.

It was shown in our paper (Molotskii and Fleurov, 1997) that the average length of the dislocation free segment depends on the magnetic field similarly to (1),

$$l(B)=l_0\left(1+\frac{B^2}{B_0^2}\right).$$
(11)
Considering a ferromagnet we should introduce the fluctuating internal local magnetic induction $`𝐁_{int}=\overline{𝐁}_{int}+\mathrm{\Delta }𝐁`$ which acts on each particular dislocation – obstacle bond. It means that we have to average equation (11) over various values and orientations of the vector $`\mathrm{\Delta }𝐁`$. This averaging for the case of an easy axis ferromagnet is carried out in Appendix), so one gets

$$\frac{\sigma _c(T)}{\sigma _{c0}}=\overline{f}(\overline{B},\mathrm{\Delta }B^2)\left[1\left(\frac{T}{T_0}\right)^{\frac{2}{3}}\right]^{\frac{3}{2}}.$$
(12)
The critical stress appears now to be a function of the average spontaneous magnetization and of its fluctuations. The temperature dependence of both these quantities is determined by equations (6), (7), and (9). As a result we now have an equation for the temperature dependence of the critical stress in the vicinity of the Curie temperature.

Figure 1 presents a comparison of the temperature dependence of the critical stress in Ni calculated by means of equation (12) with the available experimental data.

The theoretical curve is plotted using the following parameters for Ni: $`T_C=627`$K, $`M_0=510`$G, $`p=0.606`$ (Kittel, 1986), $`\beta =0.33`$, $`\gamma =1.33`$ (Kadanoff etal, 1967). The effective number, p, of Bohr magnetons in Ni is rather small and, hence, one gets a small Curie constant $`C=7\times 10^3K`$. The remaining parameters, $`T_0=2070`$K, $`B_0=0.49`$T, $`V^{}=2.64\times 10^{22}`$cm3, are determined by fitting equation (12) to the experimental data by Zackay and Hazlett (1953) for their samples with 2% deformation. One observes a reasonable general agreement between the theory and experiment. However, the measurements were carried out with a $``$ 20K interval between the points. As a result the sharp minimum in the temperature dependence of the critical stress near the Curie point, predicted by our theory, might have been overlooked. A relatively small, $``$ 15%, decrease of the critical stress in this region is indicative of a possible existence of a sharper minimum in this temperature range.

## IV Temperature Dependence of the Microhardness: Gadolinium as a test case

This section discusses the behavior of crystal microhardness near the Curie point. The results are compared with the available experimental data on gadolinium. Currently there is no microscopic theory of microhardness, and hence, no microscopic theory of the influence of a magnetic field on it. Nevertheless, the dependence of the microhardness on a magnetic field can be estimated by means of the following simple considerations. It is known that the microhardness $`H`$ varies inversely proportionally to the plasticity — the higher the plasticity, the lower the hardness. On the other hand, the plasticity is directly connected to the dislocation path length. The magnetic field dependence of $`L(B)`$ is given by equation (1). Hence, we assume that the average microhardness of a crystal as

$$H=H_0\frac{L_0}{L(B)}$$
(13)
where $`H_0`$ is the microhardness in the absence of a magnetic field.

Assuming an isotropic distribution of the orientations one gets

$$H=H_0\overline{f}(\overline{B},\mathrm{\Delta }B^2)$$
(14)
where the function $`\overline{f}(\overline{B},\mathrm{\Delta }B^2)`$ is calculated in Appendix. Similarly to the previous section, we use equations (6), (7), and (9) in order to get equation (14) for the temperature dependence of the microhardness in the vicinity of the Curie temperature.

A comparison of this theoretical result with the experimental data available for the microhardness of gadolinium is presented in Figure 2.

The Curie constant $`C`$ for gadolinium can be found using equation (8) and the values of the relevant parameters: $`n=3.02\times 10^{22}`$cm-3, $`p=7.10`$, $`M_0=2010`$G (Kittel, 1986), which lead to $`C=0.33`$K. The Curie temperature for gadolinium is $`T_C=292`$K (Nabutovskaya, 1969). The critical indices are $`\beta =0.3265`$ and $`\gamma =1.239`$ (Aliev etal, 1988).

The value of the parameter $`B_0`$ could have been possible to determine from the magnetic field dependence of the plasticity. Unfortunately, we do not know about such measurements in gadolinium. However, magneto- (Al’shits etal, 1990) and electroplastic (Okazaki etal, 1979) effects have been measured in zinc and titanium, who, similarly to gadolinium, possess hcp structures. Fitting these results to the theory (Molotskii, 2000) provides $`B_0=0.70`$T for Zn and $`B_0=0.94`$T for Ti.

The theoretical curve in Figure 2 uses the values $`B_0`$ and $`V^{}`$ as fitting parameters. The obtained parameter $`B_0=0.87`$T for gadolinium appears to be closed to its experimentally measured values for Zn and Ti. The characteristic volume is chosen $`V_0=1.85\times 10^{20}`$cm3. It corresponds to a sphere with the radius 17Å, about four interatomic spacings. The results are not very sensitive to the variation of this volume, which can be chosen smaller (up to $`2\times 10^{21}`$cm3) without damaging the overall agreement.

## V Conclusion

Accounting for the internal magnetic induction of ferromagnets and for its fluctuations leads us to the conclusion that various plastic properties of crystals should be singular near the Curie point. The critical behavior of the magnetic susceptibility is transferred to the temperature dependence of the critical stress and microhardness. It would be interesting to measure other plasticity characteristics, such as dislocation path length $`L`$ or plastic strain rate $`\dot{\epsilon }`$ near the Curie point in order to look for a possible singular behavior, which follows from the theory proposed in this paper. For example, one may expect that carrying out measurements similar to those of Zackay and Hazlett, 1953, on the temperature dependence of the critical stress in Ni, but with a smaller temperature interval a sharp minimum at $`T=627`$K may be observed.

An external magnetic field is known to smear out the phase transition, to decrease the magnetic susceptibility, and to suppress fluctuations near $`T_C`$. As shown by Dan’kov etal (1998) a 0.5T magnetic field diminishes the magnetic susceptibility near the Curie point in gadolinium by an order of magnitude. This should lead to a nearly complete suppression of the magnetization fluctuations and, hence, to a suppression of the singular temperature dependence of microhardness. As for nickel, a stronger field about 3T (Hischler and Rocker, 1966) is necessary to suppress the fluctuations, hence correspondingly higher field will be necessary to influence the singular behavior of the critical stress in nickel near the Curie point.

We may conclude that an experimental observation of a suppression of the above singularities in the temperature dependencies of plasticity characteristics in a magnetic field will be an important experimental evidence in favour of the theoretical approach to the plasticity of ferromagnets presented in this paper.

Acknowledgment

We are grateful to D. Fleurov for the help with the computation.

## Appendix A Averaging for an easy axis cubic crystal

The averaging is carried out assuming a cubic symmetry. It is certainly relevant to Ni. As for the hexagonal Gd, its anisotropy is rather weak and the result obtained here, will be also applied for this crystal. We have to average the quantity

$$f(B)=\frac{1}{1+{\displaystyle \frac{𝐁^2}{B_0^2}}},$$
(15)
where $`𝐁=\overline{𝐁}+\mathrm{\Delta }𝐁`$, over the distribution (4). Disregarding a possible anisotropy of the fluctuations of the magnetic induction in the vicinity of the dislocation – obstacle bonds this averaging can be presented as

$$f(B)\overline{f}(\overline{B},\mathrm{\Delta }B^2)=\frac{1}{(2\pi )^{1/2}V^{}\mathrm{\Delta }B^2^{3/2}}_V^{}^{\mathrm{}}dV_0^{\mathrm{}}\mathrm{\Delta }B^2d\mathrm{\Delta }B_1^1\frac{d\mathrm{cos}\theta }{2}\times $$

$$\frac{1}{1+{\displaystyle \frac{1}{B_0^2}}(\overline{B}{}_{}{}^{2}+\mathrm{\Delta }B^2+2\overline{B}\mathrm{\Delta }B\mathrm{cos}\theta )}\mathrm{exp}\left\{\frac{\mathrm{\Delta }B^2V}{2\mathrm{\Delta }B^2V^{}}\right\}$$
(16)
Now integration over $`\mathrm{\Delta }B`$ is carried out by parts and the averaging (16) is represented as

$$\overline{f}(\overline{B},\mathrm{\Delta }B^2)=$$

$$\frac{1}{2(2\pi )^{1/2}\mathrm{\Delta }B^2^{1/2}}_V^{}^{\mathrm{}}𝑑V\frac{1}{V}_{\mathrm{}}^{\mathrm{}}𝑑\mathrm{\Delta }B\frac{1}{1+{\displaystyle \frac{1}{B_0^2}}(\overline{B}+\mathrm{\Delta }B)^2}\mathrm{exp}\left\{\frac{\mathrm{\Delta }B^2V}{2\mathrm{\Delta }B^2V^{}}\right\}+$$

$$\frac{1}{2(2\pi )^{1/2}\mathrm{\Delta }B^2^{1/2}\overline{B}}_V^{}^{\mathrm{}}\frac{dV}{V}_{\mathrm{}}^{\mathrm{}}𝑑\mathrm{\Delta }B\frac{\mathrm{\Delta }B}{1+{\displaystyle \frac{1}{B_0^2}}(\overline{B}+\mathrm{\Delta }B)^2}\mathrm{exp}\left\{\frac{\mathrm{\Delta }B^2V}{2\mathrm{\Delta }B^2V^{}}\right\}$$
(17)
or

$$\overline{f}(\overline{B},\mathrm{\Delta }B^2)=\frac{1}{2\sqrt{2\pi }}_1^{\mathrm{}}\frac{dv}{v}_{\mathrm{}}^{\mathrm{}}𝑑x\frac{a^2}{a^2+(b+x)^2}(1+\frac{x}{b})\mathrm{exp}\left\{\frac{x^2v}{2}\right\}$$
(18)
where the notations $`a=B_0/\sqrt{\mathrm{\Delta }B^2}`$, $`b=\overline{B}/\sqrt{\mathrm{\Delta }B^2}`$, $`v=V/V^{}`$, and $`x=\mathrm{\Delta }B/\sqrt{\mathrm{\Delta }B^2}`$ are introduced. The identity

$$\frac{1}{a\pm i(b+x)}=_0^{\mathrm{}}𝑑p\mathrm{exp}\{p[a\pm i(b+x)]\}$$
(19)
is substituted into (18). Now integration over $`x`$ is carried out and a new integration variable $`z=p/\sqrt{v}`$ is introduced instead of $`v`$. After that, $`p`$ is substituted for $`pa`$ and one gets

$$\overline{f}(\overline{B},\mathrm{\Delta }B^2)=\frac{B_0}{\sqrt{\mathrm{\Delta }B^2}}\underset{0}{\overset{\mathrm{}}{}}dp\times $$

$$\left\{\sqrt{\frac{\pi }{2p}}e^p\mathrm{\Phi }\left(\frac{\sqrt{\mathrm{\Delta }B^2}p}{B_0\sqrt{2}}\right)\left[\mathrm{cos}\left(\frac{\overline{B}p}{B_0}\right)\frac{B_0}{4\overline{B}p}\mathrm{sin}\left(\frac{\overline{B}p}{B_0}\right)\right]+\frac{B_0}{4\overline{B}p}e^{p{\displaystyle \frac{\mathrm{\Delta }B^2p^2}{2B_0^2}}}\mathrm{sin}\left(\frac{\overline{B}p}{B_0}\right)\right\}$$
(20)
where

$$\mathrm{\Phi }(x)=\frac{2}{\sqrt{\pi }}_0^xe^{t^2}𝑑t$$
(21)
is the probability integral.

The value of $`\overline{f}(\overline{B},\mathrm{\Delta }B^2)`$ can be analytically estimated in the limit of small fluctuations, when $`\sqrt{\mathrm{\Delta }B^2}B_0`$. The integral (20) converges for small values of $`p`$ and the expansion $`\mathrm{\Phi }(x)2x/\sqrt{\pi }`$ can be used. It converges to a rather obvious result

$$\overline{f}(\overline{B},\mathrm{\Delta }B^2)=\frac{B_0^2}{B_0^2+\overline{B}^2}.$$
(22)
For large fluctuations when $`\sqrt{\mathrm{\Delta }B^2}B_0`$ the integration (20) includes a broad range of large values of $`p`$. The probability integral $`\mathrm{\Phi }(x)`$ rapidly converges to one for $`x1`$. Then we may exclude the range $`p<p_0={\displaystyle \frac{B_0\sqrt{2}}{\sqrt{\mathrm{\Delta }B^2}}}`$ whose contribution to the integral is of the order of $`\left({\displaystyle \frac{B_0}{\sqrt{\mathrm{\Delta }B^2}}}\right)^2`$ and integrate for $`p>p_0`$. Then one gets

$$\overline{f}(\overline{B},\mathrm{\Delta }B^2)=\frac{B_0}{2\sqrt{\mathrm{\Delta }B^2}}\sqrt{\frac{\pi }{2}}\left[\text{Ei}(\sqrt{2}\frac{B_0+i\overline{B}}{\sqrt{\mathrm{\Delta }B^2}})+\text{Ei}(\sqrt{2}\frac{B_0i\overline{B}}{\sqrt{\mathrm{\Delta }B^2}})\right]$$

$$\sqrt{\frac{\pi }{2}}\frac{B_0}{\sqrt{\mathrm{\Delta }B^2}}\left[\mathrm{ln}(\frac{\mathrm{\Delta }B^2}{B_0^2+\overline{B}^2})C\frac{1}{2}\mathrm{ln}2\right]$$
(23)
where $`C=0.577`$ is the Euler number. Ei($`x`$) is the integral exponential function.

REFERENCES

Aliev, Kh.K., Kamilov, I.K., and Omarov, A.M., Zh.Eksp.Teor.Fiz., 94, 153 (1988) \[Sov.Phys.JETP, 67, 2262 (1988)\].

Al’shits, V.I., Darinskaya, E.V., Getkina, I.V., and Lavrentiev, F.F., Kristallografiya, 35, 1014 (1990) \[Sov.Phys.Crystallogr., 35, 597 (1979)\]

Bolling, C.F., and Richman, R.H., Philos.Mag., 19, 247, (1969).

Dan’kov, S.Yu., Tishin, A.M., Pecharsky, V.K., and Gschneider, K.A., Jr., Phys.Rev. B 57, 3478 (1998)

Echigoya, J., Hayashi, S., and Yamamoto, Y., Phys.Stat.Sol., (a), 14, 463 (1973).

Echigoya, J., and Hayashi, S., Phys.Stat.Sol., (a), 55, 279 (1979).

Flor, H., Gudladt, H.J., and Schwink, Ch., Acta Metal., 28, 1611 (1980).

Friedel, J.,Dislocations, (Pergamon, Oxford,1964).

Gulayev, A., and Svistunova, E.L., Scr.Mater., 35, 501 (1996).

Haasen, P., in Physical Metallurgy, edited by R.W. Cahn and P.Haasen (North-Holland, Amsterdam, 1983), Part II, pp. 1341-1409.

Hischler, W., and Rocker, W., Z.Angew.Phys., 21, 386 (1966)

Kadanoff, L.P., Götze, W., Hamblen, D., Hecht, R., Lewis, E.A.S., Palciauskas, V.V., Rayl, M., Swift, J., Aspens, D., and Kane, J., Rev.Mod.Phys., 39, 395 (1967)

Kittel, C., Introduction to Solid State Physics, 6-th edition (J.Wiley & Sons, New York) 1986.

Landau, L.D., and Lifshitz, E.M., Statistical Physics, 3rd Edition, Part 1 (Pergamon Press, Oxford, 1980), Ch.14.

Ma, Shang-Keng, Modern Theory of Critical Phenomena, (Benjamin, London, 1976)

Maksimova, G.A., and Maiboroda, V.P., Metal Science and Heat Treatment, 34, 339 (1992).

Molotskii, M., Mater.Sci.Eng., A287, 248 (2000).

Molotskii, M., and Fleurov, V., Phys.Rev.Lett., 78, 2779 (1997).

Molotskii, M., and Fleurov, V., Phys.Rev., B 63, 184421 (2001)

Nabutovskaya, O.A., Fiz.Tverd.Tela (Leningrad), 11, 1434 (1969) \[Sov.Phys. Solid State, 11, 1172 (1969)\].

Nabutovskaya, O.A., Fiz.Tverd.Tela (Leningrad), 13, 3130 (1971) \[Sov.Phys. Solid State, 13, 3130 (1971)\].

Okazaki, K., Kagawa, M., and Conrad, H., Scr.Metal., 13, 473 (1979)

Retat, I., Steffens, Th., and Schwink, Ch., Phys.Stat.Sol;, (a), 92, 507 (1985).

Retat, I., Phys.Stat.Sol., (a), 99, 121 (1981).

Seeger, A., Kronmüller, H., Riegre, H., and Träuble, H., J.Appl.Phys., 35, 740 (1964)

Spichkin, Y.I., Tishin, A.M., and Gschneider, K.A., Jr., J.Magn.Magn.Mater., 204, 5 (1999).

Vonsovskii, S.V., Magnetism, (J.Wiley & Sons, New York, 1974), Vol. 2.

Wolfenden, A., Z.Metallkunde, 69, 308 (1978).

Zackay, V.F., and Hazlett, T.H., Acta Metal., 1, 624 (1953).