text:
# 1 Introduction

## 1 Introduction

Gamma ray bursts (GRBs) are amongst the most distant, energetic and enigmatic astrophysical phenomena known. Understanding GRBs is arguably the most outstanding question in astronomy, and one which may be answered by a plethora of gamma ray observatories such as INTEGRAL, SWIFT and BATSE, and corresponding infrared and optical telescopes such as REM and LT . It is well known that 99% of the energy of a supernova is emitted in the form of neutrinos, and therefore it is widely expected that GRBs are similarly a copious source of neutrinos which may be detected in future neutrino telescopes . Within this decade it is therefore likely that GRBs will become much better understood, and their exact nature and mechanisms which drive their internal engines will be revealed. For example it may turn out that a GRB results from the core collapse of a very massive supernova to a compact rotating black hole with the energy emitted in beamed relativistic fireball jets containing copious neutrino fluxes  Alternatively the GRB engine could result in the emission of beamed earth-sized cannonballs .

In this paper we are not concerned with detailed models of GRBs, but instead regard them as a high intensity, high energy neutrino beam with a cosmological baseline. We shall be interested in the time delay of the arrival of neutrinos. We show that the time delay may be used as a probe of three physical effects: (i) neutrino mass, (ii) quantum gravity, (iii) dark energy. The time delay due to neutrino mass has been noted earlier in  while effects of quantum gravity on the time of flight have been considered in  for high energy photons and in  for neutrinos. However none of the above papers considers the time delays due to massive neutrinos in the presence of quantum gravity, although  gives the dispersion relation for this case. Moreover no studies to date have calculated time delays using a formalism which correctly takes into account the time dependence of the Hubble constant due to matter and dark energy. In our study we consider the time dependence of the Hubble constant due to matter and dark energy and show that they can change the naive results by more than 100% for $`z>1`$. For large redshifts the results for time delays due to neutrino mass and quantum gravity are sensitive to the nature of dark energy. We show that the three effects may be disentangled by using low energy neutrinos to study neutrino mass, high energy neutrinos to study quantum gravity, and using large redshifts to study dark energy, leading to the results stated in the abstract.

For the determination of the neutrino mass, in principle one could compare the arrival time of the massive neutrinos with the arrival times of the photons emitted in the GRB, assuming them to be emitted at the same time. However this would be plausible only if the GRB were a point source in vacuo, which it is not. In any realistic GRB model the photons are trapped inside the fireball and are released much later – the exact amount of time delay being highly model dependent. Another strategy might be to compare the arrival times of the low energy neutrinos with that of the ultra-high energy ones, since (as we will show) the ultra-high energy neutrinos suffer negligible time delay due to mass. But again one would expect the low energy neutrinos to be produced thermally, as in supernovae, leading to a model dependent time delay.

However there are alternative strategies which could overcome these problems. To begin with, if the neutrinos are hierarchical in mass, $`m_1m_2m_3`$, then due to their mixing, neutrinos of the same energy will arrive at the detector in three bunches, corresponding to the three mass eigenstates, the arrival time only depending on their mass. We can then compare the arrival times of the different neutrino mass eigenstates and put limits on the neutrino mass. This represents a “clean strategy” and indeed similar arguments for constraining neutrino mass from time delays using supernova neutrinos have been used in the literature (see  and references therein). However we will show that with GRB neutrinos the sensitivities possible are better by orders of magnitude.

For constraining models of quantum gravity and dark energy we can use arrival times of the ultra-high energy neutrinos, which are produced over many decades of neutrino energy. We can compare the arrival times of the high energy neutrinos of different energies – which unlike the low energy thermal neutrinos are expected to be produced at almost the same time at the source – and put limits from such observations. Thus, the results in this paper do not rely on the comparison of the arrival times between photons and neutrinos, which would involve the uncertainties discussed above. However for simplicity, we shall calculate arrival time delays relative to a hypothetical low energy photon which is assumed to be emitted at the same time. This is for convenience only; in a realistic search strategy what will be important will be the comparison of time delays between neutrinos of either the same energy, or between high energy neutrinos of two different energies as discussed above.

## 2 Formalism

For a neutrino of mass $`m`$ and energy $`E(t)`$, the momentum $`p(t)`$ in the presence of the effects of quantum gravity (QG) with an effective energy scale $`E_{QG}`$ is given to lowest order as

$`p^2(t)c^2E^2(t)\left(1\xi {\displaystyle \frac{E(t)}{E_{QG}}}\right)m^2c^4`$ (1)
where $`\xi =\pm 1`$. It is worth briefly discussing the origin and reliability of the QG corrections to the dispersion relation in Eq.1. Ellis et al  have shown that quantum fluctuations in space-time foam lead to a energy dependent perturbed gravitational background metric

$`g_{ij}=\delta _{ij},g_{00}=1,g_{0i}=u_i/c`$ (2)
where $`|\stackrel{}{u}/c|E/E_{QG}`$, $`E_{QG}`$ being the scale at which the quantum gravity effects set in. Such a change in the metric violates Lorentz invariance at some high scale $`E_{QG}`$ and changes the dispersion relations of the particles. The linear energy dependence of such Lorentz invariance violating (LIV) terms in these models come mainly from the gravitational recoil effects. An identical form for the dispersion relation Eq.1 is obtained by studying massive spin-1/2 fields in the framework of loop quantum gravity where the LIV breaking correction term in Eq.1 arises due to the discrete structure of space-time at Planck scale .

In fact the last term in Eq.1 can arise in any theory which violates Lorentz invariance. The energy dependent form of the corrections may arise from a dimension 5 LIV operator of the form $`c_{\mu \nu \lambda }\overline{\psi }\gamma _\mu D_\nu D_\lambda \psi `$. The dimension 4 operator $`c_{\mu \nu }\overline{\psi }\gamma _\mu D_\nu \psi `$ would lead to an energy independent correction to the dispersion relation while the dimension 6 operator $`c\overline{\psi }\gamma _\mu D_\nu D^2\psi `$ would lead to only tiny energy independent corrections. Hence one may regard the dispersion relation in Eq.1 as a consequence of a generic type of LIV which may arise in some particular QG theories. From the point of view of measuring neutrino masses, the QG corrections to the dispersion relation represent a possible effect which may threaten to swamp the neutrino mass effect, and therefore it is important to include the largest imaginable such effects as we do here.

Assuming Eq.1 the time dependent neutrino velocity $`v(t)`$ is then given to lowest order in terms of the observed neutrino energy $`E_0`$ as 111Strictly $`E_0`$ is the neutrino energy at the time that the photons are observed on Earth $`E_0=E(t_0)`$, but since the neutrinos arrive a short time later this is to excellent approximation equal to the observed neutrino energy.

$`v(t)={\displaystyle \frac{E(t)}{p(t)}}c\left(1{\displaystyle \frac{ϵ_0^2}{2}}{\displaystyle \frac{a^2(t)}{a^2(t_0)}}+{\displaystyle \frac{3}{2}}\xi ϵ_0^{}{\displaystyle \frac{a(t_0)}{a(t)}}\right)`$ (3)
where we have expanded in $`ϵ_0^21`$ and $`ϵ_0^{}1`$ where

$`ϵ_0={\displaystyle \frac{mc^2}{E_0}},ϵ_0^{}={\displaystyle \frac{E_0}{E_{QG}}}.`$ (4)
The time-dependence of the neutrino velocity $`v(t)`$ arises from the expansion of the universe which redshifts the neutrino de Broglie wavelength $`\lambda (t)`$, and reduces their momentum and hence their velocity. The neutrino momentum $`p(t)`$ is related to the cosmic scale factor $`a(t)`$ through

$`{\displaystyle \frac{p(t)}{p(t_0)}}={\displaystyle \frac{\lambda (t_0)}{\lambda (t)}}={\displaystyle \frac{a(t_0)}{a(t)}}`$ (5)
Now suppose that low energy photons 222As discussed in the Introduction, the time delay is defined in terms of low energy photons for convenience only. What will matter in a practical search strategy is the comparison of different time delays between neutrinos of either the same energy, or between high energy neutrinos of two different energies. and high energy neutrinos are emitted from a GRB source at time $`t_e`$ and the low energy photons arrive on Earth at time $`t_0`$ while the neutrinos arrive at time $`t_\nu `$. The low energy photons will travel at the speed of light $`c`$, while the high energy neutrinos will travel with time dependent velocity $`v(t)`$ which may be smaller than $`c`$ due to the finite neutrino mass, and may be smaller or larger than $`c`$ due to the dispersive effects of quantum gravity. The co-moving distance $`\chi `$ between the GRB source and the Earth calculated in terms of the low energy photons is given by333While we have used he low energy photons for defining the co-moving distance, Eq.6 holds even for low energy massless neutrinos which travel at speed $`c`$.

$`\chi `$ $`=`$ $`{\displaystyle _{t_e}^{t_0}}{\displaystyle \frac{cdt}{a(t)}}`$ (6)
where $`a(t)`$ is the cosmic scale factor. The co-moving distance $`\chi `$ between the GRB source and the Earth calculated using the neutrinos which travel freely with time dependent velocity $`v(t)`$ is

$`\chi `$ $`=`$ $`{\displaystyle _{t_e}^{t_\nu }}{\displaystyle \frac{v(t)dt}{a(t)}}`$ (7)
$`=`$ $`{\displaystyle _{t_e}^{t_0}}{\displaystyle \frac{v(t)dt}{a(t)}}+{\displaystyle _{t_0}^{t_\nu }}{\displaystyle \frac{v(t)dt}{a(t)}}`$
The co-moving distance in Eq.7 must be equal to that calculated using the low energy photons or low energy massless neutrinos in Eq.6.

If we now assume that the neutrino velocity and the cosmic scale factor do not change much over the time scale $`t_0t_\nu `$, then equating Eq.7 to Eq.6 and using Eq.3 we find the time delay, to leading order in $`ϵ_0^2`$ and $`ϵ_0^{}`$,

$`\mathrm{\Delta }t=t_\nu t_0{\displaystyle \frac{ϵ_0^2}{2}}{\displaystyle _{t_e}^{t_0}}{\displaystyle \frac{a(t)}{a(t_0)}}𝑑t{\displaystyle \frac{3}{2}}\xi ϵ_0^{}{\displaystyle _{t_e}^{t_0}}{\displaystyle \frac{a^2(t_0)}{a^2(t)}}𝑑t`$ (8)
The time delay in Eq.8 may be expressed in terms of an integral over redshift defined as $`1+z^{}a(t_0)/a(t)`$, and Hubble constant defined as $`H\dot{a}/a`$,

$`\mathrm{\Delta }t{\displaystyle \frac{ϵ_0^2}{2}}I_2{\displaystyle \frac{3}{2}}\xi ϵ_0^{}I_1,`$ (9)
where $`I_n`$ are “redshift moments” of the inverse Hubble constant,

$`I_n={\displaystyle _0^z}{\displaystyle \frac{dz^{}}{(1+z^{})^nH(z^{})}}.`$ (10)
It is interesting to compare the expression for the time delay in Eq.9 to the result for the total time taken $`T`$ for the low energy photon or massless neutrino to travel from the GRB source to the Earth, 444The time $`T`$ represents how far back in time the GRB took place. As $`z\mathrm{}`$, $`TT_0`$ where $`T_0`$ is the age of the universe.

$`T=t_0t_e=I_1,`$ (11)
It is also interesting to compare to the result for the proper distance $`D`$ of the GRB source to the Earth at the photon arrival time, 555 $`D`$ represents the actual distance of the GRB from the Earth as measured now. As $`z\mathrm{}`$, $`DH_D`$ where $`H_D`$ is the horizon distance that light could have travelled since the beginning of the universe, which represents the size of the observable universe.

$`D=a(t_0)\chi =cI_0.`$ (12)
The upper limit of the integrals in Eq.10 $`z`$ represents the redshift of the GRB source, where $`1+z=a(t_0)/a(t_e)`$, and the Hubble constant $`H(z)`$ for a flat universe is given by

$`H^2(z)/H_0^2=\mathrm{\Omega }_M(1+z)^3+\mathrm{\Omega }_{DE}e^{I(z)}`$ (13)
where

$`I(z)=3{\displaystyle _0^z}(1+w(z^{})){\displaystyle \frac{dz^{}}{1+z^{}}}`$ (14)
where the equation of state for the dark energy $`w`$, defined as the ratio of its pressure to its density $`wp/\rho `$, in quintessence models has a redshift dependence. For constant $`w`$ we have

$`H^2(z)/H_0^2=\mathrm{\Omega }_M(1+z)^3+\mathrm{\Omega }_{DE}(1+z)^{3(1+w)}`$ (15)
and for a cosmological constant $`\mathrm{\Lambda }`$ with $`w=1`$ and $`\mathrm{\Omega }_{DE}=\mathrm{\Omega }_\mathrm{\Lambda }`$ this reduces to,

$`H^2(z)/H_0^2=\mathrm{\Omega }_M(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }.`$ (16)
The present day ratio of matter (M) density to critical density is $`\mathrm{\Omega }_M0.3`$, the present day ratio of dark energy (DE) density to critical density is $`\mathrm{\Omega }_{DE}0.7`$ and the Hubble constant is $`H_070\mathrm{k}\mathrm{m}\mathrm{s}^1\mathrm{Mpc}^1`$.

We shall consider three different cosmological models which were recently parametrised in : cosmological constant (CONST) with potential $`V\mathrm{\Lambda }^4`$; supergravity (SUGRA) inspired quintessence field (Q) model with potential $`V(Q)1/Q^{11}e^{Q^2/2}`$; inverse (INV) power law quintessence field (Q) model with potential $`V(Q)1/Q^6`$. The equation of state $`w(z)`$ for the SUGRA and INV models does not have a simple analytic form, but it may be parametrised as discussed in , and we shall use the parametrisation given there in this paper. Each of these models has $`w(z)`$ which varies with $`z`$ in such a way as to lead to a “tracking” behavior. This implies in particular that the value of the dark energy density remains of the same order as the matter density for large redshifts. By comparison the ratio of the cosmological constant energy density to the matter density falls as $`1/(1+z)^3`$ and rapidly becomes negligible at large redshifts.

In Figure 1 we calculate $`I_n`$ in Eq.10 for $`n=1,0,1,2`$ in units of the Hubble time $`t_H=1/H_014\mathrm{Gyr}`$ as a function of redshift $`z`$ for the three different cosmological models CONST, SUGRA, INV defined above. The physically relevant quantities $`D`$, $`T`$ and $`\mathrm{\Delta }t`$ are simply related to $`I_0`$, $`I_1`$, $`I_2`$ and $`I_1`$ as shown in Eq.12,11,9. Figure 1 shows that as $`n`$ increases $`I_n`$ becomes increasingly insensitive to redshift (flatter curves) and to the particular cosmological model (closer spacing between solid, dotted and dashed curves) at large redshift. Both effects can be simply understood as being due to the factor of $`(1+z^{})^n`$ in the denominator which tends to suppress the contributions from the higher redshift parts of the integration region for larger $`n`$. The difference between the curves for $`I_0`$ and those for $`I_1`$, $`I_2`$ in Figure 1 represents the error that would be made if the time delay were calculated using the naive formula

$`\mathrm{\Delta }t^{\mathrm{naive}}\left({\displaystyle \frac{ϵ_0^2}{2}}{\displaystyle \frac{3}{2}}\xi ϵ_0^{}\right){\displaystyle \frac{D}{c}}\left({\displaystyle \frac{ϵ_0^2}{2}}{\displaystyle \frac{3}{2}}\xi ϵ_0^{}\right)I_0`$ (17)
which ignores the effect of neutrino redshift as has been done up till now in the literature rather than the correct formula for $`\mathrm{\Delta }t`$ in Eq.9. For $`z=1`$ the error incurred by using the naive formula can clearly be seen to be of order 100%, with the error rapidly growing for larger redshifts. It is also clear from Figure 1 that there is some sensitivity to the nature of dark energy for large redshifts.

Finally it is interesting to consider the matter dominated limit of our time delay result corresponding to $`\mathrm{\Omega }_M=1`$, $`\mathrm{\Omega }_{DE}=0`$. In this limit the expression for the time delay in Eq.9 reduces to

$`\mathrm{\Delta }t^{\mathrm{matter}}{\displaystyle \frac{ϵ_0^2}{5H_0}}\left[1{\displaystyle \frac{1}{(1+z)^{5/2}}}\right]{\displaystyle \frac{3\xi ϵ_0^{}}{H_0}}\left[(1+z)^{1/2}1\right]`$ (18)
Time delays due to quantum gravity effects in the matter dominated limit were also considered in . However the result quoted there corresponds to the second term in Eq.17 rather than the second term in Eq.18 which correctly takes into account the expansion of the universe.

## 3 Results

The results in this section are based on the full formula for $`\mathrm{\Delta }t`$ in Eq.9, using the Hubble constant for a flat universe calculated using Eq.13 for the different dark energy models.

Figure 2 shows the time delays of the neutrinos against observed neutrino energy $`E_0`$ for a fixed GRB redshift of $`z=1`$ and assuming the cosmological constant model. In each panel the downward sloping dotted lines give the time delay due to the effect of neutrino mass which in the upper panels is chosen to be $`m=0.05`$ eV corresponding to the “atmospheric neutrino mass” defined as the square root of the atmospheric mass squared splitting , and in the lower panels we take $`m=0.005`$ eV corresponding to the “solar neutrino mass” which is the square root of the large mixing angle solar mass squared splitting . To probe neutrino mass corresponding to the solar scale the detector has to observe neutrinos with $`E_0`$ few 10 MeV with time sensitivity of a few milliseconds. The time delays corresponding to the atmospheric scale are higher and should be easier to detect in the planned neutrino telescopes.

The upward sloping dashed lines show the gravitationally induced time delay of the neutrinos. The upper panels give the time delay when the quantum gravity energy scale corresponds to the Planck scale. The lower panels are for the case where the quantum gravity corrections become important at $`10^{22}`$ GeV. The left panels have $`\xi =+1`$ corresponding to a negative time delay due to quantum gravity which tends to cancel with the positive time delay from the effect of neutrino mass. For this case the higher energy neutrinos arrive earlier. In the figure we have plotted the absolute value of the difference between their arrival times. The right panels have $`\xi =1`$ corresponding to a positive time delay due to quantum gravity which reinforces the positive time delay from the effect of neutrino mass. For both the cases we observe that the quantum gavity effects become very important for neutrinos arriving with energies greater than a few GeV. The ultra high energy neutrinos travelling cosmological distances can put severe constraints on $`E_{QG}`$. Note that the arguments in  suggest that fermions travel more slowly than $`c`$ and hence that $`\xi =1`$.

The solid lines in the figure show the time delay of the neutrinos due to the sum of the dotted and dashed lines when both the effects are present. For the case where $`\xi =+1`$ the two effects can cancel each other when from Eq.9 we have $`\mathrm{\Delta }t=0`$ which occurs at an energy

$`E_0^{}=\left({\displaystyle \frac{m^2E_{QG}}{3}}{\displaystyle \frac{I_2}{I_1}}\right)^{\frac{1}{3}}.`$ (19)
For $`\xi =1`$ even though the two effects reinforce each other, we see from Eq.9 and fig. 2 that the minima in $`\mathrm{\Delta }t`$ comes at exactly the same energy.

Figure 2 shows that time delay due to neutrino mass is important for the lower energy neutrinos with $`E_0<E_0^{}`$, while the time delay due to quantum gravity effects is important for the higher energy neutrinos with $`E_0>E_0^{}`$. This enables the two effects to be disentangled and treated separately.

### 3.1 Neutrino Mass Limits from Low Energy Neutrinos

In this section we shall focus on neutrino events with energy $`E_0<E_0^{}`$ in order to put limits on $`m`$. This implies that neutrinos with energies less than $`100`$ MeV can be effectively used to study their mass. Thermal neutrinos with energies between $`10100`$ MeV are expected to be produced in GRBs . Thermal neutrinos, coming from GRBs which are energetic enough, which are not very far away and which are probably beamed towards the Earth, should be detectable in the future km2 ice detectors like IceCube. These lowest energy neutrinos are expected to be detected mainly via $`\overline{\nu }_e+pn+e^+`$ where the positrons lead to an increase in the low PMT noise. This detection method which forms part of a Supernova Early Warning System may also be used to detect low energy neutrinos from GRB’s . However the determination of neutrino mass using the time delay techniques requires sensitivity to both arrival times as well as energy of the incoming neutrino. For a detector like IceCube since the inter-PMT distance is large the energy resolution is expected to be poor, particularly for the lower energy neutrino events. However energy resolution of the proposed megaton Super-Kamiokande type water Cerenkov detectors such as Hyper-K, UNO and TITAND should be good and can be effectively used for putting direct limits on neutrino mass. We next calculate the number of events that these detectors would observe from a GRB event.

One can estimate the number of thermal neutrinos emitted in a typical GRB by the following argument . Assume that the energy of the GRB in photons is $`E_\gamma `$ and the radius of the GRB fireball is $`R`$, then the photon temperature $`T_\gamma `$ is given by $`T_\gamma =(2E_\gamma /\sigma Vh_\gamma )^{1/4}`$, where $`V=\frac{4}{3}\pi R^3`$, $`\sigma `$ is the Stefan-Boltzmann constant and $`h_\gamma =2`$ is the number of degrees of freedom. Similarly the neutrino temperature $`T_\nu `$ is given by $`T_\nu =(2E_\nu /\sigma Vh_\nu )^{1/4}`$, where the number of degrees of freedom for the neutrinos is $`h_\nu =2\times 3\times \frac{7}{8}`$ and $`E_\nu `$ is the energy emitted in neutrinos. We shall assume $`E_\nu 100E_\gamma `$ . The total number of neutrinos emitted is given by $`N_\nu =E_\nu /E_\nu `$, where $`E_\nu =3.15T_\nu `$ is the average energy of the neutrinos, assuming a Fermi-Dirac distribution. From this we obtain

$`N_\nu =7\times 10^{57}\left({\displaystyle \frac{E_\nu }{10^{54}\mathrm{ergs}}}\right)^{3/4}\left({\displaystyle \frac{R}{100\mathrm{k}\mathrm{m}}}\right)^{3/4}.`$ (20)
The neutrino flux at Earth is

$`\mathrm{\Phi }_\nu ={\displaystyle \frac{N_\nu }{4\pi D^2}}{\displaystyle \frac{4\pi }{\mathrm{\Omega }_{\mathrm{beam}}}}`$ (21)
where $`D`$ is the distance of the GRB given by Eq.12 and $`\mathrm{\Omega }_{\mathrm{beam}}`$ is the beaming solid angle of the GRB corresponding to a cone of opening angle $`\theta _{\mathrm{beam}}`$, which we assume to be of order $`\theta _{\mathrm{beam}}1^{}`$ . The total number of events observed at the detector by $`\overline{\nu }_e`$ capture on protons is given by

$`N=N_p{\displaystyle \varphi (E)\sigma (E)𝑑E}`$ (22)
where the $`\overline{\nu }_e`$ flux $`\varphi (E)=\frac{1}{6}\mathrm{\Phi }_\nu f_{FD}(E)`$, $`f_{FD}(E)`$ being the Fermi-Dirac distribution function and we assume that the total neutrino flux $`\mathrm{\Phi }_\nu `$ is evenly distributed in all the six species. The reaction cross-section $`\sigma (E)`$ for $`\overline{\nu }_e+pn+e^+`$ is given by

$`\sigma (E)=0.94\times 10^{43}\left(EQ\right)\left((EQ)^2m_e^2\right)^{1/2}\mathrm{cm}^2`$ (23)
where $`Q=1.293`$ MeV is the neutron-proton mass difference and $`m_e`$ is the electron mass. $`N_p=\frac{2}{18}N_AM_D`$ is the number of free protons in water, where $`N_A`$ is the Avogadro’s number and $`M_D`$ is the fiducial mass of the detector. The integration in Eq.22 is done over observed neutrino energies. In Figure 3 we show the number of events expected in a one megaton water Cerenkov detector from thermal neutrinos produced in a GRB, which is beamed towards us with a beaming angle of $`1^{}`$ . The number of expected events are seen to be large for the GRBs with low redshift.

In Figure 4 we show the time delay $`\mathrm{\Delta }t`$ in seconds as a function of $`ϵ_0=mc^2/E_0`$ for redshift $`z=0.01,0.02,0.05,0.1,0.2,0.5,1,10`$, for the cosmological constant model. In this figure we assume that the neutrino energy satisfies $`E_0<E_0^{}`$ so that the effects of quantum gravity can be neglected and hence we can set $`ϵ_0^{}=0`$. For a fixed redshift the time delay increases as $`ϵ_0^2=m^2c^4/E_0^2`$ giving larger time delays for more massive neutrinos and lower observed energies, as seen in Figure 2. We are interested in measuring the smallest neutrino masses, which for a fixed $`ϵ_0`$ corresponds to the smallest observed energies $`E_0`$ and the smallest time delays. Since the neutrinos are released from the GRB in bursts on the time scale of milliseconds, the smallest time delays that are meaningful will also be milliseconds and from the figure this corresponds to $`ϵ_0=mc^2/E_010^{10}`$. For the lowest detectable energy neutrinos of $`E_010\mathrm{M}\mathrm{e}\mathrm{V}`$, $`ϵ_010^{10}`$ corresponds to neutrino masses of $`m10^3\mathrm{eV}`$. In order to discuss the limits on neutrino masses in more detail we must focus on the millisecond region of this plot.

In Figure 5 we show a blow-up of Figure 4 in the important millisecond region. We have also included the effect of different cosmological models in this figure, which can be important for high redshifts. From Figure 5 we see that time delays of $`10^3\mathrm{s}`$ correspond to $`ϵ_0=7\times 10^{10}8\times 10^{11}`$, over a range of redshift $`z=0.0110`$. For the lowest energy neutrinos of $`E_010\mathrm{M}\mathrm{e}\mathrm{V}`$ this range of $`ϵ_0`$ corresponds to a range of neutrino masses of $`m=7\times 10^38\times 10^4\mathrm{eV}`$. Although the smallest neutrino masses are associated with the highest redshifts the number of events expected in the detector falls sharply with increasing redshift. Since $`I_1`$ is not very sensitive to $`z`$ beyond $`z1`$, the fall in the number of events with distance is more acute than the rise in the time delay, and hence GRBs at lower values of redshift would be better suited for determination of the neutrino mass.

In addition to being sensitive to the absolute masses of all the three neutrino states, GRB neutrinos have the potential to probe the mass hierarchy as well. Though from solar neutrino data we know that the sign of $`\mathrm{\Delta }m_{21}^2\mathrm{\Delta }m_{}^2`$ is positive  ($`\mathrm{\Delta }m_{ij}^2=m_i^2m_j^2`$), there is still an ambiguity in the sign of $`\mathrm{\Delta }m_{32}^2\mathrm{\Delta }m_{atm}^2`$. It would be hard to determine the sign of $`\mathrm{\Delta }m_{32}^2`$ in any of the current and planned long baseline oscillation experiments and only a neutrino factory would be able to resolve this ambiguity. The delay in arrival times for the neutrinos can in principle be used to determine the neutrino mass hierarchy. Neutrinos arrive in three bunches corresponding to the three mass eigenstates. The neutrinos are detected via the electron antineutrino reaction, so for each mass eigenstate we must multiply the detection rate by the probability $`|U_{ei}|^2`$ that the mass eigenstate corresponding to the mass $`m_i`$ contains an electron antineutrino. Since $`|U_{e3}|^2`$ is small $`m_3`$ will have the smallest component of electron antineutrino and hence smallest number of events at the detector. In practice the states with mass $`m_1`$ and $`m_2`$ may be most easily detectable since $`|U_{e1}|^2`$ and $`|U_{e2}|^2`$ are not too small. The heaviest mass eigenstate will arrive last and have the largest time delay. In hierarchical models the heaviest mass is $`m_3`$ while in models with inverted mass hierarchy $`m_3`$ is the lightest. So depending on whether fewer events are detected for the neutrinos which arrive last (earliest) it would in principle be possible to conclude that the mass hierarchy is normal (inverted).

### 3.2 Quantum Gravity Limits from High Energy Neutrinos

The dispersion of velocity due to quantum gravity effects leading to delay in arrival times of the neutrinos is most significant at the higher energy end of the GRB neutrino spectrum, as evident from Figure 2. In fact for a given neutrino mass $`m`$, all neutrinos with energy $`E_0>E_0^{}`$ can be used to put limits on the scale of quantum gravity $`E_{QG}`$. Since neutrinos with a wide range of energies from $`110^9`$ GeV are expected to be produced in GRB fireball , the detection of their time delays is the most powerful method of studying/constraining models which predict dispersion relations given by Eq.1.

In this section we therefore focus on neutrinos with energy $`E_0>E_0^{}`$ in order to set limits on the quantum gravity scale. Figure 6 shows the time delay expected as a function of $`ϵ_0^{}=E_0/E_{QG}`$ for redshifts of 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1 and 10, assuming the cosmological constant model and also that $`E_0>E_0^{}`$ so that we can set $`ϵ_0=0`$. As also seen in Figure 2 the time delays in Figure 6 increase linearly with energy and hence with $`ϵ_0^{}`$. The time delay also increases significantly with the redshift $`z`$ as $`I_1`$ has a sharp $`z`$ dependence. Since we want to restrict ourselves to minimum time delays of a millisecond, which corresponds to the variability time of the GRBs, we show in Figure 7 the blow up of the region with $`\mathrm{\Delta }t10^3`$ seconds. We also show the effect of the different cosmological models on time delays on the same plot. For $`z=0.0110`$ the range of $`ϵ_0^{}`$ which would give $`\mathrm{\Delta }t10^3`$ seconds is $`ϵ_0^{}1.5\times 10^{19}2\times 10^{22}`$. For the highest energy neutrinos with $`E_010^9`$ GeV expected from the GRB, this would translate into $`E_{QG}6\times 10^{28}5\times 10^{30}`$ GeV. This should be compared with the bounds set by the time delays of the photons observed from GRBs  where the tightest limit obtained from the highest energy observed events is $`E_{QG}8.3\times 10^{16}`$ GeV . In the extreme ultra-relativistic regime the effect of neutrino mass is negligible and the dispersion of the neutrinos in the quantum space-time foam would be almost the same as that of photons, the latter having the advantage of being easier to detect. However ultra-high energy photons are subject to the GZK cut-off and beyond that for better bounds on $`E_{QG}`$ one can use ultra-high energy neutrinos. Bounds on $`E_{QG}`$ from non-observation of dispersion effects of the quantum space-time foam in neutrino oscillations experiments is slightly stronger with $`E_{QG}\stackrel{>}{}\mathrm{\hspace{0.25em}10}^{22}`$ GeV . This upper limit may be improved in the forthcoming long baseline experiments using neutrino superbeams. But it should be noted that this effect will show up in neutrino oscillations only if the dispersion due to quantum gravity is flavor dependent. However the method using time delays will be sensitive even if the effect is flavor independent. Thus the observation of time delays of GRB neutrinos is the most promising way of phenomenologically testing quantum gravity.

### 3.3 Dark Energy Limits from High Redshift Events

Different dark energy models give different values for the Hubble constant $`H(z)`$ at a given redshift. This can have significant impact on the time delays. Since the difference in $`H(z)`$ is maximum at the highest redshifts, we have to look at time delays of neutrinos coming from the farthest GRB events to probe dark energy. Since from Figure 1 we see that $`I_1`$ has a much better sensitivity to dark energy than $`I_2`$, time delays arising from quantum gravity would be a better probe for the different cosmological models. Hence it is the ultra high energy neutrinos which would be potentially more sensitive. In Figure 8 we show a further blow-up of the time delay plot due to quantum gravity effects, with $`\mathrm{\Delta }t`$ between $`5\times 10^310^2`$ seconds and $`ϵ_0^{}`$ between $`10^{21}10^{20}`$. We show the plots for $`z=`$1, 2, 4, 6 and 10. The three different line types correspond to the three different models for the dark energy that we have considered in this paper. The time delay obtained depends on the redshift as well as the cosmological model. For $`z>3`$ the time delays for different redshifts and different cosmologies begin to overlap. This ambiguity can be resolved by using smaller redshift events to determine $`E_{QG}`$, since for such events the effect of dark energy is reduced. Once $`E_{QG}`$ is determined then one can use the high redshift GRB neutrinos to probe dark energy. In order to successfully probe the correct cosmological model it is important that the redshift of the distant GRBs can be determined accurately enough from their afterglow and that the scale of $`E_{QG}`$ can be ascertained well enough by looking at the time delay for neutrinos from GRB from lower redshifts. If this can be achieved then information on $`\mathrm{\Delta }t`$ for neutrinos coming from the very distant sources can be used with information on $`z`$ and the scale of $`E_{QG}`$ to probe the correct cosmological model.

## 4 Discussion and Conclusion

In this paper we have calculated the time delays of neutrinos emitted in gamma ray bursts due to the effects of neutrino mass and quantum gravity. Our results are based on the formula for $`\mathrm{\Delta }t`$ in Eq.9, using the Hubble constant for a flat universe calculated using Eq.13 for different dark energy models. This formalism correctly takes into account the time dependence of the Hubble constant due to matter and dark energy, and can change the naive results in the literature by more than 100% for $`z>1`$. We have shown that the effects of neutrino mass, quantum gravity and dark energy may be disentangled by using low energy neutrinos to study neutrino mass, high energy neutrinos to study quantum gravity, and large redshifts to study dark energy. From low energy neutrinos one may obtain direct limits on neutrino masses of order $`10^3`$ eV, and distinguish a neutrino mass hierarchy from an inverted mass hierarchy. From ultra-high energy neutrinos the sensitivity to the scale of quantum gravity can be pushed up to $`E_{QG}5\times 10^{30}`$ GeV. By studying neutrinos from GRBs at large redshifts a cosmological constant could be distinguished from quintessence.

For convenience we have calculated all time delays with respect to a hypothetical low energy photon, assumed to be emitted at the same time from a point source as the neutrino of a given energy. We emphasise that what is important in practical search strategies is not these time delays themselves, which will be unmeasurable due to the uncertainties in the emission characteristics of low energy photons, but rather the comparison of time delays between neutrinos of either the same energy, or between high energy neutrinos of two different energies. As already mentioned, if the neutrinos are hierarchical in mass then neutrinos of the same energy will arrive at the detector in three bunches, and we can then compare the arrival times of the different neutrino mass eigenstates and put limits on the neutrino mass, as in the case of supernova neutrinos  but with sensitivities better by orders of magnitude. For constraining models of quantum gravity and dark energy one can use arrival times of the ultra-high energy neutrinos, and can compare the arrival times of the high energy neutrinos of different energies.

Acknowledgement

We particularly acknowledge very helpful discussions with Tom Marsh. We also wish to thank Juan Garcia-Bellido, Christian Kaiser, Ian McHardy, Tim Morris and Graham Shore for helpful discussions, and Francis Halzen for useful communications.


format: text

text:
# Numerical and Theoretical Studies of Noise Effects in the Kauffman Model

## 1 Introduction

The Kauffman model (or $`N`$-$`K`$ model) describes the dynamics of a network of $`N`$ Boolean spins, each controlled by $`K`$ other spins through a binary function. It was first proposed by S. A. Kauffman in 1969  as a model for cell differentiation and genetic networks. Since then, its application has been extended to many other fields in physics, biology, computational and social sciences. During the past few decades, most of the work done on Kauffman models has been dedicated to the study of the configuration space structure, the length and number of cycles, the size of basins of attraction, and the phase transition between ordered and disordered phases (for references see ). These properties are obtained by considering the deterministic dynamics of the system, which is well known by now. However, those studies have shown that some of the generic properties of the Kauffman model are far from being robust. The non-robustness of the deterministic dynamics is reflected, for example, in the fact that by slightly changing a given initial configuration of spins, the system may “jump” from one basin of attraction to a very different one. On the other hand, due to the exponential growth of the state space with $`N`$, it is often necessary to thoroughly probe the state space in order to determine a generic property of the system, such as the mean number of different basins of attraction or the mean cycle length. Actually, it has recently been shown that a systematic bias due to an under-sampling of the state space can be present in some of the results reported in the literature during the last 30 years . Therefore it is valuable to find a method which reveals the robust properties of the Kauffman model.

Real networks are always subjected to external fluctuations. Consequently, the relevant properties characterizing the network should exhibit a certain degree of robustness to external perturbations. In 1989, both Miranda *et al* and Golinelli *et al* analyzed the stochastic dynamics of the Kauffman model in the case in which an external noisy signal is present. In this work we extend the study of the stochastic dynamics of the Kauffman network with noise, by performing more accurate numerical simulations as well as analytic calculations. We focus our attention on the time it takes for two trajectories, starting out from different initial conditions, to cross. We consider two cases. First the situation in which each one of the $`N`$ spins is determined by $`K`$ other spins chosen randomly from everywhere in the system (the Kauffman net). The second case is a d-dimensional lattice in which each spin is preferentially coupled to its immediate neighbors. As we will see, both models exhibit qualitatively similar behavior.

In section 2, we introduce the Kauffman model with deterministic and stochastic dynamics. In section 3 we describe our numerical results for both Kauffman nets and lattices. In section 4, we study closely the behavior of these models in the limits of small and large noise. We summarize this work in section 5 with a brief discussion of the results.

## 2 The Kauffman model

### 2.1 Deterministic dynamics

A Kauffman model consists of $`N`$ Boolean spins $`\{S_1,S_2,\mathrm{},S_N\}`$ with $`S_i`$ being either zero or one. The value of each spin $`S_i`$ at time $`t+1`$ is determined by the values of $`K`$ other spins $`S_{i_1},S_{i_2},\mathrm{},S_{i_K}`$, which are called the controlling elements for spin $`S_i`$. (The number $`K`$ is called the *connectivity* of the system.) Once the connections in the system are established, each spin $`S_i`$ is assigned with a Boolean function $`f_i`$ of its $`K`$ controlling elements. A realization of the Kauffman model consists of the set of connections and Boolean functions assigned to every spin. The dynamics of the network is then given by

$$S_i(t+1)=f_i(S_{i_1}(t),S_{i_2}(t),\mathrm{},S_{i_K}(t))\text{for }i=1,\mathrm{},N.$$
(1)
For convenience, we will denote by $`\mathrm{\Sigma }_t`$ the state of the system at time $`t`$:

$$\mathrm{\Sigma }_t=\{S_1(t),S_2(t),\mathrm{},S_N(t)\}.$$

In different Kauffman models, the assignments of the $`K`$ controlling spins $`S_{i_1},S_{i_2},\mathrm{},S_{i_K}`$ of each spin $`S_i`$ and the dynamic rules $`f_i`$, are different. In *Kauffman nets*, the controlling elements of $`S_i`$ are assigned randomly, whereas in a *Kauffman lattice* they are chosen only among its nearest neighbors. The dynamic rules $`f_i`$ are chosen randomly in such a way that its two possible outcomes, 0 and 1, occur with probability $`\rho `$ and $`1\rho `$ respectively. If the realization of the network is time-independent, the network is called *quenched*, while if either the set of connections or the set of Boolean functions $`f_i`$ are re-assigned at every time step, the network is termed *annealed*.

Annealed models are more convenient for theoretical studies than quenched models. For example, by using the annealed approximation it has been shown analytically that Kauffman nets exhibit three different phases: frozen, critical and chaotic, depending upon the values of the parameters $`K`$ and $`\rho `$ . The critical value of the connectivity is given by $`K_c=[2\rho (1\rho )]^1`$. For $`K<K_c`$ the system is in the frozen phase, whereas if $`K>K_c`$ it is in the chaotic phase. Throughout this work we will use $`\rho =1/2`$, for which $`K_c=2`$.

But for most real cases (neural networks, genetic networks, etc.), quenched models will be more appropriate since in real networks neither the connections nor the interactions between the elements change randomly at every moment. However, it has been shown that in the limit $`N\mathrm{}`$, both the quenched and the annealed Kauffman nets are exactly equivalent with respect to the evolution of the overlap between different configurations, although not with respect to the configurations themselves . In this paper, our main focus is on quenched Kauffman nets and lattices.

Due to the finite size of the system, there are a finite number of possible configurations, to wit $`\mathrm{\Omega }=2^N`$. Therefore, starting out with an initial configuration, the system will eventually fall into a previously visited state, after which the same sequence of states repeatedly occurs again. The state space breaks up into a multitude of cycles (or attractors). The totality of points which end up in the same attractor represents its *basin of attraction*.

### 2.2 Stochastic dynamics

There are different ways of introducing noise into Kauffman models and they reveal different features of the configuration space of the model. Following Miranda and Parga , we introduced noise in the following way:

$$S_i(t+1)=\{\begin{array}{cc}f_i(S_{i_1}(t),S_{i_2}(t),\mathrm{},S_{i_K}(t))\hfill & \text{with probability }1r,\hfill \\ 1f_i(S_{i_1}(t),S_{i_2}(t),\mathrm{},S_{i_K}(t))\hfill & \text{with probability }r.\hfill \end{array}$$
(2)
In this way, every spin $`S_i`$ has a probability $`r`$ of violating the deterministic rule (1). We will say that an *$`n`$-spin flip event* has occurred at a particular time step, if $`n`$ spins violated the deterministic rule in this time step. Notice that this stochastic dynamic rule has a symmetry about the point $`r=0.5`$. For $`r>0.5`$, if we make a substitution $`f_i1f_i`$, the rule becomes identical with the case of $`1r`$. Since the Boolean functions $`f_i`$ are assigned randomly, $`f_i`$ and $`1f_i`$ are equally likely to appear in a particular realization of the model. After averaging over different realizations, the cases with probabilities $`r`$ and $`1r`$ are indeed identical. Due to this symmetry in the stochastic dynamical rule (2), we only need to consider the case $`r[0,0.5]`$. Note that for the particular value $`r=0.5`$, $`S_i(t+1)`$ is equally likely to be zero or one independently of the value of $`f_i`$.

In the presence of noise, the concept of “attractor” does not hold any more, since as the system evolves, there is a non-zero probability of “jumping” to a different attractor, and consequently every point in the state space can be reached from any initial condition. In this sense, the “boundaries” between different attractors become more diffuse as the level of noise increases . However, we shall argue that the system has a sort of effective attractor even in the presence of noise, specifically when $`r`$ is large.

One of the interesting things to study is the time it takes for two trajectories to cross. Suppose that we start with two different initial configurations, $`\mathrm{\Sigma }_0`$ and $`\stackrel{~}{\mathrm{\Sigma }}_0`$, and let them evolve according to (2), noting all the configurations produced:

$`\mathrm{\Sigma }_0`$ $``$ $`\mathrm{\Sigma }_1\mathrm{\Sigma }_2\mathrm{}\mathrm{\Sigma }_\tau `$
$`\stackrel{~}{\mathrm{\Sigma }}_0`$ $``$ $`\stackrel{~}{\mathrm{\Sigma }}_1\stackrel{~}{\mathrm{\Sigma }}_2\mathrm{}\stackrel{~}{\mathrm{\Sigma }}_\tau `$
The crossing time $`\tau `$ is then defined as the time for which either one of the trajectories coincides for the first time with a configuration previously visited by the other trajectory. For example, when $`\stackrel{~}{\mathrm{\Sigma }}_\tau `$ is equal to any of the configurations $`\mathrm{\Sigma }_0,\mathrm{\Sigma }_1,\mathrm{\Sigma }_2,\mathrm{},\mathrm{\Sigma }_\tau `$. Two important cases have to be distinguished, when $`\mathrm{\Sigma }_0`$ and $`\stackrel{~}{\mathrm{\Sigma }}_0`$ belong to the same basin of attraction, and when they belong to different basins of attraction. For those cases we will denote the crossing time by $`\tau _s`$ and $`\tau _d`$, respectively.

Miranda and Parga examined the behavior of the system by considering only the attractors with largest and next-largest basins of attraction. They then showed that for small values of $`r`$, the behavior of $`\tau _s`$ and $`\tau _d`$ are very different. At $`r=0`$, two trajectories from the same basin will cross in a time comparable with the sum of two times: first, the transient time required to enter the attractor, and second, the length of the attractor itself. Conversely, two trajectories starting out from different basins will never cross. On the other hand, they found that for sufficiently large values of $`r`$, the crossing time became independent of the starting point. It did not matter where the two trajectories start, the two basins merge into a sort of effective attractor and the trajectory bounces around within that subset of the system-states. For these larger values of $`r`$, the observed effective basin size increased with $`r`$. From this, they drew the conclusion that the disappearance of basins of attraction with the increase of $`r`$ is a sort of hierarchical process in the sense that in a finite period of time, the portion of the whole state space explored by a trajectory starting out from a given basin of attraction increases with $`r`$. Complete randomness is achieved at $`r=0.5`$, where the trajectory explores the entire state space.

As we will see, our simulation will show the same general behavior as described by Miranda and Parga. But we shall explore the behavior in more detail, showing the crossing time for the whole range of values of $`r`$ and $`K`$.

## 3 Numerical results

Kauffman nets and Kauffman lattices differ in the structure of their basins of attraction. One would expect this difference to be reflected in the response of these models to the influence of noise. For random realizations of the coupling functions $`f_i`$, what determines the basin structure is the connectivity $`K`$. Therefore, we will first analyze separately the cases with large $`K`$ (chaotic phase) and small $`K`$ (ordered phase).

We will partially follow Miranda and Parga’s approach in that we compute the average crossing time $`\tau _s`$ by using two initial configurations, $`\mathrm{\Sigma }_0`$ and $`\stackrel{~}{\mathrm{\Sigma }}_0`$, in the largest basin of attraction. For the average crossing time $`\tau _d`$, we pick one starting configuration in the largest basin and the other in the next largest one. The reasons to choose only the two largest basins of attraction will be clear in what follows.

### 3.1 Kauffman models with large $`K`$

We want first to characterize the structure of the basins of attraction. One way of doing it is by computing the distribution of basin-sizes $`W(n)`$, which is the fraction of the state space $`\mathrm{\Omega }`$ occupied by the $`n`$-th largest basin of attraction. In Fig. 1 we show $`W(n)`$ for a Kauffman net and a 1-dimensional Kauffman lattice, both with $`N=20`$ and $`K=5`$. It can be seen from this figure that both models exhibit a very similar structure in their basins of attraction in the sense that the basin sizes are similar. It is worth mentioning that in other aspects, like orbit length or transient time222The transient time is the time it takes before a trajectory enters the stable cycle., the basins of attraction can still be very different in both models.

From Fig. 1 it also can be seen that the largest and next largest basins occupy more that $`90\%`$ of the whole state space. Therefore, to a good approximation it can be assumed that the dynamics takes place mainly in these two largest basins.

Fig. 2 shows the average crossing times $`\tau _s`$ and $`\tau _d`$ as functions of $`r`$ for the Kauffman net and the $`1`$-dimensional lattice both with $`N=20`$ and $`K=5`$ (chaotic phase). Notice that these two kinds of $`N`$-$`K`$ models exhibit very similar behavior under the influence of noise.

### 3.2 Kauffman models with small $`K`$

The ordered phase is characterized by $`K=1`$ and $`K=2`$. In this section we will present the results for the minimum value of $`K`$, namely, $`K=1`$. In Fig. 3 we show $`W(n)`$ for a Kauffman net and a 1-dimensional Kauffman lattice, both with $`N=20`$ and $`K=1`$. The connections in the Kauffman net were, as usual, chosen randomly, whereas in the 1-dimensional lattice the node $`S_i`$ was connected either to $`S_{i1}`$ or to $`S_{i+1}`$ with equal probability (we use periodic boundary conditions).

From Fig. 3 it is apparent that in this case, the basin structures of the Kauffman net and the Kauffman lattice are less similar than in the chaotic phase. For the lattice, the two largest basins of attraction no longer occupy more than 90% of the whole state space, whereas in the Kauffman net they still do. However, the response to the influence of noise is mostly the same in both models, as can be seen from Fig. 4 where the crossing times $`\tau _s`$ and $`\tau _d`$ are plotted as functions of the noise intensity $`r`$.

### 3.3 Robust behavior of the crossing time

Figures 2 and 4 show that the behavior of the different Kauffman models under the influence of noise have the following general characteristics, both in the frozen and in the chaotic phases:

* For small $`r`$, $`\tau _d`$ decreases as $`1/r`$ while $`\tau _s`$ is nearly constant.
* For large $`r`$, both $`\tau _d`$ and $`\tau _s`$ increase with $`r`$ and become equal at $`r=0.5`$.
* For intermediate values of $`r`$, $`\tau _d`$ has a minimum when $`\tau _s\tau _d`$.

Let us analyze separately each one of the above characteristics.

In the limit $`r0`$, $`\tau _s`$ approaches a finite value $`\tau _s(0)`$, the mean crossing time for two trajectories in the largest basin of attraction in the absence of noise. This crossing time is roughly one half the average cycle length, plus one half the average transient time. The transient time and the cycle length are of the same order of magnitude, therefore $`\tau _s(0)`$ is expected to be approximately equal to the mean cycle length $`L`$ of the largest basin of attraction (see Fig. 5a).

On the other hand, $`\tau _d`$ diverges as $`r0`$. The numerical data (see figures 2b and 4b) suggest that in this limit, $`\tau _d`$ has the form

$$\tau _d\frac{a(K,N)}{r}+b(K,N),r0.$$
(3)
The above divergence is due to the fact that, in the absence of noise, there is a zero probability for a trajectory to jump between different attractors. Under the deterministic dynamics, every trajectory will remain within its own basin of attraction for ever. In the next section we will see that the $`1/r`$ behavior of $`\tau _d`$ is a consequence of the fact that the dynamics is governed by one-spin flip events when $`r`$ is small. For large values of $`K`$, the largest basin occupies almost the whole state space. Under these circumstances, every time a one-spin flip occurs, the trajectory in the next largest basin will have a finite probability of diverging very substantially from the path it would have followed in the absence of noise. That divergence will usually force the trajectory into the largest basin. In fact, for a fraction of order one of the noise events in the smaller basin, the noise will flip the trajectory into the largest one. Once the two trajectories are in the largest basin, they have a lifetime $`b(K,N)`$ before they cross. This lifetime is expected to be of order one of $`\tau _s(0)`$, the typical length of the largest basin of attraction. The above can actually be seen in Fig. 5, from which it is apparent that for large $`K`$, $`b(K,N)\tau _s(0)L`$.

In the opposite limit $`r0.5`$, the crossing times $`\tau _s`$ and $`\tau _d`$ become equal, which means that for high levels of noise, the barriers between different attractors become small. When $`r`$ reaches its maximum value $`0.5`$, all barriers vanish. In this case, both trajectories randomly jump from one state to another throughout the state space, and both $`\tau _s`$ and $`\tau _d`$ become equal to the time it takes for two random walks to cross. As derived in section 4.2, this crossing time is the solution to the “birthday problem”, i.e.

$$\tau _s=\tau _d2^{N/2}$$
(4)
In this way we have obtained a qualitative description of the limiting cases of figures 2 and 4. The one qualitative feature left to describe is the crossover from the small $`r`$ to the large $`r`$ behavior. As one can see from these figures, the crossover occurs when the two times $`\tau _s`$ and $`\tau _d`$ become roughly equal. This in turn happens when

$$ra(K,N)/\tau _s(0)1/N$$
(5)
Thus, the minimum in $`\tau _d`$ occurring between the two previous limit values of $`r`$, can be interpreted as the result of a sort of “competition” between the randomness in the system (coming from the presence of noise), and the barriers separating the attractors (which come from the deterministic dynamics of the system).

### 3.4 A Kauffman model with equal basin-sizes

Finally, we would like to mention that the above results are also true for Kauffman models in which all the basins of attraction have the same weight. As an example, consider a Kauffman lattice with $`N=20`$ and $`K=1`$ in which every spin is connected to itself. For $`K=1`$ there are only four Boolean functions $`f_i(S)`$: tautology $`f_i(S)=1`$, contradiction $`f_i(S)=0`$, identity $`f_i(S)=S`$ and negation $`f_i(S)=1S`$. Imagine then the very specific realization in which two of the coupling functions are identity, two are negation and all the others are either tautology or contradiction. By simple analysis, we know that for this specific model, the whole state space is composed of eight basins of attraction with equal size. Fig. 6 shows the crossing time $`\tau `$ as a function of $`r`$ for this particular model. Since in this case all the basins of attraction have the same size, the two initial conditions needed to compute $`\tau `$ were chosen randomly among the whole state space. Again, the $`\tau 1/r`$ behavior for small $`r`$ and the $`\tau 2^{N/2}`$ behavior for $`r0.5`$ are obtained. Of course we can construct many other models by choosing different coupling functions $`f_i`$ for this $`K=1`$ self-correlated case. All the Kauffman models we have explored have shown this kind of behavior under noise.

## 4 Theoretical analysis

### 4.1 Small $`r`$ limit

In this region, the main characteristic of average crossing time is that $`\tau _d1/r`$. The reason for this dependence is that the stochastic dynamics is dominated by one-spin flip events, in the following sense. The probability of a one-spin flip event ($`r`$) is much larger than the probability of a two-spin flip event ($`r^2`$). If the probability to jump to a different basin of attraction in a one-spin flip event is significantly different from zero, then the dynamics will be dominated only by this kind of events. Even if it was necessary to flip two spins to jump from one basin to another, this process can be decomposed into two one-spin flip events occurring sequentially, instead of being carried out at once in one two-spin flip event.

As we have shown in section 3, the $`1/r`$ behavior is present in a wide variety of Kauffman models. This in turn, implies that the one-spin flip events dominate the dynamics for small values of the noise. In this subsection, our goal is to derive an expression for the coefficient $`a(K,N)`$ for Kauffman nets. To do so, we will make the assumption that the dynamics takes place only in the two largest basins of attraction. Although the $`\tau _d1/r`$ behavior is generally true, as we have found, the preceding assumption is not true for all Kauffman models, especially for those with small values of $`K`$, but as we show below, it becomes more valid as $`K`$ increases.

To start the calculation of $`a(K,N)`$, let $`P_{1,2}`$ be the probability for jumping from the largest basin to the next largest basin with a one-spin flip event and $`P_{2,1}`$ be a similar probability but jumping in the opposite direction333By definition, the ratio of $`P_{1,2}`$ and $`P_{2,1}`$ is strictly the inverse of the ratio of the size of the largest basin to the next largest one in one realization.. Let us also define $`Q_1`$ and $`Q_2`$ as the probabilities of remaining in the largest basin and in the next largest one, respectively, after one-spin flip event. Simulations show that for Kauffman nets these probabilities have a slight dependence on $`N`$ but a very strong dependence on $`K`$. The result is that $`P_{1,2}+Q_1P_{2,1}+Q_2M`$, where $`M`$ is approximately constant for all values of $`K`$ and $`M>0.92`$ (see Fig. 7).

However, the dynamics depends not only on the total sum $`M`$ but also on the particular values of $`P_{1,2}`$, $`Q_1`$, $`P_{2,1}`$ and $`Q_2`$. For small values of $`K`$, both $`P_{1,2}`$ and $`P_{2,1}`$ are very small and comparable with $`(1M)`$. So, even though $`M`$ is large, the interaction between the two largest basins is weak in the sense that the probability of jumping into smaller basins is of the same order as $`P_{1,2}`$ and $`P_{2,1}`$. Therefore, for small values of $`K`$ the smaller basins play a significant role in the dynamics of the system. The above can be seen in Fig. 7, in which the probabilities $`P_{1,2}`$, $`P_{2,1}`$, $`Q_1`$, $`Q_2`$, and the sum $`M`$ are plotted as functions of $`K`$. From this figure it is apparent that the two largest basins of attraction are the dominant ones for large values of $`K`$ (say $`K5`$). The closeness of $`M`$ to $`1`$ means that a trajectory will seldom jump into a basin other than the two largest ones with only one-spin flip event. In view of this result, in some of the arguments below we will assume that $`K`$ is sufficiently large so that the dynamics takes place only in the two largest basins of attraction.

With the information about these probabilities, we can give an approximate calculation of $`\tau _d`$ as a function of $`r`$ for the case in which $`K`$ is large. We know that every spin violates the deterministic rule (1) with probability $`r`$. Therefore, the probability of a one-spin flip event is $`Nr`$ and consequently the expected time for this event to occur is $`T=1/(Nr)`$. This is true for all configurations. For sufficiently small values of $`r`$, this expectation time is much longer than the average crossing time for two configurations in the same basin. Hence, once two configurations jump into the same basin, their trajectories meet before the next spin-flip event becomes possible.

There are two cases in which the two trajectories meet after the occurrence of a one-spin flip event at time $`T`$: the configuration in the largest basin remains in it while the configuration in the next largest basin jumps into the largest one, or vice versa. The above occurs with probability $`(Q_1P_{2,1}+Q_2P_{1,2})`$. Similarly we get the probability for the crossing of the two trajectories after one-spin flip events at $`2T`$, $`3T`$, etc. The average value of $`\tau _d`$ is then:

$`\tau _d`$ $``$ $`T(Q_1P_{2,1}+Q_2P_{1,2})+2T(Q_1Q_2+P_{1,2}P_{2,1})(Q_1P_{2,1}+Q_2P_{1,2})+\mathrm{}`$ (6)
$`=`$ $`{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}mT(Q_1P_{2,1}+Q_2P_{1,2})(Q_1Q_2+P_{1,2}P_{2,1})^{m1}`$
$`=`$ $`{\displaystyle \frac{1}{Nr}}{\displaystyle \frac{Q_1P_{2,1}+Q_2P_{1,2}}{(1Q_1Q_2P_{1,2}P_{2,1})^2}}`$
$`=`$ $`{\displaystyle \frac{a(K,N)}{r}}`$
where $`a(K,N)`$ is explicitly given by

$$a(K,N)=\frac{1}{N}\frac{Q_1P_{2,1}+Q_2P_{1,2}}{(1Q_1Q_2P_{1,2}P_{2,1})^2}$$
(7)
In the above derivation, which is true for the case in which $`K`$ is large, we have ignored the time $`b(K,N)`$ it takes for two trajectories to cross after they have jumped into the same basin. Fig. 8 shows the coefficient $`a(K,N)`$ obtained from equation (7) and from simulation, for different values of $`K`$ in a Kauffman net with $`N=12`$. It can be seen that the simulation and the theoretical result agree very well for $`K5`$. It is worth emphasizing that for small values of $`K`$ and other Kauffman models where the two largest basins don’t have such dominance, the effect of other basins besides the largest and next largest ones has to be considered to perform an accurate derivation for the coefficient $`a(K,N)`$.

### 4.2 Large $`r`$ limit

When $`r`$ acquires its maximum value $`0.5`$, the barriers between different attractors vanish and the dynamics transforms into a random mapping of the state space into itself . In this limit, the crossing times $`\tau _s`$ and $`\tau _d`$ become indistinguishable and we will refer to both of them simply as $`\tau `$. To understand the limit $`r0.5`$ we will first give a simple “birthday-problem” argument to obtain the order of magnitude of $`\tau `$. Then we will proceed to a more elaborate analysis to obtain an approximate expression for $`\tau `$ as a function of $`r`$, valid for $`r`$ close to $`0.5`$.

Imagine two walkers moving at random through a space of dimension $`\mathrm{\Omega }=2^N`$. They go through a number of steps $`n=1,2,3\mathrm{}`$ As each step is completed, they have a total number of chances $`C_n`$ on landing on some place previously covered by the other walker. At the first step we could have the two walkers at identical positions. Thus, $`C_1=1`$. At the second step, $`C_2=3+C_1`$ since each can land on the original position of the other or they can both land at precisely the same place. After $`n`$ steps $`C_n=C_{n1}+2n+1=n^2`$. For large $`\mathrm{\Omega }`$ and small $`n`$, the probability of having collided with the path is then of order $`p_nC_n/\mathrm{\Omega }`$. The average time for crossing is then roughly given by the $`n`$-value for which $`p_n`$ becomes of order unity, so that $`n^2`$ is of order $`\mathrm{\Omega }`$ or $`n=O(\mathrm{\Omega }^{1/2})`$. Therefore, the crossing time $`\tau `$ satisfies $`\tau \mathrm{\Omega }^{1/2}`$.

To derive a more precise functional relation between $`\tau `$ and $`r`$ we have to compute the probability $`p_c(t)`$ for two trajectories $`\{\mathrm{\Sigma }_0,\mathrm{\Sigma }_1,\mathrm{},\mathrm{\Sigma }_\tau \}`$ and $`\{\stackrel{~}{\mathrm{\Sigma }}_0,\stackrel{~}{\mathrm{\Sigma }}_1,\mathrm{},\stackrel{~}{\mathrm{\Sigma }}_\tau \}`$ to cross at time $`\tau `$. In order to do that, we have to have first the probability $`p`$ for two configurations $`\mathrm{\Sigma }_t`$ and $`\stackrel{~}{\mathrm{\Sigma }}_{\stackrel{~}{t}}`$ to be the same (note that $`t`$ and $`\stackrel{~}{t}`$ might be different). Let $`S_i(t1)`$ be in $`\mathrm{\Sigma }_{t1}`$ and $`\stackrel{~}{S}_i(\stackrel{~}{t}1)`$ be the corresponding spin in $`\stackrel{~}{\mathrm{\Sigma }}_{\stackrel{~}{t}1}`$. Since these spins are in the same position (each in its respective configuration), the deterministic rule $`f_i`$ they obey is the same. Notice that the dynamical equation (2) can be written as

$$S_i(t+1)=\{\begin{array}{cc}f_i\hfill & \text{with probability }12r,\hfill \\ \text{evolve randomly}\hfill & \text{with probability }2r.\hfill \end{array}$$
(8)
From this expression it follows that the probability for $`S_i(t1)`$ and $`\stackrel{~}{S}_i(\stackrel{~}{t}1)`$ to evolve according to the deterministic rule $`f_i`$ is $`(12r)^2`$. Let us denote by $`p_1`$ the probability that $`S_i(t)=\stackrel{~}{S}_i(\stackrel{~}{t})`$ when both spins are updated according to the deterministic rule $`f_i`$. To calculate $`p_1`$ we follow the annealed approximation introduced by Derrida and Pomeau , which leads us to the following two possibilities:

1. The $`K`$ inputs of $`S_i(t1)`$ and $`\stackrel{~}{S}_i(\stackrel{~}{t}1)`$ are the same, which occurs with probability $`1/2^K`$. When this happens, $`S_i(t)=\stackrel{~}{S}_i(\stackrel{~}{t})`$ with probability $`1`$.
2. At least one of the inputs is different, which occurs with probability $`(11/2^K)`$. In this case, if the evolution rules $`f_i`$ are assigned in a sufficiently random way, there is a probability of 1/2 that $`S_i(t)=\stackrel{~}{S}_i(\stackrel{~}{t})`$.

From the above it follows that

$$p_1=(12r)^2\left[\frac{1}{2^K}+\left(1\frac{1}{2^K}\right)\frac{1}{2}\right]$$

On the other hand, the probability $`p_2`$ that $`S_i(t)=\stackrel{~}{S}_i(\stackrel{~}{t})`$ when the evolution rule $`f_i`$ is violated in one or both of the configurations, is simply given by

$$p_2=[1(12r)^2]\frac{1}{2}$$

Combining the values of $`p_1`$ and $`p_2`$ given above, the probability $`p`$ for both configurations $`\mathrm{\Sigma }_t`$ and $`\stackrel{~}{\mathrm{\Sigma }}_{\stackrel{~}{t}}`$ to be equal is

$`p`$ $`=`$ $`\left\{(12r)^2\left[{\displaystyle \frac{1}{2^K}}+(1{\displaystyle \frac{1}{2^K}}){\displaystyle \frac{1}{2}}\right]+[1(12r)^2]{\displaystyle \frac{1}{2}}\right\}^N`$ (9)
$`=`$ $`{\displaystyle \frac{1}{2^N}}\left\{1+{\displaystyle \frac{1}{2^K}}(12r)^2\right\}^N`$
If $`q(t)`$ is the probability that the two trajectories have not yet crossed at time $`t`$, then $`p_c(t)`$, the probability for the two trajectories to cross at time $`t`$, is given by

$$p_c(t)=q(t1)q(t)$$
(10)
The two trajectories are still separated at time $`t`$ if none of the configurations $`\{\mathrm{\Sigma }_0,\mathrm{\Sigma }_1,\mathrm{},\mathrm{\Sigma }_t\}`$ is equal to any of the configurations $`\{\stackrel{~}{\mathrm{\Sigma }}_0,\stackrel{~}{\mathrm{\Sigma }}_1,\mathrm{},\stackrel{~}{\mathrm{\Sigma }}_t\}`$. The probability for this to happen is

$$q(t)=(1p)^{(t+1)^2}$$

Substituting this value of $`q(t)`$ into equation (10) we get

$$p_c(t)=(1p)^{t^2}(1p)^{(t+1)^2}$$
(11)
Therefore, the average crossing time $`\tau =_{t=1}^{\mathrm{}}tp_c(t)`$ is given by

$`\tau `$ $`=`$ $`{\displaystyle \underset{t=1}{\overset{+\mathrm{}}{}}}t[(1p)^{t^2}(1p)^{(t+1)^2}]`$
$``$ $`{\displaystyle _0^+\mathrm{}}t{\displaystyle \frac{d(1p)^{t^2}}{dt}}dt`$
$`=`$ $`{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{\pi }{\mathrm{ln}(1p)}}}`$
Expanding the logarithm in the above equation around $`p=0`$, and retaining only the terms up to the the first order, we finally get

$$\tau \frac{\sqrt{\pi }}{2}\left[\frac{2}{1+\frac{1}{2^K}(12r)^2}\right]^{N/2}$$
(13)
It can be seen that Eqn. 13 is consistent with the “birthday-problem” argument for the case $`r=0.5`$. When $`r`$ is not exactly $`0.5`$, we do not have a simple birthday problem because the coupling between different elements, and consequently the functions $`f_i`$, still play a role in the dynamics. However, it is clear from the above equation that the problem can be viewed as a birthday problem with an effective state space $`\mathrm{\Omega }_{eff}=\left[\frac{2}{1+(12r)^2/2^K}\right]^N`$.

Fig. 9 compares the theoretical result (13) with the numerical simulation for $`K=1`$ and $`K=5`$. As can be seen, the analytic result approximates very well the numerical data in the region of $`r`$ close to $`0.5`$. In this region the annealed approximation holds because the noise breaks the correlations between the spins. However, for small values of the noise those correlations are important and cannot be neglected. We therefore do not expect agreement between the numerical and theoretical results for small values of $`r`$ since in this region the annealed appriximation is not longer valid.

## 5 Conclusions

We have considered the effect of external perturbations (noise) in the dynamics of the Kauffman model. The behavior of both, the Kauffman net and the Kauffman lattice under the influence of noise is very similar, even though these models might have a quite different structure in their basins of attraction. In this sense, the response of the Kauffman models to the effect of noise can be considered as a very robust property.

In the limit $`r0`$, the most important property is the $`1/r`$ behavior of the crossing time $`\tau _d`$, which has been always present in the Kauffman models we have studied so far. This $`1/r`$ behavior is a consequence of the fact that, for small values of the noise, the dynamics is dominated by one-spin flip events. An approximate equation relating $`\tau _d`$ and $`r`$ was obtained by taking into account the fact that, for large values of the connectivity, most of the dynamics takes place in the largest and next largest basins of attraction.

In the second limit $`r0.5`$, the barriers between different attractors disappear and the dynamics transforms into a random mapping of the state space into itself. As a consequence, $`\tau _s`$ and $`\tau _d`$ become equal. In the case in which $`r=0.5`$, the crossing time between the two trajectories can be seen as the solution of a “birthday problem” in a space of size $`\mathrm{\Omega }=2^N`$. For other values of $`r`$, but still close to $`0.5`$, the correlations between spins have to be taken into account, which have the effect of reducing the size of the region of the state space explored by the dynamics.

Between these two limit cases for the noise, there is a minimum in the value of $`\tau _d`$ as a function of $`r`$. In a loose sense, this minimum could be interpreted as the result of a “competition” between the randomness generated by the noise, which tends to homogenize the state space by diminishing the barriers across different attractors, and the deterministic dynamics, which tends to confine two trajectories within the same basin. To analyze this region it would be necessary to consider multiple-spin flip events as well as long-time step correlations.

The results and techniques presented in this work could be extended to other systems acting under the influence of noise in order to provide them with a robust characterization.

### Acknowledgements

This work was supported in part by the MRSEC Program of the National Science Foundation under award number 9808595, and by the NSF DMR 0094569. We also thank to the Santa Fe Institute of Complex Systems for partial support through the David and Lucile Packard Foundation Program in the Study of Robustness.


text:
# Origin Of Current-Induced Forces In An Atomic Gold Wire: A First Principles Study

## Abstract

We address the microscopic origin of the current-induced forces by analyzing results of first principles density functional calculations of atomic gold wires connected to two gold electrodes with different electrochemical potentials. We find that current induced forces are closely related to the chemical bonding, and arise from the rearrangement of bond charge due to the current flow. We explain the current induced bond weakening/strengthening by introducing bond charges decomposed into electrode components.

preprint: 1

The failure of interconnects in integrated circuits is often due to directionally biased diffusion of atoms caused by the presence of an electric current (electromigration)Sorbello (1998). In recent years the prospect of electronic devices operating essentially on the atomic/molecular scale has gained significant interestJoachim et al. (2000), and electromigration and current-induced conformational changes are important issues when downscaling electronic components to these sizesVentra et al. (2002). A single atom wide gold contact is among the simplest atomic scale conductors and therefore ideal for fundamental studies of the current-induced forces. It is remarkable that these systems can sustain voltages up to several volts before breakingYanson et al. (1998); Hansen et al. (2000); Nielsen et al. (2002).

Our aim here is to uncover the microscopic origin of the current-induced bond weakening/strengthening in an atomic gold wire. We analyze state-of-the-art first principles calculations of the current-induced forces while other investigations have focused on the theoretical foundation of these forcesTodorov et al. (2000). Traditionally the microscopic origin of the current induced force on an atom is related to electrostatic forces and the momentum transfer by the electron flow, the so called electron wind forceSorbello (1998). We find that the current-induced forces are rather linked to the redistribution of bond charge or overlap population in the system. By extending the concept of overlap population to a nonequilibrium system with two electrochemical potentials, we can rationalize the first principles results in terms of a simple two orbital interaction model.

We have performed our calculations within density functional theory (DFT) met ; Soler et al. (2002); Brandbyge et al. (2002). We employ the TranSIESTA programBrandbyge et al. (2002) which allows a full atomistic description of the scattering region and the electrodes. DFT has previously been used to describe the current-induced forces acting on single atomsLang (1992); Kobayashi et al. (1996), moleculesVentra et al. (2002), and an atom adsorbed on an nanotubeMingo et al. (2001). Recently, the current-induced forces in atomic wires, and the resulting “embrittlement”, was addressed by more approximate theoryTodorov et al. (2001).

We consider a simple symmetric geometry consisting of three atoms connecting electrodes in both the (100) and (111) directions, see Fig. 1. The electrode-electrode distance is in both cases chosen to be 9.6Å and the wire atoms are initially relaxed at zero bias. We choose this particular geometry since here the middle atom (2) can roughly be considered to have only two bonds namely to the left atom (1) and the right atom (3).

When the voltage is applied we observe a substantial redistribution of the electronic charge. At 0V the wires have a net negative charge of about 0.2 excess electrons. About half of this resides on atom 2. At 2V the excess electrons on atom 2 are almost halved. This charge redistribution results in an asymmetric voltage drop where the main drop is between the negative electrode and atom 2 (see Fig. 1).

In Fig. 2a we show the calculated bond forces for the (100)/(111) wires with atomic positions fixed at the relaxed valuesrel  obtained for zero bias. We find that the forces on the symmetric atoms 1 and 3 roughly follow each other, also beyond the linear bias regimeTodorov et al. (2000). We note that the order of magnitude and direction of the bond forces at 2V are similar for (100)/(111) while the detailed behavior differs: the bond forces for the (100) wire stay linear up to 1V and level off around 2V while we observe an onset of the bond forces for (111) at 1V. The magnitude of the bond forces of the order of 1-1.5nN for 1-2V is sizable compared to the force required to break a one-atom thick wire $`1.5`$nN Rubio-Bollinger et al. (2001).

Along with the bond force we also show the overlap population (OP) in Fig. 2a. The OP is a measure of the electronic charge residing in the bond. We note that except for the (111)-wire at low bias, the force and OP curves are almost mirrored. The result that an increase (decrease) in bond charge makes the bond force compressive (repulsive) is what we expect intuitively. In Fig. 2b we show the bond lengths after subsequent relaxation of the wire atoms for finite bias. The main relaxations take place in the wire bonds with smaller displacements in the wire-electrode bonds.

We aim at an qualitative understanding of the results for the dependence of the force on the voltage bias. To this end we focus on the OP. We can expressfor  the force acting on atom $`i`$ due to the valence electrons using the force operator, $`𝐅_i`$, and the density matrix $`𝐃`$,

$$\stackrel{}{F}_i=\text{Tr}[\stackrel{}{𝐅}_i𝐃]\text{where},\stackrel{}{𝐅}_i=\frac{𝐇}{\stackrel{}{R}_i}.$$
(1)
Consider the bonding between two atoms represented by orbitals $`|\varphi _1`$ and $`|\varphi _2`$ and separated by bond length $`b`$. From Eq. 1 we get the bond force,

$$F_{\text{bond}}=2\left(\frac{\varphi _2|H^{}|\varphi _1}{\varphi _2|\varphi _1}\right)O_{12},H^{}=\frac{H}{b},$$
(2)
where we assume that only the hopping element, $`\varphi _2|H|\varphi _1`$, changes with $`b`$. The bond force is proportional to the OP for the 12 bond,

$$O_{12}=2𝐒_{12}𝐃_{12}.$$
(3)
The OP is typically taken as a simple measure of the strength of a chemical bond as suggested by Mulliken for moleculesMulliken (1955). For extended systems the density matrix can be expressed in terms of the density of states (DOS) matrix, $`𝝆`$. In equilibrium we have,

$$O_{12}=2𝐒_{12}_{\mathrm{}}^{\mathrm{}}𝑑\epsilon 𝝆_{12}(\epsilon )n_F(\epsilon E_F),$$
(4)
where $`n_F`$ and $`E_F`$ is the Fermi function and energy. The contribution from states at different energy($`\epsilon `$) in the extended system to the bonding between the two orbitals is described qualitatively by the OP weighted DOS (OPWDOS) or COOP curve, $`2𝐒_{12}𝝆_{12}(\epsilon )`$, as discussed by HoffmannHoffmann (1988).

Now we consider the nonequilibrium situation where an electrical current is running though a contact connecting two reservoirs, left and right($`L`$,$`R`$), with different chemical potentials($`\mu _L`$, $`\mu _R`$). The contact contains the 12 bond. In this case the bond force will change since the density and Hamiltonian matrices deviate from equilibrium values. The appropriate density matrix for the nonequilibrium situation is constructed from scattering states, and the total spectral density matrix($`𝝆`$) is split into partial contributions corresponding to their left($`𝝆^L`$) or right($`𝝆^R`$) origin (see e.g. Brandbyge et al. (2002)),

$$𝝆(\epsilon )=𝝆^L(\epsilon )+𝝆^R(\epsilon ).$$
(5)
The density matrix for nonequilibrium is then found by filling the left and right originating states according to the respective chemical potentialsTodorov et al. (2000); Brandbyge et al. (2002),

$$𝐃=_{\mathrm{}}^{\mathrm{}}𝑑\epsilon 𝝆^L(\epsilon )n_F(\epsilon \mu _L)+𝝆^R(\epsilon )n_F(\epsilon \mu _R).$$
(6)
We split the overlap population(Eq. 4) further into its partial left and right components,

$$O_{12}=O_{12}^L+O_{12}^R,$$
(7)
where

$$O_{12}^L=2𝐒_{12}_{\mathrm{}}^{\mathrm{}}𝑑\epsilon 𝝆_{12}^L(\epsilon )n_F(\epsilon \mu _L),$$
(8)
and likewise for $`R`$. When more orbitals are involved on atom 1 and 2 the overlap populations can be broken down in orbital components. In this case each orbital overlap will contribute with different weight to the bond force according to their different force matrix elements $`H^{}`$. This complication, together with the neglect of change in $`H^{}`$ with applied voltage, makes the bond force and OP only roughly proportional, as we observe for the gold wires in Fig. 2. However, as a first approximation, we can explain the observed change in bond forces from the change in OP, which in turn means the change in left ($`\mathrm{\Delta }O^L`$) and right ($`\mathrm{\Delta }O^R`$) contributions.

In Fig. 3 we consider the change from 0V to 1V in left and right OP of the 12 and 23 bonds for the (100) case (we see a quite similar pattern for (111)). For both bonds $`\mathrm{\Delta }O^L`$ is negative (decreasing bonding) while $`\mathrm{\Delta }O^R`$ is positive (increasing bonding). We have resolved these further into the orbitals on atom 2 and display the main contributions in Fig. 3.

In order to understand the physics involved we present in Fig. 4a a simple picture. In the presence of current the electrons from the left electrode populate states above the original Fermi energy, $`E_F`$. These populate bonding or anti-bonding states depending on the position of $`E_F`$. Likewise there is an depopulation in electrons from the right. To make this picture more explicit we consider in Fig. 4b a single bond between two atoms each with a single atomic level, $`\epsilon _L`$ and $`\epsilon _R`$. The levels couple via a matrix element $`t<0`$, and an overlap $`S>0`$, and form bonding/anti-bonding orbitals, $`\epsilon _B`$ and $`\epsilon _A`$. The left and right atoms are coupled to left and right electrodes modelled by a level broadening $`\mathrm{\Gamma }_L`$ of the left level and $`\mathrm{\Gamma }_R`$ for the right. For this bond the resulting partial DOS due to the states originating from the left is,

$$\rho ^L(\epsilon )=\frac{2}{\pi }\frac{\mathrm{\Gamma }_L\left(\epsilon _R\epsilon \right)(\epsilon S+|t|)}{\left[(1S^2)(\epsilon \epsilon _A)(\epsilon \epsilon _B)\mathrm{\Gamma }_L\mathrm{\Gamma }_R\right]^2+\left[\mathrm{\Gamma }_R\left(\epsilon \epsilon _L\right)+\mathrm{\Gamma }_L\left(\epsilon \epsilon _R\right)\right]^2}.$$
(9)
We note: (i) $`\rho ^R`$ is found by exchanging $`LR`$. (ii) $`\rho ^L`$ is proportional to the coupling $`\mathrm{\Gamma }_L`$. (iii) The cross-over between bonding(positive) and anti-bonding(negative) is given by the position of $`\epsilon _R`$. From (i), (iii) we see that $`\rho _L`$ contains more anti-bonding character than $`\rho _R`$ when $`\epsilon _L>\epsilon _R`$.

We have plotted $`L`$/$`R`$ COOP curves in Fig. 4(c-e) for different fillings and $`\mathrm{\Gamma }_L=\mathrm{\Gamma }_R`$. In the top row we assume that there is no voltage drop between the atoms, $`\epsilon _L=\epsilon _R`$. Here the change in OP comes entirely from the change in filling due to the shift of $`\mu _L`$ and $`\mu _R`$ away from the equilibrium $`E_F`$ (Eq. 8). For almost filled states (Fig. 4c, top) the lowering of $`\mu _R`$ deplete states with anti-bonding character and $`O^R`$ increases. The increase in $`\mu _L`$ leads to a decrease in $`O^L`$ which does not balance the increase in $`O^R`$. The net effect is a decrease in anti-bonding character and strengthening of the bond. When the states are almost empty (Fig. 4e, top) we also observe a strengthening since now $`O^L`$ increases more than $`O^R`$ is decreasing. For a bond initially with maximal strength (Fig. 4d, top) and, in general, for $`E_F`$ located in a region where the COOP decreases, the bond weakens since both $`O^R`$ and $`O^L`$ decrease.

In the bottom row of Fig. 4(c-e) we assume that the full voltage drop takes place between the two atoms, $`\epsilon _L\epsilon _R=\mu _L\mu _R>0`$. The bonding peak becomes dominating in the right COOP and likewise for the anti-bonding peak in the left COOP. The net result is less bond strengthening in (c) and (e), and more bond weakening in (d). The largest effect of a voltage drop is seen for half filling (Fig. 4d, bottom) where a large negative $`\mathrm{\Delta }O_L`$ occurs because the bonding peak decreases and at the same time more anti-bonding states are being occupied.

In summary, for equal coupling, $`\mathrm{\Gamma }_L=\mathrm{\Gamma }_R`$, our simple model yield: (i) bond strengthening for almost empty/filled states, (ii) bond weakening for about half filled states, and (iii) a bond with a voltage drop is weaker than it would be without a drop. When the bond is not coupled equally well to left and right, $`\mathrm{\Gamma }_L\mathrm{\Gamma }_R`$, the left and right overlap contributions have to be weighted accordingly.

In the case of the gold wires we can to a good approximation assume that the $`d_{yz}`$, $`d_{zx}`$, and $`p_z`$ orbital on atom 2 only couple with the corresponding orbital on atom 1 and 3. The $`d_{yz},d_{zx}`$ states are almost filled and match the case in Fig. 4c: $`\mathrm{\Delta }O^R>0`$ while $`\mathrm{\Delta }O^L<0`$. On the other hand the $`p_z`$ states are almost empty and match the situation in Fig. 4e: $`\mathrm{\Delta }O^L>0`$ while $`\mathrm{\Delta }O^R<0`$. This is exactly the behavior we observe in Fig. 3. For the $`s`$ orbitals, which are slightly more than half-filled, we can expect a decrease in $`O^L`$, which will be especially pronounced if a voltage drop takes place in the bond (cf. Fig 4d, bottom). This is also what we observe for the 12 bond in Fig. 3.

Based on our simple model we see that the almost completely filled states($`d`$) and slightly more than half-filled states($`s`$) together yield a decrease in $`O^L`$ and increase in $`O^R`$. Since bond 12 couple more to the left electrode the total change in $`O^L`$ will be bigger than the change in $`O^R`$. This explains the bond weakening of bond 12. The opposite is the case for bond 23. Also the voltage drop taking place in bond 12 will weaken this bond with respect to bond 23.

In conclusion we have analyzed first principles calculations of the current-induced forces in an atomic gold wire. We have shown that the current-induced forces are due to a shift in the population of bonding and antibonding levels, and can be explained by examining the change in electrode decomposed bond charges(overlap population) with applied voltage. Although we have concentrated on a simple system the approach should apply to atomic scale conductors in general.

###### Acknowledgements.

We acknowledge support from the Danish Research Councils (M.B. and K.S.), European Union (SATURN project) and MCyT (J.L.M. and P.O). We also thank the $`\mathrm{\Psi }_k`$ network and CEPBA for funding collaborative visits.