format: text

text:
# Effects of columnar disorder on flux-lattice melting in high-temperature superconductors

## I Introduction

Type II superconductors tinkham ; blatter ; brandt  allow for a partial penetration of magnetic field into the bulk of the superconducting (SC) material when the applied field $`H`$ satisfies $`H_{c1}<H<H_{c2}`$. In a seminal work Abrikosov abrikosov  showed that when the ratio $`\lambda /\xi `$, where $`\lambda `$ is the magnetic field penetration depth and $`\xi `$ is the coherence length, is greater than $`1/\sqrt{2}`$ the magnetic field penetrates the SC material in the form of flux-lines (FLs). These FLs are also called vortices, since they are surrounded by circular currents. Each FL carries a quantized unit of flux $`\varphi _0=hc/2e`$ called the fluxoid. The FLs have cylindrical shape of radius $`\lambda `$ (the radius is not sharp since the magnetic field decays exponentially like $`\mathrm{exp}(r/\lambda )`$, where r is the distance from the axis) and a non-SC core of radius $`\xi `$. Due to a repulsive interaction among the Fls, they arrange themselves in a triangular lattice referred to as the vortex solid (VS). This result follows from mean-field theory.

After high-temperature superconductors were discovered in the 1980’s, it became apparent that thermal fluctuations, not included in the mean-field theory brezin , play an important role at relatively high temperatures and fields, still below $`T_c`$ and $`H_{c2}`$. These fluctuations can cause the Abrikosov lattice to melt into a disordered liquid via a first order transition (FOT), which can be roughly estimated using the Lindemann criterion known from solid state physics lindemann ; nelson3 ; pelcovits . Technologically, the melting of the FL lattice is important since the vortex liquid (VL) is not actually SC due to the dissipation caused by the FL motion when an electric current passes through the system. Pinning of FLs by naturally occurring defects in the form of vacancies, interstitials, twin and grain boundaries etc., is effective to impede FL motion in the VS phase, where the FLs form a rigid correlated network. The effectiveness of the pinning manifests itself by leading to high critical currents. In the VL phase pinning of a few vortices does not inhibit others from moving when a current is applied. Thus for practical purposes the sudden increase in resistivity occurs at the melting transition rather than when $`H=H_{c2}(T)`$ for any reasonably non-vanishing currents.

The existence of the melting transition in high-$`T_c`$ pristine materials has been established through numerous experimental exp1 ; exp2 ; cubitt ; exp3 ; exp4  and numerical numer1 ; margo ; ryu ; koshelev ; hu ; nordborg ; chin  studies. As was mentioned above, disorder in the form of points defects and sometimes more extended defects can and does occur naturally in laboratory samples. In addition artificial point defects can be induced by bombarding the sample with electrons originating from particle accelerators. Extended columnar defects in the form of linear damaged tracks piercing through the sample can be induced by heavy ion irradiation. Both naturally occurring and artificially induced defects are situated at random positions in the sample and their effective pinning strength (i.e. their interaction with FLs) can also vary from defect to defect. Thus defects play the role of quenched disorder. The adjective “quenched” refers to the immobility of these defects during experimental time scales. Introduction of disorder in terms of point defects or columnar pins affects both the properties of the solid and liquid phases and might also shift the location of the melting transition in the H-T plane blatter ; brandt . In the case of point pins, the VS phase is replaced with a Bragg glass phase giamarchi ; nattermann , characterized by quasi-long-range order. The melting transition is predicted to shift towards lower temperatures MD ; ertas ; goldschmidt . In the case of columnar pins the VS phase is replaced with a so called pinned Bose glass boseglass  where FLs are trapped by the columnar defects and the whole lattice becomes immobile. The Bose glass phase is similar to the localized phase of a two dimensional repulsive Bose gas in the presence of quenched disorder, as will be explained in more detail in the next section.

The effect of both kinds of disorders on the FLs melting has been studied experimentally in various high-temperature superconductors. Two common materials that have been extensively investigated are YBa2Cu3O7-δ (YBCO) and Bi2Sr2CaCu2O8+δ (BSCCO), both having critical temperatures ranging between 90-120 K at $`H=0`$. The main difference between these materials is their anisotropy parameter $`\gamma ^2=m_z/m_{}>1`$, where $`m_z`$ and $`m_{}`$ denote the effective masses of electrons moving along the $`c`$-axis and the $`ab`$-plane respectively. BSCCO is much more anisotropic: its anisotropy parameter $`\gamma `$ lies in the range of 50-200 compared to the range of 5-7 for YBCO blatter . This fact causes the FLs to be much “softer” or elastic. Thus in the case of BSCCO the FLs are sometimes described as a collection of loosely connected “pancakes” residing in adjacent Cu-O planes. Experimental studies on YBCO have shown a marked shift in the irreversibility line in the presence of the columnar disorder konczykowski ; civale ; samoilov ; paulius ; olsson . The irreversibility line in the H-T plane marks the onset of hysteresis effects and is located close to the melting transition on the solid phase side. For BSCCO, many experimental studies have been conducted gerhauser ; klein ; konczykowski1 ; seow ; khaykovich ; lee ; beek1 ; beek2 . The more recent ones have shown khaykovich ; lee ; beek1 ; beek2  that the melting line is not shifted when the density of columnar defects is relatively low, $`B_\varphi <<B`$, but for $`B_\varphi B`$ a shift in the position of the melting transition is observed. Here the matching field $`B_\varphi `$ is defined as $`B_\varphi =n_d\varphi _0`$ where $`n_d`$ is the density of the columnar defects and $`\varphi _0`$ is the flux quantum.

Theoretical work on columnar disorder includes Bose Glass theory boseglass . Radzihovsky radzihovsky  considered the possibility of two kinds of Bose glass phases (strongly or weakly pinned) depending on whether $`B<B_\varphi `$ or $`B>B_\varphi `$. More recently, columnar as well as point disorder were investigated by Goldschmidt goldschmidt  using replica field theory. He showed that the melting line shifts to lower temperature in the case of point disorder and to higher temperature in the case of columnar disorder.

Due to the complexity of the problem, especially in the presence of disorder, simulations have been very useful in studying the FL system. There have been many simulation studies of the vortex system in the presence of disorder. However, most simulation work has been confined to the addition of point disorder only. In particular, there is little work done on the effects of columnar disorder on the FL melting. Recent work by Wengel and Täuber wengel  concentrated on the case of high defect density region $`B_\varphi B`$ . In contrast, a recent simulation study by Sen et al. sen  uses a small density of columnar defects, each of infinite strength, but considers an extremely small magnetic field. Similarly, Nandgaonkar et al. nandgaonkar  also investigate the case of a very small magnetic field. In this paper we consider columnar defects of a finite strength ($`\eta `$) with a relatively low defect density, $`B_\varphi /B=0.2`$, and with realistic magnetic fields as used in the experiments.

## II The model

We first discuss the method implemented for YBCO:

Following Nelson nelson , we map the system of $`N`$ vortices (FLs) in a high-$`T_c`$ superconductor onto $`N`$ interacting bosons in 2-space + 1-time dimensions. Then we do a Path Integral Monte Carlo (PIMC) simulation on this system. The partition function of $`N`$ vortices can be expressed as:

$$\mathrm{\Xi }(T,L_z,N)=e^{F(T,L_z,N)/kT}=\frac{1}{N!}\underset{i=1}{\overset{N}{}}\text{D}𝐑_i(z)e^{\mathrm{}/kT}$$
(1)
where $`𝐑_i(z)`$ denotes the 2D position vector of the i’th vortex at a height $`z`$ along the $`c`$-axis, $`F`$ is the free energy as a function of temperature $`T`$, $`L_z`$ is the length of the sample along the $`z`$-direction. The London free-energy functional $`\mathrm{}`$ is given by

$$\frac{\mathrm{}}{kT}=\frac{1}{kT}_0^{L_z}𝑑z\left\{\underset{i}{}\frac{\epsilon _l}{2}\left(\frac{d𝐑_i}{dz}\right)^2+\underset{i<j}{}2\epsilon _0K_0\left(\frac{R_{ij}}{\lambda }\right)\right\},$$
(2)
where $`\epsilon _0=\varphi _0^2/`$ $`(4\pi \lambda )^2`$ is the vortex line energy per unit length, the line tension is $`\epsilon _l=\epsilon _0\mathrm{ln}(a_0/2\sqrt{\pi }\xi )/\gamma ^2`$ and $`R_{ij}(z)=|𝐑_i(z)𝐑_j(z)|`$. Here

$$a_0=\sqrt{\frac{2\varphi _0}{\sqrt{3}B}}$$
(3)
is the lattice spacing, $`B`$ is the magnetic field along the $`z`$-direction and $`\gamma `$ is the anisotropy parameter. $`K_0(R_{ij}/\lambda )`$ is the in-plane interaction potential between two FLs a distance $`R_{ij}`$ apart. $`K_0`$ denotes the modified Bessel function of the first kind. This expression for the London free-energy is an approximation that neglects the non-local interaction of real vortices and replaces it by in-plane interactions only, which is really justified if the FLs do not deviate too much from straight lines along the $`z`$-axis nordborg . With this approximation, the system of $`N`$ interacting FLs is equivalent to a system of $`N`$ bosons in $`d=2`$ dimension interacting with a pairwise potential $`K_0(R_{ij}/\lambda )`$.

The path integral representation of a system of $`N`$ bosons of mass $`m`$ each in two dimensions, interacting through a potential $`V(r)=g^2K_0(r/\lambda )`$, with $`g`$ being the strength and $`\lambda `$ being the range of the repulsive interaction, is given at finite temperature $`T_B`$ in terms of the imaginary time action

$$\frac{S}{\mathrm{}}=\frac{1}{\mathrm{}}_0^{\mathrm{}/kT_B}𝑑\tau \left\{\underset{i}{}\frac{m}{2}\left(\frac{d𝐑_i}{d\tau }\right)^2+\underset{i<j}{}g^2K_0\left(\frac{R_{ij}(\tau )}{\lambda }\right)\right\}.$$
(4)
Here $`\tau `$ is the imaginary time and $`T_B`$ is the temperature of the Bose system. We see that there is a one to one parameter mapping between the boson system and the vortex system nordborg :

$$\tau z,\mathrm{}kT,g^22\epsilon _0,\mathrm{}/kT_BL_z,m\epsilon _l,n=N/A=2/(\sqrt{3}a_0^2)=B/\varphi _0,$$

where $`n`$ is the average density of bosons (and FLs) and $`A`$ is the area of the sample. We can write the London free-energy functional in a dimensionless form as follows:

$$\frac{\mathrm{}}{kT}=_0^\beta 𝑑\tau \left\{\underset{i}{}\frac{1}{2\mathrm{\Lambda }^2}\left(\frac{d𝐑_i}{d\tau }\right)^2+\underset{i<j}{}K_0\left(\frac{R_{ij}}{\stackrel{~}{\lambda }}\right)\right\},$$
(5)
where

$$\mathrm{\Lambda }=\frac{\mathrm{}}{a_0g\sqrt{m}}=\frac{kT}{a_0\sqrt{2\epsilon _l\epsilon _0}},\beta =\frac{g^2}{kT_B}=\frac{2\epsilon _0L_z}{kT},\stackrel{~}{\lambda }=\frac{\lambda }{a_0}.$$
(6)
All lengths are measured in units of $`a_0`$ and energies in units of $`g^2`$ for bosons and $`\epsilon _0a_0`$ for FLs. We next discretize the integral along the $`z`$-axis by dividing it into $`M`$ segments:

$$\mathrm{\Xi }(\beta ,\mathrm{\Lambda },T)=\frac{1}{N!}\left(\frac{1}{2\pi \mathrm{\Lambda }^2\tau }\right)\underset{P}{}\underset{m=1}{\overset{M}{}}\underset{i=1}{\overset{N}{}}d^2𝐑_{i,m}e^{\mathrm{}[𝐑_{i,m}]/kT}$$
(7)
where $`\tau =\beta /M`$, $`m`$ labels the planes and

$$\frac{\mathrm{}[𝐑_{i,m}]}{kT}=\underset{i,m}{}\frac{\left(𝐑_{i,m+1}𝐑_{i,m}\right)^2}{2\mathrm{\Lambda }^2\tau }+\underset{i<j}{}\tau K_0\left(R_{ij,m}/\stackrel{~}{\lambda }\right).$$
(8)
We work only in the limit where $`\beta `$ is large, which amounts to taking $`L_z`$ very large and in the mapping onto $`2D`$ bosons corresponds to the ground state of the bosons at the absolute zero temperature ($`T_B=0`$). We used $`\beta =375`$ and discretized the $`z`$-axis into $`M=75`$ planes.

We used the matrix-squaring method klemm ; ceperley  to calculate the right action so that we can work with a small number of planes along the $`z`$-direction. Working with the primitive action requires the use of a large number of slicing of the $`z`$-direction and is very time consuming ceperley . The boundary conditions in the $`z`$-direction for a system of bosons are $`𝐑_i(\beta )=𝐑_j(0)`$, with all permutations $`j=P(i)`$ of the indices being allowed. This is what is meant by the summation over $`P`$ in Eq. (7). For small $`g`$ (a small repulsive interaction), which corresponds to large $`\mathrm{\Lambda }`$, we expect a Bose-Einstein condensed phase where permutations are important. This is the superfluid phase. For small $`\mathrm{\Lambda }`$, repulsion is large and permutations are rare as the Bose system is in its classical phase which is a Wigner crystal. These two phases correspond respectively to the FL liquid and solid phases. The melting line represented by $`\mathrm{\Lambda }=\mathrm{\Lambda }_m`$ is given approximately by the expression $`B_m=\mathrm{const}./T^2`$. Near $`T_c`$ this behavior is more complicated (see below).

We now discuss the method of simulations of BSCCO :

Because of the high anisotropy of the BSCCO system one can not use the simple picture given above. Here, instead, we follow a different model which can be cast in a form analogous to YBCO. We take the Lawrence-Doniach model ld  as our starting point for BSCCO. It leads to the following form of the London free-energy for inter-layer (IL) Josephson coupling feigelman ; ryu :

$$\mathrm{}_{IL}(𝐑_{i,m})=\frac{d\varphi _{0}^{}{}_{}{}^{2}}{8\pi ^3\lambda ^2}\left(1+\mathrm{ln}\left(\frac{\lambda }{d}\right)\right)\left[\frac{(|𝐑_{i,m}𝐑_{i,m+1}|)^2}{4r_{g}^{}{}_{}{}^{2}}1\right],$$
(9)
for $`|𝐑_{i,m}𝐑_{i,m+1}|<2r_g`$ , and

$$\mathrm{}_{IL}(𝐑_{i,m})=\frac{d\varphi _{0}^{}{}_{}{}^{2}}{8\pi ^3\lambda ^2}\left(1+\mathrm{ln}\left(\frac{\lambda }{d}\right)\right)\left[\frac{|𝐑_{i,m}𝐑_{i,m+1}|}{r_g}2\right],$$
(10)
for $`|𝐑_{i,m}𝐑_{i,m+1}|>2r_g`$ , where $`d`$ is the inter-layer spacing and $`r_g`$ is the healing length defined by $`r_g=\gamma d`$. For the in-plane (IP) coupling we use:

$$\mathrm{}_{IP}(𝐑_{ij,m})=\frac{d\varphi _{0}^{}{}_{}{}^{2}}{8\pi ^2\lambda ^2}K_0\left(\frac{R_{ij,m}}{\lambda }\right).$$
(11)
In principle one should add to the above interaction an electromagnetic interaction among the pancake vortices clem . The electromagnetic interaction becomes dominant in the limit of infinite anisotropy ($`\gamma \mathrm{}`$). Ryu et al. ryu  argue, (see also Ref. artemenko, ) that for the value $`\gamma =50`$ for BSCCO the Josephson coupling still dominates by one order of magnitude over the electromagnetic coupling. Their argument goes as follows: Clem clem  shows that if one has a straight array of pancake vortices along the $`z`$-axis, and one pancake is displaced a distance $`R`$ in the $`x`$-direction than the magnetic energy of the configuration increases by an amount

$$\mathrm{\Delta }E(R)=\frac{d\varphi _{0}^{}{}_{}{}^{2}}{8\pi ^2\lambda ^2}\left(𝒞+\mathrm{ln}\left(\frac{R}{2\lambda }\right)+K_0\left(\frac{R}{\lambda }\right)\right).$$
(12)
where $`𝒞`$ is Euler’s constant (=0.5772…). For large $`R`$ ($`R\lambda `$), the modified Bessel function $`K_0`$ decays exponentially and thus the energy increases like $`\mathrm{ln}(R/\lambda )`$. For small $`R`$ the Bessel function can be expanded in a power series in $`R/\lambda `$

$$K_0(R/\lambda )=\mathrm{ln}(R/2\lambda )(1+R^2/4\lambda ^2+\mathrm{})𝒞+R^2(1𝒞)/4\lambda ^2+\mathrm{},$$
(13)
and thus the magnetic energy behaves like $`R^2`$ to leading order in $`R`$ which is the same as the quadratic behavior of the Josephson energy in Eq. (9) above. The ratio of the coefficients of the quadratic terms in the magnetic and Josephson energies goes roughly like $`\gamma ^2(d/\lambda )^2`$ (where $`d/\lambda 1/100`$ for BSCCO). Thus for anisotropy $`\gamma =50`$ we get a factor of 0.25 (a somewhat more precise estimate ryu  gives a ratio of about 0.1). Thus the magnetic interaction is negligible compared to the Josephson interaction for $`\gamma =50`$. For samples with $`\gamma =200`$ these interactions are already comparable. For large values of $`R`$ the magnetic interaction increases logarithmically and the Josephson interaction increases linearly so the magnetic interaction is always negligible. The key to the estimate given above is to consider not just two pancake vortices but a whole line with one displaced pancake. This argument is valid if the deviations of the vortices from straight lines are not too large. As for the in-plane interaction Clem showed that a linear array of pancake vortices gives rise to exactly the same magnetic field at a distance $`R`$ away from it as produced by an Abrikosov vortex line. Thus Eq. (11) is consistent with the magnetic interaction of pancake vortices, again when the FLs do not deviate too much from straight lines.

Equations (9)-(11) for BSCCO can be cast in a form similar to that for the YBCO with the following substitutions:

$$\mathrm{\Lambda }=\frac{kT}{\epsilon _0a_0}\frac{r_g}{d}\sqrt{\frac{\pi }{2(1+\mathrm{ln}(\frac{\lambda }{d}))}},\tau =\frac{2\epsilon _0d}{kT}.$$
(14)
With these changes the London free-energy functional would look like

$$\frac{\mathrm{}_{IL}(𝐑_{i,m})}{kT}=\left[\frac{(|𝐑_{i,m}𝐑_{i,m+1}|)^2r_g^2}{2\mathrm{\Lambda }^2\tau }\right],$$
(15)
for $`(|𝐑_{i,m}𝐑_{i,m+1}|)<2r_g`$, and

$$\frac{\mathrm{}_{IL}(𝐑_{i,m})}{kT}=\left[\frac{(𝐑_{i,m}𝐑_{i,m+1})^2r_g^2}{2\mathrm{\Lambda }^2\tau }\right]\frac{(2r_g(|𝐑_{i,m}𝐑_{i,m+1}|))^2}{2\mathrm{\Lambda }^2\tau },$$
(16)
for $`(|𝐑_{i,m}𝐑_{i,m+1}|)>2r_g`$, when now again all lengths are measured in units of $`a_0`$. While doing simulations at a fixed magnetic field $`B`$ and temperature $`T`$ , the term $`r_g^2/2\mathrm{\Lambda }^2\tau `$ will remain constant and would drop out of $`\mathrm{\Delta }E`$ term in the Boltzmann factor. It however needs be considered during the measurement of the energy. The second term in Eq. (16) can be easily handled at the last stage of the Bisection method (see the next section for a discussion of the Bisection method).

We can make use of a reduced temperature variable to make some expressions look simpler. First, using the fact that the temperature dependence of $`\mathrm{\Lambda }`$ arises mainly through $`\epsilon _0`$ and neglecting the logarithmic corrections, one gets blatter ; margo

$$\lambda =\frac{\lambda _0}{\sqrt{1T/T_c}},\epsilon _0\frac{1}{\lambda ^2}\left(1\frac{T}{T_c}\right),$$
(17)
and hence

$$\mathrm{\Lambda }\frac{T\sqrt{B}}{1T/T_c}.$$
(18)
Defining reduced temperature as :

$$T_r=\frac{T}{1T/T_c},$$
(19)
we obtain

$$\mathrm{\Lambda }T_r\sqrt{B}.$$
(20)
This shows that the equation for the melting line is approximately

$$B_m=\mathrm{const}.(1T/T_c)^2/T^2.$$

Note that some authors anisotropy  use a temperature dependence of $`1/\sqrt{1T^2/T_c^2}`$ in $`\lambda `$, or even tinkham  $`1/\sqrt{1T^4/T_c^4}`$. All these choices coincide near $`T_c`$. The choice of temperature dependence of $`\lambda `$ is not expected to have a significant effect on the results.

## III Notes on the simulations

The technique that we use to simulate our system is called Multilevel Monte Carlo Simulation (MMC) ceperley . There are several advantages in using this technique for the simulation of the FLs over the usual metropolis Monte Carlo (MC) method. In the discrete model, we work with $`N`$ FLs with the $`z`$-axis discretized into $`M`$ planes, thus resulting in $`N`$ beads in each of the $`M`$ planes. In the usual MC method one would displace a few of these beads in a plane by small random displacements inside a two dimensional box and then would accept or reject the move based on a probability given by the Boltzmann factor.

A big disadvantage of using this technique is that it is difficult to move beads appreciably from their original positions over a number of MC steps. The reason for this is that a bead in a plane belonging to a FL finds itself in a local harmonic potential generated by the kinetic energy term involving this bead and the beads belonging to the same FL on either side of the plane. This harmonic potential becomes stronger and stronger at lower temperatures and magnetic fields. As a result, in the usual MC simulations beads keep moving around inside these local harmonic cages and end up sampling only a small part of the phase space. The other problem with the usual MC method is that there is no natural easy way of implementing FL cutting. If there are two FLs twisted around each other and if it is energetically favorable for them to reconnect each other in such a way as to lower their free energy then this step should be permitted in the MC method without regard to the question if this process occurs in reality. This is so because in the MC simulations phase space is sampled according to the probability distribution and all one needs is to generate configurations weighted by the Boltzmann factor, and the path followed in configuration space has nothing to do with any real dynamics.

These two main drawbacks of the usual MC method are easily overcome in the MMC technique. First, one moves bigger chunks of FLs encompassing beads in several planes. This way one can avoid local harmonic traps. (This is like taking an aerial route to a destination rather than going through the zigzag maze of roads.) This is much in the spirit of Fourier space Monte Carlo where one first samples modes with smaller wave numbers and then move onto higher modes.

The method of creating new FL configurations is based on the concept of the conditional probabilities. It is called the Bisection method ceperley  because one starts sampling beads by iteratively bisecting the Fls. At each stage of the division, the beads belonging to that stage are moved with some conditional probability factor $`P_i`$. It is important to make sure that the probabilities are chosen in such a way that detailed balance is satisfied at each stage of the division. One notes that the $`P_i`$’s may not be the actual Boltzmann factors for the beads to be moved at different levels. But what is required is that when all $`P_i`$’s are multiplied together, they cancel in such a way so as to leave the correct Boltzmann factor for the whole move. Thus, inherent in this algorithm is the fact that a move would finally be accepted only if it has been accepted at each stage of the Bisection method. The power of this method lies in choosing the appropriate $`P_i`$’s. If these conditional probability factors are chosen judiciously, most of the rejections would take place at the initial stages of the bisection process when not too much computational effort has been spent yet.

The cutting and reconnection of FLs is implemented naturally in a MMC method: permutation among the 3 or 4 lines chosen to be moved becomes the first among the many hierarchical steps one goes through before a move is finally accepted and the position of the beads updated accordingly. We typically moved a total of 15 to 20 beads distributed over 5 planes. Permutations were sampled by a random walk algorithm through the space of permutations nordborg (see Appendix B).

In the case of YBCO we worked with a field of $`B=4000`$ G. Working in the primitive approximation of the action would require the use of a smaller value of the dimensionless parameter $`\tau `$, which would require slicing of the $`z`$-direction into a large number of planes. To avoid this, the matrix-squaring method klemm ; ceperley  has been used to get the effective action for bigger values of $`\tau `$. For example Nordborg and Blatter nordborg  use a value of $`\tau =3.0`$ and they work with 100 planes. In this simulation a value of $`\tau =5.0`$ has been used and the $`z`$-axis has been sliced into 75 planes. Choosing a bigger value of $`\tau `$ by utilizing the matrix-squaring method makes it easier to equilibrate the system as compared to working with the primitive approximation.

For BSCCO, we did not use the matrix-squaring method because of the complications involved due to the few extra terms in free energy which contribute depending on whether $`R_{ij}`$ is smaller or bigger than $`r_g`$. Here, we used the natural inter-layer spacing $`d`$ to calculate the parameter $`\tau `$ at different temperatures and then used the MMC technique to efficiently sample the configurational space.

As mentioned previously, in the present simulations we included only the Josephson coupling. This approximation works well with YBCO but it could be less satisfactory for BSCCO because of its high anisotropy. For very anisotropic materials the electromagnetic coupling becomes important clem ; anisotropy . As discussed in Section 2, Ryu et al. ryu  estimated that for anisotropy of magnitude $`\gamma =50`$ the Josephson interaction still dominates over the electromagnetic interaction by an order of magnitude, but this will not be the case for $`\gamma =200`$ which can characterize some samples. Olson et al. olson  discuss how to include the electromagnetic interaction in a MC simulation, but they only consider the opposite limit where the Josephson coupling is totally neglected. To our knowledge there is no satisfactory treatment of both couplings included on equal footings. We carried out preliminary simulations which show that the inclusion of the electromagnetic coupling does not shift the position of the transition line much at a field of $`125`$ G, thus supporting our current conclusions. These results will be reported elsewhere sandeep .

Simulations were usually carried out for $`36`$ and $`64`$ FLs. The decay of the structure factors at the transition temperature becomes sharper when one uses a larger number of FLs. However, no appreciable shift in the transition temperature is seen while working with the smaller ($`N=36`$) or the larger system ($`N=64`$). We did not run our simulations for even larger systems as it becomes computationally very time consuming. Typical simulation times were $`3.5\times 10^6`$ MC steps. Each MC step involved moving 3-4 lines in 5 planes simultaneously. We usually averaged over 10-15 realizations of the columnar disorders, though some results have an average over as many as 20 realizations of the disorder.

Columnar disorder is modeled as an array of straight cylindrical wells of typical radius $`r_d`$=25-35 Åplaced randomly throughout the cross-sectional region of the sample and oriented along the $`z`$-direction blatter . Each columnar defect is of length $`L_z`$. The density of the columnar pins can be varied by changing their number for a given cross sectional area, and the strength is controlled by a positive dimensionless parameter $`\eta `$. If a bead happens to wander inside a columnar well, we include an extra free energy of $`\eta \frac{d\varphi _{0}^{}{}_{}{}^{2}}{8\pi ^2\lambda ^2}`$. The defect concentration was taken to be a 20 percent ratio of defects to FLs which means $`B_\varphi /B=0.2`$. The strength of disorder was set at $`\eta 0.5`$ for BSCCO. For YBCO, $`\eta =0.5`$ was found to be too large in the sense that the transition became too broad and useful information could not be extracted. Thus, for YBCO we kept $`\eta 0.3`$.

Other parameters used for YBCO are $`\gamma =5,`$ $`\lambda _0=1500`$ Å, $`\xi _{ab}(0)=15`$ Å, $`d=12`$ Å. Parameters for BSCCO are as follows: $`\gamma =50`$, $`\lambda _0=1700`$ Å, $`\xi _{ab}(0)=20`$ Å, $`d=15`$ Å.

## IV Measured Quantities

In this section we describe many different physical quantities that were monitored during the simulation. From the variation of these quantities with the temperature we can extract important conclusions about the different phases of the system.

### IV.1 Energy

In terms of the reduced temperature, the energy

$$E=T^2\frac{}{T}\mathrm{ln}(\mathrm{\Xi }(\mathrm{\Lambda },\beta ,N))$$
(21)
can be simply written as nordborg :

$$E=T_r(S_1+S_2),$$
(22)
where,

$$S_1=\underset{m,i}{}\frac{(|𝐑_{i,m}𝐑_{i,m+1}|)^2r_g^2}{2\mathrm{\Lambda }^2\tau },\mathrm{for}(|𝐑_{i,m}𝐑_{i,m+1}|)<2r_g,$$
(23)
and

$$S_1=\underset{m,i}{}\left[\frac{(𝐑_{i,m}𝐑_{i,m+1})^2r_g^2}{2\mathrm{\Lambda }^2\tau }\right]\frac{(2r_g(|𝐑_{i,m}𝐑_{i,m+1}|))^2}{2\mathrm{\Lambda }^2\tau },$$
(24)
for $`(|𝐑_{i,m}𝐑_{i,m+1}|)>2r_g`$, whereas

$$S_2=\underset{m,i>j}{}\tau K_0\left(\frac{|𝐑_{i,m}𝐑_{j,m}|}{\lambda }\right).$$
(25)
Any discontinuous jump in $`E`$ would indicate a FOT. From this discontinuous jump in energy, $`\mathrm{\Delta }E`$, we can also calculate the jump in entropy $`\mathrm{\Delta }s`$,

$$\mathrm{\Delta }s=\frac{\mathrm{\Delta }E}{T}.$$
(26)
### IV.2 Translational structure factor

The translational structure factor $`S(𝐐_{l_1,l_2})`$ is defined as,

$$S(𝐐_{l_1,l_2})=\frac{1}{MN}\underset{ij,m}{}e^{\left(i𝐐_{l_1,l_2}.(𝐑_{i,m}𝐑_{j,m})\right)}$$
(27)
where $`\mathrm{}`$ stands for the MC average, $`𝐐_{l_1,l_2}`$ is a reciprocal lattice vector and is of the form

$$𝐐_{l_1,l_2}=l_1𝐛_1+l_2𝐛_2$$
(28)
and $`l_1`$ and $`l_2`$ are some integers. $`𝐛_{1,2}`$ represent the basis vectors of the reciprocal lattice

$$𝐛_{1,2}=\frac{2\pi }{a_0\mathrm{sin}^2\theta }(𝐞_{1,2}𝐞_{2,1}\mathrm{cos}\theta )$$
(29)
where $`\theta =\pi /3`$, $`a_0`$ is the nearest neighbor distance and $`𝐞_{1,2}`$ are the unit vectors along the hexagonal unit cell such that

$$𝐞_1𝐞_2=\mathrm{cos}\theta .$$
(30)
In the simulations, $`l_1`$ and $`l_2`$ in Eq. (27) are chosen to be 1 and 0 or 0 and 1. These choices correspond to the first Bragg peak. We will normally write $`S(𝐐_{1,0})`$ as simply $`S(𝐐_1)`$. This quantity gives important information about the phase of the system. In an ordered phase where Fls sit on a triangular lattice, $`S(𝐐_1)`$ is of the order of $`N`$. In the disordered phase, it saturates to almost zero as the system size increases. In the simulations, however, there is a problem with measuring $`S(𝐐_1)`$, especially in the presence of disorder. We will return to this point at the end of the section.

### IV.3 Hexatic Structure Factor

We use Delaunay triangulation to measure the hexatic order parameter, $`\psi _6`$, which is defined as,

$$\psi _6=\underset{m=1}{\overset{M}{}}\underset{i=1}{\overset{N}{}}\frac{1}{z_{i,m}}\underset{j=1}{\overset{z_{i,m}}{}}e^{(i6\theta _{ij,m})},$$
(31)
where $`z_{i,m}`$ denotes the number of the nearest neighbors of a bead at position $`i`$, $`m`$ and it is 6 for a perfect hexagonal lattice, $`\theta _{ij,m}`$ stands for the bond angle, that is the angle that vector $`𝐑_{ij,m}`$ makes with an arbitrary axis. Just like the $`S(𝐐_1)`$, this quantity too has a large value in an ordered phase and saturates at a finite value for a system of finite size.

### IV.4 Line Entanglement

As we allow permutations of FLs, we can define a number $`N_e/N`$ as that fraction of the total number of FLs which belong to loops that are bigger than the size of a “simple” loop. A simple loop is defined as a set of $`M`$ beads connected end to end, $`M`$ being the total number of planes. Loops of size $`2M`$, $`3M`$… start proliferating at and above the transition temperature and in the corresponding 2D boson system this proliferation is related to the onset of the superfluidity.

### IV.5 Line Wandering

The transverse FL fluctuations are measured by

$$u^2(z)=(𝐑(z_0+z)𝐑(z_0))^2/a_0^2,$$
(32)
which is independent of $`z_0`$. At the transition temperature $`u^2(z)`$ undergoes a large increase for large $`z`$.

We want to emphasize one important aspect of measuring the translational structure factor. The usual way of measuring this quantity is by choosing $`𝐐_{l_1,l_2}`$ corresponding to the first Bragg peak i.e., $`𝐐_1`$ in Eq. (27). Now, it might happen that the configuration of the FLs comes close to making an almost perfect hexagonal lattice but its basis vectors are not aligned along the usual major axes of the rhombically shaped cell encompassing the system. If we use the reciprocal lattice vector $`𝐐_1`$ to measure the structure factor of such a configuration, we would end up getting a very small value for $`S(𝐐_1)`$ and might wrongly conclude that the system is in a very disordered state. This happened many times in our simulations; we got a very low value of the translational structure factor while the hexatic order was indicating a high degree of orderliness in the FL lattice. To remedy this situation the translational structure factor is measured at 60 different angles, 1 degree apart, and choose that number which gives the largest possible value of the $`S(𝐐_1(\alpha ))`$ corresponding to some angle $`\alpha `$. After implementing this technique we find that $`S(𝐐_1)`$ and $`\psi _6`$, which differ so much initially, almost follow each other.

## V Results

1. For YBCO, simulations were carried out at a magnetic field of $`B=4000`$ G. At this field we have $`2a_0\lambda _0`$. By the argument given in Ref. nordborg, , we conclude that our results should hold qualitatively for any $`B`$ such that $`a_0<\lambda `$. We checked our simulation results against those by Nordborg and Blatter nordborg . First we verified that our code was working fine by comparing our results against those given in Ref. nordborg, . Working in the limit of $`L_z\mathrm{}`$, we kept $`\tau `$ fixed as $`\mathrm{\Lambda }`$ was increased. We got a sharp transition at $`\mathrm{\Lambda }=0.0605`$. A jump was also seen in the energy as defined in the previous section and was found to be $`0.013k`$, again the same as in Ref. nordborg, .

The effect of introducing columnar disorder on the structure factor at the first Bragg peak is shown in Fig. 1. In Fig. 2 we depict the jump in the energy for a pure as well as for a disordered system.

We note that at lower strengths of the disorder, the melting line is almost unaffected and the transition takes place at $`\mathrm{\Lambda }0.0605`$. However, as we increase the strength of the disorder, the melting transition shifts towards higher values of $`\mathrm{\Lambda }`$ which means that the melting line shifts towards higher temperatures and/or higher magnetic fields. For $`\eta `$ up to $`0.1`$ no change is seen in the melting curve or in the jump in the free energy functional. However, at $`\eta `$ equal to $`0.2`$ the structure factor comes down at around $`0.0625`$ and, the jump in the energy becomes gradual. This finding is in agreement with several experimental studies where a change in the irreversibility line is seen with the introduction of columnar disorder konczykowski ; civale ; samoilov ; paulius ; olsson .

Next, we look at the FL entanglement. $`N_e`$ denotes the number of FLs which do not form a simple loop. By looking at the $`N_e/N`$ vs. $`\mathrm{\Lambda }`$ graph in Fig. 3.

it is clear that the entanglement is suppressed by as much as almost one order of magnitude just after and in the vicinity of the original transition line at $`\mathrm{\Lambda }=0.0605`$. This result is expected, as it is well known that point disorder helps in line wandering and entanglement while the columnar disorder has just the opposite effect goldschmidt  since a nearby FL is induced to align along the columnar defect, and thus its transverse fluctuations are reduced. Unfortunately, there have been few studies on columnar disorder in YBCO. Sen et al. sen  work with vanishingly small magnetic field and infinite pinning strengths, in the vicinity of the lower branch of the melting line.

Figure 4 shows the effect of columnar disorder on a system of $`64`$ FLs in YBCO.

This graph shows a transition at almost the same $`\mathrm{\Lambda }`$ as in Fig. 1. This is a confirmation that finite size effects are not important in our simulations.

2. For the BSCCO system the $`E`$ vs. $`T`$ graphs for different $`\eta `$ values are depicted in Fig. 5 (36 FLs) and Fig. 6 (64 FLs).

For the pure system ($`\eta =0`$) we see a sharp jump in the energy at exactly the point where the translational as well as hexatic structure factors sharply decline. This jump in energy is a clear signature of a first order transition. From the energy jump we can compute the jump in entropy ($`\mathrm{\Delta }s`$). The jump in energy at the transition is around $`\mathrm{\Delta }E=12.0k`$ K per vortex per layer, which gives $`\mathrm{\Delta }s=\mathrm{\Delta }E/T_m=0.19k`$ which is small compared to the experimental value of $`\mathrm{\Delta }s=0.40k`$ at $`B125`$ G. However, these values for $`\mathrm{\Delta }s`$, $`T_m`$ and $`B`$ are in qualitative agreement with the same system studied using a different model hu .

Introduction of the columnar defects of a finite strength shows some interesting effects. We put columnar pins at random positions with a concentration fixed at 20 percent of the FLs. We study the BSCCO system for pins of low strength ($`\eta =0.2`$) as well as for a higher strength ($`\eta =0.5`$). Columnar disorder of strength up to $`\eta =0.2`$ appear to have little effect on the system . This can be seen from the energy vs. temperature graph in Figs. 5 and 6, where the curve for the pure system and the one in the presence of low disorder strength ($`\eta =0.2`$) fall on top of each other, except in the transition region where the jump in case of the system with columnar pins is more gradual. Also, from Figs. 7 and 8 we see that the $`S(𝐐_1)`$ is not much affected in the presence of the columnar disorder of strength $`\eta =0.2`$. The same is true for the hexatic order depicted in Fig. 9. This means that for the lower strengths of the disorder, the vortex-vortex interaction is dominant in the range of the temperature where we carried out simulations. As a result the first order vortex-lattice melting transition (FOT) line is almost not affected.

The nature of FL melting changes, however, when the strength of the disorder is increased to $`\eta =0.5`$ as seen in Figs. 7 and 8.

We see that at the lower temperatures the order parameter is suppressed. As we increase the temperature, the order parameter starts to rise and joins the melting curve of the pure system and then falls along it at even higher temperatures. Fig. 9 indicates that the hexagonal structure factor also shows a similar behavior.

This clearly shows that the FLs start to disengage from the pinning centers as they wiggle more. At higher temperatures, columnar pins with the strength chosen here are not effective in pinning the vortex system. This becomes clear from the fact that the melting curve for the pure system and the one with columnar disorder join each other at a temperature below the FOT. Further interpretation and a possible explanation of this phenomenon will be discussed in more detail in the Conclusions section below. This convergence of the curves for the pure and disordered cases is also borne out in the $`E`$ vs. $`T`$ graphs in Fig. 5 and Fig. 6. Initially at low temperatures there is a big difference in energies of the systems with no disorder ($`\eta =0.0`$) and high disorder ($`\eta =0.5`$). However as the temperature is increased, the two curves come closer and finally merge together in the liquid phase. The jump in energy in presence of columnar defects can still be seen in the $`E`$ vs. $`T`$ graph for a disorder strength of $`\eta =0.2`$. This tends to suggest that the FOT in a BSCCO system is not affected by disorders kept at as many as 20 percent of the total number of lines used in the simulations as long as we keep the strength of the disorder less or equal to $`\eta =0.2`$. For $`\eta =0.5`$, the rise of energy with temperature becomes much more gradual and we do not see any discontinuous jump in energy at any temperature. On the other hand an abrupt change in slope of the energy vs temperature graph is observed, suggesting a discontinuity in the specific heat characterizing a second order transition.

In Fig. 10 and Fig. 11, the snapshots of FLs, projected on a plane, at temperatures less than the transition temperature and at a temperature bigger than the transition temperature are shown.

From Fig. 10 it is easily seen that at low temperatures most columnar defects have captured FLs. Also, the FLs make simple loops and have cleverly set themselves so as to make a hexagonal lattice and yet occupy as many defects as possible. The transverse fluctuations of the trapped FLs are greatly suppressed at low temperatures. At higher temperature beyond the transition point, we see that columnar defects are not occupied any more and a lot of FLs are entangled.

Inspection of snapshots like Fig. 10 gives support to the assertion of Sen et al. that the Bose glass consists of patches of ordered regions with only short range positional and orientational order. This phase is different from the Bragg-glass in systems with point pins which is characterized by quasi-long-range order.

Figure 12 shows the mean squared displacements of the FLs in the $`z`$-direction at different temperatures for a pure system.

At lower temperatures $`u^2(z)`$ saturates for large $`z`$ but in the liquid state it grows linearly koshelev . A large gap in $`u^2(z)`$ for large values of $`z`$ seen at the transition temperature signals the onset of the entanglement of the FLs.

In the presence of the columnar disorder of the strength $`\eta =0.5`$ we see from Fig. 13 that $`u^2(z)`$ at temperatures less than the transition temperature is suppressed compared to the corresponding $`u^2(z)`$ in the pure system.

This result is in agreement with the findings in Ref. goldschmidt, . Also the big gap in $`u^2(z)`$ occurring at the melting transition has moved towards a higher temperature signaling a shift in the position of the melting transition in the presence of columnar defects.

## VI Conclusions

For lower strengths of disorder ($`\eta 0.1`$ for YBCO and $`\eta 0.2`$ for BSCCO) no appreciable shift in the melting line was seen. For these strengths, a sharp drop in the translational and hexatic structure factors takes place at the transition. Also, the jump in energy at the transition is not affected much as compared with the case when there is no disorder present. This suggests that for smaller concentrations of the columnar defects the transition still remains first order. This result is in agreement with Ref. khaykovich,  where even with the introduction of the columnar disorder no shift in transition temperatures is seen as long as the concentration of columnar disorder introduced is small.

It is found that for YBCO as well as BSCCO, the melting transition shifts towards higher values of temperature and magnetic field when random disorder is introduced provided its strength exceeds a threshold which is different for YBCO and BSCCO. The size of the shift, for a given concentration of defects, depends on the strength of the disorder. For YBCO, a considerable shift in the melting line towards higher temperatures and magnetic fields was seen for $`\eta =0.2`$ and $`\eta =0.3`$. This shift was bigger for $`\eta =0.3`$ than for $`\eta =0.2`$. Similarly, for BSCCO, a large shift in the melting line towards a higher temperature was found for $`\eta =0.5`$. These results are in tune with numerous experimental findings civale ; konczykowski ; paulius ; khaykovich  as well as theoretical prediction goldschmidt  where the irreversibility line is seen to shift towards higher temperature and magnetic fields in the presence of columnar defects. The qualitative reason for this effect is that due to the interaction with the columnar defects the transverse thermal fluctuations of the FLs are reduced and thus the melting transition, as determined from the Lindemann criterion, takes place at at a higher temperature.

At $`\eta =0.5`$ the jump in energy is not discernible any more. Instead, a change of slope corresponding to a specific heat discontinuity is observed. This means that the transition is probably not of first order in the presence of columnar disorder of higher strengths, but it is rather a continuous (second order) transition.

The most dramatic outcome of this study for BSCCO is that for some values of the applied field and defect strength both the translational and hexatic structure factors start to rise at a certain temperature as the transition is approached from the lower temperature side. This is an unusual result, that to our knowledge, has not been seen in previous simulations. This fact can be explained as follows. At low temperatures, the free energy is dominated by energy effects rather than entropy considerations, and pinning effects are dominant. As a result most columnar defects capture a FL, while the rest of the FLs adjust themselves in positions such as to minimize the free energy. However, as the temperature is increased, FLs start to decouple from the defect sites, which allow the vortices pinned at interstitial positions to move themselves into a more ordered arrangement, thus increasing the structure factors compared to the situation at lower temperatures. As was mentioned before, due to its high anisotropy, FLs in BSCCO behave more like a collection of two dimensional pancakes rather than rigid rods. These pancakes are comparatively more difficult to get pinned all at once by a columnar defect. On the other hand, FLs are much stiffer in YBCO and it is easier for a columnar defect to capture a FL all along its length. It can be shown (Eq. (9.49) in Ref. blatter, ) that the depinning temperature for FLs, $`T_{dp}`$, is inversely proportional to the anisotropy parameter $`\gamma `$. Thus, we expect FL depinning to occur at comparatively smaller temperatures in BSCCO than in YBCO. Since the melting temperatures are not that different between these materials at the fields we consider, the depinning process for BSCCO occurs further below the melting temperature than in YBCO. This allows the structure factor to increase above the depinning temperature before its ultimate decline at the melting transition. For YBCO, the depinning takes place close to the melting temperature, and it is difficult to detect any rise in the structure factors, especially for small systems, because it is masked by the decrease in order due to the stronger thermal fluctuations.

In principle, The rise of the order parameters in the presence of low amount of columnar defects as the melting transition temperature is approached could be observed experimentally, if the appropriate parameters are tuned correctly. In small angle neutron scattering (SANS) cubitt ; lee  one can measure the integrated intensity over a Bragg peak of wave vector $`𝐐_1`$ which is proportional to the translational structure factor $`S(𝐐_1)`$ measured in the simulation. Another commonly used technique is muon spin rotation lee  which gives information about the width of the magnetic field distribution in the sample. This is not directly proportional to the structure factor at a given $`𝐐_1`$, but is rather given by

$$\mathrm{\Delta }B^2(T)=B^2\underset{𝐐0}{}\frac{\mathrm{exp}(𝐐^2u^2/2)}{[1+\lambda ^2(T)𝐐^2]^2},$$
(33)
where $`u^2`$ is the mean square deviation of vortices from their average position. It is possible to measure this quantity in the simulations but this was not done it in the present work. This quantity might not show the unusual rise describe above since it is dominated by $`u^2`$ which is very likely monotonically increasing with temperature yilding a monotonically decreasing line width. Thus in order to look for the effect observed in this paper we suggest using the SANS technique to look at a BSCCO sample with columnar pins described by a matching field of about 25 G and an applied field of about 125 G. It is not clear to us what is the corresponding $`\eta `$ parameter describing the pinning strength of the experimental defects. According to Blatter et al. blatter  $`\eta `$ lies in the range 0.1-1.

## VII Acknowledgements

This work is supported by the US Department of Energy (DOE), grant No. DE-FG02-98ER45686.

## Appendix A Energy sum over the images

Consider a rhombically shaped region with side $`L`$ and angle $`\theta ,`$ unit vectors are $`𝐞_1`$, $`𝐞_2`$, with $`𝐞_1𝐞_2=\mathrm{cos}\theta `$. In practice we took $`\theta =60^{}`$ but we leave the discussion here more general. The Green’s function $`G_0`$ which describe the 2D coulomb interaction between one vortex and another including all its images, as is implied by the periodic boundary conditions is given by the solution to London’s equation (see e.g. Ref. tinkham, )

$$(1\lambda ^2^2)G_0(𝐑,\lambda )=\lambda ^2\delta (𝐑),$$
(34)
with the parameter $`\lambda `$ setting the scale for the range of the interaction. The solution is given by

$$G_0(𝐑,\lambda )=\frac{2\pi \lambda ^2}{L^2\mathrm{sin}\theta }\underset{𝐐}{}\frac{\mathrm{exp}(i𝐐𝐑)}{1+\lambda ^2𝐐^2},$$
(35)
with

$$𝐑=R_1𝐞_1+R_2𝐞_2,𝐐=n_1𝐛_1+n_2𝐛_{2,}$$
(36)
where $`𝐐`$ runs over all reciprocal lattice vectors spanned by

$$𝐛_i=\frac{2\pi }{L\mathrm{sin}^2\theta }(𝐞_i𝐞_j\mathrm{cos}\theta ),$$
(37)
for $`(i,j)=(1,2),(2,1)`$ . Substituting in Eq. (35) we obtain

$$G_0(𝐑,\lambda )=\frac{\mathrm{sin}\theta }{2\pi }\underset{n_1=\mathrm{}}{\overset{\mathrm{}}{}}\underset{n_2=\mathrm{}}{\overset{\mathrm{}}{}}\frac{\mathrm{exp}[i\frac{2\pi }{L}(n_1R_1+n_2R_2)]}{L^2\mathrm{sin}^2\theta /(2\pi \lambda )^2+n_1^22n_1n_2\mathrm{cos}\theta +n_2^2}.$$
(38)
We are now going to carry out the summation over $`n_1`$ analytically. This can be done by using the formula

$$\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}f(n)=(\mathrm{residues}\mathrm{of}\pi f(z)(\mathrm{cot}\pi zi)\mathrm{at}\mathrm{poles}\mathrm{of}f(z)).$$
(39)
We have subtracted the constant $`i`$ from $`\mathrm{cot}\pi z`$ to ensure that $`|zf(z)(\mathrm{cot}\pi zi)|0`$ on the contour of integration when $`z`$ has a large negative imaginary component. The contour of integration is a square with sides parallel to the real and imaginary axes with the origin in the middle, in the limit that its size goes to infinity. In our case

$$f(z)=\frac{e^{ixz}}{(z\beta i\gamma )(z\beta +i\gamma )},$$
(40)
with appropriate values of $`\beta `$ , $`\gamma `$ , and $`x`$ . This function has two simple poles, and the residues can be easily evaluated. The final answer becomes (relabeling $`n_2`$ as $`n`$ )

$$G_0(𝐑,\lambda )=\frac{\mathrm{sin}\theta }{2}\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}\frac{\mathrm{cos}(t_2n2\pi \beta _n)\mathrm{sinh}(\gamma _nt_1)+\mathrm{cos}(t_2n)\mathrm{sinh}(\gamma _n(2\pi t_1))}{\gamma _n(\mathrm{cosh}(2\pi \gamma _n)\mathrm{cos}(2\pi \beta _n))},$$
(41)
where

$$t_1=\frac{2\pi R_1}{L},t_2=\frac{2\pi }{L}(R_1\mathrm{cos}\theta +R_2),\beta _n=n\mathrm{cos}\theta ,\gamma _n=\mathrm{sin}\theta \sqrt{n^2+L^2/(2\pi \lambda )^2}.$$
(42)
This expression is simpler than the one used by Nordborg and Blatter since it does not have different expressions for even and odd n. It also converges faster in certain regions.

## Appendix B Permutation sampling

Essentially the same method is used for permutation sampling as was used in Ref. nordborg, . The only difference is that we use permutation space of only the neighboring lines. This is so because even if we were able to get a permutation step, involving a large number of lines, accepted at the first stage of the algorithm, it would be very likely to get rejected at the following stages. So we work with only 3-5 lines. Sufficiently long segments (typically 5 planes) of a number of lines were cut. Care should be taken to make sure that chosen FLs are the nearest neighbors in the plane where the reconnection of FLs is going to take place. These points can be implemented easily with the concept of linked lists and pointers tildesley . Also, care is to be taken that even though a FL may be far away from some other FL in the rhombically shaped unit cell, it can still permute with it through one of the images of the latter. These few simple points are very important to implement the whole procedure correctly. Just for a check, we tried with the sampling procedure given in Ref. ceperley, . This gave results in good agreement with the sampling procedure given above.


text:
# Untitled Document

TAUP 2729-03

13 January, 2003

The Conformal Metric Associated

with the

U(1) Gauge of the Stueckelberg-Schrödinger Equation

O. Oron and L.P. Horwitz **Also at Department of Physics, Bar Ilan University, Ramat Gan 52900, Israel

Raymond and Beverly Sackler School of Physics and Astronomy

Tel Aviv University, Ramat Aviv 69978, Israel

Abstract: We review the relativistic classical and quantum mechanics of Stueckelberg, and introduce the compensation fields necessary for the gauge covariance of the Stueckelberg-Schrödinger equation. To achieve this, one must introduce a fifth, Lorentz scalar, compensation field, in addition to the four vector fields which compensate the action of the space-time derivatives. A generalized Lorentz force can be derived from the classical Hamilton equations associated with this evolution function. We show that the fifth (scalar) field can be eliminated through the introduction of a conformal metric on the spacetime manifold. The geodesic equation associated with this metric coincides with the Lorentz force, and is therefore dynamically equivalent. Since the generalized Maxwell equations for the five dimensional fields provide an equation relating the fifth field with the spacetime density of events, one can derive the spacetime event density associated with the Friedmann-Robertson-Walker solution of the Einstein equations. The resulting density, in the conformal coordinate space, is isotropic and homogeneous, decreasing as the square of the Robertson-Walker scale factor. Using the Einstein equations, one sees that both for the static and matter dominated models, the conformal time slice in which the events which generate the world lines are contained becomes progressively thinner as the inverse square of the scale factor, establishing a simple correspondence between the configurations predicted by the underlying Friedmann-Robertson-Walker dynamical model and the configurations in the conformal coordinates.

PACS: keywords:conformal metric, Lorentz force, generalized Maxwell theory

1. Introduction

We shall review here some of the basic foundations of a relativistically covariant classical and quantum dynamics1,2,3, in which not only space coordinates and space momenta are dynamical variables, but also the time $`t`$ and the energy $`E`$, resulting in an 8N dimensional phase space for N particles. The dynamical evolution of such a system is parametrized by a universal time $`\tau `$, essentially the time postulated by Newton, which is unaffected by motion or forces which may act on the system, and provides a universal correlation between subsystems. The existence of such a time is suggested by the early thought experiment of Einstein, for which signals emitted at some interval $`\mathrm{\Delta }\tau `$ in a frame $`F`$, according to intervals of a clock in that frame, are detected in a second inertial frame $`F^{}`$. The time of detection, recorded in terms of clocks (of the “same manufacture”) in $`F^{}`$ show that the detected interval is altered by the Lorentz time dilation. The interval in $`F^{}`$, measured on the same type of clocks ($`\tau `$) as found in the emitting frame, is then called $`\mathrm{\Delta }t`$, and participates covariantly in the Lorentz transformation. This construction would not be possible without the assumption that there are clocks of identical structure in both frames, and therefore form the basis of the assumption of a universal time. From this argument, one sees that the rate of the clocks must be the same. To be able to construct a mechanics of an N body system, in which motions of the individual parts are correlated, or to think of an evolution of the world, it was made a fundamental assumption of the theory that this invariant time $`\tau `$ is universal2. Following Stueckelberg1, consider the Hamiltonian of a free particle to be (we take $`c=1`$)

$$K=\frac{p^\mu p_\mu }{2M},$$
$`(1.1)`$
where $`p^\mu =(E,p)`$. Since $`E`$ and $`𝐩`$ are assumed independent variables, the quantity $`p^\mu p_\mu =𝐩^2E^2`$ (we use the metric $`(,+,+,+)`$) is not constrained to the constant numerical value of an a priori given mass. The Hamilton equations associated with Eq.$`(1.1)`$ are

$$\frac{dx^\mu }{d\tau }=\frac{K}{p_\mu }=\frac{p^\mu }{M},$$
$`(1.2)`$
and

$$\frac{dp_\mu }{d\tau }=\frac{K}{x^\mu }.$$
$`(1.3)`$
Dividing the space part by the time part of the equations for the $`\tau `$ derivative of $`x^\mu `$, eliminating $`d\tau `$, one finds that

$$\frac{d𝐱}{dt}=\frac{𝐩}{E},$$

which is the definition of velocity in special relativity. One finds, following a similar procedure, all of the formulas for Lorentz transformation as dynamical equations follow from the Hamilton equations generated by the free Hamiltonian $`(1.1)`$.

The quadratic form of Eq. $`(1.1)`$ makes it possible to separate variables in the two-body problem. For example, for an action-at-a-distance potential,

$$K=\frac{p_{1}^{}{}_{}{}^{\mu }p_{1}^{}{}_{\mu }{}^{}}{2M_1}+\frac{p_{2}^{}{}_{}{}^{\mu }p_{2}^{}{}_{\mu }{}^{}}{2M_2}+V(x_1x_2)$$
$`(1.4),`$
where the potential $`V`$ (of dimension mass) is a scalar function of the four-vector $`x_1x_2`$.

A model for the relativistic generalization of the (spinless) Coulomb problem was worked out classically in ref. 2, with a complete discussion and numerical solutions in ref. 4, and quantum mechanically in ref. 5, where quite general potentials as functions for the invariant $`(x_1x_2)^2`$ were studied. By the transformation

$$P^\mu =p_{1}^{}{}_{}{}^{\mu }+p_{2}^{}{}_{}{}^{\mu }X^\mu =\frac{M_1x_{1}^{}{}_{}{}^{\mu }+M_2x_{2}^{}{}_{}{}^{\mu }}{M}$$

and

$$p^\mu =\frac{M_2p_1^\mu M_1p_2^\mu }{M}x^\mu =x_1^\mu x_2^\mu ,$$

where $`M=M_1+M_2`$, one obtains

$$K=\frac{P^\mu P_\mu }{2M}+\frac{p^\mu p_\mu }{2M_{red}}+V(x),$$
$`(1.5)`$
where $`M_{red}=\frac{M_1M_2}{M_1+M_2}`$. In this form, one sees clearly that the center of mass part may be separated from the relative motion problem. By choosing a representation for the coordinates5 (called RMS coordinates) which span the spacelike part of the two subsectors of the spacelike region complementary to the light cone, one finds the usual Schrödinger spectrum for the reduced Hamiltonian. Taking the (constant) value of the generator $`K`$ at the asymptotic ionization point to be $`M/2`$ (on mass shell for both particles), for small excitation spectrum (compared to the total mass M), the center of mass energy spectrum was found to be $`Mc^2`$ plus the correct Schrödinger bound state eigenvalues, plus relativistic corrections5.

2. Stueckelberg-Schrödinger Equation

Electromagnetism may be thought of as closely associated with the gauge invariance of the nonrelativistic Schrödinger equation, in the sense that covariance of the theory under local phase transformations of the form $`\psi e^{i\mathrm{\Lambda }}\psi `$ requires the addition of compensation fields $`A_i`$ to the canonical momenta (acting as derivatives), and the field $`A_0`$ to the explicit derivative $`i\frac{}{t}`$ on the left hand side of the Schrödinger equation. By adding quadratic terms in the field strengths to the Lagrangian that produces the Schrödinger equation, one finds the inhomogeneous Maxwell equations and the equations that couple the motion of the charged particle to the fields.

In a similar way, one can demand the covariance of the Stueckelberg theory to phase transformations (functions of spacetime and $`\tau `$), and find five compensation fields, four for the spacetime derivatives, and one for the derivative in $`\tau `$ generating evolution of the system. Defining the Hamiltonian as the function which generates the $`\tau `$ derivative, i.e., putting the fifth field (as a potential) on the right hand side, one may use the Hamilton equations to derive a generalized Lorentz force. The Hamiltonian has the form

$$K=\eta ^{\mu \nu }\frac{(p_\mu ea_\mu (\xi ))(p_\nu ea_\nu (\xi ))}{2M}ea_5(\xi ),$$
$`(2.1)`$
where $`\xi ^\mu `$ are the particle coordinates in the flat Minkowski space.

The equations of motion using the Hamilton equations are then

$$\frac{d\xi ^\sigma }{d\tau }=\frac{K}{p_\sigma }=\eta ^{\sigma \nu }\frac{p_\nu ea_\nu }{M}=\frac{p^\sigma ea^\sigma }{M}$$
$`(2.2)`$
$$\frac{dp_\sigma }{d\tau }=\frac{K}{\xi ^\sigma }=e\eta ^{\mu \nu }\frac{a_\mu }{\xi ^\sigma }\frac{(p_\nu ea_\nu )}{M}+e\frac{a_5}{\xi ^\sigma }=e\dot{\xi }^\mu \frac{a_\mu }{\xi ^\sigma }+e\frac{a_5}{\xi ^\sigma }.$$
$`(2.3)`$
We now use Eq.$`(2.2)`$ to substitute for $`\dot{p}_\sigma `$, and we finally arrive at the 5D Lorentz force:

$$M\ddot{\xi }^\sigma =e\dot{\xi }^\mu f_\mu ^\sigma +ef_5^\sigma ,$$
$`(2.4)`$
where $`f_{\alpha \beta }=_\alpha a_\beta _\beta a_\alpha `$ ($`\alpha ,\beta =0,1,2,3,5`$); $`x^5=\tau `$.

We now suggest replacing the 5-potential in a flat space picture by a new Hamiltonian which contains only a 4-potential, and takes into account the fifth potential in the metric of a curved space picture . We shall designate the curved coordinates with $`\widehat{x}`$. The generator of motion is

$$K_r=g^{\mu \nu }\frac{(p_\mu ea_\mu (\widehat{x}))(p_\nu ea_\nu (\widehat{x}))}{2M}$$
$`(2.5)`$
Assuming that this functional of $`(p_\sigma ,\widehat{x}^\sigma )`$ gives Hamilton equations which are equivalent dynamically to those of the flat space we find:

$$\frac{d\widehat{x}^\sigma }{d\tau }=\frac{K_r}{p_\sigma }=g^{\sigma \nu }\frac{d\xi _\nu }{d\tau }$$
$`(2.6)`$
This equation gives us the transformation law $`d\widehat{x}^\sigma =g^{\sigma \nu }d\xi _\nu `$ and $`d\xi _\sigma =g_{\sigma \mu }d\widehat{x}^\mu `$. The second Hamilton equation gives

$$\frac{dp_\sigma }{d\tau }=\frac{K_r}{\widehat{x}^\sigma }=\frac{M}{2}\frac{g^{\mu \nu }}{\widehat{x}^\sigma }\dot{\xi }_\mu \dot{\xi }_\nu +g^{\mu \nu }\frac{a_\mu }{\widehat{x}^\sigma }\frac{(p_\nu ea_\nu )}{M}$$
$`(2.7)`$
We now replace $`\widehat{x}`$ with $`\xi `$ using the transformation law and substitute Eq.($`3`$) for $`\dot{p}_\sigma `$:

$$e\dot{\xi }^\mu \frac{a_\mu }{\xi ^\sigma }+e\frac{a_5}{\xi ^\sigma }=\frac{M}{2}\frac{g^{\mu \nu }}{\xi _\alpha }g_{\alpha \sigma }\dot{\xi }_\mu \dot{\xi }_\nu +g^{\mu \nu }e\frac{a_\mu }{\xi _\alpha }g_{\alpha \sigma }\dot{\xi }_\nu $$
$`(2.8)`$
The functionals $`K,K_r`$ are different; however, on the physical trajectories they take the same value, $`𝐊=\frac{𝐌}{\mathrm{𝟐}}`$, where K is the common numerical value. We now show that choosing a conformal metric $`g^{\mu \nu }=\mathrm{\Phi }(\widehat{x})\eta ^{\mu \nu }`$, the Lorentz force derived from $`K_r`$ is the same as the one derived from $`K`$. In this case we have:

$$g_{\mu \nu }=\frac{1}{\mathrm{\Phi }(\widehat{x})}\eta _{\mu \nu }\mathrm{\Phi }(\widehat{x})=\frac{1}{1+\frac{e}{𝐊}a_5(\widehat{x})}.$$

In this case Eq.($`2.8`$) gives

$$e\frac{a_5}{\xi ^\sigma }=\frac{M}{2}\frac{1}{\mathrm{\Phi }}\frac{\mathrm{\Phi }}{\xi _\sigma }\eta ^{\mu \nu }\dot{\xi }_\mu \dot{\xi }_\nu .$$
$`(2.9)`$
Using $`\frac{M}{2}\eta ^{\mu \nu }\dot{\xi }_\mu \dot{\xi }_\nu =\frac{𝐊}{𝚽}`$ we find that Eq.($`2.9`$) is indeed satisfied. This shows that the the Hamilton equations for the two generators are identical.

It is now interesting to examine the geodesic motion of this dynamical system, assuming the fields are static in $`\tau `$. The Lagrangian in this case is

$$L=p_\mu \dot{x}^\mu H=\frac{M}{2}g_{\mu \nu }\dot{x}^\mu \dot{x}^\nu +\dot{x}^\mu a_\mu .$$
$`(2.10)`$
We now make a small variation in $`x^\mu `$

$$x^\mu x^\mu +\delta x^\mu $$

$$\delta S=d\tau [\frac{M}{2}(\frac{g_{\mu \nu }}{x^\sigma }\dot{x}^\mu \dot{x}^\nu \delta x^\sigma +2g_{\mu \nu }\dot{x}^\mu \frac{d\delta x^\nu }{d\tau })+ea_\mu \frac{d\delta x^\mu }{d\tau }+\frac{a_\mu }{x^\sigma }\dot{x}^\mu \delta x^\sigma ]$$

From the minimal action principal we obtain, by integration by parts of the $`\tau `$ derivatives ($`\tau `$ independence of the field implies $`\frac{d}{d\tau }=\dot{x}^\mu \frac{}{x^\mu }`$)

$$0=\frac{1}{2}(\frac{g_{\mu \nu }}{x^\sigma }\dot{x}^\mu \dot{x}^\nu 2g_{\sigma \nu }\dot{x}^\mu \dot{x}^\sigma 2g_{\sigma \nu }\ddot{x}^\nu )+\frac{e}{M}(\dot{x}^\mu \frac{a_\sigma }{x^\mu }+\frac{a_\mu }{x^\sigma }\dot{x}^\mu );$$

multiplying by $`g^{\lambda \sigma }`$ we finally get

$$\ddot{x}^\lambda =\mathrm{\Gamma }_{\mu \nu }^\lambda \dot{x}^\mu \dot{x}^\nu +\frac{e}{M}\dot{x}^\mu f_\mu ^\lambda $$
$`(2.11)`$
where $`f_\mu ^\lambda =g^{\lambda \sigma }f_{\sigma \mu }`$ and $`f_{\sigma \mu }=\frac{a_\mu }{x^\sigma }\frac{a_\sigma }{x^\mu }`$.

3. The Friedmann-Robertson-Walker Universe

In the “flat space” Robertson-Walker model6 (for the spatial geometry characterized by k=0) the metric

$$ds^2=d\tau ^2\mathrm{\Phi }^2(\tau )\left(dx^2+dy^2+dz^2\right).$$
$`(3.1)`$
can be brought to the form

$$ds^2=\mathrm{\Phi }^2(t)\left(dt^2dx^2dy^2dz^2\right).$$
$`(3.2)`$
by using the transformation

$$t=\frac{d\tau }{\mathrm{\Phi }(\tau )};$$
$`(3.3)`$
$`\tau `$ is the time coordinate of a freely-falling object and therefore coincides with our notion of universal $`\tau `$. The function $`\mathrm{\Phi }(\tau )`$ is often designated by $`R`$ or $`a`$ and is the (dimensionless) spatial scale of the expanding universe. In the conformal coordinates the time-coordinate is therefore related to $`\tau `$, according to the transformation above, through

$$\frac{dt}{d\tau }=\frac{1}{\mathrm{\Phi }}.$$
$`(3.4)`$
It is interesting to use the Lorentz force in order to achieve the same result. Let us assume that $`a_5`$ depends on $`t`$ alone. In this case, the force is

$$\ddot{t}=\frac{e}{M}f_{\mathrm{\hspace{0.17em}5}}^0=\frac{e}{M}\frac{da^5}{dt}$$
$`(3.5)`$
The relation

$$\mathrm{\Phi }^2=\frac{1}{1+\frac{e}{𝐊}a^5}$$
$`(3.6)`$
then implies

$$2\frac{d\mathrm{\Phi }}{dt}\mathrm{\Phi }=\mathrm{\Phi }^4\frac{e}{𝐊}\frac{da^5}{dt},$$

i.e.,

$$\frac{da^5}{dt}=\frac{𝐊}{e}\frac{d}{dt}\left(\frac{1}{\mathrm{\Phi }^2}\right)$$

We substitute this in the force equation and multiply by $`2\dot{t}`$ to obtain

$$\frac{d\dot{t}^2}{d\tau }=2\frac{𝐊}{M}\frac{d}{d\tau }\left(\frac{1}{\mathrm{\Phi }^2}\right).$$
$`(3.7)`$
Finally, putting $`𝐊=\frac{𝐌}{\mathrm{𝟐}}`$ we arrive at the remarkable result

$$\frac{dt}{d\tau }=\frac{1}{\mathrm{\Phi }},$$

which coincides with the transformation $`(3.4)`$ from the time on the freely falling clock $`\tau `$ to the redshifted $`t`$ in the conformal form of the Robertson-Walker metric. We see that this $`t`$ corresponds to the Einstein time satisfying the dynamical Hamilton equations, and the conformal factor of the Robertson-Walker metric coincides with the conformal facter of the curved space embedding.

In this construction we have assumed the $`a_5`$ field to depend on $`t`$ alone. The generalized Maxwell equations then provide a simple connection between the Robertson-Walker scale and the event density.

The generalized Maxwell equations3 are

$$_\alpha f^{\beta \alpha }=j^\beta ,$$
$`(3.8)`$
where $`j^\beta =(j^\mu ,\rho )`$ satisfies $`_\beta j^\beta =_\mu j^\mu +_5\rho =0`$, and $`\rho `$ is the event density. In the generalized Lorentz gauge $`_\alpha a^\alpha =0`$, we have

$$_\alpha ^\alpha a^5=j^5=\rho .$$
$`(3.9)`$
Since $`a^5`$ depends on t alone $`(3.9)`$ becomes

$$_t^2a^5=\rho .$$
$`(3.10)`$
From $`(3.6)`$,

$$a_5=\frac{𝐊}{e}\left(\frac{1}{\mathrm{\Phi }^2}1\right)$$

so that from $`(3.10)`$

$$\rho =\frac{2𝐊}{e}\left[\frac{\mathrm{\Phi }_{tt}}{\mathrm{\Phi }^3}3\frac{\mathrm{\Phi }_t^2}{\mathrm{\Phi }^4}\right]$$
$`(3.11)`$
The space-time geometry is related to the density of matter $`\rho _M`$ through the Einstein equations

$$G^{\mu \nu }R^{\mu \nu }\frac{1}{2}g^{\mu \nu }R=8\pi GT^{\mu \nu },$$
$`(3.12)`$
where $`R^{\mu \nu }`$ is the Ricci tensor, $`R`$ is the scalar curvature and $`T^{\mu \nu }`$ is the energy-momentum tensor. For the perfect fluid model (isotropy implies the $`T^{\mu \nu }`$ is diagonal)

$$T^{\mu \nu }=\rho u^\mu u^\nu +P(g^{\mu \nu }+u^\mu u^\nu ).$$
$`(3.13)`$
The $`(0,0)`$ component (referring to $`\tau `$) is then

$$T^{\tau \tau }=8\pi G\rho _M.$$
$`(3.14)`$
using the affine connection derived from the metric $`(3.1)`$ one finds

$$G^{\tau \tau }=3\frac{\dot{\mathrm{\Phi }}^2}{\mathrm{\Phi }^2}=8\pi G\rho _M$$
$`(3.15)`$
and the (equal) diagonal space-space components are (for example, we write the $`x,x`$ component

$$G^{xx}=\frac{1}{\mathrm{\Phi }^2}\left[2\frac{\ddot{a}}{a}+\frac{\dot{a}^2}{a^2}\right]=8\pi GT^{xx}=\frac{8\pi GP}{\mathrm{\Phi }^2}.$$
$`(3.16)`$
Since $`T^{\tau \tau }`$ is the $`(0,0)`$ component of a tensor, it follows from $`(3.14)`$ and $`(3.4)`$ that the matter density in the conformal coordinates is given by

$$\rho _M^{}=\frac{1}{\mathrm{\Phi }^2}\rho _M.$$
$`(3.17)`$
To establish a connection between the density of events in spacetime $`\rho `$ and the density of matter(particles) $`\rho _M^{}`$, in space at a given time $`t`$, we assume that,

$$\rho _M^{}=\rho \mathrm{\Delta }t$$
$`(3.18)`$
where $`\mathrm{\Delta }t`$ is the time interval (in the conformal coordinates associated with the Stueckelberg evolution) in which the events generating the particle world lines are uniformly spread. It then follows from $`(3.18)`$ that

$$\mathrm{\Delta }t=\frac{\rho _M}{\rho \mathrm{\Phi }^2}.$$
$`(3.19)`$
We now consider two examples. For the static universe, for $`\rho _M`$ constant, it follows from Eq. $`(3.15)`$ that $`\mathrm{\Phi }`$ is given by an exponential; it then follows from $`(3.11)`$ that $`\rho `$ is constant, so that

$$\mathrm{\Delta }t\mathrm{\Phi }^2.$$
$`(3.20)`$
For the matter dominated universe, where the pressure is negligible6, one sees from $`(3.16)`$ that

$$2\mathrm{\Phi }\mathrm{\Phi }_{tt}=\mathrm{\Phi }_t^2,$$

and substituting in $`(3.11)`$, one finds after changing $`\tau `$ derivatives to $`t`$ derivatives in $`(3.15)`$ that $`\frac{\rho _M}{\rho }`$ is constant. It then follows that $`\mathrm{\Delta }t\mathrm{\Phi }^2`$ in this case as well.

This result implies that, at any given stage of development of the universe, i.e., for a given $`\tau `$, the events generating the world lines lie in an interval of the conformal time $`t`$ which becomes smaller as $`\mathrm{\Phi }`$ becomes large in the order of $`\mathrm{\Phi }^2`$ With the relation $`(3.4)`$, this corresponds, on the other hand, to a narrowing distribution, of order $`\mathrm{\Phi }^1`$ in $`\tau `$, contributing to a set of events observed at a given value of the conformal time $`t`$. In general, if one observes the configuration of a system at a given $`t`$, the events detected may have their origin at widely different values of the world time $`\tau `$ parametrizing the trajectories (world lines) of the spacetime events. It would be generally difficult to relate such configurations to the configurations in spacetime (at a given $`\tau `$, instead of at a given $`t`$) predicted by a dynamical theory. However, in this case, we see that the spreading is narrowed for large $`\mathrm{\Phi }`$, so that the set of events occurring at a given $`\tau `$ is essentially the same as the set of particles occurring at a given $`t`$. The observed configuratoins therefore become very close to those predicted by the underlying dynamical model. In the general case, the relation between $`\rho (\tau )`$ and $`\rho _M(t)`$ could be very complicated, and it may be difficult to see in the observed configurations a simple relation to the dynamical model evolving according to the world time. In the static and matter dominated Friedmann-Robertson-Walker model, the correspondence between the dynamical theory and observed configurations becomes more clear as $`\mathrm{\Phi }`$ becomes large.

4. Summary and Conclusions

We have shown that the fifth potential of the generalized Maxwell theory, obtained throught the requirement of gauge invariance of the Stueckelberg-Schrödinger equation, can be eliminated in the function generating evolution of the classical system by replacing the Minkowski metric in the kinetic term by a conformal metric. The Hamilton equations resulting from this function coincide with the geodesic associated with this metric, and with the Hamilton equations of the original form, i.e., the geodesic equations of the conformal metric describe orbits that coincide with solutions of the original Hamilton equations, as found in previous work which studied the replacement of an invariant (action-at-a-distance) potential by a conformal metric7. In this case the geodesic equations are those obtained from the conformal geodesic with the addition of a Lorentz force in standard form.

The Robertson-Walker metric can be put into conformal form. The conformal factor of the Robertson-Walker metric can then be put into correspondence with the $`a^5`$ field of the generalized Maxwell theory and therefore, through its $`t`$ derivatives (we assume no explicit $`\tau `$ dependence) with the event density. In both the static and the matter dominated models, the set of events generating the world lines of the expanding universe condense into progressively thinner slices of the conformal time. This result implies that the observed configuration of the universe at a given conformal time $`t`$, for large $`\mathrm{\Phi }`$ approximately corresponds to the configuration in spacetime predicted by the Friedmann-Robertson-Walker model at a given $`\tau `$.

References

1. E.C.G. Stueckelberg, Helv. Phys. Acta 14, 322 (1941); 14, 588 (1941).

2. L.P. Horwitz and C. Piron, Helv. Phys. Acta 46, 316 (1973).

3. M.C. Land, N. Shnerb and L.P. Horwitz, Jour. Math. Phys. 36, 3263 (1995); N. Shnerb and L.P. Horwitz, Phys. Rev. 48A, 4068 (1993).

4. M.A. Trump and W.C. Schieve, Classical Relativistic Many-Body Dynamics, Kluwer Academic, Dordrecht (1999).

5. R. Arshansky and L.P. Horwitz, Jour. Math. Phys. 30, 66,380 (1989).

6. S. Weinberg, Gravitation and Cosmology, Wiley, New York (1972); R.M. Wald, General Relativity, University of Chicago Press, Chicago (1984).

7. D. Zerzion, L.P. Horwitz and R. Arshansky, Jour. Math. Phys. 32, 1788(1991).