text:
# 1 Introduction

## 1 Introduction

The space-time structure at Planck scales $`L_P(G\mathrm{}/c^3)^{1/2}`$ will be drastically affected by the quantum gravitational effects and it is generally believed that $`L_P`$ acts as a physical cutoff for space-time intervals. The two main approaches to quantum gravity, Superstring theory and Loop quantum gravity, incorporates the fundamental length scale by considering extended structures, rather than point particles, as fundamental blocks.

The existence of a fundamental length implies that processes involving energies higher than Planck energies will be suppressed, and the ultraviolet behavior of the theory will be improved. This could arise naturally in String theories; several other models also incorporate a Planck length cut-off in a suitable manner to improve the ultra violet behaviour of the theory ?. One direct consequence of such improved behavior will be that the Feynman propagator(in momentum space) will acquire damping factor for energies larger than Planck energy. However, this propagator – which arises in the standard formulation of quantum field theory – does not take into account of the existence of any fundamental length in the space-time. On the other hand, such a fundamental length scale was introduced into the Feynman propagator in a Lorentz invariant manner by invoking the “principle of path-integral duality” ?. According to this postulate, the weightage given for a path in the path integral should be invariant under the transformation $`L_P^2/`$, where $``$ is the length of the path and the fundamental length scale $`L_P`$ is assumed to be of the order of the Planck length $`(G\mathrm{}/c^3)^{1/2}`$.\[In this paper, when we say that the fundamental length is $`L_P`$, we actually mean that it is of $`𝒪(L_P)`$.\] Padmanabhan ? has shown by rigorous evaluation of the path integral by lattice techniques that the effect of the duality principle is to modify the weightage given to a path of proper time $`s`$ from $`\mathrm{exp}(im^2s)`$ to $`\mathrm{exp}[i(m^2sL_P^2/s)]`$, where $`m`$ is associated with the mass of the particle. For example, the Feynman propagator for a free scalar field of mass $`m`$, propagating in the flat $`(3+1)`$ space-time, in Schwinger’s proper time formalism, is described by the integral

$$G_F(x,x^{})=\frac{1}{(4\pi i)^2}_0^{\mathrm{}}\frac{ds}{s^2}\mathrm{exp}(im^2s)\mathrm{exp}[i(xx^{})^2/4s].$$
(1)
The duality principle modifies the Feynman propagator to the form

$`G_F^P(x,x^{})`$ $`=`$ $`{\displaystyle \frac{1}{(4\pi i)^2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{ds}{s^2}}\mathrm{exp}(im^2s)\mathrm{exp}[i((xx^{})^2+L_P^2)/4s]`$ (2)
$`=`$ $`{\displaystyle \frac{m}{4\pi ^2}}{\displaystyle \frac{K_1(im\sqrt{x^2L_P^2iϵ})}{\sqrt{x^2L_P^2iϵ}}},`$
where $`K_1(z)`$ is the modified Bessel function of order 1. \[The metric signature we follow is $`(+,,,)`$.\] In momentum space, the modified propagator is

$$G_F^P(p)=i_0^{\mathrm{}}𝑑z\mathrm{exp}(iL_P^2/4z+i(p^2m^2+iϵ)z).$$
(3)
The presence of a fundamental length scale is a feature that is expected to arise in a quantum theory of gravity. Hence, the modification of the weightage factor as mentioned above can be interpreted as being equivalent to introducing quantum gravitational corrections into standard field theory. The ultraviolet divergences in quantum field theory arise from the singularities of the propagator functions on the light cone, and a smearing out of the light cone due to the quantum gravitational corrections, using the principle of “path integral duality”, will lead to the suppression of these divergences.

In an earlier work ?, the implications of the modified propagator to certain conventional non-perturbative quantum field theoretic results were discussed in detail. It was found that the essential feature of this prescription of path integral duality is to provide an ultra-violet cutoff at the Planck energy scales, thereby obtaining a Lorentz invariant finite results. Encouraged by this fact, we would like to estimate the renormalization factors in QED, and other radiative correction terms which cannot be obtained in the conventional QED calculations.

The standard definition of Feynman propagator for Dirac particles is

$$S(x)=(i\gamma ^\mu _\mu +m)G_F(x),$$
(4)
where, $`G_F(x)`$ is the usual Feynman propagator for the scalar particles. In our analysis of evaluating quantum gravitational corrections to the standard QED calculations, we assume that the effect of “principle of path integral duality” on the Dirac propagator is defined as

$$S^P(x)=(i\gamma ^\mu _\mu +m)G_F^P(x).$$
(5)
Hence, the effect of summing over the quantum fluctuations of the space-time structure, in the low energy scales, can be realized in any field theoretic calculations where the propagator of the scalar particles appear explicitly.

There is a similar approach in the literature by Ohanian ?, who had used a smeared propagator which is Poincare invariant, to calculate the radiative corrections in QED. The modified Feynman propagator in our case is obtained by rigorous evaluation and the duality principle introduces the fundamental length scale in a Lorentz invariant manner. We also comment on the differences between our approaches and results at the appropriate sections below.

Before proceeding to the technical aspects, it is necessary to outline certain conceptual issues related to this approach. In the conventional approaches to quantizing a field based on some classical Lagrangian, one will invariably obtain quantum corrections to the classical Lagrangian. When these corrections are local in the configuration space and contain terms which are of the same form as those in the classical Lagrangian, it is necessary to absorb them into the parameters of the original Lagrangian. This process of renormalization, has — a priori — nothing to do with divergences. Usually, however, the quantum corrections to the theory — calculated by perturbative methods — lead to divergent expressions. If the new (divergent) terms are not of the form of the original terms, then the theory cannot be interpreted perturbatively. However, if the divergent terms have the same structure as the terms in the original Lagrangian, it is possible to use the procedure of renormalization (which, a priori, has nothing to do with divergences) to give meaning to the theory. To do this properly, it is necessary to first evolve a procedure (called regularization) which allows the divergent expressions to be recast as the limit of some finite quantities. After performing the renormalization subtractions, one is free to take the required limit, leading to finite corrections.

The above approach gets modified in two essential aspects when the modified propagator is used. Firstly, when quantum gravitational corrections lead to a cutoff, the quantum corrections will not have any divergent terms; that is, the regularization is now built into the theory with Planck length acting as regulator. But renormalization of the theory is still needed and the physical and bare coupling constants will differ by a finite amount. Secondly, the regularization procedure is now fixed by our ansatz. This is important because, it is well known from standard work in quantum field theory that different regularization procedures are not equivalent. For example, dimensional regularization and momentum space techniques are not completely equivalent as regards their treatment of the symmetries of the system. Since we have no freedom in choosing a regularization, it is necessary to accept and investigate the final results arising from the ansatz. The use of modified propagator to the field theoretic calculation is a more realistic approach towards the removal of the divergences in field theory based on the relevant physics, rather than a formalistic approach based on improper mathematical manipulations. In an earlier work ?, it was shown that the physical parameters in this system like mass and coupling constant have additive finite terms, that are proportional to the powers of $`L_P`$ by calculating the effective potential for a self-interacting scalar field theory using the modified propagator.

In this paper, we evaluate the second order radiative corrections in QED using the modified propagator. The three corrections are to the vacuum polarization, electron self-energy, and vertex function. In the conventional QED calculations the divergent terms are absorbed into the physical parameters like mass, charge and spin. The modified propagator here again acts as a regulator to the ultra-violet divergences in the theory. The modified propagator introduces two kinds of regulators in the radiative correction calculations, which are logarithmic and power law. In the three radiative corrections the leading order power law regulator is $`𝒪(L_P^2m^2)`$ and this is very small when compared to the other non-divergent terms in the conventional QED results. The three renormalization factors has the logarithmic corrections which is of the $`𝒪(\mathrm{ln}(L_Pm))`$.

The following point needs to be stressed regarding the actual values of the corrections. In pure QED, the perturbative expansion is in a series in $`\alpha `$ and corrections are of the order of $`\alpha 10^2`$, $`\alpha ^210^4`$ etc in successive orders. The lowest order quantum gravitational corrections are, by and large, of order $`L_P^2m^210^{45}`$. Since $`\alpha ^{22}(L_Pm)^2`$, the first 22 order corrections of QED will dominate over the quantum gravitational corrections computed here! Much before this, electroweak corrections will start modifying the results. Thus the finite corrections computed here are only of conceptual significance — providing the lowest order corrections from quantum gravity — rather than of any operational significance. The only exception to this general situation is when quantum gravitational effects break a symmetry originally present in the theory, which — as we shall see — does happen.

In sections II and III, we evaluate the second order radiative corrections in QED, and we discuss the results and present other applications in Sec. IV.

## 2 Vacuum Polarization and Electron Self Energy

### 2.1 Vacuum Polarization

The interaction of the photon field with electron field modifies the free photon propagator. The photon propagator of momentum $`q`$, with the one loop radiation correction included, is given by

$$iD_{\mu \nu }^F(q)=iD_{\mu \nu }^F+iD_{\mu \rho }^F\frac{i\mathrm{\Pi }^{\rho \sigma }(q)}{4\pi }iD_{\sigma \nu }^F,$$

where $`\mathrm{\Pi }^{\mu \nu }`$ is the vacuum polarization tensor, and $`D_{\mu \nu }^F`$ is the free photon propagator. The free photon propagator, $`D_{\mu \nu }^F`$, in the above relation is the conventional QED propagator. We concentrate here on the corrections to $`D_{\mu \nu }`$ arising from the modification of $`\mathrm{\Pi }_{\mu \nu }`$ rather than from direct modification of $`D_{\mu \nu }`$ due to our ansatz. Using the Feynman rules of QED in the momentum space, we can write the vacuum polarization tensor as

$`\mathrm{\Pi }_{\mu \nu }(q)`$ $`=`$ $`16\pi ie^2{\displaystyle \frac{d^4k}{(2\pi )^3}\left(k_\mu (kq)_\nu +(kq)_\mu k_\nu g_{\mu \nu }(k^2qkm^2)\right)}`$ (6)
$`\times `$ $`G_F^P(k)G_F^P(kq).`$
Substituting for the propagator from Eqn.(3), and following similar calculations as in conventional QED, the vacuum polarization tensor can be separated into gauge invariant and gauge non-invariant part ?. The resultant gauge invariant part is given by

$$\mathrm{\Pi }_{\mu \nu }^1(q^2)=\frac{4e^2}{\pi }\left(q_\nu q_\mu g_{\mu \nu }q^2\right)_0^1𝑑zz(1z)K_0(\xi ),$$
(7)
where

$$\xi ^2=L_P^2\frac{m^2q^2z(1z)}{z(1z)}.$$

The series expansion of $`K_0(z)`$, about the origin, is given by

$$K_0(z)=\gamma \mathrm{ln}(z/2)\frac{z^2}{4}\left[1\gamma \mathrm{ln}(z/2)\right]+\mathrm{}.$$
(8)
Hence, Eqn.(7) takes the form

$`\mathrm{\Pi }_{\mu \nu }^1(q^2)`$ $`=`$ $`(q_\nu q_\mu g_{\mu \nu }q^2)[(Z_31){\displaystyle \frac{2e^2}{\pi }}A_1{\displaystyle \frac{e^2L_P^2}{\pi }}({\displaystyle \frac{q^2}{2}}A_1`$ (9)
$`+`$ $`{\displaystyle \frac{1}{2}}(1\gamma \mathrm{ln}(L_Pm/2))(m^2{\displaystyle \frac{q^2}{6}}))],`$
where

$$A_1=_0^1𝑑zz(1z)\mathrm{ln}\left(\frac{m^2}{m^2q^2z(1z)}\right)$$
(10)
to the lowest order of $`K_0(\xi )`$. The term $`A_1`$ is the familiar conventional QED non-divergent term and the remaining terms (of the order $`𝒪(L_P^2)`$) are the leading order quantum gravitational power law corrections/regulators to the conventional QED terms. The contribution of the quantum gravitational corrections to the vacuum polarization is extremely small as compared to the conventional non-divergent QED term, i.e. they are of the order $`10^{45}`$ which is much smaller compared to the next order radiative corrections of QED. But as we can notice the quantum gravitational corrections to the conventional QED non-divergent terms become important when $`L_Pq1`$. In this high momentum transfer limit, the propagation effects of the virtual photons will probe the small scale quantum gravitational effects and the corrections will be of the same order as the conventional QED non-divergent terms. The charge renormalization factor, $`Z_3`$, which is divergent in the usual QED calculations is now finite and is given by

$$Z_31=\frac{2e^2\mathrm{ln}(L_Pm/2)}{\pi }\left[\frac{1}{3}+\frac{6\gamma 5}{18\mathrm{ln}(L_Pm/2)}+𝒪(L_P^2m^2)\right].$$
(11)
The estimate of this factor comes out to be

$$Z_31=\frac{2e^2}{3\pi }\mathrm{ln}(L_Pm)0.1.$$
(12)
In the earlier work ?, authors have calculated charge renormalization factor using Effective action approach, which is a non-perturbative technique. In this approach, they calculated the effect of the classical electro-magnetic background on the quantum charged scalar fields, propagating in the flat space-time, using the above modified propagator. The charge renormalization factor they obtain using the effective action approach is of the form $`Z^P=\frac{q^2}{6\pi }K_0(2L_Pm).`$ On expanding $`K_0(z)`$ using Eqn.(8), their estimate of the charge renormalization factor was also of the order of $`0.1`$.

In the standard QED calculations, there does arise a gauge non-invariant part of the vacuum polarization tensor which is divergent and is renormalized to zero(for instance, see Hatfield ?). Even though, this is a standard result and can be found in textbooks we have given the relevant steps in the Appendix for the sake of completeness. By retracing the steps given in Appendix, we now obtain a finite quantity to the gauge non-invariant part of the vacuum polarization tensor using the modified propagator. The gauge non-invariant part comes out to be

$`\mathrm{\Pi }_{\mu \nu }^2(q)`$ $`=`$ $`{\displaystyle \frac{iL_P^2e_0^2}{(4\pi )^2}}g_{\mu \nu }{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{z}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz^{}}{z^{}}}{\displaystyle \frac{1}{(z+z^{})^2}}\mathrm{exp}\left({\displaystyle \frac{iL_P^2}{4}}(z^1+z^1)\right)`$ (13)
$`\times `$ $`\mathrm{exp}i\left({\displaystyle \frac{q^2zz^{}}{z+z^{}}}m_0(z+z^{})\right).`$
In Eqn.(13), transforming the variables $`z`$ and $`z^{}`$ to a new set of variables by the relation $`z=1/t`$ and $`z^{}=1/t^{}`$ leads to

$`\mathrm{\Pi }_{\mu \nu }^2(q)`$ $`=`$ $`{\displaystyle \frac{iL_P^2e_0^2}{(4\pi )^2}}g_{\mu \nu }{\displaystyle _0^{\mathrm{}}}t𝑑t{\displaystyle _0^{\mathrm{}}}t^{}𝑑t^{}{\displaystyle \frac{1}{(t+t^{})^2}}\mathrm{exp}\left({\displaystyle \frac{iL_P^2}{4}}(t+t^{})\right)`$ (14)
$`\times `$ $`\mathrm{exp}i\left({\displaystyle \frac{q^2}{t+t^{}}}m_0(t^1+t^1)\right).`$
Using the standard identity,

$$1=_0^{\mathrm{}}\frac{d\beta }{\beta }\delta \left(1\beta (z+z^{})\right)$$

and scaling $`t_it_i/\beta `$, as in the conventional QED calculations, we get

$$\mathrm{\Pi }_{\mu \nu }^2(q)=\frac{2ie_0^2}{(4\pi )^2}g_{\mu \nu }L_P^6_0^1t(1t)\frac{K_2(\xi )}{\xi ^2}𝑑t,$$
(15)
where

$$\xi ^2=L_P^2(\frac{m_0^2}{t(1t)}q^2)$$

The expansion of $`K_2(z)`$ near the origin is given by

$$K_2(z)=\frac{2}{z^2}\frac{1}{2}+\left(\frac{34\gamma }{32}+\frac{\mathrm{ln}(2)\mathrm{ln}(z)}{8}\right)z^2+\mathrm{}.$$

Substituting the series expansion of $`K_2(\xi )`$ in Eqn.(15), we obtain

$`\mathrm{\Pi }_{\mu \nu }^2(q)`$ $`=`$ $`g_{\mu \nu }{\displaystyle \frac{2ie_0^2}{(4\pi )^2}}[2L_P^2{\displaystyle _0^1}dt{\displaystyle \frac{(tt^2)^3}{(m_0^2q^2t(1t))^2}}`$ (16)
$``$ $`{\displaystyle \frac{L_P^4}{2}}{\displaystyle _0^1}dt{\displaystyle \frac{(tt^2)^2}{(m_0^2q^2t(1t))}}+\mathrm{}].`$
This clearly shows that the the gauge non-invariant part is nonzero for $`L_P0`$ and vanishes as $`L_P0`$. While the term which breaks the gauge invariance is small to be of operational significance, it does have certain conceptual importance. The following points need to be noted regarding this result: (i) Mathematically speaking, this result arises from the fact that our ansatz is equivalent to a momentum space regularization procedure in conventional QED. It is known that, momentum space regularization, in contrast to dimensional regularization, can lead to gauge breaking terms. Usually, this is considered as an argument in favor of dimensional regularization. In our approach, of course, we have no choice and the result arises automatically. In fact, it is very likely that any quantum gravitational cutoff will appear like a momentum space regulator and will break the gauge symmetry. (ii) Previously, Ohanian ? had obtained a gauge breaking term using a gravitationally smeared propagator. There is however one vital difference between our result and the one obtained by him. Note that, in our approach, the propagator reduces to that of conventional field theory when $`L_P0`$. If the procedure is to be consistent, the gauge breaking term should vanish when the limit of $`L_P0`$ is taken. This is true as regards our result in Eqn.(16) showing that this is indeed a quantum gravitational effect. However, Ohanian ? obtains a gauge breaking term which does not vanish in the corresponding limit. This suggests that, our approach does allow a consistent interpretation of the results. (iii) It is the extra term $`L_P^2(z+z^{})^2/(4zz^{})`$ which breaks the gauge invariance of the electro-magnetic field. This extra factor can be associated to the current in the charge conservation relation and hence, implying that the charge conservation is no more valid. There have been – more drastic! – attempts in the literature to break even the Lorentz invariance by introducing a coupling of the photon to charged scalar field, through gravitational couplings to the photon, etc ?. The basic aim in the process of breaking the conformal symmetry of the electro-magnetic field is allowing for the possibility of generating large scale magnetic fields within inflationary scenarios. Most of these studies have used ad-hoc interaction potential to break the conformal invariance. The breaking of gauge invariance of the electro-magnetic field in our case is from a much more deeper “principle of path integral duality”. The connection to the cosmological seed magnetic field is still under investigation.

### 2.2 Electron Self Energy

The interaction of the electron field with photon field modifies the free electron propagator. The electron propagator of momentum $`p`$, with the one-loop correction included, is given by

$$iS_F^{}(p)=iS_F(p)+iS_F(p)(i\mathrm{\Sigma }(p))iS_F(p),$$

where $`\mathrm{\Sigma }(p)`$ (a 4-spinor) is the self-energy function, and $`S_F(p)`$ is the free electron propagator. Using the Feynman rules in the momentum space, we get

$$\mathrm{\Sigma }(p)=4\pi ie^2\frac{d^4k}{(2\pi )^4}\gamma ^\mu (\gamma ^\nu p_\nu \gamma ^\nu k_\nu +m)\gamma _\mu G_F^P(k)G_F^P(pk).$$
(17)
Substituting the propagator from the Eqn.(3), the above equation reduces to

$$\mathrm{\Sigma }(p)=\frac{e^2}{4\pi ^2}_0^1𝑑z\left(2m_0\gamma ^\mu p_\mu z\right)K_0(\xi ),$$
(18)
where

$$\xi ^2=\frac{L_p^2}{z(1z)}\left[m_0^2zp^2z(1z)\right].$$

On expanding $`K_0(\xi )`$ using Eqn.(8) we obtain, the lower order terms corresponding to the conventional QED results and the higher order terms as contributions of the quantum fluctuations of the space-time to the self energy of electron. The electron wave function renormalization, and the shift in the mass are obtained by recasting the interacting field propagator to look like the free field propagator i.e. we set

$$Z_2(\gamma ^\mu p_\mu m_0\mathrm{\Sigma }(p))=\gamma ^\mu p_\mu m+\mathrm{finite}\mathrm{terms}.$$

The electron wave function renormalization factor is obtained by equating terms proportional to $`p`$ in the above equation. This results in

$$Z_2^11=\frac{e^2}{8\pi ^2}\left[\gamma +\mathrm{ln}2\mathrm{ln}(L_pm_0)+𝒪(1)\right].$$
(19)
\[Note that we have neglected the higher order terms of $`K_0(\xi )`$ as these have the dependence as $`L_P^2m_0^2`$ whose contributions are negligible\]. The last term in the above expression is a finite quantity and is of the order one and it has the dependence of the electron momentum. The shift in the mass is given by

$$\frac{\delta m}{m_0}=\frac{e^2}{8\pi ^2}\left[3\mathrm{ln}(L_Pm_0/2)+3\gamma +𝒪(1)\right]$$
(20)
The estimate of the scale factor $`Z_2`$ comes out to be

$$Z_2^11=\frac{\alpha }{\pi }\mathrm{ln}m_0L_P0.1,$$
(21)
and the fractional shift in the mass is

$$\frac{\delta m}{m}\frac{\alpha }{\pi }\mathrm{ln}m_0L_P0.1.$$
(22)
The usual infrared divergences is ignored in the calculation of the renormalization factor $`Z_2`$. Here again, we see that both the mass renormalization factor and the mass shift are also of the same order as the charge renormalization factor i.e. $`0.1`$.

## 3 Vertex Correction and Anomalous magnetic moment

### 3.1 Vertex Correction

For a free Dirac field, the current density is defined as

$$J^\mu =\overline{\psi }\gamma ^0\gamma ^\mu \psi .$$

Thus, in the low energy QED processes the current transfer is related by $`\gamma ^\mu `$, which is the low energy vertex function. The radiative corrections will modify the vertex, to the one-loop correction, as

$$ie\mathrm{\Lambda }_\mu =ie\gamma _\mu ie\mathrm{\Gamma }_\mu ,$$
(23)
where, $`\mathrm{\Gamma }_\mu `$ is the vertex function. Using the Feynman rules, we obtain

$`\mathrm{\Gamma }_\mu (p^{},p)`$ $`=`$ $`(ie_0)^2{\displaystyle \frac{d^4k}{(2\pi )^4}\left[\gamma _\mu (\gamma _\rho p^\rho \gamma _\rho k^\rho +m_0)\gamma _\mu (\gamma _\sigma p^\sigma \gamma _\sigma k^\sigma +m_0)\right]}`$ (24)
$`\times `$ $`G_F^P(k)G_F^P(p^{}k)G_F^P(pk).`$
Substituting for the propagator from the Eqn. (3), the above equation reduces to

$`\mathrm{\Gamma }_\mu (p^{},p)`$ $`=`$ $`(Z_1^11)\gamma _\mu +{\displaystyle \frac{(ie_0)^2}{(4\pi )^2}}{\displaystyle _0^1}𝑑z{\displaystyle _0^{1z}}𝑑z^{}{\displaystyle \frac{\xi }{T_1}}`$
$`\times `$ $`\left[2\gamma _\mu (p^{}p)^2(1z^{})(1z)4imz(1zz^{})(p^{}p)^\mu \sigma _{\mu \nu }\right],`$
where

$$\xi ^2=L_P^2T_1T_2,T_1=(zp^{}+z^{}p)^2zp^2z^{}p^2+(z+z^{})m_0^2,$$

and

$$T_2=(1zz^{})^1+z^1+z^1.$$

The vertex renormalization factor $`Z_1`$ is given by

$`Z_1^11`$ $`=`$ $`i{\displaystyle \frac{(ie_0)^2}{(4\pi )^2}}{\displaystyle _0^1}dz{\displaystyle _0^{1z}}dz^{}[4K_0(\xi )+2m_0^2{\displaystyle \frac{\xi }{T_1}}K_1(\xi )`$ (26)
$`\times `$ $`(2+2(z+z^{})+(z+z^{})^2)].`$
Substituting for $`K_0(\xi )`$ using Eqn.(8), and using the series expansion of $`K_1(\xi )`$ as

$$K_1(z)=\frac{1}{z}+\frac{z}{2}\left(\mathrm{ln}(z/2)+\frac{2\gamma 1}{2}\right)+\frac{z^3}{16}\left(\mathrm{ln}(z/2)\frac{54\gamma }{4}\right),$$
(27)
in the Eqn. (26), the vertex normalization factor comes out to be

$`Z_1^11`$ $`=`$ $`i{\displaystyle \frac{(ie_0)^2}{(4\pi )^2}}{\displaystyle _0^1}dz{\displaystyle _0^{1z}}dz^{}[4\mathrm{ln}(\xi /2)4\gamma +{\displaystyle \frac{2m_0^2}{T_1}}`$ (28)
$`\times `$ $`(2+2(z+z^{})+(z+z^{})^2)]`$
$``$ $`{\displaystyle \frac{\alpha }{\pi }}\mathrm{ln}(m_0L_P),`$ (29)
which is roughly of the order of $`0.1`$. The estimated renormalization factors $`Z_1`$ and $`Z_3`$ using the modified propagator are not equal unlike in conventional QED.

### 3.2 Anomalous Magnetic Moment

The triumph of QED has been the precision test of electron anomalous magnetic moment. The experimental value of the anomalous magnetic moment of an electron is in excellent agreement with the predicted perturbative calculations up to the $`4^{th}`$ order to the $`15^{th}`$ decimal place. It is therefore of interest to compute the quantum gravitational corrections to the magnetic moment of an electron using the modified propagator.

The vertex correction contribution to scattering of an electron in an external field is given by $`\overline{u}(p^{})\mathrm{\Gamma }^\mu (p^{},p)u(p)A_\mu ^c(p^{}p)`$, where $`\overline{u}`$ and $`u`$ are the spinor wave-functions. Since, our interest is in calculating the radiative and the quantum gravitational corrections to the gyro-magnetic ratio of an electron, the term involving $`\sigma _{\mu \nu }`$ in the Eqn. (3.1) is of our concern. The corresponding $``$ matrix for this process is given by ?,

$``$ $`=`$ $`\overline{u}(p^{})\gamma _\mu ^Mu(p)`$
$`=`$ $`\overline{u}(p^{})4m_0{\displaystyle \frac{(ie_0)^2}{(4\pi )^2}}{\displaystyle _0^1}𝑑z{\displaystyle _0^{1z}}𝑑z^{}z(1zz^{})(p^{}p)^\nu {\displaystyle \frac{\sigma _{\mu \nu }\xi }{T_1}}K_1(\xi )u(p).`$
In the limit of small momentum transfer, $`(p^{}p)^2m_0^2`$, one obtains by expanding $`K_1(\xi )`$ from the Eqn. (27), we get

$``$ $`=`$ $`\overline{u}(p^{}){\displaystyle \frac{\alpha }{m\pi }}(p^{}p)^\mu \sigma _{\mu \nu }{\displaystyle _0^1}𝑑z{\displaystyle _0^{1z}}𝑑z^{}{\displaystyle \frac{z(1zz^{})}{(z+z^{})^2}}u(p)`$ (31)
$``$ $`\overline{u}(p^{}){\displaystyle \frac{\alpha }{m\pi }}{\displaystyle \frac{L_P^2m_0^2}{24}}(p^{}p)^\mu \sigma _{\mu \nu }{\displaystyle _0^1}𝑑z{\displaystyle _0^{1z}}𝑑z^{}z(1zz^{})T_2`$
$`\times `$ $`\mathrm{ln}(L_Pm_0^2(z+z^{})^2T_2/4)u(p).`$
The first term in the above equation corresponds to the usual QED radiative vertex correction of the gyro-magnetic ratio (g) to the order $`e_0^2`$. The second term in the above equation is the quantum gravitational corrections to the gyro-magnetic ratio of an electron. The integral in the second term of the Eqn. (31) is convergent. The contribution of the quantum gravitational correction to the gyro-magnetic ratio to the first order in $`\alpha `$ is of the order $`L_P^2m^210^{45}`$. Obviously, this is not of practical significance.

## 4 Conclusions and Discussion

In this paper, we evaluated the quantum gravitational corrections(QGC) to three radiative corrections, in the first order of $`\alpha `$, in QED using the “principle of path integral duality”. The modified propagator is able to remove all the divergences which usually crop up in the conventional QED calculations. The main features of the modified propagator in QED are as follows: (a) The three renormalization factors($`Z_1`$, $`Z_2`$, and $`Z_3`$), and the mass shift are all of the same order, $`𝒪(\mathrm{ln}(mL_P))0.1`$. The charge renormalization factor $`Z_3`$ using the non-perturbative methods in scalar QED is also of order $`0.1`$ ?. The renormalization factors $`Z_1`$ and $`Z_3`$ are different in our case as opposed to the conventional QED calculations. (b) The modified propagator makes the gauge non-invariant part of the vacuum polarization tensor to be non-zero for $`L_P0`$ vanishes and in the limit $`L_P0`$. This breaking of the gauge symmetry is also related to also the difference between the two renormalization factors. We have briefly indicated the possible effect of this breaking of gauge invariance, to the generation of large scale magnetic fields within inflationary scenarios (c) The contribution of the QGC to the gyro-magnetic ratio is very small i.e. of the order of $`L_P^2m^2`$ $`10^{45}`$. The contribution of the quantum gravitational correction to the vacuum polarization and the electron self energy is also of the same order.

Acknowledgments

S.S. is being supported by the Council of Scientific and Industrial Research, India.

Appendix A: Evaluation of Gauge Non-invariant part

For the sake of completeness, we outline the essential steps leading to the gauge non-invariant part of the vacuum polarization tensor in standard QED and how it is regularized to zero. The vacuum polarization tensor in the conventional QED calculations in $`2n`$ dimensions is given by

$`\mathrm{\Pi }_{reg}^{\mu \nu }(q)`$ $`=`$ $`(e\mu ^{2n})^2{\displaystyle \frac{d^{2n}k}{(2\pi )^{2n}}\frac{i}{k^2m^2+iϵ}\frac{i}{(kq)^2m^2+iϵ}}`$
$`\times `$ $`2^n\left(k^\mu (kq)^\nu +k^\nu (kq)^\mu ((k^2m^2)kq)g^{\mu \nu }\right),`$
where $`\mu `$ has the dimension of mass. Using the integral representation of the propagator, i.e.

$$\frac{i}{k^2m^2+iϵ}=_0^{\mathrm{}}𝑑z\mathrm{exp}\left(iz\left(k^2m^2+iϵ\right)\right)$$
(A.2)
in Eqn.(4) and completing the squares in the exponential, we obtain

$`\mathrm{\Pi }_{reg}^{\mu \nu }(q)`$ $`=`$ $`(e\mu ^{2n})^2{\displaystyle _0^{\mathrm{}}}𝑑z_1𝑑z_2{\displaystyle \frac{d^{2n}k}{(2\pi )^{2n}}}`$
$`\times `$ $`2^n\left(k^\mu (kq)^\nu +k^\nu (kq)^\mu ((k^2m^2)kq)g^{\mu \nu }\right)`$
$`\times `$ $`\mathrm{exp}\left(i(z_1+z_2)\left(k{\displaystyle \frac{z_2q}{z_1+z_2}}\right)^2+i{\displaystyle \frac{z_1z_2q^2}{z_1+z_2}}i(m^2iϵ)(z_1+z_2)\right).`$
Shifting the variable of integration and using the relations,

$$\frac{d^{2n}p}{(2\pi )^{2n}}p^2\mathrm{exp}(iap^2)=\frac{n}{(4\pi a)^n}\frac{1}{a}\mathrm{exp}(in\pi /2),$$

$$\frac{d^{2n}p}{(2\pi )^{2n}}p^\mu p^\nu \mathrm{exp}(iap^2)=\frac{g^{\mu \nu }}{(4\pi a)^n}\frac{1}{2a}\mathrm{exp}(in\pi /2),$$

we get

$$\frac{\mathrm{\Pi }_{reg}^{\mu \nu }(q)}{(e\mu ^{2n})^2}=\frac{\mathrm{exp}(in\pi /2)}{(4\pi )^n}\left[g_{\mu \nu }Q(q^2,m^2)+(q^\mu q^\nu g^{\mu \nu }q^2)P(q^2,m^2)\right],$$
(A.4)
where

$`Q(q^2,m^2)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑z_1{\displaystyle _0^{\mathrm{}}}𝑑z_2{\displaystyle \frac{2^n}{(z_1+z_2)^n}}\left({\displaystyle \frac{n1}{(z_1+z_2)}}+im^2iq^2{\displaystyle \frac{z_1z_2}{(z_1+z_2)^2}}\right)`$ (A.5)
$`\times `$ $`\mathrm{exp}\left(iq^2{\displaystyle \frac{z_1z_2}{z_1+z_2}}i(m^2iϵ)(z_1+z_2)\right),`$
and

$`P(q^2,m^2)`$ $`=`$ $`i2^{n+1}{\displaystyle _0^{\mathrm{}}}𝑑z_1{\displaystyle _0^{\mathrm{}}}𝑑z_2{\displaystyle \frac{z_1z_2}{(z_1+z_2)^{n+2}}}`$ (A.6)
$`\times `$ $`\mathrm{exp}\left(iq^2{\displaystyle \frac{z_1z_2}{z_1+z_2}}i(m^2iϵ)(z_1+z_2)\right).`$
Global gauge invariance of the action leads to current conservation, $`_\mu j^\mu =0`$. This implies $`q_\mu \mathrm{\Pi }^{\mu \nu }(q)=0`$. The factor proportional to $`g^{\mu \nu }`$ is the gauge non-invariant part of the vacuum polarization tensor and does not seem to satisfy the gauge invariance condition.

We summarize here the standard argument \[see Hatfield ?\] to show that the first term, $`Q(q^2,m^2)`$, in Eqn.(A.4) can be shown to vanish and hence the regularization preserves the gauge symmetry. Consider only the first term of the gauge non-invariant part and define

$$I(q^2,m^2)_0^{\mathrm{}}𝑑z_1𝑑z_2\frac{1}{(z_1+z_2)^{n+1}}\mathrm{exp}\left(iq^2\frac{z_1z_2}{z_1+z_2}im^2(z_1+z_2)\right).$$
(A.7)
Rescaling the integration variables, $`z_1`$ and $`z_2`$, in the above equation by $`\beta `$, i.e. $`z_1\beta z_1`$ and $`z_2\beta z_2`$, we obtain

$`I(\beta ,q^2,m^2)`$ $`=`$ $`{\displaystyle \frac{1}{\beta ^{n1}}}{\displaystyle _0^{\mathrm{}}}𝑑z_1𝑑z_2{\displaystyle \frac{1}{(z_1+z_2)^{n+1}}}`$ (A.8)
$`\times `$ $`\mathrm{exp}\left(i\beta \left(q^2{\displaystyle \frac{z_1z_2}{z_1+z_2}}m^2(z_1+z_2)\right)\right).`$
The quantity $`I(\beta ,q^2,m^2)`$ is of course independent of the integration variables $`z_1`$ and $`z_2`$, and hence is also independent of the parameter($`\beta `$), i.e. $`I(\beta ,q^2,m^2)/\beta `$ $`=0`$. Differentiating the above expression w.r.t $`\beta `$ and regrouping terms, we get

$`0=\beta {\displaystyle \frac{I}{\beta }}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑z_1𝑑z_2{\displaystyle \frac{\beta ^{2n}}{(z_1+z_2)^n}}\left({\displaystyle \frac{n1}{(z_1+z_2)}}+im^2iq^2{\displaystyle \frac{z_1z_2}{(z_1+z_2)^2}}\right)`$ (A.9)
$`\times `$ $`\mathrm{exp}\left(i\beta \left(q^2{\displaystyle \frac{z_1z_2}{z_1+z_2}}m^2(z_1+z_2)\right)\right).`$
Rescaling the variables $`z_1`$ and $`z_2`$ by $`\beta ^1`$, i.e. $`z_iz_i/\beta `$, in the above expression, we get

$`0=\beta {\displaystyle \frac{I}{\beta }}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑z_1𝑑z_2{\displaystyle \frac{2^n}{(z_1+z_2)^n}}\left({\displaystyle \frac{n1}{(z_1+z_2)}}+im^2iq^2{\displaystyle \frac{z_1z_2}{(z_1+z_2)^2}}\right)`$ (A.10)
$`\times `$ $`\mathrm{exp}\left(i\left(q^2{\displaystyle \frac{z_1z_2}{z_1+z_2}}m^2(z_1+z_2)\right)\right)=Q(q^2,m^2).`$
Thus, the gauge non-invariant part of the vacuum polarization tensor in the standard QED vanishes. Hence, $`q_\mu \mathrm{\Pi }^{\mu \nu }(q)=0`$.

In the rest of this appendix, we retrace the above steps for the modified propagator. The gauge non-invariant part of the vacuum polarization tensor gets modified to the form,

$$\mathrm{\Pi }_{\mu \nu }^{(2)}(q^2,L_P)=(e\mu ^{2n})^2\frac{\mathrm{exp}(in\pi /2)}{(4\pi )^n}g_{\mu \nu }\mathrm{\Pi }_0^{(2)}(q^2,L_P),$$
(A.11)
where

$`\mathrm{\Pi }_0^{(2)}(q^2,L_P)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑z_1𝑑z_2{\displaystyle \frac{2^n}{(z_1+z_2)^n}}\left[{\displaystyle \frac{n1}{(z_1+z_2)}}+im^2iq^2{\displaystyle \frac{z_1z_2}{(z_1+z_2)^2}}\right]`$
$`\times `$ $`\mathrm{exp}\left[iq^2{\displaystyle \frac{z_1z_2}{z_1+z_2}}i(m^2iϵ)(z_1+z_2)+{\displaystyle \frac{iL_P^2}{4}}(z_1^1+z_2^1)\right],`$
when the modified propagator is used. We now show that the gauge non-invariant part is a finite quantity and is of the order $`L_P^2`$. Here again, we consider the first term of the gauge non-invariant part and define

$`I(q^2,m^2,L_P)`$ $``$ $`{\displaystyle _0^{\mathrm{}}}𝑑z_1𝑑z_2{\displaystyle \frac{1}{(z_1+z_2)^{n+1}}}\mathrm{exp}\left(iq^2{\displaystyle \frac{z_1z_2}{z_1+z_2}}im^2(z_1+z_2)\right)`$ (A.13)
$`\times `$ $`\mathrm{exp}\left({\displaystyle \frac{iL_P^2}{4}}(z_1^1+z_2^1)\right).`$
Rescaling the variables $`z_1`$ and $`z_2`$ by $`\beta `$, i.e $`z_i\beta z_i`$, and differentiating the resultant of the above expression w.r.t $`\beta `$ and regrouping terms, we obtain

$$\beta \frac{I(\beta ,q^2,m^2,L_P)}{\beta }Q(\beta ,q^2,m^2,L_P)R(\beta ,q^2,m^2,L_P),$$
(A.14)
where

$`Q(\beta ,q^2,m^2,L_P)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑z_1𝑑z_2{\displaystyle \frac{\beta ^{2n}}{(z_1+z_2)^n}}`$
$`\times `$ $`\left({\displaystyle \frac{n1}{(z_1+z_2)}}+im^2iq^2{\displaystyle \frac{z_1z_2}{(z_1+z_2)^2}}\right)`$
$`\times `$ $`\mathrm{exp}\left[i\beta \left(q^2{\displaystyle \frac{z_1z_2}{z_1+z_2}}m^2(z_1+z_2)\right)+{\displaystyle \frac{iL_P^2}{4\beta }}(z_1^1+z_2^1)\right],`$
and

$`R(\beta ,q^2,m^2,L_P)`$ $``$ $`{\displaystyle _0^{\mathrm{}}}𝑑z_1𝑑z_2{\displaystyle \frac{iL_P^2}{4z_1z_2}}{\displaystyle \frac{\beta ^n}{(z_1+z_2)^n}}`$
$`\times `$ $`\mathrm{exp}\left(i\beta \left(q^2{\displaystyle \frac{z_1z_2}{z_1+z_2}}m^2(z_1+z_2)\right)+{\displaystyle \frac{iL_P^2}{4\beta }}(z_1^1+z_2^1)\right).`$
The term $`I(\beta ,q^2,m^2,L_P)`$ is independent of $`\beta `$ (by the arguments stated earlier in this appendix) and hence, the partial differential $`I/\beta `$ vanishes. Hence,

$$Q(\beta ,q^2,m^2,L_P)=R(\beta ,q^2,m^2,L_P).$$

Rescaling the variables $`z_1`$ and $`z_2`$ by $`\beta ^1`$, i.e. $`z_iz_i/\beta `$, in Eqn.(A.14), we get

$$Q_{rescaled}(q^2,m^2,L_P)=R_{rescaled}(q^2,m^2,L_P)=\mathrm{\Pi }_0^{(2)}(q,L_P),$$
(A.17)
where

$`R_{rescaled}`$ $`=`$ $`{\displaystyle \frac{iL_P^2}{4}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz_1}{z_1}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz_2}{z_2}}{\displaystyle \frac{1}{(z_1+z_2)^n}}\mathrm{exp}\left({\displaystyle \frac{iL_P^2}{4}}(z_1^1+z_2^1)\right)`$ (A.18)
$`\mathrm{exp}i\left({\displaystyle \frac{q^2z_1z_2}{z_1+z_2}}m_0(z_1+z_2)\right),`$
and $`\mathrm{\Pi }_0^{(2)}(q,L_P)`$ is the quantity proportional to the gauge non-invariant part of the vacuum polarization tensor defined in Eqn.(A.11). Using our ansatz, we have shown that the gauge non-invariant part of the vacuum polarization tensor is a finite quantity and vanishes as $`L_P0`$.

References


format: text

text:
# Effect of boundary conditions on diffusion in two-dimensional granular gases

## I Introduction

The interactions among grains and between grains and the boundaries influence profoundly the macroscopic behavior of granular systems. To study such complex many body systems, numerical simulations are frequently used where one of the most important ingredients are the collision laws introduced to treat interactions . For dilute assemblies of grains one can use molecular dynamics algorithms where periodic boundary conditions are usually used. If the system is initially in a square box, a particle going out on the left re-enters the system on the right. We will show below that this kind of boundary condition modifies the general dynamics of the grains and introduces large correlations in time. This changes the diffusive behavior of the grains. In this paper, we propose an alternative approach to calculate numerically the coefficient of diffusion, accurately and with only very small finite size effects. To validate our methods, we compare our results for an elastic gas with the Langevin theory.

As granular gases are dissipative it is necessary to feed energy into the system to keep the particles agitated. To thermalize the system, we choose a random acceleration added to each grain at regular time step intervals $`dt`$. Our final goal in this paper is to study the dependence of the dynamic properties of the granular gas on the mode used to force the system. This work is a first step towards understanding the diffusion process in a binary system composed of two grain sizes. The system considered here is composed of one particle, s, of radius $`R_s`$ in a sea of particles of radius $`R_b`$. The particles are spheres constrained to move in a plane and which interact along their equators so that the system is two dimensional. The system considered here is dilute with a packing fraction of 30 %. The simulations are done with the molecular dynamics algorithms (time step driven  and event driven ).

To characterize the diffusive behavior, we focus on the mean square displacement of the s particle. It is well known that for a 2D gas, the integral of the auto-correlation function does not converge . This means that the mean square displacement does not vary linearly with time. Therefore, strictly speaking, we cannot define a diffusion coefficient in 2D. However, we show that in a limited range of time, in the stationary state, the mean square displacement can be approximated by the linear function:

$$<(\stackrel{}{r}(t+t_0)\stackrel{}{r}(t_0))^2>4Dt$$
(1)
where $`D`$ can be interpreted as a diffusion coefficient. All quantities are expressed in arbitrary units.

## II Choice of boundary conditions

In this section we show that periodic boundary conditions introduce strong correlations and therefore alter the diffusion process.

### A Periodic boundary conditions

Consistent with common practice, we have used periodic boundary conditions to simulate a system of identical spheres $`R_s=R_b=0.5`$. Initially the particles are placed randomly in a square box of length $`L`$. The number of particles is calculated for each system depending on $`L`$, $`R_s`$ and the packing fraction. Periodic boundary conditions are applied in both directions. In this case, for elastic or forced gases (section IV), we have observed a strong dependence of $`D`$ (or of the mean square displacement) on the system size.

In Fig. 1, we have plotted the mean square displacement $`<(\stackrel{}{r}(t+t_0)\stackrel{}{r}(t_0))^2>`$ (Fig. 1a) and $`C(t)𝑑t`$ (Fig. 1b), both calculated in the stationary state, as function of $`t`$. $`C(t)`$ is the normalized autocorrelation function:

$$C(t)=\frac{<\stackrel{}{\mathrm{v}}(t_0+t)\stackrel{}{\mathrm{v}}(t_0)><\stackrel{}{\mathrm{v}}(t_0)>^2}{<\stackrel{}{\mathrm{v}}(t_0)^2><\stackrel{}{\mathrm{v}}(t_0)>^2}.$$
(2)
First, we note that the mean square displacement, at large time, varies linearly with time as expected but the slope of the curve, i.e. the diffusion coefficient, increases with system size. We show, in the inset to Fig. 1a, that this dependence on $`L`$ appears already at short time, when $`<(r^2(t+t_0)\stackrel{}{r}(t_0))^2>L^2`$. This feature can be also observed in $`\underset{0}{\overset{t\mathrm{}}{}}C(t)𝑑t`$, which is proportional to the diffusion coefficient. Similarly, we observe that the relaxation time $`\tau _r`$ (i.e. $`C(\tau _r)0`$) increases with size. In summary, the bigger the system is, the longer the characteristic time $`\tau _r`$ and the larger the diffusion coefficient $`D`$ are. We recall that such dependence has been observed by Alder et al . They proposed the following law for the dependence of $`D`$ on the number of particles, $`N`$,

$$D(N)=D(\mathrm{})(12/N).$$
(3)
However their numerical simulations do not support this conjecture  since they fail to observe any saturation of $`D`$ for large systems. In addition, they found strong correlations in the velocity field characterized by the presence of vortex flow pattern at the microscopic scale. Our results confirm the lack of convergence for $`D`$ with system size. In addition, this variation of $`D`$ with $`L`$ is also observed in the case of inelastic collisions.

Another important remark is in order. If the system size is, for example, 60 (with 1400 particles of radius $`R=0.5`$), the characteristic time $`\tau _r`$ is found to be around $`20`$ which represents about 200 collisions for a particle. This means that a particle needs to undergo 200 collisions to lose completely the memory of its past. According to the Boltzmann theory this time should be limited to only a few collisions. Therefore we cannot accept this result as a valid macroscopic description of a gas. It is worth noting that the same results are found for both, the time step driven and the event driven algorithms.

We now discuss some points helpful for understanding the problem. Initially, each particle has a random velocity drawn from a Maxwellian distribution. We shift the linear and angular momenta so that the system has zero center of mass momentum and zero angular momentum relative to the center of mass. We find, however, that, although the system keeps its center of mass at rest throughout the simulation, the system is no longer isotropic, its moment of inertia becoming that of an ellipsoid. Let $`I(t)`$ be the inertia matrix of the system. Its two eigenvalues $`\lambda _n`$ and $`\lambda _p`$ are related via

$$\lambda _n+\lambda _p=m\underset{k=1}{\overset{N}{}}r_k^2(t),$$
(4)
where the sum is over all $`N`$ particles each of mass $`m`$. Following $`\lambda _n`$ and $`\lambda _p`$ in time shows that the system takes an ellipsoidal form ($`\lambda _n<\lambda _p`$). We have found, as well, an anisotropy in the diffusion tensor $`\widehat{D}(t)`$ defined from $`I(t)`$ as:

$$\widehat{D}(t)=\frac{1}{Nm}\frac{I(t+\delta t)I(t)}{\delta t}$$
(5)
As example we show, in Fig. 2 the two eigenvalues $`D_1`$ and $`D_2`$ of $`\widehat{D}`$ as functions of $`t`$ for a particular periodic system. Clearly, $`D_1`$ and $`D_2`$ are very different for all $`t`$. For all systems we studied, we found two different diffusion coefficients which depend strongly the system size. We were not able to find how these values scale with $`L`$.

In addition, we have found that, contrary to its initial condition, the system starts to rotate. This fact is put in evidence by calculating the two eigenvectors $`\stackrel{}{u}_n`$ and $`\stackrel{}{u}_p`$ of $`I(t)`$. These two (perpendicular) vectors rotate in space and, most importantly, they keep the same direction of rotation for a long time ($`\tau _r`$). We suspect that this rotation induces an anomalous temporal correlation of velocities. One should point out that this rotation phenomenon seems similar to that observed by Alder et al in their simulations with similar periodic boundary conditions.

We strongly believe that the use of periodic boundary conditions is responsible for this anomalous correlation. These boundary conditions present another inconvenience which is connected with the rotation of the system: The square geometry of the system does not permit the conservation of distances between two particles when the system is rotating. In Fig. 3 we show that after a rotation of $`\theta `$ the distance $`d_{ij}`$ between particles $`i`$ and $`j`$ can be drastically changed by the rotation if one of them goes through the boundary. Because the distance between particles is not conserved by rotation, the interaction potential used in the algorithm, which depends only on the relative positions of particles $`d_{ij}`$, is itself not invariant under rotation and so the angular momentum of the system is not conserved. Effectively the total angular momentum is fluctuating as one can see in Fig. 4. Every time a particle goes through the boundary, its angular momentum, $`l_z^i`$, changes sign. Consequently, the change in angular momentum is $`\mathrm{\Delta }L_z=2l_z^i`$. $`\mathrm{\Delta }L_z`$ is always proportional to $`L`$ (the system size) and the total number of particles, $`N`$, is proportional to $`L^2`$. However it appears that the fluctuations of $`L_z`$ get bigger with system size (see Fig. 4).

The use of periodic boundary conditions amounts to replicating the system on a square lattice. There are, therefore, several identical systems which interact through the boundaries. The rotation observed in our system is then extended to all these systems and can create some shear stress, due to frustration of rotation, between neighboring systems.

These boundary conditions can have other consequences on the dynamics of granular systems. For example, during the simulation of a cooling state the system evolves towards clusters  whose orientation depends on the type of boundary conditions . In other similar simulations  it was shown that there exist large spatial correlations between particles where the velocities stay correlated over a distance of about $`L/2`$.

### B Reflecting boundaries

In the case of reflecting boundaries the system is rotationally invariant leading to better behavior of the mean square displacement. However, in this case, the mean square displacement is limited at long time by the system size. To circumvent this problem, we proceed as follows. The test particle s is initially put at the center of system at $`t=0`$. The evolution of the position and velocity of this test particle are then followed until it reaches the boundary of the system in time $`t_w`$. We then repeat this many times collecting statistics for many test particles with different initial velocities.

The mean square displacement is calculated over 500 such trajectories and limited to time smaller than the smallest $`t_w`$. In this case, as one can see in Fig. 5, there is no dependence of the mean square displacement on the system size. Therefore, we can now trust the results of our numerical simulations.

We recall that the integral of the velocity correlation function does not converge in 2D. However, in a limited range of time (see Fig. 5), the quantity $`<(\stackrel{}{r}(t+t_0)\stackrel{}{r}(t_0))^2>`$ can be approximated by a straight line and $`D`$ calculated according to Eq. (1). Therefore, the estimate of $`D`$ with this method is an approximation.

## III Diffusion in an elastic gas

We first validate our algorithm using reflecting boundaries for an elastic gas (i.e., where the collision between particles are elastic). Then, we compare the numerical results with those given by the Langevin equation. Indeed, near equilibrium, the dynamics of s can be described approximately, by a Langevin equation:

$$\begin{array}{c}\frac{d\mathrm{v}_\mathrm{i}(t)}{dt}=\gamma \mathrm{v}_\mathrm{i}(t)+\mathrm{\Gamma }_i(t),\\ <\mathrm{\Gamma }_i(t)\mathrm{\Gamma }_j(t^{})>=q\delta _{i,j}\delta (tt^{}).\end{array}$$
(6)
where $`i`$ denotes the two direction $`x`$ and $`y`$. Integrating Eq. (6), the dependence of the mean square velocity on time is simply given by:

$$v^2(t)=v^2(0)e^{2\gamma t}+\frac{q}{\gamma }(1e^{2\gamma t}).$$
(7)
In this paper $`\mathrm{v}`$ denotes the instantaneous velocity of one particle and $`v^2`$ the mean square velocity averaged over the different s trajectories. We can rewrite Eq. (7) using the mean square velocity in the equilibrium state $`v^2(\mathrm{})=q/\gamma `$:

$$v^2(t)=v^2(\mathrm{})+(v^2(0)v^2(\mathrm{}))e^{2\gamma t}.$$
(8)
where $`1/\gamma `$ corresponds to the relaxation time.

For example, if $`R_s>>R_b`$ or, equivalently, $`m_s>>m_b`$ ($`m_{s,b}`$ is the mass of the particle of radius $`R_{s,b}`$), $`1/\gamma `$ is very large: the collision of s with a light particle $`b`$ will not affect strongly the velocity of s. So a great number of collisions is needed before s reaches its equilibrium state. Knowing the total kinetic energy of the system $`E_k^{tot}`$ (which is given by the initial velocity of each particle), we can easily calculate the square velocity in the equilibrium state $`v^2(\mathrm{})=q/\gamma `$ (using the Boltzmann distribution law for elastic gases). Using the simulations to calculate $`v^2(t)`$ for the s particles for initial conditions which are very different from the stationary state, i.e. $`v^2(0)v^2(\mathrm{})`$ and comparing with Eq. (8), we obtain the relaxation time for a big (heavy) particle (Fig. 7). Note that we are looking for agreement near the equilibrium state, where Eq. (6) is valid. Indeed the dissipation term $`\gamma `$ must depend on both velocities $`v_s^2`$ and $`v_b^2`$, as we will show.

Equation (6) also gives the mean square displacement as a function of time,

$$<(\stackrel{}{r}(t)\stackrel{}{r}(0))^2>=(v_0^2\frac{q}{\gamma })\frac{(1e^{\gamma t})^2}{\gamma ^2}+\frac{2q}{\gamma ^2}t\frac{2q}{\gamma ^3}(1e^{\gamma t}).$$
(9)
Comparing Eq. (9) and Eq. (1) at large time, the coefficient of diffusion is seen to be $`D=\frac{v^2(\mathrm{})}{2\gamma }`$. In Fig. 7, we compare the theoretical mean square displacement, Eq. (9), with the numerical one, obtained for the case $`R_s=3R_b`$. Note that in Eq. (9) all the parameters are known. Clearly, the agreement is very good. This confirms that even for a large test particle, the motion is well described by the simple Langevin equation. This observation, while reasonable, is not trivial since the limited size of our system and the radius of the particles are comparable to the mean free path.

We can now present a theoretical calculation of $`\gamma `$ which describes dissipation in the Langevin equation for all pairs ($`R_s`$,$`R_b`$). This will allow us to compare the theoretical values with the numerical ones as a function of $`R_s`$.

The value of $`\gamma `$ depends on both velocities, $`v_s`$ and $`v_b`$. To estimate theoretically the value of $`\gamma `$, we consider the deviation, due to a collision, of the particle s moving at $`v_s`$ in the $`x`$ direction. The dissipative term $`\gamma \stackrel{}{v_i}`$ appearing in Eq. (6), in the $`x`$ direction, can therefore be formally written as

$$\gamma v_s=<\frac{\stackrel{}{v^{}}_s.\stackrel{}{x}v_s}{v_s}>\omega _cv_s$$
(10)
where $`\stackrel{}{v^{}}_s`$ is the velocity after the collision and $`\omega _c`$ is the rate of collision. The symbol $`<>`$ in Eq. (10) corresponds to the average over all collisions between the s particle and the b ones. To calculate the different terms, we proceed as follows. We consider the collision of s with a particle b moving at a velocity $`\stackrel{}{v}_b`$. The collision is characterized by two angles: $`\theta `$, the angle between $`(\stackrel{}{r_s}\stackrel{}{r_b})`$ and the $`x`$ axis, and $`\phi `$ the angle between $`\stackrel{}{v_b}`$ and the $`x`$ axis. Then, for such a collision, illustrated in Fig. 6, we can calculate theoretically $`\stackrel{}{v_s^{}}(\theta ,\phi )`$, the final velocity of the s particle.

For elastic collisions, the projection of $`\stackrel{}{v_s^{}}(\theta ,\phi )`$ on the $`\stackrel{}{x}`$ direction is given by,

$$\stackrel{}{v_s^{}}(\theta ,\phi ).\stackrel{}{x}=\frac{m_sm_b}{m_s+m_b}v_scos^2(\theta )+\frac{2m_b}{m_s+m_b}v_bcos(\theta )cos(\theta \phi )+v_ssin^2(\theta ),$$
(11)
with the collision taking place only if

$$v_scos(\theta )v_bcos(\theta \phi )>0.$$
(12)
Integrating over $`\phi `$ and taking into account Eq. (12) we can write

$$<\stackrel{}{v_s^{}}(\theta ).\stackrel{}{x}>_\phi =\frac{\underset{0}{\overset{2\pi }{}}_v(v_scos(\theta )v_bcos(\theta \phi ))\stackrel{}{v_s^{}}(\theta ,\phi )\stackrel{}{x}𝑑\phi }{\underset{0}{\overset{2\pi }{}}_v(v_scos(\theta )v_bcos(\theta \phi ))},$$
(13)
where $`_v`$ is the Heavyside function. We found for Eq. (13) two solutions depending on the velocities. If $`v_s<v_b`$, we have

$$<\stackrel{}{v_s^{}}(\theta ).\stackrel{}{x}>_\phi =\frac{m_sm_b}{m_s+m_b}v_scos^2(\theta )+v_ssin^2(\theta )\frac{2m_bv_bcos(\theta )sin(\theta _p)}{(\pi \theta _p)(m_s+m_b)}$$
(14)
for all $`\theta `$ ($`0\theta \pi `$) and with $`\theta _p=\mathrm{cos}^1(v_s\mathrm{cos}(\theta )/v_b)`$. For the second case, $`v_s>v_b`$, there is a critical angle, $`\theta _c=\mathrm{cos}^1(v_b/v_s)`$, such that for $`\pi \theta _c<\theta \pi +\theta _c`$ the collision does not take place. In this case the solution of Eq. (13) is

$$\begin{array}{c}<\stackrel{}{v_s^{}}(\theta ).\stackrel{}{x}>_\phi =\frac{m_sm_b}{m_s+m_b}v_scos^2(\theta )+v_ssin^2(\theta )\text{ for}\mathrm{\hspace{0.25em}0}\theta <\theta _c,\hfill \\ <\stackrel{}{v_s^{}}(\theta ).\stackrel{}{x}>_\phi =\frac{m_sm_b}{m_s+m_b}v_scos^2(\theta )+v_ssin^2(\theta )\frac{2m_bv_bcos(\theta )sin(\theta _p)}{(\pi \theta _p)(m_s+m_b)}\text{ for}\theta _c\theta \pi \theta _c.\hfill \end{array}$$
(15)
We call $`\nu (\theta )`$ the mean relative loss of velocity, $`\nu (\theta )=\frac{<\stackrel{}{v^{}}_s(\theta ).\stackrel{}{x}>_\phi v_s}{v_s}`$, averaging only over the angle $`\phi `$. In Fig. 8 we show $`\nu (\theta )`$ for the particular case $`R_s=0.25`$ and $`R_b=0.5`$, which means that $`v_s>v_b`$. Note that in our calculation, the terms $`v_s`$ and $`v_b`$ correspond to the averaged values with respect to the appropriate Maxwellian distribution. To obtain the mean value, $`\stackrel{~}{\nu }`$, of $`\nu `$, we average by integrating numerically over $`\theta `$.

To conclude the calculation of the dissipative term, $`\gamma v_s`$, we have to estimate, using Eq. (10), the collision frequency which also depends on the velocities of the two particles. A similar calculation of $`\nu `$ can be done . In the stationary state where $`v_s^2`$ and $`v_b^2`$ are constant and the distributions of the velocities are Maxwellian one can use

$$\omega _c=\chi \sqrt{\pi }(R_s+R_b)d\sqrt{v_s^2+v_b^2},$$
(16)
where $`d`$ is the density of b particles and $`\chi `$ is a correction factor which corresponds to the local radial distribution around the s particle.

In Fig. 9, we compare, for different values of $`R_s`$, the diffusion coefficient found from the simulation with the theoretical value, $`v^2(\mathrm{})/2\gamma `$, predicted by the Langevin equation combined with our analytical calculation of $`\gamma `$. The theoretical calculation of $`\gamma `$ agrees very well with the simulation results. We recall that $`1/\gamma `$ corresponds to the characteristic time for the diffusive behaviour. It is important to notice that $`\gamma `$ can be approximated by $`\omega _c`$ only when $`R_s<<R_b`$. Effectively, the calculation for $`m_s0`$ gives $`\stackrel{~}{\nu }=1`$. Larger s particles need to suffer more than one collision to lose memory of their previous condition. For $`m_s\mathrm{}`$, $`\stackrel{~}{\nu }`$ is found equal to zero. Using these methods, we find for the elastic monodisperse case ($`R_s=R_b`$) that relaxation (decorrelation) takes place after about three collisions.

The agreement between numerical results and theoretical predictions allows us to confirm our numerical algorithm.

## IV Forced system

In a real granular system dissipation occurs through collisions, a fact that must be taken into account. Experimental mechanical properties of grains (restitution and friction coefficients) and collision laws  are used in our simulations. The collisions between grains and the walls are treated with the same inelastic properties. Due to dissipation, we need to feed energy into the system to maintain the particles agitated. To accomplish this, we choose random heating : At every time step $`\delta t`$ we give a random acceleration, $`\eta _i(t)`$, in both spatial directions to each particle. The equation of motion can now be written formally as:

$$\begin{array}{c}m\frac{d\mathrm{v}_i}{dt}=F_i^c+F_i^t,\\ <F_i^t(t)F_j^t(t^{})>=m^2\delta _{i,j}\delta (tt^{})\eta _0^2.\end{array}$$
(17)
$`F_i^c`$ is the collision force acting on a particle of mass $`m`$. We chose the random acceleration, $`F_i^t/m`$, to be independent of the mass of the particle. It is given by a Gaussian noise of variance $`\eta _0^2`$.

At long time, the loss of energy due to collisions and the gain due to $`F^t`$ balance each other such that the system reaches a steady state out of equilibrium. It can be shown  that the velocity distribution in this steady state is well described by a Maxwellian.

### A Stationary state

In the stationary state energy loss and gain balance exactly. The energy loss per unit time, $`\mathrm{\Gamma }`$, for the s particle, can be expressed as:

$$\mathrm{\Gamma }=P(m_s,m_b)\omega _cm_sv^2,$$
(18)
where $`P(m_s,m_b)`$ is the relative energy loss of particle s due to collisions. Clearly as for $`\stackrel{~}{\nu }`$, $`\mathrm{\Gamma }`$ must depend on the mass of the particle and on the two velocities $`v_s`$ and $`v_b`$. On the other hand, the gain of energy due the stochastic force is

$$\frac{1}{2}m_s[v^2(t+\delta t)v^2(t)]=m_s\eta _0^2\delta t.$$
(19)
In the steady state of the monodisperse system ($`R=R_s=R_b`$, and $`v^2(\mathrm{})=`$constant), we find, using Eq. (16), the following scaling for $`v^2(\mathrm{}):`$

$$\begin{array}{c}v^2(\mathrm{})(\eta _0^2)^{2/3}\\ v^2(\mathrm{})\tau _c\end{array}$$
(20)
We checked these two scaling laws numerically (see Fig. 10) and obtained the correct exponent $`2/3`$ for the various coefficients of restitution used in the contact laws. We have also verified the predicted dependence on $`\tau _c`$ for different values of $`R`$. The good agreement between theory and simulation indicates that we can describe the system by macroscopic continuous equations if $`\delta t<<\tau _c`$. As we explain elsewhere , the term $`P(m,m)`$ (in the monodisperse case) is independent of mass and velocity, because all particles are identical. This value of $`P(m,m)`$ was found equal approximatively to 0.145 for the mechanical properties corresponding to acetate spheres . We can then, in the case of a mono-disperse system, predict the dependence of the mean square velocity on the various parameters and, consequently, characterize the stationary state. For the bi-disperse case, the calculation is more complicated. Indeed the loss of energy depends on the two types of colliding particles and also on the different coefficients of restitution and friction introduced in the collision laws. As we show  the dependence of $`P(m_s,m_b)`$ on $`v_s/v_b`$ is not trivial.

In this paper we limit ourselves to the effect of the thermalization mode (or random force) on the diffusion coefficient. To this end, we will compare in the following section the simulation results for $`D`$ with $`v^2(\mathrm{})/2\gamma `$ from the Langevin equation.

### B Diffusion of one particle

To estimate $`D`$, we use reflecting boundaries and the same method explained in section II B. We consider here the bi-disperse case (a single particle of radius $`R_s`$ in a sea of particles of radius $`R_b`$). As we have not yet found a theoretical expression for $`P(m_s,m_b)`$ for this case, we use for the mean square velocities the values obtained from the simulations which are shown in Fig. 11a. Note that $`\eta _0^2`$ has been chosen such that the value of $`v_b^2`$ is the same as in the previous section. We see that $`v^2(\mathrm{})`$ first decreases with $`R_s`$ for $`R_s<R_b`$ but then increases when $`R_s>R_b`$. Because of dissipation and the random acceleration, the repartition of the energy with the mass is no longer proportional to $`1/m_s`$. In all cases it is possible to calculate the mean collision frequency for s with Eq. (16) and the associated $`\gamma `$ value with Eq. (10). We can then calculate the relaxation time $`\tau _r`$ for all couples ($`R_s`$,$`R_b`$) used. In Fig. 11b we show the diffusion coefficient $`D`$, and the relaxation time $`\tau _r`$ -in the inset- as functions of $`R_s`$. The behavior of $`v^2(\mathrm{})`$ strongly modifies the curve of $`\tau _r`$ and $`D`$ versus $`R_s`$ compared to the elastic case. Note that the relaxation time represented in Fig. 11 clearly increases as $`R_s`$ increases.

We have seen in the elastic case that $`D=\frac{v^2}{2\gamma }`$. In Fig. 11b, we show the numerical results for $`D`$ as a function of $`R_s`$ and the corresponding values given by $`\frac{v^2}{2\gamma }`$. One can see clearly that the external noise modifies the dynamics of the granular gas and in particular the diffusion coefficient, $`D`$. The numerical value of $`D`$ is found to be larger than that obtained by the corresponding random walk. Indeed, at short time, due to the random force, $`v^2(t)`$ is not constant. Between two collisions $`v^2(t)`$ increases linearly with $`t`$. Starting with the equation of motion of particle s (between two collisions),

$$\frac{d\mathrm{v}_\mathrm{i}(t)}{dt}=\eta _i(t),$$
(21)
and with the initial conditions $`x_i(0)`$ and $`v_i(0)`$, we can calculate mean square displacement

$$<(x_i(t)x_i(0))^2>=\underset{t_1=0}{\overset{t}{}}𝑑t_1\left(\mathrm{v}_i(0)+\underset{0}{\overset{t_1}{}}\eta _i(t_1^{})𝑑t_1^{}\right)\underset{t_2=0}{\overset{t}{}}𝑑t_2\left(\mathrm{v}_i(0)+\underset{0}{\overset{t_2}{}}\eta _i(t_2^{})𝑑t_2^{}\right).$$
(22)
In two dimensions, this yields for the interval between two collisions

$$<(r(t)r(0))^2)>=v^2(0)t^2+\frac{2\eta _0^2}{3}t^3.$$
(23)
On the other hand, in the case of a random walk (or elastic collisions) the mean square displacement at short time scales as $`t^2`$. This difference explains the disagreement between $`D`$ and $`\frac{v^2(\mathrm{})}{2\gamma }`$. As the velocity changes between two collisions the probability of collision is increasing with time too. The calculation of the coefficient of diffusion is not easy in this case, due to the correlation between the velocity and the probability of collision (see Eq. 16). For very small particles, if one approaches relaxation by the time of a new collision, i.e. $`\stackrel{~}{\nu }1`$, this calculation should be possible. Indeed we can assume that the velocities before and after a collision are not correlated and have the same distribution (molecular chaos). We can then compute the mean square displacement, knowing the dependence of the collision probability on the velocity . With this assumption we improve the estimate of $`D`$ for the smallest $`R_s`$. But for bigger particles we have seen that the velocities stay correlated over many collisions and we can no longer use molecular chaos.

## V Conclusions

We have presented here some general results about the diffusion process in an agitated granular gas. We first showed that the boundary conditions used in the simulations are of crucial importance. Indeed, periodic boundary conditions introduce artificially strong temporal correlations which alter the macroscopic properties of the gas. If we ensure that no correlations are induced by the algorithm, for example by using reflecting boundaries, the numerical results obtained for an elastic gas can be described very well by a Langevin equation. We have presented a theoretical calculation of the relaxation time which allows us to predict the diffusion coefficient in all cases studied. This was not a priori intuitive since the radius of the particles is of the order of the mean free path. Finally we have analyzed the influence of uniform heating (a random acceleration) on dissipative gases. We have shown that heating influences the dynamics at short time. This is evident through the value of the diffusion coefficient which is different from that expected from the Langevin description. We are now applying with success these results to the diffusion process in a granular mixture consisting of two type of grains (differing by mass or size) in equal proportion.

ACKNOWLEDGMENTS

This work was partially funded by the CNRS Programme International de Cooperation Scientifique PICS $`\mathrm{\#}753`$ and the Norwegian research council, NFR.