text:
# Phonon Density-of-States in MgB2

## Abstract

We report inelastic neutron scattering measurements of the phonon density–of–states in Mg11B2, which has a superconducting transition at 39.2 K. The acoustic phonons extend in energy to 36 meV, and there are highly dispersive optic branches peaking at 54, 78, 89 and 97 meV. A simple Born-von Kàrmàn model is able to reproduce the mode energies, and provides an estimate of the electron-phonon coupling of $`\lambda 0.9`$. Furthermore, the estimated boron and magnesium contributions to the isotope effect are in qualitative agreement with experiment. The data confirm that a conventional phonon mechanism, with moderately strong electron-phonon coupling, can explain the observed superconductivity.

There has already been rapid progress in our understanding of the properties of MgB2, since the first announcement of superconductivity with a critical temperature of 39K. Although some have speculated that the high value of Tc must require a novel mechanism for superconductivity, initial experimental reports suggest that a conventional phonon mechanism is sufficient. Both the observation of the isotope effect and a simple BCS form for the energy gap measured by point-contact tunneling are consistent with phonon-mediated superconductivity, with the high transition temperature made possible by the high frequency of boron vibrational modes. Initial estimates of phonon frequencies, from electronic structure calculations, predict that they are consistent with the observed critical temperature, although An et al and Kong et al predict that only a few modes have significant electron-phonon coupling.

The strength and frequency-dependence of the electron-phonon coupling is determined by measurements of both the electron-phonon spectral density, $`\alpha ^2(\omega )`$F($`\omega `$), and the bare phonon density–of–states (PDOS), F($`\omega `$). $`\alpha ^2(\omega )`$F($`\omega `$)can be derived by inversion of tunneling data, while F($`\omega `$) is derived from inelastic neutron scattering measurements. In systems where both sets of data are available, $`\alpha (\omega )`$ usuallly varies smoothly with frequency so the bare phonon density-of-states provides useful estimates of the frequency moments that determine $`T_c`$.

We report measurements of the phonon density–of–states on an isotopically enriched polycrystalline sample of Mg11B2 using time-of-flight neutron spectroscopy at the Intense Pulsed Neutron Source, Argonne National Laboratory. They confirm that boron vibrational frequencies extend as high as 100 meV. The acoustic modes peak at 36 meV, while there are dispersive optic branches at 54, 78, 89 and 97 meV, in good agreement with predictions based on electronic structure calculations. A simple Born-von Kàrmàn (BvK) model, which reproduces the energies of the main spectral features, is used to determine the logarithmic frequency moment of the PDOS, which provides a preliminary estimate of the electron-phonon coupling, $`\lambda 0.9`$ (assuming $`\mu ^{}=0.1`$). We also predict the scale of PDOS shifts with isotopic mass, providing estimates of the magnesium and boron contributions to the isotope effect, although the values depend on the definition of the PDOS used. Given the uncertainties in this analysis, the agreement with experiment is satisfactory, and allows us to conclude that the superconductivity is consistent with a conventional phonon mechanism with moderately strong electron-phonon coupling.

The MgB2 sample was synthesized from a high purity, 3mm diameter Mg rod and isotopic 11B (Eagle Picher, 98.46 atomic % 11B). The Mg rod was cut into pieces about 4mm long and mixed with the 200 mesh 11B powder. The reaction was done under moderate pressure (50 bars) of UHP Argon at 850C. At this temperature the gas-solid reaction was complete in about one hour. The sample (approximately 20 g) was contained in a 30 ml tantalum crucible with a closely fitting cover. There was surface reaction between the Ta crucible and the reactants at the synthesis temperature but the products of this reaction could be mechanically removed. X-ray diffraction showed some small, unidentified impurity peaks in the resultant powder, but the impurity concentration is small and should have no effect on the PDOS measurements.

The neutron scattering experiments were performed on the Low Resolution Medium Energy Chopper Spectrometer (LRMECS), using incident energies of 200 meV, 130 meV, and 50 meV. LRMECS has continuous detector coverage with scattering angles from 2.4 to 117.5 so that a broad range of wavevector transfers are measured simultaneously, e.g. 2 to 12 Å-1 at an energy transfer of 60 meV. This is extremely important in accurate determinations of the PDOS when the scattering is predominantly coherent, as it is in Mg11B2. The coherent (incoherent) scattering cross sections are 3.63 barns (0.08 barns) and 5.56 barns (0.21 barns) for Mg and 11B respectively. Unless the scattering is purely incoherent, it is necessary to average the measured neutron scattering over a large volume of reciprocal space in order for the resulting measurements to reflect the true phonon density–of–states.

The neutron data taken with an incident energy of 200 meV shows that the PDOS extends up to 100 meV. Above that energy, the scattering only comprises a broad multiphonon tail. Figure 1 shows data taken at 130 meV at 8K. The peak in the acoustic PDOS occurs at about 36 meV, while there are broad optic bands at 54, 78, 89 and 97 meV. The spectrometer resolution decreases with energy transfer from 5.9 meV at $`\mathrm{}\omega =40`$ meV to 3.4 meV at $`\mathrm{}\omega =90`$ meV. The large width of these features is not just a resolution effect, but arises from dispersion of the optic phonon branches. We have also performed measurements at 30K, 50K, and 100K. The only statistically significant difference between these temperatures is a small suppression of the peak at 54 meV below Tc. In a polycrystalline average, this is a very small effect, but it could still represent a significant change in an individual phonon branch. Interestingly, this is the mode that An and Pickett predict to have the strongest electron-phonon coupling, but such conclusions need to await single crystal phonon dispersion measurements.

In a multicomponent system, the inelastic neutron scattering law in polycrystalline samples is given by

$$S(Q,\omega )=\underset{i=\mathrm{Mg},\mathrm{B}}{}\sigma _i\frac{\mathrm{}Q^2}{2M_i}\mathrm{exp}(2W_i)\frac{\mathrm{G}_i(\omega )}{\omega }[n(\omega )+1]$$
(1)
where $`\sigma _i`$ and $`M_i`$ are the neutron scattering cross section and atomic mass of the $`i^{th}`$ atom and $`n(\omega )=[\mathrm{exp}(\mathrm{}\omega /k_BT)1]^1`$ is the Bose population factor. The generalized PDOS, $`\mathrm{G}(\omega )=_i\mathrm{G}_i(\omega )`$, where Gi($`\omega `$) is defined as

$$\mathrm{G}_i(\omega )=\frac{1}{3N}\underset{j𝐪}{}|𝐞_i(j,𝐪)|^2\delta [\omega \omega (j,𝐪)]$$
(2)
which, in turn, defines the Debye-Waller factor, Wi(Q) through

$$W_i(Q)=\frac{\mathrm{}Q^2}{2M_i}_0^{\mathrm{}}\frac{G_i(\omega )}{\omega }[2n(\omega )+1]$$
(3)
$`\omega (j,𝐪)`$ and $`e_i(j,𝐪)`$ are the frequencies and eigenvectors, respectively, of the phonon modes. Because the Mg and B partial contributions are weighted by $`\sigma /M`$, the measured spectrum does not represent the true PDOS, although the stronger cross section and lighter mass of the boron atoms ensure that they make the dominant contribution.

We have analyzed the data in terms of a simple BvK model, that is sufficent to describe the principal features. MgB2 has a hexagonal structure with the lattice parameters $`a=3.082`$ Åand $`c=3.515`$ Å. The lattice is composed of a honeycomb of B atoms separated by hexagonal Mg layers, leading to comparisons with graphite intercalation compounds. The nearest-neighbor distances are 3.082 and 3.515 Å for Mg-Mg pairs in the basal plane and along the $`c`$-axis, respectively, 2.501 Å for Mg-B pairs, and 1.780 and 3.515 Å for B-B pairs in the basal plane and along $`c`$-axis, respectively.

Assuming central forces, each atom-atom interaction can be represented by two parameters, the longitudidinal and transverse force constants. We have included only nearest-neighbor interactions between Mg-Mg and B-B atoms in the basal plane and along the $`c`$-axis, and between Mg-B pairs in neighboring layers. The results of the calculation, which was performed using the UNISOFT program, are shown in Fig. 1(b), where the in–plane and out–of–plane projections of the partial contributions to G($`\omega `$) from boron and magnesium modes are also plotted.

It is possible to reproduce the frequencies of the peaks in the PDOS with only six force constants. These are the longitudinal and transverse in-plane B-B interactions (90 and 7.5 N/m respectively), the longitudinal out-of-plane B-B interaction (20 N/m), the longitudinal in-plane and out-of-plane Mg-Mg interactions (both 0.6 N/m) and the longitudinal and transverse Mg-B interactions (70 and 2 N/m respectively). The remaining transverse nearest-neighbor interactions had a negligible effect on the mode frequencies.

The frequencies of the main peaks in the PDOS are well-represented by this model, but it is clear that the widths of the peaks are not, i.e. we have underestimated the dispersion of the optic branches. It would be possible to improve the agreement by the introduction of additional parameters representing longer-range interactions, but the results of such an analysis would be meaningless at the present time. Complete measurements of the phonon dispersions on single crystals would allow the model to be refined, but the present analysis is sufficient to provide good estimates of the frequency moments of the PDOS that determine the phonon contribution to Tc.

Fig. 2 shows the resulting bare PDOS, F($`\omega `$) defined as

$$\mathrm{F}(\omega )=\frac{1}{N}\underset{j𝐪}{}\delta [\omega \omega (j,𝐪)]$$
(4)
This is similar to equation 2 but does not include the vibrational eigenvectors. In a Bravais lattice, the two quantities are equivalent, but in multicomponent systems with very different atomic masses, G($`\omega `$) will tend to weight the PDOS to higher frequencies because of the smaller atomic displacements of the heavier atoms.

We now discuss the implications of this data for understanding the superconducting transition temperature. This may be estimated using the Allen–Dynes equation

$$k_BT_c=\frac{\mathrm{}\omega _{ln}}{1.2}\mathrm{exp}\left[\frac{1.04(1+\lambda )}{\lambda \mu ^{}(1+0.62\lambda )}\right]$$
(5)
where $`\lambda `$ is the electron-phonon coupling, $`\mu ^{}`$ is the Coulomb repulsion, and $`\omega _{ln}`$ is the logarithmic phonon average defined in equation 2.16 of ref. . $`\omega _{ln}`$ is an average over $`\alpha ^2(\omega )`$F($`\omega `$), so we need to make the assumption that $`\alpha (\omega )`$ is approximately energy-independent, in which case it cancels from the equation. However, we have a choice of two functions to use. It is more conventional to use the bare PDOS, F($`\omega `$), given by equation (4). However, this function fails to take into account the size of the atomic displacements of the different phonon modes; they are included in G($`\omega `$). The electron-phonon matrix elements contain the product of the mode eigenvectors and the gradient of the crystal potential, so it is likely that G($`\omega `$) represents a better approximation to $`\alpha ^2(\omega )`$F($`\omega `$). In fact, the choice of F($`\omega `$) or G($`\omega `$) only affects significantly the calculation of the isotope effect, which we discuss below.

The Allen-Dynes equation can be inverted to estimate the value of $`\lambda `$ for a given value of Tc and $`\mu ^{}`$. From integrating over the model calculation of G($`\omega `$), we estimate that $`\omega _{ln}`$ is 57.9 meV \[53.8 meV when determined from F($`\omega `$)\]. Fig. 3 shows calculations of Tc using three typical values of $`\mu ^{}`$ compared to the measured value of Mg11B2. This shows that the measured value of Tc is consistent with a phonon mechanism with values of $`\lambda `$ in the range 0.6 to 1.2. $`\lambda =0.9`$ at the most frequently cited value of $`\mu ^{}=0.1`$. $`\omega _{ln}`$ nearly coincides with the energy of the E2g mode that An and Pickett predict will dominate the electron-phonon coupling. It is impossible to determine, on the basis of our measurements, whether the entire phonon spectrum is involved in the superconductivity or only a subset that occurs close to the average energy.

Finally, we can use the BvK model to estimate the isotope effect. The dotted line in Fig. 2 shows a calculation of the PDOS in Mg10B2, illustrating the significant energy shifts in all the high frequency boron modes. Substituting 26Mg for 24Mg produces much smaller shifts. In the absence of strong coupling effects, the total isotope effect is given by the sum of the partial isotope effects, $`\beta =_{i=Mg,B}\beta _i`$, where $`\beta _i=\delta \mathrm{ln}\mathrm{T}_\mathrm{c}/\delta \mathrm{ln}M_i=\delta \mathrm{ln}\omega _{ln}/\delta \mathrm{ln}M_i`$ and should always equal 0.5 if $`\mu ^{}=0`$. More accurate predictions suggest that it can fall below this value when $`\lambda `$ becomes small, but not at the value of 0.9 that we estimate. If we use G($`\omega `$) to estimate $`\beta _B`$ and $`\beta _{Mg}`$, we obtain values of 0.42 and 0.08 respectively, whereas using F($`\omega `$), we obtain values of 0.28 and 0.22, respectively. Reported value of $`\beta _B`$ are in the range 0.28(3) to 0.30(2). Hinks et al have measured an average value of $`\beta _{Mg}=0.018(8)`$. The experimental values are slightly lower than our calculations but the difference may not be significant given the limitations of our analysis.

In conclusion, we have measured the phonon density–of–states in Mg11B2, and used the data to produce a simple BvK force-constant model of the vibrational spectra. Our analysis allows us to conclude that the measured values of the transition temperature and the isotope effect are consistent with a conventional phonon mechanism for the superconductivity with moderately strong electron-phonon coupling. It will be extremely useful to compare our data with tunneling measurements of the electron-phonon spectral density, which determine whether the whole phonon spectrum is involved in the superconductivity, or only a few specific modes.

###### Acknowledgements.

This work was supported by the U.S. DOE Office of Science under contract no. W-31-109-ENG-38.


format: text

text:
# Universality of Velocity Gradients in Forced Burgers Turbulence

## Abstract

It is demonstrated that Burgers turbulence subject to large-scale white-noise-in-time random forcing has a universal power-law tail with exponent -7/2 in the probability density function of negative velocity gradients, as predicted by E, Khanin, Mazel and Sinai (1997, Phys. Rev. Lett. 78, 1904). A particle and shock tracking numerical method gives about five decades of scaling. Using a Lagrangian approach, the -7/2 law is related to the shape of the unstable manifold associated to the global minimizer.

The universality of small-scale properties in fully developed Navier–Stokes (NS) turbulence is frequently investigated assuming that a steady state is maintained by a large-scale random force. For structure functions (moments of increments) universality with respect to the force is conjectured in the case of three-dimensional NS turbulence and proven for certain linear passive scalar models (see, e.g., Ref. rmp ). The universality of probability density functions (p.d.f.) for velocity increments and gradients is a difficult question which, so far, has been mostly addressed within the framework of the pressureless model of Burgers turbulence, usually the one-dimensional Burgers equation

$$_tu+u_xu=\nu _{xx}u+f(x,t)$$
(1)
with white-noise-in-time forcing cy95 . It is generally conjectured that, when $`\nu 0`$ and the forcing is confined to large scales, the tail of the p.d.f. of velocity gradients $`\xi `$ at large negative values follows a universal power-law $`p(\xi )|\xi |^\alpha `$. The actual value of the exponent is however a matter of controversy. Let us briefly recall some of the arguments found in the literature.

A standard approach is based on studying the inviscid limit of the Fokker–Planck equation for the p.d.f.

$$_tp_\xi \left(\xi ^2p\right)\xi p+\nu _\xi \left(_{xx}\xi |\xi p\right)=B_{\xi \xi }p,$$
(2)
where the right-hand side expresses the diffusion of probability due to the delta-correlation in time of the forcing. It was pointed out by Polyakov p95  that the inviscid limit of (2) contains anomalies due to the singular behavior of the dissipative term $`\nu _\xi \left(_{xx}\xi |\xi p\right)`$. The value $`\alpha =3`$ is obtained if anomalies are ignored gk98  or if a piecewise linear approximation is made for the solutions of the Burgers equation bm96 . An operator product expansion (OPE) method borrowed from quantum field theory has been proposed for evaluating such anomalies and an argument presented in favor of $`\alpha =5/2`$ (actually, for velocity increments and infinite systems) p95 . However, this expansion leads to a relation involving unknown coefficients which must be determined, e.g., from numerical simulations yc96b98 , and restricts the possible values to $`5/2\alpha 3`$ b01 . Anomalies cannot be understood without a complete description of the singularities of the solutions, such as shocks, and of their statistical properties. For the case of a space-periodic system (as we shall assume), a crucial observation made in Ref. ekms97  is that large negative gradients stem mainly from preshocks, that is the cubic-root singularities in the velocity preceding the formation of shocks ff83 . A simple argument was given in Ref. ekms97  for determining the fraction of space-time where the velocity gradient is less than some large negative value. This leads to $`\alpha =7/2`$ provided preshocks do not cluster. Determinations of the dissipative anomaly of (2) have been made by formal matched asymptotics eve99  and by bounded variation calculus eve00 . With the assumption that shocks are born with vanishing amplitude from isolated preshocks, the value $`\alpha =7/2`$ was obtained eve99 ; eve00 . Other attempts to derive $`\alpha =7/2`$ using also isolated preshocks have been made sevenhalf . Note that there are simpler instances, including time-periodic forcing bfk00  and decaying Burgers turbulence with smooth random initial conditions eve00 ; bf00 , which fall in the universality class $`\alpha =7/2`$, as can be shown by systematic asymptotic expansions using a Lagrangian approach. In the presence of forcing, the key issues which remained to be settled are the possible clustering of preshocks and, closely related to this, the possible birth of shocks with non-vanishing amplitude. The results presented hereafter almost completely rule out such possibilities.

Numerically solving the randomly forced Burgers equation in the limit of vanishing viscosity in such a way as to obtain clean scaling for the p.d.f. of gradients represents a significant challenge. Broadly speaking, there are two classes of methods. On the one hand, methods involving a small viscosity, either introduced explicitly (e.g. in a spectral calculation) or stemming from discretization (e.g. in a finite difference calculation). Viscosity gives rise to a power-law range with exponent $`1`$ at very large negative gradients gk98  whose presence will make the inviscid $`|\xi |^\alpha `$ range appear shallower than it actually is, unless extremely high spatial resolution is used. On the other hand, there are methods which directly capture the inviscid limit with the appropriate shock conditions such as the fast Legendre transform method of Ref. nv94  (extended to the forced case in Ref. bfk00 ). This method is very well adapted to decaying Burgers turbulence with non-smooth Brownian-type initial data vdun94  but, with spatially smooth forcing, it leads to delicate interpolation problems which have been overcome in the case of time-periodic forcing bfk00 ; with white-noise-in-time forcing, it is difficult to prevent spurious accumulations of preshocks leading to $`\alpha =3`$. To avoid such pitfalls, we develop a Lagrangian particle and shock tracking method notebennaim  which is able to cleanly separate smooth parts of the solution and is particularly effective for identifying preshocks. The main idea of the method is to consider the evolution of a set of $`N`$ massless point particles accelerated by a discrete-in-time approximation of the forcing with a uniform time step. When two of these particles intersect, they merge and create a new type of particle, a shock, characterized by its velocity (half sum of the right and left velocities of merging particles) and its amplitude. The particle-like shocks then evolve as ordinary particles, capture further intersecting particles and may merge with other shocks. In order not to run out of particles too quickly, the initial small region where particles have the least chance of being subsequently captured is determined by localization of the global minimizer (see below). The calculation is then restarted from $`t=0`$ for the same realization of forcing but with a vastly increased number of particles in that region. This method gives complete control over shocks and preshocks notefilm  and allows an accurate determination of the relevant statistical quantity while keeping a manageable number of degrees of freedom.

Fig. 1 shows the p.d.f. of the velocity gradients in log-log coordinates at negative values for a Gaussian forcing restricted to the first three Fourier modes with equal variances such that the large-scale turnover time is order unity. Quantitative information about the value of the exponent is obtained by measuring the “local scaling exponent”, i.e. the logarithmic derivative of the p.d.f. calculated here using least-square fits on half-decades. It is seen that over about five decades, the local exponent is within less than 1% of the value $`7/2`$ predicted by E et al. ekms97 . This value of the exponent was also obtained numerically (with a fewer particles) for other large-scale forcing instances with compactly supported or exponentially decreasing spectra and also for non-Gaussian forcing (e.g. with Fourier amplitudes having a Bernoulli distribution or an uniform distribution in an interval). Evidence for non-clustering of preshocks is obtained by counting the average number of shock formations per unit time. For different types of large-scale forcing, we found that the typical mean number of preshocks per turnover time is comparable to the number of forcing Fourier modes significantly excited. For such forcings, the density of preshocks is found to vary by not more than 6% when the time step varies by two orders of magnitude around $`\delta t=10^4`$, which is hardly consistent with a power-law (and even a logarithmic) divergence as $`\delta t0`$. Furthermore, we have checked that shocks are always born with vanishing amplitude (within numerical errors).

Turning now to theoretical results, let us briefly recall the construction of solutions developed by E et al. ekms00 , in terms of the dynamical system associated to the characteristics of (1) in the inviscid limit note . The force is assumed to derive from a Gaussian potential $`F(x,t)`$, delta-correlated in time, periodic of period 1 and analytic in space. A statistically stationary régime is reached by taking the initial time at $`\mathrm{}`$. The central point of the construction is the following variational characterization of the solution at an arbitrary time ($`t=0`$ chosen for convenience):

$$u(x,0)=\frac{}{x}\underset{X()}{\mathrm{min}}\left[_{\mathrm{}}^0\left[\frac{1}{2}\dot{X}^2(t)F(X(t),t)\right]𝑑t\right],$$
(3)
where the minimum is taken over all piecewise smooth (absolutely continuous) curves $`X(t)`$ with $`t(\mathrm{},0]`$ such that $`X(0)=x`$. A curve minimizing the action in (3) is called a minimizer and should be understood as a fluid particle trajectory. It obviously has to satisfy for all $`t<0`$ the Euler–Lagrange equations:

$`\dot{X}(t)`$ $`=`$ $`U(t),`$ (4)
$`\dot{U}(t)`$ $`=`$ $`f(X(t),t).`$ (5)
Except for a finite number of $`x`$-values, there exists a unique minimizer ekms00 . The locations where there are more than one minimizer correspond to shocks. The minimizers converge exponentially fast backward in time to the trajectory of the unique fluid particle which is never absorbed by a shock. This trajectory is called the global minimizer because its action is minimal at any time; it corresponds to a hyperbolic trajectory of the dynamical system (4)-(5). Associated to it, there are two curves in the phase-space $`(x,u)`$: a stable (attracting) manifold $`\mathrm{\Gamma }^{(\mathrm{s})}`$ and an unstable (repulsive) manifold $`\mathrm{\Gamma }^{(\mathrm{u})}`$. The minimizers converge backward in time to the global minimizer and, thus, the graph of the solution is made of pieces of the unstable manifold with jumps at shocks. One of these shocks, called the main shock, is singled out. It is the unique shock which has always existed in the past (whereas generic shocks are born at some finite time $`t<0`$); it may be shown that it corresponds to the position giving rise to the left-most and the right-most minimizers which approach the global one backward in time. The other shocks cut through the doublefold loops of the unstable manifold (see Fig. 2). We observe that their locations can be obtained by a Maxwell rule applied to those loops. Indeed, the difference of the two areas defined by cutting such a loop at some position $`x`$ is equal to the difference of actions of the two minimizers defined by the upper and lower branches and, thus, vanishes at the shock location.

We also observe that the structure just outlined has much in common with that appearing in the unforced Burgers equation. Indeed when $`f=0`$, the solution to the Burgers equation can be constructed from the Lagrangian manifold in the $`(x,u)`$ plane, defined as the position and the velocity of fluid particles when ignoring shocks. This manifold is parameterized by the Lagrangian coordinate $`a`$; denoting $`u_0`$ the initial velocity, we then have simply $`x=a+tu_0(a)`$ and $`u=u_0(a)`$. The actual solution with shocks is obtained by applying the standard Maxwell rule to the Lagrangian manifold. In the forced case, a parameterization of the unstable manifold (e.g. by the arclength) is now the analog of the Lagrangian coordinate. But there are two important differences: first, in the unforced case, the time evolution of the Lagrangian manifold is explicit and linear while, in the presence of a force, the Euler–Lagrange equations (4)-(5) are not, in general, explicitly solvable and the unstable manifold has a hyperbolic dynamic. Second, the smoothness of the Lagrangian manifold in the unforced case stems directly from the smoothness of the initial data, whereas in the forced case Pesin’s theory must be used to show that when the force is indefinitely differentiable in space, so is the unstable manifold ekms00 .

Using the smoothness of the unstable manifold, we now formally derive the $`7/2`$ law, by an argument mostly borrowed from the unforced case bf00 . Let $`\mathrm{\Gamma }^{(\mathrm{u})}=\{(X(s),U(s))\}`$ with $`s`$ real, be a parameterization of the unstable manifold at time $`t=0`$. It is assumed for convenience that $`s=0`$ corresponds to the global minimizer and that $`X^{}(0)>0`$, where primes denote $`s`$-derivatives. The velocity is exactly obtained by eliminating from the unstable manifold the shaded areas determined by the Maxwell rules and the parts beyond the main shock (shown as dashed lines in Fig. 2). The surviving set of parameter values (excluding shocks) is denoted $`\mathrm{\Omega }`$. Turning to the statistical description, the p.d.f. of velocity gradients may be written

$$p(\xi )=\delta \left(_xu(x,0)\xi \right).$$
(6)
Because of homogeneity, we can integrate over the space period and then change from the $`x`$ variable to the $`s`$ variable to obtain

$$p(\xi )=_0^1p(\xi )𝑑x=_\mathrm{\Omega }\delta \left(\frac{U^{}}{X^{}}\xi \right)X^{}𝑑s.$$
(7)
Note that since a finite gradient is assumed, $`x`$ cannot be at a shock position. Denoting by $`s_k`$ the parameter values where the argument of the delta function vanishes, we obtain

$$p(\xi )=\underset{k}{}\frac{X^2}{|U^{\prime \prime }\xi X^{\prime \prime }|}\delta (ss_k)ds.$$
(8)
For very large negative values of $`\xi `$, the $`s_k`$’s must be near some $`s_j`$, corresponding to a local minimum of $`X^{}`$. Taylor expansions of $`X`$ and $`U`$ in the vicinities of the $`s_j`$’s and the use of the Maxwell rule show that the $`s_j`$’s are located in space-time near preshocks satisfying $`X^{}=0`$ and $`X^{\prime \prime }=0`$ with $`X^{\prime \prime \prime }>0`$ (see Fig. 2). Proceeding as in Ref. bf00 , we finally obtain, to leading order

$`p(\xi )C|\xi |^{7/2},\xi \mathrm{},`$ (9)
$`C={\displaystyle \frac{5\sqrt{2}}{2}}{\displaystyle _{\mathrm{\Omega },X^{\prime \prime \prime }>0}}|X^{\prime \prime \prime }|^{1/2}|U^{}|^{5/2}\delta (X^{})\delta (X^{\prime \prime })𝑑s.`$ (10)
Hence, the constant involves the mean of $`|X^{\prime \prime \prime }|^{1/2}|U^{}|^{5/2}`$ at preshocks. Its evaluation requires the knowledge of the joint probability distribution of $`X^{}`$, $`X^{\prime \prime }`$, $`X^{\prime \prime \prime }`$ and $`U^{}`$. From the Euler–Lagrange equations (4)-(5), we observe that a set of ordinary differential equations with nonlinear stochastic forcing is easily obtained for $`X`$, $`U`$ and the aforementioned four variables. From these equations, using techniques similar to those developed in Ref. ekms00  (where a subset of these stochastic equations is studied), it should be possible, on the one hand, to make our derivation more rigorous (including for the non-clustering of preshocks) and, on the other hand, to obtain an upper bound for the constant $`C`$ in the $`7/2`$ law. Note that the expression for $`C`$ involves also an integral over the admissible set of parameters $`\mathrm{\Omega }`$ whose determination cannot in general be done by local analysis with ordinary differential equations. This is why only an upper bound is expected.

As noted in Ref. ekms97 , the universality with respect to the forcing of the p.d.f. of large negative velocity gradients may be extended to negative velocity increments, provided that they are not significantly influenced by shocks. Without understanding of all the mechanisms leading to small-amplitude shocks in the forced case, the issue of universality for the p.d.f.’s of velocity increments cannot be settled. A first step would be to determine numerically the distribution of shock amplitudes. Note that our technique may also be extended to the case of forcing at scales much smaller than the size of the system, a problem close to that considered by Polyakov p95 , which is left for future work.

I wish to express my gratitude to U. Frisch and K. Khanin for numerous helpful discussions and encouragements. I also thank M. Blank, G. Eyink, A. Fouxon, R. Mohayaee, E. Vanden Eijnden, M. Vergassola and V. Yakhot for useful remarks. Work in Nice was supported by the European Union, under contract HPRN-CT-2000-00162; work at CNLS was supported by the U.S. Department of Energy, under contract W-7405-ENG-36. Simulations were performed on Avalon at the Los Alamos National Laboratory.


text:
# Testing large mixing MSW solutions of the solar neutrino problem through Earth regeneration effects

## 1 Introduction

The available neutrino data from the Homestake , GALLEX-GNO , SAGE , Kamiokande  and Super-Kamiokande  experiments indicate a significant $`\nu _e`$ flux deficit with respect to the predictions of the Standard Solar Model (SSM) . The $`\nu `$ oscillation mechanism represents a possible explanation of this deficit. In particular, recent analyses  show a preference of the data for the so-called MSW  (Mikheyev-Smirnov-Wolfenstein) solutions involving large mixing, usually indicated as “large mixing angle” (LMA) and “low $`\delta m^2`$” (LOW) solutions. We point out how present and future experiments can help to discriminate such two solutions through Earth regeneration effects.

## 2 Energy dependence of the Earth regeneration effect in SK and SNO

The slightly positive indication ($`1.3\sigma `$) for an excess of nighttime to daytime events in Super-Kamiokande , if confirmed with higher statistical significance, would indicate the occurence of the Earth regeneration effect for $`{}_{}{}^{8}B`$ neutrinos. Such indication, by itself, might not be sufficient to discriminate the LMA from the LOW solution, since a slight excess is predicted in both cases. On the other hand, the Earth regeneration effect depends strongly on the neutrino energy and, in principle, one could take advantage of this feature to discriminate the two solutions. In particular, starting from the simple observation that the Earth regeneration effect is stronger at low energy for the LOW solution, and at high energy for the LMA solution, we point out that it may be useful to study the night-day asymmetry in two separate energy ranges in both the Super-Kamiokande and the SNO experiments . For definiteness, we consider the two following representative ranges for the total (measured) energy of recoiling electrons in SK and SNO,

$$\mathrm{Low}\mathrm{range}(\mathrm{L})=[5.0,7.5]\mathrm{MeV},$$
(1)
$$\mathrm{High}\mathrm{range}(\mathrm{H})=[7.5,20]\mathrm{MeV},$$
(2)
and calculate the night-day rate asymmetry in such ranges,

$$A_{H,L}=\left(\frac{ND}{N+D}\right)_{H,L}.$$
(3)
Since one expects $`A_H>A_L`$ for the LMA solution and $`A_H<A_L`$ for the LOW solution, it is useful to introduce the difference

$$\mathrm{\Delta }=A_HA_L,$$
(4)
which should change sign when passing from the LMA region ($`\mathrm{\Delta }>0`$) to the LOW region ($`\mathrm{\Delta }<0`$).

Figures 1 and 2 show the result of our calculations of $`\mathrm{\Delta }`$ (eccentricity effects removed) in SK and SNO, respectively, in the form of isolines at $`\mathrm{\Delta }\times 100=\pm 0.5`$, $`\pm 1`$ and $`\pm 2`$, superimposed to the current MSW solutions at $`90\%`$, $`95\%`$ and $`99\%`$ C.L., obtained using the latest experimental data. Such figures confirm that $`\mathrm{\Delta }`$ sign effectively changes when passing from LMA to LOW solution, thus representing a useful test to discriminate the two large mixing angle solutions.

As one can notice in the figures 1 and 2, it turns out that, in the most favorable circumstances, the value of $`\mathrm{\Delta }\times 100`$ can reach $`0.5`$ (SK) and $`1.0`$ (SNO) in the case of the LOW solution and 1.0 (SK) and 2.0 (SNO) for the LMA. This small effect could be observable in the two experiments, provided that they can reach a sensitivity to $`\mathrm{\Delta }`$ at a (sub)percent level. However, the simple preference of the SK and SNO experiments for a definite $`\mathrm{\Delta }`$ sign would already constitute a useful information, corroborating other tests envisaged to solve the LOW-LMA ambiguity .

Our choice of the low and high energy ranges in Eqs. (1, 2) is only representative, and it should be tuned to reach the best compromise between sensitivity to the energy dependence of Earth regeneration effect and statistical significance. For instance, one could also enhance the sensitivity to the energy dependence by calculating the day-night asymmetry in a larger number of bins

(D-N asymmetry spectrum); however, the statistical significance would decrease in each bin. The optimal compromise requires dedicated experimental studies.

## 3 Testing LOW solution through time variations in BOREXINO and GNO

The Earth regeneration effect occurs at relatively low values of $`\delta m^2`$ for low energy (7Be and pp) neutrinos, that can be detected by the two experiments BOREXINO (based on $`\nu _e`$ scattering) and GNO (based on $`\nu _e`$-Ga absorption). The first is a real-time experiment, able to detect a possible regeneration phenomenon as a day-night asymmetry of the event rate. The second is a radiochemical experiment with a signal extraction timescale of $`1`$ month. It can, however, detect Earth effect in the form of a seasonal modulation of the signal . Indeed, during winter, the nighttime is longer compared to summer and, furthermore, the trajectory of neutrinos probes inner layers of the Earth (in particular, the Earth core is crossed only during winter ).

In Fig. 3 we show the results of our calculation of the day-night asymmetry ($`A_{DN}\times 100`$) of the event rate for BOREXINO. In the same plot we show the current solutions of the solar neutrino problem at $`90\%`$, $`95\%`$ and $`99\%`$ C.L. (for 2 d.o.f.) respectively. This figure shows that the day-night asymmetry expected in BOREXINO in the region of the LOW solution is expected to be in the range $`1\%20\%`$.

Figure 4 shows the iso-lines of Winter-Summer rate difference $`R_{WS}`$ (in SNU) expected in GNO due to the Earth regeneration effect (eccentricity effect removed). Winter and Summer periods are defined as

$$\mathrm{Winter}[23\mathrm{september},21\mathrm{march}],$$
(5)
$$\mathrm{Summer}[22\mathrm{march},22\mathrm{september}].$$
(6)
Notice that the sensitivity in the LOW region is almost maximal since GNO is predominantly sensitive to low-energy pp $`\nu `$’s, as compared with BOREXINO, which is more sensitive to the 7Be line $`\nu `$’s. This circumstance makes these two experiments intrinsically different in testing the LOW solution. The two tests are also complementary in tracking the origin of the time variations: a possible seasonal signal in GNO, if originated by vacuum (instead of MSW) oscillations, should not produce a N-D asymmetry in BOREXINO.


text:
# Spontaneous magnetization of aluminum nanowires deposited on the NaCl(100) surface

## I Introduction

Electronic properties of nanowires adsorbed on solid surfaces depend strongly on the geometric structure of the wire. The surface enables the realization of a specific geometry but the interaction between the substrate and the adsorbed wire does not necessarily destroy the key properties the wire would have in the same geometry but as isolated. This idea has been the basis of many recent experimental and theoretical studies of nanowires on semiconducting  and on metallic  substrates. The studies of nanowires on surfaces have been motivated from the new technology point of view. Nanowires can be seen as ultimate leads when miniaturizing devices in electronics. Moreover, if they exhibit magnetization they could be used as compact magnetic storage devices. In general, the exciting quantum phenomena taking place as the size and the dimensionality of the system are reduced are the driving force in the research of nanostructures.

The manufacturing of nanowires on solid surfaces can be based on the manipulation of the surface and the adsorbate atoms with the scanning tunneling microscope (STM) . A more efficient method for fabrication may be the deposition of adsorbate atoms in the steps on metal surfaces  or the growing of stripe-like adsorbate structures either on metal  or on semiconductor  surfaces. The geometry of these structures can be monitored using STM , He scattering  or by field ion microscope . An efficient way to characterize the electronic structures of wires, e.g. their one-dimensionality, is the angle-resolved photoemission .

In this work we report state-of-the-art electronic structure calculations based on the density-functional theory (DFT) for simple metal nanowires on an insulating substrate, i.e. for Al wires on the NaCl(100) surface. The present work is largely inspired by the theoretical investigations of Ga adsorbate wires on Si(100) surface performed by Watanabe et al. . Using first-principles electronic structure methods they calculated the stable structures of Ga - dangling-bond wires on hydrogen-terminated Si(100). The resulting bandstructures show that the wires are metallic and that in some of the structures the adsorbate-induced band near the Fermi level is very flat. On this basis Watanabe et al.  discussed, using a simple tight-binding model, the possibility of the ferromagnetism of these simple metal wires. Recently, Okada and Oshiyama  returned to these Ga wires on Si(100) performing first-principles DFT calculations. They found one atomic configuration for the Ga wire exhibiting ferromagnetic ordering. We have chosen to work with Al wires because the jellium model calculations  have predicted, and the subsequent first-principles calculations  confirmed, that isolated linear Al atom chains can be magnetic for certain interatomic distances. NaCl has been selected as the substrate due to its insulating character, which is expected to be reflected as minor substrate effects in the electron structure of the Al wire, and due to the fact that the atomic distances on the NaCl(100) surface match with the interatomic distance of a magnetic isolated Al-chain. This surface is also expected to be stable, as a surface reconstruction occurs only in one of many different cleavages . An Al wire on the unreconstructed NaCl surface is a much simpler system than Ga wires on the hydrogenated reconstructed Si(100). Therefore the present calculations are expected to show clearly the essential physical phenomena for simple metal wires on insulating surfaces. Because the properties of a nanowire on a surface are largely derived from those of unsupported wires, we begin this work by calculations for isolated wires in different geometries and study their properties as functions of structural parameters.

The organization of the present article is as follows: In Sec. II we discuss the computational details. We use the Vienna ab initio simulation package (VASP)  which is a pseudo-potential-plane-wave (PPPW) code. Moreover, we have made benchmark calculations using the the full-potential linearized augmented plane wave (FLAPW) method with the WIEN97 package . FLAPW is considered more accurate but it is also much more cpu-time-consuming than the PPPW method. In Section III we present our results, first for the unsupported linear (III A) and dimerized or zig-zag wires (III B) and thereafter for wires on the NaCl(100) surface (III C). Section IV summarizes our results, and includes an outlook.

## II Computational details

We employ in our DFT electronic structure calculations the generalized gradient approximation (GGA) for the exchange-correlation effects . When describing interactions between adsorbate atoms and a substrate it has been found important go beyond the local density approximation (LDA). We do both spin-polarized and spin-compensated computations in order to find the solution of the lowest energy.

In the PPPW calculations (VASP package) we have used ultrasoft Vanderbilt pseudopotentials . In the case of unsupported Al-wires the cutoff energy of the plane-wave expansion has been 180 eV. This is sufficient to describe bulk Al. In the case of the NaCl surface we need a cutoff energy of 250 eV.

We use the supercell approximation in our calculations. This means that instead of isolated atomic chains a set of parallel chains are used. In the case of linear chains the interchain distance of 15 Å is sufficient to eliminate the interaction between the wire and its periodic replicas. The zig-zag chains are treated in a supercell with the lattice constant of 20 Å perpendicular to the chain axis. In the NaCl(100) surface calculations we have used the slab geometry with one or two layers of Na and Cl atoms. The use of one layer is found sufficient to describe the energetics of straight Al wires on the surface. The NaCl slabs are separated by a vacuum region of 20 Å. On the surface, the Al wires are separated by 3 lattice constants of NaCl. These separations were found to depress sufficiently the interactions between the slabs or the Al wires. Dipole corrections for the electrostatics of the supercell geometry are also included .

The desired periodicity along the wire and the artificial periodicity due to the supercell approximation require a Brillouin zone integration. In order to describe correctly the magnetic properties the k-point sampling in the wire directions has to be dense enough. We have used 26 k-points for this purpose. In the case of unsupported wires there are k-points only on the $`k_z`$ axis in the wire direction. When calculating Al wires on the NaCl(100) surface the k-points form a 3x26 mesh on the $`(k_x,k_z)`$ plane parallel to the surface and there are actually 26 k-points in the irreducible wedge of the Brillouin zone. The Fermi-smearing of 0.06 eV is used in the integrations.

In order to test the accuracy of the PPPW calculations, especially with respect to the predicted magnetic properties, we have performed calculations also with the FLAPW method (WIEN97 code) for the unsupported linear Al atom chain. The FLAPW method can be considered as one of the most accurate band structure methods. In the calculations, the muffin-tin sphere radii for Al is 2.10 a.u. Inside the muffin-tin spheres the maximum angular momentum $`l`$ in the radial expansion is $`l_{max}=10`$, and the largest $`l`$ value for the non-spherical part of the Hamiltonian matrix is $`l_{max,ns}=4`$. The cutoff parameters are $`RK_{max}=9`$ for the plane waves, the number of plane waves ranges up to around 4000. The Brillouin zone integrations are done using 10 special points , and a Fermi broadening of 0.001 Ry is used.

## III Results

### A Unsupported linear chain of Al atoms

First we want to test the accuracy of the pseudo-potential approach by comparing the PPPW and FLAPW results in the simple case of an unsupported linear chain of Al atoms. The total energies $`E_{tot}`$ per atom are shown in the uppermost panel of Fig. 1 as a function of the bond length. The energy minima correspond to spin-compensated solutions. According to the FLAPW method the energy minimum occurs at $`d_{nn}=`$ 2.43 Å which is in a fair agreement with the PPPW result of 2.41 Å. The equilibrium force constants $`k=\frac{d^2E_{tot}}{da^2}`$ are 104 N/m (FLAPW) and 151 N/m (PPPW). Our GGA-PPPW values are quite close to the LDA results obtained by Mehrez and Ciraci  using a PPPW method. According to Fig. 1 the difference between the FLAPW and PPPW results increases towards to the larger interatomic distances. This reflects some difficulties in the transferability of the pseudopotential for larger distances and low dimensionality, because the increase of the energy cutoff does not change the PPPW results substantially.

The stabilized-jellium model calculations have predicted that nanowires of simple metals may be magnetic . In Fig. 1 the results obtained using the spin-compensated formalism are connected with solid lines whereas solutions with magnetic moments are connected with dashed lines. At distances larger than 2.7 Å (FLAPW) or 2.6 Å (PPPW) the total energy of a magnetic solution is indeed lower than that of the spin-compensated one. The magnetic solution is due to the spin polarization of the second subband accommodating atomic 3$`p`$ electrons (see Figure 2). In the middle panel of Fig. 1 the resulting magnetic moment per atom is given as a function of the elongation. The magnetic moments calculated by both methods saturate to a value around 0.7 $`\mu _B`$ per Al atom. It is interesting to compare this value with the stabilized-jellium model prediction which gives, corresponding to the interatomic distance of 2.86 Å (bond length in bulk Al), a moment of 0.65 $`\mu _B`$ per Al atom .

The elongation force is defined as the force opposing the lengthening of the wire due to an external applied force. When negative, the elongation force would like to shorten the wire, whereas a positive elongation force without a counterbalancing force would lead to a spontaneous elongation. In the case of a linear atomic chain the elongation force is calculated as

$$F=\frac{dE_{tot}}{dr},$$
(1)
where $`r`$ is the interatomic distance. Naturally, the force vanishes at the minimum energy and is negative for larger distances. For very large distances the force should approach zero. The crucial point is that when the wire is stretched, the force reaches a minimum which determines the maximum force sustained before breaking. In both calculations, the breaking force is about 1.5 nN. This is of the same magnitude as the typical force before breaking seen in the atomic force microscope experiments  on gold nanowires. According to Fig. 1 the breaking of the linear Al chain occurs as a function of the elongation before the magnetization sets on. If the chain would be stable, the onset of the magnetization should result in a discontinuous decrease of the magnitude of the elongation force .

It is instructive to compare the one-dimensional bandstructure of the linear Al chain with the properly folded subband structure of an infinitely long stabilized jellium cylinder. This is done in Fig. 2. The Al chain is calculated using FLAPW for the interatomic distance of 2.86 Å. The splitting of the bands due to the magnetic ground state occurs at this distance. The PPPW method gives quantitatively similar results. In the stabilized-jellium calculation the density $`n`$ of the positive jellium charge corresponds to the electron density parameter $`r_s=(3/4\pi n)^{1/3}`$ = 2.07 a.u. The radius, 1.35 Å, of the positive background charge is chosen to give the linear charge density of 3 $`e`$/ 2.86 Å. The electronic structure has been solved self-consistently within the LDA. The above value of the jellium radius falls in a narrow window, in which the the second, doubly-degenerate subband, corresponding to the quantum numbers $`m`$ = 1 and $`n`$ = 1, is totally spin polarized resulting in the magnetic solution. In the atomic chain calculations the corresponding, nearly totally polarized band arises from the atomic $`p_{xy}`$ orbitals (the $`z`$ axis is parallel to the chain). In a quasi-one-dimensional system the density of states (DOS) diverges at the bottom of a subband. When the Fermi level is just above the bottom of a subband, the Fermi-level DOS will be large and the system will fulfill the Stoner criterion for magnetism . The quantum numbers of the lowest jellium subband are $`m`$ = 0 and $`n`$ = 1. The corresponding atomic chain states have the character of the atomic 3$`s`$ orbital on the lowest branch and after the band is folded at $`k_z=\pi /a`$ the character is that of the atomic 3$`p_z`$ orbitals. The contribution of the lowest subband to the magnetic moment of the wire is small in comparison with the second (nearly) totally polarized subband. The agreement between the jellium and atomic chain models is surprisingly good. This gives credence for the stabilized jellium model for Al even in the present very confined geometry.

In conclusion, in the case of the unsupported linear chain of Al atoms the PPPW calculations are able to reproduce the essential cohesive properties obtained in the all-electron FLAPW calculation. Also the magnetic properties, which are of a especial interest in the preset work, are similar in both calculations. Therefore in the following calculations for the unsupported dimerized and zig-zag Al wires as well as for the Al wires on the NaCl (100) surface we rely on the less computer intensive PPPW method.

#### 1 Unsupported dimerized and zig-zag Al wires

The next step in our study is to allow the atoms of an unsupported linear Al chain to dimerize and also to form zig-zag wires. These geometries are studied because a linear chain of atoms is susceptible to Peierls-like distortions, opening up a small band gap and lowering thereby the total energy. The notations of the geometries studied are presented in Fig. 3. In the dimerized geometry the distance $`2r`$ between every second atom remains constant, but the atom in the middle of a cell moves in the direction of $`r_{short}`$ along the wire. In zig-zag wires the nearest neighbor distance $`r`$ is kept constant while the zig-zag angle $`\theta `$ is varied. When studying the distortions of the spin-polarized chain the distance $`r`$ is 2.8 Å whereas for the spin-compensated chain the linear equilibrium value of 2.41 Å is used.

According to Table I, the energy of the Al wire increases with dimerisation, i.e., the bond length alternation for the Al chain vanishes or is at least below our error estimate of around 2 % of the interatomic distance. We have actually confirmed this energy increase also with the FLAPW calculations. The energy increases more rapidly close to the equilibrium distance corresponding to the non-magnetic wire than close to the larger distance corresponding to the magnetic wire. It is logical to expect that changes in the geometry of an wire already stretched have a smaller effect on the energy than the changes close to the equilibrium. One can also note in the case of spin-polarized wires that the magnetic moment grows slightly along dimerization.

Next we consider zig-zag wires depicted in Fig. 3. Table II gives the total energies and magnetic moments obtained. The zig-zag formation is energetically favorable at the distance $`r`$ = 2.41 Å corresponding to spin-compensated solutions as well as at $`r`$ = 2.80 Å giving magnetic solutions. The values of the equilibrium zig-zag angles $`\theta `$ are around 150o and 130o for the spin-compensated and the magnetic systems, respectively. The deformation is limited by the second-nearest-neighbor atom interactions in the atomic chain. In the magnetic chain, the longer bond length $`r`$ allows for a angle deformation larger than in the spin-compensated chain.

According to Table II the energies of the of the chains decrease when the zig-zag angle decrease from 180o. Especially for the magnetized chain the decrease is strong and it is accompanied by a decrease in the magnetic moment. The minimum values are obtained at the angle of 120o. To understand this behavior an analysis of the band structure required. We do this analysis here, in the context of the unsupported Al chains, because this analysis will be the basis to understand the magnetism of the more complicated systems of Al wires on the NaCl(100) surface. The band structure of the magnetic zig-zag wire at the equilibrium zig-zag angle of 130o is plotted in Figure 4 a). It is dramatically different from that for the linear chain in Fig. 2. As the zig-zag angle $`\theta `$ decreases, the width of the s-band decreases and the band gap at around -5 eV increases. The states of the 3$`p`$-bands that cause the magnetization undergo also significant changes. There is no longer a doubly-degenerate $`p_xp_y`$ band. Two $`sp^2`$-like bands are formed, the lower of which is occupied both by spin-up and spin-up electrons whereas the upper $`sp^2`$-like band is well above the Fermi level. A totally polarized $`p_z`$-band causes the magnetization. As one third of the Al 3$`p`$-electrons occupy the $`p_z`$ band, a magnetic moment of roughly 0.3 $`\mu _B`$ can be expected from this interpretation. This value is indeed in agreement with the actual computed result in Table II. At the zig-zag angle of 120o (see Fig. 4 b) ) the lower $`sp^2`$-like band becomes polarized at low k-values. This explains the increase of the magnetic moment at small zig-zag angles.

### B Al wire on the NaCl(100) surface

In this subsection we present our results for an Al wire in equilibrium on the NaCl(100) surface and at the end of of the subsection we analyze these results using simple model calculations. The NaCl(100) slabs used in the calculations correspond to the atomic layers cut from the perfect bulk crystals with the theoretical lattice constant of 5.68 Å. The omission of the relaxations of Na and Cl ions is justified because we are interested mainly in the effect the surface has on the wire and not vice versa. Before dealing with Al wires on the surface we have made calculations for single isolated Al atoms adsorbed at different sites on NaCl(100). The sites on the top of Na and Cl ions, on the middle of the bridge between the two ions as well as on the hollow site in the center of four ions have been considered. These calculations show that the properties of the adsorbed Al atoms (adsorption energy and the distance from the surface) do not change remarkably when three NaCl(100) layers are used instead of two or the surface supercell is increased from 3x3 to 4x4. The following results correspond to two NaCl(100) layers and the 4x4 cell or in the case of the Al wire, the interwire distance is three NaCl lattice constants.

The results for the isolated Al atom on NaCl(100) serve also as guidelines for later calculations dealing with the Al chains on the NaCl(100) surface. The optimized adsorption heights from the surface-layer atoms and the corresponding total energies are given in Table III. The isolated Al atom prefers to be on the top of a Na atom. The total energies are therefore given relative to this system. Thereafter the next favorable site is the hollow one followed by the bridge site. The adsorption of Al on the top of the Cl atom is clearly the most unfavorable case.

The results for the isolated Al on NaCl(100) in mind we start the simulations for the Al wire commensurate on the NaCl(100) surface. In these calculations there are two rigid layers of Na and Cl ions and the supercell contains two Al atoms. The equilibrium positions of Al atoms obtained reflect the Al-atom - surface interaction as well as the Al-Al interactions within the wire. It turns out that wires in which the Al atoms sit on the top of nearest Na atoms along the  direction as straight or zig-zag chains have not the lowest energies. Also, when we start (inspired by the results for unsupported wires) the structure optimization from an initial configuration, in which the Al atoms sit close to the adjacent hollow sites forming a slightly zig-zaged chain a strong structural relaxation takes place. The resulting configuration shown in Fig. 5 a) is a zig-zag wire in which the Al atoms reside in the bridges between the Na and Cl atoms. The Al atoms are closer to the Cl atoms than the Na atoms so that the zig-zag angle $`\theta =143^o`$ and the interatomic distance in the chain is $`r`$ = 2.99 Å. The distance to the surface is $`d`$= 2.90 Å. It is exciting that the Al wire in this equilibrium configuration has a magnetic moment of 0.24 $`\mu _b`$ per atom! This is just of the order of magnitude expected from Table II for unsupported zig-zag wires. Thus the interaction with the NaCl surface does not destroy the magnetic moment of the quasi one-dimensional Al wire.

The above results underline the importance of relaxing the positions of the wire atoms on the surface. This kind relaxation should be performed also when studying the early low-coverage stages of the growth of metal layers on insulating surfaces . The approximation that lateral interactions between adjacent adsorbate atoms are unimportant is maybe not valid one for all growth patterns.

In order to understand the properties of the above equilibrium wire we have made further calculations using simplified models for linear Al wires on the NaCl(100) surface. In these calculations only one rigid layer of NaCl has been used but the distance between the adjacent wires due to the supercell structure is the same as before. The Al atoms are allowed to relax only perpendicular to the surface. The geometries considered are presented and labeled in Fig. 5. An Al atom has in these geometries one (NN1), two (NN2) or four (NN4) surface ions as nearest neighbors corresponding to atop, bridge and hollow site wires. The supercells in these simulations contain 12 non-equivalent surface ions and two wire atoms.

The energies of the model wires calculated per one Al atom are shown in the lower panel of Fig. 6 as a function of the distance from the surface layer. The NN2 wire has the lowest equilibrium energy and this is used to define the energy zero. The equilibrium distance is nearly the same as for the fully-relaxed zig-zag wire discussed above. The NN4 wire has a clearly higher energy. This explains the behaviour of our above-discussed simulation starting from Al atoms near the hollow sites. The second lowest energy minimum is obtained in the NN1 geometry, in which the Al adjacent Al atoms sit on the top of both Na and Cl surface atoms. Clearly, this geometry in which both Al atoms of the supercell have the same distance from different types of surface atoms is physically not meaningful. Anyway, the fact that the minimum total energy of the NN1 wire is lower than that of the NN4 wire further supports the unstable character of the NN4 wire. In conclusion, the comparison of the results of these model wires with those for single isolated Al atom on the NaCl(100) surface reflect the balance between the Al-surface interactions and the intrawire interactions.

The upper panel of Fig. 6 shows that all the model wires are magnetic over a wide region around the equilibrium distance. It is interesting to note that the magnetism of the NN1 wire vanishes when the distance to the surface decreases. In order to understand the magnetic properties of these wires as well as those of the fully-relaxed wire on NaCl(100) we need to analyze their band structures. The equilibrium band structure of the NN2 wire along the wire direction is shown in Fig. 7. Comparing with Fig. 2 it can be seen that the 3$`p`$ bands around $`E_f`$ causing the magnetization have lost their degeneracy. However, in contrast to the unsupported zig-zag wires in Fig. 4 both of the spin-down 3$`p`$ bands are above the Fermi level. Therefore the magnetic moments of the wires attain values similar to those for unsupported linear wires. The energy bands for NN1 and NN4 geometries show similar features. In the NN2 geometry the occupied spin-up 3$`p`$ bands reach the Fermi level further away from $`k_z=0`$ than in the NN1 and NN4 geometries. This explains why the magnetic moment is largest in NN2 geometry. When the surface-wire distance is reduced, the unoccupied spin-down bands touch in the NN1 geometry the Fermi level and eventually their low-$`k_z`$ part sink below the Fermi level. This is the reason for the gradual decrease of the magnetic moment of the NN1 wire towards to shorter distances from the surface.

The band structure of the NN2 wire in Fig. 7 also shows that the Na and Cl bands do not disturb the band structure of the Al wire near the Fermi level. The Na 3$`s`$ bands are well above the Fermi level. The Cl 3$`p`$ bands form, at about 6 eV below the Fermi level, a 1.5 eV wide region which repells the Al 3s-bands. This repulsion does not affect the Al 3$`p`$ bands causing the magnetization. Addition of more NaCl layers into the supercell causes larger distortions in the Al 3$`s`$ band, but the above discussion about the origin of the magnetization will not be affected.

In case of the NN2 wire we have tested if dimerization could take place spontaneously. In the test the distance of the wire from the surface layer is kept constant corresponding to the previously determined equilibrium distance. According to Fig. 6 this is a fair approximation because the energy minima even for the very different geometries occur approximately at the same distance. Then the two Al atoms in a supercell are simultaneously displaced symmetrically in opposite directions to form a dimerized wire. Similarly as in the case of the unsupported linear wires, dimerization increases the energy of the wire. In contrast, a zig-zag deformation lowers the energy of the NN2 wire as in the unsupported situation. These notions conform the results obtained for the fully-relaxed Al wire on NaCl(100).

## IV Conclusions

We have performed density-functional computations to study electronic and magnetic properties of Al nanowires. Unsupported atomic chains of Al become magnetic along elongation in agreement with jellium results. When relaxing, straight wires deform to zig-zag wires spontaneously. The zig-zag wires exhibit also spontaneous magnetization, although the magnetic moments are lower than those of straight wires.

The results for the unsupported wires form the basis for understanding the properties of nanowires on solid surfaces. Many important features of the electronic structure are conserved when the wire is placed on the solid surface. Metallic nanowires are soft and the role of the surface is then to provide a mean to deform in a controlled way the geometry (e.g. the bond length) of the nanowire. As a concrete example we have considered commensurate Al wires on the NaCl(100) surface. A relaxation of the Al ion positions results in a zig-zag geometry and the wire conserves its magnetic moment. The resulting properties can be understood by studying the band structures of model wires on a NaCl layer.

Our calculations for Al wires on NaCl(100) are among the first studies for understanding properties of nanowires on solid surfaces. They show the importance of relaxing the ionic geometry of the wire. In order to gain deeper knowledge more geometries should be tried. For example, an Al wire at the step edge on the NaCl(100) surface would be interesting from the experimental point of view as a way to construct a straight wire. Interesting future topics could also be Al wires on surfaces of different materials, e.g. on silicon.

Studies of nanowires on crystal surfaces are of current scientific interest because such structures can be implemented in future nanotechnological devices. On the other hand, one-dimensional atomic chains on surfaces have been seen during the early stages of the growth of adsorbate layers and they have been shown to exhibit many interesting physical phenomena. We hope that our calculations have enlightened the richness of physical phenomena on these fields and will eventually encourage experiments searching for magnetism of simple metal nanowires on solid surfaces.

ACKNOWLEDGMENTS

A. Ayuela is a Marie Curie Fellow supported by the EU TMR program (Contract No. ERB4001GT954586). This research has been supported by the Academy of Finland through its Centre of Excellence Programme (2000 - 2005). We acknowledge the generous computing resources of the Center for the Scientific Computing (CSC), Espoo, Finland. The authors whish to thank G. Kresse for providing the VASP computer programs which made this work possible.


text:
# On the two types of steady hard X-ray states of GRS 1915+105.

## 1 INTRODUCTION

The X-ray source GRS 1915+105, the first Galactic object found to generate superluminal jets (Mirabel $`\&`$ Rodriguez 1994), was discovered by the GRANAT observatory as a transient in 1992 (Castro-Tirado, Brandt $`\&`$ Lund 1992). Long-term monitoring of GRS 1915+105 with the RXTE satellite (Bradt et al. 1993) has revealed a rich character of the source transient activity associated with different levels of its luminosity, which can be reduced to the sequence of varying soft steady, flaring and hard steady states qualitatively distinguished by their spectral and temporal properties (Greiner et al. 1996; Trudolyubov et al. 1999a; Muno et al. 1999). It was demonstrated that the major spectral and temporal properties of the source during these periods are similar to that of the Galactic black hole candidates in the ’canonical’ spectral states (Chen et al. 1997; Taam et al. 1997; Trudolyubov et al. 1999a).

For almost a half of the observational time the source demonstrates steady periods characterized by a relatively hard broad–band X-ray spectrum and X-ray luminosity $`(16)\times 10^{38}ergs/s`$ (assuming the source distance of $`12.5`$ kpc) – hard steady states (Muno et al. (1999); Belloni et al. (2000)). In general, the energy spectrum of GRS 1915+105 during these periods can be approximated by a sum of the weak soft thermal component with characteristic color temperature $`1`$ keV and dominant extended hard component with turnover at high energies (Fig. 1) (Trudolyubov et al. 1999a ; Muno et al. (1999)). The variation of the X-ray flux on the time scale of an individual observation does not exceed $`20\%`$. The power density spectrum (PDS) of the source is dominated by strong band-limited white noise (BLN) component accompanied by system of QPO peaks with centroid frequencies of several Hz (Fig. 2) (Morgan et al. (1997); Chen et al. (1997); Trudolyubov et al. 1999a ). The pattern of correlation between X-ray and radio emission of the source is complicated: alternating periods of high (so-called ’plateau states’) and low level radio emission associated with hard steady states have been reported (Foster et al. (1996); Fender et al. (1999); Muno et al. (1999)).

In this work we demonstrate that the source shows two distinct branches of the hard steady state (hereafter type I and type II), distinguished by the properties of its energy spectrum, high frequency X-ray variability and level of radio flux. We discuss the difference between these two types of the hard steady states in the framework of the two-phase model of the inner part of accretion flow, involving the outflow from the immediate vicinity of the compact object.

## 2 OBSERVATIONS AND RESULTS

The RXTE satellite performed regular pointed observations of GRS 1915+105 in 1996 – 2000 providing a good coverage for the detailed study of the spectral and timing properties of the source in the hard steady states. In the Table 1 we listed the groups of observations covering hard steady states and individual representative observations used for the detailed analysis.

### 2.1 Spectral and timing analysis

For the processing of the PCA and HEXTE data we used standard RXTE FTOOLS (version 5.01) tasks and methods recommended by RXTE Guest observer Facility.

For spectral analysis, we used PCA Standard 2 mode data, collected in the $`320`$ keV energy range. The PCA response matrices for individual observations were constructed using PCARMF v7.01, and background estimation was performed applying a Very-Large Events (VLE)–based model. Standard dead time correction procedure has been applied to the PCA data. Owing to the uncertainties of the response matrix, a 1$`\%`$ systematic error was added to the statistical error for each PCA energy channel. We used the HEXTE response matrices, released on October 20, 1998 and subtracted background collected in off-source observations for each cluster of detectors. In order to account for the uncertainties in the response and background determination, only data in the $`20150`$ keV energy range were used for the spectral analysis. Typical examples of the broad-band energy spectra of GRS 1915+105 during the hard steady states in 1996 – 2000 are shown in Fig. 1.

As we are interested in the general properties of the source energy spectrum during the hard steady state, we used only simple models to approximate its spectra. For the first group of observations showing a hard spectra with clear cut-off at the energies $`6090`$ keV (Fig. 1, panels a, b, c), we used only HEXTE data in the $`20150`$ keV energy range and approximated them by power law with exponential cut-off (Table 2, first part). For the remaining observations (Fig. 1, panels d, e, f) data of PCA and HEXTE instruments covering the $`3150`$ keV energy range were approximated by absorbed broken power law (Table 2, second part). For all cases there was a distinct residuals corresponding to the iron emission /absorption, thus Gaussian line profile and an absorption edge were included to refine the overall fit. An equivalent hydrogen absorbing column density was fixed at the level of $`N_HL=5\times 10^{22}`$ cm-2.

For the timing analysis in the $`230`$ keV energy range the RXTE/PCA data in the binned and event modes containing X-ray events below and above 13 keV respectively were used. We generated power density spectra (PDS) in the $`0.001512`$ Hz frequency range combining the results of the summed Fourier transforms of a short stretches of data ($`8`$ s) with $`0.001`$ s time bins for the $`0.3512`$ Hz frequency range, and a single Fourier transform on the data in $`0.125`$ s time bins for lower frequencies. The resulting spectra were logarithmically rebinned, when necessary, to reduce scatter at high frequencies, and normalized to square root of fractional variability rms. The white noise due to the Poissonian statistics corrected for the dead time effects, was subtracted (Zhang et al. (1995); Revnivtsev et al. (2000)). Typical examples of the broad-band power density spectra of GRS 1915+105 during hard steady state (in units of $`f\times (rms/mean)^2`$ (Belloni et al. (1997))) are shown in Fig. 2.

Broad-band power density spectrum (PDS) of the source is dominated by a strong band-limited white noise (BLN) component with at least two characteristic breaks (at $`0.1`$ and $``$few Hz) and a complex of relatively narrow peaks of quasi-periodic oscillations (QPO) lying near the second break in the slope of BLN continuum (Fig. 2). For some observations (Fig. 2, panels a, b, c), the additional broad high frequency component with characteristic cut-off at $`6080`$ Hz was also notable.

To quantify the properties of the source rapid aperiodic variability, PDS in the $`0.05256`$ Hz frequency range were fitted to analytic models using $`\chi ^2`$ minimization technique. We used a combination of two BLN components (approximated by zero-centered Lorentzian functions) and several QPO features (presented by the Lorentzian functions). For some observations (Fig. 2, panels a, b, c), an additional high frequency noise component was approximated by zero-centered Gaussian function with characteristic width of $`6080`$ Hz. Parameters of the PDS approximation are presented in Table 3.

### 2.2 Two types of X-ray hard steady states

As it is seen from Figs. 1 and 2, representative hard steady  state observations can be subdivided into two distinct classes with qualitatively different spectra and PDS. Further we will refer to these groups of observations as type I and type II states (see Table 2,3). Below we summarize characteristic spectral and timing properties of both groups.

Type I states. The energy spectrum is dominated by an extended power law hard component with photon index $`\alpha 1.82.3`$ and quasi- exponential cut-off at $`6090`$ keV (Fig. 1(panel a, b, c); Table 2). Broad-band PDS shows prominent noise component with average integrated rms amplitude of $`45\%`$ and characteristic cut-off frequency of $`6080`$ Hz, detectable up to $`150`$ Hz (Fig. 2, panel a, b, c; Table 3). An average total fractional rms amplitude of the rapid aperiodic variability in the $`230`$ keV energy range is about $`22\%`$ (Table 3).

Type II states. The energy spectrum exhibits clear break in power law slope near $`1316`$ keV. High energy part of the spectrum has roughly power law shape with almost constant photon index of $`\alpha 3.23.3`$, in spite of significant changes in the slope of the lower energy part (Fig. 1, panel d, e, f; Table 2). As in case of type I states, the PDS of the source is dominated by BLN and QPO components with characteristic peaks at a few Hz. Contrary to the type I states, high frequency part of the PDS in the type II state fades quickly lacking any significant variability above $`30`$ Hz. For some observations an additional noise component with characteristic cut-off at $`15`$ Hz was marginally detected (Table 3). The average total fractional rms amplitude of the rapid aperiodic variability in the $`230`$ keV energy range is $`5\%`$ lower than during the type I state observations (Table 3).

### 2.3 Correlated properties of X-ray and radio emission

To trace the simultaneous evolution of GRS 1915+105 in radio band we used publicly available data of the Green Bank Interferometer monitoring observations at 8.3 GHz (Foster et al. (1996)). Fig. 3 shows examples of the X-ray and radio lightcurves of the GRS 1915+105 during the hard steady state in $`19971998`$. We note that all periods of type I X-ray hard steady state correspond to low level of radio emission (the flux at $`8.3`$ GHz is typically lower than $`1020`$ mJy). 111Note, that the short ’spikes’ seen on the radio lightcurve correspond to the X-ray flaring states, occasionally interrupting hard steady states On the other hand, periods of type II state are characterized by much higher level of radio flux ($`50100`$ mJy). Type II state is also referred as a quasi- stable ’plateau’ radio state with flat radio spectrum (Foster et al. (1996); Fender et al. (1999)) and accompanying bright infrared emission (Bandyopadhyay et al. (1998)). To explain the properties of the source during ’plateau’ states the formation of the optically thick radio source powered by outflow of matter from the system was proposed (Fender et al. (1999)).

## 3 DISCUSSION

To explain the observational properties of black hole binaries in the hard and very high states (Tanaka & Lewin (1995)) a number of models involving the hot Comptonization region near the compact object surrounded by the optically thick accretion disk (Shakura & Sunyaev (1973)) were proposed (Chakrabarti & Titarchuk (1995)). According to these models the soft thermal component of the energy spectrum is emitted by the optically thick accretion disk, while the hot inner region is responsible for generation of the hard spectral component. It was proposed that low and high frequency noise components in the power density spectrum are associated with outer optically thick and inner optically thin parts of the accretion flow respectively (Miyamoto et al (1994)). It is often assumed that QPO phenomenon is caused by the interaction between these two distinct parts of the accretion flow occurring on the local dynamical time scale at the boundary between these regions related to the local Keplerian time (Molteni et al (1996); Titarchuk, Lapidus & Muslimov (1998)). Supporting this model, the evidence of the correlation between an observed QPO frequency and the position of the boundary between optically thick and optically thin parts of the accretion flow was reported for the flaring states of GRS 1915+105 (Trudolyubov et al. 1999b ; Chakrabarti $`\&`$ Manickam (2000)).

It looks natural to suggest that band limited white noise dominating the PDS of GRS 1915+105 may also be a product of the dynamic processes in the inner part of the accretion flow. Then the dynamic time scales on both inner and outer boundaries of this region will determine an extent of the resulting noise spectrum in the frequency domain. Given the position of high frequency cut-off in the PDS and radial dependence of the dynamic oscillation frequency, one can estimate an approximate position of the inner boundary of the optically thin region. Let us discuss individual properties of the two branches of the hard steady state on the basis of the above assumptions. We assume that characteristic frequencies of the QPO and BLN noise components of the power density spectra of the GRS 1915+105 are proportional to the Keplerian frequency at the outer and inner boundaries of the optically thin region of accretion flow, i.e. $`f=f_K/l`$, where $`f_K2200m^1r^{3/2}`$ Hz, $`m`$ – mass of the compact object in solar masses, $`r`$ – distance to the compact object in units of 3 gravitational radii.

Type I state: If the hard spectral component originates through the Comptonization of the low energy photons in the optically thin inner region, the observed cut-off energy of the spectrum indicates high electron temperature in this part of the accretion flow, $`kT_e2040`$ keV (Sunyaev & Titarchuk (1980)). Typical value of the QPO centroid frequency allows one to estimate an approximate position of the outer boundary of the optically thin part of the accretion flow: $`r_{out}\mathrm{\hspace{0.33em}170}(mf_{QPO}l)^{2/3}`$. Broad high-frequency component detected in all type I state observations is of particular interest. The position of the characteristic break of this component $`f_{cut}^{high}6080`$ Hz (Table 3) is remarkably close to the position of the reccurent stable QPO at 67 Hz observed mainly during the soft steady state (Morgan et al. (1997); Remillard & Morgan (1998)). This fact provides a possibility that the broad noise component can be related to this steady QPO and be interpreted as the effective broadening of the coherent oscillation feature. Most of the models explaining the stable 67 Hz QPO invoke some kind of oscillations near the inner edge of the accretion disk in the immediate vicinity of the compact object (Morgan et al. (1997) and references therein). These oscillations may explain the high-frequency noise component detected during the type I state observations.

Type II state: The properties of the X-ray and radio emission of GRS 1915+105 in type II state imply significantly different geometry and characteristics of the accretion flow compared to that of type I state. We suggest two possible interpretations of the observed shape of the hard spectral component in the type II state. Thermal Comptonization mechanism implies relatively low electron temperature in the inner region, $`kT_e5`$ keV. On the other hand, the observed break in the hard component can be explained as due to the transmission of the hard X-rays through dense surrounding matter with average Thomson optical depth of the order of several (assuming power law shape of the spectrum of illuminating radiation) (Sunyaev & Titarchuk (1980)). Typical value of the QPO centroid frequency during the type II states is close to that of the type I states, indicating possible similarity in the position of the outer boundary of the optically thin part of the accretion flow. As it is seen from Fig. 2, the PDS of the type II state demonstrates the lack of significant variability of the source at frequencies above $`30`$ Hz. This fact may hint at the larger inner radius of the hot comptonization region with respect to the type I state. Assuming noise generation mechanism described above, and given the characteristic cut-off frequencies of the PDS in the type I ($`f_I^{high}70`$ Hz) and type II ($`f_{II}^{high}15`$ Hz) state (Table 3), the ratio of the inner radii of the accretion flow, $`r_I^{in}`$ and $`r_{II}^{in}`$ can be estimated as $`(r_{II}^{in}/r_I^{in})(f_I^{high}/f_{II}^{high})^{2/3}3`$. To explain higher value of the inner radius in the type II state, the effective mechanism of the truncation of the accretion flow near the compact object in this state is needed. The observations at the radio and infrared wavelengths indicate the presence of a strong outflow of matter near the compact object coupled with a hot inner region (Fender et al. (1999); Eikenberry et al. (1998)). It is usually assumed that such outflows originate very close to the compact object, modifying the structure of the innermost part of the accretion flow. Finally, we may suppose that contrary to the type I state in the type II state the accretion flow extends down to some transition radius, followed by outflow region in the immediate vicinity of the compact object.

As a result, the difference between the properties of X-ray and radio emission of GRS 1915+105 during type I and type II hard steady states may be explained in terms of the different structure of the inner part of the accretion flow. Contrary to the type I X-ray state, in the type II state an inner region of the accretion flow coupled with a strong outflow of matter near the compact object, which causes the aforementioned differences in the observational properties.

Finally, our two types of the hard steady states may be compared with the previous systems of ’states’ used to characterize X-ray variability of GRS 1915+105. According to the definition, the ’C’ state of Belloni et al. (2000) (with four classes of variability $`\chi _1`$, $`\chi _2`$, $`\chi _3`$, $`\chi _4`$) corresponds to our hard steady state. It should be also noted that similar distinction for a limited sample of $`19961997`$ ’hard steady’ observations based on the level of radio flux was proposed by Muno et al. (1999). The ’radio-loud hard steady state’ and ’radio-quiet hard steady state’ of Muno et al. (1999) are directly associated with our type I and type II hard steady states.

This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center. The Green Bank Interferometer is a facility of the National Science Foundation operated by the NRAO in support of NASA High Energy Astrophysics programs.