text:
# Solving the boundary value problem for finite Kirchhoff rods

## I Introduction

The study of conformations of slender elastic rods is of substantial utility in several applications, ranging from the fields of structural mechanics and engineering to biochemistry and biology. Examples are the study of coiling and loop formation of sub-oceanic cables coyne ; zajac ; sun ; vaz , filamentary structures of biomolecules zimm ; yang ; shi ; tamar ; golds  and bacterial fibers golds2 ; klapper , the phenomenon of helix hand reversal in climbing plants goriely  and the shape and dynamics of cracking whips alain3 .

The Kirchhoff model kirchhoff  provides a powerful approach to study the statics and dynamics of elastic thin rods tamar ; olson . In this model the rod is described by a set of nine partial differential equations (the *Kirchhoff equations*) in the time and arclength of the rod. They contain the forces and torques plus a triad of vectors describing the deformations of the rod. These equations are the result of Newton’s second law for the linear and angular momentum applied to the thin rod plus a linear constitutive relationship between moments and strains. The Kirchhoff model holds true in the approximation of small curvatures of the rod, as compared to the radius of the local cross section dill . An interesting characteristic of this model, known as *Kirchhoff kinetic analogy*, is that the equations governing the static problem are formally equivalent to the Euler equations describing the motion of spinning tops in a gravity field. The Kirchhoff equations for equilibrium configurations can, therefore, be written in hamiltonian form.

The Kirchhoff equations have been solved for a number of simple situations. Shi and Hearstshi  first obtained analytical solutions of the equilibrium equations and, recently, Nizzete and Goriely nizzete  completed the study by making a classification of all kinds of equilibrium solutions. Goriely and Tabor tabor1 ; tabor2  developed a method to study the dynamical stability of these solutions and Fonseca and de Aguiar fonseca  applied this method to study the near equilibrium dynamics of non-homogeneous closed rods in viscous media. Recently, Tobias *et al.* tobias  developed the necessary and sufficient criteria for elastic stability of equilibrium configurations of closed rods.

In many cases of interest, including biological molecules, the filaments are subject to boundary conditions. Examples are the problem of multiprotein structures, such as histones and gyrase, about which long pieces of DNA wrap tobias2 , and multiprotein structures, such as the *lac repressor* complex mahade . Despite the many achievements described above, the study of the boundary value problem (BVP) associated with Kirchhoff filaments is still a big challenge. While the integration of differential equations from initial conditions is a relatively simple numerical task, the difficulties of finding solutions for given boundary conditions are well known in classical mechanics, electromagnetism and quantum mechanics. Typical examples are classical trajectories connecting two given space points in the time $`t`$, electric potentials that vanish at given surfaces and eigenvalues of the Laplace operator defined inside a finite domain (quantum billiards).

Because of the analogy with spinning tops, the case of trajectories of Hamiltonian systems is of particular interest here. The monodromy method, developed by Baranger and Davis baranger , was designed specifically to find periodic solutions of Hamiltonian systems with $`N`$ degrees of freedom. Xavier and de Aguiar xavier  extended the method to find non-periodic trajectories with any given combination of $`2N`$ position and momenta at initial and/or final times.

A widely used method to solve BVPs is the so called *shooting method* keller ; Ha . For a single second-order differential equation, the method consists in finding the proper ’velocity’ at the initial point so as to reach the desired ’target’ at the end point, similar to the shooting of a projectile. Examples of applications of this method to the Kirchhoff equations are the search for homoclinic orbits in reversible systems spence , heteroclinic orbits resembling tendril pervesions tyler  and the study of localized buckling modes of thin elastic filaments champ1 ; champ2 .

In the case of open rods some specific BVPs were solved recently. Károlyi and Domokos domokos , using symbolic dynamics, found global invariants for BVPs of elastic linkages, as natural discretization of continuous elastic beams, an old problem solved by Euler(see reference domokos ). Gottlieb and Perkins perkins  investigated spatially complex forms in a BVP governing the equilibrium of slender cables subjected to thrust, torsion and gravity. Also, the criteria of Tobias *et al.* tobias  was applied to linear segments subjected to *strong anchoring end conditions*, where not only the end points but also the tangent vector at the end points are held fixed. The dependence of DNA tertiary structure on end conditions was studied in tobias2 , where explicit expressions for equilibrium configurations were obtained for a specific case with symmetric end conditions.

Our aim in this paper is to develop a method to find equilibrium solutions of finite rods subjected to boundary conditions at their end points and with given load parameters. We emphasize that this is different from the approach in champ1 ; champ2 , where the authors use shooting methods to calculate localized buckling modes. These modes are treated as homoclinic solutions of the Kirchhoff equations, corresponding to infinite rods that become asymptotically straight in the infinite. Our objective is to find equilibrium solutions for finite rods subjected to boundary conditions at both ends.

Our method is an adaptation of the monodromy matrix method for non-periodic trajectories xavier  to the hamiltonian formulation of the static Kirchhoff equations. We work with the Kirchhoff equations directly in Euler angles, instead of using the Cartesian position and tangent vectors champ1 ; champ2 . The difficulty in working with Euler angles is that the variables that one wants to hold fixed, the spatial position of the filament end points, are not the variables appearing in the differential equations describing the rod. However, the number of differential equations to be solved is much smaller in these variables. Using a symmetry of the Hamiltonian, we end up with only two independent equations to solve.

One of the motivations of this work is its possible biological applications as, for example, the study of single DNA molecules manipulated by optical traps wuite ; meiners ; vincent , and the DNA loops between multiprotein structures (such as the *lac* repressor-operator complex) schleif ; mahadevan .

This work is organized as follows. In Sec. II we review the Kirchhoff equations and, in Sec. III, their hamiltonian formulation. In Sec. IV we describe our method for solving the BVP. The monodromy method enters as part of the solution, and proves to be a very efficient tool. In Sec. V we give numerical examples, calculating the three dimensional configuration of rods with different sets of load parameters and end positions. We also discuss the existence of solutions as function of the load parameters. In Sec. VI, motivated by the repeated sequences of base-pairs commonly found in DNA molecules, we consider rods with periodically varying Young modulus. We compare the configurations of these non-homogeneous rods against their homogeneous counterparts, fixing the same end points and mechanical parameters. In Sec. VII we summarize our conclusions.

## II The Kirchhoff Equations

The Kirchhoff model describes the dynamics of inextensible thin elastic filaments within the approximation of linear elasticity theory dill . They result from the application of Newton’s laws of mechanics to a thin rod, and consist of two equations describing the balance of linear and angular momentum plus a constitutive relationship of linear elasticity theory, relating moments to strains. The Kirchhoff model assumes that the filament is thin and weakly bent (*i.e.* its cross-section radius is much smaller than its length and its curvature at all points). In this approximation it is possible to derive a one-dimensional theory where forces and moments are averaged over the cross-sections perpendicular to the central axis of the filament.

A thin tube can be described by a smooth curve $`𝐱`$ in the 3D space parametrized by the arclength $`s`$, and whose position depends on the time: $`𝐱=𝐱(s,t)`$. A *local orthonormal basis*, (or *director basis*) $`𝐝_i=𝐝_i(s,t)`$, $`i=1,2,3`$, is defined at each point of the curve, with $`𝐝_3`$ chosen as the tangent vector, $`𝐝_3=𝐱^{}`$. In this paper we shall use primes to denote differentiation with respect to $`s`$ and dots to denote differentiation with respect to time. The two orthonormal vectors, $`𝐝_1`$ and $`𝐝_2`$, lie in the plane normal to $`𝐝_3`$, for example along the principal axes of the cross section of the rod. These vectors are chosen such that $`\{𝐝_1,𝐝_2,𝐝_3\}`$ form a right-handed orthonormal basis for all values of $`s`$ and $`t`$. The space and time evolution of the director basis along the curve are controlled by *twist* and *spin equations*

$$𝐝_i^{}=𝐤\times 𝐝_i,\dot{𝐝}_i=𝝎\times 𝐝_ii=1,2,3$$
(1)
which follow from the orthonormality relations $`𝐝_i𝐝_j=\delta _{ij}`$. The components of $`𝐤`$ and $`𝝎`$ in the director basis are defined as $`𝐤=_{i=1}^3k_i𝐝_i`$ and $`𝝎=_{i=1}^3\omega _i𝐝_i`$. $`k_1`$ and $`k_2`$ are the components of the curvature and $`k_3`$ is the twist density of the rod. The solution of the twist and spin equations determines $`𝐝_3(s,t)`$, which can be integrated to give the space curve $`𝐱(s,t)`$.

Let the material points on the rod be labeled by

$$𝐗(s,t)=𝐱(s,t)+𝐫(s,t),$$
(2)
where

$$𝐫(s,t)=x_1𝐝_1(s,t)+x_2𝐝_2(s,t)$$
(3)
gives the position of the point on the cross section $`𝒮(s)`$, perpendicular to $`𝐱^{}(s)`$, with respect to the central axis. The total force $`𝐅=𝐅(s,t)`$ and the total moment $`𝐌=𝐌(s,t)`$ (with respect to the axis of the rod) on the cross section are defined by

$$𝐅=_{𝒮(s)}𝐩_s𝑑S.$$
(4)
$$𝐌=_{𝒮(s)}𝐫\times 𝐩_s𝑑S,$$
(5)
where $`𝐩_s`$ is the contact force per unit area exerted on the cross section $`𝒮(s)`$. In terms of the director basis we write $`𝐅=_{i=1}^3f_i𝐝_i`$ and $`𝐌=_{i=1}^3M_i𝐝_i`$.

In order to derive a set of equations describing the rod as a one-dimensional object, the rod is divided into thin disks of length $`\text{d}s`$ and cross section $`𝒮(s)`$. To each of these disks the conservation laws of linear and angular momentum are applied dill . The result is

$$𝐅^{}+_{𝒮(s)}𝐟_{ext}𝑑S=_{𝒮(s)}\rho _0\ddot{𝐗}𝑑S,$$
(6)
$$𝐌^{}+𝐱^{}\times 𝐅+_{𝒮(s)}𝐫\times 𝐟_{ext}𝑑S=_{𝒮(s)}\rho _0𝐫\times \ddot{𝐗}𝑑S.$$
(7)
where $`𝐟_{ext}`$ is an external force that will not be considered in our calculations ($`𝐟_{ext}=0`$ in what follows).

In this article we are interested only in the equilibrium solutions and, therefore, we shall drop the derivatives with respect to time. Assuming that the rod has a uniform circular cross section of area $`A`$, Eqs. (6) and (7) can be simplified to yield

$$𝐅^{}=0,$$
(8)
$$𝐌^{}+𝐝_3\times 𝐅=0$$
(9)
which are a set of six equations for 9 variables: $`𝐅`$, $`𝐌`$ and $`𝐤`$ (from which we determine $`𝐝_i`$). In order to close the system of equations we need a *constitutive relation* relating the local forces and moments (stresses) to the elastic deformations of the body (strains). In the linear theory of elasticity, and for a homogeneous elastic material, the stress is proportional to the deformation. The Young’s modulus $`E`$ and the Shear modulus $`\mu `$ characterize the elastic properties of the material. Therefore, it is possible to obtain, for small deformations, a constitutive relation for the moment. For a isotropic material, in the director basis, this relation is dill :

$$𝐌=EI\left(k_1k_1^u\right)𝐝_1+EI\left(k_2k_2^u\right)𝐝_2+2\mu I\left(k_3k_3^u\right)𝐝_3,$$
(10)
where $`I`$ is the principal moment of inertia of the cross section, $`k_i`$ are the components of the strain vector and $`k_i^u`$ are the components of the twist vector in the unstressed configuration. The case $`k_i^u=0`$ corresponds to a naturally straight and untwisted rod. We shall assume $`k_i^u=0`$.

Eqs. (8), (9) and (10) can be further simplified by the introduction of scaled variables:

$$\begin{array}{c}ssL,𝐅𝐅\frac{EI}{L^2},\\ \\ 𝐌𝐌\frac{EI}{L},𝐤𝐤\frac{1}{L}.\end{array}$$
(11)
In the new variables the rod has total length $`L=1`$. The static Kirchhoff equations become

$$𝐅^{}=0,$$
(12)
$$𝐌^{}+𝐝_3\times 𝐅=0,$$
(13)
$$𝐌=k_1𝐝_1+k_2𝐝_2+\mathrm{\Gamma }k_3𝐝_3.$$
(14)
where $`\mathrm{\Gamma }=2\mu /E`$ is an elastic parameter that does not affect the equilibrium solutions.

## III Hamiltonian Formulation

In order to construct a hamiltonian formulation of the Kirchhoff equations we first note that Eqs. (12)-(14) are integrable if $`E`$ and $`\mu `$ are constant nizzete . Eq. (12) shows that the tension $`𝐅`$ is constant. Let us choose the direction of the force as the $`z`$ direction:

$$𝐅=F𝐞_Z.$$
(15)
In analogy to the spinning top, the tension $`𝐅`$ corresponds to the gravity field $`m𝐠`$. Here, $`F`$ can be considered as an external parameter and not as a first integral. Substituting Eq. (15) in Eq. (13) and projecting along $`𝐞_Z`$ we get

$$𝐌^{}𝐞_ZM_Z^{}=0$$
(16)
which does represent a first integral. By projecting the Eq. (13) along $`𝐝_3`$ we obtain another integral, $`M_3`$, since

$$𝐌^{}𝐝_3M_3^{}=0.$$
(17)
Finally, it is also possible to show that the elastic energy per unit arclength

$$H=\frac{1}{2}𝐌𝐤+𝐅𝐝_3$$
(18)
is constant, i.e., $`H^{}=0`$. Therefore $`H`$ is the last integral.

The orthonormal Cartesian basis can be connected to the director basis by Euler angles with

$$𝐝_i=\underset{j=1}{\overset{3}{}}S_{ij}𝐞_j$$
(19)
where

$$S=\left(\begin{array}{ccc}\mathrm{cos}\theta \mathrm{cos}\varphi \mathrm{cos}\psi \mathrm{sin}\varphi \mathrm{sin}\psi & \mathrm{cos}\theta \mathrm{cos}\varphi \mathrm{sin}\psi +\mathrm{sin}\varphi \mathrm{cos}\psi & \mathrm{cos}\varphi \mathrm{sin}\theta \\ \mathrm{cos}\theta \mathrm{sin}\varphi \mathrm{cos}\psi \mathrm{cos}\varphi \mathrm{sin}\psi & \mathrm{cos}\theta \mathrm{sin}\varphi \mathrm{sin}\psi +\mathrm{cos}\varphi \mathrm{cos}\psi & \mathrm{sin}\varphi \mathrm{sin}\theta \\ \mathrm{sin}\theta \mathrm{cos}\psi & \mathrm{sin}\theta \mathrm{sin}\psi & \mathrm{cos}\theta \end{array}\right).$$
(20)
The static Kirchhoff equations (12)-(14) can then be written in terms of $`\theta `$, $`\varphi `$ and $`\psi `$. We get

$$\begin{array}{c}\theta ^{\prime \prime }(\psi ^{})^2\mathrm{sin}\theta \mathrm{cos}\theta +\mathrm{\Gamma }\psi ^{}(\varphi ^{}+\psi ^{}\mathrm{cos}\theta )\mathrm{sin}\theta =F\mathrm{sin}\theta \hfill \\ \psi ^{\prime \prime }\mathrm{sin}\theta +2\psi ^{}\theta ^{}\mathrm{cos}\theta \mathrm{\Gamma }\theta ^{}(\varphi ^{}+\psi ^{}\mathrm{cos}\theta )=0\hfill \\ \psi ^{\prime \prime }\mathrm{cos}\theta =\psi ^{}\theta ^{}\mathrm{sin}\theta \varphi ^{\prime \prime }\hfill \end{array}.$$
(21)
These equations can also be derived directly from Eqs. (16)-(18). In terms of the Euler angles the Hamiltonian becomes

$$H=\frac{P_\theta ^2}{2}+\frac{P_\varphi ^2}{2\mathrm{\Gamma }}+\frac{(P_\psi P_\varphi \mathrm{cos}\theta )^2}{2\mathrm{sin}^2\theta }+F\mathrm{cos}\theta ,$$
(22)
where

$$P_\theta =\theta ^{},$$
(23)
$$P_\varphi M_3=\mathrm{\Gamma }(\varphi ^{}+\psi ^{}\mathrm{cos}\theta ),$$
(24)
$$P_\psi M_Z=\psi ^{}\mathrm{sin}^2\theta +P_\varphi \mathrm{cos}\theta .$$
(25)
Eqs. (21) correspond to Hamilton’s equations $`P_\alpha ^{}=\frac{H}{\alpha ^{}}`$ and $`\alpha ^{}=\frac{H}{P_\alpha }`$ for $`\alpha =\theta ,\psi `$ or $`\varphi `$. We see immediately that $`P_\varphi `$ and $`P_\psi `$ are constants and that $`\theta `$ is the only independent variable.

The total elastic energy of the rod can be calculated by the integration of the Eq. (22):

$$E_T=_0^1H(s)𝑑s,$$
(26)
The energy is a function of $`P_\psi `$, $`P_\varphi `$ and $`F`$. It also depends on the initial conditions $`\theta (s=0)\theta _0`$ and $`P_\theta (s=0)P_{\theta 0}`$.

The procedure to construct equilibrium solutions for given constants $`P_\psi `$ and $`P_\varphi `$ and initial condition $`(\theta _0,P_{\theta 0})`$ is as follows: first we solve the equations $`P_\theta ^{}=\frac{H}{\theta ^{}}`$ and $`\theta ^{}=\frac{H}{P_\theta }`$ to obtain $`(\theta (s),P_\theta (s))`$. Second, using Eq. (25), we obtain $`\psi (s)`$. The solutions $`\theta (s)`$ and $`\psi (s)`$ are sufficient to construct the rod by integrating the tangent vector $`𝐝_3`$:

$$𝐱(s)=_0^s𝐝_3(s^{})𝑑s^{}.$$
(27)
Explicitly,

$$x(s)=_0^s\mathrm{sin}\theta (s^{})\mathrm{cos}\psi (s^{})𝑑s^{},$$
(28)
$$y(s)=_0^s\mathrm{sin}\theta (s^{})\mathrm{sin}\psi (s^{})𝑑s^{},$$
(29)
$$z(s)=_0^s\mathrm{cos}\theta (s^{})𝑑s^{},$$
(30)
where

$$\psi (s)=\psi _0+_0^s\frac{P_\psi P_\varphi \mathrm{cos}\theta (s^{})}{\mathrm{sin}^2\theta (s^{})}𝑑s^{}.$$
(31)
Substituting equation (31) in (28-(30) and re-arranging the terms we obtain, in matrix form,

$$\left(\begin{array}{c}x(s)\\ y(s)\\ z(s)\end{array}\right)=\left(\begin{array}{ccc}\mathrm{cos}\psi _0& \mathrm{sin}\psi _0& 0\\ \mathrm{sin}\psi _0& \mathrm{cos}\psi _0& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}x_0(s)\\ y_0(s)\\ z(s)\end{array}\right).$$
(32)
where $`x_0(s)`$ and $`y_0(s)`$ are the equations (28) and (29) for $`\psi _00`$. Therefore, it suffices to find the solution with $`\psi _0=0`$. The solutions for other values of $`\psi _0`$ are simple rotations of this basic solution.

## IV The Linearized Method

In this section we present our method for finding the configuration of finite rods subject to boundary conditions in the position of its initial and final points. In fact, since the Kirchhoff equations are invariant under space translations, we can always choose the initial point be the origin. As we saw in the previous section, equilibrium solutions for the static Kirchhoff equations depend only on two initial conditions, namely, $`\theta _0`$ and $`P_{\theta 0}`$. The third initial condition, $`\psi _0`$ corresponds to a rotation of this solution around the z-axis. The problem is then that of finding a solution that starts from the origin and ends at $`z_f`$ and at a distance $`r_f\sqrt{y_f^2+x_f^2}`$ from the z-axis. Our method is based on a series of small deformations made upon an initial solution of the Kirchhoff equations which, however, does not have the desired boundary values for $`r_f`$ and $`z_f`$. We call this initial solution the trial solution. The method consists in pushing the end point of the trial solution to the desired position, step by step. The basic idea, which is a type of shooting procedure, is to find a variation in the initial conditions so as to obtain the desired variation in the end point.

In order to do so, we shall employ a variation of the Monodromy Method baranger ; xavier , originally devised to calculate periodic solutions of chaotic Hamiltonian systems. As discussed in the Introduction, working with the Kirchhoff equations in Euler angles poses an extra difficulty on the already hard problem of satisfying boundary conditions: the variables to be held fixed, $`r_f`$ and $`z_f`$, are not the ones entering the equations of motion, namely, $`\theta `$, $`P_\theta `$ and $`\psi `$. The advantage is that we can find the solutions solving only two differential equations.

The rod can be obtained from the Euler angles the equations (28), (29) and (30). The Euler angles, in their turn, obey the equations

$$\theta ^{}=P_\theta $$
(33)
$$P_\theta ^{}=\frac{(P_\psi P_\varphi \mathrm{cos}\theta )P_\varphi }{\mathrm{sin}\theta }+\frac{(P_\psi P_\varphi \mathrm{cos}\theta )^2\mathrm{cos}\theta }{\mathrm{sin}^3\theta }+F\mathrm{sin}\theta $$
(34)
and

$$\psi (s)=_0^s\frac{P_\psi P_\varphi \mathrm{cos}\theta (s^{})}{\mathrm{sin}^2\theta (s^{})}𝑑s^{}.$$
(35)
If we integrate Eqs.(33)-(35) using the initial condition provided by the trial solution and further integrate Eqs.(28)-(30) with the resulting Euler angles, we get, of course, the trial rod. Variations in these initial conditions will produce variations in the rod configuration, and, in particular, in its end point. In what follows we shall construct an explicit relation between a small variation in the initial variables $`\theta _0`$ and $`P_{\theta 0}`$ and the rod’s end point, represented by $`r_f`$ and $`z_f`$. Explicitly, we shall find the matrix $`B`$ such that

$$\left(\begin{array}{c}\delta r_f\\ \delta z_f\end{array}\right)=B\left(\begin{array}{c}\delta \theta _0\\ \delta P_{\theta 0}\end{array}\right).$$
(36)
Once $`B`$ is obtained (and if it can be inverted) we can work our way from the trial solution, whose end point is at, say, $`r_t`$ and $`z_t`$, to the desired end point at $`r_f`$ and $`z_f`$, provided we do that in a series of small steps. In each step we use the previous solution as the trial input, pushing the rod’s end point slowly towards its final destination.

Using $`r_f=\sqrt{x_f^2+y_f^2}`$, the components of the matrix $`B`$ can be written as:

$$B_{11}=\frac{r_f}{\theta _0}=\frac{x_f}{r_f}\frac{x_f}{\theta _0}+\frac{y_f}{r_f}\frac{y_f}{\theta _0}$$
(37)
$$B_{12}=\frac{r_f}{P_{\theta 0}}=\frac{x_f}{r_f}\frac{x_f}{P_{\theta 0}}+\frac{y_f}{r_f}\frac{y_f}{P_{\theta 0}}$$
(38)
$$B_{21}=\frac{z_f}{\theta _0}$$
(39)
$$B_{22}=\frac{z_f}{P_{\theta 0}}.$$
(40)
From Eqs. (28)-(30) and (35) we find

$$\delta x_f=_0^1\mathrm{cos}\theta (s)\mathrm{cos}\psi (s)\delta \theta (s)𝑑s_0^1\mathrm{sin}\theta (s)\mathrm{sin}\psi (s)\delta \psi (s)𝑑s,$$
(41)
$$\delta y_f=_0^1\mathrm{cos}\theta (s)\mathrm{sin}\psi (s)\delta \theta (s)𝑑s+_0^1\mathrm{sin}\theta (s)\mathrm{cos}\psi (s)\delta \psi (s)𝑑s,$$
(42)
$$\delta z_f=_0^1\mathrm{sin}\theta (s)\delta \theta (s)𝑑s,$$
(43)
and

$$\delta \psi (s)=_0^sA(\theta (s^{}))\delta \theta (s^{})𝑑s^{},$$
(44)
where $`A(\theta )`$ is given by

$$A(\theta )=\frac{P_\varphi }{\mathrm{sin}\theta }\frac{2(P_\psi P_\varphi \mathrm{cos}\theta )\mathrm{cos}\theta }{\mathrm{sin}^3\theta }.$$
(45)
Finally, to find the relation between the variations $`\delta \theta (s)`$ and $`\delta P_\theta (s)`$ and their values at the initial point $`s=0`$, we consider small variations of Eqs.(33) and (34) around the trial solution :

$$\begin{array}{c}\delta \theta ^{}=\delta P_\theta ,\hfill \\ \delta P_\theta ^{}=C(\theta )\delta \theta ,\hfill \end{array}$$
(46)
where $`C(\theta )`$, given by

$$C(\theta )=P_\varphi ^2\frac{(P_\psi P_\varphi \mathrm{cos}\theta )(P_\psi 4P_\varphi \mathrm{cos}\theta )}{\mathrm{sin}^2\theta }\frac{3(P_\psi P_\varphi \mathrm{cos}\theta )^2\mathrm{cos}^2\theta }{\mathrm{sin}^4\theta }+F\mathrm{cos}\theta ,$$
(47)
is computed at the trial solution.

The solution to these linear equations can be written in matrix form as

$$\left(\begin{array}{c}\delta \theta (s)\\ \delta P_\theta (s)\end{array}\right)=\left(\begin{array}{cc}M_{11}(s)& M_{12}(s)\\ M_{21}(s)& M_{22}(s)\end{array}\right)\left(\begin{array}{c}\delta \theta _0\\ \delta P_0\end{array}\right),$$
(48)
where $`M`$ is the tangent matrix, satisfying $`M(0)=1`$. In the special case where the trial solution is periodic, $`M`$ is called the monodromy matrix.

Writing $`\delta \theta (s)`$ explicitly as

$$\delta \theta (s)=M_{11}(s)\delta \theta _0+M_{12}(s)\delta P_0,$$
(49)
and using Eqs.(41)-(44) we can readily obtain

$$\begin{array}{c}\frac{x_f}{\theta _0}=_0^1\mathrm{cos}\theta (s)\mathrm{cos}\psi (s)M_{11}(s)𝑑s_0^1\mathrm{sin}\theta (s)\mathrm{sin}\psi (s)_0^sA(\theta (s^{}))M_{11}(s^{})𝑑s^{}𝑑s\hfill \\ \frac{x_f}{P_{\theta 0}}=_0^1\mathrm{cos}\theta (s)\mathrm{cos}\psi (s)M_{12}(s)𝑑s_0^1\mathrm{sin}\theta (s)\mathrm{sin}\psi (s)_0^sA(\theta (s^{}))M_{12}(s^{})𝑑s^{}𝑑s\hfill \\ \frac{y_f}{\theta _0}=_0^1\mathrm{cos}\theta (s)\mathrm{sin}\psi (s)M_{11}(s)𝑑s+_0^1\mathrm{sin}\theta (s)\mathrm{cos}\psi (s)_0^sA(\theta (s^{}))M_{11}(s^{})𝑑s^{}𝑑s\hfill \\ \frac{y_f}{P_{\theta 0}}=_0^1\mathrm{cos}\theta (s)\mathrm{sin}\psi (s)M_{12}(s)𝑑s+_0^1\mathrm{sin}\theta (s)\mathrm{cos}\psi (s)_0^sA(\theta (s^{}))M_{12}(s^{})𝑑s^{}𝑑s\hfill \\ \frac{z_f}{\theta _0}=_0^1\mathrm{sin}\theta (s)M_{11}(s)𝑑s\hfill \\ \frac{z_f}{P_{\theta 0}}=_0^1\mathrm{sin}\theta (s)M_{12}(s)𝑑s\hfill \end{array}$$
(50)
and, therefore, the matrix $`B`$,

Since we linearized the equations of motion, we have to check if the new solution, starting from $`\theta _1=\theta _0+\delta \theta _0`$ and $`P_{\theta 1}=P_0+\delta P_{\theta 0}`$ generates a rod with the chosen final point, within a given precision. If the precision is not reached, we can use the newly computed solution as a new trial one, using again Eq. (36), now with ($`\delta r_f,\delta z_f`$) corresponding to the distance between the fixed end point and the end point of the previously computed rod. The process can be repeated until the desired accuracy is obtained.

Finally we note that the elements $`M_{ij}(s)`$ can be computed by solving the linear equations (46) with proper initial conditions. Indeed, setting $`\delta \theta _0=1`$ and $`\delta P_{\theta 0}=0`$, Eq. (48) gives $`M_{11}(s)=\delta \theta (s)`$ and $`M_{21}(s)=\delta P_\theta (s)`$. If, on the other hand, we set $`\delta \theta _0=0`$ and $`\delta P_{\theta 0}=1`$ we get $`M_{12}(s)=\delta \theta (s)`$ and $`M_{22}(s)=\delta P_\theta (s)`$. Therefore, $`M_{11}(s)`$ and $`M_{21}(s)`$ are solutions of the linearized Eqs. (46) with the initial conditions $`\delta \theta _0=1`$ and $`\delta P_{\theta 0}=0`$ and $`M_{12}(s)`$ and $`M_{22}(s)`$ are the solutions of the same equations with $`\delta \theta _0=0`$ and $`\delta P_{\theta 0}=1`$.

In many cases we might want to push the rod’s end-point to a position $`𝐫_f=(r_f,z_f)`$ far from that of the initial trial solution, $`𝐫_t=(r_t,z_t)`$. To do that we can divide the line connecting $`𝐫_t`$ to $`𝐫_f`$ into $`N`$ small segments and apply the linearized method $`N`$ times, moving a small distance at each step. The number of steps required will depend on the particular configuration and possibly on the stability of the rod. In all integrations presented in this paper, we used a fourth order Runge-Kuta method with fixed step. In all cases the distance between the end-point of the trial rod and the target position was divided into $`10`$ segments and the solution converged to the desired boundary condition with a precision of $`10^6`$ in each component $`r`$ and $`z`$.

Figure 1 shows a example of the method. We have chosen the following load parameters in scaled units (see Eq. (11): $`P_\psi =1.0`$, $`P_\varphi =1.0`$ and $`F=1.0`$. The desired end point is $`r_f=0.7`$ and $`z_f=0.0`$. We plot the trial, an intermediate and the converged rods together, in order to show the process of convergence from the trial to the desired solution. The trial solution was computed integrating the Kirchhoff equations starting from $`\theta _0=0.5`$ and $`P_{\theta 0}=1.0`$, which corresponds to a rod whose final point is $`r_f0.8`$ and $`z_f0.455`$. The intermediate solution was computed from $`\theta _00.216`$ and $`P_{\theta 0}2.257`$, which corresponds to a rod whose final point is $`r_f0.73`$ and $`z_f0.15`$. Finally, the initial conditions obtained for the converged solution are $`\theta _00.307`$ and $`P_{\theta 0}2.510`$.

## V Numerical Examples

The particular trial solution used in the previous numerical example converged smoothly to the chosen final position. In some cases, however, a given trial solution does not converge to its destination no matter how many intermediate steps are used to divide the line between $`𝐫_t`$ and $`𝐫_f`$. As we shall see, this problem has to do with the existence or not of solutions for a given $`𝐫_f`$. In scaled units, it is obvious that there are no solutions if $`\sqrt{r_f^2+z_f^2}>1`$. The restrictions are actually much stronger than this simple ’length rule’, and depend on the values of $`P_\varphi `$, $`P_\psi `$ and $`F`$. It might also happen that the solution for a given $`𝐫_f`$ does exist, but that the straight line connecting $`𝐫_t`$ to $`𝐫_f`$ passes through forbidden regions, hindering the convergence. In this section we investigate the space of possible solutions and give several examples of rods computed with our method.

Each initial condition $`\theta _0`$ and $`P_{\theta 0}`$ leads to an end point $`𝐫_f`$. The easiest way to map all possible final points is to scan the space of initial conditions. Therefore, for a fixed set of parameters $`P_\psi `$, $`P_\varphi `$ and $`F`$ we calculate $`𝐫_f=𝐫_f(\theta _0,P_{\theta 0})`$ and plot the resulting figure in the $`(r_f,z_f)`$ space. Points outside this region are unreachable by the rod. Changing the parameters, such as the tension, changes the region of possible solutions, including end-points that were not previously present and excluding others.

The results in this section are presented as follows: for each fixed set of $`P_\psi `$, $`P_\varphi `$ and $`F`$ we show the $`(r_f,z_f)`$-space of possible solutions. On the same plot we draw curves of constant $`D=\sqrt{r_f^2+z_f^2}`$ and, for each $`D`$ we compute the three dimensional configuration of a few rods.

In order to compare the rods, we always adjust the value of $`\psi _0`$ such that the rod ends in $`y_f=0`$ plane. $`\psi _0`$ is determined by:

$$\mathrm{tan}\psi _0=\frac{y_f^0}{x_f^0}$$
(51)
where $`x_f^0`$ and $`y_f^0`$ are the final values of $`x`$ and $`y`$ for the rod calculated with $`\psi _0=0`$.

In all cases tested our method converged with at least six significant figures to the previously defined final values of $`r_f`$ and $`z_f`$.

### V.1 $`P_\psi =0,P_\varphi =1`$ and $`F=0.1`$

Figure 2 shows the $`(r_f,z_f)`$ map for this case. The map was generated by varying the initial condition in the intervals $`0<\theta _0<\pi `$ ($`30`$ points) and $`5.0<P_{\theta 0}<5.0`$ ($`200`$ points). For larger values of $`P_{\theta 0}`$ the total elastic energy of the rod increases and the final points tend to concentrate in the region near the origin (data not shown). The full thick line represents the curve $`D=1.0`$, which is the natural limit for the solutions. But there is a large region inside the $`D=1`$ line where no solutions exist. We shall compare it with that of other load parameters later on. We also show the lines of constant $`D`$ for $`D=0.9`$ (full line), $`D=0.7`$ (dashed line) and $`D=0.4`$ (dotted-dashed line). It is also interesting to note a forbidden region centered around $`z_f0`$ and $`r_f0.25`$. From the sequence of points crossing each other, it is evident that there are two sets of initial conditions $`(\theta _0,P_{\theta 0})`$ that generate rods with the same end point. In general they correspond to rods that are above or below the $`z`$ axis, and we shall call them the ’up’ solution and the ’down’ solution respectively.

Figure 3 shows examples of the rods whose end points are marked with circles in Fig. 2. We show the up and down solutions for each final point. Figs. 3a and 3b show three rods each for $`D=0.4`$ and $`D=0.7`$. Fig. 3c shows 5 different rods for $`D=0.9`$.

Notice that when $`P_\varphi `$ or $`P_\psi `$ are zero, the hamiltonian (22), becomes symmetric under the transformation $`FF`$ and $`\theta \pi \theta `$. In these cases the ’up’ solution for a given $`F`$ and $`\theta _0`$ is identical to the ’down’ solution for $`F`$ and $`\pi \theta _0`$. When both $`P_\varphi `$ and $`P_\psi `$ are non-zero the symmetry disappears.

### V.2 $`P_\psi =0,P_\varphi =5`$ and $`F=1`$

The $`(r_f,z_f)`$ map for these parameters is shown in Fig. 4. It has a curious pattern of thin bulbs centered around the $`z_f=0`$ axis, that degenerate for small $`r_f`$. The only rods possible in this case are those that return close to the $`z=0`$ plane.

Keeping $`P_\psi =0`$ and $`F=1`$ but increasing $`P_\varphi `$ results in even thiner bulbs. Figures 5 (a) and (b) show the up and down solutions, respectively, for $`r_f=0.7`$, $`z_f=0`$ and $`P_\varphi =1`$, $`5`$ and $`10`$. Panels (c) and (d) show the up and down solutions for the same parameters, except for $`r_f=0.9`$. Notice that large values of $`P_\varphi `$ correspond to horizontal helical rods.

### V.3 $`P_\psi =P_\varphi `$ and $`F=1`$

Figure 6 shows the $`(r_f,z_f)`$ map for $`P_\psi =P_\varphi =5`$. It resembles the map on Fig. 4, only distorted towards $`z_f=1`$. For these values of the parameters the rod admits ’vertical’ configurations, as opposed to the ’horizontal’ configurations displayed in the previous case. Figures 7 show examples of rods with $`P_\psi =P_\varphi =5`$ and $`P_\psi =P_\varphi =10`$ for $`r_f=0.39`$ and $`r_f=0.5`$.

## VI Application to Non-homogeneous DNA

As a last application of our method we shall consider the equilibrium configurations of non-homogeneous rods. We restrict ourselves to the simplest case of periodic non-homogeneities in Young’s modulus. The motivation for this study is the fact that repeated (and therefore periodic) DNA sequences form a substantial fraction of all eukaryotic genomes jenny ; brian . The calculations presented here are based on the stiffness parameters recently computed for the 32 tri-nucleotide units from DNA data gromiha . Our goal is to understand how much the equilibrium configuration of a non-homogeneous rod differs from that of the homogeneous case when the rod is subject to fixed mechanical conditions.

Repetitive DNA is formed by nucleotide sequences of varying lengths and compositions. Repeated sequences, reaching up to 100 megabasepairs of length brian , appear to have little or no functional role, and are commonly regarded as “selfish” or “junk” DNA mclister . We shall use a simple periodic formula for the (scaled) Young’s modulus that covers most of the parameter interval spanned by the tri-nucleotides given in ref. gromiha :

$$E(s)=1+\alpha \mathrm{cos}\frac{2\pi }{}s.$$
(52)
$``$ is the period of the oscillations of the Young’s modulus along the DNA and $`\alpha `$ is a parameter that depends on the specific sequence being repeated and that can not be greater than $`0.66`$.

Eqs.(22)-(25) have to be slightly modified to include the non-constant Young’s and shear moduli. We obtain

$$H=\frac{P_\theta ^2}{2E(s)}+\frac{P_\varphi ^2}{2\mathrm{\Gamma }_0\mu (s)}+\frac{(P_\psi P_\varphi \mathrm{cos}\theta )^2}{2E(s)\mathrm{sin}^2\theta }+F\mathrm{cos}\theta ,$$
(53)
with

$$P_\theta =E(s)\theta ^{},$$
(54)
$$P_\varphi M_3=\mathrm{\Gamma }_0\mu (s)(\varphi ^{}+\psi ^{}\mathrm{cos}\theta ),$$
(55)
$$P_\psi M_Z=E(s)\psi ^{}\mathrm{sin}^2\theta +P_\varphi \mathrm{cos}\theta .$$
(56)
The solutions do not depend on $`\mu (s)`$, since it does not enter in the differential equations for $`\theta `$, $`P_\theta `$ and $`\psi `$. Notice that these equations are not integrable if $`\alpha 0`$. Although $`P_\psi `$ and $`P_\varphi `$ are still constants, the elastic energy per unit arclength is not.

The method for solving the BVP for a non-homogeneous rod is the following: consider a solution extending from the origin to $`𝐫_f`$ with $`\alpha =0`$ and initial conditions $`\theta _0`$ and $`P_{\theta _0}`$. Now integrate the Kirchhoff equations from the same initial conditions but using Eq.(52) with $`\alpha 0`$. This new solution, whose end point is $`𝐫_f+\delta 𝐫`$, can be used as a trial solution for the non-homogeneous rod. Using the method of section IV we push the rod back to $`𝐫_f`$.

Figure 8 shows the up and down solutions for load parameters $`F=1`$, $`P_\psi =0`$, $`P_\varphi =10`$, $`r_f=0.9`$ and $`z_f=0`$ for the homogeneous rod (black curve) and a non-homogeneous rod with $`\alpha =0.66`$ and $`=0.1`$ (gray curve).

Finally, Figure 9 shows the effect of changing the period of oscillations in Young’s modulus. We show the up and down solutions for non-homogeneous rods with the same load parameters of the previous figure and $`\alpha =0.66`$. The curves show rods for $`=0.1`$ (gray), $`=0.5`$ (thick black) and $`=0.65`$ (thin black). The changes in the three-dimensional shape of the rods are evident for these values of load parameters. The sensitivity of the three-dimensional shape to the base-pair sequences indicated that DNA repeats may have conformational roles.

## VII Conclusions

In this work we presented a general method to solve the boundary value problem (BVP) for finite Kirchhoff filaments. The method consists in making small changes to a known trial rod that satisfies the Kirchhoff equations but not necessarily the boundary conditions. We combine a shooting technique with the method of monodromy matrix to push the end point of the trial solution to the desired position, step by step. By linearizing the Kirchhoff equations we obtain an explicit relation between a variation of the initial conditions (expressed in terms of Euler angles) and the consequent variation of the rod’s end point.

The solutions of the BVP are limited by the physical constraints of the rod, such as the moments and tension. A sketch of the allowed end points can be constructed by integrating the Kirchhoff equations for a large number of initial conditions. The regions of possible end points form complex figures reflecting the non-linear character of the equations. The regions of existence of end points may serve as a guide to find the appropriate load parameters needed for a desired solution.

We presented several examples of rods with different end positions at different distances from the origin for various sets of load parameters. In all cases the method worked very well and the BVP was solved with at least six significant figures in the values of $`r_f`$ and $`z_f`$. We also applied the method to non-homogeneous, sequence-dependent, DNAs. We modeled pieces of repeated sequences by a sinusoidal oscillation of the Young’s modulus. We showed that the tri-dimensional structure of the DNA is indeed sensitive to the presence of such sequences, a property that has been considered before hogan  but studied only in terms of the DNA intrinsic curvature manning . The effect studied here may contribute to other sequence-dependent properties that affect the three-dimensional conformations of the DNA.

###### Acknowledgements.

This work was partially supported by the Brazilian agencies FAPESP, CNPq and FINEP.


format: text

text:
# Untitled Document

Chronology Protection in Generalized Gödel Spacetime

Wung-Hong Huang

Department of Physics

National Cheng Kung University

Tainan, 70101, Taiwan

The effective action of a free scalar field propagating in the generalized Gödel spacetime is evaluated by the zeta-function regularization method. From the result we show that the renormalized stress energy tensor may be divergent at the chronology horizon. This gives a support to the chronology protection conjecture.

Classification Number: 04.62.+v; 04.20.Gz

E-mail: whhwung@mail.ncku.edu.tw

Physical Review D60 (1999) 067505

Many people have considered the problem of time travel \[1-13\]. In recent, Hawking proposed the ”chronology protection conjecture”  which states that the laws of physics will always prevent a spacetime to form the closed timelike curves (CTCs). In his argument, the conjecture may be due to the divergence of the stress tensor near the chronology horizon where CTCs are beginning to form. Several model spacetimes have been considered: the wormhole spacetime with a ”time machine” \[1-4\], the Gott,s two-string spacetime, Grant space  , and the Misner spacetime \[8,12-14\]. In a number of papers it is concluded that the time machine is quantum unstable, as the stress tensor becomes divergent near the chronology horizon. However, some authors \[7-13\]claim that the stress tensor is finite everywhere and it poses a problem for the chronology protection.

For instance, despite of the divergent behavior shown in the first calculation by Hiscock and Konkowski , a specific-choice of the vacuum state in the Misner spacetime  could render the renormalized stress tensor vanish at the chronology horizon. Boulware, Tanaka, and Hiscock  also find that for sufficiently massive fields in Gott space and Grant space, respectively, the renormalized stress tensor could remain regular at the chronology horizon. In  it also argues that the renormalized stress tensor may be smoothed out by introducing the absorption material, so that the spacetime with a time machine may be stable against vacuum fluctuation. The investigation of the Roman ring of traversable wormhole in  also find that the vacuum fluctuation can be made arbitrarily small. However, in  it has argued that, even though a self-consistent vacuum state in the Misner spacetime could be constructed, it must ultimately be unstable against a slight perturbation. Therefore, the chronology protection conjecture may be right. Some people  have further argued that the solving the problem of chronology protection is impossible within the semiclassical theory of gravitation and it requires a fully developed theory of quantum gravity. Anyway, the general proof of the chronology protection conjecture has not yet existed and it is useful to study the quantum stress tensor of a physical matter near a chronology horizon in various spacetimes with CTCs.

In this paper we will consider the generalized Gödel spacetime

$$ds^2=dt^2+dx^22\sqrt{2}\frac{\mathrm{\Omega }}{\alpha }e^{\alpha x}dydt+\frac{1}{2}(14\frac{\mathrm{\Omega }^2}{\alpha ^2})e^{2\alpha x}dy^2+\alpha ^2dz^2.$$
$`(1)`$
Note that we have introduced a constant metric $`g_{zz}=a^2`$ in the above equation as this could help us to find the renormalized stress tensor in the later calculation. The case of $`\alpha =\sqrt{2}\mathrm{\Omega }`$ is a rotating spacetime originally proposed by Gödel  . As analyzed in , there are no CTCs in the spacetimes with $`\alpha ^24\mathrm{\Omega }^2`$ . But the spacetimes with $`\alpha ^2<4\mathrm{\Omega }^2`$ admit CTCs. We will use the zeta-function regularization method  to evaluate the renormalized effective action of a free scalar field propagating in the generalized Gödel spacetime. The result shows that the action is divergent at chronology horizon , i.e. $`\alpha ^2=4\mathrm{\Omega }^2`$ . From the result we can show that the stress energy may be divergent at the chronology horizon and this thus gives a new example to support the chronology protection conjecture.

In a system with a Hamiltonian H the renormalized effective action W evaluated by the zeta-function regularization method  can be expressed as

$$W=\frac{i}{2}[\zeta ^{}(0)+\zeta (0)ln(\mu ^2)],$$
$`(2)`$
where

$$\zeta (\nu )=Tr(H)^\nu .$$
$`(3)`$
To evaluate the above zeta function exactly we shall first find the eigenvalue of the Hamiltonian H. It is easy to see that the eigenvalue $`\mathrm{}`$ of the Hamiltonian H for a massive (m) scalar field with curvature coupling ($`\xi `$) to the generalized Gödel spacetime can be easily found through the same method  as that in the old Gödel spacetime. The method is that we first assume the following form of the eigenfunction

$$\mathrm{\Psi }(t,x,y,z)=h(x)exp(i\omega t+ik_yy+ik_zz).$$
$`(4)`$
Then the function h(x) will satisfy the Whittaker’s function and the condition to have an everywhere-bounded solution is that

$$\mathrm{}=a^2k_z^2+m^2+\xi R+(\frac{\alpha }{2})^2+(4\mathrm{\Omega }^2\alpha ^2)(n+\frac{1}{2})^2[\omega 2\mathrm{\Omega }ϵ(n+\frac{1}{2})]^2,$$
$`(5)`$
where the scalar curvature $`R=2(\alpha ^2\mathrm{\Omega }^2)`$ , $`ϵ=k_z/k_z`$ , and $`n=0,1,2,\mathrm{}`$. The zeta function can then be calculated by the formula

$$\zeta (\nu )=𝑑x^4\sqrt{g}𝑑k_z𝑑\omega \underset{n=0}{\overset{\mathrm{}}{}}(\mathrm{})^\nu .$$
$`(6)`$
Note that in  we have evaluated the finite-temperature effective potential for a scalar field in the old Gödel spacetime, also with the help of zeta-function regularization method and thus the equation (6).

Now, denoting $`M^2=m^2+\xi R+(\alpha /2)^2`$ and rotating $`\omega `$ through $`\pi /2`$ in the complex plane to become $`i\omega _E`$  we have the relation

$$\zeta (\nu )=𝑑x^4\sqrt{g}𝑑k_z𝑑\omega _E[a^2k_z^2+\omega _E^2+M^2+(4\mathrm{\Omega }^2\alpha ^2)(n+\frac{1}{2})^2]^\nu $$

$$=ia\pi (\nu 1)^1𝑑x^4\sqrt{g}\underset{n=0}{\overset{\mathrm{}}{}}[M^2+(4\mathrm{\Omega }^2\alpha ^2)(n+\frac{1}{2})^2]^{1\nu }.$$
$`(7)`$
Next, using the relation

$$\underset{n=0}{\overset{\mathrm{}}{}}F[(n+1/2)^2]=\underset{n=1}{\overset{\mathrm{}}{}}F[(n/2)^2]\underset{n=1}{\overset{\mathrm{}}{}}F[n^2],$$
$`(8)`$
in which F is an arbitrary function, we have

$$\zeta (\nu )=i\pi a(\nu 1)^1𝑑x^4\sqrt{g}\underset{n=1}{\overset{\mathrm{}}{}}[M^2+(4\mathrm{\Omega }^2\alpha ^2)(n/2)^2]^{1\nu }$$

$$i\pi a(\nu 1)^1𝑑x^4\sqrt{g}\underset{n=1}{\overset{\mathrm{}}{}}[M^2+(4\mathrm{\Omega }^2\alpha ^2)n^2]^{1\nu }.$$
$`(9)`$
To proceed, let us quote a formula

$$\underset{n=1}{\overset{\mathrm{}}{}}(An^2+B)^s=\frac{1}{2}A^{1/2}B^{1/2s}\pi \mathrm{\Gamma }(s1/2)/\mathrm{\Gamma }(s)+\frac{1}{2}(A+B)^s_0^1(Ax^2+B)^s𝑑x+$$

$$i_0^{\mathrm{}}[[A(1+ix)^2+B]^s[A(1ix)^2+B]^s](e^{2\pi x}1)^1𝑑x.$$
$`(10)`$
Then, after the calculations the zeta function becomes

$$\zeta (\nu )=i\nu 𝑑x^4\sqrt{g}\frac{7}{12}a\pi \frac{1}{\sqrt{4\mathrm{\Omega }^2\alpha ^2}}[m^2+\xi R+(\alpha /2)^2]^{3/2}+O(\nu ^2)+$$

$$(finitetermswhen4\mathrm{\Omega }^2approachesto\alpha ^2).$$
$`(11)`$
Using this result and from (2) we see that

$$<T_{zz}>=\frac{2}{\sqrt{g}}\frac{\delta W}{\delta g_{zz}}=\frac{7}{12}\pi \frac{1}{\sqrt{4\mathrm{\Omega }^2\alpha ^2}}[m^2+\xi R+(\alpha /2)^2]^{3/2}+$$

$$(finitetermswhen4\mathrm{\Omega }^2approachesto\alpha ^2),$$
$`(12)`$
after letting $`a=\mu =1`$. We thus see that, at least, one component of the renormalized stress energy tensor is divergent at the chronology horizon and the investigation thus gives a new support to the chronology protection conjecture. The evaluation of all components of the stress energy tensor is more difficult and left to be found.

Finally, let us mention that the spacetime considered in this paper dose not like those will form CTCs after the cosmological evolution \[1-13\]. Thus our result dose not directly relate to the chronology protection conjecture initially suggested therein. However, as a geometry with the causality violation, such as Gödel spacetime, can be taken as a part of a globally causal spacetime . Thus our investigation dose suggest that the physical law of requiring the finite stress tensor will forbid a geometry with the causality violation to be enclosed in a region which maintains the required condition of global causality. This means that the physical law will protect a would-be advanced civilization to create a noncausal geometry to violate the chronology protection conjecture.

1. M. S. Morris, K. S. Thorne and U. turtsever, Phys. Rev. Lett. 61. 1446 (1988).
2. V. P. Frolov and I. D. Novikov, Phys. Rev. D 42, 1057 (1990).
3. V. P. Frolov, Phys. Rev. D 43, 3878 (1991).
4. S. W. Kim and K. S. Thorne, Phys. Rev. D 43, 3929 (1991).
5. J. R. Gott, Phys. Rev. Lett. 66. 1126 (1991).
6. S. Hawking, Phys. Rev. D 46, 603 (1992).
7. J. D. E. Grant, Phys. Rev. D 47, 2388 (1993).
8. D. G. Boulware, Phys. Rev. D 46, 4421 (1992); T. Tanaka and W. A. Hiscock, Phys. Rev. D 52, 4503 (1995).
9. G. Klinkhammer, Phys. Rev. D 46, 3388 (1992); S. V. Krasnikov, ibid. 54, 7322 (1996).
10. M. Lyutikov, Phys. Rev. D 49, 4041 (1994); D. Laurence, ibid. 50, 4957 (1994); M. Visser, ibid. 55, 5212 (1997).
11. S. V. Sushkov, Class. Quantum Grav. 12, 1685 (1995); ibid. 14 , 523 (1997).
12. M. J. Cassidy, Class. Quantum Grav. 14 , 3031 (1997); ibid. 15 , 2257 (1998)
13. L-X Li, Quantum Grav. 13 , 2563 (1996); L-X Li and J. R. Gott, Phys. Rev. Lett. 80. 2980 (1998).
14. W. A. Hiscock and D. A. Konkowski, Phys. Rev. D 26, 1225 (1982).
15. M. J. Reboucas and J. Tiomno, Phys. Rev. D 28, 1251 (1983) ; M. J. Reboucas and A. F. F. Teixeira, Phys. Rev. D 34, 2985 (1988).
16. K. Gödel, Rev. Mod. Phys. 21, 447 (1949).
17. N. D. Birrell and P. C. W. Davies, Quantum Field in Curved Space (Cambridge University Press, Cambridge, England, 1982).
18. B. Mashhoon, Phys. Rev. D 11, 2679 (1975); W. A. Hiscock, ibid. 26, 1497 (1982).
19. W. H. Huang, Finit Temperature Cosmological Phase Transition in a Rotational Spacetimes,Class. Quantum Grav. 8, 1471 (1991);
W. H. Huang, One-loop Effective Action on Rotational Spacetimes: $`\zeta `$ \- Function Regularization and Swinger Perturbative Expansion, Ann. Phys. 254, 69 (1997).
20. J. Toms, Phys. Rev. D 21, 2805 (1980); L. H. Ford, ibid. 21, 933 (1980).
21. M. Novello and M. C. Motta da Silva, Phys. Rev. D 49, 825 (1994)


text:
# Time Variations in the Scale of Grand Unification

## I Introduction

Recent observations of quasar absorption lines have been interpreted as indicating variation in the fine structure constant on cosmological time scales Webb:2000mn ; Murphy:2002jx ,

$$\frac{\delta \alpha }{\alpha }=(0.57\pm 0.10)\times 10^5\mathrm{at}z0.23.7.$$
(1)
In the context of grand unified theories (GUTs), such a time variation in $`\alpha `$ may be related to time variations in other gauge couplings and, in particular, to variations in the strong coupling constant $`\alpha _3`$ or, equivalently, in the QCD scale $`\mathrm{\Lambda }_{\text{QCD}}`$. Such related variations affect various observables in interesting ways Calmet:2001nu ; Langacker:2001td ; Dent:2001ga ; Flambaum:2002de ; Ichikawa:2002bt ; Chacko:2002mf ; Calmet:2002ja ; Olive:2002tz .

Most previous works have assumed that the variation in $`\alpha `$ is induced by a variation in the unified coupling at the Planck scale, $`\alpha _\mathrm{U}(M_{\mathrm{Pl}})`$ Calmet:2001nu ; Langacker:2001td ; Dent:2001ga ; Flambaum:2002de ; Ichikawa:2002bt ; Olive:2002tz . The motivation for such studies comes mainly from string theories, where the value of $`\alpha _\mathrm{U}(M_{\mathrm{Pl}})`$ is determined by the value of the dilaton and other fields. We are here interested in another possibility, namely that the variation in $`\alpha `$ is induced by a variation in the GUT scale $`M_\mathrm{U}`$, with $`M_{\mathrm{Pl}}`$ and $`\alpha _\mathrm{U}(M_{\mathrm{Pl}})`$ held fixed. In ref. Dine:2002se , it is argued that if grand unification (in the sense of unification of the gauge interactions in a four-dimensional gauge group) arises within string theory, the GUT scale is likely to be a modulus.111Of course, in string theory, coupling unification does not necessarily require a unified gauge group. Such a situation can also arise in the framework of field theories, where symmetry breaking scales are determined by vacuum expectation values (VEVs) of scalar fields, but coupling constants are not. In either case, understanding why the modulus has the requisite properties is a difficult problem, which we will not address. We will show that the relation between $`\delta \alpha /\alpha `$ and $`\delta \mathrm{\Lambda }_{\text{QCD}}/\mathrm{\Lambda }_{\text{QCD}}`$ is highly model dependent. In particular, in contrast to the case that both are induced by $`\delta \alpha _\mathrm{U}/\alpha _\mathrm{U}0`$, it may happen that $`\delta \mathrm{\Lambda }_{\text{QCD}}/\mathrm{\Lambda }_{\text{QCD}}\delta \alpha /\alpha `$ (though this is not the generic case).

We are aware of two previous works that are closely related to ours. In ref. Calmet:2002ja , the consequences of variations of $`M_\mathrm{U}`$ have been analyzed. It was assumed, however, that $`\alpha _\mathrm{U}(M_\mathrm{U})`$ remains fixed. It is then necessary that $`\alpha _\mathrm{U}(M_{\mathrm{Pl}})`$ also varies in a correlated way. We think that such a scenario is less plausible than the one that we investigate. In ref. Chacko:2002mf , the consequences for the low energy couplings of variations in the mass of particles have been analyzed. Conceptually, our work takes a similar direction. In their discussion of grand unified theories, ref. Chacko:2002mf  examine a set of models in which only the masses of a subset of GUT mass particles vary in time; in particular, only the masses of color-neutral, charged particles shift. In contrast, our basic assumption is that all heavy masses (of GUT gauge bosons, fermions and scalars) are proportional to a single breaking scale, $`M_\mathrm{U}`$, and it is the variation of this scale that we study.222To understand the difference, one may think of the consequences of a variation in the VEV of the single Higgs doublet of the Standard Model, which would shift all masses, or a shift in the VEV of a hypothetical Higgs triplet, which would shift only neutrino masses.

In this work we focus on the resulting relations between the low energy couplings. We do not concern ourselves here with the origin of the scale variation, nor with the connection to higher scale physics. Discussions of these topics can be found in Campbell:1994bf ; Chacko:2002mf . In addition, there are fundamental issues we will also not address. As pointed out in Banks:2001qc , it is difficult to understand, in the framework of local quantum field theory, how couplings could vary by so large an amount, without an enormous variation in the vacuum energy. This problem is not significantly different if it is changes in the unification scale which are responsible for this variation than if the coupling at the Planck or unification scale changes.

The plan of this paper is as follows. In section II we present the formalism for the variation of low energy couplings induced by variations in the scale of threshold corrections. In section III we obtain the relations between the variation of relevant cosmological observables and argue that, within our framework, they depend on a single parameter. In section IV we apply the formalism to the experimental data. We use it to constrain the parameter of section II. In section V we test specific models of grand unified theories with this constraint. We show that a variation in the GUT scale is unlikely to explain the claimed variation in $`\delta \alpha /\alpha `$. In section VI we explain how our formalism and results apply much more generically, to any model where there is a time variation in a scale where threshold corrections take place. We also mention some further subleading effects which might be non-negligible. We conclude in section VII.

## II Formalism

We assume that the variation in couplings is due to a variation in a single intermediate scale, $`M_\mathrm{U}`$, where GUT threshold corrections take place. The gauge couplings $`\alpha _i(Q)`$ at a scale $`Q<M_\mathrm{U}`$ are given to one loop by the RGE solution

$$\alpha _i^1(Q)=\alpha _\mathrm{U}^1(M_\mathrm{U})+\frac{b_i}{2\pi }\mathrm{ln}\left(\frac{Q}{M_\mathrm{U}}\right),$$
(2)
where $`\alpha _\mathrm{U}(M_\mathrm{U})`$ is the unified coupling at the scale $`M_\mathrm{U}`$ and $`b_i`$ is the beta function coefficient (between $`M_\mathrm{U}`$ and $`Q`$) for $`\alpha _i`$. The scale $`M_\mathrm{U}`$ is related to the mass scale of heavy degrees of freedom that are integrated out, particularly the heavy GUT gauge bosons. Since the particle content of the theory below and above $`M_\mathrm{U}`$ is different, the RGE coefficient $`b_i`$ changes at $`M_\mathrm{U}`$. In particular, one obtains

$$\alpha _\mathrm{U}^1(M_\mathrm{U})=\alpha _\mathrm{U}^1(M_{\mathrm{Pl}})+\frac{b_\mathrm{U}}{2\pi }\mathrm{ln}\left(\frac{M_\mathrm{U}}{M_{\mathrm{Pl}}}\right),$$
(3)
where $`b_\mathrm{U}`$ is the beta function coefficient for the unified group. Combining (2) and (3) we write

$$\alpha _i^1(Q)=\alpha _\mathrm{U}^1(M_{\mathrm{Pl}})+\frac{b_\mathrm{U}}{2\pi }\mathrm{ln}\left(\frac{M_\mathrm{U}}{M_{\mathrm{Pl}}}\right)+\frac{b_i}{2\pi }\mathrm{ln}\left(\frac{Q}{M_\mathrm{U}}\right).$$
(4)
Now we allow the GUT scale, $`M_\mathrm{U}`$, to change by an amount of $`\delta M_\mathrm{U}`$, while holding $`\alpha _\mathrm{U}(M_{\mathrm{Pl}})`$ and $`M_{\mathrm{Pl}}`$ fixed. The resulting variation in the low scale couplings can be calculated from (4) (a similar derivation, for $`\alpha _{\text{EM}}`$, appears in Chacko:2002mf ):

$$\frac{\delta \alpha _i(Q)}{\alpha _i(Q)}=\alpha _i(Q)\frac{\mathrm{\Delta }b_i}{2\pi }\frac{\delta M_\mathrm{U}}{M_\mathrm{U}},$$
(5)
where

$$\mathrm{\Delta }b_ib_\mathrm{U}b_i.$$
(6)
Using (5), we can now relate the variation in low scale parameters to the variation $`\delta M_\mathrm{U}`$. We focus on two parameters, the electromagnetic coupling $`\alpha `$ and the QCD scale $`\mathrm{\Lambda }_{\text{QCD}}`$:

$$\alpha ^1=\frac{5}{3}\alpha _1^1(0)+\alpha _2^1(0),$$
(7)
$$\mathrm{\Lambda }_{\text{QCD}}=M_Z^{23/27}m_b^{2/27}m_c^{2/27}\mathrm{exp}\left(\frac{2\pi }{9\alpha _3(M_Z)}\right).$$
(8)
Neglecting effects of threshold corrections below $`M_\mathrm{U}`$ (which are discussed later), we obtain the relations:

$$\frac{\delta \alpha }{\alpha }=\frac{\alpha }{2\pi }\left(\frac{5}{3}\mathrm{\Delta }b_1+\mathrm{\Delta }b_2\right)\frac{\delta M_\mathrm{U}}{M_\mathrm{U}},$$
(9)
$$\frac{\delta \mathrm{\Lambda }_{\text{QCD}}}{\mathrm{\Lambda }_{\text{QCD}}}=\frac{\mathrm{\Delta }b_3}{9}\frac{\delta M_\mathrm{U}}{M_\mathrm{U}}.$$
(10)
Finally, we use (9) and (10) to relate the variation of the two parameters:

$`{\displaystyle \frac{\delta \mathrm{\Lambda }_{\text{QCD}}}{\mathrm{\Lambda }_{\text{QCD}}}}`$ $`=`$ $`{\displaystyle \frac{2\pi }{9\alpha }}{\displaystyle \frac{\mathrm{\Delta }b_3}{\frac{5}{3}\mathrm{\Delta }b_1+\mathrm{\Delta }b_2}}{\displaystyle \frac{\delta \alpha }{\alpha }}`$ (11)
$``$ $`R{\displaystyle \frac{\delta \alpha }{\alpha }}.`$
The factor $`R\frac{95.7\mathrm{\Delta }b_3}{\frac{5}{3}\mathrm{\Delta }b_1+\mathrm{\Delta }b_2}`$ plays a crucial role in our analysis. As we show later, its value is highly GUT-model dependent.

## III Observables

In this section, we use (11) to obtain expressions for the variation of several cosmological observables. The primordial 4He abundance is given by Kolb:1990vq ; Bergstrom:1999wm

$$Y_4\frac{2(n_n/n_p)_{\text{NS}}}{1+(n_n/n_p)_{\text{NS}}},$$
(12)
where the ratio of the number density of the neutron to that of the proton at nucleosynthesis is given by

$$(n_n/n_p)_{\text{NS}}0.8\times e^{Q/T_D}.$$
(13)
Here $`T_D0.8`$ MeV is the decoupling temperature and $`Q`$ is the mass difference between the neutron and the proton. Since $`T_D`$ is a function of $`M_{\mathrm{Pl}}`$ and $`G_F`$, it is independent of the Standard Model gauge couplings and we therefore neglect its variation. The mass difference $`Q`$, on the other hand, depends on both $`\alpha `$ and $`\mathrm{\Lambda }_{\text{QCD}}`$ Gasser:1982ap :

$$Q2.05\text{MeV}+C\alpha \mathrm{\Lambda }_{\text{QCD}},$$
(14)
where $`C`$ is a dimensionless order one parameter. Assuming variations in $`\alpha `$ and $`\mathrm{\Lambda }_{\text{QCD}}`$, the value of $`Q`$ is shifted \[eq. (14)\] from its present-day value of $`Q1.29\text{ MeV}`$, and consequently $`(n_n/n_p)_{\text{NS}}`$ would be shifted \[eq. (13)\] from its ‘standard’ value of $`(n_n/n_p)_{\text{NS}}\frac{1}{7}`$. We thus obtain the change in $`Y_4`$ compared to its value calculated with the present values of the coupling constants:

$$\begin{array}{cc}\hfill \frac{\delta Y_4}{Y_4}& =\frac{1}{1+(n_n/n_p)_{\text{NS}}}\frac{2.05\text{MeV}Q}{T_D}\left(\frac{\delta \alpha }{\alpha }+\frac{\delta \mathrm{\Lambda }_{\text{QCD}}}{\mathrm{\Lambda }_{\text{QCD}}}\right)\hfill \\ & 0.8\left(1+R\right)\frac{\delta \alpha }{\alpha }.\hfill \end{array}$$
(15)
Two other relevant observables are $`X`$ and $`Y`$. The $`X`$ parameter gives the ratio between the hyperfine 21cm neutral hydrogen absorption transition to an optical resonance transition Cowie:1995sz :

$$X\alpha ^2g_p\left(\frac{m_e}{m_p}\right),$$
(16)
where $`m_e`$ is the electron mass, $`m_p`$ is the proton mass, and $`g_p`$ is the proton gyromagnetic moment. The $`Y`$ parameter determines the molecular hydrogen transition in the early universe Potekhin:1998mf ; Ivanchik:2001ji ,

$$Y\frac{m_p}{m_e}.$$
(17)
The proton mass is proportional, at first order, to $`\mathrm{\Lambda }_{\text{QCD}}`$ Langacker:2001td :

$$\frac{\delta m_p}{m_p}\frac{\delta \mathrm{\Lambda }_{\text{QCD}}}{\mathrm{\Lambda }_{\text{QCD}}}=R\frac{\delta \alpha }{\alpha }.$$
(18)
Assuming that the variations of $`m_e`$ and $`g_p`$ are negligible, we obtain:

$$\frac{\delta X}{X}2\frac{\delta \alpha }{\alpha }\frac{\delta m_p}{m_p}\left(2R\right)\frac{\delta \alpha }{\alpha },$$
(19)
$$\frac{\delta Y}{Y}R\frac{\delta \alpha }{\alpha }.$$
(20)
## IV Constraining $`R`$ with Experimental Data

Values for the observed variation of various observables were obtained through cosmological and nuclear studies. Table 1 summarizes these values (the details can be found in the quoted references). Using the formalism of the previous section, we relate the allowed range for each observable in Table 1 to an (in general, $`R`$-dependent) allowed variation in $`\alpha `$. The results are also given in Table 1.

For a given value of $`R`$, one obtains from Table 1 the implied variation of $`\alpha `$ at different $`z`$’s. The most interesting result, however, arises from requiring consistency of the three observables related to the intermediate redshift values. Note that $`\delta Y/Y`$ and $`\delta X/X`$ are consistent with each other for any value of $`R`$, because, at the $`2\sigma `$ level, they allow $`\delta \alpha /\alpha =0`$, for which the value of $`R`$ becomes irrelevant. The constraints become nontrivial when one considers the claimed signal of $`\delta \alpha /\alpha `$ \[eq. (1)\]. The requirement of consistency between the observed $`\delta Y/Y`$ at $`z2.3\text{}3`$ and the observed $`\delta \alpha /\alpha `$ at $`z0.2\text{}3.7`$ favors a range in $`R`$ that is mostly negative in sign and not very large in magnitude, $`R=𝒪(10)`$. The requirement of consistency between the observed $`\delta X/X`$ at $`z1.8`$ and the observed $`\delta \alpha /\alpha `$ at $`z0.2\text{}3.7`$ favors a range in $`R`$ that is mostly positive in sign and, again, not very large in magnitude, $`R=𝒪(+3)`$. The most interesting results follows from fitting simultaneously all three observables. The result of a $`\chi ^2`$ analysis is that the following range of $`R`$ is favored:

$$1R+6.$$
(21)
Values of $`R`$ outside this range have a probability that is $`5\%`$ to yield the observed results. (Note that even the best fit point, with $`R+2`$, has only a probability $`15\%`$ to yield the observed results.)

This result can be understood qualitatively: When $`|R|`$ is large, a value for $`\delta \alpha /\alpha `$ implies an even larger value for $`\delta \mathrm{\Lambda }_{\text{QCD}}/\mathrm{\Lambda }_{\text{QCD}}`$. This larger variation dominates the shift of nucleon masses. Since $`Y`$ is proportional to the proton mass while $`X`$ is inversely propotional to it, their relative variations must be of opposite signs, $`\delta Y\delta X0`$. The data, however, indicate that their variations are both biased to be positive. Consequently, only one of them can be consistent with a non-zero variation of $`\alpha `$. (Which one depends on the sign of $`R`$.) This situation is avoided only in case that the variation in $`\mathrm{\Lambda }_{\text{QCD}}`$ does not dominate $`\delta X/X`$ or, in other words, in case of a small $`R`$.

To summarize, eq. (21) gives the consistent range for the parameter $`R`$ in theories where the variation of the low energy coupling constants is dominated by threshold effects at a single varying scale.

## V Implications for various GUT models

In this section we apply our results to specific models. We focus on the minimal supersymmetric standard model (MSSM) embedded in a grand unified theory (GUT), in which the varying scale is the scale of breaking of the unified group, $`M_\mathrm{U}`$. In general, this is the mass scale of the GUT gauge supermultiplets, and we assume that it characterizes also the masses of all the heavy chiral supermultiplets. The one loop beta function coefficients below $`M_\mathrm{U}`$ are those of the MSSM:

$$b_1=\frac{33}{5},b_2=1,b_3=3.$$
(22)
Above the GUT scale, there is a single coefficient common to all gauge couplings, $`b_\mathrm{U}`$. Thus, the ratio $`R`$ of eq. (11) is given by

$$R=\frac{2\pi }{9\alpha }\frac{b_\mathrm{U}+3}{\frac{8}{3}b_\mathrm{U}12}.$$
(23)
Before we discuss specific GUT models, let us make some general observations:

1. For asymptotically large (positive or negative) values of $`b_\mathrm{U}`$, we obtain

$$|b_\mathrm{U}|1R\frac{\pi }{12\alpha }+36.$$
(24)
Large negative values of $`b_\mathrm{U}`$ are not relevant, as demonstrated by the very conservative lower bounds set by examining the contributions from the gauge supermultiplet and from the minimal matter representations that are necessary to accommodate the MSSM fields:

$$b_\mathrm{U}[SU(5)]>8,b_\mathrm{U}[SO(10)]>17.$$
(25)
On the other hand, realistic models often have a large, positive $`b_\mathrm{U}`$ and predict $`R=𝒪(40)`$.

2. Negative values of $`R`$ are achieved only for a very limited range of $`b_\mathrm{U}`$:

$$R<0\mathrm{for}3<b_\mathrm{U}<+\frac{9}{2}.$$
(26)
3. For a hypothetical value of $`b_\mathrm{U}=\frac{9}{2}`$, $`\alpha `$ does not vary (to the approximation that we use):

$$b_\mathrm{U}=\frac{9}{2}|R|\mathrm{}.$$
(27)
This result explains the change of sign of $`R`$ between models with $`b_\mathrm{U}=4`$ and $`b_\mathrm{U}=5`$, and the very large magnitude of $`R`$ in these models, $`R502(+574)`$ for $`b_\mathrm{U}=4(5)`$.

4. For $`b_\mathrm{U}=b_3=3`$, $`\mathrm{\Lambda }_{\text{QCD}}`$ does not vary (to the approximation that we use):

$$b_\mathrm{U}=3R=0.$$
(28)
In general, the $`𝒪(100)`$ factor in eq. (11), which is due to the exponential dependence of $`\mathrm{\Lambda }_{\text{QCD}}`$ on $`\alpha _3`$, leads to a large $`R`$. In order to overcome this factor and obtain $`|R|𝒪(1)`$, the threshold correction of $`\alpha _3`$ needs to be highly suppressed, i.e., $`\mathrm{\Delta }b_30.01\mathrm{\Delta }b_{1,2}`$. For MSSM GUT models with varying $`M_{\mathrm{GUT}}`$, this would imply $`3.5b_\mathrm{U}2.5`$. In particular, when $`b_\mathrm{U}`$ is shifted by one unit from the special value of $`3`$, one already obtains rather large values of $`|R|`$, $`R=+4.2(5.5)`$ for $`b_\mathrm{U}=4(2)`$. We learn that, if $`R`$ is to be within the range of eq. (21), the only acceptable integer values for $`b_\mathrm{U}`$ are $`3`$ and $`4`$.

We now move on to discuss several examples of specific GUT models. The $`b_\mathrm{U}`$ coefficients and the resulting $`R`$-factors of various models are given in Table 2. Each model is defined by its representation content. (Note that, for example, $`n_\mathrm{𝟓}`$ gives the number of SU(5) fundamentals plus the number of antifundamentals.) All SU(5) models have the quark and lepton fields in three generations of $`\mathrm{𝟏𝟎}+\overline{\mathrm{𝟓}}`$ and the Higgs fields related to electroweak symmetry breaking in $`\mathrm{𝟓}+\overline{\mathrm{𝟓}}`$. The minimal model has, in addition, a single $`\mathrm{𝟐𝟒}`$ to break $`SU(5)G_{\mathrm{SM}}`$. This model is, however, excluded by combining constraints from coupling unification and from proton decay Murayama:2001ur . (It may still be viable with a special flavor structure of the supersymmetric mixing matrices Bajc:2002bv .) A viable, though fine-tuned, model can be constructed by adding another $`\mathrm{𝟓}+\overline{\mathrm{𝟓}}`$ pair Hisano:1992ne ; Murayama:2001ur . To naturally induce doublet-triplet splitting, the $`\mathrm{𝟐𝟒}`$ has to be replaced by, at least, $`\mathrm{𝟓𝟎}+\overline{\mathrm{𝟓𝟎}}+\mathrm{𝟕𝟓}`$ Masiero:1982fe ; Grinstein:1982um ; Altarelli:2000fu .

All SO(10) models have the quark and lepton fields in three generations of $`\mathrm{𝟏𝟔}`$. The minimal Higgs sector has $`\mathrm{𝟒𝟓}+\mathrm{𝟏𝟔}+\overline{\mathrm{𝟏𝟔}}`$ to break $`SO(10)G_{\mathrm{SM}}`$ and a single $`\mathrm{𝟏𝟎}`$ for electroweak symmetry breaking Babu:1998wi . Note that we assume here that the breaking of $`SO(10)G_{\mathrm{SM}}`$ is done in one step, and that the VEVs of all GUT-breaking Higgs fields vary together. The investigation of a two-step breaking, $`SO(10)SU(5)G_{\mathrm{SM}}`$, with independent variation of each scale, will be presented in future work.333In case that the scale of $`SO(10)SU(5)`$ breaking varies but the scale of $`SU(5)G_{\mathrm{SM}}`$ breaking does not, the result is effectively a variation in the unified coupling constant at the $`SU(5)`$ breaking scale. This result shows how our formalism can be applied also to models of time variation in $`\alpha _\mathrm{U}`$: we should take $`\mathrm{\Delta }b_1=\mathrm{\Delta }b_2=\mathrm{\Delta }b_3`$, which gives $`R36`$. Models where doublet-triplet splitting is achieved naturally require an additional $`\mathrm{𝟏𝟎}`$ and either an additional $`\mathrm{𝟏𝟔}+\overline{\mathrm{𝟏𝟔}}`$ pair Barr:1997hq  or a more complicated Higgs sector that includes several $`\mathrm{𝟒𝟓}`$ and $`\mathrm{𝟓𝟒}`$ multiplets Babu:1993we . In each of these models, one could have a $`\mathrm{𝟏𝟐𝟔}+\overline{\mathrm{𝟏𝟐𝟔}}`$ pair in place of a $`\mathrm{𝟏𝟔}+\overline{\mathrm{𝟏𝟔}}`$ pair. Such a replacement increases $`b_\mathrm{U}`$ by $`66`$ and gives $`R=𝒪(40)`$.

Other models of interest are reviewed in, for example, ref. Mohapatra:1999vv  and lead to similar results. We conclude that most GUT models give $`|R|1`$. A remarkable exception is given by models with $`b_\mathrm{U}=3`$. This possibility is realized in the minimal $`SU(5)`$ model which, without a very special flavor structure, is, however, phenomenologically excluded. Other models with $`b_\mathrm{U}=3`$ can be constructed within the framework of larger groups. There is, however, little motivation to consider such a scenario (which would require fine-tuning for the doublet-triplet splitting) for the sole purpose of deriving the desired $`b_\mathrm{U}`$. We conclude that, within the framework of realistic GUT models that avoid fine-tuning, our mechanism predicts $`|R|1`$.

When this result is confronted with eq. (21), one is led to the conclusion that supersymmetric GUT models with a varying GUT scale are unlikely candidates to explain a variation in $`\alpha `$ of the size given in eq. (1).

## VI Generalization and Limitations

Our formalism and many of our results have a much broader applicapibility than the $`M_\mathrm{U}`$-variation in GUT models that we have focussed on. The formalism applies to any variation in a single intermediate scale where threshold corrections take place, $`\mu _{\mathrm{th}}`$. Suppose that the couplings are held fixed at some fixed high energy scale $`\mu _0`$. Define $`b_i`$ to be the one-loop beta function coefficients between $`\mu _{\mathrm{th}}`$ and a low scale $`Q`$, while $`b_i^{}`$ are the corresponding coefficients between $`\mu _0`$ and $`\mu _{\mathrm{th}}`$. We now define, as a generalization of eq. (6),

$$\mathrm{\Delta }b_ib_i^{}b_i.$$
(29)
With this definition, eq. (9) for $`\delta \alpha /\alpha `$, eq. (10) for $`\delta \mathrm{\Lambda }_{\text{QCD}}/\mathrm{\Lambda }_{\text{QCD}}`$, and eq. (11) for the ratio $`R`$, still apply.444One has to be careful about the normalization of $`\alpha _1`$ in non-GUT models. For example, if we take the conventional definition of $`U(1)_Y`$ in the Standard Model, the factors of $`\frac{5}{3}`$ in eqs. (7), (9) and (11) have to be omitted, while for eq. (22) one should use $`b_1=11`$. Similarly, the constraints on $`\delta \alpha /\alpha `$ in Table 1 and the allowed range for $`R`$ in eq. (21) hold in these more general circumstances.

An important comment is now in order: In principle, the variation $`\delta \mu _{\mathrm{th}}`$ may induce a variation in all couplings at lower scales. In the most general case, this includes a variation in the masses of all the particles which are lighter than $`\mu _{\mathrm{th}}`$. A variation of these masses could be induced by two sources:

1. The parameters that determine a particle mass, such as Yukawa couplings, may have threshold corrections at $`\mu _{\mathrm{th}}`$. In this case, a variation analogous to (5) results. To calculate this effect, one should use the (linearized) beta function of these parameters below and above the scale $`\mu _{\mathrm{th}}`$ .
2. The variation in the gauge couplings, $`\alpha _i`$, below $`\mu _{\mathrm{th}}`$, leads to variations in the low scale masses. To calculate this effect, one should use the RGE of the mass parameters below $`\mu _{\mathrm{th}}`$.

When running the gauge couplings from $`\mu _{\mathrm{th}}`$ to some low scale $`Q`$, all particles with masses $`Q<m<\mu _{\mathrm{th}}`$ are integrated out. For example, in an MSSM GUT model, with $`Qm_b`$ and $`\mu _{\mathrm{th}}=M_{\mathrm{GUT}}`$, we integrate out all supersymmetric particles at a scale $`M_{\mathrm{SUSY}}`$ as well as heavy SM particles ($`t,Z,W^\pm `$). This integration out would usually introduce additional threshold corrections, which have small effects on the parameters at $`m_b`$.

We are mainly interested in small variations of $`\alpha `$ and $`\mathrm{\Lambda }_{\text{QCD}}`$. As explained above, the variation $`\delta \mu _{\mathrm{th}}`$ induces variations in all intermediate scales below $`\mu _{\mathrm{th}}`$. But these variations induce, by threshold corrections at these scales, further variation in gauge couplings. The whole process may be pictured as a ‘chain reaction’, with a final effect that could be significant, even if the individual variations at each intermediate scale are small. We postpone the detailed study of this effect to future work.

Here, we would only like to emphasize that when $`R`$ is large, as is the case in most of the GUT models that we have investigated, the variation of $`\mathrm{\Lambda }_{\text{QCD}}`$ is the dominant effect at low energies, and it is well justified to ignore variations in couplings other than the gauge couplings. For small $`R`$, however, the modification due to threshold corrections of other parameters can become important. Thus, when we think of $`R=0`$ models, one should not conclude that $`\alpha `$ would vary with $`\mathrm{\Lambda }_{\text{QCD}}`$ remaining strictly fixed. Very likely, $`\mathrm{\Lambda }_{\text{QCD}}`$ varies due to the effects of lower thresholds and of higher loops. But the effect is expected to be small \[$`R=𝒪(1)`$\] and our qualitative conclusions would not change.

## VII Conclusions

A variation in a physical scale in which threshold corrections take place leads to variations in low scale observables. Current experimental data constrain the relations between the threshold corrections of the three different gauge couplings. In particular, consistency between the claimed variation $`\delta \alpha /\alpha `$ \[eq. (1)\] and the allowed ranges for $`\delta X/X`$ and $`\delta Y/Y`$ \[Table 1\] give

$$1R\frac{2\pi }{9\alpha }\frac{\mathrm{\Delta }b_3}{\frac{5}{3}\mathrm{\Delta }b_1+\mathrm{\Delta }b_2}+6,$$
(30)
where $`\mathrm{\Delta }b_i`$ is the difference of the one-loop beta function coefficient above and below the varying scale. In other words, the difference between the beta function coefficients of the strong coupling $`\alpha _3`$ above and below the varying scale, should be suppressed by a factor of $`𝒪(0.01)`$ compared with those of the other two couplings, $`\alpha _1`$ and $`\alpha _2`$.

We focussed our investigation on the framework where the MSSM is embedded in GUT models and the varying scale is $`M_{\mathrm{GUT}}`$. We demonstrated that the relation between the variation of $`\mathrm{\Lambda }_{\text{QCD}}`$ and that of $`\alpha `$, that is, the ratio $`R=\frac{\delta \mathrm{\Lambda }_{\text{QCD}}/\mathrm{\Lambda }_{\text{QCD}}}{\delta \alpha /\alpha }`$, is highly model dependent. For example, we can construct SU(5) models with $`R`$ as high as $`𝒪(+600)`$ and as low as $`𝒪(500)`$, and a model where $`R=0`$. (The latter is the minimal SU(5) model, and the former have, respectively, eight and seven $`\mathrm{𝟓}+\overline{\mathrm{𝟓}}`$ pairs added to the minimal model.) We argued, however, that it is difficult to obtain consistency with the requirement (30) in realistic supersymmetric GUT models in our framework.

###### Acknowledgements.

The visit of M.D. to the Weizmann Institute was supported by the Albert Einstein Minerva Center for Theoretical Physics and by the United States - Israel Binational Science Foundation (BSF). Y.N. is supported by the Israel Science Foundation founded by the Israel Academy of Sciences and Humanities.


text:
# DESY 02–147 Edinburgh 2002/09 LU-ITP 2002/018 Structure functions near the chiral limitPoster at Lattice 2002, Cambridge, USA.

## 1 INTRODUCTION

Understanding the structure of hadrons in terms of quark and gluon constituents (partons), in particular how quarks and gluons provide the binding and spin of the nucleon, is one of the outstanding problems in particle physics.

Moments of parton distribution functions, such as $`x`$ and $`g_A=\mathrm{\Delta }u\mathrm{\Delta }d`$, are benchmark calculations in lattice QCD. To compare the lattice results with experiment, one must extrapolate the data from the lowest calculated quark mass to the physical value. Up to now, most results are for quark masses of about the strange quark mass and larger. A naive, linear extrapolation of $`x`$ in the quark mass overestimates the experimental number by $``$ 40 % . While most results are for quenched QCD, where one might expect that $`x`$ is larger than the experimental value, recent unquenched results  indicate that this problem remains.

It has been argued that a linear extrapolation must fail because it omits non-analytic structure associated with chiral symmetry breaking . Indeed, chiral perturbation theory  suggests a large deviation of $`x`$ and $`g_A`$ from linearity as the pion (quark) mass tends to zero:

$`x_{NS}`$ $`=`$ $`x_{NS}^0\left(1{\displaystyle \frac{3g_A^2+1}{(4\pi f_\pi )^2}}m_\pi ^2\mathrm{ln}\left({\displaystyle \frac{m_\pi ^2}{m_\pi ^2+\mathrm{\Lambda }^2}}\right)\right)`$ (1)
$`+O(m_\pi ^2),`$
$`g_A`$ $`=`$ $`g_A^0\left(1{\displaystyle \frac{2g_A^2+1}{(4\pi f_\pi )^2}}m_\pi ^2\mathrm{ln}\left({\displaystyle \frac{m_\pi ^2}{m_\pi ^2+\mathrm{\Lambda }^2}}\right)\right)`$ (2)
$`+O(m_\pi ^2),`$
where $`\mathrm{\Lambda }`$ is a phenomenological parameter. Equation (1) fits both the lattice data and the experimental value . However, $`\mathrm{\Lambda }`$ is not yet determined by the lattice data. To constrain this parameter, and to perform an accurate extrapolation based solely on lattice results, data at smaller quark masses are crucial.

## 2 THE SIMULATION

To narrow the gap between the ‘chiral’ regime and the lowest calculated quark mass, we have started simulations at quark masses corresponding to about twice the physical pion mass. The calculations are done for Wilson fermions in the quenched approximation. Improved Wilson fermions are known to suffer from exceptional configurations, which would forbid such a calculation at currently accessible couplings.

We work at $`\beta =6.0`$. Our present data sample consists of:

| $`\beta `$ | $`\kappa `$ | Volume | Configurations |
| --- | --- | --- | --- |
| 6.0 | 0.1515 | $`16^3\mathrm{\hspace{0.17em}32}`$ | O(5000) |
| 6.0 | 0.1530 | $`16^3\mathrm{\hspace{0.17em}32}`$ | O(5000) |
| 6.0 | 0.1550 | $`16^3\mathrm{\hspace{0.17em}32}`$ | O(5000) |
| 6.0 | 0.1550 | $`24^3\mathrm{\hspace{0.17em}32}`$ | 220 |
| 6.0 | 0.1558 | $`24^3\mathrm{\hspace{0.17em}32}`$ | 220 |
| 6.0 | 0.1563 | $`24^3\mathrm{\hspace{0.17em}32}`$ | 220 |
| 6.0 | 0.1563 | $`32^3\mathrm{\hspace{0.17em}48}`$ | 50 |
| 6.0 | 0.1566 | $`32^3\mathrm{\hspace{0.17em}48}`$ | 100 |

The $`24^3\mathrm{\hspace{0.17em}32}`$ and $`16^3\mathrm{\hspace{0.17em}32}`$ lattices are from our previous runs . On all our lattices $`m_\pi L>4`$, so that finite volume effects may be expected to be small. We did not see any at $`\kappa =0.1563`$. We are currently increasing our statistics at the smaller quark masses.

## 3 HADRON MASSES

We first looked at the chiral behavior of pion, $`\rho `$ and nucleon masses. We did not find compelling evidence for non-analytic behavior in any of the three cases. In Fig. 1 we plot the nucleon mass as a function of the pion mass. Quenched chiral perturbation theory (qCPT) predicts

$$m_N=m_N^0+C_{1/2}m_\pi +C_1m_\pi ^2+C_{3/2}m_\pi ^3$$
(3)
with

$$C_{1/2}=\frac{3\pi }{2}\left(D3F\right)^2\delta 0.5.$$
(4)
Empirically one finds, to a very high precision,

$$m_N=\sqrt{m_N^{\mathrm{0\hspace{0.33em}2}}+C_1m_\pi ^2},$$
(5)
which is an analytic function in the quark mass. In Fig. 1 we fit (3) and (5) to the data. Both fits are hardly distinguishable in the range of the data. Fit (3) gives $`C_{1/2}=0.7(1)`$ and $`\chi ^2=0.8`$. Note that the coefficient $`C_{1/2}`$ comes out positive instead of negative as it should, according to qCPT. There is some preference for the single parameter analytic fit (5), which gives $`\chi ^2=0.4`$. What comes as a surprise is that the resulting physical nucleon masses (denoted by open circles in Fig. 1) still differ by more than 10 %. We conclude that pion masses of 300 MeV are not small enough to sufficiently constrain the fit function, so as to give unambiguous results.

## 4 STRUCTURE FUNCTIONS

Let us now turn to the moments of non-singlet, unpolarized and polarized structure functions. We consider $`x_{NS}`$ first. On the lattice one does not compute $`x_{NS}`$ directly but

$$R=x_{NS}m_N,$$
(6)
which can be interpreted as the fraction of the proton’s mass carried by the $`u`$ quark minus that carried by the $`d`$ quark. Because chiral perturbation theory does not provide much guidance at present quark masses, we seek a phenomenological extrapolation. In Fig. 2 we plot $`R`$ as a function of the pion mass. The data points have been properly renormalized and converted to renormalization group invariant (RGI) numbers following . We find that the data lie precisely on a straight line, all the way from $`m_\pi 1\text{GeV}`$ down to our smallest pion mass at $`m_\pi =300\text{MeV}`$. The slope comes out to be $`0.4`$, which is close to the value $`1/3`$ expected in the heavy quark limit. A linear extrapolation to the physical pion mass, $`R=R^0+R_{1/2}m_\pi `$, gives a value for $`x_{NS}^{RGI}`$ that is much closer to the experimental number than previous results. A fit of (1) to $`x_{NS}`$ still does not constrain $`\mathrm{\Lambda }`$. We find $`\mathrm{\Lambda }=360_{360}^{+120}\text{MeV}`$. Our present data cannot distinguish between the fit function (1) and $`(R^0+R_{1/2}m_\pi )/\sqrt{m_N^{\mathrm{0\hspace{0.33em}2}}+C_1m_\pi ^2}`$.

Let us now consider $`g_A`$, the axial vector coupling of the nucleon. This quantity describes the fraction of the proton’s spin carried by the $`u`$ quark minus that carried by the $`d`$ quark. In Fig. 3 we plot $`g_A`$ as a function of the pion mass. Again, the data lie on a straight line, and we do not see any sign of non-analytic behavior as suggested by (2).

## 5 CONCLUSION

It appears that at pion masses of 300 MeV we are not yet sensitive to the predictions of chiral perturbation theory. Whether the situation will improve on finer lattices has to be seen.

The numerical calculations were performed on the APE100 at NIC (Zeuthen) and on the Cray T3E at NIC (Jülich).


text:
# The NEXUS Model

( )

1 SUBATECH, Universit de Nantes – IN2P3/CNRS – EMN, Nantes, France
2 Moscow State University, Institute of Nuclear Physics, Moscow, Russia
3 Institute f. Kernphysik, Forschungszentrum Karlsruhe, Karlsruhe, Germany
4 Physics Department, New York University, New York, USA

a Invited speaker

at the International Workshop on the Physics of the Quark Gluon Plasma,

Palaiseau, France, September 4-7, 2001

Abstract

> The interpretation of experimental results at RHIC and in the future also at LHC requires very reliable and realistic models. Considerable effort has been devoted to the development of such models during the past decade, many of them being heavily used in order to analyze data.
>
> There are, however, serious inconsistencies in the above-mentioned approaches. In this paper, we will introduce a fully self-consistent formulation of the multiple-scattering scheme in the framework of a Gribov-Regge type effective theory.

With the start of the RHIC program to investigate nucleus-nucleus collisions at very high energies, there is an increasing need of computational tools in order to provide a clear interpretation of the data. The situation is not satisfactory in the sense that there exists a nice theory (QCD) but we are not able to treat nuclear collisions strictly within this framework, and on the other hand there are simple models, which can be applied easily but which have no solid theoretical basis. A good compromise is provided by effective theories, which are not derived from first principles, but which are nevertheless self-consistent and calculable. A candidate seems to be the Gribov-Regge approach, and – being formally quite similar – the eikonalized parton model. Here, however, severe inconsistencies occur, as discussed in .

As a solution of the above-mentioned problems, we present a new approach which we call “Parton-based Gribov-Regge Theory”: we have a consistent treatment for calculating (topological) cross sections and particle production considering energy conservation in both cases; in addition, we introduce hard processes in a natural way.

The basic guideline of our approach is theoretical consistency. We cannot derive everything from first principles, but we use rigorously the language of field theory to make sure not to violate basic laws of physics, which is easily done in more phenomenological treatments (see discussion above).

Let us first introduce some conventions. We denote elastic two body scattering amplitudes as $`T_{22}`$ and inelastic amplitudes corresponding to the production of some final state $`X`$ as $`T_{2X}`$ (see fig. 1).

As a direct consequence of unitarity on may write the optical theorem $`2\mathrm{I}\mathrm{m}T_{22}=_x`$$`(T_{2X})`$$`(T_{2X})^{}`$. The right hand side of this equation may be literally presented as a “cut diagram”, where the diagram on one side of the cut is $`(T_{2X})`$ and on the other side $`(T_{2X})^{}`$, as shown in fig. 2.

So the term “cut diagram” means nothing but the square of an inelastic amplitude, summed over all final states, which is equal to twice the imaginary part of the elastic amplitude.

Before coming to nuclear collisions, we need to discuss somewhat the structure of the nucleon, which may be studied in deep inelastic scattering – so essentially the scattering of a virtual photon off a nucleon, see fig. 3.

Let us assume a high virtuality photon. It is known that this photon couples to a high virtuality quark, which is emitted from a parton with a smaller virtuality, the latter on again being emitted from a parton with lower virtuality, and so on. We have a sequence of partons with lower and lower virtualities, the closer one gets to the proton. At some stage, some “soft scale” scale $`Q_0^2`$ must be reached, beyond which perturbative calculations are no longer valid. So we have some “unknown object” – indicated by “?” in fig. 3 – between the first parton and the nucleon. In order to proceed, one may estimate the squared mass of this “unknown object”, and one obtains doing simple kinematics the value $`Q_0^2/x`$, where $`x`$ is the momentum fraction of the first parton relative to the nucleon. Therefore, sea quarks or gluons, having typically small $`x`$, lead to large mass objects – which we identify with soft Pomerons, whereas valence quarks lead to small mass objects, which we simply ignore. So we have two contributions, as shown in fig. 3: a sea contribution, where the sea quark or gluon is emitted from a soft Pomeron, and a valence contribution, where the valence quark is one of the three quarks of the nucleon. The precise microscopic structure of the soft Pomeron not being known, it is parameterized in the usual way a a Regge pole.

Elementary nucleon-nucleon scattering can now be considered as a straightforward generalization of photon-nucleon scattering: one has a hard parton-parton scattering in the middle, and parton evolutions in both directions towards the nucleons. We have a hard contribution $`T_{\mathrm{hard}}`$ when the the first partons on both sides are valence quarks, a semi-hard contribution $`T_{\mathrm{semi}}`$ when at least on one side there is a sea quark (being emitted from a soft Pomeron), and finally we have a soft contribution, when there is no hard scattering at all (see fig. 4).

We have a smooth transition from soft to hard physics: at low energies the soft contribution dominates, at high energies the hard and semi-hard ones, at intermediate energies (that is where experiments are performed presently) all contributions are important.

Let us consider nucleus-nucleus ($`AB`$) scattering. In the Glauber-Gribov approach , the nucleus-nucleus scattering amplitude is defined by the sum of contributions of diagrams, corresponding to multiple elementary scattering processes between parton constituents of projectile and target nucleons. These elementary scatterings are exactly as discussed above, namely the sum of soft, semi-hard, and hard contributions: $`T_{22}=T_{\mathrm{soft}}+T_{\mathrm{semi}}+T_{\mathrm{hard}}`$. A corresponding relation holds for the inelastic amplitude $`T_{2X}`$. We use the above definition of a cut elementary diagram, which is graphically represented by a vertical dashed line, whereas the elastic amplitude is represented by an unbroken line:

$`=T_{22}`$,   $`=_X(T_{2X})`$$`(T_{2X})^{}`$.

This is very handy for treating the nuclear scattering model. We define the model via the elastic scattering amplitude $`T_{ABAB}`$ which is assumed to consist of purely parallel elementary interactions between partonic constituents, as discussed above. The amplitude is therefore a sum of terms as the one shown in fig. 5.

One has to be careful about energy conservation: all the partonic constituents (lines) leaving a nucleon (blob) have to share the momentum of the nucleon. So in the explicit formula one has an integration over momentum fractions of the partons, taking care of momentum conservation. Having defined elastic scattering, inelastic scattering and particle production is practically given, if one employs a quantum mechanically self-consistent picture. Let us now consider inelastic scattering: one has of course the same parallel structure, just some of the elementary interactions may be inelastic, some elastic. The inelastic amplitude being a sum over many terms – $`T_{ABX}=_iT_{ABX}^{(i)}`$ – has to be squared and summed over final states in order to get the inelastic cross section, which provides interference terms $`_X(T_{ABX}^{(i)})(T_{ABX}^{(j)})^{}`$, which can be conveniently expressed in terms of the cut and uncut elementary diagrams, as shown in fig. 6. So we are doing nothing more than following basic rules of quantum mechanics.

Of course a diagram with 3 inelastic elementary interactions does not interfere with the one with 300, because the final states are different. So it makes sense to define classes $`K`$ of interference terms (cut diagrams) contributing to the same final state, as all diagrams with a certain number of inelastic interactions and with fixed momentum fractions of the corresponding partonic constituents. One then sums over all terms within each class $`K`$, and obtains for the inelastic cross section

$$\sigma _{AB}(s)=d^2b\underset{K}{}\mathrm{\Omega }^{(s,b)}(K)$$

where we use the symbolic notation $`d^2b=d^2b_0d^{2A}b_A\rho (b_A)d^{2B}b_B\rho (b_B)`$ which means integration over impact parameter $`b_0`$ and in addition averaging over nuclear coordinates for projectile and target. The variable $`K`$ is characterized by $`AB`$ numbers $`m_k`$ representing the number of cut elementary diagrams for each possible pair of nucleons and all the momentum fractions $`x^+`$ and $`x^{}`$ of all these elementary interactions (so $`K`$ is a partly discrete and partly continuous variable, and $``$ is meant to represent $``$). This is the really new and very important feature of our approach: we keep explicitly the dependence on the longitudinal momenta, assuring energy conservation at any level of our calculation.

The calculation of $`\mathrm{\Omega }`$ actually very difficult and technical, but it can be done and we refer the interested reader to the literature .

The quantity $`\mathrm{\Omega }^{(s,b)}(K)`$ can now be interpreted as the probability to produce a configuration $`K`$ at given $`s`$ and $`b`$. So we have a solid basis for applying Monte Carlo techniques: one generates configurations $`K`$ according to the probability distribution $`\mathrm{\Omega }`$ and one may then calculate mean values of observables by averaging Monte Carlo samples. The problem is that $`\mathrm{\Omega }`$ represents a very high dimensional probability distribution, and it is not obvious how to deal with it. We decided to develop powerful Markov chain techniques  in order to avoid to introduce additional approximations.

To summarize: We provide a new formulation of the multiple scattering mechanism in nucleus-nucleus scattering, where the basic guideline is theoretical consistency. We avoid in this way many of the problems encountered in present day models. We also introduce the necessary numerical techniques to apply the formalism in order to perform practical calculations.

This work has been funded in part by the IN2P3/CNRS (PICS 580) and the Russian Foundation of Basic Researches (RFBR-98-02-22024).