text:
# Percolating through networks of random thresholds: Finite temperature electron tunneling in metal nanocrystal arrays

## Abstract

We investigate how temperature affects transport through large networks of nonlinear conductances with distributed thresholds. In monolayers of weakly-coupled gold nanocrystals, quenched charge disorder produces a range of local thresholds for the onset of electron tunneling. Our measurements delineate two regimes separated by a cross-over temperature $`T^{}`$. Up to $`T^{}`$ the nonlinear zero-temperature shape of the current-voltage curves survives, but with a threshold voltage for conduction that decreases linearly with temperature. Above $`T^{}`$ the threshold vanishes and the low-bias conductance increases rapidly with temperature. We develop a model that accounts for these findings and predicts $`T^{}`$.

Densely packed arrays of nanocrystals form a new class of “artificial solids” with tunable electronic, magnetic and optical properties . These properties stem from single-electron charging and quantum confinement energies on the individual particle level, mediated by the coupling to neighboring particles . However, even for the simplest such artificial solid, a close-packed monolayer of metallic nanoparticles, the effect of temperature on the overall transport properties is poorly understood.

At low temperatures, metal nanoparticle arrays show pronounced nonlinear current-voltage ($`IV`$) characteristics, with a well-defined threshold above which the applied bias has to be raised for conduction to occur. Similar threshold behavior occurs in a large family of disordered nonlinear systems, including the depinning of magnetic flux bundles in type-II superconductors  and the onset of motion in charge-density waves trapped by impurities. In metal nanoparticle arrays the disorder arises, even for perfectly ordered particle arrangements, from quenched background charges in the substrate or the matrix surrounding the particles. These quenched charges induce random shifts in the local nanoparticle potentials and lead to a distribution of electrostatic energy costs incurred by electrons as they tunnel from particle to particle. The system behaves as a network of ultrasmall, Coulomb-blockade-type tunnel junctions  with a distribution of threshold energies. For $`T=0`$, a key result was obtained a decade ago by Middleton and Wingreen (MW)  who showed through simulations and a scaling model that the effective voltage threshold $`V_t`$ is proportional to the array length and that the current-voltage ($`IV`$) characteristics exhibit power law scaling of the total current, $`I`$, with excess voltage, $`VV_t`$. These predictions are in line with experimental results on a wide range of nanocrystal arrays . More recently, the effect of structural disorder on the scaling behavior was investigated by experiments and reproduced by simulations .

For finite temperatures, on the other hand, the literature to date shows no consensus. In single-particle systems the shape of the $`IV`$ curves near threshold is exponentially sensitive to temperature . For arrays, Heath’s group reported that increasing temperature rapidly rounds out any nonlinearities . Black et al. also observed linear behavior above 40K . Other recent work found that finite temperature had a much weaker effect on the shape of the $`IV`$ characteristics and on the threshold $`V_t`$ . Theoretical approaches that could calculate temperature-dependent $`IV`$s for large arrays are currently not available.

To address this issue, we have performed systematic experiments on large, well-characterized two-dimensional arrays of gold nanocrystals. Our results show that the $`IV`$ curves retain their low-T, nonlinear scaling properties to remarkably high temperatures and that $`V_t`$ decreases only linearly with T. We can explain these findings by assuming that there is a subset of junctions, percolating across the array, for which thermal fluctuations effectively have removed the Coulomb blockade. The starting point is the picture developed by MW of charge fronts propagating through the array. However, our model goes beyond MW in three key aspects: it introduces a way to deal with the effect of temperature on a distribution of charging energies, it goes beyond square lattices and applies to arbitrary regular networks, and it can explicitly take into account capacitive coupling between neighboring islands, as is appropriate for close-packed arrays. We believe that this very general approach may also be applicable to other systems with distributions of local threshold values.

Our arrays consisted of monolayers of 1-dodecanethiol-ligated gold nanocrystals, synthesized and deposited from solution onto silicon nitride “window” substrates . Several arrays were deposited onto substrates with thicker $`\mathrm{Si}_3\mathrm{N}_4`$ layers and no windows, more suitable for measurements above 200K. Eleven arrays were studied ranging from monolayers with long-range hexagonal order, with domains of typical size equal to the electrode spacing, to monolayers with some disorder, including voids or localized regions of double layers not more than 20% in area fraction. In each sample, the nanocrystals were monodisperse to within 5% as checked by transmission electron microscopy (TEM). Different samples had mean core diameters $`4.5\mathrm{nm}<\mathrm{d}<6.5\mathrm{nm}`$, and mean spacings $`1.5\mathrm{nm}<\mathrm{s}<2.6\mathrm{nm}`$. The relevant Coulomb charging energy, $`\mathrm{\Delta }E`$, is the change in electrostatic energy of the system due to an electron tunneling from one particle to the next. $`\mathrm{\Delta }E`$ depends on the background charges on all particles involved. However, because of compensation by mobile electrons, it suffices to consider effective background charges within \[-e/2, +e/2\] and the distribution of charging energies, $`P(\mathrm{\Delta }E)`$, only extends to $`\pm \mathrm{\Delta }E_{max}`$. For close-packed monolayers, keeping only interactions between the two particles forming the junction and the surrounding 8 nearest neighbors, to good approximation $`\mathrm{\Delta }E_{max}=e^2(𝐂_{11}^1+2𝐂_{12}^1)=e^2𝐂_{11}^1(1+2\gamma )`$, with $`\gamma =𝐂_{12}^1/𝐂_{11}^1`$ where $`𝐂^1`$ is the inverse capacitance matrix. From measurements of d and s, and taking a dielectric constant $`ϵ=2`$ for the alkanethiol, the inverse capacitance matrix elements for such 10-particle system were calculated numerically . This gave $`𝐂_{11}^1`$ of order $`10^{18}\mathrm{F}^1`$ and typical values for $`\gamma 0.36`$. All arrays had a fixed width ($`2\mu `$m, or $`270`$ particles) and lengths, $`N`$, ranging from 27 - 170 particles, set by the electrode separation. For samples on window substrates, $`N`$ was measured by direct inspection with TEM after the transport measurements were completed. Otherwise, $`N`$ was estimated from the electrode separation. Samples were cooled in vacuum, and dc $`IV`$ curves were taken using Keithley 614 electrometers at voltage ramp rates $`525\mathrm{m}\mathrm{V}/\mathrm{s}`$.

At low temperatures, all arrays showed a clear voltage threshold for conduction, indicating strong Coulomb blockade behavior. Close inspection of the $`IV`$ curves (see Fig.1) reveals that the differential conductance, $`g\mathrm{d}I/\mathrm{d}V`$, for a given current $`I`$ is independent of temperature below a sample-dependent characteristic temperature in the range 100 K - 200 K. The $`IV`$ curves keep their shape and simply shift to lower voltages with increasing temperature. More quantitatively, we find that the $`IV`$s can be collapsed on top of each other by translating them by an amount $`V_{\mathrm{shift}}(T)<0`$. Defining an effective threshold $`V_\mathrm{t}(T)=V_\mathrm{t}(0)+V_{\mathrm{shift}}(T)`$, we can collapse all traces onto a master curve as shown in Fig.1b. $`|V_{\mathrm{shift}}(T^{})|=V_\mathrm{t}(0)`$ defines the temperature, $`T^{}`$, at which the array threshold first reaches zero. Once $`T>T^{}`$, the $`IV`$s exhibit non-zero differential conductance at zero bias (Fig.1a inset) prohibiting a data collapse at small bias. However, at higher bias the $`IV`$s still collapse when translated by an amount $`|V_{\mathrm{shift}}(T)|>V_\mathrm{t}(0)`$. The regime above $`T^{}`$ thus is characterized by a negative effective threshold and low-bias data falling below the master curve.

For highly ordered arrays (bottom curve in Fig.2a) MW’s $`T=0`$ scaling law, $`I(VV_t)^\zeta `$, remains valid up to $`T^{}`$ with temperature-independent exponent $`\zeta =2.25\pm 0.1`$. Details of the functional form of $`I(V)`$ in different arrays affect the value of $`T^{}`$ but are irrelevant as far as the collapse is concerned. Even for slightly disordered arrays, characterized by a somewhat s-shaped $`IV`$  in Fig.2, collapse can be achieved. For all arrays studied, $`V_\mathrm{t}(T)`$, the only adjustable parameter in the scaling procedure, drops essentially linearly with increasing temperature (Fig.2b). Both $`V_\mathrm{t}(0)`$ and the slope of this depression increase with array size $`N`$. Furthermore, for similar $`V_t(0)`$ the slope decreases with increasing charging energy.

This slow, linear depression of the voltage threshold with temperature is surprising and unexpected from the theory of few-junction networks  or from other nonlinear systems such as charge-density waves . A similar trend was first noted in 1D chains of carbon nanospheres by Bezryadin et al. , and by Ancona et al.  in small 2D arrays. Bezryadin et al. proposed that thermal fluctuations reduce the charging energy by $`k_BT`$ in all particles and that therefore $`V_t(T)`$ should decrease as $`Nk_BT/e`$. However, this simple picture cannot account for the slope values in Fig. 2b. We now outline a model that captures the linear dependence of $`V_t(T)`$, predicts the slope, and, under appropriate rescaling, collapses the data in Fig.2 onto a universal curve.

MW showed that the minimum bias voltage for charges to move all the way from one electrode to the other depends on the number of “up-steps” along any given path. Each local up-step in potential, no matter how small, requires an increment $`\mathrm{\Delta }V`$ in the externally applied bias for the charge front to advance, and thus contributes equally. ($`\mathrm{\Delta }V=\mathrm{\Delta }E_{max}/e`$ for uncoupled particles or junctions surrounded by charge neutral neighbors; for the more general case see below.) The array threshold is set by the optimal path with the smallest total number of up-steps and scales as $`V_t=\alpha N\mathrm{\Delta }V`$ in both 1D and 2D . In 2D, $`\alpha =0.338`$ for square lattices  and $`\alpha =0.227`$ for triangular lattices . Down-steps do not gain energy because all tunneling between particles is inelastic (no resonant tunneling between sharp levels). Hence, as long as the number of up-steps along a path remains fixed, a gradual reduction in the associated energy step size, induced by thermal fluctuations, will not reduce $`V_\mathrm{t}`$. We now argue that $`V_\mathrm{t}`$ can be reduced if temperature eliminates some of the up-steps completely. This occurs for the subset of junctions with step magnitudes $`|\mathrm{\Delta }E|`$ so small that they are washed out by thermal fluctuations.

We consider a scenario in which temperature has removed the nonlinear, threshold-type Coulomb blockade behavior of a fraction $`p(T)`$ of the junctions (irrespective of the sign of the step), effectively turning them into linear resistors. The effect of temperature, then, is a shift in threshold proportional to $`p(T)`$. At some temperature, $`T^{}`$, the fraction $`p(T^{})`$ will become large enough for a continuous path of junctions without discernable Coulomb blockade to form a percolating sub-network that bridges the array. At this point the overall array voltage threshold will vanish. We therefore identify $`p(T^{})`$ with the percolation threshold, $`p_c`$, for the underlying lattice. Hence we expect $`V_\mathrm{t}(T)=V_\mathrm{t}(0)(1p(T)/p_c)`$. For a triangular lattice $`p_c0.347`$. To make explicit the temperature dependence of the threshold voltage, we need to find $`p(T)`$. For a given distribution, $`P(\mathrm{\Delta }E)`$, junctions corresponding to $`|\mathrm{\Delta }E|<bk_\mathrm{B}T`$ are taken as effectively linearized. The factor $`b`$ accounts for the broadening of the Fermi-Dirac distributions in the nanocrystals. Electrons in the high-energy tail of one nanocrystal’s distribution can tunnel into the low-energy tail of accessible states of a neighboring particle, reducing the energy cost by roughly $`2k_\mathrm{B}T`$. More accurately, equating the reduction with the difference between the mean energies in the high- and low-energy tails, we calculate $`b=2.4`$ . At any given temperature, $`p(T)`$ is then given by the integral over $`P(\mathrm{\Delta }E)`$ between the limits $`\pm bk_BT`$.

Consider first the case of two particles with random charge offsets in the range $`[e/2,e/2]`$, forming a representative tunnel junction, surrounded by 8 neutral neighbors. In this case, $`\mathrm{\Delta }E_{max}=e^2𝐂_{11}^1(1\gamma )`$, and $`P(\mathrm{\Delta }E)`$ is a triangle with y-intercept $`P(0)=1/\mathrm{\Delta }E_{max}`$ (inset to Fig.3) . Integrating $`P(\mathrm{\Delta }E)`$ over the shaded area leads to $`p(T)=bx(2bx)`$ with $`x=k_BTP(0)`$. Using $`p(T^{})=p_c`$ and $`b=2.4`$ we find $`k_BT^{}0.08/P(0)`$. Note that $`p(T)`$ is proportional to $`T`$ and $`P(0)`$ as long as $`x<<1`$, which is satisfied up to temperatures far above $`T^{}`$. Allowing for random offset charges also on the 8 neighbors rounds the corners and steepens the slope of $`P(\mathrm{\Delta }E)`$. Still, the overall triangular shape remains, $`p(T)=2bx+𝒪`$$`(x^3)`$ and in the above considerations only the numerical value for $`P(0)`$ needs to be modified. Now $`\mathrm{\Delta }E_{max}=e^2𝐂_{11}^1(1+2\gamma )`$ and $`P(0)=(e^2𝐂_{11}^1)^1(11.57\gamma )/(1\gamma )^2`$ . For all arrays, independent of N, the model then predicts that the normalized threshold $`V_\mathrm{t}(T)/V_\mathrm{t}(0)`$ decreases linearly as a function of the dimensionless variable $`x=k_BTP(0)`$, with a slope $`2b/p_c`$, crossing zero near $`x=0.072`$. Within this picture, $`1/P(0)`$ can be thought of as an effective charging energy per up-step that captures all effects of the local packing geometry and the quenched charge distribution.

These predictions agree well with our data (Fig.3) . To obtain $`V_\mathrm{t}(0)`$ we extrapolated $`V_\mathrm{t}(T)`$ back to $`T=0`$ for each sample (Fig.2b). To determine $`P(0)`$ two methods were used, giving independent checks of the model. Equating $`e\mathrm{\Delta }V=1/P(0)`$ for the energy per up-step we obtained $`P(0)`$ from $`V_\mathrm{t}(0)=\alpha N\mathrm{\Delta }V`$ using $`N`$ from the electrode separations and $`\alpha =0.227`$. In addition, for the samples where TEM images of sufficient resolution were available to extract d and s, we were able to cross-check the values for $`P(0)`$ with those derived from direct numerical calculation of the $`𝐂_{ij}^1`$ .

The broad distribution of local energy differences $`\mathrm{\Delta }E`$ guarantees that, even if a fraction of steps is effectively removed, the underlying concept of current paths meandering through a random potential landscape - the premise of MW’s scaling picture - remains valid. The general shape of the $`IV`$s above $`V_\mathrm{t}(T)`$, which reflects how additional current paths open up with increasing bias beyond threshold, therefore should remain unchanged up to $`T^{}`$. This is borne out by our data. Rounding of the Coulomb blockade in the local junction $`IV`$ characteristics affects the quality of the data collapse only at the smallest $`VV_\mathrm{t}(T)`$ (Fig.2a).

The behavior for $`T>T^{}`$ is quite different. Here the low-bias conductance is dominated by a few junctions from the ($`p_c/2`$)-th percentile of $`P(\mathrm{\Delta }E)`$, since they are the key, ”bottleneck” bonds in the lattice during percolation at $`T^{}`$. Therefore, the low-bias $`IV`$s for the percolating path resemble those of a single island system with washed-out Coulomb blockade, for which the linear zero-bias conductance (Fig.1a inset) exhibits thermally activated behavior : $`g_0(T)\mathrm{exp}\left(U/k_BT\right)`$. Our model predicts $`U=bk_BT^{}=0.17\mathrm{\Delta }E_{max}`$, the largest $`\mathrm{\Delta }E`$ for the percolating path. This $`g_0`$ will dominate the low-bias array conductance as $`T`$ increases and additional, parallel paths open up that have larger $`U`$. Because of these additional paths, we expect the measured effective activation energy to be somewhat larger than $`bk_BT^{}`$. Our data indeed are compatible with simple activated behavior and a magnitude of $`U`$ 1 to 2 times $`bk_BT^{}`$. For example, for the array in Fig.1b inset $`\mathrm{\Delta }E_{max}=124`$ meV, giving $`bk_BT^{}=21`$ meV, to be compared with $`U38`$ meV from the slope of $`\mathrm{ln}(g_0)`$ vs $`T^1`$. However, because of the limited (high) temperature range available, our data cannot rule out other functional forms, such as variable range hopping conduction as seen in Ref..

Our results support a picture of thermally assisted quantum tunneling between neighboring gold nanocrystal cores through the dodecanethiol ligand. Within this picture, the ligands serve as a mechanical spacer between the particles and they set the tunneling barrier height. They do not introduce additional states into the barrier, nor is the barrier low enough to be thermally hopped over , both of which may be the case with more complex ligands . We find no evidence of a regime with positive temperature coefficient of resistance as seen in Ref. . Such behavior may require a stronger interparticle coupling and junction resistances much closer to the quantum resistance $`R_q=h/e^2`$. In the weak coupling limit, as in our arrays, the value of $`T^{}`$ as well as the collapse of the temperature-dependent threshold data onto a universal curve emerge naturally from considering the interplay of tunneling and quenched charge disorder.

We thank Qiti Guo and Benjamin Lauderdale for experimental assistance, and Eduard Antonyan, Terry Bigioni, Ilya Gruzberg, Alan Middleton, Toan Nguyen and Tom Witten for stimulating discussions. XML acknowledges support from DOE W-31-109-ENG-38. This work was supported by the UC-ANL Consortium for Nanoscience Research and by the NSF MRSEC program under DMR-0213745.


format: text

text:
# Simulations of a classical spin system with competing superexchange and double-exchange interactions

## I Introduction

Doped rare-earth manganites with general chemical formula Re1-xAxMnO3 (where Re denotes a rare-earth and A is typically Ca or Sr) have been shown to exhibit colossal magnetoresistance, and a wide range of other physical properties, which are strongly dependent upon doping exp1 . The physical mechanisms responsible for the unusual magnetoresistance properties of these materials have not been understood yet, and it has been suggested that double-exchange interactions zener  between Mn+3 and Mn+4 ions are important. However, there is much ongoing debate as to whether double-exchange alone can indeed give rise to colossal magnetoresistance or if further mechanisms have to be considered debate . It is thus important to further understand properties associated with the double-exchange theory.

Equal valence manganese ions in the parent compounds interact via a superexchange mechanism, whereas the different valence manganese ions, introduced through doping, are coupled via a double-exchange interaction exp1 ; zener ; wollan . We investigate a simple model of classical spins in which nearest-neighbor spins interact via a competing antiferromagnetic superexchange and a ferromagnetic double-exchange interaction, where the ratio of these two interactions mimic the effects of the doping $`x`$. Properties of the pure Heisenberg model (only superexchange interaction) are well-known heis  and recently a high-accuracy Monte-Carlo study obtained the paramagnet-ferromagnet transition temperature and associated static critical exponents of the pure double-exchange model alvaro . Recent studies of the phase diagram of double-exchange systems using a different approach from ours are given in Ref.arovas . We use a combination of Monte-Carlo simulations and ground-state calculations to determine the phase diagram of the model with varying ratios of the superexchange and double-exchange interaction strengths. We also use spin-dynamics techniques to study the dynamic structure factor and obtain the dispersion curve of the double-exchange model.

## II Model and Methods

The model considered here can be described by the Hamiltonian

$$=J_{AF}\underset{<𝐫,𝐫^{}>}{}𝐒_𝐫𝐒_𝐫^{}J_{DE}\underset{<𝐫,𝐫^{}>}{}\sqrt{1+𝐒_𝐫𝐒_𝐫^{}},$$
(1)
where $`𝐒_𝐫=(S_{𝐫}^{}{}_{}{}^{x},S_{𝐫}^{}{}_{}{}^{y},S_{𝐫}^{}{}_{}{}^{z})`$ is a three-dimensional classical spin of unit length at site $`𝐫`$. The first term, with $`J_{AF}<0`$, is the antiferromagnetic superexchange interaction and the second term describes a ferromagnetic double-exchange term with coupling constant $`J_{DE}>0`$ between nearest-neighbor pairs of spins. We consider $`L\times L\times L`$ simple cubic lattices with periodic boundary conditions, and divide the system into two sublattices, $`A`$ and $`B`$, with magnetizations $`𝐦_A`$ and $`𝐦_B`$, respectively. The scalar product is defined as $`𝐦_A𝐦_B=|𝐦_A||𝐦_B|\mathrm{cos}\theta y`$, where clearly $`1y1`$.

In the ground state, the condition of energy minimization leads to the expression

$$y=(J_{DE})^2/(4|J_{AF}|^2)1,$$
(2)
with $`|J_{AF}|/J_{DE}1/\sqrt{8}=0.35355\mathrm{}`$, and for $`|J_{AF}|/J_{DE}<0.35355\mathrm{}`$ the ground state is ferromagnetic. Eq.(2) admits the possibility of non-colinear phases and/or sublattice magnetizations with different magnitudes. Note that $`y=0`$ for $`|J_{AF}|/J_{DE}=1/2`$, which means that $`\theta =\pi /2`$ and/or that the magnetization in at least one sublattice is zero in this case.

For $`T>0`$ we use Monte-Carlo simulations to measure thermodynamic quantities of the model in order to determine the phase diagram. We used the Metropolis algorithm, lattice sizes $`L=12`$ and $`24`$, and typically we discarded 10000 Monte-Carlo steps for thermalization and used 10000-30000 steps in computing averages.

The dynamics of the spins are governed by the equations of motion

$$\frac{d}{dt}𝐒_𝐫=𝐒_𝐫\times (\stackrel{}{}_𝐫),$$
(3)
where the effective field

$$\stackrel{}{}_𝐫=J_{AF}\underset{𝐫^{}}{}𝐒_𝐫^{}+\frac{J_{DE}}{2}\underset{𝐫^{}}{}\frac{𝐒_𝐫^{}}{\sqrt{1+𝐒_𝐫𝐒_𝐫^{}}}$$
(4)
and the sums are over $`𝐫^{}`$ nearest-neighbor to $`𝐫`$. These equations were solved numerically using a method based on second-order Suzuki-Trotter decomposition of exponential operators krech , to a maximum time of $`t_{max}=440/J_{DE}`$.

The dynamic structure factor $`S(𝐪,\omega )`$ for momentum transfer $`𝐪`$ and frequency transfer $`\omega `$, observable in neutron scattering experiments, is given by

$$S^k(𝐪,\omega )=\underset{𝐫,𝐫^{}}{}\mathrm{exp}[i𝐪(𝐫𝐫^{})]_{\mathrm{}}^+\mathrm{}\mathrm{exp}(i\omega t)C^k(𝐫𝐫^{},t)\frac{dt}{\sqrt{2\pi }},$$
(5)
where $`C^k(𝐫𝐫^{},t)`$ is the space-displaced, time-displaced spin-spin correlation function defined, with $`k=x,y,`$ or $`z`$, as $`C^k(𝐫𝐫^{},t)=S_{𝐫}^{}{}_{}{}^{k}(t)S_{𝐫^{}}^{}{}_{}{}^{k}(0)S_{𝐫}^{}{}_{}{}^{k}(t)S_{𝐫^{}}^{}{}_{}{}^{k}(0).`$ The displacement $`𝐫`$ is in units of the lattice unit cell length $`a`$. More details on the spin-dynamics methods are given in Refs. kun ; alex .

## III Results

Our preliminary results show that $`|𝐦_A|=|𝐦_B|`$ at all temperatures for $`|J_{AF}|/J_{DE}<1/\sqrt{8}`$, decreasing smoothly from $`|𝐦_A|=|𝐦_B|=1`$ at low $`T`$, towards zero as $`T`$ increases. For $`|J_{AF}|/J_{DE}>1/\sqrt{8}`$ we find that at low temperatures $`|𝐦_A||𝐦_B|`$, as illustrated in Fig. 1 for $`|J_{AF}|/J_{DE}=0.9`$. For this range of interaction strengths our results indicate that $`|𝐦_A|=|(J_{DE})^2/(4|J_{AF}|^2)1|`$ and $`|𝐦_B|=1`$ at $`T=0`$, and Eq.(2) then implies that $`\theta =0`$ for $`|J_{AF}|/J_{DE}<1/2`$, $`\theta =\pi `$ for $`|J_{AF}|/J_{DE}>1/2`$, and the sublattice magnetizations reduce to $`|𝐦_A|=0`$ and $`|𝐦_B|=1`$ for $`|J_{AF}|/J_{DE}=1/2`$. The phase-transition temperatures are obtained from analyses of the sublattice magnetizations for $`L=24`$, hence there are some uncertainties in these estimates due to finite-size effects.

The phase diagram that we obtained from the analysis of the sublattice magnetizations, shown in Fig. 2, has five different regions: (i) ferromagnetic (F) for low $`T`$ and small values of $`|J_{AF}|/J_{DE}`$; (ii) antiferromagnetic (AF) for moderate $`T`$ and high values of $`|J_{AF}|/J_{DE}`$; (iii) paramagnetic (PM) at high $`T`$; and two low temperature phases, which we label as (iv) region I for $`1/\sqrt{8}<|J_{AF}|/J_{DE}<1/2`$, and (v) region II for $`|J_{AF}|/J_{DE}>1/2`$. Regions (i)-(iii) are characterized by $`|𝐦_A|=|𝐦_B|`$, with $`𝐦_A=𝐦_B`$ in the ferromagnetic phase and $`𝐦_A=𝐦_B`$ in the antiferromagnetic phase. In contrast, in regions I and II we find that $`|𝐦_A||𝐦_B|`$ and although our simulations suggest that $`\theta =0`$ and $`\pi `$ in regions I and II, respectively, we cannot rule out possible small deviations from the aligned and anti-aligned sublattice spin configurations in the respective regions. The question marks (?) in the phase diagram (see Fig. 2) indicate that the points at which the F-PM and AF-PM transition lines join other phase boundaries have not been determined accurately yet. For the pure double-exchange ferromagnetic model ($`J_{AF}=0`$), the critical temperature (F-PM transition point) is $`T_{c}^{}{}_{}{}^{DE}=0.74515(7)J_{DE}`$ alvaro , whereas for the pure Heisenberg model $`T_c^{SE}=1.442929(77)|J_{AF}|`$ heis  for either ferromagnetic or antiferromagnetic interaction.

Fig. 3 shows a comparison of the dispersion relations for the pure Heisenberg ferromagnet and for the pure double-exchange ferromagnetic model. The results at $`T=0`$ correspond to the linear spin-wave theory, and data at the respective critical temperatures are from our simulations.

## IV Conclusions

We studied a system of classical spins coupled via competing antiferromagnetic superexchange and ferromagnetic double-exchange interactions. As expected, when the former is much weaker than the latter, the system is ferromagnetic and in the reverse case the system is antiferromagnetic. In addition, we found two other low temperature phases, in which the sublattice magnetizations differ in magnitude and are close to being aligned or anti-aligned, depending on the ratio of the two interaction strenghts. At high temperature the system is paramagnetic. Comparing the dispersion curves we see that spin-waves in the pure double-exchange model occur at lower frequency transfer than in the pure Heisenberg ferromagnet.

Acknowledgments

This research was supported in part by NSF grant No. DMR-9727714.


text:
# 1 The main result

## 1 The main result

We announce a new proof of the uniform estimate on the curvature of solutions to the Ricci flow on a compact Kähler manifold $`M^n`$ with positive bisectional curvature.

Given a compact Kähler manifold $`M^n`$ of complex dimension $`n`$ and a Kähler metric $`\stackrel{~}{g}=\stackrel{~}{g}_{i\overline{j}}(x)dz^idz^{\overline{j}}`$ on $`M^n`$ of positive holomorphic bisectional curvature, consider the normalized Kähler-Ricci flow

$`{\displaystyle \frac{}{t}}g_{i\overline{j}}(x,t)=R_{i\overline{j}}(x,t)+g_{i\overline{j}}(x,t)`$ (1)
on $`M^n`$ with initial condition $`g_{i\overline{j}}(x,0)=\stackrel{~}{g}_{i\overline{j}}`$. By a result of Bishop-Goldberg , we can assume the Kähler class of the metric $`\stackrel{~}{g}`$ satisfies the condition:

$`[\stackrel{~}{\omega }]=[\stackrel{~}{\mathrm{\Sigma }}]=\pi c_1(M^n)`$ (2)
where $`\stackrel{~}{\omega }=(\sqrt{1}/2)\stackrel{~}{g}_{i\overline{j}}dz^idz^{\overline{j}}`$ and $`\stackrel{~}{\mathrm{\Sigma }}=(\sqrt{1}/2)\stackrel{~}{R}_{i\overline{j}}dz^idz^{\overline{j}}`$ are the Kähler form, the Ricci form of the metric $`\stackrel{~}{g}`$ respectively, while $`c_1(M)`$ denotes the first Chern class. Under the normalized initial condition (2), the first author  (see also Proposition 1.1 in ) showed that the solution $`g(x,t)=g_{i\overline{j}}(x,t)dz^idz^{\overline{j}}`$ to the normalized flow (1) exists for all time. Furthermore by the work of Mok  (and Bando  for $`n=3`$), the solution metric $`g(x,t)`$ is known to have positive bisectional curvature for any $`t>0`$. Our main result is on the uniform estimate of the curvature independent of $`t`$:

Main Theorem Let $`M^n`$ be a compact Kähler manifold of complex dimension $`n`$ and $`\stackrel{~}{g}=\stackrel{~}{g}_{i\overline{j}}(x)dz^idz^{\overline{j}}`$ be a Kähler metric on $`M^n`$ of positive holomorphic bisectional curvature satisfying condition (2). Let $`g(x,t)`$ be the solution to the normalized Kähler-Ricci flow (1) on $`M^n`$ with the initial condition $`g(x,0)=\stackrel{~}{g}`$. Then $`g(x,t)`$ is nonsingular. Namely, the curvature of $`g(x,t)`$ is uniformly bounded independent of $`t`$.

Whether the solution $`g`$ above is nonsingular or not has been one of the important open problems in the study of Hamilton’s Ricci flow on compact Kähler manifolds. For $`n=1`$, the main theorem is proved by Richard Hamilton  based on his Harnack inequality and entropy estimate for the Gaussian curvature. (see also the further improvements by B. Chow  and Hamilton ). In , Chen-Tian obtained the uniform estimate on the curvature of $`g(x,t)`$ based on a Moser-Trudinger type inequality in  for Kähler-Einstein manifolds and some new functionals. Note the existence of Kähler-Einstein metrics on $`M^n`$ follows from the solution to the Frankel conjecture by Mori  and Siu-Yau . In contrast, our proof of the main theorem does not rely on the exsitence of Kähler-Einstein metrics, but instead on first author’s Harnack inequality for the scalar curvature $`R(x,t)`$ which is a consequence of the Li-Yau-Hamilton estimate for the Kähler-Ricci flow , and a very recent local injectivity radius estimate of Perelman for the Ricci flow . Thus, this gives a crucial estimate needed in a new proof of the Frankel conjecture via Hamilton’s Ricci flow method. In next section, we’ll outline the proof of the Main Theorem. Complete details of the proof and discussions about convergence of $`g(x,t)`$ as $`t\mathrm{}`$ will appear elsewhere.

Remark Richard Hamilton and the first author  are able to prove the main theorem for $`n=2`$ by studying the geometry of singularity model from the blow up limits and certain dimension reduction argument.

## 2 Outline of the proof

Let $`g(x,t)`$ be the solution to the normalized Kähler-Ricci flow (1) on compact Kähler manifold $`M^n`$ with positive bisectional curvature whose initial metric satisfies condition (2). It is easy to see that under flow (1), the volume $`V=V(t)`$, the total scalar curvature $`_{M^n}R(x,t)𝑑\nu _t`$, and the average scalar curvature

$$r=r(t)=\frac{1}{V}_{M^n}R(x,t)𝑑\nu _t$$

of $`(M^n,g(x,t))`$ are all constant in $`t`$. In fact, $`r(t)=n`$ for all $`t`$. Clearly, it suffices to show that the scalar curvature of $`g`$ is uniformly bounded from above in $`t`$:

$$R(x,t)<C,xM^n,t[0,\mathrm{})$$

for some constant $`C>0`$ indepedent of $`t`$. To do this, let us first recall

The Harnack inequlaity for the scalar curvature (Cao ) For any $`x,yM^n`$ and $`0<t_1<t_2<\mathrm{}`$, the scalar curvature $`R(x,t)`$ of solution metric $`g(x,t)`$ satisfies the inequlaity

$`R(x,t_1)R(y,t_2){\displaystyle \frac{e^{t_2}1}{e^{t_1}1}}e^{\mathrm{\Delta }/4}.`$ (3)
Here $`\mathrm{\Delta }=\mathrm{\Delta }(x,y,t_1,t_2)=inf_\gamma _{t_1}^{t_2}|\gamma ^{}(\tau )|^2𝑑\tau `$, where the infimum is taken over all space-time curves from $`(x,t_1)`$ to $`(y,t_2)`$, and $`|\gamma ^{}(\tau )|`$ is the velocity of $`\gamma `$ at time $`\tau [t_1,t_2]`$.

Now for any $`xM^n`$ and $`t>1`$, set $`t_1=t`$ and $`t_2=t+1`$. We can find $`yM^n`$ such that $`R(y,t+1)=r(t+1)=n`$. It then follows from the Harnack inequality (3)

$`R(x,t)`$ $``$ $`R(y,t+1){\displaystyle \frac{e^{t+1}1}{e^t1}}e^{\mathrm{\Delta }(x,y,t,t+1)/4}`$ (4)
$``$ $`n(e+1)e^{\mathrm{\Delta }(x,y,t,t+1)/4}.`$
On the other hand, from (1) one can show that

$`\mathrm{\Delta }(x,y,t,t+1)e^2d^2(x,y;t)`$ (5)
where $`d(x,y;t)`$ denotes the geodesic distance between $`x`$ and $`y`$ with respect to the metric $`g_{i\overline{j}}(t)`$. In particular, by (4) and (5),

$`R(z,t)n(e+1)\mathrm{exp}({\displaystyle \frac{1}{4}}e^2)`$
for all $`zB_t(y,1)`$, the geodesic ball centered at $`y`$ of radius $`1`$ with respect to the metric $`g_{i\overline{j}}(t)`$. Then by Perelman’s no local collapsing result (Theorem 4.1 or its corollary in ), there exists a constant $`\beta >0`$ independent of $`t`$ such that the volume of the geodesic ball $`B_t(y,1)`$ has a uniform lower bound:

$`\text{V}ol(B_t(y,1))\beta .`$ (6)
From (6) one can, by a volume comparison argument of Yau (cf. ), deduce a uniform upper bound on the diameter $`d_t`$ of $`(M^n,g(x,t))`$: there is a constant $`D=D(\beta )>0`$ independent of $`t`$, such that for all $`t>1`$

$$d_tD.$$
(7)
But then, (4), (5), and (7) imply that, for any $`t>1`$ and any $`xM^n`$,

$`R(x,t)n(e+1)\mathrm{exp}(e^2D^2/4),`$
which is our desired uniform estimate.

Finally we should remark that the above argument also works for compact Kähler manifolds with nonnegative bisectional curvature.

Institute for Pure and Applied mathematics at UCLA, and Texas A&M University
E-mail: hcao@ipam.ucla.edu or cao@math.tamu.edu

Department of Mathematics, Zhongshang University, & The Institute of Mathematical Sciences, The Chinese University of Hong Kong

Email: blchen@math.cuhk.edu.hk & stszxp@zsu.edu.cn


text:
# Brane-world dark matter

## Abstract

We show that, in the context of brane-world scenarios with low tension $`\tau =f^4`$, massive brane fluctuations are natural dark matter candidates. We calculate the present abundances for both hot(warm) and cold branons in terms of the branon mass $`M`$ and the tension scale $`f`$. The results are compared with the current experimental bounds on these parameters. We also study the prospects for their detection in direct search experiments and comment on their characteristic signals in the indirect ones.

One of the most important open problems in astrophysics and cosmology is to identify the nature of dark matter. It has been known for a long time that the luminous matter observed in spiral galaxies is insufficient to explain their rotation curves. The existence of dark halos was proposed as a possible solution for the discrepancy (see Binney  and references therein), although at present, numerical simulations of the formation of such halos appear inconsistent with observations. On the other hand, different estimations of the total matter density in the universe from large scale motions, virial masses or cluster abundances, and the more recent Type Ia supernovae and CMB anisotropies observations agree in a value $`\mathrm{\Omega }_M=0.27\pm 0.04`$ PDG ; WMAP , which is much larger than the value of the total luminous mass density in the universe $`\mathrm{\Omega }_{lum}h=0.0060.002`$. In addition, the big bang nucleosynthesis (BBN) and WMAP results for the total baryonic content $`\mathrm{\Omega }_Bh^2=0.0224\pm 0.0009`$ imply that most of the matter in the universe is dark and non-baryonic (see PDG ; WMAP ; Kamion  and references therein).

A possible explanation of this puzzle is that the dominant component of dark matter consists of some non-relativistic (cold) stable and weakly interacting massive particles (WIMP) which decoupled from radiation early enough so that their relic abundances are important today. The possibility that the universe is dominated by hot dark matter seems to conflict with numerical simulations of structure formation. Thus, the only potential candidates within the known particles would be massive neutrinos. However a detailed analysis has excluded the three light and even one additional heavy fourth generation of Majorana or Dirac neutrinos. This fact has led to the search of cold dark matter candidates beyond the Standard Model (SM) Kamion .

There are two main such particles studied in the literature. On one hand axions which strictly speaking cannot be considered as WIMP’s since they are very light and produced non-thermally. On the other hand we have the lightest supersymmetric particle, which can be identified with a neutralino in most of the supersymmetric models. The latter is probably the most studied and best theoretically motivated dark matter candidate Kamion . However the large number of free parameters in supersymmetric theories make their predictions extremely model dependent. More recently, the existence of large extra dimensions has been proposed as a new setting for a possible solution to the hierarchy problem ADD . In this scenario, the SM fields are forced to live on a three-dimensional hypersurface (brane) whereas gravity is able to propagate on the higher $`D=4+N`$ dimensional bulk space. In this Brane World scenario (BWS) the fundamental scale of gravity is not the Planck scale $`M_P`$ but another scale $`M_D`$. Although in the original ADD model ADD , $`M_D`$ is supposed to be not too much larger than the electroweak scale, recently brane cosmology models have arised in which $`M_D`$ is much larger than the TeV (Langlois  and references therein). In this work we will consider the case of a general BWS without assuming any particular value for $`M_D`$.

In these models gravitons propagating through the bulk space give rise to a Kaluza-Klein (KK) tower of massive gravitons on the brane. These KK gravitons couple to the energy-momentum tensor of the SM fields $`T_{SM}^{\mu \nu }`$ and could be produced under the appropriate circumstances as real or virtual particles. Another important effect that is expected in the BWS is the presence of brane fluctuations since rigid objects do not exist in relativistic theories. In other words the brane should have some finite tension $`\tau =f^4`$. When these oscillations are taken into account two new effects appear GB . First of all, we have to introduce new fields, which for a homogeneous extra space, essentially represent the position of the brane in the bulk space $`(x^\mu ,y^\alpha \pi ^\alpha (x)/f^2)`$. The $`\pi ^\alpha (x)`$ fields are the Goldstone bosons (GB) corresponding to the spontaneous symmetry breaking (SSB) of the translational invariance produced by the presence of the brane (branons). It has been shown GB  that when these branons are properly taken into account, the coupling of the SM particles to any bulk field is exponentially suppressed by a factor $`\mathrm{exp}(M_{KK}^2M_D^2/(8\pi ^2f^4))`$, where $`M_{KK}`$ is the mass of the corresponding KK mode. As a consequence, if the tension scale $`f`$ is much smaller than the fundamental scale $`M_D`$, i.e. $`fM_D`$, the KK modes decouple from the SM particles. Therefore, for flexible enough branes, the only relevant degrees of freedom at low energies in the BWS are the SM particles and the branons. Similarly to other GB’s, branons are expected to be nearly massless and weakly interacting at low energies. Nevertheless, in general, translational invariance in the extra dimensions is not necessarily an exact symmetry and some branon mass $`M`$ is expected from such explicit symmetry breaking as shown in DoMa ; ACDM . This is similar to what happens to pions which are the GB corresponding to the SSB of the chiral symmetry of low-energy strong interactions. As gravitons do, branons couple to $`T_{SM}^{\mu \nu }`$, however in this case, the lowest order effective Lagrangian is ACDM :

$`_{Br}`$ $`=`$ $`{\displaystyle \frac{1}{2}}g^{\mu \nu }_\mu \pi ^\alpha _\nu \pi ^\alpha {\displaystyle \frac{1}{2}}M^2\pi ^\alpha \pi ^\alpha `$ (1)
$`+`$ $`{\displaystyle \frac{1}{8f^4}}(4_\mu \pi ^\alpha _\nu \pi ^\alpha M^2\pi ^\alpha \pi ^\alpha g_{\mu \nu })T_{SM}^{\mu \nu }`$
We see that branons always interact by pairs, they are stable and difficult to detect, since their interactions are suppressed by the tension scale $`f`$, and they are expected to be massive. Thus we arrive to the conclusion that the massive oscillations of the brane are natural candidates to dark matter in the BWS where $`fM_D`$. The dark matter problem has been considered in this scenario from a different point of view in Dvali  and also in models with universal extra dimensions CFM .

In order to calculate the thermal relic branon abundance, we will use the standard techniques given in Kolb  in two limiting cases, either relativistic (hot) or non-relativistic (cold) branons at decoupling. The evolution of the number density $`n_\alpha `$ of branons $`\pi ^\alpha `$, $`\alpha =1,\mathrm{},N`$ interacting with SM particles in an expanding universe is given by the Boltzmann equation:

$`{\displaystyle \frac{dn_\alpha }{dt}}=3Hn_\alpha \sigma _Av(n_\alpha ^2(n_\alpha ^{eq})^2)`$ (2)
where $`\sigma _A=_X\sigma (\pi ^\alpha \pi ^\alpha X)`$ is the total annihilation cross section of branons into SM particles $`X`$ summed over all final states. The $`3Hn_\alpha `$ term, with $`H`$ the Hubble parameter, takes into account the dilution of the number density due to the universe expansion. These are the only terms which could change their number density to the leading order. In fact, branons do not decay into other particles and since they interact always by pairs the conversions like $`\pi ^\alpha X\pi ^\alpha Y`$ do not change their number. Notice that we are considering the low-energy effective Lagrangian in (1) and assuming for simplicity that all the branons are degenerate. Accordingly, each branon species evolves independently, and in the following we will drop the $`\alpha `$ index. The total branon density will be just $`N`$ times that of a single branon. The $`\sigma _Av`$ term denotes the thermal average of the total annihilation cross section times the relative velocity. From (1), it includes, to leading order, annihilations into all the SM particle-antiparticle pairs. If the universe temperature is above the QCD phase transition ($`T>T_c`$), we consider annihilations into quark-antiquark and gluons pairs. If $`T<T_c`$ we include annihilations into light hadrons. For the sake of definiteness we have taken a critical temperature $`T_c170`$ MeV and a Higgs mass $`m_H125`$ GeV, although the final results are not very sensitive to the concrete value of these parameters.

Defining the new variable $`x=M/T`$, one finds from (2) that below the freeze-out temperature $`x_f`$ for which the annihilation rate $`\mathrm{\Gamma }_A=n_{eq}\sigma _Av`$ equals the expansion rate $`H`$, branons are decoupled from the thermal bath and their abundance remains frozen relative to the entropy density. We will denote $`g_{eff}(T)`$ and $`h_{eff}(T)`$ the effective relativistic degrees of freedom for energy and entropy densities respectively at temperature $`T`$. In order to calculate the decoupling temperature for hot relics, it is a good approximation to use the condition $`\mathrm{\Gamma }_A(x_f)H(x_f)`$. From the Friedmann equation in a radiation dominated universe we have: $`H(x)1.67g_{eff}^{1/2}M^2/(x^2M_P)`$ On the other hand, expanding $`\mathrm{\Gamma }_A(x)`$ for $`x3`$ and neglecting $`M`$, we find: $`\mathrm{\Gamma }_A^\gamma (T)16\pi ^9T^9/(297675\zeta (3)f^8)`$ for photons and $`\mathrm{\Gamma }_A^\nu (T)=\mathrm{\Gamma }_A^\gamma (T)/4`$ for neutrinos. For massive particles we cannot give closed expressions. Once we know $`x_f`$, the corresponding fraction of energy density today in the form of relics is given by:

$`\mathrm{\Omega }_{Br}h^2{\displaystyle \frac{7.8310^2}{h_{eff}(x_f)}}{\displaystyle \frac{M}{\text{eV}}}`$ (3)
For cold relics the calculation of the decoupling temperature is more complicated. The well-known result is:

$`x_f=\mathrm{ln}\left({\displaystyle \frac{0.038c(c+2)M_PM\sigma _Av}{g_{eff}^{1/2}x_f^{1/2}}}\right)`$ (4)
where $`c0.5`$ is obtained from the numerical resolution of the Boltzmann equation. The above equation can be solved iteratively. The matter density can be written as:

$`\mathrm{\Omega }_{Br}h^28.76610^{11}\text{GeV}^2{\displaystyle \frac{x_f}{g_{eff}^{1/2}}}\left({\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{c_n}{n+1}}x_f^n\right)^1`$ (5)
where we have expanded $`\sigma _Av`$ in powers of $`x^1`$ as $`\sigma _Av=_{n=0}^{\mathrm{}}c_nx^n`$. In the case of photons, the first non-vanishing coefficient is $`c_2^\gamma =68M^6/(15f^8\pi ^2)`$ and for massless neutrinos $`c_2^\nu =c_2^\gamma /4`$ ($`d`$wave annihilation) whereas for non-conformal matter we also have $`s`$ and $`p`$wave annihilation, ($`c_0,c_10`$). The corresponding expressions are more complicated and will be given elsewhere. We have performed all the expansions up to $`𝒪(x^2)`$. Coannihilation effects are absent in this case since there are no slightly heavier particles which eventually could decay into the lightest branon. Also, in order to avoid the problems of the Taylor expansion near SM thresholds, we have taken branon masses sufficiently separated from SM particles masses where the usual treatment is adequate Kolb . Such treatment is known to introduce errors of the order of 10$`\%`$ in the relic abundances.

Branons could be responsible for the observed cosmological dark matter density provided $`\mathrm{\Omega }_{Br}h^2=0.1290.095`$ at the 95$`\%`$ C.L. which corresponds to $`\mathrm{\Omega }_{Br}=0.23\pm 0.08`$ and $`h=0.790.65`$ WMAP . In FIG. 1, we have plotted these curves for hot and cold branons in the $`fM`$ plane for one single branon. For $`N`$ types of branons the corresponding abundances are simply $`N`$ times larger. In fact, the contribution of branons to $`g_{eff}`$ in (5) is negligible in the case of cold relics. For hot relics such contribution in (3) has been taken into account, although it is very small in the interesting regions. Concerning the freeze-out temperature, the results on the cold dark matter curve range from $`x_f9`$ for $`M=10^5`$ GeV to $`x_f31`$ for $`M=10^6`$ GeV. In the hot relics case, a very good approximation is given by $`\mathrm{log}(T_f/\text{GeV})(8/7)\mathrm{log}(f/\text{GeV})3.2`$.

Pure hot dark matter models are disfavored at present because relativistic matter free streams from overdense into underdense regions preventing structures from growing below the so called free-streaming scale given by Colin : $`\lambda _{FS}0.2(\mathrm{\Omega }_{Br}h^2)^{1/3}(\text{keV}/M)^{4/3}\text{Mpc}`$. In the hot branons curve in FIG. 1, the masses in the allowed region are in the range $`M=85177`$ eV, which corresponds to $`\lambda _{FS}2.41.0`$ Mpc. Such scales are much smaller than those in neutrino dark matter models and, in addition, since these branons decouple much earlier than neutrinos do, their corresponding temperatures are also lower, i.e. they could be considered rather as warm dark matter (WDM) candidates Colin  from the point of view of structure formation. However, a WDM dominated universe seems to be also disfavored by the recent observations of WMAP WMAP . On the other hand, the presence of relativistic branons during BBN could change the expansion rate of the universe, spoiling the predictions of the light elements abundances Kolb . However, the BBN limits on the number of branons contributing to $`g_{eff}`$ are not very constraining. In fact, imposing the conservative bound $`g_{eff}(T_{BBN}1\text{MeV})<12.5`$ Kolb , we find $`N<18`$ for branons which decoupled before the QCD phase transition, corresponding to $`f>60`$ GeV.

Let us clarify the main assumptions used so far in the work. We are considering models with $`fM_D`$, and also assuming that the evolution of the universe is standard up to a temperature around $`f`$. Indeed, this is the case of realistic brane cosmology models Langlois . Therefore, we have taken the conservative bound $`f>1`$ MeV, so that BBN is not affected. Moreover the effective Lagrangian (1) is only valid at low energies relative to $`f`$. We have checked that our calculations are consistent with these assumptions since the decoupling temperatures are always smaller than $`f`$ in the allowed regions in FIG. 1.

Brane fluctuations could be, not only candidates for the cosmological dark matter, but also they could make up the galactic halo and explain the local dynamics. In such case, they could be detected in direct search experiments from the energy transfer in elastic collisions with nuclei of a suitable target. The appropriate quantity to be compared with the experimental results is not the elastic branon-nucleus cross section $`\sigma `$, but the differential cross section per nucleon at zero-momentum transfer $`\sigma _n`$, which is defined by Kamion :

$`{\displaystyle \frac{d\sigma }{d|q|^2}}={\displaystyle \frac{\sigma _nA^2F^2(|q|)}{4v^2\mu ^2}}`$ (6)
where $`\mu =Mm/(M+m)`$, $`F(|q|)`$ is a nuclear form factor with the normalization $`F(0)=1`$, $`m939`$ MeV is the nucleon mass, $`v`$ is the relative velocity and $`A`$ is the mass number of the nucleus. In the limit in which the momentum transfer goes to zero, we can consider the nucleons as pointlike particles. In this case, it is possible to calculate the branon-nucleon cross section $`\sigma _n`$ from (1) just considering the nucleon as a Dirac fermion of mass $`m`$:

$`\sigma _n={\displaystyle \frac{9M^2m^2\mu ^2}{64\pi f^8}}`$ (7)
In fact, this quantity does not depend on the type of particle which couples to the branon, but only on its mass. This can be seen from (1) since in this limit, branons only couple to the $`T^{00}`$ component.

The results of our analysis are shown in FIG. 2. Lines of constant $`f`$ with 50 GeV separation are shown for reference. The area on the left of the $`\mathrm{\Omega }_{Br}h^2=`$$`0.1290.095`$ curves is excluded by branon overproduction, but the right portion is compatible with observations and will be explored in future experiments. Such region corresponds to $`f>120`$ GeV and $`M>40`$ GeV.

Another interesting possibility is to detect branons indirectly. Their annihilations in the galactic halo can give rise to pairs of photons or $`e^+e^{}`$ which could be detected by $`\gamma `$rays telescopes such as MAGIC or GLAST or antimatter detectors (AMS). Annihilation of branons trapped in the center of the sun or the earth can give rise to high-energy neutrinos which could be detectable by high-energy neutrino telescopes such as AMANDA, IceCube or ANTARES (see for example CFM ). Because annihilations of non-relativistic branons into conformal matter are $`d`$wave suppressed, the most relevant contribution will come from the secondary leptonic decays of ultra-relativistic $`Z`$ or $`W`$ bosons (in the case $`MM_Z`$). These processes will be characterized by the presence of peaks around one half of the branon mass in the leptonic or neutrino spectra. In the case of photons, softer peaks will be present at lower energies and therefore their detection will be more difficult. The hadronic decays will give rise to relatively smeared spectra at lower energies. Detailed results will be presented elsewhere.

Throughout the paper we have assumed that branons were in thermal equilibrium with radiation at some point in the history of the universe. If this is not the case, branons could still be produced non-thermally, very much in the same way as axionic dark matter Wilczek . In fact for very light branons, the energy density produced by this mechanism could be cosmologically important.

In conclusion, we have proposed branons as natural dark matter candidates in the BWS with low tension. Our results show that in a certain range of the parameters $`f`$ and $`M`$, their relic abundances could explain the missing mass problem, and that such parameters region will be explored in future direct detection experiments.

Acknowledgements: This work has been partially supported by the DGICYT (Spain) under the project numbers FPA 2000-0956 and BFM2000-1326.


text:
# Richtmyer-Meshkov instability and solid 4He melting driven by acoustic pulse

## Abstract

Recent experiments have shown remarkable dynamics of solid 4He melting and growth, driven by the normal incidence of an acoustic pulse on the solid-liquid interface. The theory of solid growth/melting, driven by the radiation pressure of the acoustic pulse, accounts well for the temperature dependence of the measured data. There is however an observed source of extra, temperature-independent, melting. We here propose that this extra melting is due to solid-liquid mixing (and consequent melting) at the interface, in a process similar to the Richtmyer-Meshkov instability: Initial undulations of the rough interface, grow when accelerated by the acoustic pressure oscillations. This model predicts a temperature-independent extra melting and its dependence on the acoustic power, which is in agreement with the measured data.

Recent experiments nomura  have shown that the radiation pressure $`P_r`$ of an acoustic pulse, incident on the solid-superfluid interface of 4He, can induce the local melting or growth of the solid. The general temperature dependence of the pulse-induced melting/growth is accounted for using the concept of acoustic radiation pressure. Nevertheless, a large deviation was found, which seems to be temperature independent. It is this extra melting which we shall describe in this paper.

We begin with the analysis of the experimental data in terms of a linear growth/melting coefficient (Eq.1-6 of nomura ). The amplitude of melting/growth is linearly related to the radiation pressure $`P_r`$ through the (temperature dependent) growth coefficient $`K(T)`$ nomura

$$h=K(T)\frac{P_r}{\rho _c}\tau $$
(1)
where $`\tau =1`$msec is the overall duration of the acoustic pulse. The radiation pressure is itself temperature dependent: $`P_r=E\left[1\frac{c_1}{c_2}+R^2\left(1+\frac{c_1}{c_2}\right)\right]`$, where the energy density of the incident wave is $`E=I/c_1`$, the power density $`I=144`$W/m2, and the reflection coefficient is

$$R=\frac{z_2z_1z_1z_2\rho _1K(T)\left(\frac{\rho _1\rho _2}{\rho _1\rho _2}\right)^2}{z_1+z_2+z_1z_2\rho _1K(T)\left(\frac{\rho _1\rho _2}{\rho _1\rho _2}\right)^2}$$
(2)
where $`z_{1,2}=\rho _{1,2}c_{1,2}`$. The growth coefficient $`K(T)`$ is given by

$$K(T)^1=AT^4+Be^{(\frac{\mathrm{\Delta }}{T})}$$
(3)
(with $`A=3\times 10^2`$,$`B=2.4\times 10^3`$,$`\mathrm{\Delta }=7.2`$K bowley ). In these expressions $`\rho _{1,2},c_{1,2}`$ are the density and sound velocities of the two phases on either side of the interface, with the index 1 corresponding to the phase on the side of the incident wave. In the experiment nomura  we have $`\rho _l=0.17`$gr/cc, $`c_l=365`$m/sec and $`\rho _s=0.19`$gr/cc, $`c_s=475`$m/sec. In Fig.1 we compare the experimental data with the calculated growth/melting amplitude $`h`$, given by Eq.(1). The overall agreement is good (dotted lines in Fig.1), but there is clearly an additional source of melting which appears for pulses incident from both the solid and liquid sides. This extra melting is not caused by the sound heating of the solid nomura0 .

The contribution from the acoustic pressure leads to growth/melting that proceeds at a constant rate for the duration of the applied pulse (1). In addition there is the inertial acceleration of the solid-liquid interface by the incident pressure oscillations. This motion should average to zero if no instability or mixing occurs at the interface. It is only due to mixing at the interface that this mechanism can contribute on average to the position of the interface at the end of the pulse. We assume here that the shape changes of the interface due to Richtmyer-Meshkov (RM)-mixing rychtmyer  do not change the average acoustic pressure and therefore can be approximately decoupled from each other. We therefore add the contribution from the acoustic pressure to the mixing contribution due to the RM-mixing. A complete description has to take into account possible interference effects of one mechanism with the other.

Let us now summarize the main results of this paper. We propose that the measured extra melting is the result of dynamic solid-liquid mixing at the interface. This mixing is due to a RM-instability rychtmyer , that typically occurs when two liquids/gases of different densities are accelerated one against the other by an incident shock-wave, causing initial perturbations of the interface to grow. In our case the interface is between two phases of the same material, but of different density, which is repeatedly accelerated to and fro by a series of rapid pressure oscillations. Under these conditions, and a rough initial interface, solid-liquid mixing should occur due to a RM-type instability rychtmyer , since we note that in the limit of weak shock-waves these behave as acoustic pulses landau . Using the linear growth theory of the RM instability we predict a temperature-independent melting, with the correct dependence on the incident acoustic power. With a single free parameter we get good quantitative account of the measured extra melting.

A RM instability will induce a layer of mixing at the interface between two materials of different densities (here the solid-liquid interface). The mixing occurs when initial perturbations of the interface between two liquids (or gases) grow after the interface is accelerated by a shock-wave rychtmyerjpg  (Fig.2). The linear single-mode theory describes the inertial growth of an initial interface undulation, due to an impacting shock-wave, gives rychtmyer

$$\frac{da}{dt}=ka_0vAa_{RM}=a_0+ka_0vA\tau $$
(4)
where $`k`$ is the wavevector of the interface undulation of initial amplitude $`a_0`$, $`v`$ is the interface velocity due to the acceleration and $`A=(\rho _2\rho _1)/(\rho _2+\rho _1)`$ is the Atwood number. The sign of $`A`$ relates to the phase-reversal process (Fig.2b) and is irrelevant here. Note that Eq.(4) is temperature-independent. Since in our case the solid-liquid interface is continuously accelerated to and fro by the pressure oscillations noziers2 , we therefore interpret the equation for the rate of growth of surface instabilities (4), as giving the rate at which a continuous flux of solid is locally mixed with the liquid and melts.

Furthermore, unlike the usual RM scenario, we have a solid crystal on one side of the interface. Due to the rapid melting-solidification at the rough (i.e. unfaceted) solid-superfluid interface, this interface deforms similar to that between two liquids andreev . On the other hand, dissipation and surface tension tend to damp out any inertial growth of instabilities. The solid-liquid surface tension, is $`\sigma 0.2`$erg/cm2 for 4He surftenshe , and acts to stabilize and damp the inertial growth of initial perturbations of the RM-type surftens . The regime where the surface tension dominates occurs whenever the damping rate $`\omega =k^{3/2}\sqrt{\sigma /(\rho _1+\rho _2)}`$ is larger than the growth rate given in Eq.(4), i.e. for $`k>(a_0vA)^2(\rho _1+\rho _2)/\sigma `$. Using the relevant values (estimating $`a_010`$Å), we find that in the experiments the surface tension dominates. This means that following a single pressure oscillation, the amplitude growth described in (4) would have been quickly damped out, and result in negligible solid-liquid mixing.

Nevertheless, in the experiment nomura , there is not a single pressure oscillation followed by inertial growth, but rather a long pulse of continuous acoustic waves. Each pressure oscillation (of $`0.1\mu `$sec duration, frequency $`9.9`$MHz nomura ) will cause the solid-liquid interface to accelerate, and existing undulations to grow in amplitude. The solid-liquid interface is therefore continuously accelerated back and forth, throughout the pulse duration. The damping due to the surface tension is therefore not relevant, and we proceed to estimate the amount of solid-liquid mixing due to this continuously driven RM-instability, using the simplest linear growth theory rychtmyer  (Eq. 4), over the entire pulse duration $`\tau `$. In this case we treat the growth rate equation (4) as a measure of the rate at which solid material is continuously mixed into the liquid, in a manner similar to the (highly exaggerated) ”spike” in Fig.2. This material then melts as it finds itself in a region which is liquid in equilibrium. The accumulated amplitude $`a_{RM}`$ represents therefore the total thickness of solid that was mixed (and melted) by this process over the duration of the acoustic pulse. In particular, it does not correspond to the final amplitude of individual interface undulations.

The interface velocity $`v`$, due to the oscillating pressure of amplitude $`\delta P=\sqrt{E\rho _{1,2}c_{1,2}^2}`$ (depending on the direction of the incident wave), is given in terms of the local material velocities at the interface, after the passing of the pressure wave noziers2 : $`v=(\rho _1v_{m,1}\rho _2v_{m,2})/(\rho _1\rho _2)`$. The local material velocities following a pressure oscillation are given by noziers2  $`v_m=\delta P/z_{1,2}`$. For power density $`I=144`$W/m2 nomura  we therefore get: $`v0.023`$m/sec.

In Fig.1 we compare the calculated combined melting/growth $`ha_{RM}`$ (1,4) with the experimental data. We find that a good agreement note  is achieved when we use an interface roughness aspect-ratio of $`ka_020`$ in Eq.(4). For a rough solid-liquid interface, which is the case in the experiments done so far nomura ; nomura0 , we would naively expect $`ka_01`$. This discrepancy indicates that the linear growth theory is not entirely adequate, especially as it usually deals with a single pressure pulse, while in our case the interface is continuously driven. The use of a constant amplitude $`a_0`$ in the rate equation (4) is a rough approximation, as the amplitude of surface undulations grows between each pressure oscillation. We therefore note $`ka_0`$ as the free parameter in our model, and takes into account our ignorance of the details of the solid-superfluid interface dynamics. Note that the agreement is better for the case of melting with an acoustic pulse from the liquid side (Fig.1).

Our model further predicts that $`a_{RM}\delta v_m\sqrt{I}`$, while the radiation-pressure driven growth (1) predicts a linear dependence on the acoustic power. In the regime where the melting is dominated by the RM-instability, i.e. at relatively high temperatures, we therefore expect to find $`a_{RM}^2I`$ behavior. The experimental data at 1.2K nomura0  does indeed show such a power law (Fig.3).

To conclude, we have analyzed recent experiments where the solid-superfluid interface is accelerated by an acoustic-pulse, and a large source of unexplained melting was observed. It was suggested that some mechanism of energy dissipation at the interface is responsible for this unexplained melting nomura . We proposed here a possible mechanism for such a dissipation; local solid-liquid mixing (and melting) occurs due to a dynamic instability of the interface, of the RM-type. The proposed model is in good agreement with the experimental data, and describes a specific mechanism by which the accelerated interface melts the solid.

It will be interesting to study this instability using simulations rychtmyer , taking into account both the material flows and the solid-liquid melting/solidification dynamics. Experiments on faceted crystals, with an atomically smooth interface, should show less melting due to the RM-like process described here. It could also be useful to repeat these experiments with a controlled initial undulation on the solid-liquid interface, of macroscopic size, and visualize its evolution when accelerated by an acoustic pulse. This initial undulation can be formed by a standing capillary wave at the interface capilary .

###### Acknowledgements.

I wish to thank R. Nomura and Y. Okuda for access to their data, many helpful discussions and their hospitality at the ULT conference in Kanazawa, Japan (2002). This research is being supported by the Louis and Anita Perlman Family Foundation and the Koshland Postdoctoral Fellowship.