text:
# The Lee-Yang theory of equilibrium and nonequilibrium phase transitions

## I Introduction

In this work we seek a mathematical understanding of phase transitions in the steady state of stochastic many-body systems. Systems at equilibrium with their environment provide examples of such steady states, and the mechanisms underlying equilibrium phase transitions are long known and understood. Experimentally, one can distinguish between two types of transition: the first-order transition at which there is phase coexistence, e.g. between a high-density solid and a low-density fluid, and the continuous transition at which fluctuations and correlations grow to such an extent as to be macroscopically observable.

From a thermodynamical perspective, one can understand first-order transitions by associating with each phase a free energy. For a given set of external parameters, the phase ‘chosen’ by the system is that with the lowest free energy, and so a phase transition occurs when the free energies of two (or more) phases are equal. The sudden changes in macroscopically measurable quantities that take place at first-order transitions are described mathematically as discontinuities in the first derivative of the free energy. Discontinuities in higher derivatives relate to continuous (higher-order) phase transitions.

The tools of equilibrium statistical physics allow one—in principle at least—to express the free energy solely in terms of microscopic interactions. More specifically, the free energy is given by the logarithm of a partition function, a quantity that normalises the steady-state probability distribution of microscopic configurations. Initially it was not universally accepted that this approach could faithfully describe phase transitions, in particular the first-order solid-fluid transition Domb96 . In order to show that the statistical mechanical approach can reproduce the correct discontinuities in the free energy at a first-order transition, Lee and Yang YL52  introduced a description of phase transitions concerning zeros of the partition function when generalised to the complex plane of an intensive thermodynamic quantity. Initially, the theory was couched in terms of zeros in the complex fugacity plane which is appropriate for fluids in the grand canonical (fluctuating particle number) ensemble. However, by mapping a lattice gas onto the Ising model LY52 , the theory was found to hold equally well for a magnet in a fixed magnetic field. Later Fisher65 ; GR69 ; GL69  it became clear that the distribution of zeros in the complex temperature plane can reveal information about phase transitions in the canonical ensemble.

In section II we present a brief, self-contained discussion of the Lee-Yang theory of equilibrium phase transitions which relates nonanalytic behaviour in the free energy to zeros of the partition function. In the remainder of this article, we move our attention to phase transitions in nonequilibrium steady states, that is, those that carry currents of mass, energy or some other quantity. We give therefore a very brief introduction to the established practice of modelling nonequilibrium physics through stochastic dynamics in section III. In particular, we will identify a candidate quantity for use in place of an equilibrium partition function when performing an Lee-Yang analysis of phase transitions. Finally, in section IV we review recent work in which the Lee-Yang theory has been successfully applied to specific nonequilibrium phase transitions.

## II Overview of Lee-Yang theory of equilibrium phase transitions

In this section we recount the relationship of zeros of a partition function to phase transitions in the system it describes. We present here a brief but self-contained tour of the key points. For more detail (and mathematical rigour) the reader is referred to textbooks, such as Reichl80  or, for a more advanced presentation, Ruelle69 , as well as early treatments of the subject YL52 ; LY52 ; Fisher65 ; GR69 ; GL69 .

For concreteness we shall describe the theory in terms of a simple model system of $`N`$ spins in thermal contact with a heat reservoir. In this model system the energy of the system can take values $`E=nϵ`$ with $`n=0,1,2,\mathrm{},M`$ and the number of microstates corresponding to the $`n^{\mathrm{th}}`$ energy level is $`g(n)`$. The canonical partition function is given by

$$Z_N(\beta )=\underset{n=0}{\overset{M}{}}g(n)\mathrm{exp}(\beta nϵ)$$
(1)
where $`\beta `$ is the inverse temperature. To investigate the zeros of this partition function, it is useful to make the change of variable $`z=\mathrm{e}^{\beta ϵ}`$. Then, the partition function is explicitly a polynomial in $`z`$ of degree $`M`$ and can be factorised as

$$Z_N(z)=\underset{n=0}{\overset{M}{}}g(n)z^n=\kappa \underset{n=1}{\overset{M}{}}\left(1\frac{z}{z_n}\right).$$
(2)
Clearly the quantities $`z_n`$ are the $`M`$ zeros of the partition function $`Z_N(z)`$; meanwhile, $`\kappa `$ is a constant which we can safely ignore in the following.

Since all the $`g(n)`$ are positive, no zero of $`Z_N(z)`$ can be real and positive: that is, the zeros $`z_n`$ will generally lie in the complex plane away from the physical values of $`z`$ which lie on the positive real axis. We now define for all complex $`z`$ except the points $`z=z_n`$ the (complex) free energy

$$h_N(z)\frac{\mathrm{ln}Z_N(z)}{N}.$$
(3)
Using (2) we rewrite this as

$$h_N(z)=\frac{1}{N}\underset{n=1}{\overset{M}{}}\mathrm{ln}\left(1\frac{z}{z_n}\right)$$
(4)
and note that a Taylor series expansion of $`h_N(z)`$ around a point $`zz_n`$ has a finite radius of convergence given by the distance to the nearest zero from $`z`$. This then implies that that $`h_N(z)`$ can be differentiated infinitely many times over any region of the complex plane that is devoid of partition function zeros. Since we identify a phase transition through a discontinuity in a derivative of the free energy, we see that such a transition can only occur at a point $`z_0`$ in the complex plane if there is at least one zero of the partition function $`Z_N(z)`$ within any arbitrarily small region around the point $`z_0`$.

Clearly this scenario is impossible if the number of zeros $`M`$ is finite, except at the isolated points $`z_n`$ where the free energy exhibits a logarithmic singularity. Since such a point cannot lie on the positive real $`z`$ axis, there is no scope for a phase transition in a finite spin system, such as the simple example (1). On the other hand, if the partition function zeros accumulate towards a point $`z_0`$ on the real axis as we increase the number of spins $`N`$ to infinity there is the possibility of a phase transition.

In order to deal with the thermodynamic limit (see Ruelle69  for rigorous considerations) we shall assume that the limit

$$h(z)=\underset{N\mathrm{}}{lim}\frac{\mathrm{ln}Z_N(z)}{N}$$
(5)
exists and we may write

$$h(z)=dz^{}\rho (z^{})\mathrm{ln}\left(1\frac{z}{z^{}}\right)$$
(6)
where $`\rho (z)`$ is the local density of zeros in the complex-$`z`$ plane in the thermodynamic limit.

Since the imaginary part of this complex free energy $`h(z)`$ is multi-valued it will at times be more convenient to work with the potential $`\phi (z)`$ defined as

$$\phi (z)\mathrm{Re}h(z)=dz^{}\rho (z^{})\mathrm{ln}\left|1\frac{z}{z^{}}\right|.$$
(7)
An expression for the density $`\rho (z)`$ in terms of the potential $`\phi (z)`$ is easily obtained once one recognises that $`\mathrm{ln}|z|`$ is the Green function for the two-dimensional Laplacian. Specifically, if $`z=x+iy`$, we have

$$^2\mathrm{ln}|z|\left(\frac{^2}{x^2}+\frac{^2}{y^2}\right)\mathrm{ln}|x+iy|=2\pi \delta (x)\delta (y).$$
(8)
Using this expression, we can take the Laplacian of both sides of (7) to find

$$\rho (z)=\frac{1}{2\pi }^2\phi (z).$$
(9)
Such an equation is familiar from electrostatics and so $`\phi (z)`$ is analogous to an electrostatic potential.

In analogy to electrostatics, as long as we can integrate $`\rho (x)`$ over bounded regions containing any singularities of $`\rho `$, the potential will be a continuous function kellogg . The significance of this statement is that we can derive a rule for locating phase boundaries given a partition function. Let us suppose that around points $`z_1`$ and $`z_2`$ in the complex plane, one has (analytic) expressions for the potential $`\phi _1(z)`$ and $`\phi _2(z)`$ such that $`\phi _1(z)\phi _2(z)`$ (i.e. not the same function over the entirety of the complex plane). In order for the potential to be continuous at all points on the complex plane, we must have a phase boundary at those values of $`z`$ for which the condition

$$\phi _1(z)=\phi _2(z)$$
(10)
holds. This is basically the definition of a transition mentioned above in the introduction. Since $`\phi _1`$ and $`\phi _2`$ are different functions, some derivatives of $`\phi (z)`$ will not exist at these values of $`z`$ and we expect the density of zeros at these points to be non-zero. (It is also possible for zeros to be present at other points in the complex plane, but we do not need to consider this possibility here.)

Typically, a solution of (10) describes a curve $`C`$ that intersects the positive real $`z`$ axis at a point $`z_0`$. Having already established that we require zeros to accumulate at the point $`z_0`$ in the thermodynamic limit for a physical phase transition to take place, we are interested in the line density of zeros per unit length of this curve. We introduce the arclength $`s`$ which measures the distance along $`C`$ from the transition point $`z_0`$.

To obtain an expression for the line density of zeros $`\mu (s)`$ along the curve $`C`$, consider a short section of $`C`$ of length $`\delta s`$ enclosed by an area $`A`$ that has two sides parallel to $`C`$ and the other two sides perpendicular to $`C`$—see figure 1. Integrating the density $`\rho (z)`$ over this area we have

$$_Adz\rho (z)=\mu (s)\delta s.$$
(11)
Meanwhile, we have from (9) and an application of the divergence theorem that

$$_Adz\rho (z)=\frac{1}{2\pi }_Adz[\phi (z)]=\frac{\delta s}{2\pi }[\phi _2(z)\phi _1(z)]\widehat{n}$$
(12)
in which the functions $`\phi _1(z)`$ and $`\phi _2(z)`$ relate to the limiting values of $`\phi (z)`$ as the point $`z`$ on the curve $`C`$ is approached from either side, and the vector $`\widehat{n}`$ is the unit vector normal to $`C`$ at that point.

Recall that $`\phi (z)`$ is the real part of the complex free energy $`h(z)`$ and therefore, away from any zeros of the partition function, satisfies the Cauchy-Riemann equations

$$_x\phi (z)=_y\psi (z)\text{and}_y\phi (z)=_x\psi (z)$$
(13)
in which $`\psi (z)`$ is the imaginary part of $`h(z)`$ and $`z=x+iy`$. Then we have that

$$\phi (z)\widehat{n}=\psi (z)\widehat{t}$$
(14)
where $`\widehat{t}`$ is the unit vector tangent to the curve $`C`$ at the point $`z`$. Thus we recognise the scalar product in (12) as the directional derivative of the imaginary part of the free energy $`h(z)`$ along the curve $`C`$. Putting this together with (11) we find that

$$\mu (s)=\frac{1}{2\pi }\frac{\mathrm{d}}{\mathrm{d}s}[\psi _2(z)\psi _1(z)]$$
(15)
in which $`\psi (z)=\mathrm{Im}[h(z)]`$ and the subscripts indicate opposite sides of a phase boundary.

Let us now assume that in the thermodynamic limit, there is a phase transition at the point $`z_0`$. Then, either side of the transition the free energy can be written as

$$f_{1,2}(\stackrel{~}{z})=h(z_0)+b_{1,2}\stackrel{~}{z}+c_{1,2}\stackrel{~}{z}^2+\mathrm{}$$
(16)
where $`\stackrel{~}{z}=zz_0`$ and $`f_1(\stackrel{~}{z})`$ is valid for $`\mathrm{Re}\stackrel{~}{z}<0`$ and $`f_2(\stackrel{~}{z})`$ valid for $`\mathrm{Re}\stackrel{~}{z}>0`$. Note that for the free energy to be real along the real $`z`$ axis, the coefficients $`b`$ and $`c`$ must also be real. From the criterion that the real part of $`h(z)`$ is continuous across a phase boundary, we find that the boundary lies along the curve

$$\stackrel{~}{y}^2=\stackrel{~}{x}^2+\frac{b_2b_1}{c_2c_1}\stackrel{~}{x}$$
(17)
where $`\stackrel{~}{x}`$ and $`\stackrel{~}{y}`$ are the real and imaginary parts of $`\stackrel{~}{z}`$ respectively. We consider the conditions under which a discontinuous (first order) or continuous (second or higher order) appear.

First-order transition
For the case where $`b_1b_2`$, and the free energy has a discontinuity in its first derivative, the curve of zeros is a hyperbola that passes smoothly through the transition point $`z_0`$. Hence the tangent to the curve of zeros is parallel to the imaginary axis at $`z_0`$ and using the rule (15) we find that

$$\mu (0)=\frac{1}{2\pi }\frac{\mathrm{d}}{\mathrm{d}\stackrel{~}{y}}[(b_2b_1)\stackrel{~}{y}+2(c_2c_1)\stackrel{~}{x}\stackrel{~}{y}]|_{\stackrel{~}{x}=\stackrel{~}{y}=0}=\frac{b_2b_1}{2\pi }.$$
(18)
Hence we see that the density of zeros at the transition point $`z_0`$ is nonzero at a first-order phase transition (i.e. one at which the first derivative of the free energy is discontinuous).

Second-order transition
If $`b_1=b_2`$ but $`c_2c_1`$ we have that the curve $`C`$ obeys the equation $`\stackrel{~}{y}=\pm \stackrel{~}{x}`$. Since the zeros of the partition function come in complex-conjugate pairs, we find that the zeros approach the point $`z_0`$ along straight lines that meet at a right-angle. If we consider the line $`\stackrel{~}{x}=\stackrel{~}{y}=s/\sqrt{2}`$, we find from (15) that

$$\mu (s)=\frac{1}{2\pi }\frac{\mathrm{d}}{\mathrm{d}s}(c_2c_1)s^2=\frac{c_2c_1}{\pi }s.$$
(19)
This reveals that at a second-order phase transition, the density of partition function zeros decreases linearly to zero at the phase transition point.

Higher-order transition
More generally, one can consider a difference in the free energies either side of the transition point to have the leading behaviour $`f_2(u)f_1(u)u^\alpha `$. Then the condition that $`\mathrm{Re}[f_2(u)f_1(u)]=0`$ suggests that the curve of zeros approaches the real axis at an angle $`\frac{\pi }{2\alpha }`$. Unless $`\alpha =1`$, the curve does not pass smoothly through the real axis. In any case, the imaginary part of the free energy difference grows as $`|u|^\alpha `$ giving rise to a density of zeros that behaves as $`s^{\alpha 1}`$ for small arclength $`s`$. This means that for the density of zeros to be finite at the transition point we must have $`\alpha 1`$.

In this section we have summarised what we refer to as the Lee-Yang theory of phase transitions. This describes how partition function zeros are related to phase transitions: the accumulation of zeros at a point along the physical (real, positive) axis of the control parameter gives the critical value and the density of zeros near to the accumulation point determines the order of the phase transition.

We derived a rule (10) for locating phase boundaries and a further rule (15) for finding the density of zeros along such boundaries. At a first-order transition there is a nonzero density of zeros at the transition point whereas the density decays as a power-law to zero at the transition point when the associated phase transition is continuous.

Although we have discussed the theory of partition function zeros with reference to the specific system described by the partition function (1) we should note that the ideas hold much more generally. Firstly, one is not restricted to the canonical ensemble: indeed in the original exposition of the theory YL52 , Lee and Yang worked in the grand canonical ensemble in which the quantity generalised to the complex plane was not a function of the temperature but rather the chemical potential. Of course, in the equilibrium theory, these intensive field-like quantities play similar mathematical roles and so there is no reason why the Lee-Yang theory shouldn’t apply to all of them. For historical reasons, zeros in the complex fugacity (or chemical potential) plane are often referred to as Lee-Yang zeros YL52 ; LY52  and those in the complex temperature plane Fisher zeros Fisher65 . Also it appears that one can just as easily generalise physical ‘field-like’ variables (such as temperature) or ‘fugacity-like’ variables (such as the quantity $`z`$ considered above) to the complex plane without altering the properties of the partition-function zero densities at first-order and continuous phase transitions described above.

Despite the apparently wide generality of the Lee-Yang theory of equilibrium phase transitions, proving its validity in the general case is a difficult task (although we note recent work in this direction BBCKK00 ). Therefore, most rigorous results, such as those discussed in Ruelle69 , tend to rely on specific properties of a particular partition function. Perhaps the most spectacular of these is that obtained by Lee and Yang in their original work. Specifically, they found that the zeros of the partition function for an Ising ferromagnet in an external field $`h`$ all lie on the circle $`|\mathrm{exp}(h)|=1`$. The significance of this result, which does not depend on the number of spatial dimensions, lattice structure and details of the spin-spin interactions, is that if there is a phase transition induced by varying the magnetic field $`h`$, it can only occur at the point $`h=0`$, and then only if the partition function zeros accumulate there in the thermodynamic limit.

## III Nonequilibrium steady states

A nonequilibrium steady state differs from its equilibrium counterpart in that it may admit a flow of, say, energy or mass. More generally, these states have a circulation of probability in the space of microscopic configurations. In the past few years, it has become customary to model nonequilibrium systems as stochastic processes, i.e. simple models defined by local dynamical rules. Extensive study of these has revealed a range of phenomena including nonequilibrium phase transitions MD99 ; Evans00 ; Hinrichsen00 ; Mukamel00 ; Schutz01 ; Stinchcombe01 ; EB02 . We present below the key elements of this approach to nonequilibrium physics with a view to understanding nonequilibrium phase transitions in the framework of the Lee-Yang theory described above. In particular, we will need to propose a quantity to use in place of the equilibrium partition function (1).

Let us begin by discussing how one might realise the dynamics of the equilibrium spin system of the previous section, for example in a computer simulation. The aim is to generate a sequence of spin configurations such that the frequency with which a particular configuration $`𝒞`$ is generated is proportional to the Boltzmann weight $`f(𝒞)=\mathrm{exp}[\beta E(𝒞)]`$ where $`E(𝒞)`$ is the energy of configuration $`𝒞`$. Usually this is achieved by choosing the next configuration $`𝒞^{}`$ from the current configuration $`𝒞`$ with a probability proportional to the transition rate $`W(𝒞𝒞^{})`$ that satisfies the detailed balance condition

$$f(𝒞)W(𝒞𝒞^{})=f(𝒞^{})W(𝒞^{}𝒞)𝒞,𝒞^{}.$$
(20)
Apart from the fact that this choice of transition rates guarantees convergence to the desired equilibrium distribution of microstates NY01 , it also has the feature that the mean current of any quantity in the steady state is zero (as it should be at equilibrium).

Since we are interested in nonequilibrium systems that support steady-state currents we must work with transition rates that do not satisfy the condition (20) or equivalent criteria (given in e.g. vanKampen92 ; Mukamel00 ; EB02 ). Clearly one has a lot of freedom in the choice of transition rates, and so in practice one often devises dynamics that seem physically reasonable given the phenomena one is trying to describe. Later, in section IV.1 we will give concrete examples in the form of driven diffusive systems, in which particle moves are biased in a particular direction to model the effect of an external field.

In order to find the set of steady-state weights associated with a particular set of prescribed transition rates, we impose the condition that the total flow of probability into a configuration $`𝒞`$ is balanced by the corresponding outflow. That is, we must have

$$\underset{𝒞^{}𝒞}{}\left[f(𝒞^{})W(𝒞^{}𝒞)f(𝒞)W(𝒞𝒞^{})\right]=0$$
(21)
for every configuration $`𝒞`$. In general there might be more than one solution to this set of equations for the rates $`f(𝒞)`$, each corresponding to a steady state that is reached with some probability that depends on the starting configuration. We shall assume for simplicity that the steady state is unique and therefore reached with certainty from every initial state. A sufficient condition for a unique steady state is that it is possible to reach each configuration from every other via a sequence of microscopic transitions EB02 .

Once the steady state weights $`f(𝒞)`$ are known from (21) one can obtain the corresponding probability distribution of microscopic configurations $`P(𝒞)`$ through a normalisation

$$Z=\underset{𝒞}{}f(𝒞)$$
(22)
such that $`P(𝒞)=f(𝒞)/Z`$. Then, one can compute averages of physical quantities and look for nonanalyticities to locate nonequilibrium phase transitions as one varies the transition rates. At equilibrium we saw that the zeros of the function (1) that normalises the Boltzmann weights encodes the phase behaviour of the system. Thus we might hope that more generally, the zeros of the normalisation (22) will provide information about nonequilibrium phase behaviour.

In section IV we review recent work that suggests that the steady-state phase behaviour of certain nonequilibrium driven diffusive and reaction-diffusion systems is correctly described by the zeros of the normalisation $`Z`$. Although we are unaware of any rigorous argument for this to be the case, some suggestive evidence is provided by an observation made in Blythe01 ; EB02  which we now outline.

First note that the equation (21) is linear in the weights $`f(𝒞)`$. This implies that these weights, and hence the normalisation, can always be written as sums of products of the transition rates $`W(𝒞𝒞^{})`$. Now, by numbering the microscopic configurations, one can construct the transition rate matrix $`𝒲`$ whose off-diagonal elements $`𝒲_{nm}`$ are equal to the transition rates from configuration $`m`$ to $`n`$ and the diagonal elements $`𝒲_{nn}`$ are negative and express the total rate of departure from configuration $`n`$. Using elementary results from matrix theory (see Blythe01 ; EB02  for the details), one can show that an expression for the normalisation $`Z`$ that is polynomial in the transition rates $`𝒲(𝒞𝒞^{})`$ is

$$Z=\underset{\lambda _i0}{}(\lambda _i)$$
(23)
in which $`\lambda _i`$ are the eigenvalues of the matrix $`𝒲`$.

With each of the eigenvalues $`\lambda _i`$ is associated an eigenvector of $`𝒲`$ describing a ‘mode’ of the stochastic process that decays exponentially with a timescale $`\tau _i=1/|\lambda _i|`$ EB02 . We have assumed that the process described by $`𝒲`$ has only one steady state, and so only one of the eigenvalues is equal to zero (since clearly the relaxation time of a steady state is infinite). Equation (23) states that the normalisation $`Z`$ can be written as a product of the remaining eigenvalues. Now, at a phase transition we expect diverging timescales: physically one encounters metastable (long-lived) states near first-order transition points and long correlation lengths and times at continuous transitions. The presence of long timescales $`\tau _i`$ implies small eigenvalues $`\lambda _i`$ of $`𝒲`$ and, from (23), that the normalisation $`Z`$ approaches zero. Hence it appears that it is appropriate to consider the zeros of $`Z`$ as given by (22) to locate nonequilibrium phase transitions.

## IV Application of Lee-Yang theory to nonequilibrium phase transitions

We shall now review recent progress in applying the Lee-Yang theory to nonequilibrium phase transitions. We consider first in section IV.1 driven diffusive systems, where most work has so far been focussed Arndt00 ; BE02 ; Jafarpour02 . An appealing feature, discussed in more detail below, is that some one-dimensional cases have been solved exactly, the normalisation (23) calculated, and nonequilibrium phase transitions analysed. (The existence of one-dimensional phase transitions contrasts with the case of one-dimensional equilibrium models which do not admit phase transitions if the interactions are short-ranged.) Then in Section IV.2 we move onto reaction-diffusion systems and directed percolation.

### IV.1 Driven diffusive systems

In their original papers on partition function zeros, Lee and Yang YL52 ; LY52  made use of the mapping between the Ising model of a magnet and the lattice gas. Essentially one associates up-spins with particles and down-spins with vacancies so that a positive interaction strength $`J`$ in the Ising model results in a particle-particle attraction in the lattice gas, and a negative interaction strength gives rise to repulsive gas particles. Note that there is implicitly at most one particle per lattice site, so there is a hard-core exclusion in the lattice gas.

As discussed in section III, one can realise the dynamics of the lattice gas through a set of transition rates that satisfy the detailed balance condition (20) with respect to the Boltzmann distribution. Thirty years after Lee and Yang’s work, Katz, Lebowitz and Spohn (KLS) KLS83 ; KLS84  introduced a driven lattice gas model in which the rate at which particles hop in the direction of an external field is enhanced and the hop-rate against the field is suppressed. This model is well-studied and many results are discussed in SZ95 . The principal effect of the driving field is to introduce anisotropy in the phase-separation that occurs below the critical temperature associated with the spontaneous magnetisation of the Ising ferromagnet ($`J>0`$) in two or more spatial dimensions and in the critical exponents that characterise this transition.

As yet, the KLS model remains unsolved for general interaction strength, although in one dimension the steady state is known for some parameters KS92 . The particular case of one dimension and zero interaction strength $`J=0`$, known as the asymmetric simple exclusion process (ASEP), had already been studied—at least at a mean-field level—by biophysicists interested in the kinetics of biopolymerisation MGP68 . The mean-field approach predicts phase transitions in the steady state as parameters controlling the rate of insertion and extraction of particles at the boundaries are varied DDM92 . The existence of these phase transitions is confirmed through an exact solution of the ASEP DDM92 ; DEHP93 ; SD93 , achieved using a powerful matrix product approach DEHP93 ; Derrida98  which has subsequently been used to solve many generalisations of the ASEP. The details of the matrix product method are not necessary for the following—suffice to say that one ends up calculating a normalisation proportional to (23) through a product of matrices, often of infinite dimension. In this way one obtains some explicit formulas for the normalisation (22) which we shall use below <sup>1</sup><sup>1</sup>1As yet, there is no comprehensive review of applications of the matrix-product method, although one should consult Derrida98 ; Hinrichsen00 ; Schutz01 ; Stinchcombe01  and references therein.

The asymmetric exclusion process with open boundaries is perhaps the simplest exactly solved nonequilibrium model that exhibits both a first-order and continuous phase transition in its steady state. Therefore it is an ideal candidate for testing the hypothesis outlined in section III that zeros of the normalisation should accumulate towards the positive real axis in the complex plane of transition rates. Before outlining the results of this analysis (the details of which are presented in BE02 ) we recall the definition of the ASEP with open boundaries.

In this system, a particle on an $`N`$-site lattice can hop one site to the right at unit rate, as long as the receiving site is empty. Meanwhile, particles are inserted onto the leftmost lattice site (if empty) at a rate $`\alpha `$ and removed from the rightmost site (if occupied) at a rate $`\beta `$—see figure 2. Along the line $`\alpha =\beta <\frac{1}{2}`$ there is a coexistence between a high- and a low-density phase, with a shock separating the two. This is indicative of a first-order transtion. Meanwhile, along the lines $`\alpha >\frac{1}{2},\beta =\frac{1}{2}`$ and $`\beta >\frac{1}{2},\alpha =\frac{1}{2}`$ there is a continuous transition (i.e. one accompanied by diverging lengthscales) to a phase in which the particle current is independent of $`\alpha `$ and $`\beta `$. The phase diagram for the model is shown in figure 3.

The steady-state normalisation (22) for the ASEP with $`N`$ sites has been calculated DEHP93  as

$$Z_N(\alpha ,\beta )=\underset{p=1}{\overset{N}{}}\frac{p(2N1p)!}{N!(Np)!}\frac{(1/\beta )^{p+1}(1/\alpha )^{p+1}}{(1/\beta )(1/\alpha )}.$$
(24)
It is a simple matter to use a computer algebra package to solve this equation for its zeros in the complex-$`\alpha `$ plane at fixed $`N`$ and $`\beta `$. (Equivalently, one could look at the complex-$`\beta `$ zeros at fixed $`\alpha `$ since $`Z_N(\alpha ,\beta )=Z_N(\beta ,\alpha )`$.) In figure 4 plots of the zeros are shown for $`\beta =1`$, where a continuous transition occurs at $`\alpha =\frac{1}{2}`$, and for $`\beta =\frac{1}{3}`$, where a first order transition occurs at $`\alpha =\frac{1}{3}`$. We immediately notice that the curve of zeros seem to intersect the real positive $`\alpha `$ axis at the correct transition point. Furthermore, the density of zeros near the first-order transition point ($`\alpha =\beta =\frac{1}{3}`$) seems to be uniform and nonzero, whereas the density of zeros near the second-order transition point ($`\alpha =\frac{1}{2},\beta =1`$) seems to decrease to zero. Both of these observations are in accord with the results known for equilibrium partition function zeros discussed in section II.

As shown in BE02 , the distribution of zeros in the thermodynamic limit can be calculated once one knows that, for large $`N`$, the normalisation behaves as $`Z_NAJ^NN^\gamma `$. In this expression, $`A,J`$ and $`\gamma `$ depend on $`\alpha `$ and $`\beta `$ and the quantity $`J`$ is the current of particles across the lattice. Thus, the complex free energy (3) in the limit $`N\mathrm{}`$ is

$$h(\alpha ,\beta )=\underset{N\mathrm{}}{lim}\frac{\mathrm{ln}Z_N(\alpha ,\beta )}{N}=\mathrm{ln}J(\alpha ,\beta ).$$
(25)
Furthermore for $`\alpha `$ smaller than the transition value $`\alpha _c`$, $`J=\alpha (1\alpha )`$ whereas for $`\alpha >\alpha _c`$, $`J=\alpha _c(1\alpha _c)`$.

Although an electostatic analogy was used in BE02  to find the zero distribution, the mathematical content is the same as that used to derive the two rules (10) and (15) in section II. Applying the first rule, which demands that the real part of the free energy be continuous across a phase boundary, we find that the zeros of $`Z_N(\alpha ,\beta )`$ should lie along the curve $`|\alpha (1\alpha )|=\alpha _c(1\alpha _c)`$. It is therefore convenient to change variable to $`\xi =\alpha (1\alpha )`$. Then, in the complex-$`\xi `$ plane, the zeros lie on the circle $`|\xi |=\xi _c=\alpha _c(1\alpha _c)`$. The density of zeros on this circle can be found by setting $`\xi =r\mathrm{e}^{i\theta }`$ and parametrising points on the circle as $`s=\xi _c\theta `$. Then, the second rule (15) gives the density of zeros $`\mu (s)`$ on the circle as

$$\mu (s)=\frac{1}{2\pi }\frac{\mathrm{d}}{\mathrm{d}s}\left[\mathrm{Im}\mathrm{ln}\xi \mathrm{Im}\mathrm{ln}\xi _c\right]=\frac{1}{2\pi \xi _c}\frac{\mathrm{d}}{\mathrm{d}\theta }\theta =\frac{1}{2\pi \xi _c}.$$
(26)
That is, in the complex-$`\xi `$ plane the zeros should become evenly distributed on a circle in the thermodynamic limit. Transforming the zeros of (24) obtained at different system sizes to the complex-$`\xi `$ plane reveals this to be the case BE02 .

Finally, one can show that near the intersection point between the curve of zeros in the complex-$`\alpha `$ plane and the positive real $`\alpha `$ axis, the zeros of (24) sit on the curve $`x=\frac{1}{2}(y^2+\frac{1}{4}\xi _c)^{1/2}`$ where $`x`$ and $`y`$ are the real and imaginary parts of $`\alpha `$. For the case $`\beta <\frac{1}{2}`$, $`\xi _c<\frac{1}{4}`$ and the transition is first order. One finds that the curve of zeros is smooth at the transition point $`\alpha =\beta `$, and that the density of zeros is $`(12\beta )/[2\pi \beta (1\beta )]`$ which is nonzero. These are precisely the properties of the equilibrium partition function zeros at a first-order transition point (see section II).

At the continuous transition point ($`\beta \frac{1}{2}`$ and $`\xi =\frac{1}{4}`$, the zeros pass through the real $`\alpha `$ axis along the line $`x=\frac{1}{2}|y|`$. Recall from section II that a transition of $`n^{\mathrm{th}}`$ order has the curve of zeros meeting the positive real axis at an angle of $`\frac{\pi }{2n}`$. Here the zeros clearly approach at an angle $`\frac{\pi }{4}`$ suggesting a second-order transition. In fact, one finds the density of zeros is $`\frac{4}{\pi }s`$ at a distance $`s`$ along this line, confirming that the transition is second-order.

In summary, we have found that the Lee-Yang theory of first-order and continuous phase transitions applies to the normalisation of the nonequilibrium asymmetric exclusion process just as it does to the partition function of equilibrium systems. Of course, this does not prove that the theory is generally applicable, and so there is some value in investigating other nonequilibrium steady states that exhibit phase transitions.

One such state is that initially studied by Arndt, Heinzel and Rittenberg AHR98 ; AHR99 . This model comprises $`M_+`$ ($`M_{}`$) positively (negatively) ‘charged’ particles on a closed ring of $`N`$ sites. When next to vacant sites, the positive particles hop to the right and negative particles to the left at unit rate, with hops in the opposite directions disallowed. Meanwhile, should two oppositely charged particles be next to one another, the following transitions can occur

$$+\stackrel{q}{}++\stackrel{1}{}+$$

where the label above the arrow indicates the rate at which the transitions occur. These dynamics are illustrated in figure 5.

Although a matrix product solution for this model is known to exist AHR99 , it is technically difficult to work out the solution in its full generality RSS00 . However in the case of either a single vacancy sasamoto00  or a single negative particle, i.e. $`M_{}=1`$ Jafarpour00  the steady-state normalisation is known explcitly and phase transitions have been identified. In the case $`M_{}=1`$ the zeros of the partition function in the complex-$`q`$ plane have been studied Jafarpour02  and, once again the continuous phase transition present in this model is correctly described by the Lee-Yang theory.

A slightly different Lee-Yang analysis, that precedes the work described above, has been used to study the nonequilbrium phase transition that occurs at a finite density of vacancies Arndt00 . The motivation behind this study arose from computer simulations of the model AHR98 ; AHR99  for the case where the numbers of positive and negative particles were equal. In these simulations, three phases were initially identified: a ‘solid’ phase for small $`q`$, a ‘fluid’ phase for large $`q`$ and a mixed phase comprising the background ‘fluid’ and a large mobile droplet in an intermediate regime $`1<q<q_c`$ where $`q_c`$ is some density-dependent quantity.

The steady-state normalisation $`Z_{N,M}`$ for the case of $`M_+=M_{}=M`$ positive and negative particles on the $`N`$-site ring is not known exactly (at least in the canonical ensemble), therefore Arndt Arndt00  chose to study the generating function

$$𝒵_N(z,q)=\underset{M}{}z^MZ_{N,M}(q)$$
(27)
in the complex-$`z`$ plane. Physically this approach is equivalent to placing the ring with nonequilibrium interactions in chemical equilibrium with a particle reservoir at fixed fugacity $`z`$. One can check that as one takes the thermodynamic limit, the relative size of the fluctuations in the density $`r=M/N`$ vanish, and so working at fixed fugacity $`z`$ is equivalent to working at fixed density $`r`$. (In fact, this is a standard ‘trick’ for dealing with closed systems in which the particle number is conserved, see e.g. DJLS93 .)

The key point here, however, is that $`z`$ is an equilibrium fugacity, and not a microscopic transition rate, and so the Lee-Yang theory of phase transitions described in section II ought to apply directly here, without reference to the discussion of section III. It was found Arndt00  that when $`q<q_c`$ the zeros of (27) in the complex-$`z`$ plane appear to lie along ellipses, intersecting the positive real axis. Conversely, for larger $`q`$ the zeros appear to describe hyperbolæ that avoid the positive real axis. After investigating numerically the density of zeros in the complex-$`z`$ plane, Arndt concluded that there is a fluid-fluid phase transition at some density-dependent point $`q_c>1`$.

Unfortunately, an exact asymptotic (i.e. large $`N`$) analysis of the generating function (27) shows that the only nonanalyticity occurs at the solid-fluid transition point $`q=1`$ RSS00 . However it was noted that physical quantities vary rapidly, but continuously, at some a point $`q>1`$. This phenomenom has been explained as an abrupt increase in a corrleation length to an anomalously large, but finite value KLMT02 . We believe that resolution with the study of Arndt lies in the fact that the system sizes considered in Arndt00  were quite small ($`N100`$), whereas it was suggested RSS00  that the distinction between this crossover behaviour and a genuine phase transition might become apparant only once $`N`$ is increased above about $`10^{70}`$ (an unfeasibly huge number!).

It would be interesting, therefore, to extend the numerical computation of zeros performed in Arndt00  to much larger systems, to see how the ellipses noted above develop. For example it might well happen that instead of approaching the positive-real fugacity axis, the zeros of (27) would terminate a short distance away from it.

### IV.2 Reaction-diffusion systems and directed percolation

A large and important class of models with stochastic dynamics is provided by reaction-diffusion systems <sup>2</sup><sup>2</sup>2As an area of ongoing research, new results are emerging continuously and there is no up-to-date review that contains them all. Nevertheless, one can consult Redner97 ; Hinrichsen00  for an introduction.. In contrast to driven diffusive systems, where the particle-particle interactions imply conservation of particles, reaction-diffusion systems are characterised by dynamics that result in a change in particle number. Moreover there are a number of such systems that have absorbing states i.e. special configurations generally devoid of particles that once entered cannot be left. Phase transitions associated with whether the system has a finite probability of not being absorbed into such a state fall within the directed percolation and related universality classes Hinrichsen00 . We shall shortly discuss the second order phase transition associated with the directed percolation university class in a little detail. Meanwhile as a simple example of a reaction-diffusion system, we review a model for which the steady state that can be solved using the matrix product approach. The approach again provides us with an explicit expression for the normalisation (22) by virtue of which we can analyse its zeros in the plane of complex reaction rates.

The system in question HSP96  has for the dynamics at neighbouring bulk sites on a one-dimensional lattice the processes

$$\begin{array}{cccccc}& \stackrel{q}{}& & & \stackrel{q^1}{}& \\ & \stackrel{q}{}& & & \stackrel{q^1}{}& \\ & \stackrel{\kappa q}{}& & & \stackrel{\kappa q^1}{}& \end{array}$$

occurring at the rates indicated, and with $``$ representing a particle and $``$ an empty lattice site. It was demonstrated HSP96  that the matrix product scheme used for the ASEP could be generalised to cater for the steady state of the present reaction-diffusion system on a lattice of $`N`$ sites with reflecting boundary conditions (i.e. which is neither periodic nor has particle input or removal at the boundaries).

In the matrix product approach, the normalisation $`Z_N`$ is given by a scalar derived from $`C^N`$ where $`C`$ is a square matrix. If $`C`$ is a finite dimensional $`m\times m`$ matrix and can be diagonalised, the resulting expression for $`Z_N`$ has the form

$$Z_N(q,\kappa )=\underset{n=1}{\overset{m}{}}a_n\lambda _n^N$$
(28)
in which $`\lambda _n`$ is the $`n^{\mathrm{th}}`$ eigenvalue of $`C`$ and the coefficients $`a_n`$ arise from the details of the way in which one obtains a scalar from the matrix $`C^N`$. A normalisation of this form leads to a complex free energy

$$h(q,\kappa )=\mathrm{ln}\left[\underset{n}{\mathrm{max}}\{\lambda _n\}\right]$$
(29)
in which the maximum means the eigenvalue with the largest absolute value. Then, there is the possibility of a phase boundary when the magnitude of the two largest eigenvalues of $`C`$ are equal. It turns out that for the closed reaction-diffusion system introduced above, $`C`$ is a $`4\times 4`$ matrix that can be diagonalised when $`q^21+\kappa `$ HSP96 . When $`q^2=1+\kappa `$, $`C`$ cannot be diagonalised on account of the largest eigenvalue being degenerate, and a phase transition is found to occur at this point. Note that this scenario contrasts with the transfer-matrix approach to one-dimensional equilibrium systems where the partition function is also written as a product of matrices. Since all elements of the transfer matrix are positive the largest eigenvalue cannot become degenerate therefore there can be no phase transition Evans00 . However there is no such restriction on the elements of $`C`$, and so eigenvalue crossing is permitted and nonequilibrium one-dimensional phase transitions can occur.

The structure of the normalisation $`Z_N(q,\kappa )`$ is even simpler in the case where particles are inserted onto and removed from the left boundary at rates $`\alpha `$ and $`\beta `$ respectively (with the right boundary remaining reflecting), and these rates satisfy the relation $`\alpha =\kappa (q^1q+\beta )`$ Jafarpour03 . In this case $`C`$ is a $`2\times 2`$ matrix and, for $`q^21+\kappa `$, is diagonalisable with eigenvalues

$$\lambda _1=1+\kappa \text{and}\lambda _2=q^2.$$
(30)
We see then immediately from (29) and an application of the rule (10) that there is a phase boundary in the complex $`q`$ plane along the circle $`|q|=\sqrt{1+\kappa }`$. Also, since one of the eigenvalues does not depend on $`q`$, the free energy $`h(q,\kappa )`$ is a constant in one of the phases which, as the analysis of the normalisation zeros for the ASEP demonstrated, implies that the density of zeros on this circle is constant in the thermodynamic limit. In turn, this implies that the phase transition is first order, as confirmed by explicitly calculating the density profile in the two phases Jafarpour03 .

We finally turn our attention to a process that comprises symmetric decoagulation (that is, $``$ and $``$ taking place with equal probability in each direction) and spontaneous decay ($``$) occurring at independent rates. Introduced as a crude model of an epidemic Harris74 , this contact process is known to exhibit a transition from a phase in which the absorbing state (empty lattice) is reached with certainty to a phase in which there is some probability that the epidemic remains active forever in the thermodynamic (infinite system size) limit Griffeath79 . This transition occurs as the ratio between the decoagulation and decay rates is increased beyond a critical value.

Although not proven, it is widely accepted on the basis of simulation and approximate methods (such as series expansions) that the transition just described is continuous and characterised by directed percolation (DP) exponents (see Hinrichsen00  for an in-depth discussion of these issues). Directed percolation itself BH57  is a geometric construction designed to model fluid flow through a random porous medium under the influence of gravity. Consider the rhombic lattice shown in figure 6 and imagine ‘pouring’ fluid in from the uppermost site. The fluid is allowed to flow along a bond pointing diagonally downwards and then only if it is open. Each bond has a probability $`p`$ of being open (and hence a probability $`1p`$ of being blocked), and the state of each bond is independent of any other and may not change with time.

The main quantity of interest in this system is the percolation probability $`P_n(p)`$ that the fluid can penetrate to a depth $`n`$ in the medium (see figure 6), averaged over all possible bond configurations. The order parameter for the model is the probability $`P_{\mathrm{}}(p)`$ that the fluid can penetrate infinitely far given a bond probability $`p`$. If $`p`$ is below some critical percolation threshold $`p_c`$, the fluid only ever penetrates a finite distance—i.e. a layer becomes dry with certainty. However, for larger $`p`$, the order parameter becomes nonzero, growing as $`P_{\mathrm{}}(p)|pp_c|^\beta `$ with $`\beta 0.276`$ near the critical point. The importance of the DP transition is that a wide range of models that have a transition into an absorbing state are expected to have that transition characterised by the DP exponents, one of which is $`\beta `$ Hinrichsen00 . Despite a huge amount of interest in directed percolation, none of the models expected to belong to its universality class has been solved exactly.

Recently an attempt has been made to shed further light on the DP transition by studying the zeros of the percolation probability $`P_n(p)`$ ADH01 ; DDH02 . At first glance, a connection between this probability and a partition function is not obvious. However, associating with each site a ‘spin’ state $`\sigma _i`$, that is up if site $`i`$ is connected to a point on the $`n^{\mathrm{th}}`$ layer and down otherwise, yields a form for $`P_n(p)`$ that has the structure of a transfer-matrix representation of a partition function for an equilibrium system with three-spin interactions BG88 .

It is not clear that being able to give $`P_n(p)`$ the appearance of a partition function necessarily implies that the Lee-Yang theory should hold. In particular, $`P_n(p)`$ has some features that set it apart from ‘standard’ equilibrium partition functions, in that it does not grow exponentially in the number of bonds $`N`$, at least in the range $`0p1`$. It is also uncertain whether a free energy can meaningfully be defined for this system, especially in the region $`0pp_c`$ in which $`lim_n\mathrm{}P_n(p)=0`$.

Given these observations, it is perhaps not surprising that the zeros of $`P_n(p)`$ calculated numerically for $`n15`$ (leading to polynomials of degree up to $`N=240`$) have a more complicated distribution than for the (exactly solved) models considered so far. As is evident from figure 7, which was presented in DDH02 , the zeros do not lie along a single curve but along a sequence of curves meeting at the critical point (as well as inside some region that encloses part of the real axis for $`p>2`$).

By considering the sequences of zeros approaching the positive real axis, one can estimate the critical point $`p_c`$ and the density of zeros as it is reached DDH02 . The latter yields a prediction for one of the DP critical exponents and one finds good agreement with the most precise estimates of both the transition point and the associated exponent obtained through other means. Hence we have further evidence for the applicability of Lee and Yang’s ideas concerning partition function zeros to a much wider range of statistical distributions than equilibrium steady states.

## V Summary and outlook

In this work we have revisited the Lee-Yang description of equilbrium phase transitions with a view to seeing whether the ideas apply to more general nonequilibrium transitions. Recently there have been a number of studies of zeros of partition-function-like quantities that arise in systems with nonequilibrium dynamics, and we have seen in our review of these works that the Lee-Yang theory, as described in section II, seems to hold quite generally.

We have argued that for dynamic models with a unique steady state, the normalisation defined as a sum over the steady-state configurational weights (22) serves as a suitable ‘partition function’ in the sense that its zeros, in the complex plane of any model parameter, should accumulate towards physical transition points in the thermodynamic limit. Furthermore, the density of zeros and angle of approach to the real axis indicate whether the transition is first-order (manifested physically through phase-coexistence) or continuous (i.e. characterised by divergent correlation lengths and times). Thus studying the zeros of the normalisation (22) provides an unambiguous classification of nonequilibrium phase transitions as do the Lee-Yang zeros in the equilibrium case.

The observation that backs up this scenario is embodied by equation (23) which reveals that the reciprocal of the steady-state normalisation is equal to the product of the characteristic relaxation times in the dynamics. Since near a phase transition one expects timescales to diverge, one also expects the normalisation to approach zero. However, a rigorous argument for this to be the case is still lacking. Moreover (23) implies a possible link between systems for which the steady state normalisation can be calculated and those for which eigenvalues of the transition matrix can in principle be calculated.

A different class of systems encompasses those whose steady state is not unique. The contact process is, in fact, an example of such a model, in which the absorbing state is reached with certainty below the critical decoagulation rate, whereas above it, and on an infinite system, a second steady state can also be reached with some nonzero probability. Since this additional steady state exists only when the lattice size becomes infinite, one must take that limit first, before taking time to infinity. Otherwise, on a finite system the steady state is simply the absorbing state and the steady-state normalisation is trivially equal to a constant. Nevertheless the work of ADH01 ; DDH02 , which we reviewed in section IV.2, indicates that the Lee-Yang theory can be relevant when one considers other properties of the nonequilibrium system such as the percolation probability.

Although we have not discussed this in great detail here, it should be noted that the Lee-Yang approach gives a method for extrapolating to the thermodynamic limit from solutions for small system sizes. In the work of DDH02 , numerical solutions for small systems were used successfully to estimate the transition point and density of zeros as it is approached. From this information one learns about the nature of the phase transition and, for example, can estimate the values of critical exponents. This is a common technique in equilibrium statistical physics (see e.g. JK02 ) and it may be the case that stochastic processes unyielding to analytical treatment could be understood this way.

###### Acknowledgements.

We thank Bernard Derrida for helpful discussions. R.A.B. acknowledges financial support under EPSRC grant GR/R53197.


format: text

text:
# Rapid Variations of Narrow Absorption Line Components in the Spectrum of P Cygni

## 1 Introduction; Long-Period DACs

P Cyg is a blue hypergiant (de Jager 1998) or a Luminous Blue Variable (LBV, Humphreys & Davidson 1994) whose observed characteristics make the star among the most peculiar ones in the Galaxy. The properties of P Cyg and its wind have been reviewed by Israelian and de Groot (1999).

Narrow absorption line components called DACs were discovered in Balmer line profiles of P Cyg by de Groot (1969). He also found them in lines of He, Fe and many other elements. Identifications of other lines were published by others. Cassatella et al. (1979) and Luud and Sapar (1980) detected DACs in UV lines of various ions, such as Cr ii, Ni ii, and others. Lamers et al. (1985) found them in Fe ii lines.

In his Balmer line observations de Groot found velocities of –215, –160 and –95 $`\mathrm{km}\mathrm{s}^1`$. The first component showed periodic $`v_R`$-variations of $`P=`$ 114 days, while the others did not present temporal variations. The observations by Luud et al. (1975) however, give contrasting results. They found that the first component was not variable, while the second was, with a period of 57 d. Since, many characteristic times of variability (we prefer this term above the too suggestive ’period’) have been mentioned: 200 d (Markova, 1986), 60 - 75 d (van Gent and Lamers, 1986), 4 to 5 months (Markova, 1993) and 200 d (Israelian et al., 1996). The apparent diversity in the communicated ’periods’ may partly be related to insufficient lengths of periods of observations; moreover, it may be a real phenomenon.

As to the interpretation of these motions, current thinking converges on a model according to which shells are periodically ejected from the star, while increasing their velocity on the way outward through the stellar wind. This hypothesis of expanding shells was proposed first by Kolka (1983), Markova and Kolka (1984) and Markova (1986). They were supported by van Gent and Lamers (1986), and are also shown in a diagram published by Israelian et al. (1996), where three or four shells appear, of which the velocity increases steadily from about 80 to 190 $`\mathrm{km}\mathrm{s}^1`$ over a period of about 1000 days.

Most of the above mentioned variations have characteristic times of the order of 2 to 7 months. There is another component, though, the period of which is not related to such presumed outgoing shells. De Groot et al. (2001) published the results of a fairly complete survey of P Cyg photometry in which they demonstrate that at each time there are two periods in the light variation. The shorter one, of 17.3 d, occurred persisitently, at least during the last 20 years albeit with variable amplitude. The longer period ranges between about 60 and 130 d, with a clear maximum near 100 d. It does not seem impossible that this latter variability is associated with the ougoing shells described above. But what to say about the shorter period ? The question arises if there is any wavelike phenomenon in the wind that could be associated with the shorter (17.3 d) light variation. Its study is the main content of this paper.

## 2 The detection of short-period absorption components

Observations of P Cyg were carried out during the week 1999, May 28–June 4, using the SOFIN Echelle Spectrograph of the 2.6-m Nordic Optical Telecope (NOT) at the ORM (La Palma). A resolving power of $`R=\lambda /\mathrm{\Delta }\lambda `$ 80 000 has been achieved with the second camera. The CCD used in this run was a $`1152\times 770`$ pixel2 EEV. We obtained echelle spectra with a slit width of 0.8 arcsec. They cover the wavelength range between 3600 and 4200 Å in 18 orders. All spectra were flat-fielded, background and scattered light subtracted. The wavelength calibration was performed with a Th–Ar lamp. The ultimate signal-to-noise of our spectra at 4000 Å is about 200.

Apparently, one week is too short for observing recurrent variability or for checking a periodicity of 17.3 d. In that respect our results have a preliminary character. The positive aspect is, though, that our observations (Fig. 1) do show the presence of rapidly changing absorption components in profiles of the higher members of the Balmer series.

The lines showed split absorption cores. These are best visible in the profiles of 28/05/99 (solid lines in Fig. 1), which have two distinct cores. Closer examination shows that these cores existed also in later spectra. The cores approached each other during the next few days and nearly merged around 02/06/99. Thereafter they reappeared. We measured the wavelength variation of these absorption cores by two methods: directly by determining the central wavelengths of the absorption component, and also by fitting the profiles by two Gaussians. This was done in the two lines, so we got 4 data points for each componentat each day. The measurements agree reasonably well; the mean scatter of the individual velocity measurements is 3 $`\mathrm{km}\mathrm{s}^1`$. Table 1 gives the average velocities measured in the two line components.

We wanted to check if these velocities are related in some way to the photometric period of 17.3 d, reported in Section 1. To that end we replotted the data in the same diagram for a period, 17.3 days later (Fig. 2). The diagram shows that it is not impossible that the observations have a repetition period of 17.3 d, and it appears moreover that there are possibly two waves that are ’swinging in counterphase’, cf. the dotted lines in Fig. 2. If this interpretation would be verified by measurements taken over a period of a month or longer, then they could be interpreted in terms of a long wave in the wind over the stellar surface having a wavelength of the order of the star’s diameter. This hypothesis will be examined in the next section.

An interesting aspect of the observations is the high acceleration that is inferred from the wavelength variation of the absorption components. It is about 10 $`\mathrm{km}\mathrm{s}^1\mathrm{d}^1`$. This value stands in contrast to the acceleration of the long-period DACs, which is about 100 times smaller, of the order of 0.1 $`\mathrm{km}\mathrm{s}^1\mathrm{d}^1`$.

## 3 A Possible Interpretation

In this section we examine the following question: if a periodic motion with $`P=17.3`$ d occured in the wind of P Cyg, what would be the physical character and the wavelength of that motion?

We base our answer on the dispersion equation for perturbations in line driven winds as developed by Abbott (1980, in his eq. (45)). That equation reads

$$\omega ^3+\omega ^2k_2F_2\omega (a^2k_{x}^{}{}_{}{}^{2}+a^2k_{z}^{}{}_{}{}^{2})a^2k_{x}^{}{}_{}{}^{2}k_zF_2=0$$
(1)
with

$$F_2=\delta g_{\mathrm{rad}}/(\delta (\delta w/\delta r))$$
(2)
where $`w`$ is the wind velocity, $`a`$ that of sound, $`k=2\pi /L`$; $`L`$ is the wavelength, $`\omega =2\pi /P`$ , where $`P`$ is the period of the disturbance and $`g_{\mathrm{rad}}`$ is the radiative acceleration. For purely horizontal or vertical waves $`k_z=0`$ or $`k_x=0`$ respectively.

For the calculation of the diagram we used the model of the wind of P Cyg given by Najarro et al. (1997), which has the main parameters $`T_{eff}`$ = 18 200 K, $`\mathrm{log}L/L_{}`$ = 5.75 and $`n_{He}/n_H=`$ 0.3. The physical conditions at the level where the short-period absorption components are formed were derived by using the observation (cf. Fig. 1) that these absorption components are formed at the level where the stellar wind velocity is between 30 and 100 $`\mathrm{km}\mathrm{s}^1`$. We took an average value of 60 $`\mathrm{km}\mathrm{s}^1`$. From Najarro et al.’s velocity law we then find that this level corresponds with $`r=2.1R_{}`$. The density and pressure at that level follow from Najarro et al.’s rate of mass loss. We assumed for that level a temperature of 12,000 K, an assumption that apeared not to be critical. That yields the continuous absorption coefficient. The value of $`\mathrm{\Gamma }_1`$ was computed using a non-LTE approach involving the 16 most abundant elements, developed by Lobel (in prep).

At that level $`R_{}=\mathrm{5.24\hspace{0.25em}10}^{12}`$ cm, $`g_{\mathrm{rad}}=22.49\mathrm{cm}\mathrm{s}^2`$ (cf. Lamers and Cassinelli, 1999) and $`\delta (\delta w/\delta r)=\mathrm{3.626\hspace{0.25em}10}^{20}\mathrm{s}^1`$. Hence $`F_2=2.15210^8\mathrm{cm}\mathrm{s}^1`$.

The two asymptotic cases are easiest to handle. With $`k_z=0`$, hence assuming strictly horizontal waves, and taking $`P=17.3`$ d we obtain for the horizontal wavelength $`L_x=\mathrm{1.982\hspace{0.25em}10}^{12}`$ cm or $`0.38R_{}`$. For purely vertical waves ($`k_x=0`$) we get $`L_z=\mathrm{1.22\hspace{0.25em}10}^{10}`$ cm or $`0.0023R_{}`$.

The waves are gravity waves since for $`P=17.3`$ d the corresponding value of $`\omega `$ falls below the Brunt-Väisälä frequency (de Jager, 2001). Since gravity waves mostly travel horizontally we tend to prefer the first solution, that for horizontal waves. An additional argument for choosing that option is that the wavelengh ( 0.38$`R_{}`$) is close to the wavelength suggested by the observations, for which we concluded in Section 2 that the waves should have lengths of the order of the stellar radius.

The conclusion of this Section is that if the short-period absorption components would have a period of 17.3 d, then they could be interpreted as horizontally propagating gravity waves with wavelengths of the order of 0.4 times the stellar radius.

## 4 Conclusions

The high S/N and high-resolution observations reported here show for the first time the existence of short-period narrow absorption components in the cores of high Balmer lines and their rapid variations. If their period was 17.3 d, which is the only known short period observed in the light of P Cyg, they should be interpreted as gravity waves with wavelengths of the order of 0.4 times the stellar radius, a conclusion that is also suggested by the observed behaviour of the absorption components, which seem to consist of two waves in counterphase. A period differing by a factor 2 or so would not change these conclusions. But, evidently, more observations are needed to check this statement and we urge observers to monitor these features.

We thank M. Friedjung, A. van Genderen, A. de Koter, O. Stahl and R. Stothers for help and useful remarks. Many thanks are due to H. Nieuwenhuijzen for advice and efficient help with part of the computations and Ilia Ilyin for help with SOFIN data reduction.


text:
# Quantum Interference to Measure Spacetime Curvature: A Proposed Experiment at the Intersection of Quantum Mechanics and General Relativity

\]physics.berkeley.edu/research/chiao/ thanks:

thanks:  (April 7, 2003)

## Abstract

An experiment in Low Earth Orbit (LEO) is proposed to measure components of the Riemann curvature tensor using atom interferometry. We show that the difference in the quantum phase $`\mathrm{\Delta }\varphi `$ of an atom that can travel along two intersecting geodesics is given by $`mR_{0i0j}/\mathrm{}`$ times the spacetime volume contained within the geodesics. Our expression for $`\mathrm{\Delta }\varphi `$ also holds for gravitational waves in the long wavelength limit.

Keywords: Riemann Curvature Tensor, Atom Interferometry

In general relativity no “gravitational force” exists Misner et al. (1973); there is only geometry. Thus, in the absence of any forces, a test particle will undergo free fall, and will travel along a geodesic determined by the local spacetime geometry. What we normally associate with the “force of gravity” on a particle in Newtonian mechanics is not a force at all: The particle is simply traveling along the “straightest” possible path in a curved spacetime.

Nowhere is this basic difference between Newtonian gravity and general relativity more apparent than for tidal forces. Consider twin astronauts in separate, but nearby, orbits around the Earth. In Newtonian gravity the astronauts accelerate towards or away from one another due to the difference in the Earth’s gravity between them. In general relativity, by contrast, because both astronauts do not feel any forces, they undergo free fall; they travel along geodesics in a spacetime curved by the Earth. Because the two geodesics are slightly different, one of the astronauts will see his twin drift toward or away from him. The rate of change of the distance $`X^\mu `$ separating the two astronauts depends on the local curvature tensor $`R_{\mu \nu \alpha \beta }`$ along $`X^\mu `$. If $`R_{\mu \nu \alpha \beta }`$ changes slowly along $`X^\mu `$, $`X^\mu `$ satisfies the geodesic deviation equation Misner et al. (1973); if $`R_{\mu \nu \alpha \beta }`$ changes appreciably along $`X^\mu `$, its evolution depends on the integral of $`R_{\mu \nu \alpha \beta }`$ along $`X^\mu `$ GLF (2002). In either case, one astronaut interprets the shift in his twin’s position as being due to a fictitious Newtonian “gravitational tidal force” acting on his twin, even though the latter feels no forces.

The twins age differently due to their differing gravitational red shifts: The twin in the lower orbit ages less than the farther twin. However, in this general relativistic twin paradox, both twins feel no forces, unlike in the special relativistic case. Their age difference cannot arise within a Newtonian framework, where time is absolute. However, at the quantum level, this age difference for interfering atoms becomes a measureable phase difference directly related to the curvature.

In this Gravity Research Foundation essay, we propose an experiment in LEO to measure components of the Riemann curvature tensor. In a quantum extension of the above example, this experiment will test quantum mechanically the tenant of general relativity that local geometry determines motion. However, instead of relying on classical test particles, we make use of a quantum phenomenon: the fringe shifts of an interference pattern caused by phase differences of quantum test particles. Consequently, this experiment will probe the intersection of general relativity and quantum mechanics Chiao (2003).

The conceptual roots of this proposed experiment are well established Stodolsky (1979); Anandan (1979); GLF (2002). Its technical roots lie in Chu’s PCC (2001) and Kasevich’s MFFSK (2002) work on atom interferometry in their measurements of the local acceleration due to gravity $`g`$ using Mach-Zehnder-type interferometers CarMly (1991); KETP (1991); RKWHB (1991); KasChu (1992); ROBSZ (1991). An accuracy based on Newtonian dynamics of $`\mathrm{\Delta }g/g10^8`$ for a single measurement, and up to $`10^{11}`$ after a two-day integration, was achieved. In these measurements, cesium atoms were thrown upwards in the Earth’s gravitational field. STIRAP (Stimulated Raman Processes) were used as beam-splitters, and the relative phase of the atom traveling along the two arms of the interferometer was measured. Because the interferometer was fixed to the Earth, this phase difference is proportional to $`g`$. We propose a similar type of experiment, but now in LEO where the interferometer is in free fall SBD (1999). The term linear in $`g`$ now disappears, and the phase difference will be proportional to the Riemann curvature tensor inside the area encircled by two geodesics.

Consider an atom interferometer in LEO at a distance $`R^i(t)`$ from the center of the Earth. We choose a local coordinate system $`X^\mu `$ fixed on the center-of-mass (CM) of the apparatus. The signature of $`g_{\mu \nu }`$ is $`(1,1,1,1)`$; in linearized gravity $`g_{\mu \nu }=\eta _{\mu \nu }+h_{\mu \nu }`$. Note also that $`|X^i||R^i|`$; we expand $`h_{\mu \nu }`$ about $`X^i=0`$.

We then release an atom traveling along the orbit of the CM with geodesic $`\gamma _1`$ (see Fig. 1). At $`t_A`$, STIRAP is used to coherently split the atomic beam. One possible geodesic for the atom is still $`\gamma _1`$; the other corresponds to a different geodesic $`\gamma _2`$. Note that the atom behaves quantum mechanically, and that the superposition principle holds. It is not possible to determine which geodesic any individual atom will take. Because the spatial projection of both $`\gamma _1`$ and $`\gamma _2`$ correspond to LEOs, they will intersect with one another again after a time $`T`$ as shown in Fig. 2. Detectors can then be used after coherent recombination at $`t_B`$ to determine the interference pattern, and thus the phase shift $`\mathrm{\Delta }\varphi `$ that the atom picks up between the two possible geodesics. The loss of fringe visibility would be a measure of decoherence of the atom. Importantly, the combined spacetime path $`\gamma =\gamma _1\gamma _2`$ is closed (see Fig. 1), and forms the boundary of a spacetime surface $`𝒟`$. If the coherence time of the atom is too short to allow the atomic beam to recombine due to drift, additional STIRAP beam splitters can be used to force recombination. However, let us assume that the coherence time is long enough to allow Fig. 2.

The Schrödinger equation for the atom is

$$i\mathrm{}\frac{\psi }{t}=\frac{\mathrm{}^2}{2m}^2\psi +i\mathrm{}N_i^i\psi mN_0\psi ,$$
(1)
where $`\psi `$ is the wavefunction, $`m`$ is the mass of the atom, and using standard methods,

$`N_0`$ $`={\displaystyle \frac{1}{4}}X^iX^j_i_jh_{00}(R(t))+{\displaystyle \frac{1}{2}}X^jX^k{\displaystyle \frac{dR^i}{dt}}_j_kh_{0i}(R(t)),`$
$`N_i`$ $`=X^j_jh_{0i}(R(t)){\displaystyle \frac{1}{2}}X^jX^k_j_kh_{0i}(R(t)).`$ (2)
Our results here hold only for stationary metrics. However, from GLF (2002) and Speliotopoulos (1995) our result for $`\mathrm{\Delta }\varphi `$ given below holds as long as the local curvature varies slowly in $`X^\mu `$. From GLF (2002), $`N_\mu =(N_0,N_i)`$ is the four-velocity field acting on the test particle induced by the tidal field as seen by an observer at the CM.

An atom forms a wave-packet that propagates along either $`\gamma _1`$ or $`\gamma _2`$. Taking the eikonal approximation $`\psi =e^{imS/\mathrm{}}\psi _0`$, where $`\psi _0`$ is the solution of eq. $`(\text{1})`$ in the absence of $`(N_0,N_i)`$,

$$\frac{S}{t}=\frac{1}{2}|S|^2+N_i^iS+N_0,0=\frac{\psi _0}{2}^2S+(_iSN_i)^i\psi _0.$$
(3)
From eq. $`(\text{2})`$, $`_iN^i=0`$. Thus, eq. $`(\text{3})`$ becomes

$$\frac{S}{t}=N_0,_iS=N_i,$$
(4)
neglecting terms of $`𝒪(N^2)`$. Solving eq. $`(\text{4})`$

$$S(X)=_0^{X^\mu }N_\mu 𝑑\stackrel{~}{X}^\mu ,$$
(5)
integrated along $`\gamma _1`$ or $`\gamma _2`$. The exponentiated form of eq. $`(\text{5})`$ is Yang’s nonintegrable phase factor Yang (1975), but with $`N_\mu `$ instead of the vector potential $`A_\mu `$. $`N_\mu `$ does pay a similar role in general relativity as the gauge field arising from local Galilean invariance in the nonrelativistic limit GLF (2002). Thus,

$$\mathrm{\Delta }\varphi =\frac{m}{\mathrm{}}_{\gamma _2}N_\mu 𝑑\stackrel{~}{X}^\mu \frac{m}{\mathrm{}}_{\gamma _1}N_\mu 𝑑\stackrel{~}{X}^\mu =\frac{m}{\mathrm{}}_𝒟R_{0i0j}(R(t))\stackrel{~}{X}^i𝑑\stackrel{~}{t}𝑑\stackrel{~}{X}^j$$
(6)
by Stokes’ theorem. Since the normal to $`𝒟`$ is a spacelike vector,

$$\mathrm{\Delta }\varphi \frac{m}{\mathrm{}}|R_{0i0j}|𝒜T$$
(7)
for $`R_{0i0j}`$ varying slowly in time. $`𝒜`$ is the area contained in the two intersecting orbits in Fig. 2, and $`T`$ is the transit time. Notice that while $`\mathrm{\Delta }\varphi `$ is proportional to $`R_{0i0j}`$ and is thus independent of the choice of frame in the nonrelativisitc limit, there is an unexpected dependence of $`\mathrm{\Delta }\varphi `$ of a test particle on its mass. For a cesium atom in a LEO of 480 km the phase difference is $`\mathrm{\Delta }\varphi _{Earth}(mGM_{}/R^3\mathrm{})𝒜T0.26𝒜T`$ rad cm-2 s-1.

The phase difference is proportional to the area inside the arms of the interferometer. Ultimately, this area is limited by the coherence time of the atoms. Using the parameters in PCC (2001) $`T=0.16`$ s and $`𝒜=0.01`$ cm2, we get $`\mathrm{\Delta }\varphi _{Earth}4\text{x}10^4`$ rad. A single-measurement fringe resolution $`0.074`$ rad has been demonstrated in Earth-based experiments PCC (2001), and after one-minute’s integration time, this resolution improves to $`0.011`$ rad. Measurement of $`\mathrm{\Delta }\varphi _{Earth}`$ requires the fringe sensitivity be increased by a factor of 27 beyond the above. However, $`\mathrm{\Delta }\varphi _{Earth}T^3`$ ($`𝒜(vT)^2`$ where $`v`$ is the velocity of the cesium atom after STIRAP), and increasing the coherence time $`T`$ from $`0.16`$ s to $`0.5`$ s will bring $`\mathrm{\Delta }\varphi _{Earth}`$ within the sensitivity of current interferometers.

That our proposed experiment is within the realm of feasibility is due to the advances in atom interferometry. These interferometers are so sensitive because beams of atoms are used instead of beams of photons. At 100 GeV/c2, the rest-mass of an atom is larger than the effective gravitational mass of a typical photon of 1 eV/c2 by eleven orders of magnitude. New interferometers based on Bose-Einstein Condensates (BECs) in development TSKSK (2000); SDECDHHRP (2000) can increase the sensitivity of atom-based interferometry significantly. In-situ phase measurements of BECs CHMW (1998); HMWC (1998); HDKTBEDJHRP (1999) offer the potential of even more sensitive measurements of curvature. Moreover, with long coherence times, vortices in BECs and superfluids may allow sensitive interference measurements of curvature.

Gravitational waves from astrophysical sources in the tens of kilohertz range would cause a phase shift four orders of magnitude smaller than $`\mathrm{\Delta }\varphi _{Earth}`$. This could in principle be detected with future improvements in atom interferometry.

###### Acknowledgements.

ADS and RYC were supported by a grant from the Office of Naval Research. We also thank Stephen A. Fulling for critically reading the manuscript.


text:
# A de Broglie-Bohm Like Model for Klein-Gordon Equation submited to Nuovo Cimento

## 1 Introduction

It is generally accepted that the Klein-Gordon equation has not a consistent de Broglie-Bohm like interpretation\[3, p. 498-503\]. However, in this paper we show that through the introduction of a hidden variable $`\mathrm{\Phi }`$ such that:

$$p_\mu =_\mu S_\mu \mathrm{\Phi },$$

where $`S`$ is the phase of the wave function, a solution to this problem can be easily found, in such way that it leads to the correct Schrödinger equation in the low-speed limit, and affords a sound interpretation of the quantum potential.

Within the theoretical framework of this work, the main assumption of the de Broglie-Bohm interpretation of quantum mechanics is not valid and a picture of particles moving under the action of electromagnetic field alone is recovered. However, the number of particles is not locally conserved. Also, given that a *hidden variable* has been introduced—on the grounds of some general electrodynamic considerations, presented in the first section—the representation of physical systems through wave functions wont be complete and, therefore, as foreseen by Einstein, Podolsky, and Rosen, quantum mechanics wont be a complete theory of motion.

## 2 On the Motion of Particles Under the Action of an Electromagnetic Field

Consider a non-stochastic ensemble of particles whose motion is described by means of a function

$$\stackrel{}{r}=\stackrel{}{r}(\stackrel{}{x},t),$$

such that $`\stackrel{}{r}(\stackrel{}{x},t)`$ represents the position, at time $`t`$, of the particle that passes through the point $`\stackrel{}{x}`$ at time *zero*—the $`\stackrel{}{x}`$-particle. This representation of motion coincides with the lagrangian representation used in hydrodynamics.

As to the function $`\stackrel{}{r}(\stackrel{}{x},t)`$, we suppose that it is invertible for any value of $`t`$. In other words, that there is a function $`\stackrel{}{x}=\stackrel{}{x}(\stackrel{}{r},t)`$, that gives us the coordinates, at time zero, of the particle that passes through the point $`\stackrel{}{r}`$ at time $`t`$. Also, we suppose that $`\stackrel{}{r}(\stackrel{}{x},t)`$ is a continuous function, altogether with its derivatives of as higher order as needed to secure the validity of our conclusions.

According to the definition of $`\stackrel{}{r}(\stackrel{}{x},t)`$, the velocity of the $`\stackrel{}{x}`$-particle at time $`t`$ is

$$\stackrel{}{u}(\stackrel{}{x},t)=\left(\frac{\stackrel{}{r}}{t}\right)_\stackrel{}{x}.$$

Writing $`\stackrel{}{x}`$ as a function of $`\stackrel{}{r}`$ and $`t`$, we get the velocity field at time $`t`$:

$$\stackrel{}{v}(\stackrel{}{r},t)=\stackrel{}{u}(\stackrel{}{x}(\stackrel{}{r},t),t).$$

As it is known from fluid kinematics, the corresponding field of acceleration is:

$$\frac{\stackrel{}{v}}{t}+(\stackrel{}{v}\stackrel{}{})\stackrel{}{v}.$$

From this, considering that

$$(\stackrel{}{v}\stackrel{}{})\stackrel{}{v}=\stackrel{}{}v^2/2\stackrel{}{v}\times (\stackrel{}{}\times \stackrel{}{v})$$

we can write the corresponding field of force:

$$\stackrel{}{f}=m\left(\stackrel{}{e}+\frac{1}{c}\stackrel{}{v}\times \stackrel{}{h}\right)$$

where

$$\stackrel{}{e}=\frac{1}{c}\frac{\stackrel{}{a}}{t}\stackrel{}{}a_0,$$
(1)
$$\stackrel{}{h}=\stackrel{}{}\times \stackrel{}{a},$$
(2)
$$a_0=c^2v^2/2,$$
(3)
$$\stackrel{}{a}=c\stackrel{}{v},$$
(4)
and $`c`$ is a constant, with dimensions of speed, that we shall take equal to the velocity of light.

Equations (1) and (2) are formally analogous to the definition of the electric and magnetic fields from the electrodynamic potentials. From them we can prove that

$$\stackrel{}{}\stackrel{}{h}=0,$$

and

$$\stackrel{}{}\times \stackrel{}{e}=\frac{1}{c}\frac{\stackrel{}{h}}{t},$$

analogous to the first pair of Maxwell equations. Furthermore, the right members of equations (3) and (4), but for a constant factor, correspond to the low-speed approximation of the components of the relativistic four-velocity, which suggest us to investigate the properties of this field, defined as:

$$v^\mu =(\frac{c}{\sqrt{1v^2/c^2}},\frac{\stackrel{}{v}}{\sqrt{1v^2/c^2}})$$

Using four-dimensional tensorial notation, the derivative of $`v_\mu `$ with respect to the proper time, along the world-line of the corresponding particle is:

$$\frac{dv_\mu }{ds}=v^\nu _\nu v_\mu $$

From the identity

$$v^\nu v_\nu =c^2$$
(5)
we prove that

$$v^\nu _\mu v_\nu =0,$$

that allows us to write

$$\frac{dv_\mu }{ds}=(_\nu v_\mu _\mu v_\nu )v_\nu .$$
(6)
If the four-force acting on the ensemble of particles has an electromagnetic origin, it is given by the expression

$$f_\mu =\frac{q}{c}F_{\mu \nu }v^\nu ,$$
(7)
where

$$F_{\mu \nu }=_\mu A_\nu _\nu A_\mu $$
(8)
is Faraday’s tensor, and $`A_\mu `$ is the electrodynamic four-potential.

From the law of motion

$$m\frac{dv_\mu }{ds}=f_\mu ,$$

and equations (6) to (8) we prove that

$$(_\mu p_\nu _\nu p_\mu )v^\nu =0,$$
(9)
where

$$p_\mu =mv_\mu +\frac{q}{c}A_\mu .$$
(10)
There is a class of solutions of equations (9) where $`p_\mu `$ is the four-gradient of a function of space-time coordinates

$$p_\mu =_\mu \varphi .$$
(11)
Therefore,

$$_\mu p_\nu _\nu p_\mu =0,$$

and (9) is obviously satisfied.

Equation (10) can be written in the form

$$mv_\mu =_\mu \varphi \frac{q}{c}A_\mu $$
(12)
If the four-potential meets the Lorentz condition,

$$^\mu A_\mu =0,$$
(13)
equation (12) is analogous to the decomposition of a classical, three-dimensional field, into an irrotational and a solenoidal part. The electromagnetic field appears thus as determining the four-dimensional vorticity of the field of *kinetic momentum* $`mv_\mu `$.

Equation (5) tells us that functions $`\varphi `$ are not arbitrary, but are subject to the condition:

$$(^\mu \varphi +\frac{q}{c}A^\mu )(_\mu \varphi +\frac{q}{c}A_\mu )=m^2c^2,$$
(14)
which is the relativistic Hamilton-Jacobi equation.\[2, Ch. VIII\]

The components of the kinetic momentum are:

$$(\frac{K}{c},\stackrel{}{p})=(\frac{1}{c}\frac{\varphi }{t}\frac{q}{c}V,\stackrel{}{}\varphi \frac{q}{c}\stackrel{}{A}),$$
(15)
where $`V`$ and $`\stackrel{}{A}`$ are the components of the electrodynamic four-potential.

In the low-speed limit, we have

$$Kmc^2+\frac{p^2}{2m}.$$

From this and (15) we get:

$$\frac{\varphi }{t}+\frac{(\stackrel{}{}\varphi \frac{q}{c}\stackrel{}{A})^2}{2m}+qV+mc^2=0,$$
(16)
that, but for the constant term $`mc^2`$, is the non-relativistic Hamilton-Jacobi equation.

Notice that

$$mv_\mu 𝑑x^\mu =\frac{q}{c}A_\mu 𝑑x^\mu ,$$

for any closed path in space time. In particular

$$\stackrel{}{p}d^3\stackrel{}{x}=\frac{q}{c}\stackrel{}{A}d^3\stackrel{}{x},$$

for any closed path of simultaneous events in three-dimensional space.

## 3 A de Broglie-Bohm Like Interpretation of Klein-Gordon Equation

Consider the Klein-Gordon equation describing an ensemble of spin-less particles with mass $`m`$ and charge $`q`$:

$$(i\mathrm{}^\mu \frac{q}{c}A^\mu )(i\mathrm{}_\mu \frac{q}{c}A_\mu )\mathrm{\Psi }=0,$$

where $`A^\mu `$ satisfies the Lorentz condition (13). In other words:

$$\mathrm{}^2^\mu _\mu \mathrm{\Psi }\frac{2i\mathrm{}}{c}A^\mu _\mu \mathrm{\Psi }+\frac{q^2}{c^2}A^\mu A_\mu \mathrm{\Psi }=m^2c^2\mathrm{\Psi }.$$
(17)
We can easily show that

$$^\mu \left(\frac{i\mathrm{}}{2}(\mathrm{\Psi }^{}_\mu \mathrm{\Psi }\mathrm{\Psi }_\mu \mathrm{\Psi }^{})\frac{q}{c}A_\mu \mathrm{\Psi }^{}\mathrm{\Psi }\right)=0.$$
(18)
By the Madelung substitution

$$\mathrm{\Psi }=\sqrt{\rho }e^{\frac{i}{\mathrm{}}S},$$
(19)
equation (20) is transformed into

$$^\mu (\rho (_\mu S\frac{e}{c}A_\mu ))=0,$$
(20)
that, in view of (12), tells us that

$$\rho =\mathrm{\Psi }^{}\mathrm{\Psi },$$
(21)
and

$$w_\mu =\frac{_\mu S}{m}\frac{q}{mc}A_\mu ,$$
(22)
might be interpreted as the density of particles in the (local) system of reference, where the particles are at rest, and the field of four-velocity, respectively. (Hence $`p_\mu =_\mu S`$.) This interpretation is not sound, however, because $`w_\mu `$ is not an implicitly time-like four-vector. Actually, from (17) and (19), it can be shown that

$$(^\mu S\frac{q}{c}A^\mu )(_\mu S\frac{q}{c}A_\mu )=m^2c^2+\mathrm{}^2\frac{^\mu _\mu \sqrt{\rho }}{\sqrt{\rho }},$$
(23)
in plain disagreement with (14).

However, considering that $`w_\mu `$ has already the required four-vorticity, we’ll suppose that there is a function $`\mathrm{\Phi }`$ such that

$$mv_\mu =_\mu S_\mu \mathrm{\Phi }\frac{q}{c}A_\mu .$$
(24)
From (14) and (23) we can show that $`\mathrm{\Phi }`$ has to satisfy the condition

$$2mv^\mu _\mu \mathrm{\Phi }+^\mu \mathrm{\Phi }_\mu \mathrm{\Phi }=\mathrm{}^2\frac{^\mu _\mu \sqrt{\rho }}{\sqrt{\rho }}.$$
(25)
We’ll show that this is a sound hypothesis by proving that it leads to the correct equations in the low-speed limit, where we suppose that

$$\frac{1}{c^2}\left(\frac{\mathrm{\Phi }}{t}\right)^2=0,\frac{\mathrm{}^2}{c^2\sqrt{\rho }}\frac{^2\sqrt{\rho }}{t^2}=0,$$

in such way that (25) can be rewritten, using classical, three-dimensional, notation, as

$$2m\frac{\mathrm{\Phi }}{t}+2m\stackrel{}{v}\stackrel{}{}\mathrm{\Phi }(\stackrel{}{}\mathrm{\Phi })^2=\mathrm{}^2\frac{\mathrm{\Delta }\sqrt{\rho }}{\sqrt{\rho }}$$
(26)
or,

$$\frac{\mathrm{\Phi }}{t}+(\stackrel{}{v}\stackrel{}{}\mathrm{\Phi })\stackrel{}{}\mathrm{\Phi }+\frac{(\stackrel{}{}\mathrm{\Phi })^2}{2m}=\frac{\mathrm{}^2}{2m}\frac{\mathrm{\Delta }\sqrt{\rho }}{\sqrt{\rho }}$$
(27)
In other words:

$$\frac{(S\frac{q}{c}\stackrel{}{A})}{m}\mathrm{\Phi }+\frac{(\stackrel{}{}\mathrm{\Phi })^2}{2m}=\frac{\mathrm{}^2}{2m}\frac{\mathrm{\Delta }\sqrt{\rho }}{\sqrt{\rho }}\frac{\mathrm{\Phi }}{t}$$
(28)
From this we can build a low-speed limit approximation of the kinetic energy:

$$Kmc^2+\frac{(S+\mathrm{\Phi }\frac{q}{c}\stackrel{}{A})^2}{2m}=mc^2+\frac{(S\frac{q}{c}\stackrel{}{A})^2}{2m}\frac{\mathrm{}^2}{2m}\frac{\mathrm{\Delta }\sqrt{\rho }}{\sqrt{\rho }}\frac{\mathrm{\Phi }}{t},$$
(29)
Finally, given that in this case,

$$K=\frac{(S+\mathrm{\Phi })}{t}qV.$$
(30)
we conclude that, in the low-speed limit:

$$\frac{S}{t}+mc^2+\frac{(S\frac{q}{c}\stackrel{}{A})^2}{2m}\frac{\mathrm{}^2}{2m}\frac{\mathrm{\Delta }\sqrt{\rho }}{\sqrt{\rho }}+qV=0,$$
(31)
that, but for the constant term $`mc^2`$, is the classical Hamilton-Jacobi equation with quantum potential. From this and (20) the corresponding Schrödinger equation is easily recovered.

In the general case, we get a picture of classical particles moving under the action of the electromagnetic field alone: there is not *quantum potential*. The field of four-velocity is given by (24), where $`S`$ is the phase of the wave function and $`\mathrm{\Phi }`$ is a *hiden variable*. $`\rho `$, $`S`$ and $`\mathrm{\Phi }`$ are solutions of a system of non-linear differential equations: (20), (23), and

$$(^\mu S^\mu \mathrm{\Phi }\frac{q}{c}A^\mu )(_\mu S_\mu \mathrm{\Phi }\frac{q}{c}A_\mu )=m^2c^2,$$
(32)
that, however, does not determine them completely. But for the restrictions imposed in the low-speed limit, the function $`\mathrm{\Phi }`$ could be equal to $`\varphi S`$ for any solution of equation (14).

The number of particles is not locally conserved. From (20) and (24) we get:

$$^\mu (\rho v_\mu )=\frac{^\mu (\rho _\mu \mathrm{\Phi })}{m}.$$
(33)
The $`\mathrm{\Phi }`$ function appears thus related to local creation/annihilation processes, which leads to the appearance of a *quantum force* in the corresponding hydrodynamic model.

Let

$$T_{\mu \nu }=m\rho v_\mu v_\nu $$
(34)
Then

$$^\nu T_{\mu \nu }=m\rho v_\nu ^\nu v_\mu +mv_\mu ^\nu (\rho v_\nu )=\frac{q}{c}\rho v^\nu F_{\mu \nu }v_\mu ^\nu (\rho _\nu \mathrm{\Phi })$$
(35)
Equations (34) and (35) describe a pressure-less fluid of charged particles under the action of the electromagnetic fluid and a force:

$$K_\mu =v_\mu ^\nu (\rho _\nu \mathrm{\Phi })$$
(36)
associated—through (33)—to local creation/annihilation processes.

## 4 The Low-Speed Limit Revisited

In the low speed limit, the equations of quantum mechanics for spin-less particles can be proved to be equivalent to:

$$m\left(\frac{\stackrel{}{u}}{t}+(\stackrel{}{u})\stackrel{}{u}\right)=q\text{E}+\frac{q}{c}\stackrel{}{u}\times \text{H}\stackrel{}{}Q$$
(37)
and

$$\frac{\rho }{t}+\stackrel{}{}(\rho \stackrel{}{u})=0,$$
(38)
where

$$\rho =\mathrm{\Psi }^{}\mathrm{\Psi },$$
(39)
$$m\stackrel{}{u}=\stackrel{}{}S\frac{q}{c}\stackrel{}{A},$$
(40)
and

$$Q=\frac{\mathrm{}^2}{2m}\frac{\mathrm{\Delta }\sqrt{\rho }}{\sqrt{\rho }},$$
(41)
is the *quantum potential*.

From (26)

$$\stackrel{}{}Q=\frac{\stackrel{}{}\mathrm{\Phi }}{t}(\stackrel{}{v}\stackrel{}{})\stackrel{}{}\varphi (\stackrel{}{}\varphi \stackrel{}{})\stackrel{}{v}\stackrel{}{}\varphi \times (\stackrel{}{}\times \stackrel{}{v})+\stackrel{}{}\frac{(\stackrel{}{}\mathrm{\Phi })^2}{2m}$$

On this basis, considering that in this limit the field of velocity is given by:

$$m\stackrel{}{v}=\stackrel{}{}S+\stackrel{}{}\mathrm{\Phi }\frac{q}{c}\stackrel{}{A},$$
(42)
and that, therefore,

$$m\stackrel{}{v}=m\stackrel{}{u}+\stackrel{}{}\mathrm{\Phi },$$

and

$$m\stackrel{}{}\times \stackrel{}{v}=\frac{q}{c}\stackrel{}{H}$$

equation (37) can be rewritten as:

$$m\left(\frac{\stackrel{}{v}}{t}+(\stackrel{}{v})\stackrel{}{v}\right)=q\text{E}+\frac{q}{c}\stackrel{}{v}\times \text{H}.$$
(43)
The continuity equation is transformed into:

$$\frac{\rho }{t}+\stackrel{}{}(\rho \stackrel{}{v})=\frac{\stackrel{}{}(\rho \stackrel{}{}\mathrm{\Phi })}{m},$$
(44)
Therefore, the number of particles is not locally conserved.

Equations (43) and (44) describe a flux of classical particles, under the action of the electromagnetic field.

## 5 Conclusions and Remarks

If the ideas exposed in this paper are proven to be valid:

1. The main assumption of the de Broglie-Bohm theory—that the local impulse of quantum particles is given by the gradient of the phase of the wave function—wont be accurate. \[See eq. (24).\]
2. However, there will be still room for a classical interpretation of quantum phenomena, in terms of particles moving along well defined trajectories, under the action of the electromagnetic field.
3. The number of particles wont be locally conserved.
4. Given that a hidden variable has been introduced—after some considerations on electrodynamics—the representation of physical systems through wave functions wont be complete and, therefore, as foreseen by Einstein, Podolsky, and Rosen, quantum mechanics wont be a complete theory of motion. Actually, but for the restrictions imposed in the low-speed limit, the function $`\mathrm{\Phi }`$ could be equal to $`\varphi S`$ for any solution of equation (14).