text:
# First order phase transitions: equivalence between bimodalities and the Yang-Lee theorem

## Abstract

First order phase transitions in finite systems can be defined through the bimodality of the distribution of the order parameter. This definition is equivalent to the one based on the inverted curvature of the thermodynamic potential. Moreover we show that it is in a one to one correspondence with the Yang Lee theorem in the thermodynamic limit. Bimodality is a necessary and sufficient condition for zeroes of the partition sum in the control intensive variable complex plane to be distributed on a line perpendicular to the real axis with a uniform density, scaling like the number of particles.

Phase transitions are universal phenomena which have been theoretically understood at the thermodynamic limit of infinite systems as anomalies in the associated equation of state (EoS). They have been classified according to the degree of non-analyticity of the thermodynamic potential at the transition point. For finite systems, it is generally believed that phase transitions cannot be defined since thermodynamic potentials of a finite system in a finite volume are analytic functions. This situation can be thought as unsatisfactory since on one hand the thermodynamic limit does not exist in nature, and on the other hand many experimental efforts are devoted to understand the link between mesoscopic and macroscopic systems, especially for spectacular properties like phase transitions wales .

A classification scheme valid for finite systems has been proposed by Grossmann grossmann  using the distribution of zeroes of the canonical partition sum in the complex temperature plane. Alternatively, it has been claimed that first order phase transitions in finite systems can be related to a negative microcanonical heat capacity gauss ; gross  or more generally to an inverted curvature of the thermodynamic potential as a function of an observable which can then be seen as an order parameter prl99 . These anomalies can also be connected to the general topology of the potential energy surface pettini . These very different approaches share the common viewpoint that at the mesoscopic or microscopic scale phase transitions are not smeared to become loosely defined smooth state changes, but can be analyzed in a rigorous way without explicitly addressing the thermodynamic limit. If this is certainly important to progress in the experimental study of phase transitions in finite systems, it is also true that the connection between these mesoscopic phenomena and the thermodynamics of the macroscopic world is not clear. This is especially annoying since the definitions of phase transitions in finite systems are not in general equivalent and different conclusions can be drawn out of the observation of a physical phenomenon depending on the theoretical framework used. In particular it can be shown deuteron  that any generic channel opening may lead to a convexity anomaly of the microcanonical entropy which does not necessarily survive to the bulk. In the same way the distribution of zeroes close to the real axis may suggest a first order phase transition even for systems where the transition temperature diverges in the thermodynamic limit isomers . In order to have a coherent picture of the physical meaning of a phase transition in a finite system it is then important to analyze in detail the conditions of persistency of the signals towards the bulk and to make a clear bridge between the various definitions.

We have recently proposed the possible bimodality of the probability distribution of observable quantities as a connection between these ideas topology . In this letter we further demonstrate that this definition is consistent with the Yang Lee theorem yanglee , $`i.e.`$ with the standard definition of first order phase transitions in the thermodynamic limit. Since our definition is a generalization of the previous ones grossmann ; gauss ; gross , the conditions of compatibility of the different signals and criteria proposed to identify first order phase transitions in finite systems can be explicitly worked out.

Let us consider that for each event $`i`$ we can measure a set of $`K`$ independent observables, $`b_k^{(i)},`$ which form a space containing one possible order parameter. We can sort events according to the results of the measurements $`\stackrel{}{b}^{(i)}\{b_k^{(i)}\}`$ and thus define a probability distribution $`P(\stackrel{}{b})`$. In the absence of a phase transition $`\mathrm{log}P(\stackrel{}{b})`$ is expected to be concave. An abnormal (e.g. bimodal) behavior of $`P(\stackrel{}{b})`$ or a convexity anomaly signals a first order phase transition topology .

This definition does not imply an a priori knowledge of the thermodynamic potential, $`i.e.`$ of the properties of the specific equilibrium reached in the experimental situation, but only requires the measurement of the event distribution. Moreover, it does not assume the a-priori knowledge of an order parameter; the order parameter can be extracted from the data themselves by looking at the best direction in the observable space to split the event cloud into two components. This topological definition is therefore easy to experimentally implement, and in fact it has already allowed to successfully recognize phase transitions in metallic clusters haberland  and in the multifragmentation of atomic nuclei lopezo .

Moreover this definition has a wider application domain than the previous ones grossmann ; gauss ; gross ; it can be used out of canonical or microcanonical equilibria, for open systems topology  and non-extensive statistics such as Tsallis equilibria next2001 . In particular it can be viewed as an extension of the definition based on the convexity anomaly of the microcanonical entropy gauss ; gross  or more generally of any thermodynamical Gibbs potential prl99 . Let us first discuss Boltzmann-Gibbs equilibrium which is obtained by maximizing the Shannon information entropy under the constraints of the various observables $`b_k`$ known in average. Then $`P(\stackrel{}{b})`$ can be written as

$$\mathrm{log}P_\stackrel{}{\lambda }(\stackrel{}{b})=\mathrm{log}\overline{W}(\stackrel{}{b})\underset{k=1}{\overset{K}{}}\lambda _kb_k\mathrm{log}Z_\stackrel{}{\lambda }$$
(1)
where $`Z_\stackrel{}{\lambda }`$ is the partition sum of this intensive ensemble controlled by the Lagrange multipliers $`\stackrel{}{\lambda }\{\lambda _k^{(i)}\}`$ and $`\overline{W}(\stackrel{}{b})=Z_0P_0(\stackrel{}{b})`$ is the partition sum of the statistical ensemble associated with fixed values $`\stackrel{}{b}`$ of all the observables. In the following we will call it the extensive ensemble using the word ”extensive” in the weak sense of ”observable” balian , meaning that observables are controlled on an event-by-event basis in the extensive ensemble. The two partition sums are related through the usual Laplace transform $`Z_\stackrel{}{\lambda }=𝑑\stackrel{}{b}\overline{W}(\stackrel{}{b})\mathrm{exp}(\stackrel{}{\lambda }\stackrel{}{b})`$. Eq.(1) clearly demonstrates that the convexity anomalies of the thermodynamical potential $`\mathrm{log}\overline{W}`$ can be traced back from $`\mathrm{log}P_\stackrel{}{\lambda }(\stackrel{}{b})`$. Indeed, the equation of states of the extensive ensemble can be related to the derivative of $`\mathrm{log}P_\stackrel{}{\lambda }(\stackrel{}{b})`$

$$\overline{\lambda }_k(\stackrel{}{b})\frac{\mathrm{log}\overline{W}(\stackrel{}{b})}{b_k}=\frac{\mathrm{log}P_\stackrel{}{\lambda }(\stackrel{}{b})}{b_k}+\lambda _k.$$
(2)
while the curvature matrix of $`\mathrm{log}\overline{W}`$ is the one of $`\mathrm{log}P_\stackrel{}{\lambda }(\stackrel{}{b})`$

$$C_{kk^{}}(\stackrel{}{b})\frac{\overline{\lambda }_k(\stackrel{}{b})}{b_k^{}}\frac{^2\mathrm{log}\overline{W}(\stackrel{}{b})}{b_kb_k^{}}=\frac{^2\mathrm{log}P_\stackrel{}{\lambda }(\stackrel{}{b})}{b_kb_k^{}}.$$
(3)
It should be noticed that these relations are valid for every set of Lagrange multipliers $`\stackrel{}{\lambda }`$ because in a finite system the probability distribution covers the entire accessible space of $`\stackrel{}{b}`$ and so contains the whole thermodynamics of the extensive ensemble. The only limitation can be a practical one for experiments or simulations: to accumulate enough statistics at every location.

Let us consider the simplified case of a unique direction $`b_kb`$ and so a unique Lagrange multiplier $`\lambda _k\lambda `$. If $`\overline{W}\left(b\right)`$ has an abnormal curvature, then the EOS $`\overline{\lambda }\left(b\right)=_b\mathrm{log}\overline{W}\left(b\right)`$ presents a back-bending. For this extensive ensemble where $`b`$ is the control parameter, coexistence can be defined as the region where $`\overline{\lambda }`$ is associated to three values of $`b`$ because of the anomalous curvature. The definition of a first order phase transition as the occurrence of a negative heat capacity gauss ; gross  is readily obtained if one identifies $`S=\mathrm{log}\overline{W}`$ as the microcanonical entropy, $`b`$ as the energy and $`\lambda =\beta `$ the inverse canonical temperature. The generalization to the EOS of any thermodynamic potential involving an extensive variable $`b`$ prl99  is also obvious. For example, a bimodal grand-canonical number of particle distribution is equivalent to a negative chemical susceptibility in the canonical ensemble, a bimodal distribution of magnetization to a negative magnetic susceptibility in the constant magnetization ensemble, a bimodal volume distribution in the isobar ensemble to a negative compressibility of the isochore ensemble, …

For $`\lambda `$ in the region of anomalous curvature the associated probability distribution presents two maxima and a minimum. In the intensive ensemble eq.(1) where $`\lambda `$ is controlled, coexistence is then signalled by the bimodality of the probability distribution and the value of $`\lambda `$ for which the two maxima have equal height is the first order transition point, analogous to the usual Maxwell construction. Therefore, the definitions of a first order phase transition through the occurrence of an abnormal curvature of the thermodynamical potential of the extensive ensemble, or through the presence of a bimodal event distribution in the associated intensive ensemble are strictly equivalent.

Let us now consider the definition grossmann  based on the zeroes of $`Z_\stackrel{}{\lambda }`$, the partition sum of the intensive ensemble eq.(1), in the complex $`\stackrel{}{\lambda }`$ plane. We shall now work out the relations between these zeroes and the concavity of $`P_\stackrel{}{\lambda }`$. Let us continue for simplicity with a single observable $`b_kb`$ and the associated Lagrange multiplier $`\lambda _k\lambda `$. Looking at the zeros of $`Z_\lambda `$ in the complex $`\lambda `$ plane, we first show that a bimodal distribution corresponds to a partition sum fulfilling the Yang Lee theorem in the double saddle point approximation topology . Then we will demonstrate the reciprocal: a distribution of zeroes which fulfills the Yang Lee theorem is associated to a bimodal probability distribution for $`b`$ .

Using equation (1) we see that the partition sum for a complex parameter $`\gamma =\lambda +i\eta `$ is nothing but the Laplace transform of the probability distribution $`P_{\lambda _0}\left(b\right)`$ for any arbitrary parameter $`\lambda _0`$ lee

$$Z_\gamma =𝑑bZ_{\lambda _0}P_{\lambda _0}\left(b\right)e^{\left(\gamma \lambda _0\right)b}𝑑bp_{\lambda _0}\left(b\right)e^{\overline{\gamma }b}$$
(4)
where $`\overline{\gamma }=\gamma \lambda _0`$ and we have defined $`p_\lambda \left(b\right)=Z_\lambda P_\lambda \left(b\right)`$. We can also write the partition sum at a complex point $`\gamma =\lambda +i\eta `$ as a Fourier transform of the probability distribution $`p_\lambda \left(b\right)`$

$$Z_\gamma =𝑑bp_\lambda \left(b\right)e^{i\eta b}$$
(5)
If $`P_\lambda \left(b\right)`$ is monomodal we can use a saddle point approximation around the maximum $`\overline{b}_\lambda `$ giving $`Z_\gamma =e^{\varphi _\gamma \left(\overline{b}_\lambda \right)}`$ , with

$$\varphi _\gamma \left(b\right)=\mathrm{log}p_\lambda \left(b\right)i\eta b+\frac{1}{2}\eta ^2\sigma ^2\left(b\right)+\frac{1}{2}\mathrm{log}\left(2\pi \sigma ^2(b)\right)$$
(6)
where $`\sigma ^2=_b^2\mathrm{log}p_{\lambda _0}\left(b\right)`$. However, if $`\overline{W}_{\lambda _0}(b)`$ has a curvature anomaly it exists a range of $`\lambda `$ for which the equation $`_b\mathrm{log}\overline{W}_{\lambda _0}=\lambda \lambda _0`$ has three solutions, one minimum and two maxima $`b_\lambda ^{(1)}`$and $`b_\lambda ^{(2)}`$ . Let us first split the probability distribution into two normal components

$$p_\lambda (b)=p_\lambda ^{(1)}(b)+p_\lambda ^{(2)}(b)$$
(7)
$`p_\lambda ^{(i)}(b)`$ being peaked at $`b_\lambda ^{(i)}.`$ The partition sum reads

$$Z_\gamma =e^{\varphi _\gamma ^{(1)}}+e^{\varphi _\gamma ^{(2)}}=2e^{\varphi _\gamma ^+}\mathrm{cosh}(\varphi _\gamma ^{})$$
(8)
where $`2\varphi _\gamma ^+=\varphi _\gamma ^{(1)}+\varphi _\gamma ^{(2)}`$ and $`2\varphi _\gamma ^{}=\varphi _\gamma ^{(1)}\varphi _\gamma ^{(2)}`$. We can now use a double saddle point approximation around the two maxima which will be valid close to thermodynamical limit lee  and we get $`\varphi _\gamma ^{\left(i\right)}=\varphi _\gamma (b_\lambda ^{(i)})`$according to eq. (6). The zeros of $`Z_\gamma `$ then correspond to

$$2\varphi _\gamma ^{}=i(2n+1)\pi $$
(9)
The imaginary part is given by $`\eta =(2n+1)\pi /\mathrm{\Delta }b`$ where $`\mathrm{\Delta }b=b_\lambda ^{(2)}b_\lambda ^{(1)}`$ is the jump in $`b`$ between the two phases. For the real part we should solve the equation $`\mathrm{}(\varphi _\gamma ^{})=0.`$ In particular, close to the real axis this equation defines a $`\lambda `$ which can be taken as $`\lambda _0`$. If the bimodal structure persists when the number of particles goes to infinity, the loci of zeros corresponds to a line perpendicular to the real axis with a uniform distribution as expected for a first order phase transition yanglee .

Let us now work out the necessary condition and show that a uniform distribution of zeroes perpendicular to the real axis with a density linearly increasing with the number of particles $`N`$ implies a bimodal probability distribution. Let us denote zeros as

$$\gamma _n=\lambda _0+i(2n+1)\frac{\pi }{N\delta }$$
(10)
such that the interval $`2\pi /\delta `$ contains $`N`$ uniformly distributed zeroes in agreement with the unit circle theorem yanglee . Since all the zeros of $`Z_\gamma `$ are periodically distributed one can define an analytic function $`f(\gamma )`$ such that

$$Z_\gamma =2\mathrm{cosh}\left(\left(\gamma \lambda _0\right)\frac{N\delta }{2}\right)f(\gamma )$$
(11)
Indeed since $`Z_\gamma `$ is analytic, $`f(\gamma )`$ could in principle diverge only on the zeroes of $`Z_\gamma `$; in these points however eq.(11) shows that $`f`$ is proportional to $`_\gamma Z_\gamma `$, $`i.e.`$ it is analytic. On the line of zeroes, relation (11) reduces to $`Z_{\lambda _0+i\eta }=2\mathrm{cos}\left(\eta N\delta /2\right)f(\lambda _0+i\eta ).`$ Using the fact that the partition sum along the line of zeroes is the Fourier transform of the reduced probability distribution at the transition point $`\lambda _0`$ (eq.(5), we can use the inverse Fourier transform to get the distribution $`p_{\lambda _0}\left(b\right)`$ at the transition point

$$p_{\lambda _0}\left(b\right)=\frac{1}{2\pi }𝑑\eta Z_{\lambda _0+i\eta }e^{i\eta b}$$
(12)
Equation (12) can then be rewritten as

$$p_{\lambda _0}\left(b\right)=g_{\lambda _0}\left(b+\frac{N\delta }{2}\right)+g_{\lambda _0}\left(b\frac{N\delta }{2}\right)$$
(13)
where

$$g_\lambda (b)=\frac{1}{2\pi }𝑑\eta f(\gamma )e^{i\eta b}$$

is a distribution. Indeed if we compute $`p_\lambda (b)`$ a little above $`\lambda _0`$, $`\lambda =\lambda _0+\epsilon `$ we get

$$p_{\lambda _0+\epsilon }\left(b\right)=e^{\epsilon N\delta /2}g_{\lambda _0+\epsilon }\left(b+\frac{N\delta }{2}\right)+e^{\epsilon N\delta /2}g_{\lambda _0+\epsilon }\left(b\frac{N\delta }{2}\right)$$
(14)
For large $`N`$ only the first term survives

$$p_{\lambda _0+\epsilon }\left(b\right)e^{\epsilon N\delta /2}g_{\lambda _0+\epsilon }\left(b+\frac{N\delta }{2}\right)p_{\lambda _0+\epsilon }^{\left(1\right)}\left(b\right)$$
(15)
which is the first term in the distribution at the transition point eq.(12). If we conversely compute $`p_\lambda (b)`$ a little below the transition point we get

$$p_{\lambda _0\epsilon }\left(b\right)e^{\epsilon N\delta /2}g_{\lambda _0\epsilon }\left(b\frac{N\delta }{2}\right)p_{\lambda _0\epsilon }^{\left(2\right)}\left(b\right)$$
(16)
which is the second term in the transition point distribution eq.(12). At the transition point $`p_{\lambda _0}\left(b\right)`$ is the sum of two shifted (identical) distributions, $`p_{\lambda _0}^{\left(1\right)}`$ and $`p_{\lambda _0}^{\left(2\right)}`$. If $`g_\lambda `$ is a normal ($`i.e.`$ monomodal) distribution the central limit theorem guarantees that for a large number of particles its width will scale as $`\sqrt{N}`$, $`i.e.`$ will grow slower than the distance between the two peaks that scales as $`N`$. This implies that (for not too small $`N`$) $`p_{\lambda _0}\left(b\right)`$ is bimodal. This stays true if $`g_\lambda `$ itself is bimodal (or multimodal), with the only difference that the total distribution $`p_{\lambda _0}\left(b\right)`$ will then represent the coexistence between more than two phases. Let consider the simplest case of a bimodal structure for $`p_{\lambda _0}\left(b\right)`$, close to the thermodynamical limit; a little before the transition point only the lowest peak of the $`b`$ distribution is present, while a little above only the second one remains. At the transition both are present and the most probable $`b`$ jumps from $`b_<=\overline{b}+N\delta /2`$ to $`b_>=\overline{b}N\delta /2`$ where $`\overline{b}`$ is the maximum of $`g_{\lambda _0}\left(b\right).`$ $`\delta `$ represents the discontinuity of the variable $`b`$ per particle. If $`b`$ is the energy $`\delta `$ is the latent heat per particle.

The important difference between finite systems and the thermodynamic limit is that in the latter case the probability distribution is bimodal only at the transition point $`\lambda _0`$. On the other hand in finite systems, in an interval $`\mathrm{\Delta }\lambda `$ of non-zero measure around the transition point $`\lambda _0`$ the two phases coexist, $`i.e.`$ the distribution is bimodal, each peak being associated with a phase having a finite probability of occurrence. Eq.(14) shows that this interval is the larger the smaller is $`\delta `$, $`i.e.`$ the lower is the density of zeroes. This extension of the phase transition point to a region in which the intensive ensemble presents a bimodal distribution makes a direct measurement of phase coexistence possible lopezo , contrary to the common belief that phase transitions would only be loosely defined out of the thermodynamic limit. Moreover, the fact that the distribution is non zero between the two maxima $`b_<`$ and $`b_>`$ in the case of small systems (see eq.(14)) implies that microstates can be accessed that do not exist at the thermodynamic limit in the intensive ensemble. These states are specific of the coexistence region and can lead to spectacular phenomena as negative heat capacity gauss  and negative compressibility prl99  in the extensive ensemble.

This demonstration can be extended to account for more important finite size effects lee  when zeroes $`\gamma _n`$ are still periodic but not yet uniformly distributed on a straight line lee ; grossmann . Indeed, it is important to remark that the distribution of zeroes is periodic in the imaginary direction before the thermodynamic limit. If for example we consider the grancanonical ensemble, the Lagrange parameter $`\lambda \alpha \beta \mu `$ represents the logarithm of the fugacity $`z=e^\alpha `$ and the observable is the number of particles $`bN`$. The partition sum $`Z_\alpha e^{\alpha N}`$ in the complex $`\alpha `$ plane is periodic for any arbitrary value of $`N`$. In the more general case, this will not be true for very small systems since the extensive variables $`b_k`$ are not in general proportional to the number of constituents of the system. However, if the forces are short-ranged and $`N`$ is sufficiently high, $`b_kN`$, which guarantees the periodicity of $`Z_\lambda `$ in the complex $`\lambda `$ plane. This periodicity is a constraint on the distribution of zeroes.

To take into account finite size distortions, we can introduce a transformation

$$m(\gamma )=\gamma \mathrm{\Delta }m(\gamma )$$
(17)
which is the identity on the real axis and which maps the zeroes $`\gamma _n`$ of the partition sum onto a uniform density perpendicular to the real axis

$$m(\gamma _n)=\lambda _0+i(2n+1)\frac{\pi }{N\delta }$$
(18)
As before we can introduce an analytic function $`\overline{f}(\gamma )`$ such that

$$Z_\gamma =2\mathrm{cosh}\left(\left(m\left(\gamma \right)\lambda _0\right)\frac{N\delta }{2}\right)\overline{f}(\gamma )$$
(19)
Using this expression in the inverse Fourier transform eq.(12) the probability distribution at the transition point reads

$$p_{\lambda _0}(b)=\overline{g}_{\lambda _0}^+(b+\frac{N\delta }{2})+\overline{g}_{\lambda _0}^{}(b\frac{N\delta }{2})$$
(20)
where

$$\overline{g}_{\lambda _0}^\pm (b)=N_0𝑑\eta \overline{f}(\lambda _0+i\eta )e^{\mathrm{\Delta }m(\lambda _0+i\eta )N\delta /2}e^{i\eta b}$$
(21)
If a phase transition is present at the thermodynamic limit, according to the Yang-Lee theorem $`\mathrm{\Delta }m`$ goes to zero faster than $`N`$ such that the two $`\overline{g}_{\lambda _0}^\pm `$ converge to the same distribution $`g_{\lambda _0}`$, and $`p_{\lambda _0}(b)`$ tends to eq.(13). Following the discussion of eq.(20), this means that a deformed distribution of zeroes that converges to a straight line of equally spaced zeroes perpendicular to the real axis gives rise to a bimodal probability distribution function.

In finite systems, if $`\mathrm{\Delta }m`$ is not zero, the sum (20) may or may not present a concavity anomaly depending on the actual properties of $`\mathrm{\Delta }m`$. This is true irrespectively of the fact that the straight distribution of zeroes eq.(10) is asymptotically reached ($`i.e.`$ a first order phase transition exists in the bulk) or not. In this case, the physical phenomenon could be classified differently with the zeroes of the partition sum grossmann  or with the concavity criteria gross ; prl99 . Indeed, we have demonstrated their equivalence only for large systems.

In conclusion, our topological definition of phase transitions through the bimodality of the probability distribution function, generalizes the definition based on the microcanonical entropy gauss ; gross  and is equivalent to the definition based on convexity anomalies of generic thermodynamic potentials prl99  for any number of particles. We have demonstrated that a bimodality is equivalent to the Yang Lee definition of phase transitions close to and at the thermodynamical limit. A uniform distribution of zeroes perpendicular to the real axis of a Lagrange multiplier $`\lambda `$ at the position $`\lambda _0`$ implies that the associated distribution of events is bimodal for this $`\lambda _0`$ , and conversely a bimodal distribution at a given $`\lambda _0`$ generates a uniform distribution of zeroes with a real part $`\lambda _0`$ and a regularly spaced imaginary part. This means that the occurrence of bimodalities can be considered as a valid generalization of the concept of phase transition to systems of any size. This topological definition is extremely powerful since it gives access to the order parameters defined as the observables for which the distribution is bimodal. We have demonstrated that this definition agrees with the definition based on the zeroes of $`Z`$ grossmann  when the bimodal distribution can be approximated by two gaussians, i.e. in the double saddle point approximation. This means that the different definitions are coherent for $`N`$ big enough that fluctuations around the two phases are at the gaussian level, but they can differ for smaller systems. The compatibility of the different definitions at the thermodynamic limit means that the great effort done by the different research groups to recognize and classify phase transitions in finite systems grossmann ; gross ; noi  is not only rigorous but also consistent with macroscopic thermodynamics.


format: text

text:
# Formation of clean dimers during gas-source growth of Si(001)

## I Introduction

Gas-source growth of Si(001) using hydrogen-based precursors (such as SiH4, silane, and Si2H6, disilane) is of great scientific and technological interestM.J.Bronikowski et al. (1993); Y.Wang et al. (1994); J.H.G.Owen et al. (1997a, b) — in particular, hydrogen can act as an effective surfactant, and has been shown to reduce roughness and intermixing during growth of Ge/Si alloys and pure Ge on Si(001)N.Ohtani et al. (1994). Understanding the reactions that occur and the intermediate structures that are formed during this growth will enable greater control of surface and interfaces during growth. STM observations of the growth of Si(001) from disilane, both at room temperature following annealsM.J.Bronikowski et al. (1993); Y.Wang et al. (1994) and at elevated temperatureJ.H.G.Owen et al. (1997a, b), along with careful electronic structure calculationsJ.H.G.Owen et al. (1997a, b); D.R.Bowler and C.M.Goringe (1996) have mapped out the growth pathway. A key observation in this pathway is that the islands formed during growth are clean, while the substrate remains covered with a certain amount of hydrogenJ.H.G.Owen et al. (1997a). The fundamental building block in gas-source growth is the clean ad-dimer (as opposed to solid-source growth, where fast-moving ad-atoms are keyB.S.Swartzentruber (1996)); yet, the mechanism to form such clean dimers from the hydrogenated dimers that occur naturally during gas-source growth is unknown. In particular, they are observed to form at 450K while desorption from the monohydride phase occurs at 790K, indicating that their formation must be completely different to the desorption of hydrogen from the monohydride phase. In this paper, we present a first-principles investigation of the mechanism for formation of clean ad-dimers from hydrogenated dimers, with the aim of explaining how these form at a comparatively low temperature.

Disilane (which is used in preference to silane as it decomposes more easily) adsorbs on Si(001) as SiH3 (which soon breaks down to form SiH2) or SiH2M.J.Bronikowski et al. (1993), sometimes with accompanying hydrogen. These SiH2 groupsD.R.Bowler and C.M.Goringe (1996) start to diffuse at 400–500KJ.H.G.Owen et al. (1997a). When two groups are on adjacent dimer rows, they react to form a hydrogenated ad-dimer (that is, an ad-dimer with both dangling bonds saturated with hydrogen, illustrated in Fig. 1 (a)) over the trench between the dimer rowsJ.H.G.Owen et al. (1997a). This then decomposes to form clean ad-dimers and hydrogen on the surface at around 450KM.J.Bronikowski et al. (1993); Y.Wang et al. (1994); J.H.G.Owen et al. (1997a), via a pathway to be investigated in this paper. A hydrogenated ad-dimer (which is the starting point) is illustrated in Fig. 1 (a), along with a partially hydrogenated ad-dimer (the result of the first part of the pathway) in Fig. 1 (b) and a clean ad-dimer (the final point) in Fig. 1 (c). Once formed, the clean ad-dimers diffuseC.M.Goringe and D.R.Bowler (1997) and form a square feature, which is believed to be the precursor to dimer stringsJ.H.G.Owen et al. (1997c), followed by short strings of dimersY.Wang et al. (1994) which later increase to form larger islandsY.Wang et al. (1994); J.H.G.Owen et al. (1997b).

The calculations to be presented are based on density functional theory (DFT) in the generalized gradient approximation (GGA), with a plane wave basis set and pseudopotentials. We have searched for possible pathways both by applying constraints to specific atoms (for instance constraining a hydrogen to lie in a given plane) and by using the nudged elastic band technique (NEB)H.Jónsson et al. (1998); Henkelman et al. (2000a) which allows accurate determination of reaction barriers given an initial approximation to a pathway. One key result is that the dehydrogenation proceeds via a meta-stable intermediate state (this is discussed fully in Section III and illustrated in Fig. 2).

The rest of the paper is organised as follows: the next section gives details of the computational techniques used; this is followed by a detailed discussion of the structure of the meta-stable state which plays a key role in the dehydrogenation; the diffusion pathways are then presented, looking at the mechanism for both hydrogens, followed by a conclusion section.

## II Computational Details

The theory underlying DFTP.Hohenberg and W.Kohn (1964); W.Kohn and L.J.Sham (1965) and their application to electronic structure calculations have been extensively reviewedR.O.Jones and O.Gunnarsson (1989), as has the use of pseudopotential and plane-wave techniquesM.C.Payne et al. (1992). The calculations in this paper were performed using the VASP codeG.Kresse and J.Furthmüller (1996), using the standard ultra-soft pseudopotentialsD.Vanderbilt (1990) that form part of the code. The approximation we use for exchange-correlation energy is the generalised-gradient approximation (GGA) due to Perdew and Wang (PW91) Y.Wang and J.P.Perdew (1991); J.P.Perdew et al. (1992). We chose the GGA rather than the local density approximation (LDA) rather deliberately. As the barriers that we will be calculating are sensitive to bonding and stretched bonds, and the GGA is known to be rather more accurate in these situations (LDA generally overbinds), we considered its use to be essential for this work.

We use periodic boundary conditions, as is standard for DFT calculations with plane-waves, and we therefore used a periodic slab for the surface, with a vacuum layer between the slabs. Our simulations were performed within a unit cell two dimers long and two dimer rows wide, with five layers of Si (the bottom of which was terminated in hydrogen and constrained to remain fixed in bulk-like positions). The vacuum gap of 6.9 Å is equivalent to five atomic layers of Si, and provides sufficient isolation between vertical periodic images. We used a plane wave cutoff of 150 eV and a $`2\times 2\times 1`$ Monkhorst-Pack k-point mesh. All these parameters were tested, and found to converge energy differences to better than 0.01 eV. The system contains an even number of electrons, but has various saddle points which might involve unpaired electrons, so we checked the effect of performing spin-polarised calculations for these points. The effect was found to be negligible (both on energies and geometries) and so was not used in the calculations.

To investigate the diffusion pathways, we used two techniques: first, constraining the diffusing hydrogen to lie in a particular plane, and calculating static energies for different locations of the hydrogen; second, the Nudged Elastic Band methodH.Jónsson et al. (1998); Henkelman et al. (2000a). This second method requires the simultaneous relaxation of a number of images of the system, which can be done in parallel. However, this has the potential to become extremely computationally intensive, which is why we chose to use the smallest realistic unit cell (with four dimers in the surface).

The initial exploration of the system used the first method (static calculations, constraining the hydrogen). It was using this method that the meta-stable state (discussed in Sec. III) was found, and it is unlikely that it would have been found using the NEB without significant effort (for instance performing simulated annealing on the initial images of the system) or using a more complicated technique such as the dimer methodHenkelman and Jónsson (1999). The diffusion barriers presented in the paper were all calculated using a variant of the original technique which actively seeks the saddle point, the climbing image NEBHenkelman et al. (2000b), with 8 images relaxed in the chain.

While we have calculated the diffusion barriers, we have not calculated attempt frequencies, which have been assumed to be $`10^{13}\mathrm{sec}^1`$, typical for diffusion processes. While DFT-GGA is sufficiently accurate to calculate reaction barriers to within 0.1 eV, it is not able to predict attempt frequencies accuratelyP.E.Blöchl et al. (1990); for instance, in previous work on solid-source growth of Si(001), it was shown that a factor of five in the attempt frequency was required to understand the results, but was not accurately predictedA.P.Smith and H.Jónsson (1996).

## III The Meta-stable Intermediate

The lowest energy diffusion pathway, and the only one which has an energy barrier which is in line with the temperature at which the dehydrogenation is observed to occur, proceeds via a meta-stable intermediate. This is an unusual and rather important structure, and will be discussed in detail in this section. It is illustrated in Fig. 2, and should be contrasted with the hydrogenated ad-dimer illustrated in Fig. 1 (a).

The atom labelled ‘A’ in Fig. 2 is one of the substrate dimers, to which the hydrogenated (and clean) ad-dimer is bonded; ‘B’ is a second layer atom in the substrate, to which the substrate dimer is normally bonded (for instance in Fig. 1 (a)); and ‘C’ is the ad-dimer atom itself, which is now clean (having started hydrogenated). The hydrogen is now bonded to ‘A’, which has broken its bond to ‘B’, while ‘C’ has formed a bond to ‘B’ (not easily seen, owing to the geometry).

In terms of the bonding of the atoms, the atoms A, B and C are all saturated as they are in Fig. 1 (a) — the bonding has merely cycled around (so that the A–B bond is now an A–H bond and a B–C bond, while the C–H bond is now a C–B bond). It is this saturation that gives the structure its stability. While some of the bond angles are rather strained (in particular, the bonds associated with B and C make $`60^{}`$ angles) the bond lengths are all close to equilibrium, and there are no further broken bonds, leading to an energy difference of 0.57 eV relative to the starting point, but no more. It is interesting to note that there are other structures where $`60^{}`$ bond angles are found during growth of Si(001), which also exhibit stability which might seem counter-intuitiveD.R.Bowler and C.M.Goringe (1996).

In terms of the formation of this structure, as we shall see in the next section, there is not a large barrier. The H never has to move a long distance from either A or C, leading to relatively strong bonds being present at all times; the second-layer atom B moves up slightly; and while the substrate dimer atom A and the ad-dimer atom C do move up and down respectively, they do this gradually while maintaining their bonding. It is this relatively small perturbation on the overall structure, and the ease with which it is reached, which allows the formation of this state, and gives it its importance.

## IV Diffusion Pathways

In this section, we describe the diffusion pathways that we have explored with DFT calculations. For simplicity, and because it is likely to be physically realistic, we allow the hydrogens on the ends of the dimer to diffuse off independently — i.e. we consider the diffusion off one end of the ad-dimer while the other hydrogen remains on the ad-dimer. Then we allow the remaining hydrogen to diffuse off the now partially-hydrogenated ad-dimer onto the substrate. In order to avoid the complications of spin and half-filled bands, we maintain both hydrogens in the unit cell at all times (the first hydrogen to diffuse off stays on the substrate, illustrated in Fig. 1(b)). The three stable points of the process (fully hydrogenated ad-dimer, partially hydrogenated ad-dimer with a hydrogen on the substrate and clean ad-dimer with both hydrogens on the substrate) are illustrated in Fig. 1. The atomic positions during the diffusion pathways are presented below in an aggregated form (due to space constraints): only the position of atoms which move significantly are shown. All of the atomic structures at each step in all the diffusion pathways are available elsewhereDeH .

The experimental data that we are comparing against comes from two separate experiments: first, where the Si(001) surface was exposed to a dose of disilane (Si2H6), annealed at different temperatures for different times, and then observed at room temperature in STMM.J.Bronikowski et al. (1993); Y.Wang et al. (1994); second, where an elevated-temperature STM was used to observe the results of dosing with disilane in real time at different temperaturesJ.H.G.Owen et al. (1997a, b). The results of both these types of experiment are identical: around 450K, clean, non-rotated dimers are formed over the trench between dimer rows. In other words, the monohydride dimers lose their hydrogen to the substrate in a matter of minutes at this temperature (for instance, an anneal to 470K for two minutes led to the dehydrogenation of all ad-dimersY.Wang et al. (1994)). Assuming an attempt frequency of $`10^{13}\mathrm{sec}^1`$ and a successful dehydrogenation rate of 1/60 Hz, we obtain a barrier of 1.28 eV. This changes by about 0.03 eV if the rate is doubled or halved, giving us a good estimate of the likely reaction barrier.

### IV.1 The First Hydrogen

There are two diffusion paths considered for the first hydrogen diffusing off the ad-dimer: a direct diffusion path; and diffusion via the meta-stable considered in Sec. III and shown in Fig. 2. We will discuss these separately, starting with the direct diffusion, and then contrast their results.

The diffusion barrier for direct diffusion is shown in Figure 3, with open circles and long dashes. The barrier is 1.93 eV, which is extremely high; the reason for this can be seen from the atomic positions, which are illustrated in Figure 4. At the saddle point, the ad-dimer bond is extended greatly (from 2.51 Å at the start to 2.82 Å) while the bond from the hydrogen to the ad-dimer is more extended (from 1.51 Å at the start to 1.85 Å). Inspecting the charge density, it is clear that the ad-dimer remains bonded (though weakly) and that the H has made a weak bond to the substrate dimer (which is 2.32 Å away) as well as maintaining a slightly weakened bond to the ad-dimer. It is this lengthening and weakening of bonds at the saddle point that causes the high barrier. Assuming Arrhenius behaviour and an attempt frequency of $`10^{13}`$ Hz, we find a hopping rate of $`10^{10}\mathrm{sec}^1`$ at 450K, which is many orders of magnitude below the observed rate.

The diffusion barrier into and out of the meta-stable state is also shown in Fig. 3, with open squares and a dot-dashed line. The barrier from the start to the meta-stable state is 1.24 eV, while the barrier from the meta-stable state to the end is 0.56 eV (and the reverse path, from the meta-stable state to the start is 0.66 eV).

The pathway from the starting position to the meta-stable state (shown in Fig. 5) involves considerable rearrangement: first, the hydrogen inserts into the bond between the ad-dimer (labelled ‘C’ in Fig. 2) and the substrate dimer (labelled ‘A’ in Fig. 2); at the saddle point, the hydrogen is 1.65 Å from the ad-dimer, and 1.94 Å from the substrate dimer (compared to an equilibrium distance of 1.51 Å), while the distance between the ad-dimer and the substrate dimer is 3.10 Å (compared to 2.48 Å at the start); second, after the hydrogen has transferred to the substrate dimer, the ad-dimer bonds to a second layer atom in the substrate (labelled ‘B’ in Fig. 2); third, the substrate dimer bonds back to the ad-dimer and breaks its bond to the second layer atom in the substrate. The first part of this rearrangement is the area where most of the energy change happens: the energy actually falls by about 0.6 eV during the second and third parts of the rearrangement.

The pathway from the meta-stable state to the end position (shown in Fig. 6) is much simpler, involving only the movement of the hydrogen from one end of the substrate dimer to the other, while the substrate dimer ‘C’ reforms its bond to the second layer atom ‘B’. At the saddle point, the hydrogen is 1.72 Å from the substrate atom and 2.05 Å from the end atom.

The barrier of 1.24 eV from the starting point to the meta-stable state fits extremely well with the observed temperature behaviour: at 450K with an attempt frequency of $`10^{13}\mathrm{sec}^1`$, it would correspond to a hopping rate of 0.044 Hz, or one hop every 23 seconds. But this is only into the meta-stable state, and there are two low energy paths out of that. The hopping rate from the meta-stable state to the end state is $`4\times 10^6`$ Hz, while from the meta-stable state to the start state is $`2\times 10^5`$ Hz, so that only 10% of meta-stable states would return to the starting point. We also expect that the equilibrium populations of the start and end states would differ, since the end state is 0.53 eV lower than the start (roughly, at 450K, we would expect a relative population about $`10^6`$ times higher in the lower state). There is also the question of whether the hydrogen could return, via the meta-stable state, from the end to the start. The barrier from the end point back to the meta-stable state it is 1.66 eV, making it extremely unlikely that the hydrogen would return to the meta-stable state (and even if it did, it would be ten times more likely to drop back to the end state than to return to the start state). Clearly, it is the low barrier from the starting state to the meta-stable state that allows the first part of the dehydrogenation of the ad-dimer to proceed, and the energy difference between the start and end points, as well as the high barrier out of the end state that makes the reaction effectively irreversible.

### IV.2 The Second Hydrogen

Once the first hydrogen has diffused off the ad-dimer, we retain it on the substrate, as shown in Figure 1(b). This is computationally convenient (as it maintains a filled set of bands) but also physically reasonable: hydrogen does not begin diffusing along the dimer rows on Si(001) at an appreciable rate until about 550K with a barrier of 1.68 eVJ.H.G.Owen et al. (1996).

As with the first hydrogen, the second hydrogen can diffuse either directly, or via a meta-stable state, which is exactly equivalent to the meta-stable state for the first hydrogen (shown in Fig. 2) and therefore not illustrated here. As before, we will discuss these results separately, starting with the direct diffusion.

The diffusion barrier for direct diffusion is shown in Figure 7, plotted with open circles and dashes. The shape is much broader than for the first hydrogen’s direct path, with a lower barrier of 1.59 eV. The reason for this can be seen in part in the atomic positions, which are shown in Figure 8. This is a little more confusing than previous plots, as the positions both of the hydrogen and the ad-dimer atom to which it is bonded have been plotted. As the hydrogen moves across towards the substrate dimer (the end point), the bond between the silicon atoms in the ad-dimer breaks, with the atom that the hydrogen is bonded to following the hydrogen as it diffuses. At the saddle point, the hydrogen is 1.67 Å from the ad-dimer and 2.03 Å from the substrate dimer, while the distance between ad-dimer atoms is 4.48 Å. Beyond this point, the hydrogen transfers to the substrate, and the ad-dimer reforms slowly.

The bond in the clean or partially clean ad-dimer is not as strong as the other bonds to the substrate, which explains why the energy cost for breaking it is relatively small, and why this pathway is followed in contrast to the pathway for the first hydrogen. Even with the reduced barrier, the hopping rate at 450K will be $`4\times 10^6\mathrm{sec}^1`$, which is still far too low to be consistent with the experimental observations.

The diffusion barrier into and out of the meta-stable state is also shown in Figure 7, plotted with open squares and dot-dashed lines. The barrier from the start to the meta-stable state is 1.14 eV, while the barrier from the meta-stable state to the end is 0.58 eV (and the reverse path, from the meta-stable state to the start is also 0.58 eV). The atomic positions are almost identical to those for the first diffusion (shown in Figs. 5 and 6) and are not shown (though these figures, and many other pieces of supplementary information such as animations of the processes can be found elsewhereDeH ).

The barrier of 1.14 eV is 0.1 eV lower than the barrier for the first hydrogen, suggesting that once the first hydrogen has diffused off the ad-dimer, the second will follow slightly more quickly; it is still in excellent agreement with observed experimental behaviour. The barriers from the meta-stable state to the start and end states are now identical, meaning that 50% of meta-stable states will return to the starting state. However, the end state is 0.48 eV lower in energy than the start, so that (as before with the first hydrogen) we would expect the population in thermal equilibrium at 450K to be about $`10^6`$ times higher in the end state than the start state. The barrier from the end state back to the meta-stable state is 1.63 eV, which again makes the reaction effectively irreversible. Of course, the clean ad-dimer can also diffuse away along the trench between dimer rows, with a barrier of 1.15 eVC.M.Goringe and D.R.Bowler (1997), which would make the reforming of the hydrogenated ad-dimer impossible.

## V Conclusions

We have presented ab initio calculations, modelling the diffusion of hydrogen off a hydrogenated ad-dimer, which is a key stage in gas-source growth of Si(001). We have shown that the diffusion proceeds via a meta-stable intermediate, and that the energy barriers calculated (1.24 eV for the first hydrogen and 1.14 eV for the second hydrogen) are in excellent agreement with temperatures at which these features are observed in experiment.

We have used the climbing image nudged elastic band method to find the diffusion barriers, and have found it to be extremely effective, particularly for the direct diffusion which was difficult to model simply by picking a single constraint. However, the problem of exploring phase space is still a difficult one, as the existence of the meta-stable state (which was discovered through application of a single constraint) shows. There are techniques for exploring energy surfaces, such as the dimer methodHenkelman and Jónsson (1999) and variants on hyperdynamicsA.Voter (1997), but there is still a large amount of work to be done in this field.

###### Acknowledgements.

We gratefully acknowledge useful discussions with Dr J.H.G.Owen, Dr. C.M.Goringe and Prof. M.J.Gillan. DRB thanks the Royal Society for funding through a University Research Fellowship. Calculations were performed at the HiPerSPACE Centre, UCL which is partly funded by JREI grant JR98UCGI. We also acknowledge the Jonsson group in University of Washington for making their implementation of the NEB and climbing image NEB within VASP publicly available.


text:
# An efficient quantum algorithm for the one-dimensional Burgers equation

(August 14, 2002)

## Abstract

We analyze one-dimensional classical and quantum microscopic lattice-gas models governed by a lattice Boltzmann equation at the mesoscopic scale, achieved by ensemble averaging over microscopic realizations. The models are governed by the Burgers equation at the macroscopic scale, achieved by taking the limit where the grid size and time step both approach zero and by performing a perturbative Chapman-Enskog expansion. The quantum algorithm exploiting superposition and entanglement is more efficient than the classical one because the quantum algorithm requires less memory. Furthermore, its viscosity can be made arbitrarily small.

Quantum Computing, Quantum Lattice Gas, Burgers equation  preprint: Draft Revision 1.2

Here we present the simplest example where a quantum computer is demonstrably more efficient at numerically predicting the time-dependent solutions of an important nonlinear one-dimensional partial differential equation, the classical Burgers equation

$$_tu+u_xu=\nu _x^2u,$$
(1)
a simplified model of turbulence and shock formation with flow field $`u(x,t)`$, sound speed $`c_s`$, and kinematic viscosity $`\nu `$. We shall derive (1) as the general effective field theory describing the large-scale behavior of microscopic one-dimensional lattice gas models with two particles per site. The occupation probabilites for the two particles at position $`x`$ and at time $`t`$ are denoted by $`p_+(x,t)`$ and $`p_{}(x,t)`$, respectively. The mesoscopic kinetic transport dynamics is governed by the lattice Boltzmann equation

$$p_\pm (x\pm \delta x,t+\delta t)=p_\pm (x,t)\pm \mathrm{\Omega }(p_+,p_{}),$$
(2)
where $`\mathrm{\Omega }`$ denotes the collision term, a nonlinear function of $`p_+`$ and $`p_{}`$. The particular functional form of $`\mathrm{\Omega }`$ depends on the model type, either classical or quantum mechanical. For any $`\mathrm{\Omega }`$, the model (2) always conserves particle number density defined as $`\rho p_++p_{}`$. The propagation speed of particles is the ratio of the lattice cell size to the time step interval, $`c=\frac{\delta x}{\delta t}`$. In these models, the macroscopic flow field is

$$u(x,t)=c\left(\rho (x,t)1\right).$$
(3)
The total classical computational complexity is the product of the lattice size $`N_x`$, the ensemble of size $`N`$, and local resources of size $`\varrho `$

$$C_{\text{cl}}NN_x\varrho .$$
(4)
The quantity $`\varrho `$ is the amount of resources needed to encode $`\mathrm{\Omega }`$ for the local microscopic collisions. Since we are modeling the time evolution of a classical system, the value of $`u(x,t)`$ must be known everywhere at every time step whence the need for continual measurement in any quantum model. Therefore, the memory load factor for the lattice is still $`N_x`$ and not $`\mathrm{log}_2N_x`$ in the quantum case yepez-ijmpc00a . Furthermore, because of wave function collapse by Von Neuman projective measurement, ensemble averaging is required just as in the classical case, so a lower bound for the quantum computational complexity $`C_{\text{qu}}`$ is

$$C_{\text{qu}}NN_x\mathrm{log}_2\varrho .$$
(5)
For some quantum computing technologies, such as spatial nuclear magnetic resonance quantum computing yepez-pravia-cpc2001 ; yepez-pravia-pre2002 , $`N`$ may be counted as order unity overhead. Relating $`\varrho `$ to the minimum number of bits per node needed to encode (2), as a parsimonious demonstration of (4) and (5), it is known (1) can be modeled classically with 3 bits per node whereas we prove quantum mechanically only 2 qubits per node is required.

Multiplying (1) by $`u`$ and integrating over all space, with periodic boundaries, gives a relation for energy conservation where the time-rate of change of the turbulent kinetic energy density $`_t(\frac{u^2}{2})`$ is balanced by the viscous dissipation $`\epsilon \nu \left(\frac{u}{x}\right)^2\frac{u_{}^3}{}`$, where $``$ is the characteristic scale of the largest feature in the flow field and $`u_{}`$ is the standard deviation of the turbulent kinetic energy or eddy velocity at that scale.

The flow velocity, the kinematic viscosity, and the viscous dissipation quantities have the dimensions: $`\left[u\right]=\frac{L}{T}`$, $`\left[\nu \right]=\frac{L^2}{T}`$, and $`\left[\epsilon \right]=\frac{L^2}{T^3}`$. The dissipation scale $`\lambda \left(\frac{\nu ^3}{\epsilon }\right)^{\frac{1}{4}}`$ is the smallest spatial scale where macroscopic effective field theory (1) is physically applicable and the smallest physical velocity at the dissipation scale is the dissipation-scale velocity $`u_\lambda \left(\nu \epsilon \right)^{\frac{1}{4}}=\frac{\nu }{\lambda }`$.

The Reynolds number characterizing the fluid’s nonlinearity is $`\text{Re}\frac{}{\lambda }\frac{u_{}}{u_\lambda }=\frac{u_{}}{\nu }`$, and using $`\epsilon =\frac{u_{}^3}{}`$ to eliminate $``$, we also have $`\text{Re}=\frac{u_{}^4}{\nu \epsilon }=\left(\frac{u_{}^2}{u_\lambda ^2}\right)^2=\left(\frac{}{\lambda }\right)^{\frac{4}{3}}`$. The computational complexity of numerically modeling fluidic behavior at the macroscopic scale can be expressed as a function of Re. First, the number of grid points $`N_x`$ to sufficiently resolve the flow field down to the dissipation scale is $`N_x=\frac{}{\lambda }=\text{Re}^{\frac{3}{4}}`$. Second, we get $`N`$ by physically limiting $`|\delta p_\pm |`$. Since, the occupation probabilities are measured by ensemble averaging over $`N`$ independent microscopic realizations, $`|\delta p_\pm |\frac{1}{\sqrt{N}}`$, due to either classically stochastic shot noise or quantum mechanically stochastic projective measurement. Following Orszag and Yakhot orszag-86 , the value of the statistical fluctuation $`\delta u(x,t)`$ of the numerical flow field must be much much less than $`u_\lambda (x,t)`$ of the macroscopic effective field theory (1). Using (3), we have $`u_\lambda c(\delta p_++\delta p_{})\frac{c}{\sqrt{N}}`$. Hence, the ensemble size $`N\frac{c^2}{u_\lambda ^2}\frac{\text{Re}^{\frac{1}{2}}}{M^2}`$, where the Mach number is $`Mu_{}/c`$.

It is convenient to treat the occupation probabilities as a two-component field

$$|p=\left(\begin{array}{c}p_+\\ p_{}\end{array}\right).$$
(6)
We expand $`|p`$ about equilibrium its value denoted $`|d`$ so that $`|p=|d+|\delta p+𝒪(\epsilon ^2)`$, where $`\epsilon \frac{1}{N_x}`$ is the Knudsen number. The equilibrium condition $`\mathrm{\Omega }|_{p=d}=0`$ leads to a tractable polynomial equation for $`d_\pm `$, whence the linearized finite-difference equation

$$|p(x\pm \delta x,t+\delta t)|p(x,t)=J|\delta p(x,t),$$
(7)
where the Jacobian of the collision term is

$$J\left(\begin{array}{cc}\frac{\mathrm{\Omega }}{p_+}& \frac{\mathrm{\Omega }}{p_{}}\\ \frac{\mathrm{\Omega }}{p_+}& \frac{\mathrm{\Omega }}{p_{}}\end{array}\right)|_{p=d}=\left(\begin{array}{cc}J_+& J_{}\\ J_+& J_{}\end{array}\right).$$
(8)
The left and right eigenvectors of $`J`$ are

$`\xi _1|=\left(\begin{array}{cc}1& 1\end{array}\right)`$ $`|\xi _1=`$ $`{\displaystyle \frac{1}{J_{}J_+}}\left(\begin{array}{c}J_{}\\ J_+\end{array}\right)`$ (9)
$`\xi _2|={\displaystyle \frac{1}{J_+J_{}}}\left(\begin{array}{cc}J_+& J_{}\end{array}\right)`$ $`|\xi _2=`$ $`\left(\begin{array}{c}1\\ 1\end{array}\right)`$ (10)
with associated eigenvalues $`\lambda _1=0`$ and $`\lambda _2=J_+J_{}`$. $`\xi _i|\xi _j=\delta _{ij}`$. $`J`$ may be rewritten as

$$J=\lambda _2|\xi _2\xi _2|=\left(\begin{array}{cc}J_+& J_{}\\ J_+& J_{}\end{array}\right)$$
(11)
$`J`$ is singular. Nevertheless, its generalized inverse is

$$J_{\text{gen}}^1=\frac{1}{\lambda _2}|\xi _2\xi _2|=\frac{1}{J_+J_{}}\left(\begin{array}{cc}J_+& J_{}\\ J_+& J_{}\end{array}\right).$$
(12)
Now we invoke the continuum limit where $`\delta x0`$ and $`\delta t0`$ so $`|p`$ is a continuous and differentiable two-component field. We obtain a first order equation by Taylor expanding (7) in $`x`$ and $`t`$ and keeping terms only terms first order in $`\epsilon `$:

$$\sigma _z\delta x_x|d=J|\delta p+𝒪(\epsilon ^2),$$
(13)
where $`\sigma _z=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)`$. Then using (12) we have

$$|\delta p=\frac{1}{J_+J_{}}\sigma _z\delta x_x|d+𝒪(\epsilon ^2).$$
(14)
Taking the difference of the respective components gives

$$\delta p_+\delta p_{}=\frac{1}{J_+J_{}}\delta x_x\rho +𝒪(\epsilon ^2).$$
(15)
Similarly from (7), we obtain the second order equation:

$`\delta t_t|d`$ $`+`$ $`\sigma _z\delta x_x\left(|d+|\delta p\right)`$ (16)
$`+`$ $`{\displaystyle \frac{\delta x^2}{2}}_x^2|d+𝒪(\epsilon ^3)=\left(\begin{array}{c}\mathrm{\Omega }\\ \mathrm{\Omega }\end{array}\right).`$
We now take the sum of the respective components:

$`\delta t_t\rho `$ $`+`$ $`\delta x_x\left(d_+d_{}+\delta p_+\delta p_{}\right)`$ (17)
$`+`$ $`{\displaystyle \frac{\delta x^2}{2}}_x^2\rho +𝒪(\epsilon ^3)=0.`$
Inserting (15) into the above equation gives the general effective field theory for any one-dimensional two-particle-per-site lattice gas conserving particle number

$`_t\rho `$ $`+`$ $`c_x\left(d_+d_{}\right)+{\displaystyle \frac{\delta x^2}{\delta t}}{\displaystyle \frac{_x(J_+J_{})}{(J_+J_{})^2}}_x\rho `$ (18)
$`+`$ $`{\displaystyle \frac{\delta x^2}{2\delta t}}\left({\displaystyle \frac{2}{J_+J_{}}}+1\right)_x^2\rho +𝒪(\epsilon ^3)=0.`$
The Brieger-Bonomi model of the Burgers equation has one bit per site with a 3-bit stencil where the center bit is updated by a stochastic Masters equation brieger-jsp92 . This nonlocal classical model is equivalent to a 3-bit per site local lattice-gas model. We briefly discuss the classical Boghosian-Levermore lattice-gas model boghosian-87 ; boghosian-89  before proceeding to our quantum model. The classical collision term is

$$\mathrm{\Omega }=\frac{1}{2}(p_{}p_+)+\frac{\alpha }{2}\left(p_++p_{}2p_+p_{}\right),$$
(19)
where $`\alpha `$ is occupation probability of an additional random bit used to bias the collision 111Actually, in the Boghosian-Levermore model, 59 additional bits per site were used to generate a single random bit in a strictly reversible fashion.. Taking the equilibrium to be of the form $`d_\pm =\frac{\rho }{2}\pm A`$, then $`\mathrm{\Omega }|_{p=d}=0`$ leads to a quadratic equation in $`A`$. We take the negative root $`A=\alpha \frac{\rho }{2}\left(1\frac{\rho }{2}\right)`$, and hence

$$d_+d_{}=\alpha \rho (1\frac{\rho }{2}).$$
(20)
With $`J_\pm =\frac{\mathrm{\Omega }}{p_\pm }|_{p=d}`$, we find

$$J_+J_{}=1+2\alpha ^2\rho (2\rho )1+𝒪(\alpha ^2).$$
(21)
Then substituting the two results (20) and (21) into (18) yields the desired effective field theory

$$_t\rho +c\alpha (1\rho )_x\rho \frac{\delta x^2}{2\delta t}_x^2\rho +𝒪(\epsilon ^3,\epsilon \alpha ^2)=0,$$
(22)
which is the nonlinear Burgers equation for $`u=c(\rho 1)`$ with variable sound speed $`c_s=c\alpha `$ and fixed kinematic viscosity $`\nu =\frac{\delta x^2}{2\delta t}`$.

In the special case when $`\alpha =0`$, (22) reduces to the diffusion equation. However, there exists an abnormal case when $`\alpha =1`$ where the random bit is fixed to the value of 1. The 3-bit model reduces to a 2-bit one that does not model the Burgers equation. Its collision term is $`\mathrm{\Omega }=(1p_+)p_{}`$ and the equilibrium occupations are $`d_+=\rho `$ and $`d_{}=0`$. The components of the Jacobian matrix are $`J_+=0`$ and $`J_{}=1+\rho `$. Hence $`d_+d_{}=\rho `$ and $`J_+J_{}=1\rho `$, so (18) becomes

$$_t\rho +c_x\rho +\frac{\delta x^2}{\delta t}\frac{\left(_x\rho \right)^2}{(1\rho )^2}+\frac{\delta x^2}{\delta t}\left(\frac{\rho +1}{\rho 1}\right)_x^2\rho +𝒪(\epsilon ^3)=0.$$
(23)
See figure 1 for numerical solutions of (22) and (23).

Now we consider the quantum algorithm for the Burgers equation. Initially we encode the qubits with their respective occupation probabilities

$$|q_\pm =\sqrt{p_\pm }|1+\sqrt{1p_\pm }|0,$$
(24)
ignoring the possibility of an internal phase angle. The initial ket $`|\psi =|q_+|q_{}=\sqrt{p_+p_{}}|11+\sqrt{p_+(1p_{})}|10+\sqrt{(1p_+)p_{}}|01+\sqrt{(1p_+)(1p_{})}|00`$ is transformed by application of a unitary matrix:

$$|\psi ^{}=\widehat{U}|\psi .$$
(25)
The identity matrix and $`\widehat{n}=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)`$ denote the single qubit number operator, the multi-qubit number operators $`\widehat{n}_1=\widehat{n}\mathrm{𝟏}`$ and $`\widehat{n}_2=\mathrm{𝟏}\widehat{n}`$ are used to determine the new probabilities of the respective updated qubits

$`p_+^{}`$ $``$ $`\psi ^{}|\widehat{n}_1|\psi ^{}=p_++\mathrm{\Omega }_{\text{QLG}}(p_+,p_{})`$ (26)
$`p_{}^{}`$ $``$ $`\psi ^{}|\widehat{n}_2|\psi ^{}=p_{}\mathrm{\Omega }_{\text{QLG}}(p_+,p_{}).`$
(26) implicitly determines the functional form of the mesoscopic collision term $`\mathrm{\Omega }_{\text{QLG}}(p_+,p_{})`$ associated with the microscopic operator $`\widehat{U}`$. We use a conservative collision operator as our 2-qubit quantum gate entangling the qubits using only the microscopic states $`|01`$ and $`|10`$

$$\widehat{U}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& e^{i\xi }\mathrm{cos}\theta & e^{i\zeta }\mathrm{sin}\theta & 0\\ 0& e^{i\zeta }\mathrm{sin}\theta & e^{i\xi }\mathrm{cos}\theta & 0\\ 0& 0& 0& 1\end{array}\right).$$
(27)
Inserting (27) into (25), and then substituting the resulting $`|\psi ^{}`$ into (26), we find the collision term $`\mathrm{\Omega }_{\text{QLG}}`$ is yepez-jstatphy01

$`\mathrm{\Omega }_{\text{QLG}}`$ $``$ $`\mathrm{sin}^2\theta (p_{}p_+)`$
$`+`$ $`\mathrm{sin}2\theta \mathrm{cos}(\zeta \xi )\sqrt{p_+(1p_+)p_{}(1p_{})}.`$
The equilibrium condition $`\mathrm{\Omega }_{\text{QLG}}|_{p=d}=0`$ becomes:

$$\frac{d_+}{1d_+}\frac{d_{}}{1d_{}}=2\mathrm{cot}\theta \mathrm{cos}(\zeta \xi )\sqrt{\frac{d_+}{1d_+}\frac{d_{}}{1d_{}}},$$
(29)
which is a legitimate statement of detailed-balance of collisions at the mesoscopic scale since the quantum model’s evolution operator is unitary. We take the equilibrium occupation probabilities to have the following form:

$$d_+=\frac{1}{\gamma z+1}\text{and}d_{}=\frac{1}{\frac{z}{\gamma }+1}.$$
(30)
Substituting (30) into (29) gives a quadratic equation in $`\gamma `$ that has the solution $`\gamma =\sqrt{\alpha ^2+1}+\alpha `$ or $`\frac{1}{\gamma }=\sqrt{\alpha ^2+1}\alpha `$, where $`\alpha \mathrm{cot}\theta \mathrm{cos}(\zeta \xi )`$. Next, substituting (30) into the total number density, $`\rho =d_++d_{}`$, we obtain a quadratic equation in $`z`$

$$\rho z^2+\left(\gamma +\frac{1}{\gamma }\right)(\rho 1)z+\rho 2=0.$$
(31)
Substituting the positive root solution of (31) into (30), we find after much algebraic manipulation

$$d_+=\frac{1+\gamma ^2+(1\gamma ^2)\rho \gamma \sqrt{\left(\frac{1}{\gamma }+\gamma \right)^2(\rho 1)^2+4(\rho 2)\rho }}{2(1\gamma ^2)}.$$
(32)
Then substituting $`\gamma =\sqrt{\alpha ^2+1}+\alpha `$ into (32) gives the result

$$d_\pm =\frac{\rho }{2}\frac{1}{2\alpha }\left(\sqrt{1+\alpha ^2}\sqrt{1+\alpha ^2(\rho 1)^2}\right).$$
(33)
This implies that

$$d_+d_{}=\frac{1}{\alpha }\left(\sqrt{1+\alpha ^2}\sqrt{1+\alpha ^2(\rho 1)^2}\right).$$
(34)
Again, we compute the components of $`J`$:

$$J_\pm =\frac{\mathrm{\Omega }_{\text{QLG}}}{p_\pm }=\mathrm{sin}^2\theta \left(1\alpha \frac{(2d_\pm 1)d_{}(1d_{})}{\sqrt{d_+(1d_+)d_{}(1d_{})}}\right).$$
(35)
And this implies

$$J_+J_{}=2\mathrm{sin}^2\theta (1+\alpha ^2f),$$
(36)
where the factor $`f=f(\alpha ,\rho )`$ is too complicated an expression to write out here but has the important property that $`f(\alpha ,\rho )=1+𝒪(\alpha )`$. Finally, substituting the two results (34) and (36) into (18) gives the effective field theory

$$_t\rho +c\mathrm{cot}\theta \mathrm{cos}(\zeta \xi )(1\rho )_x\rho =\mathrm{cot}^2\theta \frac{\delta x^2}{\delta t}_x^2\rho +𝒪(\epsilon ^3,\epsilon \alpha ^2),$$
(37)
which is the nonlinear Burgers equation for $`u=c(\rho 1)`$ with independently tunable sound speed $`c_s=c\mathrm{cot}\theta \mathrm{cos}(\zeta \xi )`$ and kinematic viscosity $`\nu =\mathrm{cot}^2\theta \frac{\delta x^2}{\delta t}`$ by appropriately choosing the Euler angles in (27). In our quantum case, the trigonometric term $`\mathrm{cot}\theta \mathrm{cos}(\zeta \xi )`$ plays the role of the expectation value $`\alpha `$ of the additional random bit required in the classical model. Figures 1 and 2 show the time evolution of the 2-qubit quantum algorithm versus the 3-bit classical Boghosian-Levermore algorithm both carried out on a $`N_x=256`$ lattice. The vertical axis is the particle number density $`\rho =p_++p_{}`$ plotted in the range of $`0.5\rho 1.5`$. The time step is in the upper left corner of each snapshot. The viscosity of the quantum model is close to zero.

The quantum algorithm is unconditionaly stable, obeys detailed-balance, requires less memory than its classical counterpart, and can achieve arbitrarily high Reynolds numbers. Having a variable transport coefficient that can be made small, it is consistent with the inviscid Burgers equations when the Euler angle $`\theta \frac{\pi }{2}`$. It is possible to generalize this type of quantum algorithm to three-dimensions to efficiently handle the important application of computational fluid dynamics.

I would like to acknowledge Professor Boghosian for helpful discussions about his classical lattice-gas model and Owen Cote for helpful discussions about turbulence.


text:
# Power-law temporal auto-correlations in day-long records of human physical activity and their alteration with disease

## Abstract

We investigate long-duration time series of human physical activity under three different conditions: healthy individuals in (i) a constant routine protocol and (ii) in regular daily routine, and (iii) individuals diagnosed with multiple chemical sensitivities. We find that in all cases human physical activity displays power law decaying temporal auto-correlations. Moreover, we find that under regular daily routine, time correlations of physical activity are significantly different during diurnal and nocturnal periods but that no difference exists under constant routine conditions. Finally, we find significantly different auto-correlations for diurnal records of patients with multiple chemical sensitivities.

Healthy free-running physiologic systems have complex self-regulating mechanisms which process inputs with a broad range of characteristics Shlesinger , and may generate signals that have scale-invariant dynamics Musha ; Peng95 . Physical activity is a physiologic signal of great interest due to its (i) impact on other rhythms such as heart rate, and (ii) significance for certain psychiatric or psychosomatic conditions Teicher95 ; psy . Indeed, a number of researchers has used actigraph recordings to classify types of motion—such as sitting, moving or standing children —with the goal of developing ambulatory methods for the diagnosing of attention-deficit hyperactivity disorder in children adhd . In such applications, the advantage of using actigraph recordings is that it may replace more limiting, cumbersome, and intrusive techniques such as videotaping monitoring . Surprisingly, up to now there has been no attempt to quantify the time correlations of the fluctuations in the levels of physical activity in “free-running” ambulatory settings or their alteration under special conditions, such as constrained activity or psychiatric disease Teicher .

In this Letter, we study long-duration time series of human physical activity under two different conditions: a constant routine (CR) protocol where physical activity and postural changes are kept to a minimum Constant\_rout , and (ii) regular daily routine (RDR). We find that physical activity displays power law decaying temporal correlations for these conditions. Moreover, we find that during diurnal periods the time correlations of physical activity are significantly different for regular daily routine and CR protocols. In contrast, we find no differences between diurnal and nocturnal periods for the CR protocol while we find significant differences between diurnal and nocturnal periods for the regular daily routine. Further, we find similarities between the CR patterns of human physical activity and those of nocturnal periods for regular daily routine note-cr .

Our finding of long-range time-correlations in activity of healthy individuals under regular daily activity and the change in those correlation properties under controlled conditions, suggests that possibility that the fractal index characterizing the long-range auto-correlations of the time series of physical activity may be an additional probe Teicher95  of the abnormal patterns of physical activity of patients with psychiatric disorders such as depression, bipolar disorder, or schizophrenia DSM , and with psychosomatic illnesses associated with exaggerated fatigue. We test this possibility with data collected for individuals diagnosed with multiple chemical sensitivities (MCS). This is a disease mcs\_criteria  characterized by a positive chemical exposure history and multi-system symptoms, such as fatigue, headache, sleep disturbance, and myalgia, which are elicited after exposure to various chemical compounds buchwald94 . The clinical picture of MCS is shown buchwald94  to largely overlap with that of other psychosomatic illnesses with severe fatigue symptoms holmes88 ; fukuda94 . Interestingly, we find a significant difference in the degree of diurnal auto-correlations in physical activity for MCS patients and healthy controls leading a regular daily routine.

We analyze datasets from seven healthy, non-smoking male subjects Tokyo  (ages: 21–30 yr). For six of the subjects, we obtained two datasets per subject, the first under CR conditions, and the second under regular daily routine conditions Constant\_rout ; na00 ; Amaral01 . We recorded na00  body acceleration using a portable, waist-worn, long-term ambulatory monitoring device (Amx720, Nihon-Koden Wellness Corporation, Japan). The monitoring device has two shock sensors that measure body acceleration at the waist in the vertical and horizontal directions with a sensitivity of 0.08$`g`$. The acceleration was recorded after being full-wave rectified and integrated over 8 s intervals na00 . In the present study, we consider only body acceleration in the vertical direction.

Figure 1 displays vertical acceleration at the waist for subject #2 for both CR protocol and regular daily routine, Vertical acceleration at the waist is a good proxy for physical activity because it correlates to entire body motion. The figure suggests a decrease in the number of episodes of sustained physical activity during CR but the quantitative nature of such time organization is not clear from the graph. Thus, we apply the detrended fluctuations analysis (DFA) method to quantify long-range time-correlations in the physical activity time series Peng95 ; Ashkenazy01 . The DFA method is defined as follows: One first integrates the physical activity time series. One then divides the time series into “boxes” of length $`n`$ and perform, in each box, a least-squares polynomial fit of order $`k`$ fit to the integrated signal. Next, one calculates in each box the root-mean-square deviations of the integrated signal from the linear fit. The basic idea is that the polynomial fit represents the local trend in each box while the magnitude of the fluctuations around the trend signals the degree of the time correlations in the signal.

This procedure is repeated for different box sizes (time scales) $`n`$. For fractal signals one finds a power-law relation between the average magnitude of the fluctuations $`F(n)`$ and the number of points $`n`$

$$F(n)n^\alpha ,$$
(1)
where the scaling exponent $`\alpha `$ quantifies the degree of the correlations. Uncorrelated time series yield $`\alpha =0.5`$, long-range anti-correlations result in $`\alpha <0.5`$, and long-range positive correlations results in $`\alpha >0.5`$.

In Fig. 2, we show the results of the DFA analysis for subject #2 in our database. We consider separately diurnal and nocturnal periods because of recent reports indicating altered physiologic control during sleep Ivanov ; Bunde\_sleep . It is visually apparent that the data for diurnal and nocturnal periods and for both the CR protocol and the regular daily routine follow straight lines in the double-logarithmic plots. This fact confirms our expectation that Eq. (1) holds. Surprisingly, we find different exponents during diurnal periods for CR and regular daily routine while we find similar exponent values for the nocturnal periods. This finding suggest a difference in the type of voluntary physical activity note  during diurnal regular daily routine.

We test this finding for all other subjects in the database (Fig. 3 and Table 1). Our analysis reveal a value of $`\alpha `$ systematically larger for regular daily routine than for CR during the daytime. We test this possibility with two statistical test nr : the Student t-test, which tests differences in mean values under the assumption of Gaussianity of the random variables, and the Kolmogorov-Smirnov (KS) test, which is a non-parametric test for distributions. We find that the difference in the average values of $`\alpha `$ for the regular daily routine and CR during diurnal periods is significant at the $`p<0.001`$ level for both Student’s t-test and the KS test. We also find that the difference in the average values of $`\alpha `$ for the regular daily routine and CR during nocturnal periods is not significant note\_magnitude .

To investigate the medical significance of our results for the characterization of psychiatric and psychosomatic diseases, we next analyze data collected over a seven day period for four individuals (ages: 34–47 yr) diagnosed with MCS mcs-subjects . We find that the auto-correlations of the diurnal physical activity time series are characterized by a DFA exponent $`\alpha =0.93\pm 0.001`$ (Table 1 and Fig. 3) which we find to be significantly different from the results for healthy subjects at the $`p<10^5`$ level for both the Student’s t-test and the KS test.

A key finding of this study is that voluntary physical activity has a power-law decaying time correlations whose characteristic index $`\alpha `$ changes when patterns of physical activity are controlled or disturbed. For diurnal periods, we measure a smaller value of $`\alpha `$ for CR than for regular daily activity. The picture emerging from our results raises the possibility that conditions such as hyperactivity, depression and other psychiatric and psychosomatic disorders will modify the “pull” of daily routines on physical activity and lead to different values of $`\alpha `$. Hence, it is plausible that our finding of fractal correlations in the activity levels of healthy mature individuals provides a new window into the study of psychiatric and maturation disorders which is not hindered by complex (and more falsifiable) procedures such as classifying motion types. Indeed, we test this possibility for MCS patients and confirm that correlations in activity are altered by disease. Further, we have also found Kyoko  changes in correlation for patients with a severe fatigue illness called chronic fatigue syndrome holmes88 ; fukuda94 .

Our findings are also of importance because of their implications on the question of the influence of fractal stimuli, such as physical activity, on the properties of other rhythms such as heart rate. Specifically, it has been recently shown that time series of healthy human interbeat intervals belong to a special class of complex signals that display multifractal properties Ivanov99 ; multi . An explanation is that the neuroautonomic control mechanisms —in the presence of even weak external noise— endogenously generate multifractal dynamics.

A recent study showed that the heart rate variability of healthy adults displays similar multifractal properties for CR and regular daily routine conditions Amaral01 . In contrast, sympathetic or parasympathetic blockades were found to lead to a significant loss of complexity of HRV Amaral01 . It would thus appear that—in contrast to neuroautonomic control—voluntary physical activity may not be an important factor in the generation of the fractal Meesmann  or multifractal properties of heartbeat dynamics Amaral01 . However, an important assumption for validating the above conclusion is that voluntary physical activity during the CR protocol Constant\_rout  is significantly different from that during regular daily routine. Indeed, if during CR conditions there is only a decrease of the amplitude in the physical activity, then one would need to show that such a decrease during the CR protocol lowers the amplitude of the physical activity below a threshold value in order to justify the endogenous character of the multifractal properties of heartbeat dynamics. Our finding that the fractal time organization of physical activity during CR conditions is significantly different from the time organization during diurnal regular daily routine strengthens the conclusion that voluntary physical activity is not responsible for the multifractal character of heartbeat dynamics.

We thank P. Ch. Ivanov, A. L. Goldberger, and H. E. Stanley for stimulating discussions. L.A.N.A. thanks the hospitality of the Universidade Federal do Rio Grande do Norte and support from NIH/NCRR (P41 RR13622) Physionet . D.J.B.S., L.R.S. and L.S.L. acknowledge support from the Brazilian Agencies CNPq and FINEP. Y.Y. acknowledges a MECSST Grant-in-Aid for Scientific Research and support from Japan Science and Technology Corporation.


text:
# Born–Infeld theory of electroweak and gravitational fields: Possible correction to Newton and Coulomb laws

In this note one suggests a possibility of direct observation of the $`\theta `$-parameter, introduced in the Born–Infeld theory of electroweak and gravitational fields, developed in quant-ph/0202024. Namely, one may treat $`\theta `$ as a universal constant, responsible for correction to the Coulomb and Newton laws, allowing direct interaction between electrical charges and masses.

1. In the note  a Born–Infeld (BI)  theory of electroweak and gravitational fields with the action invariant under $`\theta `$-rotations of the Lorentz group, imposed on pairs ($`\gamma _a,\pi _a`$), was developed. The following equations (slightly modified*** Here definition for the energy is $`E=c\overline{\psi }p^0\psi 𝑑\mathrm{\Sigma }_0`$, and for the electrical charge is $`Q=q\overline{\psi }\gamma ^0\psi 𝑑\mathrm{\Sigma }_0`$, where $`\psi `$ is a spinorial wave function of a particle, and $`q`$ is electrical charge of electron (absolute value). from (6.5) and (6.6) in ) were obtained:

$`E^{}`$ $`=`$ $`\mathrm{cosh}\theta E+{\displaystyle \frac{\mathrm{}c}{\mathrm{}q}}\mathrm{sinh}\theta Q;`$ (1)
$`Q^{}`$ $`=`$ $`\mathrm{cosh}\theta Q+{\displaystyle \frac{\mathrm{}q}{\mathrm{}c}}\mathrm{sinh}\theta E.`$ (2)
Here

$`\mathrm{}^2`$ $`=`$ $`{\displaystyle \frac{10}{3\alpha }}(1+2\mathrm{sin}^2\theta _W)L_p^2,`$ (3)
where $`L_p`$ is the Planck’s length, $`\alpha `$ is the fine structure constant for electron, and $`\theta _W`$ is the Weinberg mixing angle . Equations (1) and (2) mean that if ($`E`$, $`Q`$) are possible solutions of field equations following from the theory for the energy and the electrical charge of, say, an electron, then ($`E^{}`$, $`Q^{}`$) are also possible solutions for the energy and the electrical charge of the same electron.

2. Consider, now, nonrelativistic case, so that $`Emc^2`$, and $`E^{}m^{}c^2`$. Suppose, that $`\theta 1`$, so that one may neglect terms of order $`\theta ^2`$. Substituting $`e`$ and $`e^{}`$ instead of $`Q`$ and $`Q^{}`$ respectively, one obtains from (1) and (2),

$`m^{}`$ $`=`$ $`m+\kappa ^{\frac{1}{2}}\theta e;`$ (4)
$`e^{}`$ $`=`$ $`e+\kappa ^{\frac{1}{2}}\theta m.`$ (5)
Here

$$\kappa =\frac{q^2\mathrm{}^2c^2}{\mathrm{}^2}.$$
(6)
One may interpret $`\mathrm{\Delta }m=\kappa ^{\frac{1}{2}}\theta e`$ and $`\mathrm{\Delta }e=\kappa ^{\frac{1}{2}}\theta m`$ as contributions of electrical charge, $`e`$, and mass, $`m`$, to the observable mass, $`m^{}`$, and electrical charge, $`e^{}`$, respectively.

3. Keeping only leading terms of the action’s expansion in $`\mathrm{}`$ , one would obtain standard equations for Coulomb’s electrostatics and Newton’s gravity in nonrelativistic case. Consider two electrical charges, $`(e,m)`$ and $`(Q,M)`$. Then, observable quantities are $`(e^{},m^{})`$ and $`(Q^{},M^{})`$. Here $`Q^{}`$ and $`M^{}`$ are expressed through $`Q`$ and $`M`$ similarly to (5) and (4). Then, one obtains for the Newton and Coulomb potential energy of two particles, separated at distance $`r`$,

$$V(r)=\frac{e^{}Q^{}km^{}M^{}}{r}.$$
(7)
Here $`k`$ is the Newton constant of gravitational interaction. Substituting ‘unprimed’ variables, one obtains for the potential energy,

$$V(r)=\frac{eQkmM+\sigma (eM+mQ)}{r},$$
(8)
where

$$\sigma =\theta \left(\sqrt{\kappa }\frac{k}{\sqrt{\kappa }}\right).$$
(9)
The new effect (proportional to $`\sigma `$) may be responsible for interaction between cosmic charged particles (protons, electrons, etc.) and the Earth’s ‘charge’, $`Q_E=\sigma M_E`$. Here $`M_E`$ is the Earth’s mass. Assuming that $`\sigma <0`$, one obtains that Earth’s mass contributes negative electric charge, which explains, why primary cosmic rays consist mainly from positively charged particles . The same assumption explains the fairweather electric field (100 Volts per metre at surface of the Earth) . Due to the same effect rotation of the Earth leads to electric currents, which may be responsible for Earth’s magnetic field.