text:
# 1 Introduction

## 1 Introduction

Born-Infeld (BI) electrodynamics , was proposed in 1934. as a theory in which the energy of electrically charged point particle is finite, in contrast to the Maxwell electrodynamics. The Born-Infeld action is built similarly to the action of relativistic point particle and it introduces dimensional parameter, $`\beta `$, - the “maximal field strength”. It is usually written in one of the following forms

$`S_{BI}=`$ $`\beta ^2{\displaystyle d^4x\left(\sqrt{det\left(g_{\mu \nu }+\frac{1}{\beta }F_{\mu \nu }\right)}\sqrt{detg_{\mu \nu }}\right)}`$ (1)
$`=`$ $`\beta ^2{\displaystyle d^4x\sqrt{g}\left(\sqrt{1+\frac{1}{2\beta ^2}F_{\mu \nu }F^{\mu \nu }\frac{1}{16\beta ^4}(F_{\mu \nu }F_{}^{}{}_{}{}^{\mu \nu })^2}1\right)},`$ (2)
where $``$ is the Hodge-dual. This action has many interesting properties , among them duality symmetry, physical propagation, absence of birefringence, etc.

Actions of the BI-type arise in string/M theory in two main contexts. BI action represents the non-derivative part of the effective open string action. As it was shown in , the bosonic field partition function for the open string in an external field reduces to the BI lagrangian in the string theory limit. On the other hand, BI action is related to D-branes. This comes from the result that the effective action for the open strings ending on D-branes, after the integration of string degrees of freedom , is Dirac-Born-Infeld (DBI) action:

$$S_{DBI}=d^{p+1}x\sqrt{det(\eta _{\mu \nu }+F_{\mu \nu }+_\mu y^i_\nu y^i)},$$
(3)
where $`F_{\mu \nu }`$ is the field strength and $`y^i`$’s are scalar fields. BI action is obtained from (3) for $`y^i=0`$. Conversely, DBI action can be related to BI action in higher dimensions by dimensional reduction.

The generalization of BI electrodynamics to nonabelian gauge theory is not unique. In the general case, if $`F_{\mu \nu }`$ is the field strength of the nonabelian gauge group $`𝒢`$ and $`F_{\mu \nu }=F_{\mu \nu }^aT_a`$ ($`T_a`$ are the generators of $`𝒢`$, $`[T_a,T_b]=if_{abc}T_c`$), the ”determinant” form of the action (1) is not equal to the “square-root” form (2). Different definitions of nonabelian Born-Infeld (NBI) lagrangians are possible, regarding the way of tracing the group indices. The symmetrized trace version of Tseytlin  is often regarded as the one which describes the non-derivative approximation of string theory; however, there are other proposals . Usually NBI lagrangians cannot be put in the closed form in the component fields $`F_{\mu \nu }^a`$.

Following Gal’tsov and Kerner , in this paper we will analyze the simplest version of NBI action in which the trace over the group indices is done under the square-root sign. Gal’tsov and Kerner found particle-like finite energy solutions for the NBI action for the $`SU(2)`$ gauge group. Motivated by this result and by the fact that the dyonic solutions are of interest in the brane theory, we analyze a more general class of solutions. We also discuss the existence of pure monopole solutions.

## 2 Action and field equations

The initial point of our analysis is the following nonabelian Born-Infeld action in Minkowski space:

$$S=\frac{1}{4\pi }d^4x(1),$$
(4)
where $``$ is defined as

$$=\sqrt{1+\frac{1}{2}F_{\mu \nu }^aF^{\mu \nu a}\frac{1}{16}F_{\mu \nu }^aF_{}^{}{}_{}{}^{\mu \nu a}}.$$
(5)
We put that the maximal field strength equals unity, $`\beta =1`$. Lorentz inices $`\mu `$, $`\nu `$ run from 0 to 3 and we will often split them into the temporal part 0 and the spatial part, $`i`$, $`j`$ = 1, 2, 3. The signature which we use is $`(,+,+,+)`$. $`F_{\mu \nu }^a`$ are the field strengths of the $`SU(2)`$ gauge group,

$$F_{\mu \nu }^a=_\mu A_\nu ^a_\nu A_\mu ^a+ϵ^{abc}A_\mu ^bA_\nu ^c,$$
(6)
with $`a`$, $`b`$ = 1, 2, 3. The equations of motion which follow from the NBI action (4) are

$$D_\mu P^{\mu \nu }=0,$$
(7)
where $`P_{\mu \nu }`$ are the ”displacements” defined by

$$P_{\mu \nu }^a=\frac{}{F^{\mu \nu a}}=\frac{F^{\mu \nu a}GF_{}^{}{}_{}{}^{\mu \nu a}}{}.$$
(8)
The quantity $`F^{}`$ denotes the Hodge-dual of $`F`$

$$F_{}^{}{}_{}{}^{\mu \nu }=\frac{1}{2}ϵ^{\mu \nu \rho \sigma }F_{\rho \sigma },$$
(9)
and we use the shorthand notation

$$G=\frac{1}{4}F_{\mu \nu }^aF_{}^{}{}_{}{}^{\mu \nu a}.$$
(10)
The equations of motion (7) can be complemented with the Bianchi identities

$$D_\mu F_{}^{}{}_{}{}^{\mu \nu }=0.$$
(11)
It is important to note that NBI theory has the duality symmetry as BI:

$$F^{\mu \nu }P_{}^{}{}_{}{}^{\mu \nu },P^{\mu \nu }F_{}^{}{}_{}{}^{\mu \nu }.$$
(12)
Duality invariance can be seen from the vacuum equations (7) and (11), too. It can be used to generate new vacuum solutions from the given ones.

The ansatz for the gauge potentials of  was the “monopole ansatz”,

$$A_0^a=0,A_i^a=ϵ_{aik}\frac{1w(r)}{r}\frac{x^k}{r},$$
(13)
and it describes the purely magnetic configurations. Electric and magnetic fields are defined by:

$$E_i^a=F_{i0}^a,B_i^a=\frac{1}{2}ϵ_{ijk}F_{jk}^a.$$
(14)
We will generalize the ansatz (13) – in fact, we will consider the general spherically symmetric static potential of the $`SU(2)`$ group (Witten’s ansatz, ). It is given via three real functions $`a_0(r)`$, $`a_1(r)`$ and $`w(r)`$ of the radial coordinate $`r`$. The components of the gauge potential read:

$$A_0^a=a_0(r)\frac{x^a}{r},$$
(15)
$$A_i^a=a_1(r)\frac{x^ax^i}{r^2}+ϵ_{aik}\frac{1w(r)}{r}\frac{x^k}{r}.$$
(16)
Here $`x^a`$, $`x^i`$ and $`x^k`$ ($`a`$, $`i`$, $`k`$ = 1, 2, 3) are the Cartesian coordinates. The field strengths for this ansatz are

$$E_i^a=a_0^{}\frac{x_ix_a}{r^2}\frac{a_0w}{r}\frac{x_ix_a\delta _{ia}r^2}{r^2},$$
(17)
$$B_i^a=2\delta _{ia}\frac{1w}{r^2}+\frac{(1w)^2}{r^2}\frac{x_ix_a}{r^2}+\left(\frac{1w}{r^2}\right)^{}\frac{x_ix_a\delta _{ia}r^2}{r}+\frac{a_1w}{r^2}ϵ_{iak}x_k$$
(18)
and prime denotes the derivative $`\frac{d}{dr}`$ . The square root $``$ becomes

$$=\sqrt{1+\frac{(1w^2)^2}{r^4}+2\frac{w_{}^{}{}_{}{}^{2}}{r^2}+2\frac{a_1^2w^2}{r^2}2\frac{a_0^2w^2}{r^2}a_{0}^{}{}_{}{}^{2}\frac{[a_0(1w^2)]_{}^{}{}_{}{}^{2}}{r^4}}.$$
(19)
In order to find the equations for $`a_0(r)`$, $`a_1(r)`$ and $`w(r)`$ we can consider the condition of extremality of the action or introduce the ansatz (16) into (7)–(11). After the integration of angular variables, the action is proportional to the lagrangian $`L`$,

$$L=_0^{\mathrm{}}r^2(1)𝑑r.$$
(20)
Varying the unknown functions $`a_0`$, $`a_1`$ and $`w`$, we obtain the set of the equations:

$$w^2a_1=0,$$
(21)
$$(1w^2)\left(\frac{[a_0(1w^2)]^{}}{r^2}\right)^{}=\frac{2w^2a_0}{}\left(\frac{r^2a_0^{}}{}\right)^{},$$
(22)
$$wa_0\left(\frac{[a_0(1w^2)]^{}}{r^2}\right)^{}=\frac{2w(1w^2)}{r^2}\left(\frac{2w^{}}{}\right)^{}\frac{wa_0^2}{}+\frac{wa_1^2}{}.$$
(23)
## 3 NBI dyons

The system of equations (2123) is a complicated nonlinear system. We will search for particular solutions of this system with finite energy. The energy of the static configurations is equal to the negative value of the lagrangian,

$$M=_0^{\mathrm{}}r^2(1)𝑑r.$$
(24)
The convergence of this integral on both boundaries imposes restrictions on the asymptotic behavior of the functions $`a_0`$, $`a_1`$ and $`w`$, which we will discuss later.

Let us first consider the simplest equation, (21): it implies that either $`w(r)=0`$ or $`a_1(r)=0`$. But one can see rather easily that the configuration $`w(r)=0`$, $`a_1(r)0`$ is gauge equivalent to the configuration $`w(r)=0`$, $`a_1(r)=0`$. Indeed, for $`w(r)=0`$ we obtain that (21) and (23) are identically fulfilled, leaving $`a_1(r)`$ undetermined. This means that $`a_1(r)`$ represents the gauge freedom. The value of $`a_1(r)`$ does not influence the values of the field strengths in the case $`w(r)=0`$, as can be seen from (1718). Therefore, we will always assume that $`a_1(r)=0`$ and denote $`a_0(r)=a(r)`$ in the following, keeping the indexed notation like $`a_0`$, $`a_1`$, $`w_0`$ etc. for the coefficients in the power series expansions.

The solutions with $`a(r)=0`$, $`w(r)0`$ were discussed by Gal’tsov and Kerner in detail. In this case, the equations of motion reduce to

$$\left(\frac{w^{}}{}\right)^{}=\frac{w(1w^2)}{r^2},$$
(25)
and the square root $``$ to the expression

$$=\sqrt{1+\frac{(1w^2)^2}{r^4}+2\frac{w_{}^{}{}_{}{}^{2}}{r^2}}.$$
(26)
The simplest solution of (25), $`w(r)=\pm 1`$, is the pure gauge. $`w(r)=0`$ is also a solution, and it has the form of the Dirac monopole: this is embedded $`U(1)`$ monopole. Its energy is finite:

$$M_e=\frac{\pi ^{3/2}}{3\mathrm{\Gamma }(3/4)^2}1.2360.$$
(27)
There is also an infinite discrete set of finite energy solutions $`w_n(r)`$, $`n`$, the so-called Gal’tsov-Kerner glueballs. These solutions can be found numerically using the condition that function $`w(r)`$ with the allowed asymptotic forms at $`r0`$ and $`r\mathrm{}`$ is smooth in the intermediate region. The asymptotic expansions are:

$`r0`$ $`:`$ $`w(r)=1br^2+O(r^4),`$
$`r\mathrm{}`$ $`:`$ $`w(r)=\pm 1+{\displaystyle \frac{c}{r}}+O\left({\displaystyle \frac{1}{r^2}}\right).`$ (28)
Let us note that the solutions behaving at infinity as $`w(r)0`$ are excluded, thus leaving only the configurations with no magnetic charge. Solutions $`w_n(r)`$ behave as magnetic dipoles and have energies which tend to the energy $`M_e`$ of the embedded monopole as $`n\mathrm{}`$ .

The other simple possibility, $`w(r)=0`$, $`a(r)0`$, is also nontrivial. The equations of motion in this case reduce to

$$\left(\frac{a^{}}{r^2}\right)^{}=\left(\frac{r^2a^{}}{}\right)^{},$$
(29)
where now we have

$$=\sqrt{\frac{(1+r^4)(1a_{}^{}{}_{}{}^{2})}{r^4}}.$$
(30)
This equation can be solved explicitly and its solution is a two-parameter family

$$a(r;C,\alpha )=C\pm _0^r\sqrt{\frac{\alpha 1}{\alpha +r^4}}𝑑r,$$
(31)
where $`C`$ and $`\alpha `$ are the integration constants and $`\alpha >1`$. As the energy and the field strenghts do not depend on $`C`$ and the equations are invariant under $`a(r)a(r)`$, we will take $`C=0`$ and the + sign in front of the square root. The explicit form of the solution is given in terms of the elliptic integral ,

$$a(r;\alpha )=\frac{1}{2}(\alpha 1)^{1/2}\alpha ^{1/4}F(\mathrm{arccos}\frac{\sqrt{\alpha }r^2}{\sqrt{\alpha }+r^2},\frac{1}{2}).$$
(32)
The function $`a(r;\alpha )`$ is shown in the Figure 1 for different values of $`\alpha `$. The limiting value of the parameter, $`\alpha =1`$, gives $`a(r)=`$ const, a configuration which is gauge equivalent to the embedded monopole $`w(r)=0`$, $`a(r)=0`$. The energy of the solution (31) is

$$M(\alpha )=\frac{\pi ^{3/2}}{\mathrm{\Gamma }(3/4)^2}\frac{1}{2\alpha ^{1/4}}\left(1\frac{\alpha }{3}\right).$$
(33)
It is unbounded below with the maximum $`M_e`$ at $`\alpha =1`$. We observe that the existence of the electric field decreases the total energy.

We call the solution (31) dyon , as in the asymptotic region, $`r\mathrm{}`$, the behavior of the electric and magnetic fields is given by

$$E_i^a\sqrt{\alpha 1}\frac{x_ix_a}{r^4},B_i^a\frac{x_ix_a}{r^4},$$
(34)
and describes the field strenghts of point-like sources. The “electric charge” of the source is proportional to $`\sqrt{\alpha 1}`$, while the “magnetic charge” is 1.

Let us discuss the duals of the forementioned solutions. One defines the splitting of the displacement tensor in terms of the vectors $`D_i^a`$ and $`H_i^a`$ as

$$P_{i0}^a=D_i^a,P_{ij}^a=ϵ_{ijk}H_k^a.$$
(35)
The duality transformation (12) can then be reexpressed as

$$E_i^aH_i^a=\frac{B_i^aGE_i^a}{},B_i^aD_i^a=\frac{E_i^a+GB_i^a}{}.$$
(36)
In the case of Gal’tsov-Kerner glueballs we have $`E_i^a=0`$, $`G=E_i^aB_i^a=0`$, so the duality transforms

$$E_i^a\frac{B_i^a}{},B_i^a\frac{E_i^a}{}=0.$$
(37)
This means that from the magnetic dipole solution we obtain purely electric solution, which behaves as a dipole since $`1`$ asymptotically.

In the case of a dyon we see that $`Gr^4`$ and $`1`$ at infinity. The leading behavior of the transformed configuration will be

$$B_i^a\frac{x_ix_a}{r^4},E_i^a\sqrt{\alpha 1}\frac{x_ix_a}{r^4},$$
(38)
i. e., the electric and magnetic charges interchange.

## 4 General case

We now turn to the analysis of the general case, $`w(r)0`$, $`a(r)0`$. The first condition that we want to impose on our solutions is finiteness of the energy. This condition restricts the possible behavior of the functions $`w(r)`$, $`a(r)`$ at the boundaries of the integral (24). If we expand $`a(r)`$ and $`w(r)`$ in the power series around $`r=0`$, we conclude that (24) converges if there are no poles in the series, i. e. if they are of the form:

$$a(r)=\underset{0}{\overset{\mathrm{}}{}}a_nr^n,w(r)=\underset{0}{\overset{\mathrm{}}{}}w_nr^n.$$
(39)
When we analyze the other boundary, $`r\mathrm{}`$, we obtain the similar asymptotics:

$$a(r)=\underset{0}{\overset{\mathrm{}}{}}A_nr^n,w(r)=\underset{0}{\overset{\mathrm{}}{}}W_nr^n,$$
(40)
but now the convergence imposes $`W_0A_0=0`$.

In order to analyze the relations among the coefficients in (3940) further, we will assume that the equations of motion are satisfied order by order in $`r`$ (or respectively in $`1/r`$).

Case $`r0`$. From the equation (22) we obtain that $`w_0`$ must be $`\pm 1`$ or $`0`$. As the equations are invariant to the transformation $`w(r)w(r)`$ and to $`a(r)a(r)`$, we will discuss only $`w_0=0`$ and $`w_0=1`$. The similar situation will repeat in the next case.

For $`w_0=0`$ the equation (22) gives: $`w_1=w_2=w_3=\mathrm{}=0`$; the function $`w(r)`$ vanishes. At the same time, from (21) we get $`a_2=a_3=a_4=0`$, $`a_5=\frac{a_1(a_1^21)}{10}`$, $`a_6=a_7=a_8=0`$, $`a_9=\frac{a_1(a_1^21)^2}{24}`$, etc. We also obtain that $`|a_1|<1`$. Thus, both expansions show that this case corresponds to the dyon solution (31) with $`a_1=\sqrt{\frac{\alpha 1}{\alpha }}`$.

For $`w_0=1`$ we get

$`a(r)`$ $`=`$ $`a_1r+a_3r^3+O(r^5),`$
$`w(r)`$ $`=`$ $`1+w_2r^2+w_4r^4+O(r^5),`$ (41)
where

$$a_3=\frac{8a_1^3w_2+8a_1w_2^32a_1w_2}{10a_1^220w_2^25},$$

$$w_4=\frac{6w_2^2+a_1^4(2+20w_2^2)+16w_2^4(7+22w_2^2)a_1^2(1+42w_2^2+408w_2^4)}{20(1a_1^2+4w_2^2)(12a_1^2+4w_2^2)},$$

etc., and $`|a_1|<\frac{1}{\sqrt{3}}`$. We will analyze this asymptotics in the following, let us just note here that for $`a_1=0`$ it is the one obtained in .

Case $`r\mathrm{}`$. We consider separately the possibilities $`W_0=0`$ and $`A_0=0`$.

If $`W_0=0`$, the assumption that the equations of motion are satisfied order by order in $`1/r`$ leads to $`W_1=W_2=W_3=\mathrm{}=0`$. For coefficients of $`a(r)`$ we get $`A_2=A_3=A_4=0`$ , $`A_5=\frac{A_1(A_1^2+1)}{10}`$ , $`A_6=A_7=A_8=0`$, $`A_9=\frac{A_1(A_1^2+1)^2}{24}`$, etc. Again, we obtain the power series expansion of the dyon (31), in this case around infinity.

For the second possibility, $`A_0=0`$, the solutions behave asymptotically as

$`a(r)`$ $`=`$ $`{\displaystyle \frac{A_2}{r^2}}+{\displaystyle \frac{A_3}{r^3}}+{\displaystyle \frac{A_4}{r^4}}+O\left({\displaystyle \frac{1}{r^5}}\right)`$
$`w(r)`$ $`=`$ $`1+{\displaystyle \frac{W_1}{r}}+{\displaystyle \frac{W_2}{r^2}}+{\displaystyle \frac{W_3}{r^3}}+{\displaystyle \frac{W_4}{r^4}}+O\left({\displaystyle \frac{1}{r^5}}\right),`$ (42)
where the following relations are fulfilled

$$A_3=A_2W_1,A_4=\frac{18A_2W_1^2A_2^3}{20},W_2=\frac{6W_1^2A_2^2}{8},$$

$$W_3=\frac{22W_1^39A_2^2W_1}{40},W_4=\frac{17A_2^4540A_2^2W_1^2+772W_1^4}{1920}.$$

From this analysis we see that, in order to find new solutions, we need to join the asymptotics (41) and (42) smoothly. Our first attempt was to do the numerical integration from $`r=0`$ to the right or from $`r=\mathrm{}`$ to the left, with the initial conditions defined appropriately. Doing this, we obtain the generic solution of a typical form shown in the Figure 2. The coupling of $`a(r)`$ and $`w(r)`$ induces the oscillations of $`w(r)`$ which reduce its initial value 1 for $`r=0`$ to 0 for $`r=\mathrm{}`$.

This solution is interesting, as it has the behavior of ’t Hooft-Polyakov monopole. However, it is numerically unstable: if we keep the same values of $`w_2`$, $`a_1`$ and $`r_i`$ but decrease the integration step $`h`$, we obtain the functions given in the Figure 3.

The oscillations of $`w(r)`$ increase to a larger region of $`r`$, while the asymptotic value of $`a(r)`$ changes. We conclude that the solutions of this type are nonanalytic. Indeed, from the previous discussion of asymptotics we know that $`w(r)0`$ as $`r\mathrm{}`$ is compatible only with $`w(r)=0`$. Further numerical analysis of energy confirms this conclusion: the values of energy differ for orders of magnitude for different integration steps and therefore signal that the energy diverges. We see that in the NBI case, as in the pure Yang-Mills theory, $`w(r)=0`$ and $`w(r)=1`$ are separated by an infinite energy barrier and it is impossible to find the solution of finite energy which interpolates between them.

The second possibility for numerical investigation is to start the integration from both sides $`r=0`$ and $`r=\mathrm{}`$ with the given asymptotics, and try to join smoothly the solution in the intermediate region by varying the parameters $`w_2`$, $`a_1`$, $`W_1`$ and $`A_2`$. A numerical programme which handles this type of boundary conditions  was made, and proved to be correct and very efficient in the simple case of small $`a_1`$, $`A_2`$ (glueballs). However, no new solutions were found using this programme for a wide range of initial parameters. This might be a consequence of some weaknesses of the implemented variational procedure (Newton-Raphson), due to the high dimensionality of the parameter space. We are, however, inclined to interpret this as a strong numerical evidence that there are no further finite-energy solutions of the system (2123).

## 5 Conclusions

The set of equations (2123), which represent the equations of motion for the static spherically symmetric configurations of $`SU(2)`$ NBI action (4), is analyzed. The asymptotic analysis shows that, if one imposes finiteness of energy, there are only three possible types of solutions: glueballs, dyons and solutions of the form (39)–(41).

Dyon solutions are of importance in the brane-theory, as they represent strings ending on three-brane . The name dyon, introduced after , is used in the generalized sense: there is no Higgs field to determine the unbroken $`U(1)`$ group. As in the case of Julia-Zee dyon, the electric charge of this solution is continuous while the magnetic charge is 1. However, the hope that the components $`A_0^a`$ of the vector potential (given via the function $`a(r)`$) can, through the nonlinear interaction, take the role of Higgs and counterbalance the magnetic field to produce the monopole of the ’t Hooft-Polyakov type failed. Instead of the exponential decay, $`a(r)`$ induces the oscillations of $`w(r)`$ with infinite energy. This could be expected from the fact that the change of the action from Yang-Mills to NBI does not change the topology of the fields which are included, necessary for the existence of monopole . The solutions of the NBI models with Higgs fields were discussed in .

Finally, let us add that, although the solutions of the third mentioned type are allowed by the energy considerations, we have a strong numerical indication that they do not exist. This problem might deserve further numerical analysis.

Acknowledgments. M. B. wishes to thank Prof. F. W. Hehl for his hospitality during the stay at the Institute of Theoretical Physics of the University of Cologne, when this work was initiated. We also want to thank Profs. A. Belić and D. Galtsov for the useful comments concerning the numerical integration. The work was supported in part by the DAAD grant A/00/17208. The numerics was done at the IPCF of the Institute of Physics, Belgrade.


format: text

text:
# Direct Monte Carlo simulations of accretion discs

## 1. Introduction

The Direct Simulation Monte Carlo method was originally developed to treat rarefied gas flow (see Bird 1994; Matsuda, Mizutani & Boffin 2001). In DSMC, microscopic molecular motions and their mutual collisions are calculated. Macroscopic quantities such as density, velocity and temperature of fluids are calculated by taking some mean of microscopic quantities. In DSMC, molecular collisions are handled in a stochastic manner, hence the Monte Carlo denomination of the method. The collision model in DSMC is such that the molecules are assumed to be hard billiard balls. Since we solve the molecular motions, molecular viscosity and molecular thermal conduction are automatically included in the system. This ’physical’ viscosity may mimic a turbulent viscosity ($`\alpha `$ viscosity) assumed in the standard theory of accretion disk.

## 2. Application of DSMC to a three-dimensional accretion disc

We consider a binary system with mass ratio $`q=0.5`$. We work in the rotational frame. We restrict our computational region to a square, the center of which is the primary and just including the L1 point from which gas particles are injected constantly into the computational region. The region is divided into $`100\times 100\times 20`$ cells. The cell size is such that $`dx=dy=2dz`$. We start injection of gas particles at a constant rate per time step and follow the time evolution of the system. The number of particles in the computational region increases with time and reaches a steady state. After a steady state is reached, we further calculate to take a time average. Typically we computed up to t=60-100 in dimensionless unit, that is about 10-15 orbital periods.

In accretion discs, radiative cooling must occur, and some mechanism of cooling is essential in our calculation. We here consider a simple model of cooling: we reduce the thermal speed of the molecules by a factor $`f=0.97`$ at each time step.

## 3. Results

Figure 1 displays the time averaged iso-density surface, on which temperature distribution is pasted, and velocity vectors in the rotational plane of our simulation. Two cases of injected particles number, $`n_{inj}=6`$ and $`n_{inj}=30`$, are compared. We can observe the gas stream from the L1 point, the accretion disc and spiral-shaped shocks. In the case of $`n_{inj}=6`$, we find that the accretion disc is very circular. The L1 stream hitting at the disc is deflected from its direction downward by the disc flow. Because of the collision between the L1 stream and the disc flow, the iso-density surface is elevated and the temperature on the elevated area is raised (a hot line or heated wall). In the case of $`n_{inj}=30`$, the L1 stream is so strong that the direction of the disc flow is deflected and the disc shape is elongated. We may observe a pair of spiral shocks in this case.

## References

Bird, G.A. 1994, Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Oxford Science Publications)

Matsuda, T., Mizutani, H. & Boffin, H.M.J. 2001, in IAU Symp. 208, Astrophysical Supercomputing using Particle Simulations, ed. J. Makino & P. Hut.


text:
# Dynamics of Quantum Phase Transition in an Array of Josephson Junctions

## Abstract

We study the dynamics of the Mott insulator-superfluid quantum phase transition in a periodic 1D array of Josephson junctions. We show that crossing the critical point at a finite rate with a quench time $`\tau _Q`$ induces finite quantum fluctuations of the current around the loop proportional to $`\tau _Q^{1/6}`$. This scaling could be experimentally verified with in array of weakly coupled Bose-Einstein condensates or superconducting grains.

Thermal fluctuations critically slow down close to a continuous symmetry breaking thermodynamic phase transition. Both the relaxation time $`\tau `$ and the correlation length $`\xi `$ diverge when a dimensionless parameter $`ϵ`$, which measures the distance from the critical point, tends to zero: $`\xi |ϵ|^\nu `$ and $`\tau |ϵ|^{z\nu }`$, with $`z`$ and $`\nu `$ critical exponents. In the infinitesimally slow transition the symmetry broken phase is entered in a state of equilibrium with a fully ordered complex order parameter. However, in a quench occuring at a finite rate transition the critical slowing down ($`\tau \mathrm{}`$ when $`ϵ0`$) implies that the system goes out of equilibrium some time before the transition. As a result the complex order parameter which emerges after the transition will assume different random phases in different domains of space. This leads to formation of topological defects such as vortices wherever the circulation of the phase around a closed loop happens to be nonzero. This process is known as Kibble-Zurek mechanism (KZM) . The size of the domains follows from a simple argument : Close to the transition one can linearize the time dependence of $`ϵ`$ as $`ϵ(t)t/\tau _Q`$. The transition rate is $`r(t)=\dot{ϵ}/ϵ=1/|t|`$. The system goes out of equilibrium at the time $`\widehat{t}`$ before the transition when $`r(\widehat{t})=\tau ^1(\widehat{t})`$. After that time the state of the system essentially does not change until $`+\widehat{t}`$, when the rate $`r(t)`$ becomes again equal to the relaxation rate, $`r(+\widehat{t})=\tau ^1(+\widehat{t})`$. At $`\widehat{t}=\tau _Q^{z\nu /(1+z\nu )}`$ fluctuations of the order parameter with wavelengths longer than $`\widehat{\xi }=\xi (\widehat{t})=\tau _Q^{\nu /(1+z\nu )}`$ begin to grow exponentially while short wavelength fluctuations remain unchanged. The size of the correlated domains is given by $`\widehat{\xi }`$ and the density of vortices by $`1/\widehat{\xi }^2`$.

Theoretical and experimental studies of KZM have concentrated so far on thermal continuous phase transitions with the dynamics of the order parameter governed by an effective irreversible time-dependent Ginzurg-Landau theory. Most attention has been devoted to the normal-superfluid transition in $`{}_{}{}^{3}He`$ , superconductors , and, more recently, dilute Bose-Einstein condensates . To date, truly microscopic quantum approaches have been too complicated to extract useful predictions.

A quench-induced quantum phase transition (QPT) at temperature $`T=0`$ must be treated in a microscopic way. It is a common wisdom that some properties of a quantum transition can be obtained by an exact map from a thermodynamic transition : the correlation length in the ground state of the quantum system scales like $`\xi |ϵ|^\nu `$ and the gap between the ground state and the first excited state like $`dE|ϵ|^{z\nu }`$. To study the dynamics of a quantum transition, which drives the system out of its ground state, we also need information about correlations in excited states, which is not provided by this map. As we will see below, the essence of the KZM, which is the competition between the transition rate and the timescale on which the system can react, is applicable to quantum transitions. However, the quantum scenario and in particular interpretation of its results differ from the thermal case. The main reason is the reversibility of the quantum dynamics, as opposed to the dynamical irreversibility of the thermal critical dynamics.

In this Letter, we study the appearance of a nonzero current while the system is undergoing a (zero-temperature) quantum phase transition. We develop a microscopic, dynamical theory and suggest experiments to test our predictions.

Josephson junction arrays. A prototype system displaying a continuous QPT is an array of mesoscopic Josephson junctions (JJ). Superconducting JJ networks are now available thanks to important advances in techniques for ultra-small superconducting grains although existence of gauge fields may complicate the analysis . On the other hand, manipulation of optically trapped Bose-Einstein condensates promises the observation of QPT in neutral quantum fluids . In such systems, a superconducting/superfluid-insulator phase transition is driven by the competition between two physical magnitudes: the Josephson coupling energy $`E_J`$, which governs the tunneling through the intra-well barriers, and the on-site interparticle interaction energy $`E_C`$. When they are comparable, there is a competition between long-range order (which is favoured by $`E_J`$) and localization (induced by $`E_C`$). When the latter prevails, no net current can flow through the junctions.

An effective Hamiltonian describing the quantum dynamics of an array of JJ’s is given by the Quantum Phase Model (QPM) :

$$\widehat{H}=\underset{l=1}{\overset{N_s}{}}\frac{E_C}{2}\widehat{n}_l^2\underset{k,l}{}E_J\mathrm{cos}(\widehat{\varphi }_k\widehat{\varphi }_l)$$
(1)
with the last sum running over the nearest-neighbour sites. For our purposes, it is enough to consider the case of a diagonal charging energy matrix $`E_c`$, independent of the site indices. The phase $`\widehat{\varphi }_l`$ and the number of atoms $`\widehat{n}_l`$ in each site of the junction are (with some important caveats ) non-commuting conjugate observables $`[\widehat{n}_l,\widehat{\varphi }_l]=i`$, and, in the $`\varphi `$-representation, $`\widehat{n}_lni\frac{}{\varphi _l};\widehat{\varphi }_l\varphi _l`$. Therefore, $`\varphi _l`$ and $`n_l`$ play the role of coordinate and conjugate momentum, and satisfy the Heisenberg uncertainty relation The insulator phase is characterized by large quantum phase fluctuations in the ground state, which destroy the long-range order among sites. The 1D Hamiltonian Eq.(1) exhibits a continuous Mott phase transition at the critical value $`G_cE_J/E_C=0.617`$ . In what follows we study the microscopic dynamics of quantum phase transitions in a one-dimensional chain with periodic boundary conditions.

The dynamics is governed by the QPM Hamiltonian Eq.(1), and satisfies a Schrödinger equation

$$i\frac{}{t}\mathrm{\Psi }=\frac{1}{2}\underset{l=1}{\overset{N_s}{}}\frac{^2}{\varphi _l^2}\mathrm{\Psi }G\underset{k,l}{}\mathrm{cos}(\varphi _k\varphi _l)\mathrm{\Psi }.$$
(2)
with $`G=E_J/E_c`$ and the time rescaled as $`tE_ct`$. This equation is valid for a large number of atoms per site, $`n1`$, and for $`Gn^2`$ . The ground state below the transition ($`G<G_c`$) is close to the Fock state $`|n,n,n,\mathrm{}`$ with all sites occupied by the same number of atoms $`n`$. In the phase representation this state is described by a uniform wavefunction, $`\mathrm{\Psi }(\varphi _l)=\mathrm{const}`$, where all phase differences between nearest neighbor sites have the maximal dispersion of $`\mathrm{\Delta }\varphi 1`$. Above the transition $`(G>G_c)`$ the ground state is a number squeezed state which continuously tends to a Fock state when $`GG_c^+`$, and to a coherent state for $`Gn^2`$. Therefore, when $`GG_c`$ one can describe the low energy part of the spectrum of Eq.(2) in a harmonic approximation

$$i\frac{}{t}\mathrm{\Psi }=\frac{1}{2}\underset{l=1}{\overset{N_s}{}}\frac{^2}{\varphi _l^2}\mathrm{\Psi }+\frac{G}{2}\underset{k,l}{}(\varphi _k\varphi _l)^2\mathrm{\Psi }.$$
(3)
This equation is diagonalized by normal modes numbered by momentum $`\mu \{N_s+1,\mathrm{},+N_s\}`$, $`\mathrm{\Psi }=_\mu \mathrm{\Psi }_\mu (\mathrm{\Phi }_\mu )`$. There is one zero mode $`\mathrm{\Phi }_0_l\varphi _l`$. All other modes $`\mathrm{\Phi }_\mu `$ have nonzero frequencies $`\sqrt{\gamma _\mu G}`$\[$`G^{1/2}(\frac{2\pi |\mu |}{N_s})`$ for small $`\mu `$\], which scale as $`G^{1/2}`$,

$$i\frac{}{t}\mathrm{\Psi }_\mu =\frac{1}{2}\frac{^2}{\mathrm{\Phi }_\mu ^2}\mathrm{\Psi }_\mu +\frac{\gamma _\mu G}{2}\mathrm{\Phi }_\mu ^2\mathrm{\Psi }_\mu .$$
(4)
To implement quench with a quench timescale $`\tau _Q`$ we linearly ramp the control parameter $`G`$ in Eq.(2) as

$$G(t)=\frac{t}{\tau _Q}.$$
(5)
The time $`t`$ runs from $`0`$ to $`\tau _QG_{\mathrm{max}}`$ with $`G_{\mathrm{max}}G_c`$.

Adiabatic transition. In the limit $`\tau _Q\mathrm{}`$ the transition is adiabatic and the system state adiabatically follows its ground state from the Fock state $`|n,\mathrm{},n`$ at $`G=0`$ to a coherent state with all the nearest neighbor phase differences close to zero, $`\mathrm{\Delta }\varphi 0`$, for $`G_{\mathrm{max}}G_c`$. The ground state at $`G_{\mathrm{max}}G_c`$ is well described by the ground state of the harmonic oscillators (4). The dispersion of phases in any harmonic oscillator ground state at $`G_{\mathrm{max}}`$ is proportional to $`G_{\mathrm{max}}^{1/4}`$. The final dispersion of phase differences in the adiabatic transition is

$$\mathrm{\Delta }\varphi _{\tau _Q\mathrm{}}G_{\mathrm{max}}^{1/4}.$$
(6)
Instantaneous transition. In the opposite limit of an istantaneous transition, when $`\tau _Q0`$, the system is not able to adjust its quantum state to the changing hopping rate $`G(t)`$, and it remains in the initial Fock state till $`G_{\mathrm{max}}`$. The final dispersion of the phase differences in the uniform Fock state, $`\mathrm{\Psi }(\varphi _l)=\mathrm{const}`$, is

$$\mathrm{\Delta }\varphi _{\tau _Q0}\mathrm{\hspace{0.33em}1}.$$
(7)
For a generic quench the evolution of the system is approximately adiabatic when $`G`$ is far from $`G_c`$, and impulse when $`G`$ is close to $`G_c`$; compare Fig.1 for a 3-site periodic lattice. This is similar to the thermal case . As $`G`$ increases starting from $`0`$ the system initially follows its ground state, which is an incoherent superposition of Fock states, until a $`\widehat{G}_{}<G_c`$ when the transition ceases to be adiabatic. Between $`\widehat{G}_{}`$ and a certain $`\widehat{G}>G_c`$ the evolution is impulse: to a first approximation the state of the system does not change. The system arrives at $`\widehat{G}`$ in an incoherent superposition of Fock states with a dispersion $`\mathrm{\Delta }\varphi 1`$ as in the initial Fock state. After $`\widehat{G}`$ the transition again becomes adiabatic.

Fast transitions with $`(E_c)\tau _Q1`$. Let us estimate $`\widehat{G}`$ for a dimensionless $`\tau _Q1`$ when we expect that $`\widehat{G}G_c`$ and we can use the harmonic approximation (3,4) for $`G\widehat{G}`$. For the linear quench (5) the periods of the oscillators decay like $`G(t)^{1/2}=\sqrt{\tau _Q/tG_{\mathrm{max}}}`$. The transition timescale is $`G/\dot{G}=t`$. This rate becomes comparable to the periods of the oscillators at the time $`\widehat{t}`$ when $`G^{1/2}G/\dot{G}`$: compare Fig.1 for a 3-site periodic lattice. A solution of this equation gives

$`\widehat{t}_{\tau _Q1}\tau _Q^{1/3}(\mathrm{or}E_c^{2/3}\tau _Q^{1/3}),`$ (8)
$`\widehat{G}_{\tau _Q1}G(\widehat{t}_{\tau _Q1})\tau _Q^{2/3}.`$ (9)
As anticipated, when $`\tau _Q1`$ we have $`\widehat{G}_{\tau _Q1}G_c`$.

We can use the fact that $`\widehat{G}_{\tau _Q1}`$ to estimate the final dispersion $`\mathrm{\Delta }\varphi `$ at $`G_{\mathrm{max}}`$. At $`\widehat{G}_{\tau _Q1}`$ the dispersion is $`\mathrm{\Delta }\varphi 1`$. For $`G>\widehat{G}_{\tau _Q1}`$ the evolution of the harmonic oscillators is adiabatic. There is no mixing between the eigenstates of any oscillator $`\mu `$ (4). The dispersion of the phase $`\mathrm{\Phi }_\mu `$ in any given eigenstate of the oscillator (4) scales as $`G^{1/4}`$. In the adiabatic evolution for $`G(t)>\widehat{G}_{\tau _Q1}`$ the dispersion of $`\mathrm{\Phi }_\mu `$ shrinks like $`[\widehat{G}_{\tau _Q1}/G(t)]^{1/4}`$. $`\mathrm{\Delta }\varphi `$ shrinks in the same way. The final dispersion at $`G_{\mathrm{max}}`$ is $`\mathrm{\Delta }\varphi [\widehat{G}_{\tau _Q1}/G_{\mathrm{max}}]^{1/4}`$ or

$$\mathrm{\Delta }\varphi _{\tau _Q1}\tau _Q^{1/6}G_{\mathrm{max}}^{1/4}(\mathrm{or}E_c^{1/6}\tau _Q^{1/6}G_{\mathrm{max}}^{1/4}).$$
(10)
In a 1D periodic lattice this $`\mathrm{\Delta }\varphi _{\tau _Q1}`$ translates into a dispersion of the angular momentum $`L=i_l(a_l^{}a_{l+1}a_{l+1}^{}a_l)`$ through the formula

$$\mathrm{\Delta }L\sqrt{N_s}n\mathrm{\Delta }\varphi $$
(11)
This $`\mathrm{\Delta }L`$ may translate into dispersion of the winding number $`\sqrt{N_s}\mathrm{\Delta }\varphi `$ when the atoms after the quench are forced to condense. Eq.(10) is valid when the predicted $`\widehat{G}_{\tau _Q1}\tau _Q^{2/3}`$ is much less than the final $`G_{\mathrm{max}}`$. This condition is satisfied when $`\tau _QG_{\mathrm{max}}^{2/3}`$. In fact for the slowest $`\mu =\pm 1`$ modes to become adiabatic before the maximal possible $`G_{\mathrm{max}}=n^2`$ we need $`\tau _Qn/N_s`$.

Figures 2 and 3 show results of numerical simulations of the 3-site periodic lattice with a total of $`300`$ atoms. In Fig.2 we show that $`\mathrm{\Delta }\varphi `$ remains the same as in the initial Fock state for $`G<\widehat{G}_{\tau _Q1}`$. In Figure 3 we verify the scaling $`\widehat{t}\tau _Q^{1/3}`$. We also performed similar simulations for a total of $`120`$ (commensurate) and $`120+1`$ atoms (non-commensurate density of atoms). Plots like those in Fig.2,3 for $`120`$ and $`121`$ atoms are impossible to distinguish which shows that an extra noncommensurate atom does not make any difference for fast transitions.

Slow transitions with $`(E_c)\tau _Q1`$. So far we concentrated on $`\tau _Q1`$ with $`\widehat{G}`$ well in the harmonic regime. In slower transitions $`\widehat{G}`$ is close to $`G_c`$ where the harmonic approximation cannot be applied. The slow transitions ($`\tau _Q1`$) are probing the critical behaviour close to $`G_c`$.

In the perioric $`3`$-site array there is no phase transition but, rather, a crossover: the gap $`dE`$ is minimal at the crossover point $`G_c0.5`$ but it does not vanish there, compare Fig.1. For $`\tau _Q1`$ the transition is adiabatic and the system follows its ground state. The phase difference dispersion at $`G_{\mathrm{max}}G_c`$ is the dispersion in the ground state of the harmonic oscillators (6). A thermal crossover transition was studied in .

In a large $`2`$D array with a commensurate density of atoms there is a phase transition and not a crossover . The gap $`dE`$ vanishes at $`G=G_c`$. Close to the critical point, $`GG_c`$, the energy gap scales as $`dE(GG_c)^\nu `$, where $`\nu =2/3`$ is an exact renormalization group critical exponent. Note that in the harmonic regime for $`GG_c`$ we have $`dEG^{1/2}`$ which is consistent with a mean-field value of $`\nu =1/2`$. With a linearized $`G(t)G_c=\frac{t}{\tau _Q}`$ the gap $`\mathrm{\Delta }E(GG_c)^{2/3}=(t/\tau _Q)^{2/3}`$ becomes equal to the quench rate $`r(t)=\dot{G}/(GG_c)=1/t`$ at $`\widehat{G}G_c+1/\tau _Q^{3/5}`$. Before $`\widehat{G}`$ the system is in an incoherent superposition of Fock states. After $`\widehat{G}`$ the evolution is adiabatic, and with increasing $`G`$ the phase dispersion shrinks together with the phase widths of the system eigenstates. The dispersion measured at $`G_{\mathrm{max}}G_c`$, where we can use the harmonic approximation, is

$$\mathrm{\Delta }\varphi _{\tau _Q1}^{2\mathrm{D}}\frac{\widehat{G}^{1/4}}{G_{\mathrm{max}}^{1/4}}\frac{(G_c+\frac{1}{(E_c^{3/5})\tau _Q^{3/5}})^{1/4}}{G_{\mathrm{max}}^{1/4}}.$$
(12)
This formula is consistent with Eq.(6) because $`G_c=𝒪(1)`$. Here we again use the fact that the phase dispersion of harmonic oscillator eigenstates shrinks like $`G^{1/4}`$. The critical behaviour is realized in the $`(E_c^{3/5})\tau _Q^{3/5}`$ term.

Concluding remarks. The state of the system after the transition does not have definite angular momentum or definite phase differences between lattice sites. That would be the case in the thermal KZM, where $`\mathrm{\Delta }\varphi `$ would describe the dispersion in an ensemble of different possible classical outcomes. In a quantum phase transition at zero temperature the state of the system after the transition is in a coherent superposition of states with different $`\mathrm{\Delta }\varphi `$. Either a measurement or decoherence  are needed to convert this coherent superposition into a mixtures of states, each with definite current.

Another difference with respect to the thermal KZM is due to the reversibility of the quantum dynamics. In the thermal case the characteristic lenghtscale $`\widehat{\xi }`$ is frozen after the symmetry breaking transition is completed. $`\widehat{\xi }`$ is a permanent record of the transition rate $`\tau _Q`$: manipulations with $`ϵ`$ do not change the winding number as long as the system remains in the symmetry broken phase with $`ϵ>0`$. In contrast, even after the diabatic quantum transition is completed one can change $`\mathrm{\Delta }\varphi _{\tau _Q}`$. Adiabatic variations of $`G`$ away from $`G_{\mathrm{max}}`$ are accompanied by changes in $`\mathrm{\Delta }\varphi `$, $`\mathrm{\Delta }\varphi =\mathrm{\Delta }\varphi _{\tau _Q}(G_{\mathrm{max}}/G)^{1/4}`$.

In conclusion, we have predicted the phase dispersion after a diabatic insulator-superfluid quantum phase transition in an array of Josephson junctions. This theory is a quantum counterpart of the Kibble-Zurek mechanism for topological defect formation in classical thermal phase transitions. Our predictions can be tested experimentally in 1D or 2D superconducting JJ arrays. The possibility of tuning the Josephson coupling energy in a superconducting JJ has been demonstrated recently . 1D (or 3D ) JJ arrays realized with Bose-Einstein condensates can be, perhaps even more easily, accurately tailored . For the case of a 1D array our Eq.(10) predicts a dispersion of the phases $`\mathrm{\Delta }\varphi G_{\mathrm{max}}^{1/4}\tau _Q^{1/6}`$ that can be directly measured by the interference techniques which have been used in the experiments .

We acknowledge useful discussions with James Anglin and Diego Dalvit. This work was partially supported by the DOE and the Cofinanziamento MURST.


text:
# I Introduction

## I Introduction

The problem of the cosmological constant and the vacuum energy associated with it is of high interest these days. There are many questions related to it at the quantum level, all of which are related to quantum gravity. Why there exists the critical mass density and why the cosmological constant has this value? Trying to answer these questions and others were recently the subject of many publications \[1-18\].

In this paper it is shown that the cosmological relativity theory  predicts the value $`\mathrm{\Lambda }=2.28\times 10^{35}`$s<sup>-2</sup> for the cosmological constant. This value of $`\mathrm{\Lambda }`$ is in excellent agreement with the measurements recently obtained by the High-Z Supernova Team and the Supernova Cosmological Project \[20-26\].

## II The Cosmological Constant

In 1922 Friedmann solved the Einstein gravitational field equations and obtained nonstatic cosmological solutions presenting an expanding universe . Einstein, who thought at that time that the universe should be static and unchanged forever, suggested a modification to his original field equations by adding to them the so-called cosmological term which can stop the expansion. The field equations with the added term are:

$$R_{\mu \nu }\frac{1}{2}g_{\mu \nu }R+\mathrm{\Lambda }g_{\mu \nu }=\kappa T_{\mu \nu },$$
$`(1),`$
where $`\mathrm{\Lambda }`$ is the cosmological constant, the value of which is supposed to be determined by experiment. In Eq. (1) $`R_{\mu \nu }`$ and $`R`$ are the Ricci tensor and scalar, respectively, $`\kappa =8\pi G`$, where $`G`$ is Newton’s constant and the speed of light is taken as unity.

Soon after that Hubble  found experimentally that the distant galaxies are receding from us, and the farther the galaxy the bigger its velocity as determined by its redshift.

After Hubble’s discovery that the universe is expanding, the role of the cosmological constant to allow static homogeneous solutions to Einstein’s equations in the presence of matter, was looked upon as unnecessary. For a long time the cosmological term was considered to be of no interest in cosmological physical problems.

## III The Friedmann Universe

For a homogeneous and isotropic universe with the line element

$$ds^2=dt^2a^2(t)R_0^2\left[\frac{dr^2}{1kr^2}+r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2\right)\right],$$
$`(2)`$
where $`k`$ is the curvature parameter ($`k=1,0,1`$) and $`a(t)=R(t)/R_0`$ is the scale factor, with the energy-momentum tensor

$$T_{\mu \nu }=(\rho +p)u_\mu u_\nu pg_{\mu \nu },$$
$`(3)`$
Einstein’s equations (1) reduce to the two Friedmann equations

$$H^2\left(\frac{\dot{a}}{a}\right)^2=\frac{\kappa }{3}\rho +\frac{\mathrm{\Lambda }}{3}\frac{k}{a^2R_0^2},$$
$`(4)`$
$$\frac{\ddot{a}}{a}=\frac{\kappa }{6}\left(\rho +3p\right)+\frac{\mathrm{\Lambda }}{3}.$$
$`(5)`$
In Eqs. (4) and (5) $`H`$ is Hubble’s parameter, $`\rho `$ is the mass density and $`p`$ is the pressure. These equations admit a static solution ($`\dot{a}=0`$) with $`k>0`$ and $`\mathrm{\Lambda }>0`$.

From the Friedmann equation (4) it then follows that for any value of the Hubble parameter $`H`$ there exists a critical mass density $`\rho _c=3H_0^2/\kappa `$ at which the spatial geometry is flat ($`k=0`$). One usually measures the total mass density in terms of the critical density $`\rho _c`$ by means of the density parameter $`\mathrm{\Omega }=\rho /\rho _c`$.

In general, the mass density $`\rho `$ includes contributions from various distinct components. From the point of view of cosmology, the relevant aspect of each component is how its contribution to the total energy density evolves as the universe expands. A positive $`\mathrm{\Lambda }`$ causes acceleration to the universe expansion, whereas a negative $`\mathrm{\Lambda }`$ and ordinary matter tend to decelerate it. Moreover, the relative contributions of the components to the energy density change with time. For $`\mathrm{\Omega }_\mathrm{\Lambda }<0`$, the universe will always recollapse to a Big Crunch. For $`\mathrm{\Omega }_\mathrm{\Lambda }>0`$ the universe will expand forever unless there is sufficient matter to cause recollapse before $`\mathrm{\Omega }_\mathrm{\Lambda }`$ becomes dynamically important. For $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$ we have the familiar situation in which $`0<\mathrm{\Omega }_M1`$ universes expand forever and $`\mathrm{\Omega }_M>1`$ universes recollapse. (For more details see Ref. 19.)

## IV The Supernovae Experiments

Recently two groups (the Supernovae Cosmology Project and the High-Z Supernova Team) presented evidence that the expansion of the universe is accelerating \[20-26\]. These teams have measured the distances to cosmological supernovae by using the fact that the intrinsic luminosity of Type Ia supernovae is closely correlated to their decline rate from maximum brightness, which can be independently measured. These measurements, combined with redshift data for the supernovae, led to the prediction of an accelerating universe. Both teams obtained

$$\mathrm{\Omega }_M0.3,\mathrm{\Omega }_\mathrm{\Lambda }0.7,$$
$`(6)`$
and strongly ruled out the traditional ($`\mathrm{\Omega }_M`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$)=(1, 0) universe. This value of the density parameter $`\mathrm{\Omega }_\mathrm{\Lambda }`$ corresponds to a cosmological constant that is small but, nevertheless, nonzero and positive,

$$\mathrm{\Lambda }10^{52}\text{m}^210^{35}\text{s}^2.$$
$`(7)`$
## V The Cosmological Relativity Theory

In Ref. 19 a four-dimensional cosmological relativity theory that unifies space and velocity was presented. The theory predicts that the universe accelerates and hence it is equivalent to having a positive value for $`\mathrm{\Lambda }`$ in it. As is well known, in the traditional work of Friedmann when added to it a cosmological constant, the field equations obtained are highly complicated and no solutions were obtained so far.

Cosmological relativity theory, on the other hand, yields exact solutions and describes the universe as having a three-phase evolution with a decelerating expansion followed by a constant and an accelerating expansion, and it predicts that the universe is now in the latter phase. In the framework of this theory the zero-zero component of Einstein’s equations is written as

$$R_0^0\frac{1}{2}\delta _0^0R=\kappa \rho _{eff}=\kappa \left(\rho \rho _c\right),$$
$`(8)`$
where $`\rho _c=3/\kappa \tau ^23H_0^2/\kappa `$ is the critical mass density and $`\tau `$ is Hubble’s time in the zero-gravity limit.

Comparing Eq. (8) with the zero-zero componentof Eq. (1), one obtains the expression for the cosmological constant in cosmological relativity theory,

$$\mathrm{\Lambda }=\kappa \rho _c=3/\tau ^23H_0^2.$$
$`(9)`$
Assuming that Hubble’s constant $`H_0=85`$km/s-Mpc, then

$$\mathrm{\Lambda }=2.28\times 10^{35}\text{s}^2.$$
$`(10)`$
This result is in excellent agreement with the recent supernovae experimental results.

## VI Conclusions

We have seen how the cosmological constant can be determined in a natural way without even adding it explicitely to the Einstein field equations. Rather, the introduction of the effective mass density $`\rho _{eff}=\rho \rho _c`$ is enough to ensure that the universe expands in the same way using Einstein’s field equations with a cosmological constant. But there is a big difference now: The theory determines the numerical value of the cosmological constant, and experiments confirm it.


text:
# 1 Introduction

## 1 Introduction

FUSE observations towards HD 34078 have been undertaken to study the small scale structure in the spatial distribution of H2 (see the paper by Boissé et al. in this volume). The first FUSE spectrum turned out to be exceptionally rich since absorptions from the 19 first levels of H2 have been identified. After a brief presentation of the line of sight, we describe the results obtained for H2 and HD. A preliminary interpretation of the observed excitation is then performed using Photon Dominated Region (PDR) and magnetohydrodynamic shock models.

## 2 Description of the line of sight to HD 34078

HD 34078, also known as AE Aurigae, is a O9.5V spectral type star of visual magnitude mV= 5.95, located at 500 pc from the sun. It is a runaway star with a transverse velocity of 100 km$``$s-1 and a radial one of 59 km$``$s-1. The color excess has been estimated to be EB-V= 0.52, corresponding to AV=1.6 for RV=3.1. Several observations of HD 34078 have already been performed: HI and CO column densities were measured by Mc Lachlan and Nandy  from IUE observations while CH, CN and C2 have been detected by Federman et al. . The density is estimated to be in the range 200 - 500 cm-3. In addition, CH+ and CH, observed by M. Allen , give a Doppler parameter, b $``$ 3 km$``$s-1.

## 3 Observations

The first observation has been performed in january 2000, with a total integration time of 7145s. Data have been processed and calibrated in a standard way using the FUSE pipeline software. In the best portions of the LiF1A spectrum, the S/N ratio is 30 per 15 mÅ pixel; the spectral resolution is about 17000. Fig. 1 displays an extract of the LiF1A spectrum over a 10 Å interval around 1050 Å , showing the numerous absorption lines present in the spectrum.

Column densities have been obtained using both a curve of growth analysis and the fitting program ”Owens” written by M. Lemoine. Since high resolution CH observations performed in the visible have revealed no marked velocity structure (Federman, private communication), we assumed that there is only one component along the line of sight. Most of the observed lines fall on the flat part of the curve of growth and then poorly constrain the b Doppler parameter. We have thus adopted the value derived by Allen  for the lowest levels of H2 (up to J=4): b=3 km$``$s-1. However, using ten transitions arising from the J=9 state, we have been able to build a curve of growth for this level, from which we get an upper limit on b(J=9) of 6 km$``$s-1. For H2 transitions from J=4 to J=9, we have increased the b parameter from 3 to 5 km$``$s-1; an uncertainty of 1 km$``$ s-1 on the b parameter has been considered. Error bars are given at three sigma, meaning that values outside of the error bars can be excluded. Transitions arising from levels higher that J=9 are optically thin and thus, corresponding column densities are independent on any assumption on b; the same is true for the J=0 and 1 levels for which absorption lines are damped.

## 4 Results

We clearly detect absorption lines arising from vibrationally excited H2 levels up to v=1 and J=5 (8000 K). Other molecules and atoms such as HD, CO, C I, Ar etc… are detected as well. Absorption lines from levels with still higher excitation energies (e.g. v=4, J=2) may also be present but, due to the numerous blends present, it is difficult to establish unambiguously their identification. The table displays the column densities derived for the 19 first energy levels sorted by increasing energy. We derive a total molecular hydrogen column density of:

$$N(H_2)=6.4\times 10^{20}cm^2$$

Seven lines from HD J=0 are clearly detected but some of them are blended and it appears difficult to fit all lines with a single N(HD) value. We get the approximate value, N(HD J=0) $``$ 1015 cm-2; the accuracy is too low to derive any D/H ratio. The weak HD L7-0R(1) transition arising from the J=1 level is tentatively detected at 1022 Å giving N(HD J=1)) = 1.0 (-0.4 +0.8) $`\times `$ 1014 cm-2. A previous detection of excited HD was obtained using Copernicus by Wright and Morton  towards $`\zeta `$ Oph.

## 5 Discussion

Figure 2 displays the excitation diagram of H2. The 3 first H2 levels are consistent with a temperature of 80 K, and a thermal equilibrium ortho-to-para ratio. For higher excited levels, the distribution can be artificially decomposed in two thermal ones with excitation temperatures of 400K (for levels between J=3 and J=7) and 1000K (beyond J=7). Using first the PDR model of Le Bourlot et al.  and adopting a density of 500 cm-3 (constrained by the observations of Federman et al ) together with a radiation field of 2 times the ISRF, it is possible to reproduce the column density of the two first H2 levels and other molecules, except CH+. However, this model considerably underestimates column densities for levels J $``$ 3. A more energetic mechanism is required to account for the latter. We have run several C shock models in the presence of a radiation field (reference ). Reasonable agreement is obtained for levels J=3 to J=7 with a shock velocity of 25 km$``$s-1 and a magnetic field of 10 $`\mu `$G. But again, the large amount of H2 observed at J $``$ 8 is well above the prediction of these models. One plausible scenario could involve H2 located in a bow shock located around HD 34078. Since it is a runaway star, it is now well detached from the cloud where it formed, unlike HD 37903 towards which Meyer et al.  found a lot of H2 excited to high energies by fluorescence. However, due to the large velocity of HD 34078, a strong bow shock is expected to be present at the interface of the stellar wind and of the ambient ISM (Mc Low et al. ). Gas trapped in this dense layer might contribute to the observed absorption \- a scenario that we need to explore in more details.


text:
# VLHC predictions for 𝐻→𝜏𝜏 in weak boson fusion

## Abstract

Higgs production in weak boson fusion with subsequent decay $`H\tau \tau e\mu /p_T`$ provides a means to measure Higgs Yukawa couplings and Higgs interactions to weak bosons. The potential precision of cross section measurements at a VLHC is investigated.

preprint: CALT-68-2356preprint: MADPH-01-1246

Precision measurements of Higgs boson properties are an important goal of any post-LHC accelerator. The LHC is expected to yield a determination of the Higgs boson mass at the $`10^3`$ level ATLAS ; CMS . A combination of observations of the Higgs boson, in various decay channels ($`H\gamma \gamma ,WW,ZZ,\tau \tau `$), in inclusive searches and in the weak boson fusion (WBF) channel, is expected to yield measurements of various Higgs partial widths with accuracies of order 10–20% zknr . A linear collider (LC) can improve these measurements by up to one order of magnitude, achieving a clean separation of individual couplings lc .

In contrast, very little is known about the capabilities of a higher energy hadron collider for improving our knowledge of the Higgs sector. In this note we analyze a particular weak boson fusion process, $`qqqqH,H\tau \tau `$ with subsequent leptonic tau decays, within the SM, as a case study for possible VLHC reach in precision measurements of Higgs couplings. Unlike gluon fusion, WBF has small NLO QCD corrections, thus promising small systematic uncertainties. Second, the $`H\tau \tau `$ decay modes prz\_tautau  are key processes in the coupling determination at the LHC. We closely follow this previous LHC analysis and determine signal rates and cross sections for the main physics backgrounds in $`pp`$ collisions at $`E_{\mathrm{cm}}=50`$, 100, and 200 TeV. These cross sections allow for a first estimate of the statistical accuracy of coupling measurements at a VLHC.

The signal process is $`qqqqH`$, which is mediated by $`t`$-channel $`W`$ and $`Z`$ exchange. We consider $`H\tau \tau `$ decays with subsequent leptonic decays of both taus, i.e. the final state is $`jj\mathrm{}^+\mathrm{}_{}^{}{}_{}{}^{}/p_T`$ where, for simplicity, we only consider the $`e^\pm \mu ^{}`$ combinations.In a previous LHC analysis prz\_tautau  it was found that the physics backgrounds, $`Zjj`$ production with $`Z\tau ^+\tau ^{}`$, dominate overall backgrounds. Reducible backgrounds are mostly from leptonic decays of $`W`$ pairs, e.g. from top-decays or $`WWjj`$ events. These can be distinguished from tau pairs by the angular correlation of charged lepton momenta and the missing momentum vector. Thus, we only consider QCD $`Zjj`$ and electroweak $`Zjj`$ production, i.e. $`\tau ^+\tau ^{}jj`$ production at order $`\alpha _s^2\alpha ^2`$ and $`\alpha ^4`$, respectively. These backgrounds and the signal are generated as in Ref. prz\_tautau  and include off-shell $`Z`$’s and photon-$`Z`$ interference.

As compared to the LHC analysis we impose a higher $`p_T`$ of the two forward “tagging jets” and a larger dijet invariant mass of these two candidate quark jets. Specifically the signal is required to have two jets and two charged leptons with

$`p_{T_j}30\mathrm{GeV},|\eta _j|5.0,\mathrm{}R_{jj}0.6,`$
$`p_{T\mathrm{}_1}20\mathrm{GeV},p_{T\mathrm{}_2}10\mathrm{GeV},|\eta _{\mathrm{}}|3,\mathrm{}R_j\mathrm{}0.6,`$ (1)
the two jets must be separated by at least 4.2 units of pseudorapidity, reside in opposite detector hemispheres and the two charged leptons must lie between the two jet definition cones of radius 0.6. The neutrinos must lead to missing transverse momentum of 50 GeV or more and we require

$$M_{jj}>1000\mathrm{GeV}$$
(2)
for the invariant mass of the two tagging jets. All other cuts are as in Ref. prz\_tautau .

The surviving $`e\mu /p_T+2`$ jet events allow a reconstruction of the original tau-pair invariant mass by making use of the large boost of the individual taus and their decay products tautaumass : the charged leptons trace the original tau directions and the transverse momenta of $`e`$, $`\mu `$ and missing neutrinos allow to solve for the $`\tau ^+`$ and $`\tau ^{}`$ energies. The reconstruction relies on good missing transverse momentum resolution. We assume a performance of VLHC detectors equivalent to ATLAS, i.e. energy smearing and $`/p_T`$ resolution is handled as in Ref. prz\_tautau  and reproduces ATLAS simulations for $`H/A\tau \tau `$ events cavalli .

The resulting reconstructed $`M_{\tau \tau }`$-spectra are shown in Figure 1 for VLHC energies of 50 and 100 TeV. The Higgs peak is clearly visible above the physics backgrounds, even for Higgs masses as large as 155 GeV. The statistical significance of the signal can be estimated by considering the signal and background rates in a reasonable mass window, which we take as $`m_H15\mathrm{GeV}<M_{\tau \tau }<m_H+15\mathrm{GeV}`$.

Resulting cross sections are given in Table 1. It is quite remarkable that we may expect signal to background ratios well above 1 over the entire intermediate mass range. This is combined with large signal rates, of order hundreds to a thousand, for an integrated luminosity of 100 fb-1. Given the expected number of signal events, $`S`$, and background events, $`B`$, we can estimate the statistical error with which the signal cross section can be estimated as $`\mathrm{\Delta }\sigma _H/\sigma _H=\sqrt{S+B}/S`$. The corresponding statistical errors are listed in the last column of Table 1 and reside in the 5–10% range for a 50 TeV machine and Higgs masses between 120 and 150 GeV, and improve to the 3–6% range for a 200 TeV machine.

The numbers presented in the left half of Table 1 do not include efficiency factors for jets or for lepton identification and isolation. In previous studies these factors amounted to overall efficiency factors of 0.67 for signal and backgrounds prz\_tautau . (Note that geometric acceptance is included in our simulation.) Other reduction factors for the signal will arise from parton shower corrections: including QCD radiation, the selection of the 2 quark jets is not unique. Misidentified tagging jets then tend to fall below the $`M_{jj}`$ threshold of Eq. (2). Since QCD radiation is relatively suppressed for the $`t`$-channel color singlet exchange of the signal, this effect is modest wbfexp  and we may assume an overall efficiency of $`ϵ_{\mathrm{WBF}}=0.5`$ for WBF processes compared to $`ϵ=0.67`$ for the QCD $`Zjj`$ background.

The cross sections of Table 1 do not yet make use of another characteristic of WBF events: they are mediated by the $`t`$-channel exchange of a colorless object. Destructive interference between initial and final state radiation in such processes leads to suppressed gluon radiation in the central region, between the two quark jets bjgap . QCD background processes, which are dominated by $`t`$-channel gluon exchange, predominantly radiate in the central region. This difference can be exploited by a veto on central jets, in the region

$$\eta _{j,min}^{\mathrm{tag}}+0.6<\eta _j^{\mathrm{veto}}<\eta _{j,max}^{\mathrm{tag}}0.6.$$
(3)
A detailed study on feasible $`p_T`$ thresholds at a VLHC is still needed. For the LHC, an parton level analysis indicates that a veto of central jets of $`p_T>20`$ GeV will suppress the signal by a factor $`P_{\mathrm{surv}}=0.9`$, while analogous factors for the EW and the QCD $`Zjj`$ backgrounds are 0.75 and 0.3, respectively DR\_thesis . Due to the increased $`M_{jj}`$ values and the larger rapidity distance between the quark jets of the signal, these numbers may improve at a VLHC. However, higher $`p_T`$ thresholds of vetoed jets will counteract such improvements. We take the LHC values as an educated guess for a VLHC as well. Including these efficiency factors and survival probabilities, signal and background rates are given in the right half of Table 1. Signal to background rates improve further, due to the background suppression of the central jet veto. However, the statistical accuracy of the cross section determination suffers somewhat, due to the more realistic assumptions on reconstruction efficiencies. Nevertheless, statistical errors of 5% or better are possible at a VLHC with a mere 100 fb-1 of data. This can be improved significantly by including other tau decay channels like hadronic decays of one tau and leptonic decays into $`e^+e^{}/p_T,\mu ^+\mu ^{}/p_T`$ prz\_tautau ; zknr .

A point of concern are Higgs backgrounds. $`HWWe\mu /p_T`$ events become the dominant reducible $`WW`$ background for Higgs masses above 125 GeV at the LHC, and contribute about a third of the overall background at $`m_H=150`$ GeV prz\_tautau . The same should be expected at a VLHC. In addition, gluon fusion becomes a more pronounced source of $`Hjj`$ events at the higher VLHC energies oleari , and thus the clean separation of gluon fusion and WBF cross sections will require additional efforts. However, we do not expect these complications to significantly alter the previous optimistic findings.

We conclude that a VLHC can study $`H\tau \tau `$ decays with high precision in WBF events. Integrated luminosities of a few hundred fb-1 provide high statistics and clean samples which allow to measure $`H\tau \tau `$ Yukawa couplings with statistical errors at the $`10^2`$ level. At the same time the small QCD NLO corrections to WBF promise small systematic errors. Detailed studies of WBF channels at a VLHC are required to assess experimental systematic errors.

###### Acknowledgements.

This research was supported in part by the University of Wisconsin Research Committee with funds granted by the Wisconsin Alumni Research Foundation and in part by the U. S. Department of Energy under grant No. DE-FG02-95ER40896. S.S. was supported by DOE grant DE-FG03-92-ER-40701.


text:
# CERN-TH/2001-269 HU-EP-01/37 October 2001 The charm quark’s mass in quenched QCD presented at the XIX International Symposium on Lattice Field Theory “Lattice 2001”, August 19 – 24, 2001, Berlin, Germany

## 1 INTRODUCTION

Quark masses are among the fundamental parameters of the Standard Model, and yet not directly accessible to experimental measurement. Their determination in terms of experimental quantities requires a good quantitative control of the non-perturbative dynamics of QCD, which may be achieved by numerical simulations of the lattice regularized theory. Extending previous work on the $`\mathrm{\Lambda }`$-parameter and the strange quark mass in quenched QCD  we report here on a precise determination of the charm quark mass, taking the experimental value of the $`D_s`$ meson mass as essential new input. For previous work on the charm quark mass we refer to two recent papers  and the references given there.

## 2 STRATEGY

We use O($`a`$) improved Wilson quarks with the Wilson gauge action. All necessary renormalization constants and improvement coefficients are known with non-perturbative accuracy in the relevant range of bare couplings. We refer to  for a review and further references.

### 2.1 Setting the parameters

In quenched QCD it is convenient to measure dimensionful quantities in units of the scale $`r_0`$ , since $`r_0/a`$ has been determined very precisely for a large range of bare couplings . To obtain the scale in physical units we always use $`r_0=0.5\mathrm{fm}`$, which is roughly equivalent to setting the scale with $`F_K=160\mathrm{MeV}`$ . The renormalization group invariant (RGI) strange quark mass has been obtained in ,

$$r_0M_s=0.348(13)M_s=138(6)\mathrm{MeV},$$
(1)
taking $`m_K`$ from experiment and the mass ratio

$$M_s/\widehat{M}=24.5\pm 1.5,\widehat{M}=\frac{1}{2}(M_u+M_d),$$
(2)
from chiral perturbation theory . Up to intrinsic O($`a^2`$) ambiguities, the connection between $`M_s`$ and the bare strange quark mass can be established by combining the results of refs. , so that no tuning is required to satisfy the renormalization condition (1). The bare charm quark mass may then be fixed by matching directly the experimental result $`m_{D_s}=1969\mathrm{MeV}`$ , which translates to $`r_0m_{D_s}=4.98`$. This is justified since electroweak effects on $`D`$ meson masses are expected to be small on the scale of the expected statistical errors.

### 2.2 Some technical details

The simulations were done at four values of $`\beta `$, $`6.0\beta 6.45`$, which correspond to lattice spacings in the range $`0.050.1\mathrm{fm}`$. Spatial lattice volumes ranged from $`16^3`$ to $`32^3`$, and the Euclidean time extent from $`40`$ to $`80`$. The linear extent of the spatial volume was never less than $`1.5\mathrm{fm}`$ in physical units, so that finite volume effects on the $`D_s`$ meson mass can be safely neglected. Quark propagators were computed for two mass parameters in the strange quark region, and three values around the expected charm quark mass. To compute meson masses we used correlation functions derived from the Schrödinger functional (SF) . In the pseudoscalar channel these are given by

$$f_\mathrm{X}(x_0)=\frac{1}{2}a^6\underset{𝐲,𝐳}{}\overline{c}(x)\mathrm{\Gamma }_\mathrm{X}s(x)\overline{\zeta }_s(𝐲)\gamma _5\zeta _c(𝐳),$$
(3)
where $`\mathrm{\Gamma }_\mathrm{X}`$ stands for $`\gamma _0\gamma _5,\gamma _5`$ for $`\mathrm{X}=\mathrm{A},\mathrm{P}`$ respectively. The pseudoscalar meson mass is then obtained by looking for a plateau in the effective masses,

$$am_{\mathrm{eff}}^\mathrm{X}(x_0+\frac{1}{2}a)=\mathrm{ln}\left\{f_\mathrm{X}(x_0)/f_\mathrm{X}(x_0+a)\right\},$$
(4)
as a function of Euclidean time $`x_0`$. Autocorrelation times were found to be completely negligible, so that the produced O(100) configurations at each $`\beta `$-value are statistically independent. An example for an effective mass plot is given in fig. 1. Excited states are expected to contribute significantly when either $`x_0`$ or $`Tx_0`$ becomes small. It is possible to estimate these contributions and thereby quantify the resulting deviation from the true plateau value. Requiring this effect to be smaller than $`0.5\%`$, the plateau region in fig. 1 was identified with the interval $`25x_0/a35`$; the meson mass was then estimated by averaging the effective mass in this interval.

### 2.3 Bare and renormalized quark masses

Given the pseudoscalar masses at the chosen parameters, an interpolation in some bare charm quark mass is required to match the experimental $`D_s`$ meson mass. Besides the bare mass $`m_0`$ we also consider two alternative definitions of the bare charm quark mass, which derive from the PCAC relation. First we use the heavy–light axial current and density to define,

$$m_{cs}(x_0)=\frac{\stackrel{~}{}_0f_\mathrm{A}(x_0)+c_\mathrm{A}a_0^{}_0f_\mathrm{P}(x_0)}{2f_\mathrm{P}(x_0)},$$
(5)
which upon renormalization yields the average charm and strange quark mass. With our parameter choices yet another bare charm quark mass can be obtained from the PCAC relation involving a hypothetical strange quark, assumed to be mass degenerate with the charm quark. We denote this PCAC mass as $`m_{cc}`$ but one should bear in mind that its definition does not involve flavour singlet operators.

The renormalized and O($`a`$) improved charm quark mass $`m_\mathrm{R}`$ can now be obtained in various ways. Choosing the SF scheme we first obtain

$$m_\mathrm{R}=Z_\mathrm{P}^1Zm_\mathrm{q}(1+b_\mathrm{m}am_\mathrm{q}),$$
(6)
where $`m_\mathrm{q}`$ is the subtracted bare charm quark mass. From the PCAC mass $`m_{cs}`$, one gets

$$m_\mathrm{R}+m_{\mathrm{R},s}=Z_\mathrm{P}^1Z_\mathrm{A}\left\{1+(b_\mathrm{A}b_\mathrm{P})a\overline{m}_\mathrm{q}\right\}2m_{cs},$$
(7)
where $`\overline{m}_\mathrm{q}`$ denotes the average of the bare subtracted charm and strange quark masses, and we have assumed $`x_0=T/2`$, which is kept approximately constant in physical units. Finally, we consider

$$m_\mathrm{R}=Z_\mathrm{P}^1Z_\mathrm{A}\left\{1+(b_\mathrm{A}b_\mathrm{P})am_\mathrm{q}\right\}m_{cc},$$
(8)
and convert to the RGI quark mass using the flavour-independent ratio

$$M/m_\mathrm{R}=1.157(15).$$
(9)
## 3 RESULTS AND CONCLUSIONS

After subtraction of the strange quark mass (1) in eq. (7) we have three definitions of $`r_0M_c`$ which should coincide up to terms of O($`a^2`$). In each case we extrapolate to the continuum with an ansatz of the form

$$r_0M_c=c_0+c_1(a/r_0)^2,$$
(10)
excluding the data point at the coarsest lattice spacing ($`\beta =6.0`$). As can be seen in fig. 2, the results are nicely compatible with each other in the continuum limit, although the differences at finite cutoff are quite large. As our best (preliminary) estimate we quote

$$r_0M_c=4.26(16)M_c=1684(64)\mathrm{MeV}.$$
(11)
Using the four-loop anomalous dimension and $`\beta `$-function in the $`\overline{\mathrm{MS}}`$ scheme , and $`r_0\mathrm{\Lambda }_{\overline{\mathrm{MS}}}=0.586(48)`$ , we obtain

$$m_c(m_c)=1314(40)(20)(7)\mathrm{MeV},$$
(12)
where the first error stems from the RGI quark mass, the second is induced by the error on $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$, while the last number is the difference between three- and four-loop perturbative evolution down to the charm scale.

We conclude that precise results can be obtained for charmed observables using the same techniques as for the light quarks, provided that the continuum extrapolation is performed. We emphasize that we have used the quenched approximation, which is known to lead to inconsistencies in the hadron spectrum at the $`10\%`$ level. However, the present result constitutes a significant improvement over previous lattice estimates and will be a useful reference point for more realistic calculations involving dynamical quarks. For the time being we speculate that the mass ratio $`M_c/M_s=4.26(16)/0.348(13)12.2\pm 1.0`$ may be less sensitive to quenched scale ambiguities than the masses themselves.

This work is part of the ALPHA collaboration research programme, and partially supported by the European Community under the grant HPRN-CT-2000-00145 Hadrons/Lattice QCD. Simulations were carried out on machines of the APE100 and APE1000 series at DESY-Zeuthen. We thank the staff at the computer centre for their help, and P. Ball, R. Sommer, H. Wittig and U. Wolff for useful discussions.


text:
# Reconstruction of dynamical equations for traffic flow

## I Introduction

By increasing the number of licensed vehicles on our roads it becomes more and more necessary to reduce appearances of traffic congestion. Getting higher capacities of highways means looking for optimized flow rates of cars. In order to do so, it is necessary to investigate the complex system of traffic flow and to understand its regularities. In a subsequent step, intelligent traffic control systems may use these laws to influence the traffic flow and thus to increase the highway’s capacity utilization.

For our investigation a large amount of traffic flow data were collected at the highway near Köln-Nord (Germany) over more than one week. For each car crossing an induction loop detector the following data were recorded: (i) the time, when the car crossed the detector, with an resolution of 1 sec, (ii) the type of car (passenger car, truck), (iii) the lane number, (iv) the velocity, truncated to an 8-bit integer (0…255 km/h), accuracy of approximately 3 %, (v) the length of the car in meters (8-bit integer, not calibrated) and (vi) the distance to the car driving ahead (an integer in the range 0…999 m, not calibrated). The lanes are labeled from A to C, where lane A is the right driving lane where mostly trucks are found (because of a German law), and lanes B and C are the fast drivers’ lanes.

In the following we present a new method to derive from the measured data dynamical equations for the traffic situation. We start with the common presentation of the fundamental diagram. Next, we evaluate the deterministic and the stochastic content of of the traffic dynamics by means of a Langevin equation. Finally, we give an interpretation of the determined Langevin equation.

## II Fundamental diagrams

Here and in the following we calculate the car density $`k`$ \[km-1\] and the flux state $`q`$ \[h-1\] of each car using the measured data: the velocity $`v`$ \[km/h\], the length of the car $`l`$ \[m\] and the distance $`d`$ \[m\] to the car ahead:

$`k`$ $`=`$ $`{\displaystyle \frac{1000}{l+d}},`$ (1)
$`q`$ $`=`$ $`kv.`$ (2)
To avoid an overload of the presentation we restrict the diagrams and the following calculations to two cases. Firstly, only the traffic of a single lane C, secondly the cumulative traffic of all three lanes A, B and C are considered. Furthermore we calculate for each density state $`k`$ the mean flux state

$$q_k=\frac{1}{n(k)}q(k).$$
(3)
The obtained fundamental diagrams (qi;ki) and ($`<`$q$`>`$;k) are shown in Fig. 1. We note that for this presentation no significant difference in the traffic dynamics of one lane and the cumulative dynamics of three lanes can be detected. In both diagrams we find a maximum flux in free traffic flow of $`q_{\mathrm{max}}=8000`$ cars/h, according to the results of . For the flux out of traffic jams we find in both cases $`q_{\mathrm{out}}=5000`$ cars/h. So we have a ratio of $`\eta :=q_{\mathrm{max}}^{\mathrm{free}}/q_{\mathrm{out}}=1.6`$, which meets the value $`\eta 1.5`$ found in  quite well.

## III Langevin equation for the traffic flow

In order to grasp the underlying dynamics of the traffic flow (one lane and all three lanes), we utilize a new method to analyze the traffic data more extensively. In particular, the iterative dynamics of the traffic state $`𝐱_N`$, given by the velocity $`𝐯_N`$ and the flux $`𝐪_N`$ of the $`N`$-th car, as a function of traffic state $`𝐱_{N1}`$ of the $`N1`$-th car is investigated. Note, that other state variables could have been choosen as well.

For the iterative dynamics of the traffic state $`𝐱_N`$ (which may be taken in a generalized case as a $`r`$-dimensional variable) we propose a description by a stationary Langevin equation, taking into account a combination of deterministic and random (noisy) forces cf. :

$$x_{i,N+1}=x_{i,N}+h_i(𝐱_N)+\underset{j=1}{\overset{r}{}}g_{ij}(𝐱_N)\mathrm{\Gamma }_{j,N},$$
(4)
where the indices $`i`$ and $`j`$ denote components of the multidimensional variables, $`N`$ the time, and $`\mathrm{\Gamma }_{j,N}`$ are $`r\times N`$ independent Gaussian noise variables with zero mean and with variance 2, i.e.,

$$\mathrm{\Gamma }_{j,N}=0,\mathrm{\Gamma }_{j,M}\mathrm{\Gamma }_{k,N}=2\delta _{j,k}\delta _{M,N}.$$
(5)
The central part of our following work is that it is possible to determine the functions $`h_i`$ and $`g_{ij}`$ directly from empirical data. Taking (4) as the Ito presentation of the stochastic process the following relation to the Kramers-Moyal coefficients can be given,

$`D_i^{(1)}(𝐱)`$ $`=`$ $`h_i(𝐱)`$ (6)
$`D_{ij}^{(2)}(𝐱)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{r}{}}}g_{ik}(𝐱)g_{jk}(𝐱),`$ (7)
where $`D^{(1)}`$ and $`D^{(2)}`$ are called drift and diffusion coefficient. These coefficients can be evaluated by the conditional moments

$`D_i^{(1)}(𝐱)=x_{i,N+1}x_i|_{𝐱_N=𝐱},`$ (8)
$`D_{ij}^{(2)}(𝐱)=(x_{i,N+1}x_i)(x_{j,N+1}x_j)|_{𝐱_N=𝐱}.`$ (9)
Recently it has been shown that with the analogous definition of these Kramers-Moyal coefficients it is possible to reconstruct from time continuous dynamics the underling stochastic differential equation .

Before presenting our results on the dynamics we want to comment on the validity of this ansatz to describe the traffic flow by the Langevin equation (4). This ansatz implies that the dynamics is in the class of Markovian processes, i.e. the system does not have a memory. This can be tested by conditional probabilities

$$p(𝐱_𝐍|𝐱_{𝐍\mathrm{𝟏}},\mathrm{},𝐱_{𝐍𝐦})=𝐩(𝐱_𝐍|𝐱_{𝐍\mathrm{𝟏}})$$
(10)
or by the necessary condition of the Chapman-Kolmogorov equation

$$p(𝐱_𝐍|𝐱_{𝐍𝐫})=\underset{𝐱_{𝐍𝐬}}{}𝐩(𝐱_𝐍|𝐱_{𝐍𝐬})𝐩(𝐱_{𝐍𝐬}|𝐱_{𝐍𝐫}),$$
(11)
where $`r>s`$. From our data, conditional probabilities have been evaluated and the validity of the Chapman-Kolmogorov equation was found for the iterative dynamics of both quantities, the velocity and the flux. If this Markovian property holds, the inherent noise of the dynamics ($`\mathrm{\Gamma }`$) can be taken as $`\delta `$-correlated. It should be noted, that even in the case where the noise is not $`\delta `$-correlated, the deterministic part of the dynamics can be reconstructed from given data (8), as we found by analysing numerically generated test data, .

As expressed by (9), the knowledge of the conditional probabilities $`p(𝐱_𝐍|𝐱_{𝐍\mathrm{𝟏}})`$ provides the basis to estimate the Kramers-Moyal coefficients from the traffic data. First we consider the simplified case of the onedimensional dynamics of the velocity only. The results for and $`D^{(2)}`$(v) for the traffic of one lane and for the cumulative traffic of all three lanes are shown in Fig. 2. The one dimensional deterministic dynamics can also be expressed by the potential $`\mathrm{\Phi }_D`$(v), defined as - $`\frac{\delta \mathrm{\Phi }_D}{\delta v}`$ = D(1)(v). The corresponding potentials are shown in Fig. 3. From these results three noticable velocities $`v_137`$ km/h, $`v_2`$ 75 km/h and $`v_3`$ 107 km/h appear, which allows to identify the so called congested flow for $`vv_1`$, correlated flow for $`v_1<vv_2`$ and the free flow for $`v>v_2`$ , respectively. Note these velocities can be defined as fixed points ($`D^{(1)}(v)=0`$) of the deterministic part of the cumulative traffic dynamics (see Fig. 2b).

For the traffic dynamics of lane C we find in the congested and in the correlated regime metastable traffic states, the deterministic drift term $`D^{(1)}`$ gets zero over finite intervals. This corresponds to the plateau structure in the potential, see Fig. 3a. A clearly different behaviour is found for the cumulative traffic dynamics of all three lanes, see Figs. 2b and 3b. The different flow regimes are seperated by two fixed points at $`v_1`$ and $`v_2`$. The slope of these fixed points defines the stability, thus the congested and the correlated flow regimes are separated by a stable fixed point, whereas the correlated and the free flow regime are separated by an instable fixed point. For the free flow, in both cases of one lane or three lane traffic a stable fixed point is found at $`v_3`$, correspondingly the drift potential has its local minimum. Because of a speed limit of 100 km/h given on the inspected highway we see a great increase of the potential for $`v>v_3`$: the faster a car is driving, the stronger the attraction is to the potential’s minimum  .

To get an understanding of the real traffic dynamics grasped by these drift coefficients or drift potentials, the additional noise has to be taken into account. In Fig. 2c and d the corresponding magnitude of the noise are expressed by the evalutaed diffusion coefficients $`D^{(2)}`$. The noise will now cause transitions between different flow states. For the traffics dynamics of one lane, the noise will effect larger fluctuations as it is the case for the cumulative traffic dynamics, which has two clear minima in the potential. A further interesting detail is that the magnitude of $`D^{(2)}`$ has a minimum around the stable fixed point at $`v_3`$. This indicates a pronounced stability of this traffic situation.

Next we present the results of a higher dimensional analysis by taking $`𝐱`$ with the components $`x_1=v`$ and $`x_2=q`$. Now also the drift coefficient $`D^{(1)}`$ becomes a vector depending on $`v`$ and $`q`$, as shown in Fig. 4. These results were obtained by binning the velocity flux data into a $`25\times 25`$ matrix, corresponding to a binning of the velocity into intervals of $`\mathrm{\Delta }v=5`$ km/h.

The solid lines in Fig. 4 show the states where we have no drift of the flux component: $`D_q^{(1)}=0`$. On these lines we find only a velocity drift with a constant flux. In accordance to Fig. 2a we find in Fig. 4a for slow velocities mainly no velocity drift, in Fig. 4b (like in Fig. 2b) there seem to exist stable velocity drift states at the same velocities ($`v_i37`$km/h,$`107`$ km/h). The topology of the instable fixed point gets now a saddle point which is attractive for larger and smaller flux values but instable in the direction of larger and smaller velocity values. Again clear differences of the dynamics of on lane and three lanes is found.

## IV Discussion and Conclusion

By the investigation of traffic flow data as an iterative stochastic process we were able to calculate from the given data the 1-dimensional drift and diffusion coefficients and thus to find the deterministic and stochastic part of the corresponding Langevin equation. We are able to find stable, metastable and unstable states (fixed points) in the deterministic part of free, correlated and congested traffic flow. For a fully description of the whole dynamics also the diffusion coefficient has to be taken into account, which provides transition probabilities between the different (meta-)stable states. Without this noisy part of the traffic flow dynamics, the stable states would never be left, i.e. a congestion would stay forever if once prepared.

To see the dependency of velocity and flow from each other, the investigations were expanded to a higher dimensional analysis. In this case we find the deterministic and stochastic part of the 2-dimensional Langevin equation. Now we are able to identify stable velocity and flux states of the deterministic part.

Interestingly, the results found in this study are in agreement with empirical investigations that have identified three phases of traffic , together with transitions connecting theses phases. Especially the transition from free flow to correlated flow (or synchronized flow in the terminology of ) has similarities.

Finally we want to point out, that we presented here a new method to analyse traffic data with respect to a derivation of dynamical equations from pure data analysis. Furthermore we could show that our method provides more insight into the traffic dynamics than the presentations of diagrams like the fundamental diagram. A clear difference in the dynamics of one lane and the dynamics of cumulative three lanes was found. Our analysis provides evidence of the presence of fixed points, which are of practical importance if a control of a traffic should be achieved. At last one may conclude that this method will be helpful to perform a more thorough comparison between traffic flow models and empirical data.

Acknowledgement: Helpful discussions with Ch. Renner, St. Lück and M. Siefert are acknowledged. We also would like to thank the Landschaftsverband Rheinland and the Northrhine-Westfalia Ministry for Economy and Transport for providing the data used in this study.