text:
# 1 Introduction

## 1 Introduction

Existence of magnetic fields on astronomically large scales remains one of unsolved cosmological mysteries. It is known from observations that there are magnetic fields in galaxies with the strength about micro-gauss which are homogeneous over galactic size $`l_{gal}`$ ( a few) kpc, and maybe, even more puzzling, intergalactic magnetic fields with the strength approximately three orders of magnitude weaker but homogeneous on much larger intergalactic scale of hundreds kpc (for review see e.g. refs. ). Though in stellar processes much stronger magnetic fields can be generated, their coherence scale is negligible in comparison with $`l_{gal}`$. There are two competing ideas of explanation of the phenomenon . The first one is based on traditional astrophysics and pursues the possibility that strong magnetic fields generated in multiple stellar catastrophes could be ejected into interstellar space and by field line reconnection might create magnetic fields coherent on galactic size (for a recent review see ref. ). However there are no compelling quantitative arguments in favor of this hypothesis. Moreover, the energy density of galactic magnetic fields are of the same magnitude as the energy density of cosmic microwave background radiation (CMBR) so its magnitude is about $`10^{10}`$ of the total galactic mass/energy. Such a huge contribution is difficult to explain by conventional mechanisms. It is even more difficult to explain in this way intergalactic magnetic fields, if they exist.

Other possible mechanisms of generation of galactic magnetic fields are based on physical processes in the early universe - for the review see e.g. refs. . One can also find an extensive list of literature in ref. . Basically there are three different mechanisms discussed:
1. Breaking of conformal invariance of electromagnetic interaction at inflationary stage. The latter could be realized either through new non-minimal (and possibly non gauge invariant) coupling of electromagnetic field to curvature , or in dilaton electrodynamics , or by conformal anomaly in the trace of the stress tensor induced by quantum corrections to Maxwell electrodynamics .
2. First order phase transitions in the early universe  producing bubbles of new phase inside the old one. A different mechanism but also related to phase transitions is connected with topological defects, in particular, cosmic strings . A recent discussion and a model of generation of large scale magnetic field can be found in ref. .
3. Creation of stochastic inhomogeneities in cosmological charge asymmetry, either electric , or e.g. leptonic  at large scales which produce turbulent electric currents and, in turn, magnetic fields.

In this work we will consider a new mechanism somewhat similar to those mentioned in point 3 above but in contrast to them this mechanism does not demand any new physics and could be realized at relatively late stage of the universe evolution, namely, at red-shifts $`z10^310^4`$, near hydrogen recombination. The mechanism suggested here is intermediate in time between the early universe ones and the astrophysical mechanism which took place practically in contemporary universe after galaxies were formed and stellar explosions took place.

The basic features of the suggested model of magnetic field generation are the following. We show that, despite low Reynolds number (3), non-zero, though small, vorticity can be generated in the cosmic electron-photon fluid due to inhomogeneity of the latter and especially due to different spectrum of inhomogeneities of electronic and photonic components. Such difference could be created at relatively late stage of cosmological evolution (near hydrogen recombination) even from initially adiabatic density perturbations. The motion of cosmic plasma would create electric currents because of different velocities of electrons and protons and since the vorticity of the motion would be non-vanishing the same would be true for the currents, $`\times 𝐉0`$. Such currents are known to generate magnetic fields. An attractive feature of the proposed mechanism is that no new physics has to be invoked. The mechanism operates in the standard cosmological model in the frameworks of the usual Maxwell electrodynamics. In what follows we will estimate the magnitude of the generated magnetic field and show that it can be strong enough and have sufficiently large coherence scale so that after dynamo amplification (for the review of the latter see refs. ) it can explain the observed fields in the galaxies and, possibly, even intergalactic ones.

## 2 Hydrodynamics of cosmic plasma at recombination

Fluid motion under pressure forces is governed by the hydrodynamical equation (see e.g. the book ):

$`\rho \left(_tv_i+v_k_kv_i\right)=_ip+_k\left[\eta \left(_kv_i+_iv_k{\displaystyle \frac{2}{3}}\delta _{ik}_jv_j\right)+_i\left(\zeta _jv_j\right)\right]`$ (1)
where $`v`$ is the velocity of the fluid element, $`\rho `$ and $`p`$ respectively are the energy and pressure densities of the fluid, and $`\eta `$ and $`\zeta `$ are the first and second viscosity coefficients. In the case of constant viscosity coefficients this equation is reduced to the well known Navier-Stokes equation. The coefficient $`\eta `$ is related to the mean free path of particles in fluid as

$`\eta /\rho \nu =l_f`$ (2)
In what follows we disregard second viscosity $`\zeta `$.

The character of the solution to eq. (1) crucially depends upon the value of the Reynolds number

$`R_\lambda ={\displaystyle \frac{v\lambda }{\nu }}`$ (3)
where $`\lambda `$ is the wavelength of the velocity perturbation. If $`R1`$, then the fluid motion would become turbulent and non-zero vorticity would be created by spontaneously generated turbulent eddies. In the opposite case of low $`R`$ the motion is smooth and in the case of incompressible fluid with homogeneous $`\rho `$ and viscosity coefficients $`\eta `$ and $`\zeta `$ the fluid velocity would have potential character with vanishing vorticity, $`\times 𝐯=0`$. However, the above assumptions of homogeneity is not precise and some vorticity can be generated even with low $`R`$.

Let us consider now the cosmological epoch before the hydrogen recombination, when the plasma consists of three components: photons, electrons, and baryons mainly constituted by protons (for simplicity, we shall neglect small amount of 4He nuclei). For temperatures $`T1`$ eV the energy density of $`(e,p,\gamma )`$-plasma is dominated by photons. Indeed, at the present day the CMBR energy density is $`\rho _\gamma 0.26`$ eV/cm3, i.e. approximately $`10^3`$ smaller than the energy density of baryonic matter $`\rho _b`$, and thus $`\rho _b`$ and $`\rho _\gamma `$ become comparable at red-shift $`z10^3`$ or at $`T0.23`$ eV.111 Throughout this paper, we use the following values of the cosmological parameters: the CMBR temperature $`T_\gamma ^{(0)}=2.725`$ K, the baryon-to-photon ratio $`\beta =n_B/n_\gamma =6\times 10^{10}`$, the Hubble constant $`h=0.71`$, the matter density of the universe $`\mathrm{\Omega }_mh^2=0.135`$, and the baryon density $`\mathrm{\Omega }_bh^2=0.0224`$.  The contribution of dark matter (DM) into energy density is not important for plasma hydrodynamics as far as the DM does not interact with plasma. For the moment, we also neglect cosmological expansion considering characteristic times smaller than the Hubble time $`H^1`$ at the appropriate epoch:

$`H^1={\displaystyle \frac{27\mathrm{kpc}}{T_{\mathrm{eV}}^{3/2}[T_{\mathrm{eV}}+0.76]^{1/2}}},`$ (4)
where $`T_{\mathrm{eV}}=(T/1\mathrm{eV})`$. At the present scale, it corresponds to $`(1+z)H^1=110T_{\mathrm{eV}}^{1/2}[T_{\mathrm{eV}}+0.76]^{1/2}`$ Mpc. This result is obtained for cosmological relativistic matter consisting from photons and neutrinos (contributing 68% with respect to photons); the cosmological constant is not essential. For labeling the cosmological epoch we will interchangeably use the temperature $`T`$ or red-shift $`1+z=4260T_{\mathrm{eV}}`$.

As far as the cosmological plasma was dominated by photons, the viscosity coefficient $`\nu `$ is determined by the photon mean free path:

$`\nu =l_\gamma ={\displaystyle \frac{1}{\sigma _Tn_eX_e}}{\displaystyle \frac{30\mathrm{pc}}{X_e(T)T_{\mathrm{eV}}^3}},`$ (5)
where $`\sigma _T=8\pi \alpha ^2/3m_e^2=6.65\times 10^{25}`$ cm2 is the Thomson cross-section, $`n_e=\beta n_\gamma =610^{10}0.24T^3`$ is the electron number density, and $`X_e`$ is a fraction of the free electrons: $`X_e(z)`$ is practically 1 for $`z>1500`$, and sharply decreases for smaller $`z`$’s, reaching values $`10^5`$ at $`z<1000`$. For $`T1`$ eV the mean free path (5) is much smaller than the Hubble horizon:

$`{\displaystyle \frac{l_\gamma }{H^1}}{\displaystyle \frac{1.1\times 10^3}{X_e(T)T_{\mathrm{eV}}}}\left[1+0.76T_{\mathrm{eV}}^1\right]^{1/2}`$ (6)
This ratio becomes comparable to unity only at temperatures $`T0.3`$ eV, when $`X_e1`$. At the present day scale the photon mean free path (5) would be $`l_\gamma ^{(0)}130\mathrm{kpc}X_e^1T_{\mathrm{eV}}^2`$.

Let us first estimate the Reynolds number of the fluid motion created by the pressure gradient in eq. (1). To this end we assume that the liquid is quasi incompressible and homogeneous, so that the second term in the r.h.s. of this equation can be neglected. This is approximately correct and the obtained magnitude of fluid velocity is sufficiently accurate. In this approximation eq. (1) reduces to a much simpler one:

$`_t𝐯+\left(𝐯\right)𝐯\nu \mathrm{\Delta }𝐯={\displaystyle \frac{p}{\rho }}`$ (7)
A comment worth making at this stage. The complete system of equations includes also continuity equation which connects the time variation of energy density with the hydrodynamical flux (see below eq. (19) and the Poisson equation for gravitational potential induced by density inhomogeneities. We will however neglect the gravitational force and the back reaction of the fluid motion on the density perturbation. This approximation would give a reasonable estimate of the fluid velocity for the time intervals when acoustic oscillations are not yet developed, i.e. for $`t<\lambda /v_s`$, where $`\lambda `$ is the wave length of the perturbation and $`v_s`$ is the speed of sound (in the case under consideration $`v_s^2=1/3`$). In fact the wave length should be larger than the photon mean free path, to avoid diffusion damping, and the characteristic time interval, as we see below, should be somewhat larger than $`\lambda `$. So we may hope that our estimates of the velocity are reasonable enough. Neglecting gravitational forces, especially those induced by dark matter would result in a smaller magnitude of the fluid velocity, so the real effect should be somewhat larger.

For small velocities (or sufficiently small wavelengths) we may neglect the second term in the l.h.s. with respect to the third one. In this approximation the equation becomes linear and can be easily solved for the Fourier transformed quantities. Assuming that the parameters are time-independent (though it is not necessary) we obtain:

$`𝐯_𝐤={\displaystyle \frac{i𝐤}{3k^2\nu }}\delta _k\left[1\mathrm{exp}(\nu k^2t)\right]`$ (8)
where $`\delta _k=(\delta \rho /\rho )_k`$ is the Fourier transform of relative density perturbations, $`\delta \rho /\rho `$; its natural value is $`10^4`$, though it might be much larger at small scales. The coefficient $`1/3`$ comes from equation of state of relativistic fluid, $`p=\rho /3`$.

Therefore, for the Reynolds number we obtain:

$`R_k={\displaystyle \frac{\delta _k}{3\left(k\nu \right)^2}}\left[1\mathrm{exp}\left(\nu k^2t\right)\right].`$ (9)
If $`\delta _k`$ is weakly dependent on $`k`$, then $`R_k`$ is a monotonically rising function of the wavelength $`\lambda =2\pi /k`$. For $`t\lambda ^2/\nu `$ it takes the value

$`R_k={\displaystyle \frac{t}{3\nu }}\delta _k`$ (10)
Therefore, considering the perturbation with the wavelength $`\lambda =2\pi /k`$ which enters the horizon at the time $`tH^1`$, we find that the maximum value of the Reynolds number is given by $`R_k^{\mathrm{max}}(H^1/3l_\gamma )\delta _k`$. For later times, the Reynolds number for the comoving wavelength decreases approximately as $`T`$ until the density perturbation is completely damped by the photon diffusion. Thus, for $`T1`$ eV, using eq. (6), we obtain $`R_k^{\mathrm{max}}300\delta _kT_{\mathrm{eV}}`$, and so for the development of turbulence one needs $`\delta _kT_{\mathrm{eV}}>0.1`$. This condition can be satisfied either at high temperatures, $`T100`$ eV, or for large density perturbations, $`\delta _k10^4`$. This is not the case for the adiabatic fluctuations with a nearly flat spectrum, as is generically predicted in inflationary scenarios. However, one cannot exclude a situation that at smaller wavelengths, corresponding to the present day scales $`100`$ kpc, there are substantial isocurvature fluctuations, e.g. the baryon density perturbations $`(\delta \rho /\rho )_b`$ are larger than the dark matter fluctuations $`(\delta \rho /\rho )_{\mathrm{DM}}`$.

The baryon density perturbation of the wavelength $`\lambda `$, after many oscillations, will be damped by the Silk effect  at the effective timescale $`t_{\mathrm{eff}}\lambda ^2/\nu `$. However, because of a large Reynolds number, it could happen that during this time the photon diffusion dragging electrons would also produce non-zero vorticity.222 The case of large Reynolds numbers, that might be created at earlier epoch, $`T1`$ MeV, due to large leptonic asymmetries, was considered in ref. .  In this case the generation of magnetic fields might be strongly amplified. We shall keep in mind an interesting possibility of large $`\delta _k`$ at smaller scales, however now we shall mostly concentrate on a more plausible possibility of small Reynolds numbers and non-turbulent flow of cosmic fluid in $`T1`$ eV range of temperatures.

At first sight vorticity in laminar fluid motion is not generated in the approximation given by eq. (7). However, this is not quite so. If both $`p`$ and $`\rho `$ are different functions of space points, pressure gradient may create motion with non-vanishing $`\times 𝐯`$. Indeed, from eq. (7) follows

$`_t𝛀\nu \mathrm{\Delta }𝛀=\times \left({\displaystyle \frac{p}{\rho }}\right)`$ (11)
where $`𝛀=\times 𝐯`$ and we assume that velocity is small so that the term quadratic in $`v`$ was neglected. If the r.h.s. is non-vanishing, then $`𝛀`$ would be non-zero too. However usually pressure density is proportional to the energy density, $`p=w\rho `$, with a constant coefficient $`w`$ and hence $`\times \left(p/\rho \right)=0`$.

In the case under consideration, electrons are strongly coupled to photons. Their effective mean free path is

$`l_e={\displaystyle \frac{\sqrt{3m_e/T}}{\sigma _{\mathrm{Th}}n_\gamma }}`$ (12)
The factor in the numerator takes into account large thermal momentum of electrons, $`p\sqrt{3m_eT}`$, so that they need many collisions, with momentum transfer $`\mathrm{\Delta }p=T`$, for a significant change of their momentum. Still despite the presence of this factor, the mean free path of electrons is much smaller than that of photons simply because there are more than billion photons per a single electron in the plasma. So electrons are practically frozen in the plasma, while photons may diffuse to much larger distance. One should keep in mind that the plasma must be locally electrically neutral and thus redistribution of electron inhomogeneities is accompanied by the same redistribution of baryonic ones. This makes the mean free path of charged species of matter even smaller. On the other hand, a motion of homogeneous component of electron number density does not require dragging of baryons.

For the wavelengths larger than mean free path of charged particles we may consider plasma as a fluid where charged particles are strongly coupled to photons. So for the estimate of the rotational velocity we will use the hydrodynamic equation (11). As we have already mentioned, vorticity, $`\times 𝐯`$ could be non-zero only if the fluid is non-homogeneous and strictly speaking we should use eq. (1) with all parameters depending upon space points. However, we believe that for a simple estimate eq. (11) may be sufficient. Since before recombination the interaction rate between radiation and charged particles is very high, plasma/liquid should be in local thermal equilibrium with all constituents having the same temperature $`T(x)`$. If $`T`$ would be the only parameter which determines the state of the medium, then vorticity would not be generated because we would have in our disposal only $`T`$ and it is impossible to construct non-vanishing $`\times 𝐯`$ from the gradient of only one scalar function. However, distributions of charged particles depend upon one more function, their chemical potential:

$`f=\mathrm{exp}\left[{\displaystyle \frac{E}{T(x)}}+\xi (x)\right]`$ (13)
where the dimensionless chemical potential $`\xi `$ can be readily expressed through particle number density $`n_en_B=\beta (x)n_\gamma `$ with $`\beta (x)=610^{10}+\delta \beta (x)`$:

$`\xi (x)=\mathrm{ln}\beta (x)+\mathrm{const}`$ (14)
Hence we will find that the source term in the vorticity equation (11) is equal to:

$`S_kϵ_{ijk}_j\left({\displaystyle \frac{_ip}{\rho }}\right)=ϵ_{ijk}{\displaystyle \frac{_i\rho _\gamma }{3\rho _{\mathrm{tot}}}}{\displaystyle \frac{_j\beta }{\beta }}{\displaystyle \frac{\rho _b}{\rho _{\mathrm{tot}}}}`$ (15)
An essential feature here is that the spatial distribution of charged particles does not repeat the distribution of photons and hence the vectors $`\rho _\gamma `$ and $`\beta `$ are not collinear. This could occur if for lambda corresponding to subgalactic scales, there exist baryon isocurvature fluctuations and thus $`\rho (x)`$ and $`\beta (x)`$ have different profiles. As we have mentioned above, different mean free paths of photons and charged particles would maintain such non-collinearity of the order of unity at the scales $`\lambda l_\gamma `$. Moreover, even in the case of adiabatic perturbations a shift in the distribution of photons and charged particles could also be created because of acoustic oscillations that proceeded with different phases of radiation and matter densities. At the scales $`\lambda l_\gamma `$ perturbations in the the plasma temperature would be erased by the diffusion damping , while for $`\lambda l_\gamma `$ diffusion processes are not efficient and one would expect self-similar perturbation leading to collinearity of $`\rho _\gamma `$ and $`\beta `$. On the other hand, when $`\lambda `$ entered under horizon acoustic oscillations begun which destroyed the self-similarity. Thus the expected wavelengths of vorticity perturbations should be between $`l_\gamma <\lambda <H^1`$.

We assume that there is no additional suppression of the source term (15) and by the order of magnitude its amplitude corresponding to the wavelength $`\lambda `$ can be evaluated as

$`S{\displaystyle \frac{16\pi ^2}{3\lambda ^2}}\left({\displaystyle \frac{\delta T}{T}}\right)_\lambda \left({\displaystyle \frac{\delta \beta }{\beta }}\right)_\lambda \left({\displaystyle \frac{\rho _b}{\rho _{tot}}}\right)10\left({\displaystyle \frac{\delta T}{T}}\right)_\lambda \left({\displaystyle \frac{\delta \beta }{\beta }}\right)_\lambda T_{\mathrm{eV}}^1`$ (16)
Taking $`\delta T/T\delta \beta /\beta 310^5`$ we obtain $`S10^8T_{\mathrm{eV}}^1`$.

Now equation (11) can be solved in the same way as eq. (7) and we find:

$`|𝛀|={\displaystyle \frac{0.27}{l_\gamma T_{\mathrm{eV}}}}\left({\displaystyle \frac{\delta T}{T}}\right)_\lambda \left({\displaystyle \frac{\delta \beta }{\beta }}\right)_\lambda \left[1\mathrm{exp}\left({\displaystyle \frac{4\pi ^2l_\gamma t}{\lambda ^2}}\right)\right]`$ (17)
with $`l_\gamma `$ given by eq. (5).

Vorticity could also be generated even if perturbations in plasma are determined by a single scalar function, for example, by $`T(t,𝐱)`$ because if might be proportional to the product $`_iT(t,𝐱)_jT(t^{},𝐱)`$. These two gradients generally are not collinear if taken at different time moments $`t`$ and $`t^{}`$. To see that, let us start from the Boltzmann equation for the distribution function $`f(t,𝐱,E,𝐩)`$ of photons:

$$\left(\frac{}{t}+𝐕H𝐩\frac{}{𝐩}+𝐅\frac{}{𝐩}\right)f(t,𝐱,E,𝐩)=I_{\mathrm{coll}}\left[f_a\right],$$
(18)
where $`𝐕=𝐩/E`$ is the particle velocity (not to be confused with the velocity $`𝐯`$ of macroscopic motion of the medium), for photons $`V=1`$, while $`v1`$, $`E`$ and $`𝐩`$ are respectively the particle energy and spatial momentum, $`H`$ is the universe expansion rate, $`𝐅`$ is an external force acting on particles in question (the latter is assumed to be absent), and $`I[f_a,f_b,\mathrm{}]`$ is the collision integral depending on the distributions $`f_a`$ of all participating particles (for the definition of the $`I_{\mathrm{coll}}`$ see e.g. eq. (43) of review ).

At temperatures in eV-range only the Thomson scattering of photons on electrons is essential, so the collision integral is dominated by the elastic term. Integrating both parts of eq. (18) over $`d^3p/(2\pi )^3`$ we arrive to the continuity equation:

$`\dot{n}(𝐱)+𝐉=0`$ (19)
where $`𝐉`$ is the photon flux given by

$`𝐉𝐯n={\displaystyle \frac{d^3p}{(2\pi )^3}\frac{𝐩}{E}f}`$ (20)
and $`𝐯`$ is the average macroscopic velocity of the photon plasma. Using the standard arguments one can derive from eq. (19) the diffusion equation:

$`\dot{n}=D\mathrm{\Delta }n`$ (21)
where $`Dl_\gamma /3`$ is the diffusion coefficient. We will use this equation below to determine time evolution of the photon temperature $`T`$.

If the elastic reaction rate $`\mathrm{\Gamma }_{\mathrm{el}}=\sigma _{\mathrm{Th}}n_eX_e=1/l_\gamma `$ is sufficiently large, local thermal equilibrium would be established and the photon distribution would be approximately given by

$`ff_0=\mathrm{exp}\left(E/T+\xi \right)`$ (22)
where the temperature and effective chemical potential could be functions of time and space coordinates: $`T=T(t,𝐱)`$ and $`\xi =\xi (t,𝐱)`$, and the photon mean free path is given by eqs. (5,6). Evidently $`f_0`$ annihilates the collision integral. We can find correction to this distribution, $`f=f_0+f_1`$, substituting this expression into kinetic equation (18) and approximating the collision integral in the usual way as $`\mathrm{\Gamma }_{\mathrm{el}}f_1`$:

$`\left(K+\mathrm{\Gamma }_{\mathrm{el}}\right)f_1=Kf_0`$ (23)
where $`K`$ is the differential operator, $`K=_t+\left(𝐕\right)`$. The solution of this equation is straightforward:

$`f_1(t,𝐱,E,𝐕)={\displaystyle _0^t}𝑑\tau _1\mathrm{exp}\left[{\displaystyle _0^{\tau _1}}𝑑\tau _2\mathrm{\Gamma }_{el}(t\tau _2,𝐱𝐕\tau _2)\right]Kf_0(t\tau _1,𝐱𝐕\tau _1)`$ (24)
Using this result we can calculate the average macroscopic velocity of the plasma. The calculations are especially simple if elastic scattering rate is high and the integrals are dominated by small values of $`\tau _1`$. In this case we obtain:

$`v_j(t,𝐱)={\displaystyle \frac{d^3pV_jf_1(t,𝐱,E,𝐕)}{d^3pf_0(t,𝐱,E)}}`$ (25)
and its vorticity, $`\mathrm{\Omega }_i=ϵ_{ijk}_jv_k`$:

$`\mathrm{\Omega }_i6ϵ_{ijl}l_\gamma ^2\left({\displaystyle \frac{_jT}{T}}\right)\left({\displaystyle \frac{_l_tT}{T}}\right)`$ (26)
To estimate time derivatives of the temperature we will use the diffusion equation (21), from which we find $`_tT=D\mathrm{\Delta }T`$ and finally obtain for vorticity at the scale $`\lambda `$:

$`|𝛀|_\lambda 2\left({\displaystyle \frac{\delta T}{T}}\right)_\lambda ^2l_\gamma ^3k^4310^3\left({\displaystyle \frac{\delta T}{T}}\right)_\lambda ^2{\displaystyle \frac{l_\gamma ^3}{\lambda ^4}}`$ (27)
Since the photon diffusion erases temperature fluctuations at the scales $`\lambda <l_\gamma `$ the vorticity has the maximum near $`\lambda l_\gamma `$. This magnitude of vorticity is considerably larger than found previously (17) and we will rely on it in the estimates of magnetic field presented in the following section.

To avoid confusion let us mention that at the moment when perturbation with a given wave length enters the horizon, $`\lambda l_\gamma `$. Later $`(\lambda /l_\gamma )`$ scales as $`a(t)^2`$ and till $`\lambda `$ remains larger than $`l_\gamma `$ the amplitude of perturbations does not decrease significantly and only when $`\lambda l_\gamma `$ the photon diffusion damps density perturbations.

## 3 Electrodynamics of cosmic plasma with helical flows and magnetic field generation

Generation and evolution of magnetic field strongly depends upon the electric conductivity of plasma which can be estimated as follows. Equation of motion of a charged particle in external electric field $`𝐄`$ has the usual form

$`m\dot{𝐕}_d=e𝐄`$ (28)
So the drift velocity gained during the time $`\mathrm{\Delta }t`$ is equal to $`V_d=eE\mathrm{\Delta }t/m`$. The charged particles (electrons) keep on to be accelerated approximately during time between collisions, $`\mathrm{\Delta }t=l_e/V_T`$ where $`V_T=\sqrt{3T/m}`$ is the thermal velocity and the electron mean free path $`l_e`$ is given by eq. (12). After collision the particle loses (forgets) its previous velocity and the process repeats. This is true for sufficiently weak fields when the drift velocity is small in comparison with thermal velocity, otherwise run-away charge carriers would be produced and the conductivity would be much larger.

The induced electric current is $`𝐉=en_e𝐕`$, where $`n_e`$ is the number density of charge carriers. Comparing this with the definition of current conductivity $`\kappa =J/E`$ we find:

$`\kappa ={\displaystyle \frac{3}{2\alpha }}{\displaystyle \frac{n_e}{n_\gamma }}{\displaystyle \frac{m_e^2}{T}}`$ (29)
The conductivity is very high, so the generation of magnetic field by the source currents, created by the cosmological inhomogeneities, is governed by the well know equation of magnetic hydrodynamics:

$`_t𝐁=\times \left(𝐯\times 𝐁\right)+{\displaystyle \frac{1}{\kappa }}\times 𝐉`$ (30)
The electric current $`J`$ induced by the relaxation of the density inhomogeneities would contain two components: electronic and protonic. However the first one is surely dominant because it is much easier to drift electrons than heavier protons. This is why a non-zero current can be induced. Of course motion of electrons should not produce any excess of electric charge but this can be realized because the current is created in the dominant homogeneous part of the charge particle distribution.

The solution of eq. (30) can be roughly written as

$`B{\displaystyle _0^t}𝑑t_1\left({\displaystyle \frac{2\pi J}{\lambda \kappa }}\right)e^{2\pi vt_1/\lambda }`$ (31)
The exponential factor under the integral, which presents pregalactic dynamo effect is normally rather weak, the exponent is about

$`{\displaystyle \frac{2\pi vt}{\lambda }}500T_{\mathrm{eV}}\left({\displaystyle \frac{\delta T}{T}}\right)_\lambda `$ (32)
where eqs. (8) and (6) have been used and the density fluctuations were taken as $`\delta _\lambda =4\left(\delta T/T\right)_\lambda `$. With $`\left(\delta T/T\right)_\lambda 310^5`$ and $`T`$ eV the exponent is about $`0.015T_{\mathrm{eV}}<1`$ and pregalactic dynamo is not important. It may be significant if the density perturbations at relatively small scales are much larger than their accepted canonical values. Exponential enhancement can be large at higher temperatures, $`T>10^3`$ eV, which correspond to the present day scales below kpc. Equipartition between magnetic field and CMBR can be expected on such scales and after ’Brownian” line reconnection and relatively mild galactic dynamo the observed galactic magnetic fields may be created.

An estimate of magnetic field without pregalactic dynamo enhancement can be easily done if the helical source current is known, $`\times 𝐉=en_e𝛀`$. With $`𝛀`$ given by eq. (27) and $`B`$ by eq. (31) we obtain:

$`{\displaystyle \frac{B_0}{T^2}}=0.2410^3\left(4\pi \alpha \right)^{3/2}\left({\displaystyle \frac{t}{\lambda }}\right)\left({\displaystyle \frac{l_\gamma }{\lambda }}\right)^3\left({\displaystyle \frac{T}{m_e}}\right)^210^8T_{\mathrm{eV}}^3`$ (33)
where we took the wavelength equal to the photon mean free path, $`\lambda =l_\gamma `$.

If we take into account that linear compression of pregalactic medium in the process of galaxy formation is approximately $`r10^2`$, the seed field in a galaxy after its formation would be $`r^2B_0`$, i.e. 4 orders of magnitude larger than that given by eq. (33) and, for $`T=1`$ eV, a relatively mild galactic dynamo, about $`10^4`$, is necessary to obtain the observed galactic magnetic field of a few micro-Gauss at the scale $`l_B(100/r)`$ kpc $`=1`$ kpc. The seed magnetic fields formed earlier (at higher $`T`$) would have larger magnitude ($`T^3`$) but their characteristic scale would be smaller by factor $`1/T^2`$. Chaotic line reconnection could create magnetic field at larger, galactic scale $`l_{gal}`$, but the magnitude of this field would be suppressed by Brownian motion law - it would drop by the factor $`(l_B/l_{gal})^{3/2}`$. It is interesting that according to our results all scales give comparable contributions at $`l_{gal}`$. This effect may lead to an enhancement of the field but it is difficult to evaluate the latter. Let us also note that magnetic fields generated by the discussed mechanism at the cluster scale, 10 Mpc, should be not larger than $`10^8`$ $`\mu `$G if no additional amplification took place.

## 4 Discussion and conclusion

We have shown in this work that photon diffusion produced by temperature fluctuations of CMBR could create vortical flows in the primeval plasma in the second order, i.e. proportional to $`(\delta T/T)^2`$. Because of different mobility of electrons and protons (or helium ions) in the plasma the helical macroscopic motion of the latter would create helical electric currents which would in turn generate magnetic fields. The amplitude of such fields generated at the moment when plasma temperature was equal to $`T`$ is given by eq. (33). The characteristic scale of the field should be equal to the photon mean free path in the plasma at temperature $`T`$ which at the present day is about 100 kpc $`(\mathrm{eV}/T)^2`$. This scale is quite close to the galactic size. Adiabatic compression in the process of galaxy formation by the factor $`r100`$ would lead to an enhancement of magnetic field by $`r^2=10^4`$. Thus to meet the observational data this seed field should be amplified by the galactic dynamo only by 4-5 orders of magnitude. Even if the magnitude of the seed field presented above is somewhat overestimated there is a huge reserve in the amplification by the galactic dynamo which could be as large as 15 orders of magnitude .

A nice feature of this mechanism is that it does not demand any new physics for the realization. Of course, density/temperature perturbations at large scales were, most probably, created during inflationary epoch and in this sense new physics is invoked but independently of the mechanism we know from observations that $`\delta T/T0`$ and, based on these data, we may estimate the seed magnetic fields using good old physics.

It is interesting that, according to the existing indications to intergalactic magnetic fields, their strength is 3-4 orders of magnitude weaker than the strength of galactic fields. This fact (if it is a fact) hints that intergalactic and galactic magnetic fields might have common origin, but galactic fields are larger due to mentioned above adiabatic compression which enhanced the field by 4 orders of magnitude. If it is so, the mechanism discussed in this work may be irrelevant because some galactic dynamo is needed to amplify the galactic seed field up to the observed magnitude. On the other hand, such dynamo seemingly does not operate on intergalactic scales.

Larger density perturbations would be helpful for generation of larger magnetic field for which dynamo might be unnecessary. Though much bigger $`\delta T`$ is not formally excluded at the scale about 100 kpc, but to have them at the level $`(\delta T/T)^210^4`$ seems to be too much. A natural idea is to turn to a later stage, to onset of structure formation when $`\delta \rho /\rho `$ becomes larger than $`10^2`$. With such density perturbations strong enough magnetic fields may be generated without dynamo amplification. However after recombination the number density of charge carriers drops roughly by 5 orders of magnitude. Correspondingly $`l_\gamma `$ rises by the same amount and the strength of the seed field would be 5 orders of magnitude smaller if density perturbations and the temperature of formation remained the same. However both became very much different. Density perturbations rose as scale-factor, $`(\delta \rho /\rho )^2(T_{eq}/T)^2`$, where $`T_{eq}1`$ eV is the temperature when radiation domination changed into matter domination and density perturbations started to rise. Since, $`B/T^2T^3`$, according to eq. (33), the net effect of going to smaller $`T`$ is a decrease of $`B/T^2`$ which would be difficult to cure even by later reionization. Still, as argued in ref. , magnetic field generation, driven by anisotropic and inhomogeneous radiation pressure (and in this sense similar to our mechanism) at the epoch of reionization, could end up with the field of about $`810^6\mu \mathrm{G}`$. This result is 8 orders of magnitude larger than that found in the earlier papers  and quite close to ours (33), though these two mechanisms operated during very different cosmological epochs and were based on different physical phenomena.

After this paper was sent to astro-ph we became aware of the work  ”Generation of Cosmic Magnetic Fields at Recombination” where a much weaker effect was found. We think that the difference can be attributed to the following effects. We considered earlier period when the photon mean free path was much smaller than the horizon. It gives a factor about $`10^3`$ in fluid velocity, eq.(8). Moreover, since in our case the electrons are tightly bound to photons the electron-photon fluid moves as a whole (while protons and ions are at rest) and the electric current induced by macroscopic motion/oscillations of plasma is noticeably larger.

Acknowledgments

We thank C. Hogan and P. Naselsky for critical comments. This work was partially supported by the MIUR research grant ”Astroparticle Physics”.


format: text

text:
# On the Coalescence of Nanoscale Metal Clusters

## Abstract

We study the coalescence of nanoscale metal clusters in an inert-gas atmosphere using constant-energy molecular dynamics. The coalescence proceeds via atomic diffusion with the release of surface energy raising the temperature. If the temperature exceeds the melting point of the coalesced cluster, a molten droplet forms. If the temperature falls between the melting point of the larger cluster and those of the smaller clusters, a metastable molten droplet forms and freezes.

The ability to predict and control the size and structure of nanoparticles is of great importance for applications, as the arrangement of atoms is a key determinant of the optical, electronic and thermal properties of such particles. However, the structure of an atomic cluster often differs from that of the corresponding bulk material Ino (1966). As the number of surface atoms is comparable to the number of interior atoms in a cluster, the surface energy plays an important role in determining the overall structure. The delicate balance between surface and internal energies often produces a complex dependence of equilibrium structure upon cluster size Raoult et al. (1989); Cleveland et al. (1997).

In an inert-gas aggregation (IGA) source, where metal clusters are produced by evaporation and then condensation from a vapor in an inert gas atmosphere, experiments have shown that the final cluster structure depends on the growth rate Reinhard et al. (1997). Molecular dynamics simulations of silver clusters support this finding, with studies of both the growth of clusters from small seeds Baletto et al. (2000) and the freezing of clusters from molten droplets Baletto et al. (2002) finding a dependence of structure on kinetic factors. These examples illustrate that it is important to understand cluster kinetics if one wishes to predict cluster structure.

The growth of clusters by coalescence is a kinetic process which will also effect structure. Indeed, the aggregation of atomic clusters has generated much interest due to the variety of different large-scale structures that can form, and the relevance of this to the functionality of devices built from clusters. For example, the deposition of clusters onto surfaces can result in highly ramified fractal structures Brechignac et al (2001), highly symmetric droplets Goldby et al (1996), as well as many intermediate structures Kaiser et al (2002). A recent proposal Schmelzer et al. (2002) uses cluster deposition to form wire-like cluster chains resembling nanowires with conduction of electrical current between neighboring particles. Understanding cluster coalescence is therefore crucial for reliable design and fabrication of such devices.

Further, the distribution of cluster sizes produced in IGA sources suggests that coalescence is an important growth mechanism Granqvist and Buhrman (1976, 1976). Studies by electron diffraction of lead clusters grown in an IGA source have found structures that could not be fitted by standard model structures but which might be accounted for by coalesced structures Hyslop et al. (2001). To predict the structures of metal clusters produced in IGA sources we need to study the coalescence of free metal clusters.

The collision and coalescence of small clusters, containing tens of atoms, has been studied extensively using molecular dynamics. For example, an early study showed that the decrease in surface area during fusion can lead to an increase in kinetic energy Schmidt et al. (1991). Simulations of the coalescence of larger clusters have been predominantly conducted at constant temperature Zhu and Averback (1996); Lewis et al. (1997), which probably best represents the coalescence of supported clusters with the substrate acting as a thermostat Lewis et al. (1997). However, as we will argue below, the coalescence of nanometer-sized free clusters in an IGA source is best described by constant energy molecular dynamics. We are aware of only one such study, which considered the fusion of both pairs of liquid and pairs of ”glassy” solid silicon clusters Zachariah and Carrier (1999), and demonstrated that the kinetic energy generated during coalescence can considerably shorten coalescence times by enhancing atomic diffusion. Recent kinetic Monte Carlo simulations have shown, however, that the coalescence of facetted clusters can be considerably slower than surface-roughened particles Combe et al. (2000).

Our purpose is to study the coalescence of facetted solid lead clusters using molecular dynamics under conditions similar to those found in an IGA source. We examine the coalescence of lead clusters using the empirical glue potential developed by Lim, Ong and Ercolessi Lim et al. (1992). Recently it was shown that many-body effects in the glue potential produce unusual surface-reconstructed icosahedra in simulations of freezing Hendy and Hall (2001). Further, these surface-reconstructed icosahedra are thought to be globally optimum structures for the glue-potential for lead over a wide range of sizes Hendy and Hall (2001); Doye and Hendy (2003); Hendy and Doye (2002). It is of interest then to study the coalescence of such icosahedra under conditions that might occur in an IGA source.

An estimate of the mean time, $`\tau _g`$, for a cluster of radius $`R`$ between collisions with an inert-gas atmosphere of atomic mass $`m_g`$, pressure $`P`$ and temperature $`T`$ in the source can be made as follows Baletto et al. (2002):

$$\tau _g\frac{1}{PR^2}\sqrt{\frac{m_gkT}{8\pi }}.$$

For a 2 nm radius lead cluster in a helium atmosphere at 5 mbar and 500 K, the time between collisions is approximately 1 ns. Thus over periods of a nanosecond or so, a constant-energy simulation of cluster coalescence will be a reasonable approximation to conditions in an IGA source. We will use constant-energy molecular dynamics here.

Two 565-atom surface-reconstructed icosahedral clusters are initially prepared by construction as detailed in Hendy and Doye Hendy and Doye (2002). Each of these clusters consist of a 147 Mackay icosahedral core with two reconstructed outer shells which differ from the usual Mackay icosahedral surface termination. Each cluster is then equilibrated in isolation during a constant-temperature simulation for $`10^5`$ time-steps at some common initial temperature $`T_i`$. The clusters are then placed in a common simulation cell so that the edges are aligned and the corresponding edge atoms are 0.5 nm apart(approximately next-nearest neighbor distance but within the glue potential cut-off). In this work, we will not examine the effects of relative orientation, initial structure or impact parameter of the clusters in detail. Further simulations reveal that changing relative orientations does not alter the conclusions reached here, although we have not investigated the case of non-zero impact parameter.

The constant-energy simulation then proceeds at a total energy equal to the sum of the final energies from each isolated equilibration run. Here we use a 3.75 fs time step which ensures energy is conserved to within 0.01 eV/atom over a 10 ns simulation. Figure 1 shows a series of snapshots of the initial stages of coalescence between two 565-atom clusters each equilibrated at 300 K. Coalescence typically begins tens of picoseconds into the simulation when an atom on one of the aligned edges makes contact with the adjacent edge of the other cluster. The first snapshot in figure 1 shows the situation at this initial moment of contact, with the formation of an initial neck shown in subsequent frames. Neck growth is essentially complete by the final snapshot in figure 1. Figure 2 shows the evolution of temperature and the aspect ratio 111Aligning the clusters along the x-axis, we define the aspect ratio to be $`(R_y+R_z)/2R_x`$ where $`R_\alpha `$ is the radius of gyration along axis $`\alpha =x,y,z`$.

As the neck grows, the surface area of the merging clusters decreases, leading to the release of considerable surface energy. This energy is converted to heat, as can be seen in Figure 2 by the large increase in temperature during neck growth. After initial contact and the period of rapid neck growth, the aspect ratio evolves more slowly. In Figure 2, we can see that the cluster is essentially spherical in shape after a period of 0.5 ns after initial contact, as the aspect ratio is very close to one by this stage. Examination of cluster structure reveals that it is static after this time. We conclude that coalescence is complete within 1 ns. The final structure of the 1130-atom cluster is that of a surface-reconstructed icosahedron (see Figure 2).

The amount of surface energy released as heat will clearly play a role in enhancing atomic diffusion during coalescence. In fact, if the initial temperature of the two clusters is close enough to the melting point, the final coalesced cluster may be molten. Such a case is shown in Figure 3, which illustrates the evolution of temperature and aspect ratio of two clusters initially at 420 K. Here the coalescence is far more rapid than in the previous example, occurring in tens of picoseconds rather than over 0.5 ns. The clusters form a neck and then melt, with the coalesced cluster rapidly becoming spherical. Note the initial spike in temperature prior to melting, followed by a drop in temperature due to the heat of melting. The final cluster is molten at a temperature of approximately 480 K. In the IGA source, this coalesced cluster will cool slowly and freeze before reaching equilibrium with the inert gas atmosphere.

It is instructive to consider simple models of coalescence and melting. Consider the coalescence of a spherical particle of radius $`R_1`$ with a spherical particle of radius $`R_2R_1`$. Conservation of energy tells us that the increase in thermal energy of the coalescing cluster must balance the reduction in surface energy. Thus, the increase in temperature of the aggregate can be estimated as:

$$\mathrm{\Delta }T=\frac{3\sigma }{\rho c_v}\frac{1}{R_2}\frac{(1+(R_1/R_2)^2)(1+(R_1/R_2)^3)^{2/3}}{(1+(R_1/R_2)^3)}$$
(1)
where $`\sigma `$ is the surface tension, $`c_v`$ is the heat capacity and $`\rho `$ is the density of the two particles. In general, $`\sigma `$, $`c_v`$ and $`\rho `$ will depend on cluster size. However, we will ignore this dependence here and use the bulk values for these parameters. For the simulation presented in figure 2, $`R_1=R_2=R(565)=1.6`$ nm and equation (1) predicts $`\mathrm{\Delta }T=160`$ K, which is close to the value attained in the simulation of 150 K. This increase in temperature is maximized when the ratio of change in surface area to volume ratio is maximized i.e. when $`R_1=0.852R_2`$ for a given $`R_2`$. If $`T_i+\mathrm{\Delta }T`$ exceeds the cluster melting point then we would expect the cluster to melt during coalescence.

However, it is well-known that the melting point, $`T_M(R)`$, of a metal cluster depends strongly on the radius, $`R`$. A simple model for this dependence can be written as Pawlow (1909):

$$T_M(R)=T_B\left(1\frac{L}{R}\right)$$
(2)
where $`L`$ is some length which depends on bulk material parameters. Obviously, this model fails for clusters of radius $`RL`$ and does not account for the irregular variation in melting point observed in small clusters Kusche et al. (1999), but experimental results for lead Ben David et al. (1995), where $`L0.5`$ nm, suggest it is adequate for $`R2`$ nm. By simulation we have identified the melting point of the initial 565-atom clusters to be $`T_M(565)=400\pm 10`$ K and that of the coalesced 1130-atom clusters to be $`T_M(1130)=480\pm 10`$ K. Equation (2) predicts the melting point for the 565-atom clusters ($`R=1.6`$ nm) to be 410 K and for the 1130-atom clusters ($`R=2.0`$ nm) to be 450 K. This rough agreement between (2) and simulation will be sufficient for our purposes here.

Clearly one can always arrange an initial temperature $`T_i`$, such that $`T_i+\mathrm{\Delta }T>T_M((R_1^3+R_2^3)^{1/3})`$. However, for clusters of equal size $`R_1=R_2`$, comparison of equation (1) and equation (2) suggests that there is a critical cluster size $`R=1`$ nm, below which coalescence will always cause melting (i.e. since at $`T_i=0`$, $`\mathrm{\Delta }T>T_M((R_1^3+R_2^3)^{1/3})`$). However, we have not tested this using simulations here as the glue potential Lim et al. (1992) may be unsuitable for such small clusters Doye and Hendy (2003).

Figure 4 shows the final temperatures of a series of pairs of 565-atom clusters that have undergone coalescence. The initial temperature of the pairs has been varied from 250 K to 400 K. The solid line is the value of $`\mathrm{\Delta }T`$ predicted by equation (1). We see from the figure, that equation (1) is adequate (considering the approximations involved) until an initial temperature of $`T_i=340K`$. Above this initial temperature the final state of the coalesced clusters is molten. The highest final temperature reached in this series of simulations is $`T_i+\mathrm{\Delta }T=478K`$. Thus when the melting point of the final cluster is exceeded during coalescence, a molten cluster results. The drop in the final temperature is due to the latent heat of melting. Note, that as the initial temperatures of the simulations correspond to the total energy, figure 4 is essentially a caloric curve for the coalesced clusters. The coexistence of both liquid and solid clusters between final temperatures of 425 K and 475 K corresponds to the van der Waals loop in the microcanonical ensemble Wales and Berry (1994).

In figure 5 we illustrate the coalescence of two 565-atom icosahedra initially prepared at 330 K. Coalescence occurs rapidly, with the aspect ratio reaching one 0.2 ns after initial contact. There is a corresponding increase in temperature of approximately $`\mathrm{\Delta }T90K`$ at this time. Inspection of the cluster structure at this time reveals that it is highly disordered and is consistent with that of a molten structure 222In fact, it is difficult to be certain whether this transient state is liquid or a disordered glassy structure.. The smaller than expected increase in temperature during coalescence also suggests the cluster is molten. However this structure is evidently unstable as at $`t`$ 0.5 ns, the cluster structure undergoes a transition to form a more regular icosahedral structure. There is a corresponding increase in temperature at this time. This route to coalescence via an unstable molten intermediate state was found to be typical of coalescing clusters at initial temperatures of 320-330 K. It is interesting to note that the temperature, $`T=420`$ K, of the unstable droplet falls between $`T_M(565)`$ and $`T_M(1130)`$. This suggests a crude interpretation of this process, namely that the two individual 565-atom clusters melt when $`T`$ exceeds $`T_M(565)`$ during coalescence, but that the resulting 1130 atom molten droplet is now unstable as $`T<T_M(1130)`$. Finally, it seems reasonable to characterize the liquid state seen immediately after coalescence as metastable: while for small clusters dynamical coexistence between the solid and liquid state has been observed Honeycutt and Andersen (1987), the energy barriers between the liquid and solid states scale as $`N^{2/3}`$ Lynden-Bell and Wales (1994), hence repeated melting and freezing of the 1130-atom cluster will not be observed.

To summarize, we have identified three regimes for the coalescence of free solid metal clusters: 1. solid-solid coalescence by diffusion. This will produce highly-defective clusters, perhaps with large numbers of stacking faults (see Ref. Lewis et al. (1997)). 2. complete melting of solid clusters as they coalesce to form a single liquid droplet. In an IGA source, the final structures will depend on the cooling rate due to equilibration with the inert gas Baletto et al. (2002). 3. complete melting of solid clusters coalescing to form an unstable liquid droplet which subsequently solidifies. The solid structures produced by this mechanism will resemble those produced by a rapid quench.


text:
# Acknowledgments

## Acknowledgments

This work is supported in part by the U.S. Department of Energy under contracts DE-FG02-91ER40677 (RGL).


text:
# Isoscaling in Peripheral Nuclear Collisions around the Fermi Energy and a Signal of Chemical Separation from its Excitation Energy Dependence

## Abstract

The isoscaling is investigated using the fragment yield data from fully reconstructed quasi-projectiles observed in peripheral collisions of 28Si with 124,112Sn at projectile energies 30 and 50 MeV/nucleon. The excitation energy dependence of the isoscaling parameter $`\beta ^{}`$ is observed which is independent of beam energy. For a given quasi-projectile produced in reactions with different targets no isoscaling is observed. The isoscaling thus reflects the level of N/Z-equilibration in reactions with different targets represented by the initial quasi-projectile samples. The excitation energy dependence of the isoscaling parameter $`\beta ^{}`$, corrected for the trivial 1/T temperature dependence, does not follow the trend of the homogeneous system above 4 MeV/nucleon thus possibly signaling the onset of separation into isospin asymmetric dilute and isospin symmetric dense phase.

The isotopic composition of nuclear reaction products provides important information on the reaction dynamics and the possible occurrence of a phase transition in the asymmetric nuclear matter Muller and Serot (1995); Kolomietz et al. (2001), which is supposed to lead to separation into a symmetric dense phase and asymmetric dilute phase. It has been discussed in literature Li and Ko (1997); Baran et al. (2001); Colonna et al. (2002); Marqueron and Chomaz (2003) to what extent such a phase transition is generated by fluctuations of density or concentration, typically suggesting a coupling of different instability modes. The N/Z ( neutron to proton ratio ) degree of freedom and its equilibration was studied experimentally in detailed measurements of the isotopic distributions of emitted fragments Yennello et al. (1994); Johnston et al. (1996); Ramakrishnan et al. (1998); Laforest et al. (1999); Veselsky et al. (2000a); Rami et al. (2000); Xu et al. (2000). Isotopically resolved data in the region Z=2–8 have revealed systematic trends, which however are substantially affected by the decay of the excited primary fragments. It has recently been shown Tsang et al. (2001a) that the effect of sequential decay of primary fragments can be overcome by comparing the yields of fragments from two similar reactions. For statistical fragment production in two reactions with different isospin asymmetry, but at the same temperature, the ratio R21(N,Z) of the yields of a given fragment (N,Z) exhibits an exponential dependence on N and Z. This scaling behavior is termed isoscaling Tsang et al. (2001a) and has been observed in a variety of reactions under the conditions of statistical emission and equal temperature Tsang et al. (2001b, c); Botvina et al. (2002); Souliotis et al. (2003).

In this study, we present an isoscaling analysis of a recent fragment yield data Laforest et al. (1999); Veselsky et al. (2000a, b, 2001) obtained by charged-particle calorimetry of hot thermally equilibrated quasi-projectiles from the reactions 28Si+124,112Sn at 30 and 50 MeV/nucleon. The angular coverage and granularity of the forward-angle multidetector array FAUST Gimeno-Nogues et al. (1997) allowed full reconstruction of the hot quasiprojectiles from mid-peripheral binary collisions Veselsky et al. (2000a). Excellent description of the fragment observables was obtained using the model of deep-inelastic transfer ( DIT ) Tassan-Got and Stefan (1991) for the early stage of collisions and the statistical multifragmentation model ( SMM ) Bondorf et al. (1995) for de-excitation. The large N/Z range of quasi-projectiles allowed unique studies of the dependence of fragment observables on quasi-projectile N/Z. An increasingly inhomogeneous N/Z-distribution between light charged particles and heavier fragments was observed with increasing isospin asymmetry Veselsky et al. (2000b). An exponential scaling of the isobaric ratio Y(3H)/Y(3He) with the quasi-projectile N/Z was observed and the temperature was extracted from the logarithmic slope Veselsky et al. (2001). In the present work we examine the isoscaling phenomenon in unique conditions of highly-selective data. Typically, the investigations of isoscaling focused on yields of light fragments with Z=2-8 originating from de-excitation of massive hot systems produced using reactions of mass symmetric projectile and target at intermediate energies Tsang et al. (2001a, b, c) or by reactions of high-energy light particle with massive target nucleus Botvina et al. (2002). In a recent article Souliotis et al. (2003), we report the investigations of isoscaling using the heavy residue data from the reactions of 25 MeV/nucleon 86Kr projectiles with 124Sn,112Sn and 64Ni, 58Ni targets. Here we present the investigation of the isoscaling phenomenon on the full sample of fragments emitted by the hot thermally equilibrated quasi-projectiles with mass A=20-30.

A thermally equilibrated system undergoing statistical decay can be, within grand-canonical approach, characterized by a yield of fragments with neutron and proton numbers N and Z Albergo et al. (1985); Randrup and Koonin (1981):

$$Y(N,Z)=F(N,Z)\mathrm{exp}[B(N,Z)/T]\mathrm{exp}(N\mu _n/T+Z\mu _p/T)$$
(1)
where $`F(N,Z)`$ represents contribution due to the secondary decay from particle stable and unstable states to the ground state; $`\mu _n`$ and $`\mu _p`$ are the free neutron and proton chemical potentials; $`B(N,Z)`$ is the ground state binding energy of the fragment, and $`T`$ is the temperature.

The ratio of the isotope yields from two different systems, having similar excitation energies and similar masses, but differing only in N/Z, cancels out the effect of secondary decay and provides information about the excited primary fragments Tsang et al. (2001a). Within the grand-canonical approximation ( Eq. 1 ), the ratio $`Y_2(N,Z)/Y_1(N,Z)`$ assumes the form

$$R_{21}(N,Z)=Y_2(N,Z)/Y_1(N,Z)=C\mathrm{exp}(\alpha N+\beta Z)$$
(2)
with $`\alpha `$ = $`\mathrm{\Delta }\mu _n`$/T and $`\beta `$ = $`\mathrm{\Delta }\mu _p`$/T, with $`\mathrm{\Delta }\mu _n`$ and $`\mathrm{\Delta }\mu _p`$ being the differences in the free neutron and proton chemical potentials of the fragmenting systems. C is an overall normalization constant. Alternatively, as introduced by Botvina et al. Botvina et al. (2002), the isoscaling can be expressed as

$$R_{21}(N,Z)=Y_2(N,Z)/Y_1(N,Z)=C\mathrm{exp}(\alpha ^{}A+\beta ^{}(NZ))$$
(3)
thus introducing the scaling parameters which can be related to the isoscalar and isovector components of free nucleon chemical potential since $`\alpha ^{}`$ = $`\mathrm{\Delta }(\mu _n+\mu _p)`$/2T and $`\beta ^{}`$ = $`\mathrm{\Delta }(\mu _n\mu _p)`$/2T. In this variant of isoscaling, the mass dependence is rather weak and the information is concentrated in the isoscaling parameter $`\beta ^{}`$.

In Figs. 1 and 2, we present the isoscaling data from statistical decay of hot quasiprojectiles from the reactions 28Si+124,112Sn at projectile energy 50 and 30 MeV/nucleon, respectively. Observed charged particles with Z$``$5 have been isotopically resolved and the total observed charge was close to the charge of the projectile ( Z=12-15 ) Laforest et al. (1999); Veselsky et al. (2000a). The observed fragment data provide full information ( with the exception of emitted neutrons ) on the de-excitation of thermally equilibrated hot quasi-projectiles with known mass ( A=20-30 ), charge, velocity and excitation energy Veselsky et al. (2000a). In Figs. 1 and 2, the isoscaling plots are presented not only for the full data but also for five bins of excitation energy. The isoscaling occurs at both projectile energies, both for the inclusive data and for the sub-sets with different excitation energies. The slope depends on the excitation energy in a similar way as the slope of the dependence of isobaric ratio Y(3H)/Y(3He) on the quasi-projectile N/Z observed in Veselsky et al. (2001). On the other hand, the isoscaling parameters for individual excitation energy bins do not depend on the projectile energy and the difference of the slopes of the inclusive data is given by the different excitation energy distributions of hot quasi-projectiles at two projectile energies.

In Fig. 3 are presented isoscaling plots for three selected quasiprojectiles 24,26,28Al ( representing the central part of quasiprojectile mass and charge distribution ) in the reactions of 28Si+124,112Sn targets at projectile energy 50 MeV/nucleon. Isoscaling behavior is not observed which means that no memory of the entrance channel interaction leading to production of a given quasi-projectile remains and the isoscaling analysis is not sensitive to the statistical fluctuations in de-excitation. The slight hint of isoscaling in the case of 28Al may be due to low statistics or due to larger uncertainty in the number of unobserved neutrons of the reconstructed neutron-rich quasiprojectile 28Al caused by more intense neutron emission than for proton-rich isotopes 24,26Al. A behavior similar in appearance to Fig. 3 was reported by Tsang et al. Tsang et al. (2001b) using the yields of fragments with Z=3-8 from the reaction of 16O with 197Au, 208Pb targets at projectile energy 20 MeV/nucleon Gelbke et al. (1978), measured at 40. There, the absence of isoscaling behavior suggests that the sample of original prefragments is independent of the target isospin asymmetry, what for such inclusive data suggests a geometric reaction scenario such as breakup or abrasion and the observed fragments are most probably the spectator parts of the original projectile, ejected to a large laboratory angle by the recoil in the projectile breakup.

As suggested by Botvina et al. Botvina et al. (2002), the isoscaling parameter trends are similar to the observed double isotope ratio temperature. The temperature for the current data was studied in detail in Veselsky et al. (2001) and we use the results of that work for comparison with the isoscaling parameters. In Fig. 4 we present the observable $`\beta ^{}T`$, canceling out the trivial 1/T-dependence of the isoscaling parameter. In Fig. 4a we use the double isotope ratio temperature Albergo et al. (1985) from the isotopic ratios Y(2H)/Y(3H) and Y(3He)/Y(4He) and in Fig. 4b we use the temperature obtained from the slope of the dependence of Y(3H)/Y(3He) on the quasi-projectile N/Z ( mirror nucleus temperature ) Veselsky et al. (2001). Grossly, the $`\beta ^{}T`$ should reflect the level of N/Z equilibration Souliotis et al. (2003), which is known for the current data Veselsky et al. (2000a) where full N/Z-equilibration was not reached. The horizontal lines represent the estimate of the zero temperature values of $`\beta ^{}T`$ of reconstructed quasi-projectiles using the known N/Z equilibration Veselsky et al. (2000a). We estimate the zero temperature proton and neutron chemical potentials by the proton and neutron separation energies ( $`\mu _{p,n}S_{p,n}`$ ). Similar approximation was used to extract the mirror nucleus temperature Veselsky et al. (2001). The zero temperature estimates of $`\beta ^{}T`$ for all excitation energy bins are consistent with horizontal lines which is a result of the selection of quasiprojectiles with Z=12-15. Any experimental deviations from the zero temperature values of $`\beta ^{}T`$ can be understood as a non-trivial dependence representing the details of de-excitation. The thick and thin lines represent the $`\beta ^{}T`$ estimates with and without the correction for neutron emission from backtracing of DIT/SMM simulations Veselsky et al. (2000a). The dashed/full lines represent the projectile energies 30 and 50 MeV/nucleon, respectively. The open circles/squares represent the experimental data at 30 and 50 MeV/nucleon while the full squares represent the combined data from both energies. The experimental dependences of $`\beta ^{}T`$ on excitation energies indicate an initial decrease at low excitation energies ( consistent with expansion of the homogeneous excited source ), further they exhibit a turning-point at 4 MeV/nucleon followed by increase toward the zero temperature value at 6-7 MeV/nucleon. The trend is similar at both projectile energies independent of the method of temperature extraction and quality of the statistical sample. The existence of such a turning-point can be understood in the context of work Veselsky et al. (2001), where the viability of the approximation $`\mu _{p,n}S_{p,n}`$ at high temperatures was supported by a possible counterbalance of two effects, namely expansion of hot nucleus ( leading to a decrease of the absolute values of chemical potentials compared to the zero temperature values as indicated by a solid curve ) and chemical separation into an isospin symmetric heavy fraction and an isospin asymmetric nucleon gas, leading to stronger dependence of $`\mu _n\mu _p`$ of the nucleon gas ( dilute phase ) on N/Z of the system Muller and Serot (1995); Kolomietz et al. (2001) thus resulting in higher sensitivity of the observed yield ratios. Such a conclusion is supported by the inhomogeneous isospin distribution observed for the current data Veselsky et al. (2000b) where the isospin asymmetry of light charged particles exhibits strong sensitivity on N/Z of the quasiprojectile while the isospin asymmetry of heavier fragments is rather insensitive. The estimate of $`\beta ^{}T`$ for homogeneous system was obtained assuming $`\rho ^{2/3}`$-dependence and using the estimate of free volume ( eq. 52 of ref. Bondorf et al. (1995) ) for each excitation energy bin from the successful DIT/SMM simulations. Since the isoscaling parameter depends directly on the free nucleon chemical potentials, the turning-point at 4 MeV/nucleon can be understood as a signal of the onset of chemical separation which reverts the decrease of the free nucleon chemical potential consistent with expansion of the homogeneous system. This is also consistent with the recent theoretical study of thermodynamical properties of finite nuclei of Sanzhur et al. Sanzhur  where the hot nucleus ( at a given initial excitation energy ) can expand up to the turning point where the central part begins to contract again while the increasing number of nucleons form the dilute component surrounding the dense central part. The values of the temperature characterizing the turning point in the central density for a given excitation energy reproduce well the caloric curve Veselsky et al. (2001) for the data investigated here. The existence of the turning-point explains the difference of caloric curves ( around 4 MeV/nucleon ) obtained in Veselsky et al. (2000a) using double isotope ratio and mirror nucleus temperatures, former free of assumptions on chemical potentials ( the same applying to Fig. 4a ) while latter using the approximation $`\mu _{p,n}S_{p,n}`$ ( thus artificially lowering the values of $`\beta ^{}T`$ around the turning-point ). When estimating the density of the homogeneous system at the turning point ( assuming $`\rho ^{2/3}`$-dependence of $`\beta ^{}T`$ ), it is consistent with theoretical estimates of the position and shape of the spinodal contour, typically at almost constant total density 0.6$`\rho _0`$ in the wide range of asymmetries Marqueron and Chomaz (2003). Thus, inside the spinodal region the homogeneous nuclear medium is quickly replaced by the inhomogeneous system.

In summary, the isoscaling is investigated using the fragment yield data from fully reconstructed quasi-projectiles observed in peripheral collisions of 28Si with 124,112Sn at projectile energies 30 and 50 MeV/nucleon. The excitation energy dependence of the isoscaling parameter $`\beta ^{}`$ is observed which is independent of the beam energy. For a given reconstructed quasi-projectile produced in reactions with different targets no isoscaling is observed. Thus, the magnitude of the isoscaling is determined by the difference of initial quasi-projectile samples produced in reactions with different targets leading to different level of N/Z-equilibration. The excitation energy dependence of the observable $`\beta ^{}T`$, corrected for the 1/T temperature dependence, exhibits a turning-point at 4 MeV/nucleon which can be interpreted as a signal of the onset of separation into isospin asymmetric dilute and isospin symmetric dense phase. The initial trend of expansion of the homogeneous system is reverted by the transition into an inhomogeneous system with two phases having different proton and neutron concentrations. A novel experimental information on the equation of state of the asymmetric nuclear matter is thus obtained.

We are thankful to L. Tassan-Got for his DIT code and to A. Botvina for his SMM code. This work was supported in part by the Robert A. Welch Foundation through grant No. A-1266, the Department of Energy through grant No. DE-FG03-93ER40773, and the Slovak Scientific Grant Agency through grant VEGA-2/1132/21.


text:
# High Magnetic Field Phase Diagram of PrOs4Sb12

## Abstract

The magnetic phase diagram of PrOs4Sb12 has been investigated by specific heat measurements between 8 and 32 T. A new Schottky anomaly due to excitations between two lowest crystalline-electric-field (CEF) singlets, has been found for both $`H(100)`$ and $`H(110)`$ above the field where the field-induced ordered phase (FIOP) is suppressed. The constructed $`HT`$ phase diagram shows weak magnetic anisotropy and implies a crossing of the two CEF levels at about 8 - 9 T for both field directions. These results provide an unambiguous evidence for the $`\mathrm{\Gamma }_1`$ singlet being the CEF ground state and suggest the level crossing (involving lowest CEF levels) as the driving mechanism of FIOP.

A recent discovery of superconductivityBauer  in PrOs4Sb12 is an important milestone in the field of strongly correlated electrons for a number reasons. Firstly, a large discontinuity in the specific heat at the superconducting transition temperature $`T_c`$ unambiguously establishes that Pr-based materials can support heavy fermion (hf) states. Secondly, proposed crystalline electric field (CEF) schemes suggest non-magnetic ground states and exclude a conventional Kondo effect, believed to be the source of hf behavior in Ce- and U-based metals. The only other serious microscopic mechanism considered as the origin of the hf behavior is a controversial quadrupolar Kondo effect. Cox  Finally, two superconducting transitions observedVollmer  in this compound imply an anisotropic order parameter. Therefore, PrOs4Sb12 is a candidate for a novel mechanism of superconducting pairing that is due to neither phonons nor magnetic interactions.

Correct interpretations of these new exciting phenomena hinge crucially on the proper accounting for the CEF configuration of Pr. The consensus is that the non-magnetic (non-Kramers) $`\mathrm{\Gamma }_3`$ doublet is the ground state, and $`\mathrm{\Gamma }_5`$ is the first excited triplet. However, available experimental data do not rule out the possibility of the $`\mathrm{\Gamma }_1`$ singlet as the ground state. Among the strongest arguments for both configurations were those provided by specific heat in zero and small magnetic fields. In particular, the zero field Schottky-like anomaly at 3.1 K can be related to the $`\mathrm{\Gamma }_3`$ \- $`\mathrm{\Gamma }_5`$ model, assuming these two levels are split by 7.5 K, or $`\mathrm{\Gamma }_1`$ \- $`\mathrm{\Gamma }_5`$ model with the splitting of 8.4 K. Inelastic neutron scattering measurementsMaple  yield the lowest-energy peak at 8 K. Entropy changes in small fields and as a function of temperature were also used to argue for either of the modelsBauer ; Vollmer ; Maple ; Aoki . The difficulty with interpreting these low temperature, low field results is related to strong hybridization of 4f and conduction electrons, as evident from the large electronic specific heat coefficient and the size of the discontinuity in the specific heat ($`C`$) at $`T_c`$. Therefore, we have performed specific heat measurements in high magnetic fields, in which a coupling between conduction and f electrons is expected to be suppressed, to get insight into the CEF configuration of Pr in PrOs4Sb12.

The crystal used in our investigation was from the same batch of high quality samples for which a host of other measurements were reported.Aoki ; dHvA  The specific heat study to 14 T was performed in a superconducting magnet. The measurements above 14 T were carried out in a resistive Bitter magnet. The calorimetry employed the relaxation method, in which a small sample of about 1.5 mg was used to optimize measuring conditions. Figs. 1 and 2 show specific heat in fields ranging from 10 to 32 T. The field was applied along the nominal (100) direction. Our estimate of the probable error in the sample alignment is about $`5^o`$. All the data presented are after subtracting the phonon background approximated by a Debye temperature ($`\mathrm{\Theta }_D`$) of 165 K (after Vollmer et al.)Vollmer ). This value of $`\mathrm{\Theta }_D`$ is quite controversial. In fact, Bauer et al.Bauer , Aoki et al. Aoki , and Maple et al.Maple  suggested $`\mathrm{\Theta }_D`$ to be 304, 320, and 259 K, respectively. Changing $`\mathrm{\Theta }_D`$ to 320 K does not alter the main conclusions of our study, but has some quantitative impact on the results, as discussed further. A relatively high low-temperature limit was chosen to avoid complications associated with a nuclear contribution of Pr. This nuclear contribution is strongly enhanced by coupling with orbital moments of f-electrons.nuclear  The specific heat, at temperatures where the nuclear degrees of freedom dominate, is difficult to measure by a conventional relaxation method because of an additional time scale(s) entering the experiment, the nuclear spin-lattice relaxation time $`T_1`$. AT  Strongly non-exponential temperature decays observed at the lowest temperatures (e.g., below 0.5 K in the field of 10 T and below 1.5 K in the field of 32 T) indicate the importance of nuclear degrees of freedom and cannot be analyzed using the so called $`\tau _2`$ correction. Therefore, these lowest temperature points carry large uncertainty and are not analyzed in detail.

Figure 1 shows the low temperature specific heat for 10 and 12 T. Pronounced peaks in the vicinity of 1 K correspond to a transition to a field induced ordered phase (FIOP), which was first observed by Aoki et al.Aoki  for fields larger than 4 T also applied along the (100) direction. The highest field used in this latter study was 8 T. This long range order was confirmed by other investigations, including the specific heat of Vollmer et al.Vollmer  and magnetization study of Tayama et al.Tayama  It has been proposed that the order parameter corresponding to FIOP is antiferro-quadrupolar type. This type of order is consistent with large anomalies in the specific heat and a very small value of the ordered (antiferromagnetic) moment (about 0.025 $`\mu `$B at 0.25 K)Kohgi . Combining our results with those of Aoki et al. Aoki  and Vollmer et al. Vollmer  shows that the transition temperature ($`T_n`$) identified as a position of the peak in the specific heat ($`C`$) reaches a maximum value somewhere near 9 T. Similarly, $`C(T_n)`$ has a non-monotonic field dependence, attaining the highest value between 8 and 10 T. Notice a large reduction of $`C(T_n)`$ and small but clearly discernible decrease of $`T_n`$ between 10 and 12 T. It is interesting to note a relatively large discrepancy for $`C`$ at $`T_n`$ between all three reports for overlapping fields. This discrepancy is probably due to the anisotropic response of the specific heat in magnetic fields near $`T_n`$, discussed below, and some misalignment of the sample with respect to the (100) direction. An increase of the field from 12 to 13 T (Fig. 2) results in a small reduction of $`T_n`$ and large suppression of $`C(T_n)`$. Moreover, a shoulder appears on the high temperature side of the FIOP anomaly. The specific heat value at this shoulder is about 3400 mJ/K mol. This shoulder evolves into a broad maximum for a slightly higher field of 13.5 T. Again, the value of $`C`$ at the maximum is between 3400 and 3500 mJ/K mol. In this field we no longer detect a signature corresponding to FIOP. Thus, our results strongly imply a disappearance of FIOP before $`T_n`$ reaches 0. A similarly broad peak exists at all fields studied up to at least 32 T. The temperature of this broad maximum increases with the field strength (Fig. 2).

The magnitude of the maximum, in fields of 13 T and larger, is essentially field independent and ranges from 3300 to 3500 mJ/K mol. These values are within about 10% of the maximum value for a Schottky anomaly of a two level system with identical degeneracies. Gopal  The uncertainty of our specific heat measurements in these fields (and at temperatures where nuclear contribution is small) is about 10%. In addition, there is the aforementioned uncertainty associated with a subtraction of the phonon term. Increasing the Debye temperature from 165 K, used in our subtraction, to the other extremal value proposed, 320 K, raises the estimate of the electronic part of $`C`$ by about 290 mJ/K mol at 3.5 K. Thus, the extracted values at the maximum are well within the realistic error bar of the theoretical 3650 mJ/K mol. Since the 32 T field is large enough to significantly split any degenerate levels, the observed Schottky anomaly has to represent excitations between two singlets. The temperature of the maximum, $`T_m`$, is related to the energy separation of the two levels $`\delta `$, $`T_m=0.417\delta `$. An extrapolation of $`T_m`$ to $`T=0`$ (Fig. 3) determines the field at which the two levels cross. Our best estimate of this crossing field, using the lowest fields, is between 8 and 9 T.

These results can be used to infer new information pertaining to a plausible crystal field configuration of Pr. Pr can be modeled by the following single-site HamiltonianKohgi .

$$=_{\mathrm{CEF}}g_J\mu _B\mathrm{𝐉𝐇}𝒥𝐉^{}𝐉\underset{i}{}Q_iO_i^{}O_i,$$
(1)
where $`_{\mathrm{CEF}}`$, $`𝐉`$ and $`O_i`$ represent the CEF Hamiltonian for the cubic $`T_h`$ symmetry, the total angular momentum and the $`i`$-th quadrupole moment (see Ref.Kohgi  for detailed definition). Note that the conventionally used $`_{\mathrm{CEF}}`$Bauer  with the $`O_h`$ symmetryLea  is only approximate because it ignores the the tetrahedral symmetry of Sb ions.Takegahara . This latter CEF Hamiltonian describes the Zeeman effect roughly but cannot account for the neutron diffraction results in the FIOP regionKohgi . In particular, the lower $`T_h`$ symmetry is needed to explain the appearance of a small antiferromagnetic moment along the (010) direction, and not along (100), when the field is applied parallel to the (001) direction. (Directions (100) and (010) are not equivalent). Using the CEF parameters proposed by KohgiKohgi  for the $`\mathrm{\Gamma }_1`$ \- $`\mathrm{\Gamma }_5`$ CEF configuration, $`T_m`$ (with $`Q_i=0`$) and the $`O_{yz}`$-type quadrupolar ordering temperature $`T_n`$ are calculatedcomment . Note that this model consistently reproduces the inelastic neutron dataMaple , characteristic maxima appearing in the magnetic susceptibilityBauer ; Aoki , and the zero-field specific heatBauer ; Vollmer ; Aoki . The results of our calculations are presented in the inset to Fig. 3.

The measured phase diagram (Fig. 3) and the theoretical one expected for the $`\mathrm{\Gamma }_1`$ \- $`\mathrm{\Gamma }_5`$ model (inset to Fig. 3) are strikingly similar. In particular, in both diagrams, the crossing field is very close to the one at which the transition temperature of FIOP becomes maximum. Some overestimation of $`T_n`$ can be attributed to the mean-field nature of the calculations. On the other hand, there is less defined correlation between these two fields in the alternative $`\mathrm{\Gamma }_3`$ \- $`\mathrm{\Gamma }_5`$ model. Instead, it has been argued for this latter model that FIOP is related to the crossing of the upper $`\mathrm{\Gamma }_3`$ and lower $`\mathrm{\Gamma }_5`$ levels.Vollmer  Thus, the observed correlation between the two characteristic fields constitutes our first argument for the $`\mathrm{\Gamma }_1`$ singlet being the lowest CEF level.

More conclusive arguments regarding the CEF configuration can be obtained from the study of the anisotropy of the Zeeman effect. The Zeeman effect is expected to be strongly anisotropic for the $`\mathrm{\Gamma }_3`$ doublet and almost isotropic for the $`\mathrm{\Gamma }_1`$ singlet,Tayama  for the range of CEF parameters that are consistent with the susceptibilityBauer . In particular, there is no crossing involving the lowest CEF levels when the field is applied along the (110) direction in the $`\mathrm{\Gamma }_3`$ \- $`\mathrm{\Gamma }_5`$ model. We emphasize that this prediction is quite general; it does not depend on whether one uses the the CEF Hamiltonians and parametrs proposed by Bauer et al.Bauer  or Kohgi et al.Kohgi . Figure 4 shows the measured high field phase diagram for this field direction. For this direction we observe a decrease of $`T_n`$ values with respect to the (100) direction for the corresponding fields, consistent with the previous magnetization measurementsTayama . On the other hand, within the uncertainty of the measurement, there is no change in the position of the Schottky anomaly at 13 and 14 T, as expected for the $`\mathrm{\Gamma }_1`$ CEF ground state and inconsistent with the $`\mathrm{\Gamma }_3`$ scenario. Moreover, for the (110) orientation we clearly observe the Schottky anomaly already at 12 T. This lower field limit for the Schottky maximum is probably due to competition between the two types of anomalies, discussed below, and lower values of $`T_n`$ for the (110) direction. A straight line fit for the three $`T_m`$ points results in the crossing field value of 8.5 $`\pm `$ 0.5 T. This value agrees with our estimate for the (100) direction.

Existence of the crossing field for the (110) direction provides an unambiguous evidence for the $`\mathrm{\Gamma }_1`$ \- $`\mathrm{\Gamma }_5`$ model, in fact the strongest evidence to date. We stress that a small misalignment of the sample with respect to the field in either of the measurements cannot explain essentially identical crossing fields for both directions. In fact, the measured difference in $`T_n`$ values for (100) and (110) directions provides an additional check of the alignment. Similarly to the (100) direction, there seems to be a close correlation between the crossing field and the field corresponding to $`T_n`$ maximum.

Figures 3 and 4 imply a strong competition between the field-induced order and Schottky peak. The FIOP transition in the specific heat abruptly disappears before $`T_n`$ reaches zero. Precise magnetizationTayama  and de Haas van Alphen measurementsdHvA , on the other hand, were able to map $`T_n`$ as a function of the magnetic field all the way to $`T_n0`$. This apparent contradiction can be explained by a very small entropy available for the FIOP transition above 13 and 12 T for fields parallel to the (100) and (110) directions, respectively. Specific heat, being a bulk measurement, can be less sensitive than magnetization techniques in this situation. A strong competition is to be expected in the $`\mathrm{\Gamma }_1`$ \- $`\mathrm{\Gamma }_5`$ scenario. The ground state pseudo-doublet formed at the level crossing carries both magnetic and quadrupolar moments. Since a quadrupolar moment operator does not commute with a dipolar one, the quadrupolar interactions leading to FIOP compete with the magnetic Zeeman effect.

Summarizing, we find an evidence for the level crossing occurring around 8-9 T. This field value is very close to the field at which both $`T_n`$ and $`C(T_n)`$ become maximum, suggesting that the crossing of the lowest CEF levels is the driving mechanism of FIOP. In fields above 10 T there is a strong competition between FIOP and the Schottky peak. Our high field specific heat results provide compelling arguments for the $`\mathrm{\Gamma }_1`$ singlet being the CEF ground state of Pr in this system. This results call for a new model of electronic mass enhancement, instead of the ordinary quadrupolar Kondo effect, which has been discussed mainly in the context of the $`\mathrm{\Gamma }_3`$ non-Kramers doublet ground state.

This work has been supported by the U.S. Department of Energy, Grant No. DE-FG02-99ER45748, National Science Foundation, DMR-0104240, the National High Magnetic Field Laboratory, and the Grant-in-Aid for Scientific Research from MEXT of Japan. We thank G.R. Stewart, P. Kumar, K. Ingersent, and R. Shiina for stimulating discussions.


text:
# A Hybrid Monte Carlo Method for Surface Growth Simulations

## Abstract

We introduce an algorithm for treating growth on surfaces which combines important features of continuum methods (such as the level-set method) and Kinetic Monte Carlo (KMC) simulations. We treat the motion of adatoms in continuum theory, but attach them to islands one atom at a time. The technique is borrowed from the Dielectric Breakdown Model. Our method allows us to give a realistic account of fluctuations in island shape, which is lacking in deterministic continuum treatments and which is an important physical effect. Our method should be most important for problems close to equilibrium where KMC becomes impractically slow.

Epitaxial growth on surfaces is of central importance both for applications and as a very interesting example of statistical processes out of equilibrium. We may idealize the process as the introduction of new atoms (adatoms) onto a crystal surface with flux, $`F`$; the adatoms then diffuse, with diffusion coefficient, $`D`$, nucleate islands or attach to existing islands. During the early stages of growth, the submonolayer case, the island size and shape distribution is a matter of substantial practical and theoretical interest.

The island growth process is commonly modeled by kinetic Monte Carlo (KMC) or continuum models. In KMC, internal noise processes are automatically represented within the model and each adatom is represented individually. Therefore, when there are many adatoms (e.g. close to equilbrium) such simulations slow down considerably. A deterministic continuum model which represents the adatoms as a continuous fluid does not have this problem, and should be much faster. There has been considerable work in the development of such models for epitaxial growth (see Gyure et al. (1998); Chen et al. (2001) and references therein). In some cases they have been quite successful, but sometimes they do not reproduce experimental results. One reason for such problems is that deterministic continuum models neglect important fluctuations. In this paper we present a method of dealing with some fluctuations without giving up the advantages of a continuum treatment. We call this approach Hybrid Monte Carlo (HMC). The most important use of this method will be in cases where fluctuations are important, but which would be difficult to treat with KMC because of the presence of a large number of adatoms. A computation in a similar spirit has been given in Schulze and E (2001); Schulze et al. (2003). However, the present approach differs in a number of important ways.

There are several sources of fluctuations in the growth process. The one we consider here is the fact that when atoms attach to islands they do so one at a time – that is, there is shot noise in the island growth process. This is important because island growth limited by diffusion is intrinsically unstable. The surface of the island will grow fingers due to the analogue of the well- known Mullins-Sekerka instability of metallurgy Mullins and Sekerka (1963). In the context of thin film growth the instability was discussed in detail by Bales and Zangwill Bales and Zangwill (1990). The reason for unstable growth is easy to see: if a finger on the edge of an island starts to grow it will project out onto the terrace and be fed by more adatoms than the portions behind. The finger will grow longer, and be fed by still more flux, etc. Edge diffusion and other restructuring processes smooth out the fingers, and the final shape depends on the competition between unstable growth and smoothing.

An extreme case of the diffusively unstable growth is represented by the Diffusion-Limited Aggregation (DLA) model of Witten and Sander Witten and Sander (1981). In this model all smoothing processes are neglected and growth takes place so slowly that one random walking adatom at a time is considered. DLA clusters are sprawling fractal objects with many branches which resemble some cases of island growth Hwang et al. (1991); Tsui et al. (1996). There is a variation of the model called the Dielectric Breakdown Model (DBM) Niemeyer et al. (1984) in which random walkers are not used. Instead the Laplace equation is solved outside the aggregate for a field, $`\rho `$, which represents the probability density of walkers, and the growth algorithm is to add one particle at a time with probability proportional to $`\rho /n`$ at the surface. Thus the DLA limit can be successfully treated by a model in which the adatoms are a continuum. However, methods such as the level-set method Chen et al. (2001) cannot go to this limit and thus cannot produce dendritic islands which are seen in experiment.

The unstable modes for the interface of the islands are present in level-set models, of course. However, the reason why the DBM limit is not achieved is that the growth process is represented by deterministically advancing a continuum interface according to the flux of the adatom fluid into the surface (see Eq. (4) below). In this algorithm the amplitude of the perturbations to a smooth interface are not correctly represented: they are given either by the initial conditions or by computer roundoff errors. In the experiment, however, there is a mechanism for feeding the instability: each adatom attaches not as a spread-out advance of the interface, but as an atom. In the case of DLA the result is that noise in the shape is present at all scales Somfai et al. (1999) and does not average out.

We should mention that there are other fluctuations which we will not treat. For example, there are density fluctuations in the adatom fluid which are important in the nucleation of new islands. There is a method Ratsch et al. (2002) to treat this within level-set theory that we could employ. We will not consider such processes in this work, but rather look at the shape fluctuations of existing islands.

The HMC algorithm goes as follows: we treat the islands on the surface as crystals containing discrete atoms which occupy the sites of a lattice. To illustrate the method we use a square lattice here. On the other hand, the adatoms are treated as a continuum whose surface density is $`\rho `$, and which is governed by:

$$_t\rho =D^2\rho +F$$
(1)
In practice we solve this equation numerically on a discrete square grid which is commensurate with the crystal.

Eq. (1) is solved with periodic boundary conditions on the edge of the system. The important physics of growth in incorporated into the boundary condition at the surface of the island. We put:

$$D\frac{\rho }{n}=k_\pm (\rho \rho _o)$$
(2)
Here $`k_{}`$ and $`k_+`$ are the respective attachment rates on the upper and lower edge of the island boundary.

In the case of irreversible growth only the first term in brackets in Eq. (2) would be present. The other term, $`\rho _o`$, accounts for the detachment of atoms from the island, and depends on the position on the island boundary. Physically, the boundary condition must allow for faster detachment at corners, say, than at flat surfaces. We represent this in a way that allows us to compare directly with the bond-counting version of KMC:

$$\rho _o=\mathrm{exp}(nE/k_BT)$$
(3)
Here $`n`$ is the number of nearest neighbor bonds that must be broken to *completely detach* (see below) the atom in question and add it to the adatom sea. The value of $`n`$ depends on the environment of the detaching atom. We imagine that atoms can break bonds by moving along lattice directions. Also, $`E`$ is the bond energy, and $`T`$ the temperature. Of course, we can easily incorporate bonding to more distant neighbors.

In a pure continuum model, the velocity of the island growth would be determined by mass conservation:

$$v_n=a^2D\left[\frac{\rho }{n}\right]$$
(4)
where $`a`$ is the lattice constant and $`[]`$ denotes the jump across the island boundary. Note that we can interpret Eq. (2) in a way which is familiar in studies of crystal growth Langer (1980). By combining Eqs. (2,4) we find an equation for the density at the surface:

$$\rho \overline{\rho }=(\rho _o\overline{\rho })+\alpha v_n$$
(5)
where $`\alpha =1/(a^2[k_++k_{}])`$, and $`\overline{\rho }`$ is the equilibrium density near a flat surface. That is, we are including both local equilibrium and kinetic terms in our boundary condition. The difference $`\rho _o\overline{\rho }`$ is a measure of the number of dangling bonds on the surface, and thus of the curvature (upon coarse-graining). That is, the first term in Eq. (5) is related to the familiar Gibbs-Thompson boundary condition of crystal growth, and the second is a kinetic term.

In HMC we implement Eq. (4) in a way that *includes fluctuations*. Consider first a case where attachment is the only important process. Then we solve Eq. (1, 2) and compute the total flux onto the island boundary using Eq. (2, 4). When the total flux exceeds one atom then an adatom is attached to the boundary at random with the probability proportional to $`v_n`$ (exactly as in the DBM model.) In the case where detachment is also present we consider the surface to be partitioned into the part where the net flux is inward (growth), and attach atoms with probability density $`v_n`$, and outward (detachment) and remove island atoms with probability density $`v_n`$.

We have implemented HMC and compared the results to a KMC code on the same square lattice using nearest neighbor bonding. We use the hopping rate of an adatom to set the unit of time, and the lattice constant to be unity. We have two independent parameters namely $`D/F`$, and $`ϵ=E/k_BT`$ We can also add edge diffusion, but here we have not done so. The KMC code is written using the method described in Newman and Barkema (1999) which takes advantage of the fact that there are only a few independent jump probabilities for a bond-counting model. We find that, as expected, at large $`ϵ`$, the KMC code is much faster. However, for $`ϵ=1.5`$ and with about 300 adatoms in a 40x40 system the speeds are comparable. We note that in situations such as heteroepitaxy Lam et al. (2002) there are a great many independent probabilities to jump and KMC is slower.

The interpretation of the number of bonds, $`n`$, is a bit delicate. In Fig. (1) we show some examples of what we mean by ‘complete detachment’. For our square lattice, to detach an atom on a $`[10]`$ surface we need to break 3 bonds. However, for a $`[11]`$ surface, in order to detach from the surface, an atom must first break two bonds, and then subsequently, one more – see Fig. (1). Thus, for a $`[11]`$ surface we set $`n=3`$ as well because this corresponds to the product of the probabilities of the two processes. In effect, we have coarse-grained. For other possible surface environments it is not difficult to tabulate the correct value of $`n`$.

This procedure may seem arbitrary, but we can justify it by its results. To this end, we show that HMC gives the correct shape for an island in equilibrium with adatoms. This shape is known exactly Rottman and Wortis (1981) from a mapping to the 2d Ising model. We did simulations with $`F=0`$ starting with a square island, and ran our code for a long enough time that the system seemed to be in equilibrium. In Fig. (2) we show some results superimposed on the exact results for various temperatures. The simulation results are ensemble averages; that is, we did 20 independent simulations and averaged the density of the island after shifting the center of mass of each one to be at the origin. The dotted line in Figure 2 is the contour line where the averaged island density is 1/2.

We now present some more results to demonstrate the technique. If growth dominates detachment we should have DLA-like structures. Whether this occurs depends on $`T`$ and $`D/F`$ Barabasi and Stanley (1995). For low temperature, Fig. (3), we see the transition in a very clear way. Fig. (4) shows a higher $`T`$ case. Of course, all of these effects can be seen in KMC simulations. In other simulations (not shown), we have demonstrated that edge diffusion also smoothes out dendritic shapes, as expected. The virtue of our method will be to treat systems near equilibrium where the dynamics of fluctuations are of interest and where there are a great number of adatoms, as in Fig. (4).

Another example is the thermal broadening of steps due to repeated attachment and detachment of adatoms. There is an extensive literature in this area. Jeong and Williams (1999), and the theoretical expectation Jeong and Williams (1999); Khare and Einstein (1998) is that the thermal width, $`w`$, of a step should depend on the rate-limiting mechanism for step motion. In our model, without surface diffusion, this will be the either detachment from the step or diffusion on the terrace. In these cases $`w^2t^{1/2},t^{1/3}`$, respectively. We show, in Fig. (5), $`w^2(t)`$ for a pair of steps one atomic layer high which cross a $`100\times 100`$ terrace. We see indications of both of the expected behaviors at late times. The early time behavior is *kinetic* roughening Kardar et al. (1986), $`w^2t^{2/3}`$. This is because we started with no adatoms, and, initially, the steps were retreating.

The HMC technique should be most useful in situations where there is a large separation of time scales between the diffusion of the adatoms and fluctuations of island boundaries. Other examples for which it might be used are in studies of homogeneous nucleation of large islands near equilibrium on surfaces Theis and Tromp (1996). Another example is in the catalysis of the reaction CO+O $``$ CO2 on a noble metal. The diffusion of CO is very fast in this case, and makes a realistic KMC simulation very difficult. Evans and collaborators Suchorski et al. (2001) have used an approximation where $`\rho _{CO}`$ is constant where ever there is no O on the surface. Our approach (with suitable changes in boundary conditions to allow for the reaction) should give a better account of the kinetics.

Finally, we should mention the case of heteroepitaxy. There have been a number of KMC simulations of this very important process Lam et al. (2002); Orr et al. (1992); Khor and Sarma (2000). This is a very difficult problem to attack numerically. The time-consuming step in such calculations is not the dynamics of the adatoms, but rather the elasticity computation. Nevertheless, our approach to the adatoms offers a number of advantages, and we intend to pursue this in a future publication.

We acknowledge useful conversations with B. Orr, R. Calflisch, and J. Evans. PS is supported by NSF grant DMS 0207402. GR acknowledges a travel grant from the Michigan Center for Theoretical Physics.