text:
# QUARK-GLUON PLASMA IN THE EARLY UNIVERSE

## 1 Introduction

Who has not wondered what happens when matter gets squeezed to higher and higher densities or heated to higher and higher temperatures? At first molecules would dissociate into atoms, then atoms would ionize, then atomic nuclei would undergo a liquid-gas type phase transition to a vapor of protons and neutrons. At high baryon density there will be so many protons and neutrons per unit volume that their constituent quarks will no longer know to which nucleon they belong. This is expected to happen when nucleons are about a factor of two closer to each other than they are in atomic nuclei, that is, at a density of about 8-10 times that in a nucleus. (Normal nuclear density is 0.15 nucleons per cubic fermi - equivalent to $`2\times 10^{14}`$ g/cm3.) At this point quarks will be free to move over distances much greater than the confinement distance of 1 fm. This type of matter is a cold quark gas. It is possible, but not yet known, that it exists in the cores of neutron stars.

At high temperature there will be many particle-antiparticle pairs produced. The lightest hadrons, the pions, will be the first to be produced in abundance, quickly followed by all the hadrons listed in the Particle Data Tables . Hagedorn realized already in the 1960’s that the level density of hadrons grows exponentially $`dN/dm\mathrm{exp}(m/T_0)`$ where $`m`$ is the hadron mass . The Hagedorn temperature $`T_0`$ has the numerical value 160 MeV. A thermal ensemble of hadrons treated as noninteracting particles therefore has a limiting temperature equal to $`T_0`$ because of this exponential growth. In a sense it is the end of the hadronic world; higher temperatures bring a quark-gluon plasma. Intuitively, this transition comes about because the hadron density is so great that the constituent quarks and gluons are confused about which hadron they ought to be confined to.

The main interest in a cosmological hadron to quark-gluon plasma phase transition arises from its potential to influence big bang nucleosynthesis, discussed in section 2. Whether or not QCD with its known set of parameters undergoes a first order transition or something smoother is still not settled. This issue is discussed in section 3 using perturbation theory and lattice gauge theory. Assuming that there is a first order phase transition one needs nucleation theory to understand how the transition proceeds, and this topic is discussed in section 4. The cosmological QCD phase transition is then analyzed in section 5. Conclusions, and paths for future work, are given in section 6.

## 2 Inhomogeneous Big Bang Nucleosynthesis

A cosmological first order phase transition at $`T160180`$ MeV could leave spatial inhomogeneities in the baryon to entropy ratio and in the ratio of protons to neutrons. If these inhomogeneities survive to $`T0.11`$ MeV they could influence nucleosynthesis. This was first pointed out and analyzed by Witten , by Applegate, Hogan, and Scherer , and by Alcock, Fuller, and Mathews . In thermal and chemical equilibrium one might expect that the baryon density in the quark-gluon phase is higher than in the hadron phase. This is called the baryon density contrast. Assuming a critical temperature of $`160<T_c<180`$ MeV, Olive and I computed this baryon density contrast to be 1.5 to 2.5 when hadronic interactions were neglected and 5 to 7 when they were included . One would expect the last regions of space to undergo the phase conversion to contain more baryons per unit volume than the first regions to phase convert because of the lack of time for baryons to diffuse. After phase completion the neutrons will diffuse more rapidly than protons because they are electrically neutral and therefore do not Coulomb scatter on electrons. This leads to isospin inhomegeneities, at least temporarily.

A detailed calculation of inhomogeneous nucleosynthesis with a comparison to observed abundances of the light elements was done by Kurki-Suonio, Matzner, Olive, and Schramm . They considered baryon density contrasts ranging from 1 to 100 and a fraction of matter in the high density regions ranging from 1/64 to 1/4. The average separation of the high density regions $`l`$ was left as a free parameter as was the average baryon to photon ratio of the universe. Differential diffusion of protons and neutrons was accounted for and then a standard nucleosynthesis code was run. By fitting the observed abundances of 4He, D, 3He, and 7Li they concluded that the baryon to photon ratio must lie between $`2\times 10^{10}`$ and $`7\times 10^{10}`$ (or $`2\times 10^{10}`$ if certain constraints on 7Li were relaxed). They also concluded that $`l<150`$ m at the time of nucleosynthesis. At the completion of the QCD phase transition this upper limit would have been about 1 m. Early estimates of this distance were uncertain but promising . See also the reviews listed in . A better estimate of this scale is the purpose of sections 4 and 5.

Recently the inhomogeneous nucleosynthesis calculation was redone with technical improvements and updated estimates of the cosmic abundances of the relevant light elements by Kainulainen, Kurki-Suonio, and Sihvola . Their results are shown in figure 1.

The high density matter was distributed in spheres. The inhomogeneities are ineffective in influencing nucleosynthesis unless the high density regions were separated by more than about 150 m at $`T=1`$ MeV.

## 3 QCD Equation of State

The QCD equation of state can, in principle, be computed from the functional integral formula for the partition function .

$$Z=[d\overline{\psi }][d\psi ][dA_\mu ^a][d\overline{C}_a][dC_a]\mathrm{exp}\left\{_0^\beta 𝑑\tau _Vd^3x_{eff}\right\}$$
(1)
Here $`\psi `$ is a quark field, $`A_\mu ^a`$ is the gluon field, and $`C_a`$ is the ghost field. The Lagrangian is

$`_{eff}`$ $`=`$ $`{\displaystyle \frac{1}{\rho }}\left(^\mu A_\mu ^a\right)^2+\mu _f\overline{\psi }_f\gamma ^0\psi _f+_\mu \overline{C}_a^\mu C_a+gf_{abc}\overline{C}_a_\mu A_b^\mu C_c,`$ (2)
$``$ $`=`$ $`{\displaystyle \frac{1}{4}}F_a^{\mu \nu }F_{\mu \nu }^a\overline{\psi }_f\left(i\overline{)}m_fg\overline{)}A^aG^a\right)\psi _f,`$ (3)
$`F_a^{\mu \nu }`$ $`=`$ $`^\mu A_a^\nu ^\nu A_a^\mu gf_{abc}A_b^\mu A_c^\nu .`$ (4)
The subscript $`f`$ refers to quark flavor. The $`G^a`$ are the generators of the color group $`SU(N_c)`$. The quark fields are Grassmann variables antiperiodic in imaginary time $`\tau =it`$, $`0<\tau <\beta =1/T`$, the gluon fields are periodic in this time, and the ghost fields are Grassmann variables periodic in imaginary time. (The complex ghost field may be thought of as subtracting 2 of the 4 degrees of freedom implied by $`A^\mu `$ to yield the 2 physical transverse degrees of freedom of the massless gluon.) The Lagrangian is written in covariant gauge with $`\rho `$ the gauge parameter. Physical observables are independent of $`\rho `$.

There are two well-known ways to evaluate this functional integral. One is perturbation theory in the coupling constant $`g`$; the other is numerical estimation using Monte Carlo techniques on a discrete lattice. Before reviewing results obtained in these approaches, let us try to predict what the result might be. Suppose that all quarks were massless and that all chemical potentials were zero. Then there is only one parameter in the QCD Lagrangian, the coupling $`g`$ which is dimensionless. If this coupling is really a constant then dimensional analysis implies that the pressure must be some function of g times the temperature to the fourth power: $`P=f(g)T^4`$. Then the entropy density is $`s=4f(g)T^3`$, and the equation of state may be written as $`ϵ=3P`$, the same as a free gas of massless particles. This is a trivial equation of state without a phase transition. In contrast the quantum version of the theory displays dimensional transmutation. The dimensionless coupling constant $`g`$ gets replaced by a renormalization group running coupling. At finite temperature and at one-loop order

$$g^2\overline{g}^2=\frac{24\pi ^2}{(11N2N_f)\mathrm{ln}\left(T/\mathrm{\Lambda }_{QCD}\right)}$$
(5)
for an $`SU(N_c)`$ gauge group ($`N_c=3`$ for QCD) and $`N_f`$ flavors of quark. The scale parameter $`\mathrm{\Lambda }_{QCD}`$, which has dimensions of energy, enters as a result of the necessity to regulate and renormalize, that is, to render correlation functions finite at some arbitrarily chosen energy or momentum point. As a result the equation of state with $`P=f(\overline{g}(T))T^4`$ can now become complicated and not predictable by dimensional analysis alone. This allows nontrivial thermodynamics and the possibility of a phase transition. This is on account of quantum mechanics, not classical interacting field theory.

The partition function for QCD was calculated within perturbation theory to order $`g^3`$ in the late 1970’s by myself and then to order $`g^5`$ in the mid 1990’s by Arnold and Zhai , by Kastening and Zhai , and by Braaten and Nieto . For $`N_f`$ flavors of massless quarks and zero chemical potentials the pressure is

$`P`$ $`=`$ $`{\displaystyle \frac{8\pi ^2}{45}}T^4[1+{\displaystyle \frac{21}{32}}N_f{\displaystyle \frac{15}{4}}(1+{\displaystyle \frac{5}{12}}N_f){\displaystyle \frac{\alpha _s(M)}{\pi }}`$ (6)
$`+`$ $`30\left(1+{\displaystyle \frac{1}{6}}N_f\right)^{3/2}\left({\displaystyle \frac{\alpha _s(M)}{\pi }}\right)^{3/2}+F_4(N_f,\mathrm{ln}(M/2\pi T))\left({\displaystyle \frac{\alpha _s(M)}{\pi }}\right)^2`$
$`+`$ $`F_5(N_f,\mathrm{ln}(M/2\pi T))\left({\displaystyle \frac{\alpha _s(M)}{\pi }}\right)^{5/2}+𝒪\left(\alpha _s^3\mathrm{ln}\alpha _s\right)].`$
Here $`\alpha _s=g^2/4\pi `$, and $`F_4`$ and $`F_5`$ are lengthy untransparent functions of the number of flavors and of the renormalization scale $`M`$. The latter should be chosen judiciously so as to minimize the contribution of higher order terms in the expansion. For sake of illustration consider $`N_f=3`$, $`M=2\pi T`$, $`T=500`$ MeV, and $`\alpha _s0.2`$. Then

$$P=\frac{19\pi ^2}{36}T^4\left[10.18+0.30+0.060.37\mathrm{}\right]$$
(7)
representing the contributions of $`\alpha _s`$ to the power of 0, 1, 3/2, 2, and 5/2. At so low a temperature, only about three times the expected $`T_c`$, this expansion is apparently not very quickly converging. Although very useful for very high temperature because asymptotic freedom forces the running coupling to zero logarithmically, perturbation theory cannot reliably predict the location or order of a possible phase transition.

By its very nature perturbation theory using free particle basis states is incapable of describing a phase transition. An alternative is to discretize space and (imaginary) time and estimate the functional integral expression for $`Z`$ numerically with Monte Carlo techniques . One works with a spatial volume $`V=(N_sa)^3`$ and a time dimension $`\beta =N_\tau a`$, where $`a`$ is the lattice spacing and $`N_s`$ and $`N_\tau `$ are the number of sites in each direction. The thermodynamic limit is $`a0`$, $`N_s\mathrm{}`$, $`N_\tau \mathrm{}`$, subject to $`\beta `$ held fixed and $`V\mathrm{}`$. The difficulties associated in reaching for the thermodynamic limit include the enormous computer power and storage necessary for the huge number of degrees of freedom inherent in a quantum field theory, and the theoretical difficulty with describing fermions on a lattice. A typical lattice for pure gauge theory with no quarks is $`64^3\times 8`$, and with quarks is $`16^3\times 4`$. These are not very large lattices, especially in the latter case. Nevertheless computer power and storage is increasing with time, and practitioners are constantly improving the numerical algorithms. For example, improved lattice actions, where derivatives are better approximated via higher order finite differences and short distance behavior can be supplemented with perturbation theory on account of asymptotic freedom, have helped significantly. Here I have relied on a recent review by Karsch .

Pure quarkless $`SU(2)`$ gauge theory is expected to undergo a second order phase transition reflecting color deconfinement at high temperature whereas $`SU(N_c)`$ with $`N_c>2`$ is expected to undergo a first order phase transition. The $`SU(3)`$ equation of state computed in lattice gauge theory is shown in figure 2.

One physical quantity is required to set the overall energy scale. This is chosen to be the string tension $`\sigma `$ at zero temperature; its numerical value is taken to be the same as in our physical world. Then the critical temperature is about 270 MeV and the latent heat is about 0.5 GeV/fm3. Below $`T_c`$ is a gas of (presumably) weakly interacting very massive glueballs. The lightest one is 1700 MeV. Above $`T_c`$ all 16 degrees of freedom of massless (interacting) gluons are liberated. These computational results seem to be stable with time and therefore well established and accepted.

A compendium of results from different groups for QCD with two degenerate flavors of quarks is shown in figure 3.

In the limit that the quarks are massless a second order transition is expected, with no proper thermodynamic transition when they are massive. Even with massive quarks one sees a sharp rise in the energy density in a narrow window of temperature. The middle of this range can be used to define a crossover temperature still called, somewhat incorrectly, $`T_c`$. In order to set the physical energy scale the mass of the lightest vector meson, the $`\rho `$-meson, is calculated and set equal to its real world value of 770 MeV. The mass of the lightest pseudoscalar meson, the $`\pi `$-meson, is calculated and compared to the real world value of 140 MeV. The pseudoscalar mass is a monotonically increasing function of the quark mass so it can be used to determine what quark mass ought to be used in the lattice calculation to represent as closely as possible the real world. The $`T_c`$ is plotted against the ratio of the pseudoscalar to vector meson masses, in units of the vector mass. Taking the limit of zero quark mass is notoriously difficult. The actual range of computed pion masses starts at 1.2 GeV and decreases to 300 MeV but no lower. There is a lot of scatter among the results from different groups using different methods but the trend is towards $`T_c=170`$ MeV in the zero quark mass limit.

QCD with three flavors of massless quarks is expected to undergo a first order phase transition. The question is: What is the order of the transition for the real world with up and down quark masses of order 5-7 MeV and a strange quark mass of order 110-150 MeV? Some results are shown in figure 4.

The lattice is $`16^3\times 4`$ in size. The quark masses are actually temperature-dependent (the usual situation for numerical reasons) with $`m_{\mathrm{light}}=0.4T`$ and $`m_{\mathrm{heavy}}=T`$. Results for two light flavors, two light and one heavy flavor, and three light flavors are shown in this figure. The energy density $`ϵ`$, in units of $`T^4`$, is plotted against $`T/T_c`$. In all three cases there is a big jump in the energy density over a narrow window of temperature. The 2+1 flavor case is almost indistinguishable from the 3 flavor case in the critical region. The statistical error bars are very small. More significant are the systematic errors that are much more difficult to estimate. Due to the relatively small lattice size and still rather large light quark mass it is difficult to say whether or not the real world has a first order phase transition. Even if not, the large jump in energy density may make the real situation practically the same as if it did.

## 4 Relativistic Nucleation Theory

The dynamics of first-order phase transitions has fascinated scientists at least since the time of Maxwell and Van der Waals. Much work on the classical theory of nucleation of gases and liquids was carried out in the early part of the 20th century. There were and still are many important applications, such as cloud and bubble chambers, freezing of liquids, and precipitation in the atmosphere. The goal of nucleation theory is to compute the probability that a bubble or droplet of one phase appears in a system initially in the other phase near the critical temperature. Homogeneous nucleation theory applies when the system is pure; inhomogeneous nucleation theory applies when impurities cause the formation of bubbles or droplets. For the application we have in mind, namely the early universe, it seems that homogeneous nucleation theory is appropriate. In the everyday world it is usually the opposite; dust or ions in the atmosphere are much more efficient in producing precipitation. Nucleation theory is applicable for first-order phase transitions when the matter is not dramatically supercooled or superheated. If substantial supercooling or superheating is present, or if the phase transition is second-order, then the relevant dynamics is spinodal decomposition.

Suppose that a system is cooled below its critical temperature. Then there exists a critical sized droplet (or bubble, depending on whether the energy density in the lower temperature phase is greater or less than the higher temperature phase). If a droplet that forms because of statistical fluctuations is too small, its surface free energy is relatively large and the cost in total free energy is positive. The droplet will evaporate. If the droplet is large, its surface free energy is unimportant, and the droplet will accrete molecules and grow. A droplet of critical size is metastable, it is balanced between evaporation and accretion. The classical theory of Becker and Döring , which is nicely reviewed by McDonald , says that the probability per unit time per unit volume to nucleate the dense liquid phase from a dilute gas is given by

$$I=a(i_{})\left(\frac{ϵ^{\prime \prime }(i_{})}{2\pi T}\right)^{1/2}\overline{n}(1)e^{ϵ(i_{})/T},$$
(8)
where $`ϵ(i_{})`$ is the formation energy of a critical sized droplet consisting of $`i_{}`$ molecules, prime denotes differentiation with respect to the number of molecules $`i`$, $`\overline{n}(1)`$ is the density of single molecules and $`a(i_{})`$ is the accretion rate of single molecules on a critical droplet. Usually the accretion rate is taken to be

$$a(i_{})=\frac{1}{2}\overline{n}(1)\overline{v}4\pi r_{}^2s,$$
(9)
which is the flux of particles ($`\overline{v}`$ is the mean speed of gas molecules) striking the surface of the critical droplet times a ‘sticking fraction’ $`s`$ less than one. The first term in this formula is a dynamical factor influencing the growth rate, the second term characterizes fluctuations about the critical droplet, and the product of the third and fourth terms gives the quasi-equilibrium number density of critical sized droplets.

The modern theory of nucleation was pioneered by Langer . Langer’s theory is based in a more fundamental way on the microscopic interactions of atoms and molecules. It can also be applied close to a critical point where in fact most of the current interest in the condensed matter community has been. Langer’s modern theory of nucleation yields the following formula for the rate:

$$I=\frac{\kappa }{2\pi }\mathrm{\Omega }_0e^{\mathrm{\Delta }F_{}/T}$$
(10)
where $`\mathrm{\Delta }F_{}`$ is the change in the free energy of the system due to the formation of the critical droplet. $`\mathrm{\Omega }_0`$ is a statistical prefactor which measures the available phase space volume.

$$\mathrm{\Omega }_0=2\left(\frac{\sigma }{3T}\right)^{3/2}\left(\frac{r_{}}{\xi _g}\right)^4$$
(11)
Here $`\xi _g`$ is the correlation length in the gaseous phase. $`\kappa `$ is a dynamical prefactor which determines the exponential growth rate of critical droplets which are perturbed from their quasi-equilibrium radius $`r_{}`$.

$$\kappa =\frac{d}{dt}\mathrm{ln}[r(t)r_{}]$$
(12)
The basic structure is the same as in the classical theory, but the prefactors are different. The dynamical prefactor has been calculated by Langer and Turski  and by Kawasaki  for a liquid-gas phase transition near the critical point, where the gas is not dilute, to be

$$\kappa =\frac{2\lambda \sigma T}{\mathrm{}^2n_{\mathrm{}}^2r_{}^3}.$$
(13)
This involves the thermal conductivity $`\lambda `$, the surface free energy $`\sigma `$, the latent heat per molecule $`\mathrm{}`$ and the density of molecules in the liquid phase $`n_{\mathrm{}}`$. The interesting physics in this expression is the appearance of the thermal conductivity. In order for the droplet to grow beyond the critical size, latent heat must be conducted away from the surface into the gas.

For a relativistic system of particles or quantum fields which has no net conserved charge, such as baryon number, the thermal conductivity vanishes. The reason is that there is no rest frame defined by the baryon density to refer to heat transport. Hence the Langer formula cannot be applied to such systems. The extension to such systems was made by Csernai and me . The change in free energy due to the appearance of a bubble of hadronic matter in quark-gluon plasma is

$$\mathrm{\Delta }F=\frac{4\pi }{3}r^3\left[P_\mathrm{q}(T)P_\mathrm{h}(T)\right]+4\pi r^2\sigma ,$$
(14)
where $`r`$ is the radius. The critical sized bubble has radius

$$r_{}=\frac{2\sigma }{P_\mathrm{h}(T)P_\mathrm{q}(T)}$$
(15)
which leads to

$$\mathrm{\Delta }F_{}=\frac{4\pi }{3}\sigma r_{}^2.$$
(16)
The expression for the statistical prefactor is unchanged, but the formula for the dynamical prefactor, or growth rate, is proportional to the shear $`\eta _q`$ and bulk $`\zeta _q`$ viscosities in the quark-gluon plasma instead of the thermal conductivity.

$$\kappa =\frac{4\sigma (3\zeta _q+4\eta _q)}{3(\mathrm{\Delta }w)^2r_{}^3}$$
(17)
This is inversely proportional to the square of the enthalpy ($`w=ϵ+P)`$ difference between the two phases. The nucleation rate is

$$I=\frac{4}{\pi }\left(\frac{\sigma }{3T}\right)^{3/2}\frac{\sigma (3\zeta _q+4\eta _q)r_{}}{3(\mathrm{\Delta }w)^2\xi _q^4}e^{\mathrm{\Delta }F_{}/T}.$$
(18)
For numerical purposes we use an MIT bag model type equation of state with

$`P_q`$ $`=`$ $`\left(45.5+14.25\right){\displaystyle \frac{\pi ^2}{90}}T^4B`$
$`P_h`$ $`=`$ $`\left(5.5+14.25\right){\displaystyle \frac{\pi ^2}{90}}T^4.`$ (19)
The 45.5 approximates the effective number of degrees of freedom arising from massless gluons and up and down quarks and a strange quark mass comparable to the temperature. The 5.5 approximates the hadronic equation of state near $`T_c`$ arising from a multitude of massive hadrons. The 14.25 arises from photons, neutrinos, electrons, and muons common to both phases. The bag constant $`B`$ is chosen to give $`T_c=160`$ MeV. Furthermore Csernai and I  estimated that $`\sigma =50`$ MeV/fm3, $`\xi _q=0.7`$ fm, and $`\eta _q=18T^3`$ . Generally for relativistic systems the bulk viscosity is small in comparison to the shear viscosity and so we neglect it.

## 5 Cosmological QCD Phase Transition

Given the nucleation rate one would like to know the (volume) fraction of space $`h(t)`$ which has been converted from the quark-gluon plasma to the hadronic gas at the proper time $`t`$ in the early universe. This requires kinetic equations that use the nucleation rate $`I`$ as an input. Here we use a rate equation first proposed by Csernai and me . The nucleation rate $`I`$ is the probability to form a bubble of critical size per unit time per unit volume. If the system cools to $`T_c`$ at time $`t_c`$ then at some later time $`t`$ the fraction of space which has been converted to the hadronic phase is

$$h(t)=_{t_c}^t𝑑t^{}I(T(t^{}))[1h(t^{})]V(t^{},t).$$
(20)
$`V(t^{},t)`$ is the volume of a hadronic bubble at time $`t`$ which was nucleated at the earlier time $`t^{}`$; this takes into account bubble growth. The factor $`1h(t^{})`$ takes into account the fact that new bubbles can only be nucleated in the fraction of space not already occupied by the hadronic gas. This conservative approach neglects any spatial variation in the temperature. However, it is an improvement over the formula proposed by Guth and Weinberg  for first order cosmological phase transitions in that it allows for completion of the transition: the Guth and Weinberg formula is only valid during the early stage of the transition.

Next we need a dynamical equation which couples the time evolution of the temperature to the fraction of space converted to the hadronic phase. The following analysis is very similar to the one used to study the electroweak phase transition in the early universe by Carrington and me . We use Einstein’s equations as applied to the early universe, neglecting curvature. The evolution of the energy density is

$$\frac{dϵ}{dR}=\frac{3w}{R},$$
(21)
where $`R`$ is the scale factor at time $`t`$. This assumes kinetic but not phase equilibrium, and is basically a statement of energy conservation. We express the energy density as

$$ϵ=hϵ_h(T)+[1h]ϵ_q(T),$$
(22)
where $`ϵ_h`$ and $`ϵ_q`$ are the energy densities in the two phases at the temperature $`T`$. There is a similar equation for the enthalpy $`w`$. The time dependence of the scale factor is determined by the equation of motion

$$\frac{1}{R}\frac{dR}{dt}=\sqrt{\frac{8\pi Gϵ}{3}}.$$
(23)
This expression can be used to relate the time to the scale factor using the normalization $`R(t_c)=1`$.

We also need to know how fast a bubble expands once it is created. This is a subtle issue since by definition a critical size bubble is metastable and will not grow without a perturbation. After applying a perturbation, a critical size bubble begins to grow. As the radius increases, the surface curvature decreases, and an asymptotic interfacial velocity is approached. The asymptotic radial growth velocity will be referred to as $`v(T)`$. The expected qualitative behavior of $`v(T)`$ is that the closer $`T`$ is to $`T_c`$ the slower the bubbles grow. At $`T_c`$ there is no motivation for bubbles to grow at all since one phase is as good as the other. The bubble growth velocity was studied by Miller and Pantano . Their hydrodynamical results may be parameterized by the simple formula

$$v\gamma =3\left(1\frac{T}{T_c}\right)^{3/2}$$
(24)
which indeed has the expected behavior. A simple illustrative model for bubble growth is then

$$V(t^{},t)=\frac{4\pi }{3}\left(r_{}(T(t^{}))+_t^{}^t𝑑t^{\prime \prime }v(T(t^{\prime \prime }))\right)^3.$$
(25)
This expression can also be written in terms of $`R,R^{},R^{\prime \prime }`$ instead of $`t,t^{},t^{\prime \prime }`$.

We now have a complete set of coupled integro-differential equations which must be solved numerically. These equation take into account bubble nucleation and growth, energy conservation, and Einstein’s equations. They make no assumption about entropy conservation.

Figure 5 shows the temperature as a function of the scale factor.

For practical purposes, nucleation begins near the bottom of the cooling line. Thereafter, nucleation and growth of bubbles releases latent heat, which causes the temperature to rise. The increasing temperature shuts off nucleation, and the phase transition continues due to the growth of already nucleated bubbles. The temperature can never quite reach $`T_c`$; if it did, bubble growth would cease and the transition would never complete. This is a result of the equations of motion and is not an imposition.

Figure 6 shows the average bubble density

$$n(R(t))=_{t_c}^t𝑑t^{}I(T(t^{}))[1h(t^{})]$$
(26)
as a function of the scale factor. The bubble density rises rapidly just before $`R`$ reaches 1.007 and reaches its asymptotic value just after 1.007.

Figure 7 shows the nucleation rate as a function of scale factor. The rate has a very sharp maximum between 1.0070 and 1.0071. The turn on and turn off of the nucleation rate corresponds precisely with the fall and rise of the temperature shown in figure 5.

Figure 8 shows the fraction of space $`h`$ which has made the conversion to the hadronic phase. When $`h=1`$ the transition is complete and the temperature will begin to fall again. This occurs when $`R_f1.4464`$, to be compared with the value one would obtain from an ideal Maxwell construction $`R_{\mathrm{Maxwell}}=(239/79)^{1/3}=1.44630\mathrm{}`$. In fact the whole curve $`h(R)`$ is very close to the ideal Maxwell construction, apart from its delayed start apparent in the figure. The interested student will work out the Maxwell formula from the equations given here.

Figure 9 shows the average bubble radius as a function of scale factor.

$$\frac{4\pi }{3}\overline{r}^3n=h$$
(27)
It grows with time and with the scale factor, of course. At the end of the phase transition it is of the order of 1 cm. This is also the order of magnitude of the distance between the final quark-gluon plasma regions. Unfortunately, this is two orders of magnitude too small to affect nucleosynthesis. This result is rather robust against reasonable variations in any of the input parameters.

## 6 Conclusion

At temperatures exceeding approximately 160-180 MeV matter takes the form of quark-gluon plasma rather than hadronic gas. In principle QCD can be used to calculate the partition function and all other properties of interest. Perturbation theory has been employed to calculate the equation of state at asymptotically high temperatures to order $`g^5\alpha _s^{5/2}`$, but the results diverge at low temperature and cannot predict the properties of a phase transition. For this purpose numerical Monte Carlo simulations of lattice gauge theory are necessary. For the real world with 2 light flavors and one moderately light flavor of quark there is a big increase of the energy density in a narrow window of temperature centered on 170 MeV. Despite twenty years of effort it is still not known whether there is a true thermodynamic phase transition of first or second order or just a very rapid crossover. This is due to the great difficulty of computation with light dynamical quarks on the lattice.

Nucleosynthesis is affected by remnant inhomogeneities in the baryon to entropy ratio and isospin if the high baryon density regions immediately following a QCD phase transition are separated by at least 1 m. A set of dynamical equations can be written and solved for the evolution of the universe through such a phase transition all the way to completion. The evolution of the temperature and hadronic volume fraction as functions of time and scale factor are hardly any different than the results of an idealized Maxwell construction. The information not available in the latter construction is the length scale of the inhomogeneities, that is, bubble sizes and so on. The characteristic distance between the last regions of quark-gluon plasma seem to be on the order of 1 cm, too small to affect nucleosynthesis. However, there are qualifications and improvements can be made. For example, when the fraction of space occupied by bubbles exceeds about 50%, interactions among the bubbles probably cannot be neglected. It is unlikely, though, that further improvements in the dynamics would qualitatively change the current picture of the transition. Indeed, crude estimates of the effects of bubble fusion on the dynamics of the QCD transition in heavy ion collisions indicate that the transition completes only a little faster, and that the average bubble size is greater . At least that is in the right direction to be interesting.

The Relativistic Heavy Ion Collider (RHIC) is an accelerator at Brookhaven National Laboratory designed and constructed explicitly to search for quark-gluon plasma and the nature of a possible phase transition. Operations began in the summer of 2000, and a wealth of data is expected soon. Essentially RHIC is attempting to reproduce the first microsecond of the early universe, to produce matter with temperatures exceeding well over 200 MeV. Many of us look forward to exploitation of this little big bang in the laboratory to cosmology .

## Acknowledgements

This work was supported by the US Department of Energy under grant DE-FG02-87ER40328.


format: text

text:
# Dynamic instability of a rotating Bose-Einstein condensate

## Abstract

We consider a Bose-Einstein condensate subject to a rotating harmonic potential, in connection with recent experiments leading to the formation of vortices. We use the classical hydrodynamic approximation to the non-linear Schrödinger equation to determine almost analytically the evolution of the condensate. We predict that this evolution can exhibit dynamical instabilities, for the stirring procedure previously demonstrated at ENS and for a new stirring procedure that we put forward. These instabilities take place within the range of stirring frequency and amplitude for which vortices are produced experimentally. They provide therefore an initiating mechanism for vortex nucleation.

Quantized vortices in superfluid helium II, in particular the issue of vortex nucleation in a rotating container, have long been the subject of intense work . With the recent production of gaseous Bose-Einstein condensates  the subject has gained a renewed interest. On the experimental side, three groups have succeeded in obtaining vortices in atomic condensates, with two different techniques: a phase imprinting technique at JILA  and the equivalent of the helium rotating bucket experiment at ENS  and at MIT . At ENS a rotating laser beam superposed onto the magnetic trap holding the atoms creates a harmonic rotating potential with adjustable anisotropy $`ϵ`$ and rotation frequency $`\mathrm{\Omega }`$. For a well chosen range of variation for $`\mathrm{\Omega }`$ one or several vortices are nucleated, and then detected as holes in the density profile of the gas after ballistic expansion  or by a measurement of the angular momentum of the condensate . A striking feature of the ENS experimental results is that, for a very weak anisotropy $`ϵ`$, nucleation of vortices takes place in a narrow interval of rotation frequencies $`[\mathrm{\Omega }_{\text{min}},\mathrm{\Omega }_{\text{max}}]`$ around $`0.7\omega _{}`$, where $`\omega _{}`$ is the mean oscillation frequency of the atoms in the $`xy`$ plane, whatever the number of atoms or the oscillation frequency $`\omega _z`$ along $`z`$ in the experiment .

While the JILA experiment is well understood theoretically  the situation is more involved for the ENS experiment. Several theoretical articles, inspired by the case of superfluid helium, have tried to predict the value of the lower vortex nucleation frequency $`\mathrm{\Omega }_{\text{min}}`$ from purely thermodynamic arguments . The proposed values for $`\mathrm{\Omega }_{\text{min}}`$ are significantly different from the observed value of $`0.7\omega _{}`$, or depend on the trap aspect ratio $`\omega _z/\omega _{}`$ or on the atom number, in contradiction with the observations at ENS. Also thermodynamical reasonings are not able to predict the upper vortex nucleation frequency $`\mathrm{\Omega }_{\text{max}}`$, which is also close to $`0.7\omega _{}`$ for low anisotropy $`ϵ`$.

In this paper we consider the time dependent problem of a condensate subject to a harmonic stirring potential. We use the classical hydrodynamic approximation to the time dependent Gross-Pitaevskii equation (GPE), an approximation well justified for the ENS parameters . We are then able to reformulate the partial differential hydrodynamic equations in terms of ordinary differential equations, which allows an almost analytical solution . Our main result is the discovery of dynamical instabilities in the evolution of the condensate for a certain range of the rotation frequency and of the trap anisotropy. These instabilities will invalidate the classical hydrodynamic approximation after some evolution time. We have checked with a numerical solution of the Gross-Pitaevskii equation that vortices then enter the condensate.

The existence of such a dynamical instability explains why in earlier numerical work the time dependent Gross-Pitaevskii equation was found to nucleate vortices . Furthermore the instability range that we predict is very close to the experimentally observed range of vortex nucleation, for various stirring procedures. For the stirring procedure of  we recover the “universal” numerical value $`0.7`$ for $`\mathrm{\Omega }/\omega _{}`$ leading to vortex nucleation for low anisotropies $`ϵ`$. We provide a simple physical interpretation of this value: for $`\mathrm{\Omega }/\omega _{}=1/\sqrt{2}0.7`$ the harmonic stirring potential resonantly excites a quadrupole mode of the condensate , which induces large oscillations of the condensate and eventually a dynamical instability sets in. We also investigate a new excitation procedure to nucleate vortices, that has recently been implemented at ENS : the rotation of the stirring potential is set up very slowly. The gas then follows adiabatically a branch of steady state until the branch becomes dynamically unstable. The corresponding lower rotation frequency $`\mathrm{\Omega }`$ that we predict for vortex nucleation is also very close to the experimental value.

In our model atoms are trapped in a harmonic potential rotating at the instantaneous frequency $`\mathrm{\Omega }(t)`$ around $`z`$ axis. For convenience all the calculations of this paper are done in a rotating frame where the trap axes are fixed. The trapping potential then reads:

$$U(\stackrel{}{r},t)=\frac{1}{2}m\omega _{}^2\left\{[1ϵ(t)]x^2+[1+ϵ(t)]y^2+\left(\frac{\omega _z}{\omega _{}}\right)^2z^2\right\}$$
(1)
where $`m`$ is the mass of an atom, $`ϵ(t)`$ is the trap anisotropy at time $`t`$. The parameters $`\omega _{}`$ and $`\omega _z`$ are the oscillation frequencies of the atoms in transverse and axial directions for vanishing anisotropy of the stirring potential. Within the mean field approximation, the time evolution of the condensate field or macroscopic wavefunction $`\psi (\stackrel{}{r},t)`$ can be described by the time dependent Gross-Pitaevskii equation :

$$i\mathrm{}\frac{\psi }{t}=\left[\frac{\mathrm{}^2}{2m}\stackrel{}{}^2+U(\stackrel{}{r},t)+g|\psi |^2\mathrm{\Omega }(t)\widehat{L}_z\right]\psi ,$$
(2)
where $`g=4\pi \mathrm{}^2a/m`$ is the coupling constant, proportional to the $`s`$-wave scattering length $`a`$ of the atoms, here taken to be positive, and where the inertial term proportional to the angular momentum operator $`\widehat{L}_z`$ along $`z`$-axis accounts for the frame rotation. The condensate field $`\psi `$ can be written in terms of density $`\rho `$ and phase $`S`$,

$$\psi (\stackrel{}{r},t)=\sqrt{\rho (\stackrel{}{r},t)}e^{iS(\stackrel{}{r},t)/\mathrm{}}.$$
(3)
The equation obtained from the GPE for $`\rho `$ is just the continuity equation. The equation for $`S`$ contains the so-called quantum pressure term $`\mathrm{}^2\stackrel{}{}^2\sqrt{\rho }/2m\sqrt{\rho }`$ that we neglect here as compared to the mean-field term $`\rho g`$ in the Thomas-Fermi approximation. We obtain:

$`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`\text{div}\left[\rho \left({\displaystyle \frac{\stackrel{}{}S}{m}}\stackrel{}{\mathrm{\Omega }}(t)\times \stackrel{}{r}\right)\right]`$ (4)
$`{\displaystyle \frac{S}{t}}`$ $`=`$ $`{\displaystyle \frac{(\stackrel{}{}S)^2}{2m}}+U(\stackrel{}{r},t)+g\rho (\stackrel{}{\mathrm{\Omega }}(t)\times \stackrel{}{r})\stackrel{}{}S.`$ (5)
A very fortunate feature of the harmonic trap is that these superfluid hydrodynamic equations can be solved exactly for a condensate initially at equilibrium in the non-rotating trap with the following quadratic ansatz for the condensate density and phase :

$`\rho _c(\stackrel{}{r},t)`$ $`=`$ $`\rho _0(t)+{\displaystyle \frac{m\omega _{}^2}{g}}{\displaystyle \underset{i,j=1}{\overset{3}{}}}x_iA_{ij}(t)x_j,`$ (6)
$`S_c(\stackrel{}{r},t)`$ $`=`$ $`s_0(t)+m\omega _{}{\displaystyle \underset{i,j=1}{\overset{3}{}}}x_iB_{ij}(t)x_j,`$ (7)
where $`x_1`$, $`x_2`$ and $`x_3`$ are the coordinates along $`x`$, $`y`$ and $`z`$ axes respectively. The time dependent dimensionless coefficients $`A_{ij}`$ and $`B_{ij}`$ form $`3\times 3`$ symmetric matrices $`A`$ and $`B`$ which from Eqs.(4,5) obey the evolution equations:

$`\omega _{}^1{\displaystyle \frac{dA}{dt}}`$ $`=`$ $`2A\text{Tr}B2\{A,B\}+{\displaystyle \frac{\mathrm{\Omega }}{\omega _{}}}[R,A],`$ (8)
$`\omega _{}^1{\displaystyle \frac{dB}{dt}}`$ $`=`$ $`2B^2WA+{\displaystyle \frac{\mathrm{\Omega }}{\omega _{}}}[R,B]`$ (9)
where $`\{,\}`$ stands for the anti-commutator, $`[,]`$ stands for the commutator of two matrices, the matrix $`W`$ is diagonal, with components $`W_{11}=(1ϵ(t))/2`$, $`W_{22}=(1+ϵ(t))/2`$, and $`W_{33}=(\omega _z/\omega _{})^2/2`$, and the matrix $`R`$, originating from the vectorial product in $`\widehat{L}_z`$, has vanishing elements except for $`R_{12}=R_{21}=1`$ . Note that these equations do not depend on the number of atoms nor on the coupling constant $`g`$.

In a first stage, it is very important to study steady state solutions of the above equations. We restrict to solutions that have the same symmetry as the initial state: even parity along $`z`$, this parity being preserved by time evolution. We then find an unique class of solutions, reproducing the results of : the condensate phase varies as $`S(\stackrel{}{r})=m\omega _{}\beta xy`$ where $`\beta `$ is a real root of

$$\beta ^3+\left(12\frac{\mathrm{\Omega }^2}{\omega _{}^2}\right)\beta \frac{\mathrm{\Omega }}{\omega _{}}ϵ=0.$$
(10)
The steady state matrix $`A`$ is diagonal with elements given in . We have plotted in figure 1 the values of $`\beta `$ as function of the rotation frequency $`\mathrm{\Omega }`$ for a fixed anisotropy $`ϵ`$. For $`\mathrm{\Omega }`$ between zero and a bifurcation value $`\mathrm{\Omega }_{\text{bif}}(ϵ)`$ depending on $`ϵ`$ there is a single branch of solution for $`\beta `$. This branch is supplemented by two extra branches when $`\mathrm{\Omega }>\mathrm{\Omega }_{\text{bif}}(ϵ)`$.

We now turn back to the time dependent problem. Clearly a condensate with a vortex cannot be described within the quadratic ansatz (6,7) as the phase $`S_c`$ corresponds to an irrotational velocity flow. The actual scenario for the vortex nucleation that we put forward is the following: initially very small deviations $`\delta \rho (\stackrel{}{r},t)`$ of the condensate density and $`\delta S(\stackrel{}{r},t)`$ of the condensate phase from the quadratic shapes $`\rho _c`$ and $`S_c`$ may grow exponentially fast in the course of time evolution, eventually leading the condensate to a structure very different from Eqs.(6,7). This may happen when a dynamical instability is present.

To reveal such an instability we obtain from the evolution equations (4, 5) linearized equations of motion for initially small deviations $`\delta \rho `$ and $`\delta S`$ from $`\rho _c`$ and $`S_c`$:

$`{\displaystyle \frac{D\delta \rho }{Dt}}`$ $`=`$ $`\text{div}\left(\rho _c{\displaystyle \frac{\stackrel{}{}\delta S}{m}}\right)\delta \rho {\displaystyle \frac{\stackrel{}{}^2S_c}{m}},`$ (11)
$`{\displaystyle \frac{D\delta S}{Dt}}`$ $`=`$ $`g\delta \rho .`$ (12)
In these equations, we have introduced the convective derivative $`\frac{D}{Dt}\frac{}{t}+\stackrel{}{v}_c(\stackrel{}{r},t)\stackrel{}{}`$ where $`\stackrel{}{v}_c=\stackrel{}{}S_c/m\stackrel{}{\mathrm{\Omega }}\times \stackrel{}{r}`$ is the velocity field of the condensate in the rotating frame. A polynomial ansatz for $`\delta S`$ and $`\delta \rho `$ of an arbitrary total degree $`n`$ in the coordinates $`x`$, $`y`$ and $`z`$ solves these linear equations exactly . This is another nice consequence of the harmonicity of the trap. In practice, we calculate the evolution operator $`𝒰_n(t)`$ mapping the coefficients of the polynomials at time zero onto their values after a time evolution $`t`$. Dynamical instability takes place when one or several eigenvalues of $`𝒰_n`$ grow exponentially fast with time $`t`$. Note that after rescaling of the variables, Eqs.(11,12) become independent of the number of atoms and of the coupling constant $`g`$, in a way similar to Eqs.(8,9).

Now we perform a linear stability analysis for two different stirring procedures of the condensate.

Procedure I: The ellipticity of the stirring potential $`ϵ`$ is kept fixed and the rotation frequency $`\mathrm{\Omega }(t)`$ of the stirrer is very slowly ramped up from zero to its final value. The condensate, initially in steady state with a vanishing parameter $`\beta `$, adiabatically follows the upper branch of steady states with $`\beta 0`$, see figure 1. It is then sufficient to determine the dynamic stability of the upper branch of steady states. This greatly simplifies the calculation as one just has to identify eigenmodes of Eqs.(11,12) evolving in time as $`\mathrm{exp}(\lambda t)`$ where the eigenvalue $`\lambda `$ is a complex number. Dynamical instability takes place when $`\lambda `$ can have a strictly positive real part. As shown in figure 1 for $`ϵ=0.022`$ we find that the upper branch for $`\beta `$ is stable for modes of degree $`n=2`$ but presents two instability intervals for the modes of degree $`n=3`$, a very narrow interval around $`\mathrm{\Omega }=0.58\omega _{}`$ and a broader interval starting at $`\mathrm{\Omega }=0.778\omega _{}`$.

We have investigated in a systematic way the instability range of the upper branch of steady states, by varying the anisotropy $`ϵ`$, the rotation frequency $`\mathrm{\Omega }`$ and the degree $`n`$. The instability domain in the $`\mathrm{\Omega }ϵ`$ plane for modes of degree 3 is mainly made of a crescent, and the inclusion of higher degree modes ($`n=4,5`$) add extra crescents from above, see figure 2a. Each crescent has on the $`ϵ=0`$ axis (i) a broad basis at $`\mathrm{\Omega }>\omega _{}/\sqrt{2}`$, with a non-zero instability exponent, and (ii) a very narrow edge at $`\mathrm{\Omega }=\omega _{}/\sqrt{n}<\omega _{}/\sqrt{2}`$ with a vanishing instability exponent. We show in figure 2b that the maximal instability exponent for $`n=3`$ has a remarkably weak dependence on $`\omega _z/\omega _{}`$.

Procedure II: This is the original experimental scenario of , where the stirring potential is rotated at a fixed frequency and the ellipticity of the stirrer is turned on from zero to its final value $`ϵ_f`$ abruptly. In this case we cannot rely on adiabatic following for the condensate density and phase, so that we solve the time dependent equations (8,9) for $`\rho _c(t)`$ and $`S_c(t)`$. Then we perform a linear stability analysis as discussed above: we evolve a generic polynomial ansatz of degree $`n`$ for the fluctuations $`\delta \rho `$ and $`\delta S`$ according to Eqs.(11,12), which allows to construct the evolution operator $`𝒰_n(t)`$ and to calculate $`Z_{\text{max}}(t)`$, the eigenvalue of $`𝒰_n(t)`$ with the largest modulus. Then we define the mean instability exponent Re $`\lambda `$ as the mean slope of $`\mathrm{ln}|Z_{\text{max}}(t)|`$ as function of time.

This reveals that within certain range of rotation frequency the system becomes dynamically unstable, see the solid line in figure 3. In the limit of a low anisotropy $`ϵ`$ the instability sets in when the rotation frequency $`\mathrm{\Omega }`$ is close to the value $`0.7\omega _{}`$: in the lab frame, the stirring potential of frequency $`2\mathrm{\Omega }`$ is then resonant with a quadrupole mode  of the condensate of frequency $`\sqrt{2}\omega _{}`$, and induces large amplitude oscillations of the condensate, resulting in a dynamical instability. More precisely, the condensate described by the quadratic ansatz $`\rho _c,S_c`$ has an angular momentum oscillating in time around a non-zero value $`\overline{L}_z`$. This value $`\overline{L}_z`$ is a peaked function of the rotation frequency $`\mathrm{\Omega }`$, shown in dashed line in figure 3. The peak of $`\overline{L}_z`$ is not exactly located at $`\mathrm{\Omega }=0.7\omega _{}`$ because of non-linear effects in Eqs.(8, 9). The peak structure of the instability exponent in figure 3 is alike the peak structure of $`\overline{L}_z`$, with a narrower width as dynamical instability of the vortex free solution $`\rho _c,S_c`$ sets in for the higher values of $`\overline{L}_z`$ only. For values of $`\mathrm{\Omega }`$ significantly above or below 0.7$`\omega _{}`$ the stirrer is out of resonance with the quadrupole mode and induces only small and stable oscillations of the condensate. For larger values of $`ϵ`$, the instability interval in $`\mathrm{\Omega }`$ broadens. We have also checked that the instability interval depends weakly on $`\omega _z/\omega _{}`$.

What is the connection between the dynamical instabilities found here and the nucleation of vortices ? To obtain a theoretical answer to this question, one has to go beyond a linear stability analysis to determine the evolution of the condensate in the long run: for a few values of the rotation frequency $`\mathrm{\Omega }`$ and for procedures I and II, we have checked by a numerical integration of the time dependent GPE in three dimensions, that vortices are indeed nucleated in the predicted instability domains: after some evolution time, the angular momentum in the numerical solution suddenly becomes larger than the classical hydrodynamic prediction, as vortices enter in the condensate. An experimental answer to this question for the stirring procedure I has been provided recently at ENS : the first clear evidence of a vortex appears for a rotation frequency $`\mathrm{\Omega }=0.77\omega _{}`$, very close to our prediction $`\mathrm{\Omega }=0.778\omega _{}`$. The agreement is also very good for procedure II as shown in figure 3: our instability domain in $`\mathrm{\Omega }`$ coincides with the experimental vortex nucleation interval within a few percent.

In summary, the dynamical instabilities that we have identified provide an initiating mechanism for the production of vortices in a condensate stirred by a harmonic potential, in excellent agreement with the experimental results at ENS.

We thank S. Rica, V. Hakim, G. Shlyapnikov, F. Chevy, K. Madison and J. Dalibard for helpful discussions. We acknowledge financial support from Ministère de la Recherche et de la Technologie. LKB is a unité de recherche de l’Ecole normale supérieure et de l’Université Pierre et Marie Curie, associée au CNRS.


text:
# Pythagoras’ Theorem on a 2D-Lattice from a “Natural” Dirac Operator and Connes’ Distance Formula

## I Introduction

Lattice Dirac operator is a long-standing problem embarrassing lattice field theorists. No-Go theorem  makes the implementation of chiral fermions on lattices almost impossible. Recent years, some breakthroughs have been achieved, e.g. the rediscovery of Ginsparg-Wilson relation , the devices of domainwall and overlap Dirac operators . On the other hand, lattices can be considered as a simplest realization of noncommutative geometry (NCG) which has drawn more and more attention of theoretical physicists due to its applications in standard model of particle physics , lattice field theory , and string/M-theory . NCG provides a powerful candidate of mathematical framework for geometrical understanding of fundamental physical laws. In Alain Connes’s version of NCG, a generalized Dirac operator plays a central role in introducing the metric structure onto a noncommutative space . Then intuitively, it becomes an interesting question whether Connes’ NCG idea could brighten the problem of lattice Dirac operator. In fact, some groups have explored this question. G. Bimonte et al first pointed out that naïve Dirac operator is not able to induce the conventional distance on a 4D lattice by Connes’ construction and Wilson-Dirac operator gives an even worse result . Starting from this observation, E. Atzmon computed this “anomalous distance” induced by naïve Dirac operator for an 1D lattice precisely , which gives

$$d(0,2n1)=2n,d(0,2n)=2\sqrt{n(n+1)}(n𝒩)$$

Hence, we can say that NCG provides another criterion for the rational choices of Dirac operator on lattices. A. Dimakis et al discovered an 1D operator whose Connes’ distance coincides with the usual linear distance, when the number of lattice sites is finite . This serial work indicates that it is a highly nontrivial problem to devise a proper Dirac operator whose Connes’ distance is a desired one.

In this paper, we construct a “natural” Dirac operator on a 2D lattice. The induced metric of this operator endows conventional Euclidean geometry on this lattice, i.e. the Pythagoras’ Theorem holds for the Connes’ distance. The restoration of metric relies on that this operator has a so-called “local eigenvalue property”, such that the norm $`[D,f]`$ is completely solvable (diagonalizable, integrable). The Dirac operator devised by us can be regarded as a generalization of Dimakis’ operator in 1D case; it can be generalized into any high dimensional lattices easily.

The paper is organized as follows. In section II, we give a brief introduction to Connes’ NCG and his distance formula. The “natural” Dirac operator is defined in section III. In section IV, we explore the “local eigenvalue property” of our Dirac operator. In section V, we prove that our Dirac operator implies the Pythagoras’ Theorem on a 2D Lattice. In section VI, we generalize the “natural” Dirac operator into any higher dimensional lattices and make some open discussions.

## II Noncommutative Geometry and Connes’ Distance Formula

The object of a Connes’ NCG is a triple $`(A,H,D)`$ called a K-cycle, in which $`A`$ is an involutive algebra, represented faithfully and unitarily as a subalgebra of bounded operators on a Hilbert space $`H`$, and $`D`$ is a self-adjoint operator on $`H`$, which is called generalized Dirac operator, with compact resolvent so that $`[D,\widehat{a}]`$ is bounded for all $`a`$ in $`A`$. Here we use $`\widehat{a}`$ to denote the imagine of $`a`$ on $`H`$; without introducing any misunderstanding below, we will just omit the hat on $`a`$. A K-cycle is required to satisfy some axioms such that it recovers the ordinary spin geometry on a differential manifold when $`A`$ is taken to be the algebra of smooth functions over this spin manifold, $`H`$ is the space of $`L^2`$-spinors and $`D`$ is the classical Dirac operator.

From a K-cycle, we can define a metric $`d_D(,)`$ on the state space of $`A`$ denoted as $`S(A)`$, which we have referred as induced metric or Connes’ distance.

$$d_D(\varphi ,\psi )=sup\{|\varphi (a)\psi (a)|/aA,[D,a]1\}$$
(1)
for any $`\varphi ,\psi S(A)`$. One can check that this definition satisfies all of the three axioms for a metric easily. Once $`A`$ is commutative, pure states correspond to characters which can be interpreted as points, with $`A`$ being the algebra of functions over these points (Gel’fand-Naimark Theorem); and Eq.(1) can, under this circumstance, be rewritten as

$$d_D(p,q)=sup\{|f(p)f(q)|/fA,[D,f]1\},p,qS(A)$$
(2)
Noticing the fact that

$$a^2=sup\{(a\psi ,a\psi )/\psi H,\psi =1\}=sup\{\lambda (a)a^{}a\psi =\lambda (a)\psi \},aB(H)$$

where $`B(H)`$ is the algebra of bounded operators on $`H`$, we can express the inequality constraint in Eq.(2) by using eigenvalues

$$d_D(p,q)=sup\{|f(p)f(q)|/fA,\lambda ([D,f]^{}[D,f])1\},p,qS(A)$$
(3)
Notations
$`\sigma _i,i=1,2,3`$ are Pauli matrices defined in the ordinary way. If $`S`$ is a finite set, then $`|S|`$ is the number of the elements in $`S`$. Let $`i=1,2,3`$, define $`\gamma `$-matrices as

$$\gamma ^i=\left(\begin{array}{cc}& \sigma _i\\ \sigma _i& \end{array}\right);\gamma ^4=\left(\begin{array}{cc}& i\\ i& \end{array}\right)$$

The vacant matrix elements are understood as zeros and this will be taken as a convention all through this paper. Introduce $`\gamma _\pm ^1=\frac{1}{2}(\gamma ^1\pm i\gamma ^2),\gamma _\pm ^2=\frac{1}{2}(\gamma ^3\pm i\gamma ^4)`$, whose explicit matrix representations are

$$\gamma _+^1=\left(\begin{array}{cccc}& & & 1\\ & & 0& \\ & 1& & \\ 0& & & \end{array}\right);\gamma _{}^1=\left(\begin{array}{cccc}& & & 0\\ & & 1& \\ & 0& & \\ 1& & & \end{array}\right);$$

$$\gamma _+^2=\left(\begin{array}{cccc}& & 0& \\ & & & 1\\ 1& & & \\ & 0& & \end{array}\right);\gamma _{}^2=\left(\begin{array}{cccc}& & 1& \\ & & & 0\\ 0& & & \\ & 1& & \end{array}\right)$$

satisfying Clifford algebra relations

$$\{\gamma _\pm ^i,\gamma _\pm ^j\}=0,\{\gamma _\pm ^i,\gamma _{}^j\}=\delta ^{ij},i,j=1,2$$

$`𝒵`$ is referred to the set of integrals and $`𝒞`$ is for complex numbers. We adopt the convention $`a=1`$ for the lattice constant $`a`$ in this paper.

## III “Natural” Dirac Operator on a 2D Lattice

First we give a detailed re-formulation of Dimakis’ operator on the 1D lattice. Let $`L_1=\{x|x𝒵\}`$, $`A_0=\{f:L_1𝒞\},H=l^2(L_1)𝒞^2=\{\psi \}`$. The operator developed by Dimakis et al acting on $`H`$ can be written in the form

$$D=\sigma _+^++\sigma _{}^{}=\left(\begin{array}{cc}& ^+\\ ^{}& \end{array}\right)$$

where $`\sigma _\pm =(\sigma _1\pm i\sigma _2)/2`$, $`(T^\pm \psi )(x):=\psi (x\pm 1)`$, $`(^\pm \psi )(x)=((T^\pm \mathrm{𝟏})\psi )(x)=\psi (x\pm 1)\psi (x)`$ and $`\psi (x)=(\psi _1,\psi _2)^T(x)`$. Refine $`A_0`$ to be $`A=\{fA_0/[D,f]<\mathrm{}\}`$.

We point out that $`D`$ has a “local eigenvalue property”, showing below

$$[D,f]=\left(\begin{array}{cc}& (^+f)T^+\\ (^{}f)T^{}& \end{array}\right)=\left(\begin{array}{cc}& T^+(^{}f)\\ T^{}(^+f)& \end{array}\right)$$

$$[D,f]^{}=[D,\overline{f}]$$

Introduce a “Hamiltonian”$`H(df):=[D,f]^{}[D,f]`$

$$=\left(\begin{array}{cc}|^+f|^2& \\ & |^{}f|^2\end{array}\right)$$

Consider the eigenvalue equation

$$H(df)\psi =\lambda (df)\psi $$

$$|^+f|^2(x)\psi _1(x)=\lambda \psi _1(x),x𝒵$$

$$|^{}f|^2(x)\psi _2(x)=\lambda \psi _2(x),x𝒵$$

So that the constraint in Eq.(3) implies

$$\lambda (x,1)=|^+f|^2(x)1,x𝒵$$
(4)
$$\lambda (x,2)=|^{}f|^2(x)1,x𝒵$$
(5)
Now we observe that each eigenvalue of $`H(df)`$ is just related to one link of $`L_1`$, to which we refer as “local eigenvalue”. As a comparison, the eigenvalue equation for 1D naïve Dirac operator $`D_N:=(T^+T^{})/2`$ is

$$\frac{1}{4}(|^+f|^2(x)+|^{}f|^2(x))\psi (x)+\frac{1}{4}((^+\overline{f})(x)(^+f)(x+1)\psi (x+2)+(^{}\overline{f})(x)(^{}f)(x1)\psi (x2))=\lambda \psi (x)$$

for all x in $`𝒵`$, with $`A_N=A`$, $`H_N=l^2(L_1)`$, which possesses no “local eigenvalues” evidently.

Now we induce $`d_D(,)`$ on $`L_1`$. Let $`fA`$ subjected to (4)(5), i.e. $`[D,f]1`$, so that $`|f((x+m)f(x)|m,m=1,2,3,\mathrm{}`$; hence, $`d_D(x+m,x)`$ has an upper bound $`m`$. Let $`f_0(x)=x`$, then $`f_0`$ saturates this upper bound. Thus, we have proved $`d_D(m,n)=|mn|,m,nL_1`$.

The 2D lattice being considered as a set is parametrized as $`L_2=\{x=(m,n)m,n𝒵\}`$.

$$A_0:=\{f:L_2𝒞\},H:=l^2(L_2)𝒞^4=\{\psi :\psi =(\psi _1,\psi _2,\psi _3,\psi _4)^T\}$$

$$D=\underset{i=1}{\overset{2}{}}\underset{s=\pm }{}\gamma _s^i_i^s=\left(\begin{array}{cccc}& & _2^{}& _1^+\\ & & _1^{}& _2^+\\ _2^+& _1^+& & \\ _1^{}& _2^{}& & \end{array}\right)$$
(6)
$$(T_i^\pm \psi )(x)=\psi (x\pm \widehat{i})$$

$$(_i^\pm \psi )(x)=((T_i^\pm \mathrm{𝟏})\psi )(x)=\psi (x\pm \widehat{i})\psi (x)$$

where $`\widehat{1}=(1,0),\widehat{2}=(0,1)`$. Refine $`A_0`$ to be $`A:=\{fA_0/[D,f]<\mathrm{}\}`$. By symmetries of $`D`$, we have

$$d_D((m,n),(m^{},n^{}))=d_D((0,0),(m^{}m,n^{}n))$$

so that to consider $`d_D((0,0),(m,n))`$ is enough. What’s more,

$$d_D((0,0),(m,n))=d_D((0,0),(n,m)),d_D((0,0),(m,n))=d_D((0,0),(n,m)),$$

so we just need to consider $`mn0`$.

To end this section , we adopt the method in  to prove

$$d_D((0,0),(i,0))=i,i=1,2,\mathrm{}$$
(7)
###### Lemma 1

Let $`fA`$, $`[D,f]1`$ and $`\stackrel{~}{f}(m,n):=f(m,0)`$, then $`[D,\stackrel{~}{f}]1`$.

Proof:
By the definition of $`D`$ in (6), there is

$$[D,\stackrel{~}{f}]\psi =((_1^+\stackrel{~}{f})(T_1^+\psi _4),(_1^{}\stackrel{~}{f})(T_1^{}\psi _3),(_1^+\stackrel{~}{f})(T_1^+\psi _2),(_1^{}\stackrel{~}{f})(T_1^{}\psi _1))^T$$

Thus, notice the definition of $`\stackrel{~}{f}`$,

$$([D,\stackrel{~}{f}]\psi ,[D,\stackrel{~}{f}]\psi )=\underset{m,n}{}(|_1^+\stackrel{~}{f}|^2(m,n)(\underset{i=1}{\overset{4}{}}|\psi _i^{}|^2(m,n)))$$

$$=\underset{m}{}(|_1^+f|^2(m,0)(\underset{n}{}\underset{i=1}{\overset{4}{}}|\psi _i^{}|^2(m,n)))$$

$$[D,\stackrel{~}{f}]=sup\{|_1^+f|(m,0)/m𝒵\}<\mathrm{}$$
(8)
Define $`\widehat{H}:=\{\psi H/\psi _i=0,i=1,2,3\}`$, then

$$[D,f]_{\widehat{H}}[D,f]1$$
(9)
and

$$([D,f]\psi ,[D,f]\psi )|_{\widehat{H}}=\underset{m,n}{}|_1^+f|^2(m,n)|\psi _4|^2(m+1,n)+|_2^+f|^2(m,n)|\psi _4|^2(m,n+1)$$

For any $`\psi H`$, we define $`\widehat{\psi }\widehat{H}`$ by that $`\widehat{\psi }_4(m,n)=0,n0`$ and that $`\widehat{\psi }_4(m,0)`$ satisfy
$`|\widehat{\psi }_4(m,0)|^2=_n_{i=1}^4|\psi _i^{}|^2(m,n)`$, for all $`m`$. Therefore,

$$([D,\stackrel{~}{f}]\psi ,[D,\stackrel{~}{f}]\psi )=\underset{m}{}(|_1^+f|^2(m,0)|\widehat{\psi }_4(m,0)|^2([D,f]\widehat{\psi },[D,f]\widehat{\psi })|_{\widehat{H}}[D,f]_{\widehat{H}}$$

Notice (9), and there are

$$[D,\stackrel{~}{f}][D,f]_{\widehat{H}}[D,f]1$$

$`\mathrm{}`$
Following this lemma, we can just consider $`f`$ with $`_2^+f=0`$ and reach the conclusion $`d_D((0,0),(i,0))i,i=1,2,\mathrm{}`$; $`f_0(m,n)=m`$ saturates this upper bound. Hence, Eq.(7) holds.

## IV Local Eigenvalues

We have to consider the eigenvalue problem in Eq.(3), when discussing distance $`d_D((0,0),(m,n)),mn1`$. Fortunately, our “natural” Dirac operator has the “local eigenvalue property” also. Here we give a very detailed calculation.

Let $`D(f):=[D,f]`$

$$=\left(\begin{array}{cccc}& & (_2^{}f)T_2^{}& (_1^+f)T_1^+\\ & & (_1^{}f)T_1^{}& (_2^+f)T_2^+\\ (_2^+f)T_2^+& (_1^+f)T_1^+& & \\ (_1^{}f)T_1^{}& (_2^{}f)T_2^{}& & \end{array}\right)$$

One can check $`[D,f]^{}=D(\overline{f})`$. Define the “Hamiltonian” $`H(df):=[D,f]^{}[D,f]=D(\overline{f})D(f)`$

$$=\left(\begin{array}{cc}|_1^+f|^2+|_2^{}f|^2& (T_2^{}(_2^+\overline{f}_1^+f)T_1^+(_1^{}\overline{f}_2^{}f))T_1^+T_2^{}\\ (T_1^{}(_1^+\overline{f}_2^+f)T_2^+(_2^{}\overline{f}_1^{}f))T_1^{}T_2^+& |_1^{}f|^2+|_2^+f|^2\end{array}\right)$$

$$\left(\begin{array}{cc}|_1^+f|^2+|_2^+f|^2& (T_2^+(_2^{}\overline{f}_1^+f)T_1^+(_1^{}\overline{f}_2^+f))T_1^+T_2^+\\ (T_1^{}(_1^+\overline{f}_2^{}f)T_2^{}(_2^+\overline{f}_1^{}f))T_1^{}T_2^{}& |_1^{}f|^2+|_2^{}f|^2\end{array}\right)$$

with the eigenvalue equation

$$H(df)\psi =\lambda (df)\psi $$
(10)
Fortunately, Eq.(10) can be reduced to a collection of equation sets, with each equation set being related to a fundamental plaque $`\{(m,n),(m+1,n),(m,n1),(m+1,n1)\}L_2`$

$$\left(\begin{array}{cc}\rho _1^2+\rho _4^2\lambda & \overline{\mathrm{\Delta }_4}\mathrm{\Delta }_3\overline{\mathrm{\Delta }_1}\mathrm{\Delta }_2\\ \mathrm{\Delta }_4\overline{\mathrm{\Delta }_3}\mathrm{\Delta }_1\overline{\mathrm{\Delta }_2}& \rho _2^2+\rho _3^2\lambda \end{array}\right)\left(\begin{array}{c}\psi _1(m,n)\\ \psi _2(m+1,n1)\end{array}\right)=0$$
(11)
$$\left(\begin{array}{cc}\rho _3^2+\rho _4^2\lambda & \overline{\mathrm{\Delta }_4}\mathrm{\Delta }_1+\overline{\mathrm{\Delta }_3}\mathrm{\Delta }_2\\ \mathrm{\Delta }_4\overline{\mathrm{\Delta }_1}+\mathrm{\Delta }_3\overline{\mathrm{\Delta }_2}& \rho _1^2+\rho _2^2\lambda \end{array}\right)\left(\begin{array}{c}\psi _3(m,n1)\\ \psi _4(m+1,n)\end{array}\right)=0$$
(12)
in which $`\mathrm{\Delta }_1(m,n):=(_1^+f)(m,n)`$, $`\mathrm{\Delta }_2(m,n):=(_2^+f)(m+1,n1)`$, $`\mathrm{\Delta }_3(m,n):=(_1^+f)(m,n1)`$, $`\mathrm{\Delta }_4(m,n):=(_2^+f)(m,n1)`$ and, without misleading, we have omitted the arguments $`(m,n)`$. Let $`\mathrm{\Delta }_i:=\rho _ie^{i\theta _i},i=1,2,3,4`$, $`\theta :=\theta _1+\theta _3\theta _2\theta _4`$, and

$$A:=\rho _1^2+\rho _4^2,B:=\rho _2^2+\rho _3^2,C:=\overline{\mathrm{\Delta }_4}\mathrm{\Delta }_3\overline{\mathrm{\Delta }_1}\mathrm{\Delta }_2$$

$$A^{}:=\rho _3^2+\rho _4^2,B^{}:=\rho _1^2+\rho _2^2,C^{}:=\overline{\mathrm{\Delta }_4}\mathrm{\Delta }_1+\overline{\mathrm{\Delta }_3}\mathrm{\Delta }_2$$

then the secular equation for Eq.(11) is

$$(\lambda A)(\lambda B)C\overline{C}=\lambda ^2(A+B)\lambda +ABC\overline{C}=0$$

in which

$$C\overline{C}=\rho _3^2\rho _4^2+\rho _1^2\rho _2^22\rho _1\rho _2\rho _3\rho _4\mathrm{cos}\theta $$

The solutions are

$$\lambda _\pm =\frac{1}{2}(A+B\pm \sqrt{(AB)^2+4C\overline{C}})$$
(13)
Similarly for Eq.(12),

$$\lambda _\pm ^{}=\frac{1}{2}(A^{}+B^{}\pm \sqrt{(A^{}B^{})^2+4C^{}\overline{C^{}}})$$
(14)
in which

$$C^{}\overline{C^{}}=\rho _1^2\rho _4^2+\rho _2^2\rho _3^22\rho _1\rho _2\rho _3\rho _4\mathrm{cos}\theta $$

Therefore, $`[D,f]1(m,n)`$,

$$\lambda _+(m,n)1,\lambda _+^{}(m,n)1$$
(15)
Equivalently, Eqs.(13)(14) and inequality (15) $``$

$$1+\rho _1^2\rho _3^2+\rho _2^2\rho _4^2+2\rho _1\rho _2\rho _3\rho _4\mathrm{cos}\theta \rho _1^2+\rho _2^2+\rho _3^2+\rho _4^2$$
(16)
$$2\rho _1^2+\rho _2^2+\rho _3^2+\rho _4^2$$
(17)
There is another constraint, the closedness condition

$$\mathrm{\Delta }_1+\mathrm{\Delta }_4=\mathrm{\Delta }_2+\mathrm{\Delta }_3$$
(18)
Hence we get (16)(17)(18) as the specific expressions of $`[D,f]1`$ on $`L_2`$.

To end this section, we give Eq.(7) a second proof using (15). In fact, just notice

$$\rho _1^2\rho _1^2+\rho _4^2\lambda _+1$$

and we have $`\rho _i1,i=1,2,3,4`$. Consequently, $`d_D((0,0),(i,0))`$ has the upper bound $`i`$.

## V Pythagoras’ Theorem for $`d_D(,)`$ on $`L_2`$

Let $`mn1`$ below and we claim that

###### Theorem 1

(Pythagoras’ Theorem on 2D lattice)

$$d_D((0,0),(p,q))=\sqrt{p^2+q^2},pq1$$
(19)
The proof needs only inequalities (17)(18) and we will treat $`q=1`$ and $`q=2,3,4,\mathrm{}`$ respectively. We need an important inequality in mathematics analysis

$$\underset{i=1}{\overset{n}{}}a_i\underset{i=1}{\overset{n}{}}b_in\underset{i=1}{\overset{n}{}}a_ib_i$$
(20)
where the equality holds if $`a_1=a_2=\mathrm{}=a_n`$ or $`b_1=b_2=\mathrm{}=b_n`$.

The general philosophy of the proof can be illustrated as following. First we define the concept of canonical path, which means a subset of $`L_2`$ generated by one point denoted as start and the operations $`T_1^+,T_2^+`$. Then it must have a end if $`T_1^+,T_2^+`$ just act finite times. Second we pick out a rectangle subset $`L(p,q)`$ of $`L_2`$ which is defined as $`\{(m,n)/m=0,1,\mathrm{},p,n=0,1,..,q\}`$, and define the link set on $`L(p,q)`$, $`B(p,q):=\{(m,n,s)/s=1,2;s=1:m=0,1,\mathrm{},p1,n=0,1,\mathrm{},q;s=2:m=0,1,\mathrm{},p,n=0,1,\mathrm{},q1\}`$. For any $`f|_{L(p,q)}`$ defined on $`L(p,q)`$, there is a function $`\mathrm{\Delta }f`$ defined on $`B(p,q)`$ correspondingly, with $`\mathrm{\Delta }f(m,n,1)=(_1^+f)(m,n),\mathrm{\Delta }f(m,n,2)=(_2^+f)(m,n)`$. There is another set induced by $`L(p,q)`$, the set of all canonical paths starting at $`(0,0)`$ and ending at $`(p,q)`$, written as $`\mathrm{\Gamma }(p,q)`$. Any element $`\gamma \mathrm{\Gamma }(p,q)`$ can be considered as $`p+q`$ sequential links, namely a subset of $`B(p,q)`$. As for the first step of proof, we want to show for any $`\gamma \mathrm{\Gamma }(p,q)`$, $`|_{l\gamma }\mathrm{\Delta }f(l)|\sqrt{p^2+q^2}`$, if $`f`$ subjects to (17)(18). Adopting reduction of absurdity, we suppose there is a $`\gamma _0`$ such that $`|_{l\gamma _0}\mathrm{\Delta }f(l)|>\sqrt{p^2+q^2}`$. Then by virtue of the closedness condition (18), $`|_{l\gamma }\mathrm{\Delta }f(l)|>\sqrt{p^2+q^2},\gamma \mathrm{\Gamma }(p,q)`$, which implies

$$\underset{l\gamma }{}|\mathrm{\Delta }f(l)|^2>\frac{p^2+q^2}{p+q}$$

by (20). Therefore

$$\underset{\gamma \mathrm{\Gamma }(p,q)}{}(\underset{l\gamma }{}\mathrm{\Delta }f(l)^2)>R(p,q)\frac{p^2+q^2}{p+q}$$
(21)
where we introduce $`R(p,q)=_{\gamma \mathrm{\Gamma }(p,q)}1=|\mathrm{\Gamma }(p,q)|`$, the number of all canonical paths starting from $`(0,0)`$ and ending $`(p,q)`$, which is equal to the binomial coefficient $`\left(\begin{array}{c}p+q\\ p\end{array}\right)`$. We change the summation for paths in $`\mathrm{\Gamma }(p,q)`$ in (21) to the summation in $`B(p,q)`$,

$$\underset{lB(p,q)}{}(L(l)\mathrm{\Delta }f(l)^2)>R(p,q)\frac{p^2+q^2}{p+q}$$
(22)
where $`L(l)`$ is the weight function on $`B(p,q)`$ introduced by the changing of indice of summations. The geometric interpretation of $`L(l)`$ is the number of canonical paths in $`\mathrm{\Gamma }(p,q)`$ passing through link $`l`$, which is expressed as

$$L(m,n,1)=R(m,n)R(pm1,qn),L(m,n,2)=R(m,n)R(pm,qn1)$$

We induce the third set from $`L(p,q)`$, the set $`Q(p,q)`$ of all fundamental plaques as subsets in $`L(p,q)`$, which can be expressed as $`\{(m,n)/m=0,1,\mathrm{},p1,n=0,1,\mathrm{},q1\}`$. Rewrite (17) as

$$\underset{i=1}{\overset{4}{}}|\mathrm{\Delta }f|^2(p,i)2$$
(23)
which can be considered as a constraint on the four links in one fundamental plaque $`p`$. Now sum (23) for all plaques in $`Q(p,q)`$

$$\underset{pQ(p,q)}{}\underset{i=1}{\overset{4}{}}|\mathrm{\Delta }f|^2(p,i)2|Q(p,q)|=2pq$$

Again we change the summation for plaques to the summation for links and introduce another weight function $`S(l)`$,

$$\underset{lB(p,q)}{}S(l)|\mathrm{\Delta }f|^2(l)2|Q(p,q)|$$
(24)
where $`S(m,0,1)=S(m,q,1)=S(0,n,2)=S(p,n,2)=1,S(other)=2`$, namely to links $`l`$ shared by two plaques $`S(l)=2`$, else $`S(l)=1`$. Using (20) again and noticing $`L(l)0`$, we get

$$2|Q(p,q)|\underset{lB(p,q)}{}S(l)|\mathrm{\Delta }f|^2(l)\underset{lB(p,q)}{}(\frac{S(l)}{L(l)})(L(l)|\mathrm{\Delta }f|^2(l))$$

$$\frac{1}{|B(p,q)|}\underset{lB(p,q)}{}(\frac{S(l)}{L(l)})\underset{lB(p,q)}{}(L(l)|\mathrm{\Delta }f|^2(l))$$

$$\frac{2|Q(p,q)||B(p,q)|}{_{lB(p,q)}\frac{S(l)}{L(l)}}\underset{lB(p,q)}{}(L(l)|\mathrm{\Delta }f|^2(l))$$
(25)
where $`|B(p,q)|=2pq+p+q`$. The contradiction lays between inequalities (22) and (25), if only

###### Lemma 2

$$\frac{2|Q(p,q)||B(p,q)|}{_{lB(p,q)}\frac{S(l)}{L(l)}}R(p,q)\frac{p^2+q^2}{p+q}$$
(26)
As we declared that the description above is just general philosophy of proof, since the operating on $`L(p,q)`$ is complicated, we will just prove $`q=1`$ case in Lemma 2 and leave $`q2`$ cases to a revised $`L(p,q)`$ construction.

As for $`q=1`$, (26)$``$

$$2p(1+3p)(1+p^2)\underset{lB(p,1)}{}\frac{S(l)}{L(l)}$$
(27)
Using the definitions of $`L(l)`$ and $`S(l)`$, we have $`_{lB(p,1)}\frac{S(l)}{L(l)}=2(p+_{n=1}^p\frac{1}{n})`$. Then (27) reduces to show

$$3p^2p^3+(1+p)(\underset{n=1}{\overset{p}{}}\frac{1}{n})$$
(28)
However, (28) is easily checked by $`p=1`$, $`p=2`$, and $`p=3,4,\mathrm{}`$ respectively. Accordingly, we reach the conclusion $`d_D((0,0),(p,1)),p=1,2,3,\mathrm{}`$ has an upper bound $`\sqrt{1+p^2}`$.

We define $`\stackrel{~}{L}(p,q)`$ to be a “folding ruler”, $`\{(m,0),(m,1),(p1,n),(p,n)/m=0,1,..,p,n=2,\mathrm{},q\}`$, for $`pq2`$, and apply our general philosophy of proof to $`\stackrel{~}{L}(p,q)`$. Namely, we induce $`\stackrel{~}{B}(p,q),\stackrel{~}{P}(p,q),\stackrel{~}{Q}(p,q)`$ from $`\stackrel{~}{L}(p,q)`$, and introduce weight functions $`\stackrel{~}{L}(l),\stackrel{~}{S}(l)`$ on $`\stackrel{~}{B}(p,q)`$. If

###### Lemma 3

$$\frac{2|\stackrel{~}{Q}(p,q)||\stackrel{~}{B}(p,q)|}{_{l\stackrel{~}{B}(p,q)}\frac{\stackrel{~}{S}(l)}{\stackrel{~}{L}(l)}}\stackrel{~}{R}(p,q)\frac{p^2+q^2}{p+q}$$
(29)
holds, then we find a contradiction and conclude that $`d_D((0,0),(p,q))`$ has an upper bound $`\sqrt{p^2+q^2}`$ for all $`pq2`$. And here, $`|\stackrel{~}{Q}(p,q)|=p+q1`$, $`|\stackrel{~}{B}(p,q)|=3p+3q2`$, $`\stackrel{~}{R}(p,q)=1+pq`$; $`\stackrel{~}{S}(m,0,2)=\stackrel{~}{S}(p1,n,1)=2,m=1,2,\mathrm{},p1,n=1,\mathrm{},q1,\stackrel{~}{S}(other)=1`$; $`\stackrel{~}{L}(m,0,1)=1+q(pm1),\stackrel{~}{L}(m,0,2)=q,m=0,1,\mathrm{}p1`$; $`\stackrel{~}{L}(p,0,2)=1`$; $`\stackrel{~}{L}(m,1,1)=q(m+1),m=0,1,\mathrm{}p2`$; $`\stackrel{~}{L}(p1,n,2)=p(qn),\stackrel{~}{L}(p,n,2)=1+pn,n=1,\mathrm{},q1`$; $`\stackrel{~}{L}(p1,n,1)=p,n=1,\mathrm{},q`$.

Thus,

$$\underset{l\stackrel{~}{B}(p,q)}{}\frac{\stackrel{~}{S}(l)}{\stackrel{~}{L}(l)}=2+2(\frac{p}{q}+\frac{q}{p})+\delta 6$$

in which

$$\delta =\frac{1}{q}(\underset{n=1}{\overset{p1}{}}\frac{1}{n}1)+\frac{1}{p}(\underset{n=1}{\overset{q1}{}}\frac{1}{n}1)+\underset{n=1}{\overset{p1}{}}\frac{1}{1+nq}+\underset{n=1}{\overset{q1}{}}\frac{1}{1+np}0$$

We claim Lemma 3 is implied by

$$(p+q)|\stackrel{~}{Q}(p,q)||\stackrel{~}{B}(p,q)|3\stackrel{~}{R}(p,q)(p^2+q^2)$$

$$(p+q)(p+q1)(3(p+q)2)3(1+pq)(p^2+q^2),pq2$$
(30)
We prove (30) by induction. First, fixing one $`q=2,3,\mathrm{}`$, we show that (30) holds when $`p=q`$. In this case, (30) is reduced to

$$12q^3+2q3q^4+13q^2,q2$$

The above inequality can be checked to be valid when $`q=2`$, $`q=2`$, and $`q3`$. Second, with this fixed $`q`$, we suppose that the statement is valid in the case $`p=p_0`$, and show that it will also hold in $`p=p_0+1`$. Without misunderstanding, we drop the subscript “$`0`$”. It is sufficient to check

$$9pq+9q^29(q1)p^2+3q^3+7p+4q+3$$

Since $`9pq9(q1)p^2,pq2`$, we just check

$$9q^23q^3+7p+4q+3$$

which is valid when $`q=2`$ and $`q3`$. So the upper bound $`\sqrt{p^2+q^2}`$ follows.

We choose $`f_{(p,q)}(m,n)=\frac{pm+qn}{\sqrt{p^2+q^2}}`$ which can be checked easily to saturate this upper bound. Therefore, theorem 1 follows.

## VI Discussions

Our “natural” Dirac operator is able to be generalized easily into any dimension greater than two. Let $`\mathrm{\Gamma }^i,i=1,2,\mathrm{},2d`$ be the generators of $`Cl(E^{2d})`$ which satisfy

$$\{\mathrm{\Gamma }^i,\mathrm{\Gamma }^j\}=2\delta ^{ij},i,j=1,2,\mathrm{},2d$$
(31)
Define $`\mathrm{\Gamma }_\pm ^k=(\mathrm{\Gamma }^{2k1}\pm i\mathrm{\Gamma }^{2k})/2,k=1,2,\mathrm{},d`$, and (31) changes the form into

$$\{\mathrm{\Gamma }_\pm ^m,\mathrm{\Gamma }_\pm ^n\}=0,\{\mathrm{\Gamma }_\pm ^m,\mathrm{\Gamma }_{}^n\}=\delta ^{mn},m,n=1,2,\mathrm{},d$$

Now let

$$D_{(k)}=\underset{k=1}{\overset{d}{}}\underset{s=\pm }{}\mathrm{\Gamma }_s^k_k^s$$

and one can check the “square root property”

$$(D_{(k)})^2=\underset{k=1}{\overset{d}{}}_k^+_k^{}$$

B. Iochum et al gave a confirmative answer to the question whether there exists a Dirac operator which gives a desired metric on a finite space in . For infinite case, as we have mentioned, the sequential works  together with ours show that this problem is highly nontrivial. Besides, work on the implications of our “natural” Dirac operator on physics is in proceeding.

Acknowledgements
This work was supported by Climb-Up (Pan Deng) Project of Department of Science and Technology in China, Chinese National Science Foundation and Doctoral Programme Foundation of Institution of Higher Education in China. One of the authors J.D. is grateful to Dr. Hua Yang in Princeton for pointing out one mistake in the proof, to Dr. L-G. Jin in Peking University and Dr. H-L. Zhu in Rutgers University for their careful reading on this manuscript.