text:
# Dynamics of Turing patterns under spatio-temporal forcing

## Abstract

We study, both theoretically and experimentally, the dynamical response of Turing patterns to a spatio-temporal forcing in the form of a travelling wave modulation of a control parameter. We show that from strictly spatial resonance, it is possible to induce new, generic dynamical behaviors, including temporally-modulated travelling waves and localized travelling soliton-like solutions. The latter make contact with the soliton solutions of P. Coullet Phys. Rev. Lett. 56, 724 (1986) and provide a general framework which includes them. The stability diagram for the different propagating modes in the Lengyel-Epstein model is determined numerically. Direct observations of the predicted solutions in experiments carried out with light modulations in the photosensitive CDIMA reaction are also reported.

The study of pattern dynamics under external forcing provides a powerful tool to deeply probe their inherently nonlinear mechanisms under non-equilibrium conditions. A great deal of attention has been focused on resonances or locking of spatially structured states, either stationary or oscillatory, under temporal (spatially uniform) 1 ; 2 ; 3 ; 4 ; 5  or spatial (steady) modulations 6 ; 7 ; 8 ; 9 ; 10 . Steady patterns in reaction-diffusion systems typically arise from the celebrated Turing mechanism 11 ; 12 . According to it, inhomogeneous distributions of chemical concentrations self-organize spontaneously out of a non-structured medium as a result of a competition between autocatalytic reaction steps and the differential diffusivities of an activator (smaller) and an inhibitor-like (larger diffusion) species. Turing patterns are endowed with an intrinsic wavelength, depending only on the kinetic and diffusion parameters, but lack an intrinsic frequency, in contrast to oscillatory chemical systems 12 . Genuine Turing patterns were first observed in quasi two-dimensional gel reactors (pre-loaded with appropriate chemical indicators) in the CIMA 13  and the CDIMA 16  reactions, and appear as patterned distributions of iodide. The CDIMA reaction has the interesting feature of being photosensitive 14 .

It seems thus timely to search for generic behavior in the unexplored perspective of spatio-temporal forcing of pattern forming systems. Specifically, we aim at studying the dynamical response of photosensitive Turing patterns to the simplest external spatio-temporal forcing consisting of a travelling-wave modulation of the control parameter associated to the illumination. Through the mechanism of pure spatial resonance an external frequency will thus be imposed in an otherwise nonoscillatory system. As a consequence, new nontrivial dynamical modes are expected to arise which allow to connect the two trivial limiting cases, namely: a travelling pattern locked to the forcing at low velocities and a standing pattern resulting from the time-averaging of the forcing at large velocities. Analytical and numerical results will be reported featuring the simplest of these spatio-temporal behaviors. Experiments conducted with the CDIMA reaction are also provided which fully confirm our theoretical predictions. Beyond the particular chemical context that motivates our study, such solutions, and the conditions of their appearance, are sufficiently generic to be applicable to a rather general class of pattern forming systems including for instance Rayleigh-Bénard convection. In a sense, our study may be viewed as a development of the work of Coullet 7 ; 15  on commensurate-incommensurate transitions in nonequilibrium (spatially) forced systems.

The experimental system under study is modelled within the Lengyel-Epstein scheme 16 , once modified to include the effect of illumination, as 14 :

$`_tu`$ $`=`$ $`acu4{\displaystyle \frac{uv}{1+u^2}}\varphi +{\displaystyle \frac{^2u}{x^2}},`$ (1)
$`_tv`$ $`=`$ $`\sigma (cu{\displaystyle \frac{uv}{1+u^2}}+\varphi +d{\displaystyle \frac{^2v}{x^2}}).`$ (2)
Here $`u`$ and $`v`$ are the dimensionless concentrations of two of the chemical species; $`a`$, $`c`$, $`\sigma `$, and $`d`$ denote dimensionless parameters of the chemical system. The effect of external illumination is introduced through the $`\varphi `$-terms. This contribution can be decomposed into the mean value $`\varphi _0`$ and a modulation part: $`\varphi (x,t)=\varphi _0+ϵ\mathrm{cos}(k_fx+\omega t)`$. For purely homogeneous illumination, $`ϵ=0`$, the equations admit a solution which in the following will be referred to as base state: $`u_0=(a5\varphi _0)/(5c)`$, $`v_0=a(1+u_0^2)/(5u_0)`$.

All our numerical results are obtained through integration of the model reaction-diffusion equations (1),(2) with periodic boundary conditions by means of a pseudo-spectral method with a linear-implicit time-stepping. From here on we fix the parameters to the following: $`a=16,c=0.6,d=1.07,\text{and}\sigma =301`$. These values were chosen to reproduce the experimental conditions referred to below. The large $`\sigma `$ value, corresponding to a strong diffusion contrast between the two species, guarantees that we are far from the oscillatory regime of this chemical system. Consequently, only the Turing bifurcation will play a role. The remaining parameters concern the forcing term. One of them, $`\varphi _0`$, will serve as the parameter to locate the position of the Turing bifurcation in the homogeneous problem. For the given parameters it occurs at $`\varphi _02.3`$ (the base state being unstable to Turing patterns below $`\varphi _02.3`$) and the critical wavenumber is $`k_c1.07`$.

We begin our analysis with the case of exact 1:1 spatial resonance $`k_c=k_f`$. Choosing the length $`L`$ of the periodic domain to be $`10\lambda `$, where $`\lambda =2\pi /k_c`$ is the critical wavelength at instability, we fix $`L=58.72`$ and the amplitude $`ϵ=0.1`$.

As is known from the analysis of time-independent forcing, the spatial 1:1 resonance yields an imperfect bifurcation to Turing patterns 7 . Accordingly, the base state ceases to be a stationary solution and is modified into a non-homogeneous state for every value of the bifurcation parameter $`\varphi _0`$. For $`\varphi _0`$ in the stable region (above 2.3) this state is a travelling wave (TW). The TW locks to the forcing wave, adopting the same wave-number and frequency, with only a constant phase difference. This is the trivial state to be expected for slow driving. For large $`\omega `$ the amplitude of the TW approaches 0. The TW’s exist to the right of the solid curve in fig. 1.

Crossing the solid curve in fig. 1 the TW state undergoes a first instability into a state with temporally modulated amplitude (MTW, triangles in fig. 1). This is the signature of a Hopf bifurcation which introduces a new frequency $`f_H`$ (see fig. 2a). Note that the modulation of the amplitude occurs uniformly in the entire system. We found that for large $`\omega `$ the Hopf frequency converges to the frequency $`\omega `$ of the forcing wave.

The solutions and the transition described above can be rationalized in terms of an amplitude equation. Following standard envelope techniques near threshold 7 ; Cross , and with the forcing being invariant under the transformation $`tt+T`$, $`xx\omega T/k_f`$, for sufficiently small $`ϵ`$ and $`\omega `$, the slowly varying modulations of the travelling mode, in the case of perfect 1:1 resonance ($`k_f=k_c`$) will be given by the amplitude equation

$$\dot{A}=\mu A|A|^2A+ϵ\mathrm{exp}(i\omega t)+^2A/x^2.$$
(3)
Using polar coordinates, $`A=R\mathrm{exp}i\mathrm{\Theta }`$, we look for homogeneous solutions with $`\mathrm{\Theta }=\mathrm{\Theta }_0\omega t`$. As in the steady case for $`ϵ0`$ there is a non-zero solution for every $`\mu `$, the dimensionless distance to threshold. Its amplitude approaches $`0`$ with increasing $`\omega `$ in accordance with our observations for the Lengyel-Epstein model. $`\mathrm{\Theta }_0`$ is the phase shift between the forcing wave and the resulting pattern.

We further determine the stability of this solution with respect to homogeneous perturbations. Directly from the amplitude equation:

$`\dot{R}`$ $`=`$ $`\mu RR^3+ϵ\text{cos}(\psi +\mathrm{\Theta }_0),`$ (4)
$`\dot{\psi }`$ $`=`$ $`{\displaystyle \frac{ϵ}{R}}\text{sin}(\psi +\mathrm{\Theta }_0)+\omega ,`$ (5)
where we have defined: $`\psi =\mathrm{\Theta }+\omega t\mathrm{\Theta }_0`$. Linearization about the locked solution ($`\psi =0`$, $`R=Q=const`$) yields the following eigenvalues: $`\lambda =\mu 2Q^2\pm \sqrt{}(Q^4\omega ^2)`$. Corresponding to fig. 1, for large $`\omega `$, the marginal curve approaches the line $`\mu =0`$, and the imaginary part of the eigenvalue, converges to the driving frequency, as also observed numerically above footnote2 .

We now address the more generic case of inexact 1:1 resonance footnote3 , introducing a slight wavenumber misfit, $`k_fk_c`$. To allow for a continuous variation of the misfit in a finite system, we will fix the integer wavenumber ($`n=kL/2\pi `$) of the forcing to $`n=10`$ and change smoothly the length $`L`$ of the periodic domain. For example, for $`L=65`$ the 10th wave number corresponds to $`k_f0.97`$. Fig. 3 depicts the complex behavior that was found changing $`k_f`$ from 0.9 to 1.26 ($`L=70,\mathrm{},50`$). The average illumination $`\varphi _0`$ was fixed to 2.25 during all of the simulations. For the purely homogeneous forcing (i.e. $`\varphi ^{}=0`$) this value corresponds to a slightly unstable base state.

In fig. 3 the stability domain in $`\omega `$ and $`k_f`$ is given for four different states. The TW state, the solution locked to the forcing, is the only stable solution for approximately $`\omega 1`$. Above a roughly horizontal transition curve close to $`\omega =1.1`$ the TW states are unstable. For $`k_f=k_c`$ and close to this point the transition occurs at $`\omega 1.15`$ (cf. fig. 1) and results in the MTW.

The domain of stability of the MTW is given by the solid curve in fig. 3. Outside of this domain the attractors are either the TW state or one of two different new localized states which we call soliton waves (SOW)footnote1 . A typical space-time plot is shown in fig. 2b. Spatial plots of these states show that they resemble the soliton-like solutions for the case of non-travelling forcing 7  but they exhibit the following new properties: the soliton, i.e., the localized suppression of the amplitude moves relatively to the underlying travelling pattern with a relative velocity which may be positive or negative. Furthermore, as the soliton travels along the pattern the phase of the background pattern moves in either direction. In our case the integer wavenumber of the background pattern is either $`n=9`$ or $`n=11`$ and the corresponding states are represented in fig. 3 with circles and diamonds, respectively.

Beyond the range of forcing wavenumbers shown in fig. 3 there exist further soliton states in accordance with the prediction for non-travelling forcing 7 . Here we will describe only the two states that are adjacent to the MTW state in fig. 3. The states represented by diamonds (resp. circles) carry a soliton that moves to the right (resp. left). The approximate range of stability for these states is given by the dashed curve (resp. dotted curve).

Note that the stability domains of the soliton states and the MTW state overlap. Remarkably, at sufficiently large $`\omega `$ the soliton states are even stable for $`k_f=k_c`$, i.e., exact 1:1 resonance. This contrasts to the case of static forcing 7  where only the presence of a misfit can lead to solitons.

The above localized states (SOW) can be described as quasi-periodic in the sense that their behavior can be reproduced as superposition of two modes, namely, the forcing (travelling) mode with wavelength $`L/10`$ and an adjacent mode with wavelength $`L/9`$ or $`L/11`$, typically the closer one to the characteristic wavelength of the Turing instability. The excited Turing mode is standing so that the superposition of both yields a localized envelope (soliton) moving to the right ($`L/11`$) or to the left ($`L/9`$). Furthermore, one finds the velocity of the soliton to be proportional to the quotient of $`\omega `$ and the difference of the wavenumbers of the two modes, which agrees with numerical results of the full problem.

Note that the phase diagram depicted in fig. 3 exhibits large regions of multistability, with the corresponding hysteresis. It is worth remarking that the stability boundaries here reported may well be affected by to finite-size effects, since possible long wave-length, phase instabilities may be prevented by the finite size of our simulations. A full envelope and phase-diffusion description of this scenario in an infinite system deserves a detailed study and will be addressed elsewhere.

A link between the SOW’s and the solitons for steady forcing is apparent for smaller $`ϵ`$. We decreased the value of $`ϵ`$ to 0.003 and determined the domain of stability for the soliton 11 state for $`\omega `$ values close to 0 (fig. 4). It appears that even for vanishing $`\omega `$ the soliton states persist and the domain of stability forms a “tongue” with a finite range at $`\omega =0`$. For larger $`ϵ`$ and $`\omega =0`$ the SOW’s do not exist since the driving term dominates the equations and forces a homogeneous pattern. Therefore it seems that the relative suppression of the forcing term by a effective time-averaging due to a fast travelling forcing may account for the existence of soliton states for larger $`ϵ`$.

To prove the feasibility of controlled spatio-temporal forcing and to check the validity and robustness of our predictions in a real system undergoing a Turing instability, we have carried out experiments on the (photosensitive) CDIMA chemical reaction, using the projection of travelling light patterns as the controlled forcing. Experiments were performed in a thermostatized continuously fed unstirred one-feeding-chamber reactor at 4 $`\pm `$ 0.5 C. Structures appear in an agarose gel layer (2% agarose, thickness 0.3 mm, diameter 20 mm). The gel layer was separated from the feeding chamber by an Anapore membrane (Whatman, pore size 0.2 mm) and a nitrocellulose membrane (Schleicher & Schnell, pore size 0.45 mm). Under the chosen set of reagent concentrations (see caption of fig.5), the system at dark spontaneously yields disordered stripe patterns with a wavelength of 0.54 mm. In a typical experiment, parallel light stripes with a characteristic wavelength of precisely 0.54 mm were focused on the gel layer and were moved in the horizontal direction with well-controlled and constant velocity $`v`$. For very small values, trivial travelling stripes, following adiabatically the imposed pattern, were found, as expected. As the passing velocity was increased they readily transformed into the modulated striped mode, as predicted by the theory (see fig. 5a). Furthermore, when considering a slight misfit (see caption of fig.5), a localized structure was observed, propagating in the opposite direction to that of the stripes (see fig. 5b).This again constitutes a neat confirmation of our theoretical findings. Further experiments are presently being conducted to look for other modes of dynamical responses to such a spatio- temporal modulation and will be published elsewhere.

We have described a generic mechanism to induce new pattern dynamics through spatial resonance. The phenomenon is claimed to be generic for systems undergoing a Turing instability. For instance we have already observed the same type of response in the 1d Swift-Hohenberg equation. Within the framework of chemical Turing patterns, the consideration of the 2d case is readily accessible and points out to even richer phenomena.

S.R. is supported by the NSF under contract DMR-0100903. Financial support from DGI (Spain) under projects BXX2000-0638-C02 and BFM2000-0348, and also from European Commission under network HPRN-CT-2002-00312 is acknowledged.


format: text

text:
# Statistical Test of Anarchy

## Abstract

“Anarchy” is the hypothesis that there is no fundamental distinction among the three flavors of neutrinos. It describes the mixing angles as random variables, drawn from well defined probability distributions dictated by the group Haar measure. We perform a Kolmogorov–Smirnov (KS) statistical test to verify whether anarchy is consistent with all neutrino data, including the new result presented by KamLAND. We find a KS probability for Nature’s choice of mixing angles equal to 64%, quite consistent with the anarchical hypothesis. In turn, assuming that anarchy is indeed correct, we compute lower bounds on $`|U_{e3}|^2`$, the remaining unknown “angle” of the leptonic mixing matrix.

preprint: FERMILAB-Pub-03/002-Tpreprint: UCB-PTH-02/62preprint: LBNL-51969

All fermions in the Standard Model of particle physics (SM) seem to come in threes. The three copies of each fundamental matter particle have in common all properties except one – the mass. It is common to say that there are three families, generations, or flavors of each matter particle in the SM. Currently we do not know the reason behind the number three, nor why the matter particles should “repeat” at all. Therefore, it is important to look for any information that may shed light into the origin of flavor.

Within the SM, it has been known for quite some time that different quark flavors can mix quantum mechanically, and that the weak interactions can turn one flavor into another. The “amount” of mixing is summarized by the so-called Cabibbo–Kobayashi–Maskawa (CKM) unitary matrix. The CKM matrix, in turn, can be parameterized by three mixing angles $`\theta _{12},\theta _{13},\theta _{23}`$ and one complex phase $`\delta `$ (throughout, we use the “PDG parameterization” PDG  for the mixing matrices). A non-vanishing phase $`\delta `$ indicates that SM processes can violate CP invariance, distinguishing matter from anti-matter in a subtle manner. With the beautiful data from the $`B`$-factory experiments, we have been able to confirm the CKM framework, and measure all angles and the CP-odd phase with $`O(10)\%`$ accuracy.

A noteworthy feature of the CKM matrix is that it is rather well approximated by the unit matrix, meaning that the quark mixing angles are all small. This fact, combined with the fact that the quark masses are quite distinct (the ratio of the lightest to heaviest quark mass is $`O(10^5)`$), is interpreted as evidence for the existence of some underlying symmetry or physical mechanism that differentiates the quark families and hence explains the hierarchy in the quark masses and the small mixing angles.

In the SM, all neutrinos are exactly massless. This being the case, one can always choose a basis where the Maki–Nakagawa–Sakata (MNS) unitary matrix, the leptonic analog of the CKM matrix, is the unit matrix without loss of generality. This means that there are no SM processes through which one lepton flavor can turn into another. This hypothesis has been indeed confirmed by all experimental searches for charged lepton flavor violation to date PDG .

If neutrinos have masses, and these masses are distinct, there is no reason to expect that the MNS matrix is trivial, and lepton flavor transitions are observable in principle. In this case, the most sensitive probes for lepton flavor transitions are neutrino oscillation processes, through which a neutrino produced in a well-defined flavor state $`\nu _\alpha `$ is detected in a different flavor state $`\nu _\beta `$ after propagating over a macroscopic distance $`L`$. The transition probabilities depend on the mixing angles and the CP-odd phase of the MNS matrix, plus the difference of the neutrino masses-squared, $`\mathrm{\Delta }m_{ij}^2m_i^2m_j^2`$.

Since 1998, there is compelling evidence that neutrino flavor transitions do occur when the neutrinos traverse macroscopic distances. Atmospheric atm , solar solar , and, very recently, reactor neutrino experiments KamLAND  have all observed data consistent with the neutrino oscillation hypothesis. In light of all the experimental evidence, it appears that neutrinos have masses, and that leptonic flavors mix.

There are two striking features regarding the values of the oscillation parameters which are extracted from the current neutrino data. One is that the neutrino masses are extremely small. Neutrino oscillation experiments have determined that the neutrino mass-squared differences are 111We define the neutrino mass eigenvalues such that $`m_2^2>m_1^2`$, and $`\mathrm{\Delta }m_{12}^2<|\mathrm{\Delta }m_{13,23}^2|`$. $`|\mathrm{\Delta }m_{23}^2|=(27)\times 10^3`$ eV2 atm  and $`\mathrm{\Delta }m_{12}^2=(420)\times 10^5`$ eV2 KamLAND . These results, combined with direct searches for neutrino masses PDG , yield that the heaviest neutrino mass is less than $`O(1)`$ eV, over six orders of magnitude smaller than the smallest charged fermion mass of which we know (the electron mass). The other is that, of the mixing angles, two ($`\theta _{12},\theta _{23}`$) are known with some precision, and are both large: $`\mathrm{sin}^22\theta _{23}0.9`$ atm  and $`\mathrm{sin}^22\theta _{12}0.4`$ KamLAND .

Assuming a three family mixing scenario, there are two more parameters in the MNS mixing matrix that are still unknown: $`\theta _{13}`$ and $`\delta `$. In particular, if $`\delta 0`$ neutrino oscillation processes need not conserve CP. Leptogenesis models leptogenesis , on the other hand, try to relate the existence of matter but no anti-matter in the Universe to the CP violation present in the neutrino sector, making its observation of the utmost interest. CP-violating effects parameterized by the CP-odd phase $`\delta `$ of the MNS matrix can be probed in accelerator-based long-baseline neutrino oscillation experiments if, for example, one compares the flavor transformation probabilities of neutrinos and anti-neutrinos (written here in vacuum),

$`P(\nu _e\nu _\mu )P(\overline{\nu }_e\overline{\nu }_\mu )=16s_{12}c_{12}s_{13}c_{13}^2s_{23}c_{23}`$
$`\mathrm{sin}\delta \mathrm{sin}{\displaystyle \frac{\mathrm{\Delta }m_{12}^2L}{4E}}\mathrm{sin}{\displaystyle \frac{\mathrm{\Delta }m_{13}^2L}{4E}}\mathrm{sin}{\displaystyle \frac{\mathrm{\Delta }m_{23}^2L}{4E}},`$ (1)
where $`s_{ij}=\mathrm{sin}\theta _{ij}`$, $`c_{ij}=\mathrm{cos}\theta _{ij}`$. It is well known that the observation of CP violation in neutrino oscillations is possible only if $`\theta _{12}`$ and $`\mathrm{\Delta }m_{12}^2`$ are “large enough” (and the atmospheric parameters are also large, as has been established by the atmospheric data). The KamLAND result has shown that this is the case. The remaining question, therefore, is whether $`\theta _{13}`$ is also large enough to render the experimental search for CP violation possible. The only information we currently have is that $`\theta _{13}`$ is relatively small: $`\mathrm{sin}^2\theta _{13}0.05`$, constrained by the CHOOZ experiment CHOOZ .

The purpose of this letter is two-fold. First, we examine if the current data “requires” new symmetry principles in order to control the structure of the MNS matrix, analogous to the situation in the quark sector. Saying that there is no symmetry principle behind the MNS matrix means there is no fundamental distinction among the three flavors of neutrinos. If this is the case, the MNS matrix is distributed (statistically) according to the bi-invariant Haar measure of group theory, which dictates the probability distribution of the mixing angles. The hypothesis here is that Nature has chosen one point according to this probability distribution. This is the concept of “anarchy” in neutrinos anarchy1 ; anarchy2 . We would like to examine if the data are consistent with anarchy by performing a Kolmogorov–Smirnov (KS) statistical test. We find that they are perfectly consistent.

Second, given the empirical success of anarchy, we study what it has to say about $`\theta _{13}`$. Anarchy prefers large values for $`\theta _{13}`$, meaning that a small $`\theta _{13}`$ would be inconsistent with the anarchical hypothesis. By turning this argument around, we can place a lower limit on $`\theta _{13}`$ at various confidence levels, again using the KS test.

Consider the following situation: there is a model that “predicts” that a certain quantity is described by a probability distribution. For example, one may construct a model that predicts that a given quantity $`x`$ may have any value from 0 to 1, with equal probability. This means that the probability density $`f(x)`$ is 222The probability density function is defined in such a way the probability that $`x`$ has a value between $`x_0`$ and $`x_0+\mathrm{d}x`$ is given by $`f(x_0)\mathrm{d}x`$.

$$f(x)=\{\begin{array}{cc}1\hfill & \mathrm{if}x[0,1],\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}$$
(2)
Let us assume that the value of $`x`$ is known: $`x=x_0`$. The question to be addressed is how well does the result $`x=x_0`$ agree with the model presented above (that the probability density for $`x`$ is given by Eq. (2))? This question can be answered using the KS test. Given that we have drawn the specific value $`x=x_0`$, we would like to test the hypothesis $`_f`$ that the probability distribution associated with the random variable $`x`$ is $`f(x)`$.

In order to do so we define the distribution function 333The distribution function is defined in such a way the the probability that $`x[a,b]`$ is $`F(b)F(a)`$. $`F(x)_{\mathrm{}}^xf(x^{})dx^{}`$. For Eq. (2),

$$F(x)=\{\begin{array}{cc}0\hfill & \mathrm{if}x0,\hfill \\ x\hfill & \mathrm{if}x[0,1],\hfill \\ 1\hfill & \mathrm{if}x1.\hfill \end{array}$$
(3)
We then compare $`F(x)`$ with the best possible guess for a distribution function $`F_{\mathrm{guess}}(x)`$ that can be obtained given that $`x=x_0`$ has been “drawn,” namely,

$$F_{\mathrm{guess}}(x)=\theta (xx_0).$$
(4)
Note that it is very easy to generalize this to $`N`$ random drawings of $`x`$, which yield, say, $`x_0,x_1,\mathrm{},x_{N1}`$ KS\_ref .

The (two-sided) KS statistic (“$`D`$-function”) is defined by KS\_ref

$$D=\mathrm{sup}_x[|F_{\mathrm{guess}}(x)F(x)|].$$
(5)
In the example we have been discussing, $`D_0=x_0`$ if $`x_00.5`$ or $`D_0=1x_0`$ if $`x_00.5`$ (note that the two expressions agree at $`x_0=0.5`$, and we assume that $`x_0[0,1]`$). If the hypothesis $`_f`$ is correct, the probability that a larger value of $`D`$ (i.e. a “worse fit”) would be computed from a different random drawing of $`x`$ is KS\_ref

$$P(DD_0)=2(1D_0),$$
(6)
which is, in the example we have been discussing,

$$P(x_0)=\{\begin{array}{cc}2x_0\hfill & \mathrm{if}x_00.5,\hfill \\ 2(1x_0)\hfill & \mathrm{if}x_00.5.\hfill \end{array}$$
(7)
The smaller the value of $`P(x_0)`$, the less likely it is that $`_f`$ is correct. In this context, we allow statements such as $`_f`$ is only allowed at the $`[1P(x_0)]`$ confidence level.

We wish to apply the test described above to the MNS and CKM mixing matrices for leptons and quarks, respectively. Our model is that the mixing matrices are random variables drawn from a “flat” distribution of unitary $`3\times 3`$ matrices. Following the PDG convention, we define the three mixing angles as in Table 1.

Within this convention, the hypothesis is that the marginalized probability density function is given by (see anarchy2  for a detailed discussion of this point)

$$f(U(3))\mathrm{d}(\mathrm{phases})=f(\mathrm{cos}^4\theta _{13},\mathrm{sin}^2\theta _{12},\mathrm{sin}^2\theta _{12})=1,$$
(8)
where we have integrated over all (both physical and unphysical) complex phases. The mixing angles are defined such that $`\theta _{ij}[0,\pi /2],i,j`$. The probability distribution is flat in $`\mathrm{sin}^2\theta _{12}`$, $`\mathrm{sin}^2\theta _{23}`$, and $`\mathrm{cos}^4\theta _{13}`$. It is clear that $`f=1`$ is correctly normalized,

$$f\mathrm{d}(\mathrm{cos}^4\theta _{13})\mathrm{d}(\mathrm{sin}^2\theta _{12})\mathrm{d}(\mathrm{sin}^2\theta _{23})=1,$$
(9)
as it should be.

Since anarchy implies that the three mixing angles are distributed as uncorrelated random variables according to Eq. (9), we are allowed to perform a separate KS test for each of the three mixing angles. The three distinct $`D`$-functions are (from Eq. (5) and the line that succeeds it),

$`D_{\theta _{13}^0}`$ $`=`$ $`(1\mathrm{sin}^2\theta _{13}^0)^2,`$ (10)
$`D_{\theta _{12}^0}`$ $`=`$ $`1\mathrm{sin}^2\theta _{12}^0,`$ (11)
$`D_{\theta _{23}^0}`$ $`=`$ $`1\mathrm{sin}^2\theta _{23}^0.`$ (12)
The superscript $`0`$ refers to the randomly picked value (i.e., the physical value, “drawn” by Nature) of the corresponding mixing angle. We have assumed that $`\mathrm{sin}^2\theta _{12,23}^0<1/2`$, $`\mathrm{cos}^4\theta _{13}^0>1/2`$. The generalization for all values of $`\theta _{ij}^0`$ is trivial and does not add to our discussion 444Of course, all angles in the quark sector satisfy $`\mathrm{sin}^2\theta _{ij}1/2`$. The same is true of $`|U_{e3}|^2`$, while the preferred values of $`\mathrm{sin}^2\theta _{\mathrm{sol}}`$ are smaller than 1/2 at the three sigma level. We don’t know whether $`\mathrm{sin}^2\theta _{\mathrm{atm}}`$ is less than or greater than 1/2. The experimental information we do have is such, however, that $`\theta _{\mathrm{atm}}`$ and $`\pi /2\theta _{\mathrm{atm}}`$ cannot be discriminated. These “degenerate” solutions lead to the same $`D_{\theta _{23}^0}`$..

Because the three “random variables” are not correlated, we calculate the probability that a different random draw would yield a worse result by computing the area where the product of three one-variable probabilities

$`ϵ^0`$ $`=`$ $`P(\theta _{12}^0)\times P(\theta _{13}^0)\times P(\theta _{23}^0)`$ (13)
$`=`$ $`8\mathrm{sin}^2\theta _{12}^0(2\mathrm{sin}^2\theta _{13}^0\mathrm{sin}^4\theta _{13}^0)\mathrm{sin}^2\theta _{23}^0`$
is worse than the data. Here, $`P(\theta _{ij}^0)=2(1D_{\theta _{ij}^0})`$, as in Eq. (6). Therefore

$$P(KS)=\mathrm{d}(\mathrm{cos}^4\theta _{13})\mathrm{d}(\mathrm{sin}^2\theta _{12})\mathrm{d}(\mathrm{sin}^2\theta _{23})\theta (ϵϵ^0),$$
(14)
where $`ϵ`$ is given by the same expression as $`ϵ^0`$, evaluated at the observed values of the mixing parameters, and $`\theta (ϵϵ^0)`$ is the usual step function. We obtain

$$P(KS)=ϵ\left(1\mathrm{log}ϵ+\frac{1}{2}\mathrm{log}^2ϵ\right).$$
(15)
By using the best fit values $`\mathrm{sin}^2\theta _{12}=0.3`$ Fogli:2002pb  and $`\mathrm{sin}^2\theta _{23}=0.5`$ atm  for the MNS matrix, we find

$$ϵ=2.4(\mathrm{sin}^2\theta _{13}\frac{1}{2}\mathrm{sin}^4\theta _{13}).$$
(16)
Given the bound $`\mathrm{sin}^2\theta _{13}0.05`$, the anarchical hypothesis is consistent with the current data, with probability $`64\%`$.

One can also check whether anarchy works in the quark sector. Using the values tabulated in Table 1, one obtains a probability smaller than $`6\times 10^6`$, implying that the hypothesis that the CKM matrix is a random unitary $`3\times 3`$ matrix is safely discarded. Hence, a fundamental distinction among the three flavors of quarks seems to be required.

Once we have established as consistent the hypothesis that the MNS matrix is a matrix drawn from a random sample of unitary $`3\times 3`$ matrices, we now turn the argument around, and try to place a lower limit on $`\theta _{13}`$. What we require is that $`P(KS)>1P_0`$, where $`P_0`$ is defined to be the confidence level of the limit.

Fig. 1 depicts $`P(KS)`$ for the MNS matrix as a function of $`\mathrm{sin}^2\theta _{13}|U_{e3}|^2`$ within the three sigma bounds allowed experimentally for $`\mathrm{sin}^2\theta _{12}\mathrm{sin}^2\theta _{\mathrm{sol}}`$ and $`\mathrm{sin}^2\theta _{23}\mathrm{sin}^2\theta _{\mathrm{atm}}`$, as tabulated in Table 1. For the best fit values of $`\mathrm{sin}^2\theta _{\mathrm{atm}}`$ and $`\mathrm{sin}^2\theta _{\mathrm{sol}}`$, one is able to “rule out” $`|U_{e3}|^2<0.0007`$ at the two sigma level and $`|U_{e3}|^2<0.00002`$ (which is, curiously, the upper bound for $`|V_{ub}|^2`$) at the three sigma level. Fig. 1 also depicts $`P(KS)`$ as a function of $`\mathrm{sin}^2\theta _{13}|V_{ub}|^2`$ for the CKM matrix within the 90% experimentally allowed ranges defined in Table 1.

Note, however, that these bounds are obtained a posteriori, and turn out to be rather weak. This is due to the fact that the observed values of the angles $`\theta _{12}`$ and $`\theta _{23}`$ agree “too well” with the anarchical hypothesis, hence allowing a larger-than-usual fluctuation for $`\theta _{13}`$. We believe that $`P(KS)`$ is a good tool for testing the anarchical hypothesis against the data, but not as useful a tool for studying what values of $`|U_{e3}|^2`$ are preferred by the hypothesis. For this reason, we choose to make use of another method of obtaining lower bounds on $`|U_{e3}|^2`$ assuming anarchy. This method can be thought of as yielding an a priori prediction for $`|U_{e3}|^2`$, which, we believe, is more appropriate, and will be described promptly.

Let us compute the marginalized probability density function of $`\theta _{13}`$ only. From Eq. (9),

$$f\mathrm{d}(\mathrm{sin}^2\theta _{12})\mathrm{d}(\mathrm{sin}^2\theta _{23})g(\mathrm{cos}^4\theta _{13})=1.$$
(17)
This probability distribution ($`g`$) is rather similar to $`f`$, (see Eq. (9)) but is to be interpreted in a slightly different way. Remember that $`f`$ is the probability distribution function derived for $`U(3)`$ marginalized over all CP-odd phases, including: three unphysical phases that can be removed by redefining fermionic SM fields, two “Majorana phases” which may or may not be physical, depending on whether or not the neutrinos are Majorana fermions, and one “Dirac phase.” Marginalyzing over unphysical phases is the only reasonable procedure to follow, and we have chosen to also marginalize over the physical phase(s) in order to address whether the current information we have on mixing angles fits the anarchical hypothesis. This choice only makes sense if the anarchical hypothesis for the “whole” MNS matrix is consistent, which is the case as we have observed above. Similarly, $`g`$ is the probability distribution function derived for $`U(3)`$ marginalized over all CP-odd phases and the “solar” and “atmospheric” mixing angles.

The single variable probability is

$$P(KS)_{1d}=P(\theta _{13}^0)=4(\mathrm{sin}^2\theta _{13}\frac{1}{2}\mathrm{sin}^4\theta _{13}).$$
(18)
Fig. 2 depicts $`P(KS)_{1d}`$ as a function of $`\mathrm{sin}^2\theta _{13}|U_{e3}|^2`$. Again, we can interpret $`P(KS)_{1d}`$ as the probability that a “worse fit” is obtained assuming that $`\theta _{13}`$ is a random variable drawn from the probability distribution Eq. (17). These bounds are a priori predictions of the anarchy unlike the previous ones, and should be taken more seriously. We exclude $`|U_{e3}|^2<0.011`$ (0.0007) at the two (three) sigma confidence level. Note that the CKM-equivalent $`P(KS)_{1d}=P(V_{ub})`$ is less than $`10^4`$ (this can be easily read off from Fig. 2), again indicating that the anarchical hypothesis in the quark sector can be safely discarded.

Finally, It is worth recalling that anarchy predicts a flat probability distribution for the CP-violating phase $`\delta `$ anarchy2 , and hence the distribution in $`\mathrm{sin}\delta `$ is $`1/|\mathrm{cos}\delta |`$, peaked at $`\mathrm{sin}\delta =\pm 1`$. If anarchy is correct, chances are that the observation of CP violation in long-baseline oscillation experiments is indeed within reach!

We now summarize our results, with more discussions to follow. We have statistically tested the hypothesis that the MNS matrix is a matrix drawn from a random “flat” sample of unitary $`3\times 3`$ matrices. According to the KS test performed, this “anarchical hypothesis” is consistent with the data. The anarchical hypothesis fails the KS test when it is performed with the CKM matrix. Our result is different from other attempts to statistically “test” anarchy. For example, the authors of no\_anarchy  have claimed that the neutrino sector prefers the existence of some symmetry behind neutrino masses and mixing angles to completely random entries. We have not attempted to perform such a “compartive test,” which is, at least, hard to interpret in a well defined way. We do not believe that such tests are capable of indicating whether one hypothesis is favored with respect to the other. Our test has a well defined statistical interpretation, and directly probes whether anarchy in the neutrino sector is a good hypothesis.

Having checked that anarchy is consistent with our current understanding of the MNS matrix, we were able to use the anarchical hypothesis to “predict” the value of the still unobserved mixing angle $`\theta _{13}`$. At the two sigma level, anarchy requires that $`\mathrm{sin}^2\theta _{13}>0.011`$, for example (bound obtained from $`P(KS)_{1d}`$, see Fig. 2). If there is indeed no structure in the leptonic mixing matrix, it seems very likely that one should be able to observe CP-violation in long-baseline neutrino oscillation experiments, as not only are all angles large, but the CP-odd parameter $`\mathrm{sin}\delta `$ is also “predicted” to be large.

We have nothing to say about the value of the neutrino masses. The hypothesis we tested is that the MNS matrix is “random,” independent of whether the masses are degenerate, partially degenerate or hierarchical anarchy2 . Even in the case of non-LMA solutions to the solar neutrino puzzle (currently ruled out at 99.95% C.L. KamLAND ), one can obtain random mixing matrices min\_model . Incidently, it is interesting to note that the the neutrino masses seem to be “less hierarchical” than the charged fermion masses. Assuming that the neutrino masses are not degenerate, it turns out that $`m_3/m_2\sqrt{\mathrm{\Delta }m_{23}^2/\mathrm{\Delta }m_{12}^2}=313`$, not too far away from unity (of course, we do not know $`m_2/m_1\mathrm{}`$). This is consistent with random mass matrices generated via the seesaw mechanism anarchy1 .

We would like to underline important assumptions and limitations of our result. By hypothesis, the probability distributions for the mixing angles are uncorrelated. Our discriminatory procedure does not include information regarding whether the different variables are more likely to be correlated than not. Given the minimal statistics (provided by the fact that we live in only one Universe), adding this sort of information would not lead to different conclusions, although one should start to worry if, say, it turns out that $`\mathrm{sin}^22\theta _{23}=\mathrm{sin}^22\theta _{12}=1`$. One should also be warned that the KS test performed here need not be the most powerful test for the anarchical hypothesis, statistically speaking KS\_ref .

Finally, we emphasize what our result does not imply. Although the anarchical hypothesis is consistent with the data, neutrino mass models which rely on flavor symmetries and nontrivial “textures” are not disfavored in any well defined way. Some are perfectly justified by top-down arguments, including, say, grand unification of matter fields. We would like to point out, however, that the “burden of proof” is with the models that assume that there is structure in the leptonic mixing matrix. The anarchical hypothesis may be viewed as the simplest of flavor models – a model of flavor without flavor. In light of our long experience with quark masses and mixing angles, it is remarkable that, in the neutrino sector, one can do without new symmetry principles in order to appreciate the entries of the MNS mixing matrix.

Note Added – After the first version of the this manuscript became publicly available, a preprint discussing our results Espinosa  appeared. All of the comments contained there apply to this version of our manuscript as well. While we appreciate most of the arguments contained in Espinosa , we disagree with its author in a few key points. Most importantly, we do not agree with the claim that $`\mathrm{sin}^2\theta _{12},\mathrm{sin}^2\theta _{23},\mathrm{cos}^4\theta _{13}[0,1]`$ are “angular variables.” Note that the Haar measure is flat in these variables rather than in the mixing angles, which are convention dependent. Therefore we stand by our claim that a KS test can be used to test the anarchical hypothesis. Furthermore, Espinosa  contains an alternative statistical test of the anarchical hypothesis (a Kuiper’s test, see Espinosa  for details), and the result that maximal mixing is “preferred” is obtained, in qualitative agreement with the results presented here. The author of Espinosa , however, dismisses the result of the Kuiper’s test, claiming that the statistical sample is too small. While we appreciate that some are uneasy about the statistics of very small data samples, namely one chosen by Mother Nature, we point out that the dismissal of such results is not mathematically justified.

###### Acknowledgements.

We thank Paolo Creminelli for very useful questions and discussions, and José Ramon Espinosa for corresponce regarding Espinosa . The work of AdG is supported by the US Department of Energy Contract DE-AC02-76CHO3000. The work of HM is supported in part by the Director, Office of Science, Office of High Energy and Nuclear Physics, Division of High Energy Physics of the U.S. Department of Energy under Contract DE-AC03-76SF00098 and in part by the National Science Foundation under grant PHY-00-98840.


text:
# Universal interface width distributions at the depinning threshold

## Abstract

We compute the probability distribution of the interface width at the depinning threshold, using recent powerful algorithms. It confirms the universality classes found previously. In all cases, the distribution is surprisingly well approximated by a generalized Gaussian theory of independant modes which decay with a characteristic propagator $`G(q)=1/q^{d+2\zeta }`$; $`\zeta `$, the roughness exponent, is computed independently. A functional renormalization analysis explains this result and allows to compute the small deviations, i.e. a universal kurtosis ratio, in agreement with numerics. We stress the importance of the Gaussian theory to interpret numerical data and experiments.

The scaling properties of driven elastic interfaces in random media play an important role in a wide variety of physical situations, ranging from stochastic surface growth to domain walls in disordered magnetic materials, the spreading of fluids on rough substrates, and the dynamics of cracks Stanley.Fractal\_concepts.95 ; Kardar.dynamics.98 . These problems share many features with critical phenomena and provide a challenge for theoretical approaches to disordered systems and non-equilibrium phenomena Nattermann.pinning.97 ; Fisher.Collective.98 ; Narayan.Pinning.93 ; LeDoussalWiese2002a ; Wiese.2loop.01 ; Vannimenus.interface.01 .

Here we study interfaces described by a scalar height function $`h(𝐱)`$, where $`𝐱`$ is the $`d`$-dimensional internal coordinate. We measure the deviation from the mean position as $`u(𝐱)=h(𝐱)h`$, where $`\mathrm{}`$ stands for the spatial average over all $`𝐱`$ of a given interface (cf. Fig. 1). The mean square width of a *single* interface, $`w^2(\{u(𝐱)\})=u^2`$, can be used to characterize its roughness, and explore universal properties: After averaging over the ensemble of interfaces, $`w^2`$ grows with the lateral extension $`L`$ of the system as

$$\overline{w^2}L^{2\zeta }\text{for}L\mathrm{},$$
(1)
where $`\zeta `$ is the roughness exponent.

An interesting property is that, for positive $`\zeta `$, $`w^2`$ fluctuates even in the thermodynamic limit Racz.Random.94 ; Racz.Growth.94 ; Plischke.Curvature.94 . This means that the long-range geometric features of the interface are not characterized by the roughness exponent alone, but require the complete probability distribution $`P(w^2)`$. $`P(w^2)`$ has been computed for several linear stochastic growth equations without disorder as the Edwards-Wilkinson model, the Mullins-Herring model, and the $`1`$-d KPZ model Plischke.Curvature.94 ; Racz.Growth.94 . In these models, the probability distribution $`P(w^2)`$ can be rescaled into a form independent of system size and of microscopic details

$$P(w^2)=(1/\overline{w^2})\mathrm{\Phi }(w^2/\overline{w^2})\text{for}L\mathrm{}.$$
(2)
Although $`\overline{w^2}`$ may contain a non-universal scale, the function $`\mathrm{\Phi }(z)`$ is universal. It has been argued that the shape of $`P(w^2)`$ can thus be used as a sensitive tool, distinct from $`\zeta `$, to distinguish between different universality classes Racz.Random.94 ; Racz.Growth.94 ; Plischke.Curvature.94 ; Antal.1overf.01 . Furthermore, $`\mathrm{\Phi }(z)`$ is expected to converge to a $`\delta `$-function above the upper critical dimension $`d_{\text{uc}}`$. This has motivated attempts to determine $`d_{\text{uc}}`$ for e.g. the KPZ equation Marinari.KPZ.02 . Probability distributions of order parameters have received much attention for related models such as polymers, spin glasses, and random diffusion bouchaud.georges . The quantity we study here, $`P(w^2)`$, is the distribution of the lowest order observable which tests the whole function $`h(x)`$ for $`0<x<L`$. It appears as a fundamental quantity in disordered systems.

The aim of this Letter is to compute the width distribution (WD) $`\mathrm{\Phi }(z)`$ for elastic interfaces driven in random media, exactly at the depinning threshold, numerically and from field theory. As in the linear problems treated earlier, we confirm the existence of universal properties in various dimensions $`d`$ and with several functional forms of the elasticity. The surprising finding is that in all cases $`P(w^2)`$ (i.e. its shape) is extremely well approximated by a simple generalized Gaussian approximation (GA), without any fit parameter, and depends only on $`\zeta `$, which is determined independently. This suggests that the complicated morphology of interfaces (cracks, domain walls, etc.) may be rendered by a simple ansatz of independent modes with a characteristic decay. This may have important consequences for the analysis of numerical and experimental data. Our numerical results are then understood within a functional renormalization group calculation, detailed in a companion paper technical .

We consider the zero temperature equation of motion of an interface given by

$$_th(𝐱,t)=\frac{E}{h(𝐱,t)}=f+\eta (𝐱,h(𝐱,t))\frac{E_{\text{el}}}{h(𝐱,t)},$$
(3)
where the functional $`E(\{h,𝐱\})`$ represents the total energy incorporating potential energy due to the driving force $`f`$, the short-range correlated disorder force $`\eta (𝐱,h)`$, as well as its internal convex elastic energy $`E_{\text{el}}`$. Equation (3) is non-linear, and has not been solved exactly. We are interested in the depinning limit ($`f=f_c`$) where the velocity of the elastic manifold goes to zero. We use periodic boundary conditions, and recall that the WD $`\mathrm{\Phi }(z)`$, although independent of small scale details, does depend on the boundary condition at large scale.

For our numerical study we use a very efficient algorithm Krauth.zeta.01 ; Rosso.MonteCarlo.02 ; Rosso.contact.02  which directly determines critical forces $`f_c`$ as well as the critical interface $`h_c(𝐱)`$ for a wide range of models. In particular we calculate the WD for interfaces of dimensions $`d=1`$ and $`d=2`$, where the elastic energy has the harmonic form $`E_{\mathrm{el}}(\{h,𝐱\})(h)^2`$. We have also tested the universality of $`\mathrm{\Phi }(z)`$ in $`d=1`$, by means of a directed polymer model with an anharmonic quartic elasticity, and for a lattice model with hard local constraint, which have the same $`\zeta =0.63`$ Krauth.zeta.01 . As expected, $`\mathrm{\Phi }(z)`$ is always size independent and the WD associated to non-harmonic models can be distinguished from the one resulting from an harmonic elasticity. The harmonic models, in fact, have an exponent $`\zeta =1.2`$, and thus belong to a different universality class.

For our field theory calculation we use the functional renormalization group method (FRG) originally developped to one loop to describe the model with harmonic elasticity and correct the predictions of dimensional reduction Nattermann.pinning.97 ; Narayan.Pinning.93 . Recently a renormalized field theory was constructed to 2-loop order Wiese.2loop.01  which overcomes the deficiencies of the 1-loop analysis; notably to distinguish between statics and driven dynamics, and to account for the large values of the roughness exponent $`\zeta `$ measured e.g. in Rosso.MonteCarlo.02 ; Rosso.contact.02 ; Rosso.manifolds.02  as compared to an earlier conjecture Narayan.Pinning.93  $`\zeta =(4d)/3`$. We find that the FRG both suggests the GA as a lowest order approximation in $`ϵ=4d`$ and allows to define and compute universal ratios which probe high cumulants of $`P(w^2)`$ and deviations from the GA, and are thus more sensitive to details of the universality class. The simplest of them is the generalized kurtosis

$$R=\frac{_{x,y}\overline{\left(u(x)^2u(y)^2\right)}^\mathrm{c}}{2_{x,y}\left(\overline{u(x)u(y)}\right)^2},$$
(4)
where the subscript $`\mathrm{c}`$ indicates the connected expectation value. $`R`$ is found to be small but non-zero. This directly proves that the correct description of interfaces must go beyond the independent-mode picture.

To introduce the Gaussian approximation in the most elementary way, we first recall Racz.Random.94  the simple periodic random walk of size $`L`$, with a Fourier expansion

$$u(x)=\underset{n=1}{\overset{\mathrm{}}{}}a_n\mathrm{cos}\left(\frac{2\pi n}{L}x\right)+b_n\mathrm{sin}\left(\frac{2\pi n}{L}x\right).$$
(5)
The standard Gaussian probability measure associated with $`u`$, i.e. $`𝒫[u]\mathrm{exp}[\frac{1}{2}_0^Ldx(\frac{u(x)}{x})^2]`$ gives

$$𝒫[u]\mathrm{exp}\left[\underset{n=1}{\overset{\mathrm{}}{}}\frac{(\pi n)^2}{L}(a_n^2+b_n^2)\right].$$
(6)
The probability distribution $`P(w^2)`$

$$P(w^2)=𝒟[u]\delta (w^2u^2)𝒫[u]$$
(7)
is obtained from the generating function of its moments

$$W(\lambda )=_0^{\mathrm{}}dw^2P(w^2)\mathrm{e}^{\lambda w^2}.$$
(8)
Writing Eq. (8) as an integral over $`a_n`$ and $`b_n`$, we obtain

$`Z(\lambda )={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle da_ndb_n\mathrm{e}^{\frac{(\pi n)^2}{L}(a_n^2+b_n^2)}\mathrm{e}^{\frac{\lambda }{2}(a_n^2+b_n^2)}}`$
$`W(\lambda )={\displaystyle \frac{Z(\lambda )}{Z(0)}}={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left(1+{\displaystyle \frac{\lambda }{2}}{\displaystyle \frac{L}{(\pi n)^2}}\right)^1.`$ (9)
For the random walk Eq. (6), $`P(w^2)`$ can be obtained exactly by inverse Laplace transform of Eq. (Universal interface width distributions at the depinning threshold):

$$P(w^2)=\frac{4\pi ^2}{L}\underset{n>0}{}n^2(1)^{n+1}\mathrm{e}^{2w^2(\pi n)^2/L}.$$
(10)
Using $`\overline{w^2}=\frac{\mathrm{d}W}{\mathrm{d}\lambda }|_{\lambda =0}=\frac{L}{12}`$, Eq. (10) can be written in a scaling form

$`\mathrm{\Phi }(z)`$ $`=`$ $`\overline{w^2}P(w^2),z=w^2/\overline{w^2}`$ (11)
$`=`$ $`{\displaystyle \frac{\pi ^2}{3}}{\displaystyle \underset{n>0}{}}n^2(1)^{n+1}\mathrm{e}^{\frac{\pi ^2}{6}zn^2}.`$
The size dependence thus appears only through the average width $`\overline{w^2}`$. We can generalize Eq. (6), where each mode $`a_n,b_n`$ has a weight $`n^2`$, to an arbitrary function of independent Fourier modes

$$𝒫_{\text{gauss}}[u]\mathrm{exp}\left[\frac{L}{4}\underset{n>0}{}(a_n^2+b_n^2)G^1\left(\frac{2\pi n}{L}\right)\right].$$
(12)
which, in real space, corresponds to

$$𝒫_{\text{gauss}}[u]\mathrm{exp}\left[\frac{1}{2}_0^L_0^Ldxdyu(x)G_{xy}^1u(y)\right].$$
(13)
The function $`G_{xy}=\overline{u(x)u(y)}`$ is the exact disorder-averaged 2-point function and can be computed from numerical data. This allows to obtain $`𝒫_{\text{gauss}}`$ even for a finite system. In the thermodynamic limit, $`𝒫_{\text{gauss}}[u]`$ is obtained from the behavior of $`G_{xy}=G_{xy}`$ for large $`|xy|`$ (small $`q`$), where $`G(q)C/q^{d+2\zeta }`$. This means that a single observable, $`\zeta `$, fixes $`𝒫_{\text{gauss}}[u]`$ on large scales.

We again determine the generating function for the moments, but for arbitrary $`\zeta `$ and $`d`$:

$$W(\lambda )=\underset{q0}{}\left(1+2\stackrel{~}{\lambda }G(q)\right)^{1/2},$$
(14)
where $`\stackrel{~}{\lambda }=\lambda /L`$, $`q=2\pi n/L`$, $`n^d`$. Due to the symmetry $`qq`$, no fractional power appears in Eq. (14), as in Eq. (Universal interface width distributions at the depinning threshold), where the exponent $`1`$ stems from the double sum over the $`a_n`$ and $`b_n`$. An explicit sum over poles allows to obtain $`\mathrm{\Phi }_{\text{gauss}}(z)`$ for all $`\zeta `$ and $`d`$ with excellent precision. All GA interfaces $`\{u(x)\}`$ can be directly sampled by Monte Carlo methods. For details, including the extension to open boundary conditions, see Antal.1overf.01 ; RossoKrauth03

In Fig. 2 we compare, for different models, the exact scaling function $`\mathrm{\Phi }(z)`$ to $`\mathrm{\Phi }_{\text{gauss}}(z)`$, using $`G(n)=C/n^{d+2\zeta }`$. The roughness exponent was previously obtained using both field theory Wiese.2loop.01  and numerical methods Rosso.MonteCarlo.02 . The agreement between $`\mathrm{\Phi }`$ and $`\mathrm{\Phi }_{\text{gauss}}`$ is clearly spectacular. The scatter of the data, visible in Fig. 2, is mostly due to the finite width of histogram bins.

Tiny—yet significant—differences between $`\mathrm{\Phi }`$ and $`\mathrm{\Phi }_{\text{gauss}}`$ are best resolved in the integrated probability distributions, which need no discretization. The difference between the integrated distributions is

$$\mathrm{\Delta }H(z)=_0^zdt(\mathrm{\Phi }_{\text{gauss}}(t)\mathrm{\Phi }(t)),$$
(15)
where $`H(z)=_0^zdt\mathrm{\Phi }(t)`$ is the fraction of samples with a renormalized width below $`z`$.

In Fig. 3, we show $`\mathrm{\Delta }H(z)`$ obtained from $`N=\mathrm{2\hspace{0.17em}10}^5`$ independent samples. Statistical fluctuations in this quantity are of order $`1/\sqrt{N}`$ and the signal-to-noise ratio would be smaller than $`1`$ if $`N`$ was an order of magnitude smaller. The absence of systematic finite-size effects shows that the asymptotic regime of large interfaces has been reached and thus to conclude that the exact distribution for large systems is not Gaussian.

We now discuss the field theoretical calculation. To lowest order in perturbation theory, we show that the generalized Gaussian approximation appears naturally. This is instructive since it identifies the diagrams which are obtained by assuming the theory to be Gaussian, albeit non-trivial, since it involves a non-trivial roughness exponent $`\zeta `$. Using dynamical field theoretic methods Wiese.2loop.01 , one starts again from the Laplace transform $`W(\lambda )`$ and expands in powers of the correlator of the pinning force $`\mathrm{\Delta }(u)`$. To lowest order one finds technical  that $`\mathrm{log}W(\lambda )`$ is the sum of all connected 1-loop diagrams. The loop with $`N`$ disorder vertices and $`N`$ insertions of $`w^2`$ is

$$\frac{1}{2N}\underset{q}{}\left(\frac{2\stackrel{~}{\lambda }\mathrm{\Delta }(0)}{(q^2)^2}\right)^N,$$
(16)
where the sums over $`q`$ thus run over a $`d`$-dimensional lattice with spacing $`\frac{2\pi }{L}`$, and the 0-mode is excluded, as appropriate for periodic BC. Resumming (16) over $`q`$ would give Eq. (14) with $`G(q)1/q^4`$, i.e. the dimensional reduction (Larkin) result. In fact, the FRG tells us that the calculation should be performed with the running disorder $`\mathrm{\Delta }(0)\mathrm{\Delta }_l(0)=\mathrm{e}^{(ϵ2\zeta )l}\stackrel{~}{\mathrm{\Delta }}^{}(0)`$ where $`\stackrel{~}{\mathrm{\Delta }}^{}(0)`$ is the (non-universal) value of the fixed point Wiese.2loop.01 . For the present case of periodic boundary conditions and momentum infrared cutoff, one can replace $`l\mathrm{log}(1/q)`$, and finally obtains Eq. (14) with $`G(q)=C/q^{d+2\zeta }`$. This calculation is valid to dominant order in $`ϵ=4d`$, i.e. near $`d=4`$. If the same class of diagrams are resummed in any $`d`$ it leads to the GA, as we now illustrate considering e.g. the second connected cumulant of the WD. This cumulant is not connected w.r.t. $`h`$, and thus there is an exact relation:

$$\overline{(w^2)^2}^c=\overline{(w^2)^2}\overline{(w^2)}^2=2\left(1+R\right)_{x,y}G_{xy}^2.$$
(17)
The first term results from Wick’s theorem and would be the full result if the measure were Gaussian. Analogous formulae exist for higher cumulants, and if the measure of $`h`$ is purely Gaussian can be resummed into Eq. (14). Even though the GA is not exact, the deviations, given by the last term in (17) are expected to be small; indeed they are of order $`ϵ^4`$. Thus the GA is already exact to the two lowest leading orders $`ϵ^2`$ and $`ϵ^3`$, which explains why it is so accurate even in low dimension.

The calculation of the deviations using the field theory is delicate technical . The kurtosis $`R`$ in Eq. (4) which characterizes the importance of non-Gaussian effects is found to be $`R=0.13ϵ^2`$ to lowest order for small $`ϵ=4d`$. It is easy to see that this strongly overestimates $`R`$ in low dimensions. Another method is to work in fixed dimension and to truncate to one loop, yielding $`R=0.036`$ ($`d=3`$), $`R=0.048`$ ($`d=2`$), $`R=0.01`$ ($`d=1`$), which in view of the numerical results below seems to underestimate $`R`$. The small values obtained in low dimensions arise from kinematic constraints in the diagrams, presumably a genuine effect indicating large corrections from higher orders in $`ϵ`$ to the above $`O(ϵ^2)`$ even in $`d=3`$. Note that the sign of the result indicates a distribution more peaked than a Gaussian and is in agreement with the only other known exact result cecile  for the (random field) statics in $`d=0`$, $`R=0.080865..`$.

We have computed from Eq. (4) the generalized kurtosis function, in a model-independent way. We have checked on the $`1d`$ harmonic model that, as $`\mathrm{\Phi }(z)`$, $`R`$ is not affected by finite size effects and – using $`10^6`$ samples – we find: $`R=0.054\pm 0.002`$ ($`1d`$ harmonic $`L=256`$); $`R=0.067\pm 0.002`$ ($`1d`$ anharmonic $`L=64`$); $`R=0.053\pm 0.002`$ ($`2d`$ harmonic $`L=32`$). As proven by FRG calculations, $`R`$ is small but definitely different from zero. Direct information on the non-Gaussian effects can also be obtained from the Fourier transforms of the interfaces $`u(x)`$ in Eq. (5). In Fourier space, for $`d=1`$, the expression of $`R`$ is:

$$R=\frac{_{n_1,n_2}\overline{(a_{n_1}^2+b_{n_1}^2)(a_{n_2}^2+b_{n_2}^2)}^c}{2_n\overline{(a_n^2+b_n^2)}^2}.$$
(18)
We remark that $`R`$ detects correlations in the disorder-averaged fourth moments $`\overline{|u(q_1)|^2|u(q_2)|^2}`$, which cannot be expressed simply through the second moments $`\overline{|u(q)|^2}`$.

To summarize, we have computed both numerically and within field theory the width distribution of critical configurations at depinning, with consistent results. The shapes of the distributions are strongly dominated by the value of $`\zeta `$. On the other hand, it will be difficult to distinguish different universality classes from the forms of $`\mathrm{\Phi }(z)`$, if their roughness exponents are similar. Other universal quantities, such as $`R`$ defined here, directly involve the non-Gaussian part of the distribution. Their precise determination however requires more work, both numerically and within field theory. Also, since the WD is so tied up to $`\zeta `$, finite size effects in both quantities are connected. Finite-size effects will need to be well understood in order to resolve open issues Buldyrev.anisotropic.93 ; Rosso.manifolds.02 ; LeDoussalWiese2002a  concerning $`d_{\mathrm{uc}}`$ for the anisotropic depinning class. It would be interesting to carry similar calculations on other pure, and disordered models, in particular for the equivalent static system.

We thank Z. Rácz for stimulating discussions. K.J.W. is supported by Deutsche Forschungsgemeinschft under grant Wi 1932/1-1.


text:
# Molecular electronics and first-principles methods

## I Introduction

Nanoscale electronics constitutes the backbone of present and future technological advances or, in other words, of Nanotechnology. Ultimately, the functionality of electronic devices will rely on the conduction properties of nanoscopic regions composed of a number of atoms that can range typically from several thousands down to a single one. Probably, the most promising research field within nanoscale electronics is what is known as molecular electronics. The main idea behind molecular electronics is the possibility that functional units can be built out of very stable and well-characterized moleculesJoachim et al. (2000) such as fullerenesDresselhaus et al. (1996) or carbon nanotubesSaito et al. (1998); Dekker (1999). The simplest molecular device consists of two large metallic electrodes, several nanometers apart, joined by a molecule or molecules anchored to themReed et al. (1997). We will name this system molecular bridge. Alternatively, the electrodes can be simply connected by a chain of atoms of the same element as the electrodesMuller et al. (1992); Agrait et al. (1993). These bridges, generically known as atomic contacts or metallic nanocontacts, are not expected to be of any practical application, but constitute an excellent benchmark for us to learn about the world of electrical transport at the atomic scale.

The design and fabrication of electronic devices at the molecular and atomic scale poses new challenges for theorists who must develop and implement new techniques to address the upcoming problems. The basics to calculate the zero-bias, zero-temperature conductance $`G`$ of a molecular bridge or metallic nanocontact were established by Landauer long before the concept of molecular electronics was commonplace. In Landauer’s formalism $`G`$ is simply related to the quantum mechanical transmission probability $`T`$ of electrons at the Fermi level to go from one electrode to the otherDatta (1995)(we assume spin degeneracy):

$$G=G_0T(E_\mathrm{F})=\frac{2e^2}{h}T(E_\mathrm{F}).$$
(1)
The transmission can be easily estimated on generic considerations for metallic nanocontactsCuevas et al. (1998); Hasmy et al. (2001), but it is much more difficult to do so for molecular bridges. The reason is simple: We do not know a priori where the electrode Fermi energy $`E_\mathrm{F}`$ lies with respect to the molecular levels. The positioning of $`E_\mathrm{F}`$ depends on two closely-related factors: (i) the coupling or hybridization of the molecular levels with the free-electron levels in the electrodes and (ii) the charge transfer between molecule and electrodes. The conductance is thus strongly dependent on the particular molecule, the detailed atomic arrangement of the electrodes where the molecule binds, and the chemical nature of the various elements at play. In order to give an answer to this problem we have to rely on first-principles or ab-initio calculations. To calculate the atomic arrangement of the electrodes or the way the molecule binds to the electrodes is a major problem in itself that deserves a whole study on its own and will not be given further consideration in these notes. Still, even if we ignore this important issue, to implement Landauer formalism requires to know the electronic structure of a finite region embedded in an infinite system with no particular symmetry or periodicity. This highly non-trivial problem is what we discuss in what follows.

## II The basics

The main advantage of the numerical implementation that we discuss here with respect to similar proposals that have appeared in the literatureLang (1995); Hirose and Tsukada (1995); Brandbyge et al. (2002); Taylor et al. (2001) is the use of a standard quantum chemistry code such as GAUSSIAN98Frisch et al.  (see Refs. Palacios et al., 2001a; Palacios et al., 2002; Palacios et al., 2001b; Damle et al., 2001). The GAUSSIAN98 code provides a versatile method to perfom first-principles or ab-initio calculations of clusters, incorporating the major advancements in the field in terms of functionals, basis sets, pseudopotentials, etc.. The procedure goes as follows. A standard electronic structure calculation of the region that includes the molecule and a significant part of the electrodes is performed \[see Fig. 1(a)\]. This calculation can be performed within any mean-field-like approach: Hartree-Fock or density functional (DF) theory in any of its multiple approximations. As far as transport is concerned, the self-consistent hamiltonian $`\widehat{H}`$ (or Fock matrix $`\widehat{F}`$) of this central cluster or supermolecule contains the relevant information. The retarded(advanced) Green’s functions associated with $`\widehat{F}`$

$$\left[(E\pm i\delta )\widehat{F}\right]\widehat{G}^{(\pm )}=\widehat{I}$$
(2)
needs to incorporate the rest of the infinite electrodes:

$$\left[(E\pm i\delta )\widehat{F}\widehat{\mathrm{\Sigma }}^{(\pm )}\right]\widehat{G}^{(\pm )}=\widehat{I}$$
(3)
In this expression

$$\widehat{\mathrm{\Sigma }}^{(\pm )}=\widehat{\mathrm{\Sigma }}_\mathrm{R}^{(\pm )}+\widehat{\mathrm{\Sigma }}_\mathrm{L}^{(\pm )},$$
(4)
where $`\widehat{\mathrm{\Sigma }}_\mathrm{R}`$($`\widehat{\mathrm{\Sigma }}_\mathrm{L}`$) denotes a self-energy operator that accounts for the part of the right(left) semi-infinite electrode which has not been included in the calculation, and $`\widehat{I}`$ is the unity matrix. In a non-orthogonal basis the Green function takes the form

$$G_{NO}^{(\pm )}=S\left[(E\pm i\delta )SF_{NO}\mathrm{\Sigma }_{NO}^{(\pm )}\right]^1S$$
(5)
where S is the overlap matrix. To simplify the notation we have dropped the energy argument in the Green function and the selfenergy. The added self-energy matrices can only be explicitly calculated in ideal situations. We choose to describe the bulk electrode with a Bethe lattice tight-binding model with the coordination and parameters appropriate for the electrodes included in the cluster. Figure 1(c) depicts schematically a Bethe lattice of coordination three in two dimensions. The advantage of choosing a Bethe lattice resides in that it reproduces fairly well the bulk density of states of most commonly used metallic electrodes. For each atom in the outer planes of the cluster, we choose to add a branch of the Bethe lattice in the direction of any missing bulk atom (including those missing in the same plane). In Fig. 1(b) the directions in which Caylee tree branches are added are indicated by bigger atoms which represent the first atom of the branch in that direction. Assuming that the most important structural details of the electrode are included in the central cluster, the Bethe lattices should have no other relevance than that of introducing a generic reservoir.

The above definition of the Green function is not a standard one. As many other ab-initio implementations, GAUSSIAN98 makes use of non-orthogonal localized orbital basis sets. The appearance of $`S`$ in the definition of the Green function guarantees that the total charge of the system $`N`$ can be obtained through the standard expression used in quantum chemistry:

$$N=\mathrm{Tr}[PS],$$
(6)
where the trace runs over all the atomic orbitals of the cluster and where $`P`$ is the density matrix defined in a standard way in terms of the retarded Green function as

$$P=\frac{1}{\pi }_{\mathrm{}}^{E_\mathrm{F}}\mathrm{Im}\left[S^1G_{\mathrm{NO}}^{(+)}S^1\right]$$
(7)
$`E_\mathrm{F}`$ is the Fermi energy which can be set by imposing overall charge neutrality in the cluster. The integral in Eq. (7) can be efficiently calculated along a contour in the complex planeTaylor et al. (2001); Brandbyge et al. (2002); Xue et al. . As explained in the next section, we force GAUSSIAN98 to evaluate $`F`$ using the density matrix defined in Eq. 7 instead of the standard one obtained from filling up the molecular orbitals obtained from diagonalizing the Fock matrixSzabo and Ostlund (1989). It is also interesting to note that the applicability of this approach, which we name the Gaussian Embedded Cluster Method (GECM) can be easily implemented in any quantum-chemistry or ab-initio code as long as it is based on a localized orbital expansion of the self-consistent, single-electron wavefunctions.

The transmission probability that appears in Eq. 1 can be calculated through the general expressionDatta (1995)

$$T(E)=\mathrm{Tr}[\widehat{\mathrm{\Gamma }}_\mathrm{L}\widehat{G}^{(+)}\widehat{\mathrm{\Gamma }}_\mathrm{R}\widehat{G}^{()}],$$
(8)
where Tr denotes the trace over all the orbitals of the cluster and where the operators $`\widehat{\mathrm{\Gamma }}_{\mathrm{R}(\mathrm{L})}`$ are given by $`i(\widehat{\mathrm{\Sigma }}_{\mathrm{R}(\mathrm{L})}\widehat{\mathrm{\Sigma }}_{\mathrm{R}(\mathrm{L})}^{})`$. Finally, in order to single out the contribution of individual channels to the current, one can diagonalize the transmission matrix. While the size of the matrix $`[\widehat{\mathrm{\Gamma }}_L\widehat{G}^{()}\widehat{\mathrm{\Gamma }}_R\widehat{G}^{(+)}]`$ can be as large as desired, the number of eigenvalues with a significant contribution will be typically much smaller, being determined by the narrowest part of the nanocontact or the molecule. The symmetry of each channel can be identified by looking at its weight on the atomic orbitals of the cluster.

Notice that, although the above expression provide us with the transmission as a function of energy, only the value of $`T`$ at the Fermi energy, $`T(E_\mathrm{F})`$, is strictly correct (within the approximations inherent to DF theory). Away from $`E_\mathrm{F}`$, $`T(E)`$ is only indicative. For instance, if a gate voltage is applied and the system gets charged or discharged, $`T`$ needs to be recalculated for the new Fermi energy. To obtain the current $`I`$ at non-zero temperature and finite bias voltage $`V`$ one integrates the transmission probability:

$$I=\frac{2e}{h}_{\mathrm{}}^{\mathrm{}}T(E,V)[f_\mathrm{L}(EeV/2)f_\mathrm{R}(E+eV/2)]dE,$$
(9)
where $`f_\mathrm{L}`$ and $`f_\mathrm{R}`$ are the Fermi distribution functions of the left and right electrodes, respectively. Notice, once again, that this innocent-looking expression is not a trivial generalization of Eq. 1 extended to finite bias voltages and temperature. The transmission coefficient that appears in the integrand $`T(E,V)`$ depends on temperature, but, most importantly, depends on $`V`$. To calculate this transmission coefficient, knowledge of the Fock matrix in the presence of a bias voltage is required which, in turn, depends on the density matrix out of equilibrium. Based on standard non-equilibrium Green function techniques, it has been repeatedly shown in the literatureDatta (1995) that

$$P=\frac{i}{2\pi }_{\mathrm{}}^{\mathrm{}}dE\left[S^1G^<S^1\right]$$
(10)
where

$$\widehat{G}^<(E)=i\widehat{G}^{()}(E)\left[f_\mathrm{L}(E)\widehat{\mathrm{\Gamma }}_\mathrm{L}(E)+f_\mathrm{R}(E)\widehat{\mathrm{\Gamma }}_R(E)\right]\widehat{G}^{(+)}(E).$$
(11)
There are several technical and conceptual issues regarding the computation of the above equations that need some discussion. We refer the reader to the literature for a detailed account of most of themTaylor et al. (2001); Palacios et al. (2002); Brandbyge et al. (2002). Here we will simply mention two of them. Maybe the most important conceptual issue is that related to the fact that using DF theory for non-equilibrium situations is not justified: DF theory is a ground-state theory. We will make use of it simply because we do not know a better way to deal with out-of-equilibrium problems in an operative way. Second, there is technical issue that, in our view, has not received due attention in the literature and that we find it worth a detailed discussionLouis et al. : The out-of-equilibrium electrostatics. In most self-consistent approachesLang (1995); Taylor et al. (2001); Damle et al. (2001); Brandbyge et al. (2002) one solves the Poisson equation with boundary conditions appropriate for the electrode geometries. This is equivalent to imposing an external electric field and it can only be done for simple geometries such as infinite planes. The Keldysh formalism, however, does not necessarily require to deal with the Poisson equation if a significant part of the metallic electrodes has already been included in the central cluster. Thus, one can consider realistic electrode geometries if necessary. One simply imposes an electrochemical potential difference between electrodes $`\mu _{ec}^\mathrm{L}\mu _{ec}^\mathrm{R}=eV`$ which is actually what a battery does. There are two contributions to $`\mu _{ec}`$:

$$\mu _{ec}=\mu _{\mathrm{chemical}}+\mu _{\mathrm{electrostatic}},$$
(12)
where we define

$`\mu _{\mathrm{chemical}}`$ $`=`$ $`\underset{\delta N0}{lim}\delta E_{\mathrm{core}}/\delta N`$
$`\mu _{\mathrm{electrostatic}}`$ $`=`$ $`\underset{\delta N0}{lim}\delta E_{ee}/\delta N.`$
$`E_{\mathrm{core}}`$ is the standard energy contribution to the total energy from the cores of the atoms and $`E_{ee}`$ is the contribution from the electron-electron interaction111Note that this definition is not a standard one in the sense that one usually thinks of external fields when referring to the electrostatic potential contribution to the electrochemical potential. Instead we will make use of this term when referring to any contribution coming from electron-electron interactions. Strictly speaking, however, one should not use the term electrostatic potential since there are exchange and correlation contributions included in DF theory.. Specifically

$`\delta E_{\mathrm{core}}`$ $`=`$ $`\mathrm{tr}[\delta Ph_{core}]`$
$`\delta E_{\mathrm{ee}}`$ $`=`$ $`\mathrm{tr}[\delta Ph_{ee}],`$
where $`h_{core}`$ and $`h_{ee}`$ are the core and interaction matrices composing the Fock matrix. In the selfconsistent procedure one electrode gets charged and the other discharged (by the same amount in case of mirror symmetry). Once self-consistency has been achieved, an electrostatic potential difference $`V`$ between atoms one or two layers inside opposite electrodes develops while their chemical potentials are similar, i.e., they remain neutral. On the other hand, the electrostatic potential drop will be smaller than $`V`$ between atoms at the surface of opposite electrodes since they carry the charges. We show this process schematically in Fig. 2.

## III Example: Au nanocontacts

### III.1 Clean nanocontacts

Calculations that illustrate the methodology described above have been carried out for the Au nanocontact shown in the inset of Fig. 3: Two (111) planes containing 19 atoms each plus a four-atom chain. Interatomic bulk distances have been taken for the whole cluster (2.88Å), although one should keep in mind that an analysis of the structural stability of the nanocontact model is an important issue when one is interested in the interpretation of experimental data. The results shown next correspond to zero temperature and for the DF calculations we have used the Becke’s three-parameter hybrid functional using the Lee, Yang and Parr correlation functional (B3LYP)Becke (1993) together with the semilocal shape consistent pseudopotential (SCPP) and minimal basis sets of Christiansen et al.Pacios and Christiansen (1985); Hurley et al. (1986). Fig. 3 shows the average on-site energies of the $`5d6s6p`$ orbitals for all the atoms of the nanocontact at zero and 5 V bias. This magnitude reflects the local self-consistent electrostatic potential on each atom. The results represented by white dots correspond to zero bias where one can see the characteristic potential profile due to the parallel surfaces reflected on the chain atoms lying across the vacuum between surfaces. The black dots correspond to a 5 V bias. One can see how an electrostatic potential drop of almost 5 V has developed between plates without having imposed an external field, just an electrochemical potential difference. The major drop in the electrostatic potential occurs at the chain-plane contacts. In the case of zero bias the results are symmetric with respect to the geometric center of symmetry, as expected. Instead, a similar symmetry is absent for 5 V. Namely, whereas the potential drop between the left electrode and the first atom in the chain is 1.92 eV, it is only of 0.85 eV between the chain end and the right electrode. This feature seems common to all previously reported resultsBrandbyge et al. (1999, 2002) and reflects bulk band structure features. The transmission $`T`$ for zero and 5 V is shown in Fig. 4. $`T`$ oscillates around the Fermi energy with an upper limit value of one for both cases. This limit is well documented theoretical and experimentally for Au nanocontactsAgrait et al.  since there is only one $`s`$-like channel around the Fermi energy for Au chains. The oscillations are due to scattering at the electrode-chain contact. This explains why conductance histograms exhibit the lowest peak slightly below $`G=2e^2/h`$Rego et al. . It is worth noting that although the total current calculated by integrating $`T(E,V=0)`$ in a window \[-$`V`$/2, $`V`$/2\] does not differ much from that obtained with the full non-equilibrium approach, the differential conductance is significantly different. For instance, the gap below $`0.5`$ eV is partially filled at finite bias.

### III.2 H2 molecular bridges in Au nanocontacts

We finally examine the possible consequences of adding hydrogen (H2) to a Au nanocontact. A similar experimental analysis has been recently reported for Pt with conclusive results: The presence of H2 molecules alters the conductance histograms of clean nanocontactsSmit et al. (2002). It is thus important to determine the extent of this influence when interpreting conductance histograms under low-vacumm conditions or in the presence of gases. We considered here two electrodes in the form of pyramids in the (001) direction with a H2 molecule anchored between them (see inset in Fig. 5). The orientation of the molecule has been chosen on the basis of the report in Ref. Smit et al., 2002 for Pt, although other alternatives have been recently proposed222Y. García, J. J. Palacios, E. Louis, J. A. Vergés, E. SanFabián, and A. J. Pérez-Jiménez, unpublished.. We have checked that this orientation is stable for a range of distances $`d`$ between outer planes which are fixed in the relaxation. Finally we have computed the conductance for $`d=`$8 Å(shown in Fig. 5). As one can see the conductance at the Fermi energy is greatly reduced from the quantum of conductance $`G_0`$ due to the presence of the H2 molecule. A detailed analysis of this problem will be published elsewhere333Y. Garcia, J. J. Palacios, A. J. P ’erez-Jiménez, E. Louis, E. SanFabián, and J. A. Vergés, unpublished.

Financial support by the spanish MCYT (grants BQU2001-0883, PB96-0085, and MAT2002-04429-C03) and the Universidad de Alicante is gratefully acknowledged.