text:
# Kinematics of the Pencil Nebula (RCW 37) and its association with the young Vela supernova remnant RX J0852.0-4622

## 1 Introduction

RX J0852.0-4622 is a young nearby supernova remnant (SNR) recently discovered (Aschenbach, 1998; Iyudin et al., 1998) near the southeastern perimeter of the well known old Vela SNR. RX J0852.0-4622 has generated much interest since the distance and age could be as low as $`200\mathrm{pc}`$ and $`700\mathrm{yr}`$ respectively and thus it could have been generated by the nearest supernova explosion in recent human history.

The SNR was discovered in Rosat hard X-ray data (shown as contours in Fig. 1) and at these energies has a shell-like morphology. There is no obvious optical counterpart to the main body of the SNR but Redman et al. (2000) showed that a fragment of X-ray emission (labelled D/D’ by Aschenbach et al. 1995) coincides closely with a bright optical nebula, RCW 37. The X-ray fragment is located just beyond the main circular body of the remnant and is clearly visible in hard X-ray images. The main X-ray shell is not complete and there is a break in the emission in a direction coincident with that of the X-ray fragment and with RCW 37 (see Fig. 1). Redman et al. (2000) suggested that RCW 37 is physically associated with RX J0852.0-4622 and represents a venting of hot gas from the interior of the remnant to beyond the roughly circular shell as delimited in the X-ray.

RCW 37 (NGC 2736) was discovered in the 1840s by Sir John Herschel and is a bright optical nebula known to amateur astronomers as the Pencil Nebula (Fig. 2). The unusual, intricate morphology of the nebula and its large size and brightness make it surprising that there have been few studies of this object. Its apparent location at the dust-obscured eastern side of the six degree diameter old Vela SNR, away from the photogenic western filaments may be one reason. Blair et al. (1995) used the Ultraviolet Spectrometer on the Voyager 2 spacecraft to investigate O vi emission from RCW 37 and inferred that shock speeds of $`160300\mathrm{km}\mathrm{s}^1`$ are present. This is consistent with the limited line profiles obtained by Redman et al. (2000) at the edge of the nebula. The morphology of the nebula is striking and led Redman et al. (2000) to suggest that its funnel-like appearance could be taken to indicate that there is indeed a collimated flow of hot gas taking place from the remnant interior. Alternatively, a pre-existing cloud could be being shocked by escaping hot gas. The kinematics were not sufficient to establish clearly the dynamics of the nebula.

In this paper, new optical forbidden line profiles from RCW 37 from an extended slit position across the filamentary bulk of RCW 37 are presented in order to determine the structure and dynamics of the nebula. The kinematics are discussed and used to describe a model for the origin of RCW 37 that places it in the context of RX J0852.0-4622 and the main Vela SNR. It is argued that RX J0852.0-4622 is embedded within the older larger Vela SNR.

## 2 Observations and Results

Contours of the smoothed hard ($`E>1.3\mathrm{keV}`$) RASS discovery data of RX J0852.0-4622 are shown in Fig. 1 overlain on an ESO archive IIIaJ (green) image of the region. Fig. 2 is an ESO archive green image of RCW 37. Contours of smoothed Chandra ACIS data of D/D’ obtained by Plucinsky et al. (2002) are shown in Fig. 3 overlain on an ESO archive green image of RCW 37.

Spatially resolved, longslit echelle spectra of the \[O iii\] 5007Å and \[S ii\] 6716 & 6731 Å emission lines were obtained with the Manchester echelle spectrometer (MES; Meaburn et al. 1984) combined with the f/7.9 Cassegrain focus of the Anglo–Australian telescope (AAT). The observations were made on 2001 February 13-14 during reasonable seeing. A 90 Å wide interference filter was used to isolate the echelle order containing the emission lines. A Tektronix CCD with 1024$`\times `$1024 24 $`\mu `$m square pixels was the detector. A slit width of 150 $`\mu `$m ($`11\mathrm{km}\mathrm{s}^1`$ and 1″) was used. The CCD was binned by a factor of two in the spatial direction to give 512 $`\times `$ 0″.32 pixels along the slit length. Each integration was debiased and wavelength calibrated to better than $`0.5\mathrm{km}\mathrm{s}^1`$ using Thorium–Argon spectra obtained between exposures. Data were obtained from five slightly overlapping east-west slit lengths to produce an effective slit length of about 700″. The slit position is marked on Fig. 2. The data from each slit was reduced and the five were mosaiced together.

### 2.1 Kinematics

A negative gray-scale representation of a position-velocity (PV) array of \[O iii\] 5007Å line profiles is displayed in Figs 4 and 5. Data from the full effective slit length is displayed in Fig. 4 while Fig. 5 is the section of the data from the eastern side of the effective slit. The bright vertical band is the spectrum of a star intersected by the slit which can be used as a reference point between the two figures. The following features can readily be seen. Throughout the data at zero heliocentric radial velocity there is faint background \[O iii\] 5007Å emission. The bright eastern edge of the nebula begins at an offset of around 40″. The data here exhibit velocity knots and small velocity loops that correspond to the complex filaments readily seen in Fig. 2. With increasing offsets, the line profiles begin to split. At offsets of greater than about 300″ the profiles are no longer split and only a single faint component at positive velocities (up to $`120\mathrm{km}\mathrm{s}^1`$) is present. The positive velocity component then approaches the systemic with increasing offset and we identify the bright feature at around 700″ as being the western edge seen in Fig. 2. There is a negative velocity component from offsets of 550″ onwards that we attribute to the irregular diffuse emission present towards the western side of the nebula and not associated with RCW 37.

## 3 Discussion

### 3.1 Dynamics of the nebula

The position-velocity arrays show a coherent velocity structure across the nebula. An almost complete classical velocity ellipse is clearly apparent. The morphology of the nebula is reminiscent of an incomplete funnel (or tube) but could also be simply a curved sheet of emission. This is discussed below but, broadly speaking, the velocity data indicate gas that is expanding radially at about $`100\mathrm{km}\mathrm{s}^1`$. At the edge of the nebula, the motion is into the plane of the sky while towards the middle there is a velocity splitting between the near and far side of the emitting gas that increases towards the middle of the nebula. Eventually, the red-shifted emission disappears while that the blue-shifted emission returns to the systemic. At this point there is background emission (seen in Fig. 2) that contributes at negative velocities.

### 3.2 Funnel or Wavy sheet?

In the IIIaJ image in Fig. 2 RCW 37 exhibits a curved morphology but there is also a parallel line of emission to the west. This led Redman et al. (2000) to suggest that together with the eastern side, the structure may form a tube/funnel for hot gas to be vented from the remnant. Just such a structure is present in the Crab SNR (e.g. Fesen & Staker 1993) and also in the more evolved remnant DEM 34a (Meaburn, 1987). The previous kinematic data were not sufficient to test this possibility. The new kinematic data across RCW 37 are nearly exactly those expected of a funnel of circular cross-section undergoing radial expansion. However, for a funnel, there should be a similar velocity splitting from the western edge towards the central of the putative funnel as for the eastern edge and so the ionization of the funnel would have to be incomplete. Another difficulty with this interpretation is that such a funnel would be pointing to the centre of the young SNR. Also the soft RASS X-ray maps reveal a more extensive bow shaped feature and it is only the hard RASS X-ray component of this that is adjacent to RCW 37 (Aschenbach, 1998).

A plausible alternative interpretation is that the nebula gas is in the form of a thin ‘wavy sheet’ (Hester, 1987) of emission which at the eastern side overlaps along the line of sight. At the western side the single thin sheet curves to become edge on, forming the western edge. Interestingly, whether the single sheet curves into or out of the plane of sky is ambiguous kinematically. The whole system appears to be undergoing a bulk expansion generating double-peaked line profiles at the eastern side where the sheet overlaps itself.

### 3.3 Association with RX J0852.0-4622

Redman et al. (2000) showed that RCW 37 and D/D’ are spatially coincident but the RASS data (see Fig. 1) were not of high enough resolution to allow further investigation and the unlikely possibility of a chance alignment could not be ruled out. The new Chandra data of D/D’ displayed in Fig. 3 show morphological similarities to the optical data (see Plucinsky et al. 2002 for an excellent image). There is an edge to the X-ray emission from D/D’ that runs parallel to the bright optical edge of RCW 37 and the X-ray peak is close to the peak of the optical emission. A chance alignment can be ruled out in the light of this new data. Plucinsky et al. (2002) argue that the morphology of D/D’ makes it unlikely that the object is a discrete ‘bullet’ of ejecta and favour a shock break-out model where a change in density ahead of the old Vela SNR boundary has led to a localised distortion in the shock-front.

More controversial is whether, as suggested by Redman et al. (2000), RCW 37 and D/D’ are physically associated with RX J0852.0-4622. They argued that RCW 37 and D/D’ represent a venting of hot gas from the interior of the remnant to beyond the roughly circular shell as delimited in the X-ray. This view differs from that of Slane et al. (2001a) and Plucinsky et al. (2002) who argue that the X-ray spectrum of D/D’ shows enhanced abundances of O and Ne and is thermal in contrast to the spectrum of the main remnant body which is non-thermal. This led Slane et al. (2001a) to suggest that RCW 37 and D/D’ are not associated with RX J0852.0-4622. However, as Plucinsky et al. (2002) note, non-equilibrium ionisation models do not yield such enhanced abundances. Non-equilibrium models have to be used when the ionisation timescale is longer than the dynamical timescale as is likely to be the case for younger remnants. This is clearly important for RCW 37 and D/D’ since the dynamical timescale is either 12,000 years or 1000 years depending to which SNR it is attributing. Plucinsky et al. (2002) calculate an ionisation timescale range of $`7.0\times 10^{10}<n_\mathrm{e}t<7.0\times 10^{11}\mathrm{cm}^3\mathrm{s}`$ and find an upper limit of $`4.0\times 10^{11}\mathrm{cm}^3\mathrm{s}`$ by assuming an ambient density of $`1\mathrm{cm}^3`$ and a timescale of 12’000 years. However, this density is likely to be an underestimate since Redman et al. (2000) used \[S ii\] 6716 & 6731 Å line ratios to measure the local electron density within RCW 37, finding values of $`n_\mathrm{e}10^2\mathrm{cm}^3`$. The ambient density will have been up to a quarter of this value if there is little post-shock cooling but less if cooling and compression has occurred. Clearly, an ambient medium density of $`10\mathrm{cm}^3`$ and a timescale of $`1000\mathrm{yr}`$ or less would also be consistent with the Plucinsky et al. (2002) data. It seems reasonable to suggest therefore that the differing spectral properties between the main shell and D/D’ are due to non-equilibrium ionisation effects. A reason for this, that D/D’ has impacted the old wall of the Vela SNR is discussed below.

Combi et al. (1999), using the data of Duncan et al. (1996), mapped radio emission from RX J0852.0-4622 and have shown that a patch of non-thermal radio emission from the remnant beyond the X-ray shell also coincides with D/D’ and RCW 37 (see their figure 1). Interestingly, they identify a further three patches of radio emission beyond the main X-ray remnant boundary. This could indicate that the X-ray shell is fragmented and interior gas is escaping through several points. This would strengthen the case that RCW 37 and the emitting gas giving rise to the X-ray source D/D’ are indeed associated with RX J0852.0-4622. However, the region is very confused and Duncan & Green (2000) argue that two of the features that Combi et al. (1999) identify are part of an (unrelated) extended structure and that the third (feature C, in the Combi et al. 1999 nomenclature) is isolated and unconnected to RX J0852.0-4622. Duncan & Green (2000) suggest that none of the Combi et al. (1999) features are associated with RX J0852.0-4622. Follow-up observations will likely clarify this important issue.

### 3.4 RCW 37, RX J0852.0-4622 and the old Vela SNR

The velocity structure and X-ray properties described above may be explained in the following simple way. We suggest that RX J0852.0-4622 exploded within the main Vela SNR. The Vela SNR is approximately $`11,000\mathrm{yrs}`$ old and by this stage it will have swept up the interstellar medium into a cold dense shell at its boundary. The unusual optical morphology, and X-ray properties of RCW 37 and D/D’ are simply due to a localised venting of hot gas from RX J0852.0-4622 impacting the old wall of the Vela SNR. This creates an X-ray ‘hot spot’ and a curved sheet of optical emission tracing the inside edge of the old shell wall. Dubner et al. (1998) have mapped the whole Vela SNR in H i and find it to be a located in a thin shell closely correlated with the X-ray emission of Aschenbach et al. (1995). This is the wall with which we conjecture the blast wave is interacting. The hot gas may even have caused localised expansion of the wall at this point - the eastern boundary of the Vela SNR as mapped in soft X-rays by Aschenbach et al. (1995) is distorted at the position of RCW 37 and D/D’.

The differing X-ray spectral properties of feature D/D’ from the interior of the RX J0852.0-4622 led Slane et al. (2001a) to conclude that the two objects are unrelated to each other (see above). However, in the context of the above scenario, differing X-ray properties are expected since the X-ray gas in fragment D/D’ has only recently impacted the dense neutral wall of the older Vela SNR resulting in non-equilibrium ionisation. The new Chandra data are indicative of a shock break-out (Plucinsky et al., 2002) which in our model would be caused by the impact of the X-ray emitting gas with the neutral wall. Interestingly, to the north of the main body of RX J0852.0-4622 there is another much fainter X-ray fragment than D/D’. It coincides with a knot of non-thermal radio emission (feature C in Combi et al. 1999) but there is no obvious optical emission associated with it. We would expect that the X-ray spectra of this gas is more akin to that of RX J0852.0-4622 than to D/D’ since it is probably not interacting with a pre-existing obstacle.

Duncan & Green (2000) plot radio surface brightness of known galactic supernovae against diameter and find that unless RX J0852.0-4622 is located at around $`1\mathrm{kpc}`$, it is extremely faint. Duncan & Green (2000) note the unusual radio properties of RX J0852.0-4622 and, as concluded by others, suggest that the remnant must be expanding into a hot low density medium. If, as argued above, RX J0852.0-4622 is located within Vela SNR then the hot interior of the older remnant can provide just such an environment. The radio faintness of RX J0852.0-4622 could have important implications in the investigation of supernova remnants in extreme environments (such as the radio SNR in the starburst galaxy M82 by Pedlar and co-workers, see e.g. McDonald et al. 2001; Pedlar et al. 1999; Wills et al. 1997)

### 3.5 The distance to RX J0852.0-4622

Estimates of the distance to RX J0852.0-4622 are very uncertain, if it is not assumed to be physically associated with the old Vela SNR. Duncan & Green (2000) have summarised some of the difficulties. Aschenbach (1998) estimate that the lack of absorption of the X-ray flux gives an upper limit of $`1\mathrm{kpc}`$. Comparison with the surface brightness of SN1006 gives a lower limit of $`200\mathrm{pc}`$. Iyudin et al. (1998) date the remnant at $`700\mathrm{yr}`$ giving a distance of $`200\mathrm{pc}`$ but this is based on rather uncertain $`{}_{}{}^{44}\mathrm{Ti}`$ yields and an estimated shock velocity. The calcium decay products of $`{}_{}{}^{44}\mathrm{Ti}`$ were detected by Tsunemi et al. 2000, giving an age of between 630 and 970 $`\mathrm{yr}`$. Burgess & Zuber (2000) re-examined nitrate abundance data from Antartic ice cores obtained by Rood et al. (1979) which had peaks that seemed to correspond to historical supernovae. They suggested that an unidentified peak at $`700`$ could have been caused by RX J0852.0-4622 though of course, there is no indication of the location in the sky of the event. Note however that the claim by Rood et al. (1979) that SN events can be detected in this way has been questioned in subsequent studies (e.g. Herron 1982, Risbo et al. 1981). H i studies cannot constrain the distance to RX J0852.0-4622 because of the highly confused nature of this region which is close to the galactic plane (Duncan & Green, 2000).

Slane et al. (2001b) use the non-thermal spectrum of portions of RX J0852.0-4622 with a power-law model of the emission to compare derived H i column densities with other parts of the Vela SNR. Slane et al. (2001b) acknowledge the uncertainties in then scaling the differences in column density to obtain a distance but argue that the higher column density of RX J0852.0-4622 implies that it is located somewhere between the back of the Vela SNR and the Vela molecular ridge at $`12\mathrm{Kpc}`$.

The lack of an obvious compact central object within the remnant also seemed at first to indicate that RX J0852.0-4622 must be located at a large distance (Mereghetti, 2001). However, an X-ray source (but with no associated optical counterpart) has recently been discovered (Pavlov et al., 2001). It appears that the central object might be rather similar to that of Cas A - a radio quiet young isolated neutron star. However, it is worth noting that there is a 64 ms pulsar with coordinates within the supernova remnant boundary. Due to the uncertainties involved in the distance of RX J0852.0-4622 an association cannot be ruled out (Redman & Meaburn, in prep). The various distance estimates are not consistent with each other and follow-up work is clearly required to remove the uncertainties inherent in the different techniques.

Our optical study offers a clear distance constraint if RX J0852.0-4622 was generated within the Vela SNR. Cha et al. (1999) have used optical absorption lines towards a significant sample of OB stars in the direction of the Vela SNR to constrain the distance as $`250\pm 30\mathrm{pc}`$ with a conservative upper limit of $`390\pm 100\mathrm{p}\mathrm{c}`$. Clearly, if indeed RX J0852.0-4622 is located within the older Vela SNR then the distance is constrained to be $`250\mathrm{pc}`$. We note that if subsequent studies firmly indicate that RX J0852.0-4622 lies well beyond the old Vela SNR, then our model can be ruled out.

## 4 Conclusions

New \[O iii\] 5007Å line profiles of RCW 37 have been presented that show the kinematics of this nebula for the first time. A partial velocity ellipse is discovered in the pv array of line profiles. The kinematics and morphology could suggest that the structure of RCW 37 is that of a thin wavy sheet of optical emission that overlaps itself towards the eastern edge and is undergoing a systematic expansion. The western edge curves towards the line of sight but does not appear to form a complete tube or funnel of emission. We compared this feature with those found towards other SNRs. The evidence that the RCW 37 optical nebula and associated X-ray feature, D/D’ are in fact part of RX J0852.0-4622 has been discussed and we conclude that it is likely that this is the case. A simple explanation for the origin of the morphology of RCW 37 is that RX J0852.0-4622 has occured within the older, larger Vela SNR and that a portion of the supernova ejecta from RX J0852.0-4622 has impacted the pre-existing cold dense wall of the Vela SNR. The thin sheet of optical emission then traces out the inside edge of this shocked wall while the X-ray emission marks shock-heated gas. This model predicts that the distance to RX J0852.0-4622 will be similar to that of the main Vela SNR which has been recently measured to lie at $`250\pm 30\mathrm{pc}`$.

## Acknowledgements

We thank the staff of the Anglo-Australian observatory for their excellent assistance during the observations. We have made use of the Chandra Data Archive, Smithsonian Astrophysical Observatory, Cambridge MA, USA; the UKST archive, Royal Observatory Edinburgh, Scotland UK; and the ROSAT Data Archive of the Max-Plank Institut für extraterrestrische Physik at Garching, Germany. MPR and DJH are supported by PPARC.


format: text

text:
# Spin Stiffness of Stacked Triangular Antiferromagnets

## I Introduction

The magnetic ordering of geometrically frustrated antiferromagnets differs substantially from the conventional magnetic ordering in nonfrustrated magnetsKawamura (1998, 2001). Indeed, the nature of the phase transition in the case of stacked triangular antiferromagnets remains controversialGaulin (1994); Collins and Petrenko (1997); Kawamura (1992); Mailhot et al. (1994); Antonenko and Sokolov (1994); Antonenko et al. (1995); Tissier et al. (2000); Pelissetto et al. (2001a, b); Loison and Schotte (2000). The geometry of the stacked triangular lattice has triangles as the elementary units and this arrangement inhibits an anti-parallel alignment of the spins in each triangular layer. Consequently, the system is said to be geometrically frustrated. This frustration leads to a comprise where the spins on each triangle adopt a non-collinear spin ordering at low temperatures. The spins form a planar configuration in which nearest neighbours are oriented at an angle of 120 degrees with respect to one another. The ground state can be described by a matrix like order parameter giving the orientation of each spin on the elementary triangles and forms an SO(3) parameter spaceDombre and Read (1989). This unusual symmetry of the order parameter and the appearance of ’chiral’ degrees of freedom which correspond to the ground state having two possible realizations, left and right handed, lead KawamuraKawamura (1985) to conjecture the existence of a new chiral universality class . The chiral degrees of freedom are believed to be responsible for the novel critical behaviour but they are not decoupled from the spin degrees of freedom and the two quantities order simultaneously. While recent field-theoretic renormalization group studies of this system using an expansion up to six loops in fixed dimension $`d=3`$ indicate the existence of a stable fixed point that corresponds to the proposed chiral universality classPelissetto et al. (2001a, b), similar analyses using a three loop perturbation technique as well as an epsilon expansion approach to the same order, did not find a stable fixed point and hence exclude the possibility of a continuous phase transition for this frustrated systemAntonenko and Sokolov (1994); Antonenko et al. (1995). Non perturbative RG approaches find that the phase transition is possibly a very weak first order transition with effective critical exponentsTissier et al. (2000).

In the present work we use both a standard heat bath Monte Carlo method as well as a recently developed broad histogram methodde Oliveira et al. (1996) to study the classical isotropic antiferromagnet on this geometry. In particular, we study the spin stiffness which provides a direct measure of the correlation length exponent $`\nu `$. The spin stiffness is a convenient quantity to study since it does not require knowledge of the order parameter but does measure the rigidity of the ordered phase against fluctuations. Our results confirm the picture of a continuous transition belonging to a new chiral universality class.

## II Model and Methods

The model is described by the following Hamiltonian

$$H=\underset{i<j}{}J\stackrel{}{S}_i\stackrel{}{S}_j\underset{k<l}{}J^{}\stackrel{}{S}_k\stackrel{}{S}_l.$$
(1)
where $`\stackrel{}{S}_i`$ is a classical three component vector of unit length located at the sites $`i`$ of a hexagonal lattice. The first sum is over nearest neighbours in the triangular planes which interact with an antiferromagnetic coupling $`J<0`$ and the second sum is over inter-plane nearest neighbours which are taken to have a ferromagnetic coupling $`J^{}>0`$ with $`|J^{}|=|J|=1`$. Hence all energies and temperatures are measured in units of $`|J|`$.

We study the response of the system to a virtual twist of the spin system. The spin stiffness , or helicity modulusFisher et al. (1973) , measures the increase in free energy associated with twisting the order parameter in spin space by imposing a gradient of the twist angle about some axis $`\widehat{n}`$ in spin space along some direction $`\widehat{u}`$ in the lattice. The spin stiffness can be written as a second derivative of the free energy with respect to the strength of the gradient and can be calculated as an equilibrium response functionSouthern and Young (1993) . Finite size scaling theory predicts that the spin stiffness should vanish at the critical point with an exponent related to the correlation length exponent.

We calculate the diagonal elements of the spin stiffness tensor corresponding to twists about three principal directions in spin space. If we divide the lattice sites into three equivalent sublattices $`A,B`$ and $`C`$ corresponding to the vertices of the elementary triangles, then a chirality vector can be defined to characterize the non-collinear ordering of the spins. The chirality is defined locally for each upward (downward) triangle by the following expression

$$\stackrel{}{K}_{\mathrm{}}=\stackrel{}{S}_A\times \stackrel{}{S}_B+\stackrel{}{S}_B\times \stackrel{}{S}_C+\stackrel{}{S}_C\times \stackrel{}{S}_A$$
(2)
In the ground state, the chirality is uniform and perpendicular to the spin planes. This symmetry of the order parameter suggests that the average chirality direction ($`\widehat{K}`$) be chosen for one of the principal axes and the other two directions ($`\widehat{}_1,\widehat{}_2`$) are chosen to be in the spin plane perpendicular to the average chirality vector such that the three vectors form an orthonormal triad. The spin stiffness component $`\rho _\alpha `$ at temperature $`T`$ can be written as Southern and Young (1993)

$$\rho _\alpha =\frac{1}{N}\underset{i<j}{}J_{ij}(\widehat{e}_{ij}\widehat{u})^2S_i^\beta S_j^\beta +S_i^\gamma S_j^\gamma \frac{1}{NT}\left(\underset{i<j}{}J_{ij}(\widehat{e}_{ij}\widehat{u})\left[S_i^\beta S_j^\gamma S_i^\gamma S_j^\beta \right]\right)^2$$
(3)
where $`\alpha ,\beta ,\gamma =\widehat{K},\widehat{}_1,\widehat{}_2`$ and the indices are taken in cyclic order. The twist is taken to be along the $`\widehat{u}`$ direction in the lattice and $`\widehat{e}_{ij}`$ is a unit vector directed along the nearest neighbour bond from site $`i`$ to $`j`$ . The angular brackets indicate a thermal average in the canonical ensemble. Since the ground state is a planar spin arrangement, the stiffnesses satisfy a perpendicular axis theorem $`\rho _{\widehat{K}}=\rho _{\widehat{}_1}+\rho _{\widehat{}_2}`$ at zero temperature. Deviations from this relationship are a measure of fluctuations of spins from the planar order.

We perform numerical simulations using both a conventional Monte Carlo (MC) heat bath method and the more recent broad histogram method (BHM) introduced by Oliveira et. al. The latter method is based on the microcanonical ensemble approach to statistical sampling at fixed energy and allows an accurate estimate of the energy density of statesde Oliveira et al. (1996); Muñoz and Herrmann (1999); de Oliveira (1998) $`g(E)`$ . By knowing the density of states $`g(E)`$ and the microcanonical averages of various quantities $`<Q>_E`$, their temperature dependence can be determined by using the following expression for the canonical averages

$$<Q>_T=\frac{_E<Q>_Eg(E)e^{E/T}}{_Eg(E)e^{E/T}}$$
(4)
In the conventional heat bath method temperature is tuned as an external parameter and number of temperature points is limited by number of computer runs. The BHM method allows us to probe the system in a continuous range of $`T`$ but requires a large number of energy bins for large system sizes. We simulate spin systems of size $`N=L^3`$ with $`L=24,30,42,60`$ and $`66`$ for the heat bath method and only up to $`L=60`$ for the BHM method. Periodic boundary conditions are employed for both methods. We find excellent agreement between these two numerical methods.

## III Results

The broad histogram method(BHM) is based on the relation

$`g(E)<N_{up}(E)>=g(E+\mathrm{\Delta }E)<N_{dn}(E+\mathrm{\Delta }E)>`$ (5)
where $`<N_{up}(E)>,<N_{dn}(E)>`$ are microcanonical averages which measure the number of moves which increase (decrease) the energy by the amount $`\mathrm{\Delta }E`$. Once these microcanonical averages are known, the microcanonical temperature $`T_m(E)`$ can be determined from

$`1/T_m(E)`$ $``$ $`{\displaystyle \frac{d\mathrm{ln}g(E)}{dE}}`$ (6)
$``$ $`{\displaystyle \frac{1}{\mathrm{\Delta }E}}\mathrm{ln}{\displaystyle \frac{<N_{up}(E)>}{<N_{dn}(E+\mathrm{\Delta }E)>}}`$
and we can then integrate this expression in some range of energy to obtain the energy density of states $`\mathrm{ln}g(E)`$ as a function of $`E`$. In our case the energy is a continuous variable and we divide the energy axis into bins of a fixed size $`\mathrm{\Delta }E=1.8`$ such that $`\mathrm{\Delta }E<<E`$, where $`E`$ is total energy of interest. We employ a simple microcanonical dynamics to sample phase space and the energy density of states $`g(E)`$ (up to a multiplicative constant) is determined using the BHM relation above. One microcanonical sweep consists of a random sweep through the lattice sites and generating a new configuration of the spins by restricting the choice of a new random orientation of the spin at site $`i`$ with respect to the local field of the nearest neighbours such that the total energy of the system remains within the energy interval $`E,E+\mathrm{\Delta }E`$. At any given value of $`E`$ , 75 microcanonical sweeps were performed and 25 sample measurements were taken of various thermodynamic quantities such as the energy, specific heat and spin stiffness. Before sampling the next energy interval, 40 initial microcanonical sweeps were performed to avoid correlations. This procedure was repeated using different seeds for random numbers and errors were determined using the standard deviations for these separate measurements.

Figure 1 shows our results for $`\mathrm{ln}g(e)`$ as a function of $`e=E/N`$ in the case of a $`42\times 42\times 42`$ lattice. The units are arbitrary since we integrate equation (6) starting from $`e=2.1`$ and not the ground state value $`e_0=2.5`$ .The number of energy bins used for this energy range was 61740. For general values of $`L`$, the number of energy bins required to study this same range with the same fixed size of energy bin is $`5L^3/6`$ and is thus of the same order as the number of sites. When the energy density of states is combined with the microcanonical averages $`<Q>_E`$ for various thermodynamic quantities, we can then plot them as continuous functions of $`T`$ using equation (4). Figure 2 shows the energy per site and the specific heat obtained using the BHM method for various linear sizes $`L`$. The energy displays strong finite size effects near the temperature where the specific heat has a maximum. The figures clearly indicate that a phase transition occurs near $`T0.96`$ in agreement with previous MC studies.Kawamura (1998)

We have used both the BHM method and a Monte Carlo heat bath method at fixed values of $`T`$ to calculate the spin stiffness. In the heat bath method, we discard the first 1000 sweeps and perform 45000 MC steps in each run. Figure 3(a) shows both our heat bath results, indicated by points, and the BHM results, indicated by lines, for the three stiffnesses for various lattice sizes $`L`$ as a function of the temperature $`T`$. The relation $`\rho _{\widehat{K}}=\rho _{\widehat{}_1}+\rho _{\widehat{}_2}`$ is well satisfied for all values of $`T<0.95`$ indicating that there is a relatively small deviation from the planar spin configuration. All three stiffnesses are nonzero at low $`T`$ and vanish near $`T.96`$ which corresponds to the specific heat divergence in figure 2. Figure 3(b) shows the heat bath data for $`\rho _{\widehat{K}}`$ on an enlarged temperature scale. The stiffnesses clearly show large finite size effects and approach zero near $`T.96`$. The points labelled infinity are obtained by plotting $`\rho _{\widehat{K}}`$ versus $`1/L`$ at various values of $`T`$ and extrapolating to the large $`L`$ limit as shown in figure 4.

These finite size effects can be used to determine the correlation length exponent $`\nu `$ directly. Finite size scaling considerations for $`\rho (T,L)`$ predict

$`\rho (T,L)={\displaystyle \frac{1}{L}}f(L/\xi )={\displaystyle \frac{1}{L}}f(L^{1/\nu }|t|)`$ (7)
where $`t`$ is the reduced temperature. This form suggests that we can plot $`L\rho (T,L)`$ versus $`T`$ to identify $`T_c`$ as the temperature where the curves for different values of $`L`$ intersect. Figure 5 shows our heat bath results for $`L\rho _K`$ as a function of $`T`$ for lattice sizes $`L=24,30,42,60,66`$. Linear interpolations of neigbouring temperature points indicate that the curves intersect at a value of $`T_c=0.958\pm 0.002`$. We have also used our BHM results in the same temperature range and we obtain the same estimate for $`T_c`$.

In the limit as $`L\mathrm{}`$, the scaling form predicts $`\rho |t|^\nu `$. Using the values of the stiffness obtained by extrapolating to large values of $`L`$ as in figure 4 and then plotting these versus $`|t|`$ on a ln-ln scale, we can obtain an estimate of $`\nu `$. Figure 6 shows our results for $`\rho _K`$ which yields the value $`\nu =.589\pm .007`$. This value agrees very well with previous Monte Carlo estimatesKawamura (1998) but is slightly larger than the value found by the recent six loop renormalization group calculations in three dimensions.Pelissetto et al. (2001a, b)

Figure 7 shows a finite size scaling plot of our stiffness results using the values of $`T_c`$ and $`\nu `$ quoted above. The data obtained from both the heat bath MC method for sizes $`L=24,30,42,60,66`$ and the BHM method for $`L=24,30,42`$ collapse very well to a universal function for temperatures below $`T_c`$. The value $`\nu =.589\pm .007`$ is certainly very different from the value $`\nu =0.7113`$ which describes the three dimensional Heisenberg universality class.Campostrini et al. (2002)

## IV Summary

We have calculated the spin stiffness of the isotropic Heisenberg antiferromagnet on the stacked triangular geometry using both a MC heat bath and BHM method. The spin stiffness has the advantage that it measures the rigidity of the ordered phase in response to a virtual twist without having to specify the order parameter. The results obtained from both numerical approaches agree and predict a continuous phase transition which belongs to the new chiral universality class proposed by Kawamura.

###### Acknowledgements.

This work was supported by the Natural Sciences and Research Council of Canada and the HPC facility at the University of Manitoba.


text:
# Angle-resolved photoemission spectroscopy of Na-doped Ca2CuO2Cl2 single crystals: Fingerprints of a magnetic insulator in a heavily underdoped superconductor

\[

## Abstract

Electronic evolution from an antiferromagnet to a high-$`T_\mathrm{c}`$ superconductor is revealed by angle-resolved photoemission experiments on tetragonal Ca1.9Na0.1CuO2Cl2 single crystals, which were successfully grown for the first time under high pressures. In this underdoped superconductor, we found clear fingerprints of the parent insulator: a shadow band and a large pseudo-gap. These observations are most likely described by a “chemical potential shift”, which contrasts clearly with the prevailing wisdom of the “pinned chemical potential” learned from the prototype La2-xSrxCuO4, demonstrating that the route to a high-$`T_\mathrm{c}`$ superconductor is not unique.

PACS numbers: 74.25.Jb, 71.18.+y, 74.72.Jt, 79.60.-i

\]

Elucidating the mechanism of high-$`T_c`$ superconductivity (HTS) has been one of the most attractive challenges in condensed matter physics today. A key ingredient for HTS is the two dimensional CuO2 plane, which, upon carrier doping, switches from an antiferromagnetic insulator to HTS. An important step towards solving the HTS puzzle is to understand the electronic evolution across the insulator to metal transition. Angle-resolved photoemission spectroscopy (ARPES) is one of the most powerful tools to do this, as it can measure the momentum dependent electron excitation spectrum directly. So far, because of their superior surface quality associated with an easy cleavage plane, most ARPES studies have been conducted on Bi2Sr2Can-1CunOy ($`n`$ = 1, 2). In Bi cuprates, however, it is not easy to access the composition in the vicinity of the insulator-metal transition. This has been preventing us from attacking the fundamental issue of electronic evolution from parent insulator to HTS. At present, most information regarding this issue stems from La2-xSrxCuO4, where two electronic components and pinning of the chemical potential have been observed . It is of great interest to study whether the results from La2-xSrxCuO4 are universal.

Alkaline-earth copper oxychlorides, A2CuO2Cl2 (A = Ca, Sr), have been known as an ideal parent insulator. They have a simple tetragonal K2NiF4-type structure with undistorted and single CuO2 planes. Their single crystals can be cleaved very easily, and, indeed, the best ARPES data on insulators have been obtained on these systems ; e.g., a $`d`$-wave-like dispersion of the Mott gap was disclosed . Due to the success in the parent insulator, single crystals of oxychlorides doped with mobile carriers have been envisaged to explore the evolution of the Fermi surface from the $`d`$-wave-like Mott gap . Hole-doped Ca2CuO2Cl2 can be realized by substituting Na+ for Ca2+. Ca2-xNaxCuO2Cl2 shows superconductivity ($`T_{\mathrm{c}}^{}{}_{}{}^{\mathrm{max}}`$ = 28 K at $`x`$ 0.2) but forms only under high pressures of several GPa . Previous attempts to grow single crystals of Ca2-xNaxCuO2Cl2 have been unsuccessful, mainly because of the difficulty in the crystal growth under high pressures. Recently, we have succeeded in growing Ca2-xNaxCuO2Cl2 single crystals for the first time using a specially designed high pressure furnace .

In this Letter, we report the first ARPES data on superconducting, Ca1.9Na0.1CuO2Cl2 single crystals, which demonstrate that the band dispersion of the doped oxychloride is a consequence of the valence band of the parent insulator shifting to the chemical potential. The resulting fingerprints of the parent insulator manifest itself in the way of a shadow band and a large pseudo-gap. These results are remarkably different from the prevailing picture proposed for the prototypical HTS material, La2-xSrxCuO4, where the chemical potential remains fixed while new states are created around it by doping. This indicates that the electronic evolution across the insulator to metal transition is not unique among HTS’s.

The Ca2-xNaxCuO2Cl2 crystals used in this study were prepared by a flux growth . A powder mixture of Ca2CuO2Cl2 \+ 0.2NaCl + 0.2NaClO4 was heated and slowly cooled under a pressure of 4 GPa. Thin plate like crystals, shown in the inset of Fig. 1, were obtained from the solidified melt. They were easily cleaved like mica, which was the key to the success of the present ARPES study. Their Na content, $`x`$, was estimated to be 0.1 from x-ray diffraction and electron-probe-micro-analysis measurements . Superconductivity in these crystals was confirmed by both resistivity (Fig. 1) and magnetization measurements. The transition temperature of $`T_\mathrm{c}`$ = 13 K is in good agreement with $`x`$ 0.1 powder results . Since the onset of superconductivity in Ca2-xNaxCuO2Cl2 is reported to occur at $`x0.07`$ , the samples with $`x`$ = 0.1 are located just beyond the insulator-metal transition, corresponding to a heavily underdoped superconductor.

The ARPES measurements were performed at beamline 5-4 of the Stanford Synchrotron Radiation Laboratory. The data were collected using a SES-200 hemispherical analyzer with a monochromated synchrotron beam (16.5 $``$ 25.5 eV photons) or a He lamp (21.2 eV photons). The angular and energy resolution was 0.25 and $``$16 meV, respectively. At the end of the ARPES measurements, low-energy electron diffraction (LEED) was measured at 20 K. The four-fold symmetry of the LEED pattern without any super-lattice diffraction (the inset of Fig. 1, bottom right) confirms that the crystal structure is tetragonal and that there is no surface reconstruction also in the Na-doped Ca2CuO2Cl2 compound.

We begin with an examination of the raw ARPES data at 20 K along the high symmetry directions as shown in Figure 2. The spectra in Fig. 2(a) show that a peak is moving up towards the Fermi level $`E_\mathrm{F}`$ on going from (0,0) to ($`\pi `$,$`\pi `$). Slightly before ($`\pi /2`$,$`\pi /2`$), the peak reaches $`E_\mathrm{F}`$ and a sharp cutoff is clearly seen in the spectra. This is evidence for a Fermi level crossing, indicating that the Na-doping successfully changes the oxychloride into a metal. However, taking a further look at Fig. 2(a), we notice something unusual is happening around ($`\pi `$/2,$`\pi `$/2). After the Fermi level crossing, a broad feature remains and is pulling away from $`E_\mathrm{F}`$, though losing weight rapidly. The presence of a band dispersing back near ($`\pi `$/2,$`\pi `$/2), a shadow band, is strongly reminiscent of that in the parent insulator, where a band folding occurs at ($`\pi `$/2,$`\pi `$/2) because of the unit cell doubling due to an antiferromagnetic ordering. Because no change in the symmetry of the crystal structure was confirmed by the LEED measurements \[inset of Fig. 1\], the most likely source of the band folding is antiferromagnetic correlations still remaining in the carrier-doped metallic samples. Due to the presence of a structural modulation, the origin of the shadow band previously reported in Bi2Sr2CaCuOy  has been controversial. In this context, the present observation unambiguously shows that, even without any lattice distortion, the shadow band can be present. In addition to the magnetic shadow band, we notice that the fine structure of the spectra near ($`\pi /2`$,$`\pi /2`$) is more complicated than the simple Fermi cutoff. As indicated by tick-marks in Fig. 2(a), the spectra very close to the Fermi crossing consists of not a single component but two components: a sharp quasiparticle peak and a broad hump, which are quite similar to the “peak-dip-hump” structure observed in other HTS materials . We will discuss this peak-dip-hump like structure in more detail later.

We observe unusual features also around ($`\pi `$,0): almost complete loss of spectral weight at $`E_\mathrm{F}`$ over an extremely large energy scale of a tenth of an eV. The ARPES spectra from (0,0) to ($`\pi `$,0) to ($`\pi `$,$`\pi `$) are displayed in Fig. 2(b). A broad peak, which is very weakly dependent on the momentum near ($`\pi `$,0), is apparent. In the case of the hole-like Fermi surface expected from band theory, schematically shown in Fig. 2, the dispersion should reach $`E_\mathrm{F}`$ on going from ($`\pi `$,0) to ($`\pi `$,$`\pi `$). In contrast to the unambiguous Fermi level crossing in the diagonal direction \[Fig. 2(a)\], such behavior is not observed in Fig. 2(b). Instead, the broad feature simply appears to lose weight and vanishes completely on its approach to $`E_\mathrm{F}`$. The lack of a well-defined Fermi surface near ($`\pi `$,0) has been ascribed to a pseudo-gap. Remarkably, in comparison with other HTS’s, the energy scale of the pseudo-gap in Ca1.9Na0.1CuO2Cl2 is extremely large and distinct as a HTS material. The leading edge midpoint of the spectra at ($`\pi `$,0), which has often been used as a measure of the pseudo-gap, gives a value of 130 meV, which should be compared with 30 meV for Bi2Sr2Ca1-xDyxCu2Oy at a similar doping level ($`x`$ = 0.1) . The value of the leading edge midpoint along ($`\pi `$,0) to ($`\pi `$,$`\pi `$) reaches its minimum of 65 meV around ($`\pi `$,$`\pi `$/3). Note that the observed pseudo-gap magnitude is almost identical with the energy scale of the $`d`$-wave-like dispersion in the parent insulator .

To better visualize the band dispersion, the intensity of the second derivative of the spectra in Figs. 2(a) and 2(b) is mapped on the momentum vs. energy plane in Fig. 2(c). Around $`E_\mathrm{F}`$, we see more clearly the band dispersing back near ($`\pi `$/2,$`\pi `$/2) (the shadow band) as well as the missing weight near ($`\pi `$,0) (the extremely large pseudo-gap). We argue that these anomalous features seen in the Ca1.9Na0.1CuO2Cl2 crystal originate from the parent insulator. To make our point, it may be instructive to compare the Ca1.9Na0.1CuO2Cl2 with the parent insulator. In Fig. 2(c), we also reproduce the band dispersion of the parent insulator. Immediately we notice that these two sets of dispersion are nearly identical and overlap with each other if we shift the dispersion of the parent insulator by $``$700 meV. This indicates that, upon doping, the chemical potential simply drops to the top of the valence band, resulting in the shadow band near ($`\pi `$/2,$`\pi `$/2) and the large pseudo-gap around ($`\pi `$,0). This provides very clear evidence that the “large” pseudo-gap in the heavily underdoped region is intimately linked with the $`d`$-wave-like dispersion of the Mott gap.

In the rigid band shift scenario, the dispersion closest to $`E_\mathrm{F}`$ in the parent insulator \[i.e., near ($`\pi `$/2,$`\pi `$/2)\] touches $`E_\mathrm{F}`$ to form the Fermi surface by carrier-doping, switching its state into a metal (superconductor). As seen in Fig 2(a), the Fermi crossing is indeed observed near ($`\pi `$/2,$`\pi `$/2) in Ca1.9Na0.1CuO2Cl2, with additional low-energy structure in the spectra. It is found that the peak-dip-hump like structure in Ca1.9Na0.1CuO2Cl2 is robust at temperature far above $`T_\mathrm{c}`$ = 13 K \[Fig. 3(a)\]. The peak-dip-hump like structure observed in Ca1.9Na0.1CuO2Cl2 differs from that in Bi2Sr2CaCu2Oy in two ways . Here it is observed in the nodal direction and above $`T_\mathrm{c}`$, while in Bi2Sr2CaCu2Oy it is observed at ($`\pi `$,0) and below $`T_\mathrm{c}`$. There have been several proposals for the origin of the peak-dip-hump structure in the spectral lineshape of Bi2Sr2CaCu2Oy, many resulting from a coupling of the quasiparticles with some collective modes . The present results are quite different from the drastic change in the ($`\pi `$,0) spectra across $`T_\mathrm{c}`$ in Bi2Sr2CaCu2Oy , and further experiments such as the doping dependence study will be helpful in elucidating the origin of this feature.

Now, let us examine the overall distribution of the low-energy excitations in momentum space. To distinguish different energy scales of contributions, we adopted two methods. One is the gray scale intensity plot obtained by integrating the spectral-weight over a 100 meV window below $`E_\mathrm{F}`$ as shown in Figure 3(b). The majority of the intensity lies about ($`\pi `$/2,$`\pi `$/2), mainly because the large pseudo-gap removes most of the low-energy spectral weight around ($`\pi `$,0). This intensity plot with a large energy-integration window gives a rather broad structure around ($`\pi `$/2,$`\pi `$/2), which extends not only towards (0,0) but also towards ($`\pi `$,$`\pi `$). This reflects the presence of the shadow band from ($`\pi `$/2,$`\pi `$/2) to ($`\pi `$,$`\pi `$). The observed distribution of spectral weight is consistent with the most naive picture of a rigid band shift. In this scenario, we would expect the insulating bands with a maximum at ($`\pi `$/2,$`\pi `$/2) to be simply cut by the Fermi function, creating a small Fermi surface pocket around ($`\pi `$/2,$`\pi `$/2). If we extract the contribution only from the vicinity of $`E_\mathrm{F}`$, however, the map looks drastically different. We roughly estimated this low energy excitations by plotting the $`k`$-points where the sharp peak is observed. As shown by the dots in Fig. 3(b), the sharp peaks form a Fermi surface in the form of an “arc” rather than a “pocket”.

We briefly consider here a few theories which appear particularly relevant to this system. In a theory where umklapp scattering is believed to create gaps about ($`\pi `$,0) in both the spin and charge channels, Fermi arcs have been predicted which appear consistent with the present data in Fig. 3(b) . However, the shadow band demonstrated in Fig. 2 seems to escape this formalism. Alternatively, the RVB flux phase , or the more general arguments made by Chakravarty et al, naturally contains shadow bands. In this case the Fermi surface is predicted to appear as four hole pockets centered about ($`\pm \pi `$/2,$`\pm \pi `$/2). However, only the Fermi surface segment associated with the main band contains additional electronic structure (coherent-like peak), while that with the shadow band does not. This may provide a natural starting point for the pocket to evolve into the large Fermi surface without shadow Fermi surface.

Experimentally, the picture that emerges from the present data differs significantly from the electronic evolution recently reported for La2-xSrxCuO4 . Figure 4 demonstrates the comparison of the doping dependence of ARPES spectra between Ca2-xNaxCuO2Cl2 and La2-xSrxCuO4  at ($`\pi `$/2,$`\pi `$/2) and ($`\pi `$,0). In the undoped compounds, both La2CuO4 and Ca2CuO2Cl2, the broad peaks far below $`E_\mathrm{F}`$ represent the insulating gap with the $`d`$-wave-like dispersion . In the case of La2-xSrxCuO4 \[Fig. 4(b)\], the features of the insulator remain at the same energy upon doping while new states are created inside the gap as seen in the ($`\pi `$,0) spectrum . This duality indicates that the chemical potential is fixed upon doping. As summarized in Fig. 4(a), the situation in the doped oxychloride is strikingly different. The evolution of the electronic states in this system can be understood roughly by a shift of the chemical potential to the top of the valence band. Intimately linked with the difference in the electronic evolution, the distribution of the spectral-weight in momentum space contrasts substantially between the two systems. As illustrated in Fig. 4, the predominance of the low energy spectral weight in metallic La1.9Sr0.1CuO4 occurs near ($`\pi `$,0) and (0,$`\pi `$, while it occurs near ($`\pi `$/2,$`\pi `$/2) in the metallic Ca1.9Na0.1CuO2Cl2.

The duality of La2-xSrxCuO4, together with the 1D nature of the ($`\pi `$,0) Fermi surface and the lack of the chemical potential shift, has been linked with stripe formation ; a phase separation into 1D charge and spin stripes. The rigid band shift of the oxychlorides suggests that antiferromagnetic fluctuations, originating from the parent insulator, instead dominate the metallic region for oxychlorides. Why is the oxychloride system so distinct from La2-xSrxCuO4? Structurally, the main difference between the two lies in the apical site (Cl vs. O) and the distortion of the CuO2 planes (undistorted vs. distorted). These slight structural differences might help stabilize various different phases, resulting in enhanced antiferromagnetic correlations in the oxychlorides. In any event, the doped oxychloride system provides a new paradigm in the electronic evolution not previously observed in the cuprates. Surprisingly, however, the appearance of HTS upon doping does not crucially depend on which path leads the parent insulator to a metal.

We would like to thank A. Ino, A. Fujimori, T. Tohyama, S. Maekawa, N. Nagaosa, S. Uchida, and K. Kitazawa for fruitful discussions. This work was partly supported by a Grant-in-Aid for Scientific Research from the MEXT, Japan. ARPES experiments were carried out at SSRL, a national user facility operated by Stanford University on behalf of the U.S. DOE.


text:
# International Lattice Data Grid

## 1 INTRODUCTION

We envisage that the International Lattice Data Grid (ILDG) would be a ‘grid of grids’, integrating the emerging national data grids being developed to manage and exchange lattice QCD data, particularly gauge configurations . Through co-ordination, we could capitalise on these grid developments by sharing expertise and tools, and avoiding duplication of effort. By taking the further step of creating a global grid infrastructure, we would give added international ‘backbone’ to our work.

Since access policies could be determined locally, each collaboration would still be free to decide which, if any, of its configurations to share and with whom. Our proposal is complementary to the Lattice Data Group , which is compiling a reference set of lattice QCD results. Clearly, if combined with ILDG, this would add value to the global grid environment.

The operation of the ILDG would require locally managed, but interoperable, data grids conforming to a set of standards. Data formats, labelling, access policies, management and costs would be handled locally. To achieve this, we will need an international forum to define the standards for interfacing between the grids and to support the federation of databases. We will organise a workshop in Edinburgh on 19/20 December 2002 to begin the technical discussions .

We describe below how an XML schema could be developed for labelling QCD data. Conformance to XML would allow different schemas to be adopted by different data grids, with the additional overhead of developing translation tools to enable each to access the others’ data through a single interface.

## 2 XML SCHEMA

To share resources, either computing or data, one must be able to describe them in a standard language. To share configurations, one must have metadata associated with them, containing all the information about how each configuration was generated and for which physics parameter values. Most collaborations already have some metadata describing their configurations, but, to share configurations automatically, a complete description of the data from the ‘action’ to the ‘algorithm’ is necessary. The UKQCD collaboration is designing an XML schema to describe lattice data, which will be extendable, in principle, to cover any lattice data.

XML is an acronym for eXtensible Markup Language. It is similar to HTML in that it uses ‘$`<>`$’ tags. Unlike HTML, it is for structuring data. XML is license free, platform independent, and well supported. The commercial web services industry uses XML as standard . In fact, XML is a family of technologies. One of these is XML schema. This defines the rules and structure of an XML instance document. It is designed to allow machines to carry out rules made by people.

XML schemas are designed to evolve, to be updated and extended. XSLT provides a method of translating one XML instance document conforming to one schema, into another instance document conforming to different schema. Thus, each collaboration could have its own schema describing the form of its own (meta)data. However, this would involve some duplication of effort, so there would be some benefit if all collaborations in the ILDG used the same schema. UKQCD is designing a schema for our QCDgrid which could form the prototype for such a standard. Our progress can be followed on the UKQCD web pages . Comments are welcome.

## 3 METADATA CATALOGUE

Metadata are mapped into file names by a metadata catalogue. Standard database tools may be used to search the catalogue and, thereby, to access a file through a ‘high-level’, physics-based description of it. We are developing a graphical user interface, initially to access and manage files, and eventually to build and execute programs using them. Again, the catalogue and browser could be made available to all the ILDG collaborations.

We are keen to encourage the inclusion of both state-of-the-art and legacy data on the ILDG. Obviously, the former may have considerable value and access to it might be restricted, at least for an initial period. The ILDG would facilitate the phased release of data under the full control of their originator.

Each collaboration would publish the data it wished to make available to the ILDG on its local node. By browsing the metadata catalogue, someone wishing to use the data could discover what is available and what usage conditions apply. Access control can be achieved using software, more correctly known as middleware, being developed by the HEP experimental community to meet the challenge of LHC computing, e.g. UKQCD is using middleware developed by the EU Datagrid Project , which is itself built on GLOBUS .

## 4 CONCLUSIONS

Compared with a centrally managed respository for gauge configurations, the ILDG concept offers

* more flexible data access control (not ‘all or nothing’);
* the ability to handle different internal data formats easily;
* distributed management, control and operating costs;
* organic growth of both content and functionality; and
* a basis for greater information and resource sharing in future.

We believe that the development of an International Lattice Data Grid would enhance lattice activity worldwide. The collaboration needed to define the ILDG standards could avoid duplication of effort and rapidly promote best practice. If successful as a vehicle for data exchange, ILDG could expand to support open-source software, algorithm development, and eventually even shared computers. A global grid infrastructure would provide a continually growing environment for our science, encouraging wider participation and facilitating physics collaborations. As a symbol of international co-operation, it might even prove useful in making cases for joint machine funding.


text:
# Classical Implications of the Minimal Length Uncertainty Relation

## I Introduction

In this note, we continue our investigation Chang:2001kn ; Benczik:2002tt  of the phenomenological implications of “stringy” commutation relations Kempf:1995su  which embody the minimal length uncertainty relation $`\mathrm{\Delta }x(\mathrm{}/2)(\mathrm{\Delta }p^1+\beta \mathrm{\Delta }p)`$ of perturbative string theory gross . In particular, we study the classical limit of the “deformed” commutation relations

$`[\widehat{x}_i,\widehat{p}_j]`$ $`=`$ $`i\mathrm{}\left\{(1+\beta \widehat{p}^2)\delta _{ij}+\beta ^{}\widehat{p}_i\widehat{p}_j\right\},`$ (1)
$`[\widehat{p}_i,\widehat{p}_j]`$ $`=`$ $`0,`$ (2)
$`[\widehat{x}_i,\widehat{x}_j]`$ $`=`$ $`i\mathrm{}{\displaystyle \frac{(2\beta \beta ^{})+(2\beta +\beta ^{})\beta \widehat{p}^2}{(1+\beta \widehat{p}^2)}}\left(\widehat{p}_i\widehat{x}_j\widehat{p}_j\widehat{x}_i\right),`$ (3)
leading to the following “deformed” Poisson brackets,

$`\{x_i,p_j\}`$ $`=`$ $`(1+\beta p^2)\delta _{ij}+\beta ^{}p_ip_j,`$ (4)
$`\{p_i,p_j\}`$ $`=`$ $`0,`$ (5)
$`\{x_i,x_j\}`$ $`=`$ $`{\displaystyle \frac{(2\beta \beta ^{})+(2\beta +\beta ^{})\beta p^2}{(1+\beta p^2)}}\left(p_ix_jp_jx_i\right).`$ (6)
The classical Poisson bracket is required to possess the same properties as the quantum mechanical commutator, namely, it must be anti-symmetric, bilinear, and satisfy the Leibniz rules and the Jacobi Identity. These requirements allow us to derive the general form of our Poisson bracket for any functions of the coordinates and momenta as Benczik:2002tt

$$\{F,G\}=\left(\frac{F}{x_i}\frac{G}{p_j}\frac{F}{p_i}\frac{G}{x_j}\right)\{x_i,p_j\}+\frac{F}{x_i}\frac{G}{x_j}\{x_i,x_j\},$$
(7)
where repeated indices are summed. Thus, the time evolutions of the coordinates and momenta in our “deformed” classical mechanics are governed by

$`\dot{x}_i`$ $`=`$ $`\{x_i,H\}=\{x_i,p_j\}{\displaystyle \frac{H}{p_j}}+\{x_i,x_j\}{\displaystyle \frac{H}{x_j}},`$ (8)
$`\dot{p}_i`$ $`=`$ $`\{p_i,H\}=\{x_i,p_j\}{\displaystyle \frac{H}{x_j}}.`$ (9)
In Ref. Benczik:2002tt , we analyzed the motion of objects in central force potentials subject to these equations and found that orbits in $`r^2`$ and $`1/r`$ potentials no longer close on themselves when $`\beta `$ and/or $`\beta ^{}`$ are non-zero. This allowed us to place a stringent limit on the value of the minimal length from the observed precession of the perihelion of Mercury,

$$\mathrm{}\sqrt{\beta }<2.3\times 10^{68}\mathrm{m},$$
(10)
which was 33 orders of magnitude below the Planck length.

The natural question to ask next is whether there exist other “deformations” of classical mechanics due to $`\beta `$ and/or $`\beta ^{}`$ that are either 1) observable even for such a small value of $`\mathrm{}\sqrt{\beta }`$, or 2) lead to an even more stringent limit due to their absence. In the following, we will look at the “deformation” of Kepler’s third law. We will find that while the deformation is not observable on the scale of the solar system, it may be observable at galactic scales. In fact, such an effect may have already been seen in the rotation curves of gas and stars in spiral galaxies.

## II Aspects of the Deformed Classical Mechanics

Before we investigate the deformation of Kepler’s third law, it is worthwhile to discuss some peculiarities of this novel classical dynamics.

First, it is obvious that the Poisson bracket $`\{x_i,x_j\}`$ in Eq. (6) is not invariant under translations since it is proportional to $`(\stackrel{}{x}\times \stackrel{}{p})_{ij}`$. For the special case of $`\beta ^{}=2\beta `$, the Poisson bracket $`\{x_i,x_j\}`$ is zero to order $`𝒪(\beta )`$ and we can recover approximate translational invariance Kempf:1995su .

Second, it is not obvious how to extend the fundamental Poisson brackets to multi-particle systems, nor how to define appropriate ‘canonical transformations’ that leave the fundamental Poisson bracket invariant, even in the case of one-particle systems . For instance, consider the motion of 2 bodies of equal mass in 1D. Let us assume that the dynamics will be described by two sets of ‘canonical’ variables, $`(x_1,p_1)`$ and $`(x_2,p_2)`$, which satisfy

$`\{x_i,p_j\}`$ $`=`$ $`\left(1+\beta p_i^2\right)\delta _{ij},`$ (11)
$`\{p_1,p_2\}`$ $`=`$ $`\{x_1,x_2\}=\mathrm{\hspace{0.33em}0}.`$ (12)
Compare this expression to the multi-dimensional case, Eq. (6). The right hand side of the first line is not $`(1+\beta p^2)\delta _{ij}`$, where $`p^2=_ip_i^2`$, but $`(1+\beta p_i^2)\delta _{ij}`$. This allows us to assume $`\{x_1,x_2\}=0`$. If we change the variables naively to those of the center of mass frame,

$$\begin{array}{cc}x_{\mathrm{cm}}=\frac{x_1+x_2}{2},\hfill & p_{\mathrm{cm}}=p_1+p_2,\hfill \\ x_{\mathrm{rel}}=x_1x_2,\hfill & p_{\mathrm{rel}}=\frac{p_1p_2}{2},\hfill \end{array}$$
(13)
then we find that the sets $`(x_{\mathrm{cm}},p_{\mathrm{cm}})`$ and $`(x_{\mathrm{rel}},p_{\mathrm{rel}})`$ are no longer ‘canonical’. Of course, in this deformed dynamics momentum is no longer equal to mass$``$velocity in general, so it is not surprising that the ‘canonical’ momentum of the center of mass, if it exists, is not equal to the sum of momenta of the individual masses.

Third, consider 1D motion with the Hamiltonian given by

$$H=\frac{p^2}{2m}+V(x).$$
(14)
We will not consider any deformations of the Hamiltonian itself though it is conceivable that some sort of modification may be required for the consistency of the theory. The equations of motion read

$`\dot{x}`$ $`=`$ $`\{x,H\}=(1+\beta p^2){\displaystyle \frac{p}{m}},`$ (15)
$`\dot{p}`$ $`=`$ $`\{p,H\}=(1+\beta p^2)\left({\displaystyle \frac{V}{x}}\right).`$ (16)
As mentioned above, the momentum $`p`$ is no longer equal to $`m\dot{x}`$. From these equations, we can derive

$$(1+\beta p^2)(1+3\beta p^2)F=m\ddot{x},$$
(17)
where

$$F\frac{V}{x}.$$
(18)
Thus we obtain a deformation of Newton’s second law $`F=m\ddot{x}`$. Now, notice that if the force $`F`$ is gravitational and proportional to the mass $`m`$, the acceleration $`\ddot{x}`$ is not mass-independent as usual because of the residual $`m`$-dependence through the momentum $`p`$. Therefore, in Eq. (17) the equivalence principle is dynamically violated.

It should be noted that the equivalence principle is already known to be violated in the context of perturbative string theory mende . Fundamental strings, due to their extended nature, are subject to tidal forces and do not follow geodesics. However, in contrast to the current discussion where the systems under consideration are of macroscopic dimensions, the violation discussed in Ref. mende  was microscopic in nature. Whether this is another example of UV/IR correspondence in string theory remains to be seen.

## III Motion in Central Force Potentials

Let us now consider the motion of an object in a general central force potential where the Hamiltonian is given by

$$H=\frac{p^2}{2m}+V(r),r=\sqrt{x_ix_i}.$$
(19)
We will apply the results of this section later to the case when $`V(r)`$ is the gravitational potential and derive a “deformed” version of Kepler’s third law.

As shown in Ref. Benczik:2002tt , the rotational symmetry of the Hamiltonian Eq. (19) leads to the conservation of the ‘deformed’ angular momentum

$$L_{ij}=\frac{x_ip_jx_jp_i}{(1+\beta p^2)},$$
(20)
which in turn implies that the motion will be confined to a 2-dimensional plane. Expressing the position of the object in the plane in polar coordinates $`(r,\varphi )`$, the equations of motion can be cast into the form

$`\dot{r}`$ $`=`$ $`{\displaystyle \frac{1}{m}}[\mathrm{\hspace{0.17em}1}+(\beta +\beta ^{})p^2]p_r,`$ (21)
$`\dot{p}_r`$ $`=`$ $`{\displaystyle \frac{1}{mr}}[\mathrm{\hspace{0.17em}1}+(\beta +\beta ^{})p^2]\left\{(p^2p_r^2)\left(mr{\displaystyle \frac{V}{r}}\right)\left({\displaystyle \frac{1\beta p^2+2\beta p_r^2}{1+\beta p^2}}\right)\right\},`$ (22)
$`\dot{\varphi }`$ $`=`$ $`{\displaystyle \frac{L}{mr^2}}\left\{[\mathrm{\hspace{0.17em}1}+(\beta +\beta ^{})p^2](1+\beta p^2)+[(2\beta \beta ^{})+(2\beta +\beta ^{})\beta p^2]\left(mr{\displaystyle \frac{V}{r}}\right)\right\},`$ (23)
where

$`p_r`$ $``$ $`{\displaystyle \frac{(px)}{r}}=\sqrt{p^2{\displaystyle \frac{L^2(1+\beta p^2)^2}{r^2}}},`$ (24)
$`L^2`$ $``$ $`{\displaystyle \frac{1}{2}}L_{ij}L_{ji}={\displaystyle \frac{p^2r^2(px)^2}{(1+\beta p^2)^2}}.`$ (25)
We restrict our attention to circular orbits and impose the constraints $`\dot{r}=0`$ and $`\dot{p}_r=0`$. These conditions lead to

$$L=\frac{pr}{(1+\beta p^2)},$$
(26)
and

$$\left(mr\frac{V}{r}\right)=p^2\left(\frac{1+\beta p^2}{1\beta p^2}\right).$$
(27)
The radius of the orbit $`r`$, and the magnitude of the momentum $`p=\sqrt{p^2}`$ are now constants of motion. So is the angular velocity $`\omega \dot{\varphi }`$ which can be expressed as

$`\omega `$ $`=`$ $`{\displaystyle \frac{p}{mr}}\left\{[\mathrm{\hspace{0.17em}1}+(\beta +\beta ^{})p^2]+{\displaystyle \frac{(2\beta \beta ^{})+(2\beta +\beta ^{})\beta p^2}{(1+\beta p^2)}}\left(mr{\displaystyle \frac{V}{r}}\right)\right\}`$ (28)
$`=`$ $`{\displaystyle \frac{p}{mr}}\left\{[\mathrm{\hspace{0.17em}1}+(\beta +\beta ^{})p^2]+{\displaystyle \frac{(2\beta \beta ^{})+(2\beta +\beta ^{})\beta p^2}{(1\beta p^2)}}p^2\right\}`$ (29)
$`=`$ $`\left({\displaystyle \frac{p}{mr}}\right){\displaystyle \frac{(1+\beta p^2)^2}{(1\beta p^2)}}.`$ (30)
Solving Eq. (27) for $`p`$ and substituting into Eq. (30) yields

$$m^2v^2=m^2r^2\omega ^2=\left(mr\frac{V}{r}\right)f(\epsilon ),$$
(31)
where

$`f(\epsilon )`$ $``$ $`{\displaystyle \frac{1}{2}}\left[(1+4\epsilon \epsilon ^2)+(1+\epsilon )\sqrt{1+6\epsilon +\epsilon ^2}\right],`$ (32)
$`\epsilon `$ $``$ $`\beta mr{\displaystyle \frac{V}{r}}.`$ (33)
The function $`f(\epsilon )`$ is a monotonically increasing function of $`\epsilon `$ which behaves as

$$f(\epsilon )=\{\begin{array}{cc}1+4\epsilon +𝒪(\epsilon ^2)\hfill & \text{when }\epsilon 1,\hfill \\ 4\epsilon +𝒪(\epsilon ^1)\hfill & \text{when }\epsilon 1.\hfill \end{array}$$
(34)
The graph of $`f(\epsilon )`$ is shown in Fig. 1.

Note that $`\beta ^{}`$ does not appear in our expression for $`v^2`$. This means that the orbital period $`T=2\pi r/v`$ is not affected by $`\beta ^{}`$. In particular, it is independent of whether approximate translational invariance ($`\beta ^{}=2\beta `$) is realized or not. Note also that $`\epsilon `$, and thus the deformation factor $`f(\epsilon )`$, depends not only on $`r`$ but also on $`m`$. As was the case with Eq. (17), this last fact can be interpreted as the dynamical violation of the equivalence principle when our formalism is applied to motion in gravitational fields.

## IV Kepler’s Third Law

For a spherically symmetric mass distribution, the gravitational force on a mass $`m`$ at a distance $`r`$ from the center is given by

$$\frac{V}{r}=m\frac{GM(r)}{r^2},$$
(35)
where $`G`$ is Newton’s gravitational constant, and $`M(r)`$ is the total mass enclosed within the radius $`r`$. Thus, Eq. (31) leads to the relation

$$v^2=\frac{GM(r)}{r}f(\epsilon ),$$
(36)
with

$$\epsilon =\beta m^2\frac{GM(r)}{r}.$$
(37)
In terms of the orbital period $`T=2\pi r/v`$, we obtain

$$\frac{r^3}{T^2}=\frac{GM(r)}{4\pi ^2}f(\epsilon ).$$
(38)
This is the ‘deformed’ version of Kepler’s third law.

Let us estimate how large this deformation can be for the planets in our solar system. Using the data from Ref. AstroData , we calculate the deformation factor $`f(\epsilon )`$ for the planets setting $`\beta `$ equal to the upper bound provided in Eq. (10). The results are shown in Table 1.

We also list the values of

$$\left(\frac{a^3}{T^2}\frac{4\pi ^2}{GM_{}}\right)1,$$
(39)
which show the apparent violation of Kepler’s third law. Note that for elliptical orbits, the radius $`r`$ must be replaced by the semi-major axis $`a`$. (We have not included Mercury or Pluto in our list since the eccentricities of these planets are rather large and the circular orbit approximation is poor.)

As can be seen from Table 1, the expected size of the deviation due to our deformed dynamics is orders of magnitude below the apparent deviation which can be explained by classical Newtonian dynamics as a result of perturbations of the planets on each other. Even in the case of Jupiter, for which the predicted deviation is the largest, the deviation due to the deformation is an order of magnitude smaller than the apparent deviation.

Therefore, we can conclude that as long as the constraint given in Eq. (10) is fulfilled, any deformation of Kepler’s third law due to non-zero $`\beta `$ is completely hidden under apparent violations due to more conventional effects.

## V Rotation Curves of Spiral Galaxies

The deformation of Kepler’s third law will affect the $`r`$-dependence of $`v^2`$. So let us consider whether our deformed dynamics could provide a non-dark matter solution to the halo dark matter problem.

The halo dark matter problem pertains to the rotation curves of interstellar gas in spiral galaxies. ‘Rotation curve’ refers to the plot of $`v^2`$ as a function of the distance $`r`$ from the galactic center. If gravity due to the visible stars is the only force responsible for the centripetal acceleration, then the observed distribution of stars and Kepler’s third law tell us that $`v^2`$ should fall off asymptotically as $`1/r`$ as $`r`$ is increased. However, the actual observed rotation curves approach constant values instead, suggesting the existence of dark matter halos to account for the missing mass halo1 .

In our deformed dynamics, the correction term $`f(\epsilon )`$ does increase the value of $`v^2`$ for a given mass distribution. Unfortunately, the deformation parameter $`\epsilon =\beta m^2GM(r)/r`$ itself falls off as $`1/r`$ with $`r`$, so the correction vanishes for large $`r`$ and cannot explain the plateau of the rotation curves. What our formalism does predict, however, is that the size of the correction will be different for different masses due to the $`m`$-dependence of $`\epsilon `$: larger masses will rotate at higher speeds at the same radius $`r`$, in clear violation of the equivalence principle.

Curiously, it has been recently observed that the rotation curves of stars and interstellar gas do not in general agree when measured separately halo1 ; halo2 . Could our ‘deformed’ dynamics account for this difference? For most galaxies it turns out that the gas rotates at higher speeds than the massive stars. This may have a natural explanation by positing that the gas interacts differently with the dark matter and magnetic fields as compared to the stars. However, for two exceptional galaxies reported in Ref. halo2 , the stars are rotating at much higher speeds than the gas. For the galaxy NGC7331, the speeds of the stars and gas at distance $`R_e/4`$ from the center of the galaxy ($`R_e`$ is the radius which encloses half of the visible light from the galaxy) are reported to be

$`v_{\mathrm{gas}}`$ $`=`$ $`79\pm 20\mathrm{km}/\mathrm{s},`$ (40)
$`v_{\mathrm{star}}`$ $`=`$ $`186\pm 10\mathrm{km}/\mathrm{s}.`$ (41)
If we interpret this galaxy as a case where the mechanisms that would normally accelerate the gas over the stars are absent and the effect of ‘deformed’ dynamics is manifest, we can ask what value of $`\beta `$ would account for this difference. Since

$`v^2`$ $`=`$ $`{\displaystyle \frac{GM(r)}{r}}\left\{1+4\beta m^2{\displaystyle \frac{GM(r)}{r}}+\mathrm{}\right\}`$ (42)
$`=`$ $`{\displaystyle \frac{GM(r)}{r}}\left(1+4\beta m^2v^2+\mathrm{}\right),`$ (43)
then

$$\frac{v_{\mathrm{star}}^2}{v_{\mathrm{gas}}^2}1+4\beta m_{\mathrm{star}}^2v_{\mathrm{star}}^2,$$
(44)
where we have neglected the correction term for the gas. If we assume the solar mass

$$M_{}=2\times 10^{30}\mathrm{kg},$$
(45)
as a typical mass of a star in this galaxy, we find

$$\mathrm{}\sqrt{\beta }10^{69}\mathrm{m},$$
(46)
which is well within the limit given by Eq. (10) and thus compatible with observations on the scale of the solar system.

## VI Conclusion and Discussion

In this note, we have investigated possible observable consequences of the “deformed” Kepler’s third law as implied by the classical limit of the “stringy” commutation relations, Eq. (3), which were in turn based on the minimal length uncertainty relation of perturbative string theory gross .

Due to the stringent limit on $`\beta `$, Eq. (10), obtained from the precession of the perihelion of Mercury in a previous publication Benczik:2002tt , the predicted sizes of the deviations for the planets are too small to be observed. However, the effect of the deformation can be amplified at galactic scales due to the large momenta involved and lead to significant changes in the rotation curves of stars in spiral galaxies.

It would be interesting to contrast our deformed classical mechanics with another modification of Newton’s theory, MOND (Modified Newtonian Dynamics) MOND , which apparently succeeds in predicting a flat rotation curve by modifying Newton’s second law to

$$Fm\frac{|\ddot{x}|}{a_0}\ddot{x},$$
(47)
at very small accelerations ($`\ddot{x}a_010^{10}\mathrm{m}/\mathrm{s}^2`$). Our deformed dynamics predicts different rotation curves for different masses, so it is apparently quite different from MOND which does not violate the equivalence principle. Nevertheless, it is interesting to ask whether there exists a deformation of the usual Hamiltonian dynamics which incorporates some kind of UV/IR relation and is in the same universality class as MOND. We hope to address this question in the future.

###### Acknowledgements.

We would like to thank Vishnu Jejjala, Achim Kempf, Yukinori Nagatani, Makoto Sakamoto, Hidenori Sonoda, and Joseph Slawny for helpful discussions, and Alan Chamberlin for his help with the planetary data. This research is supported in part by a grant from the U.S. Department of Energy, DE–FG05–92ER40709, Task A.


text:
# Gauge Invariant Evaluation of Nuclear Polarization with Collective Model

## Abstract

The nuclear-polarization (NP) energies with the collective model commonly employed in the NP calculations for hydrogenlike heavy ions are found to have serious gauge violations when the ladder and cross diagrams only are taken into account. Using the equivalence of charge-current density with a schematic microscopic model, the NP energy shifts with the collective model are gauge invariantly evaluated for the $`1s_{1/2}`$ states in $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+ and $`{}_{92}{}^{}{}_{}{}^{238}`$U91+.

High-precision Lamb-shift measurement on high-Z hydrogenlike atoms BE95  spurred a renewed interest in the quantum electrodynamical (QED) calculation of electronic atoms. Comparison of theoretical results with experimental data allows sensitive tests of QED in strong electromagnetic fields SA90 ; MO98 . In this context, the study of the nuclear-polarization (NP) effect becomes important because the NP effect, as a non-QED effect which depends on the model used to describe the nuclear dynamics, sets a limit to any high-precision test of QED.

A relativistic field-theoretical treatment of NP calculation was presented by Plunien et al. PL89 ; PL95  utilizing the concept of effective photon propagators with nuclear-polarization insertions. In these studies, only the Coulomb interaction was considered based on the argument that the relative magnitude of transverse interaction is of the order of $`(v/c)^2`$ and the velocity $`v`$ associated with nuclear dynamics is mainly nonrelativistic.

The effect of the transverse interaction was studied in the Feynman gauge by Yamanaka et al. YA01  with the same collective model used in PL89 ; PL95 ; NE96  for nuclear excitations. They found that the transverse contribution is several times larger than the Coulomb contribution in heavy electronic atoms before the contributions of the positive and negative energy states cancel. However, due to the nearly complete cancellation between them, the transverse effects become small and the net effect is destructive to the Coulomb contribution in both $`1s_{1/2}`$ states of $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+ and $`{}_{92}{}^{}{}_{}{}^{238}`$U91+. As a result, the total NP energy almost vanishes in $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+.

Recently, the NP effects for hydrogenlike and muonic $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+ were calculated in both the Feynman and Coulomb gauges, using a microscopic random phase approximation (RPA) to describe nuclear excitations HHT02 ; HHT02a . It was found that, in the hydrogenlike atom, the NP effects due to the ladder and cross diagrams have serious gauge dependence and inclusion of the seagull diagram is indispensable to restore the gauge invariance HHT02 . In contrast, the magnitude of the seagull collection is a few percent effect in the muonic atom, although it improves the gauge invariance HHT02a .

In the present paper, we report that the nuclear collective model employed for hydrogenlike ions in NE96 ; YA01 ; PL89 ; PL95  also leads to a large violation of gauge invariance as far as the ladder and cross diagrams only are considered. Then it is shown, based on the equivalence of the transition density of the collective model and a microscopic nuclear model with a schematic interaction between nucleons, that the seagull corrections should also be calculated with the collective model in order to obtain gauge invariant NP results. The resulting gauge invariant NP energy shifts are given for the $`1s_{1/2}`$ states in $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+ and $`{}_{92}{}^{}{}_{}{}^{238}`$U91+.

For spherical nuclei, the Hamiltonian of the small amplitude vibration with multipolarity $`L`$ is written as

$`H_L={\displaystyle \frac{1}{2}}({\displaystyle \frac{1}{D_L}}{\displaystyle \underset{M}{}}\widehat{\pi }_{LM}^{}\widehat{\pi }_{LM}+C_L{\displaystyle \underset{M}{}}\widehat{\alpha }_{LM}^{}\widehat{\alpha }_{LM}),`$ (1)
where $`\widehat{\pi }_{LM}`$ are the canonically conjugate momenta to the collective coordinates $`\widehat{\alpha }_{LM}`$. The lowest vibrational modes are expected to have density variations with no radial nodes, which may be referred to as shape oscillations. The corresponding charge density operator with the multipolarity $`L`$ is written as

$`\widehat{\rho }_L(t,𝒓)`$ $`=\rho _L(r){\displaystyle \underset{M}{}}Y_{LM}^{}\widehat{\alpha }_{LM}(t)`$ (2)
to the lowest order of $`\widehat{\alpha }_{LM}^{}(t)`$.

The liquid drop model of Bohr (BM)Bohr52  is a simple model of such a shape oscillation obtained by considering deformation of the nuclear radius parameter while leaving the surface diffuseness independent of angle:

$`R(\mathrm{\Omega })=R_0\left[1+{\displaystyle \underset{LM}{}}\alpha _{LM}Y_{LM}^{}(\mathrm{\Omega })\right],`$ (3)
where $`R_0`$ is the nuclear radius parameter of the ground state. The radial charge density of BM is given by

$`\rho _L(r)=R_0{\displaystyle \frac{d}{dr}}\varrho _0(𝒓),`$ (4)
where $`\varrho _0(𝒓)`$ is a charge distribution with spherical symmetry.

If we assume that under distortion, an element of mass moves from $`𝒓_0`$ to $`𝒓`$ without alteration of the volume it occupies, i.e., the nucleus is composed of an inhomogeneous incompressible fluid, a harmonic vibration of an originally spherical surface $`r=r_0`$ in the nucleus is given by

$`r(\mathrm{\Omega })=r_0\left[1+{\displaystyle \underset{LM}{}}\left({\displaystyle \frac{r_0}{R_0}}\right)^{L2}\alpha _{LM}Y_{LM}^{}(\mathrm{\Omega })\right].`$ (5)
For this model we obtain

$`\rho _L(r)={\displaystyle \frac{1}{R_0^{L2}}}r^{L1}{\displaystyle \frac{d}{dr}}\varrho _0(𝒓).`$ (6)
This version will be hereafter referred to as the Tassie Model (TM) Tassie . In Eqs.(4) and (6), $`\varrho _0(𝒓)`$ is usually taken to be equal to the ground-state charge distribution.

In either case, the motion of nuclear matter is assumed to be incompressible and irrotational, hence the velocity field $`𝒗(t,𝒓)`$ is given by a velocity potential as $`𝒗(t,𝒓)=\mathbf{}\mathrm{\Phi }(t,𝒓)`$. This implies the nuclear current defined by $`𝑱(𝒓)=\varrho _0(r)𝒗(𝒓)`$ yields the transition multipole density of current operator

$`\widehat{𝑱}_L(t,𝒓)`$ $`=J_{LL1}(r){\displaystyle \underset{M}{}}𝒀_{LL1M}^{}\widehat{\alpha }_{LM}(t).`$ (7)
Note that the $`J_{LL+1}(r)`$ part does not appear in the transition density of current operator given by (7).

Therefore, in this kind of collective model, the continuity equation of charge gives

$`[i\mathrm{\Delta }E_L\rho _L(r)+\sqrt{{\displaystyle \frac{L}{2L+1}}}({\displaystyle \frac{d}{dr}}{\displaystyle \frac{L1}{r}})J_{LL1}(r)]=0,`$ (8)
where $`\mathrm{\Delta }E_L`$ is the excitation energy of the surface oscillation. Hence the transition density of current is given by

$`J_{LL1}(r)=i\mathrm{\Delta }E_L\sqrt{{\displaystyle \frac{2L+1}{L}}}r^{L1}{\displaystyle _r^{\mathrm{}}}x^{1L}\rho _L(x)𝑑x`$ (9)
in terms of the transition density of charge. If we assume the uniform charge distribution $`\varrho _0(r)=\varrho _0\mathrm{\Theta }(R_0r)`$, we obtain, for both BM and TM,

$`\rho _L(r)`$ $`=J_fr^LY_LJ_i{\displaystyle \frac{1}{R_0^{L+2}}}\delta (R_0r),`$ (10)
$`J_{LL1}(r)`$ $`=J_fr^LY_LJ_ii\mathrm{\Delta }E_L\sqrt{{\displaystyle \frac{2L+1}{L}}}{\displaystyle \frac{r^{L1}}{R_0^{2L+1}}}\mathrm{\Theta }(R_0r).`$ (11)
The transition densities given by (10) and (11) have been employed in the previous NP calculations for $`L1`$ PL89 ; PL95 ; NE96 ; YA01 . It should be mentioned that, although the surface oscillation applies to the case of the multipolarity $`L2`$, Eqs. (10) and (11) with $`L=1`$ give the transition densities of the giant dipole resonance given by the Goldhaber-Teller model describing the relative motion of neutrons and protons DW66 . For the monopole vibration, it is also possible to construct corresponding charge and current densities YA01 ; PL89 .

In general, the charge conservation relation between the charge and current densities is necessary but not sufficient for the gauge invariance of the NP calculation. Unfortunately, it is practically impossible to construct a model that incorporates gauge invariance explicitly in terms of the collective variables of the model. However, it is possible to evaluate the NP energy gauge invariantly with the above collective model as is shown below.

The NP calculations with the collective model assume that a single giant resonance with spin multipolarity $`L`$ saturates the energy-weighted $`B(EL)`$ strength for each isospin. In this respect, let us recall the fact that the transition densities of charge to the sum-rule saturated levels are given in terms of the ground-state charge density DF73 . This can be seen as follows. For a pair of single-particle operators $`g(𝒓)=g(r)Y_{LM}(\mathrm{\Omega })`$ and $`f(𝒓)=f(r)Y_{LM}(\mathrm{\Omega })`$, the energy-weighted sum rule can be generalized to

$`{\displaystyle \frac{1}{2J_i+1}}{\displaystyle \underset{n}{}}(E_nE_i)[J_ng(r)Y_LJ_i^{}J_nf(r)Y_LJ_i]`$
$`={\displaystyle \frac{2L+1}{4\pi }}{\displaystyle \frac{h^2}{2M}}{\displaystyle r^2𝑑r\varrho _0(r)[g^{}(r)f^{}(r)+\frac{L(L+1)}{r^2}g(r)f(r)]},`$ (12)
where $`\varrho _0(r)`$ is the charge distribution of the ground state normalized as $`r^2𝑑r\varrho _0(r)=Z`$ BM75 . When a single excited state $`|J_fM_f`$ saturates the $`B(EL)`$ strength, $`|J_fM_fr^LY_{LM}|J_iM_i`$, the transition density of charge to this state is derived from the sum-rule relation (12) model independently and given by

$`\varrho _{fi}(r)={\displaystyle \frac{1}{2L+1}}{\displaystyle \frac{J_fr^LY_LJ_i}{J_i|r^{2L2}|J_i}}r^{L1}{\displaystyle \frac{d}{dr}}\varrho _0(r).`$ (13)
If the charge distribution of the ground state is assumed to be a uniform distribution with a radius $`R_0`$, this becomes

$`\varrho _{fi}(r)=J_fr^LY_LJ_i{\displaystyle \frac{1}{R_0^{L+2}}}\delta (rR_0),`$ (14)
which is equal to the matrix element of the charge density operator of the collective model given by (10).

On the other hand, it is well known that the schematic RPA with a separable interaction

$`V_S(𝒓_i,𝒓_j)=\kappa _L{\displaystyle \underset{M}{}}r_i^LY_{LM}(\mathrm{\Omega }_i)r_j^LY_{LM}^{}(\mathrm{\Omega }_j).`$ (15)
for particle-hole excitations $`|mi^1`$ with a degenerate particle-hole excitation energy $`ϵ`$ gives a collective state $`|LM`$, which exhausts the energy-weighted sum rule for the single particle operator $`r^LY_{LM}`$:

$`\mathrm{\Delta }E_L|LM|r^LY_{LM}|0|^2=ϵ{\displaystyle \underset{mi}{}}|m|r^LY_{LM}|i|^2.`$ (16)
where $`\mathrm{\Delta }E_L`$ is the excitation energy of $`|LM`$ RS80  . If the ground state is assumed to be a filled major shell of the harmonic oscillator potential:

$`H_{HO}={\displaystyle \frac{1}{2M_N}}𝒑^2+{\displaystyle \frac{M_N\omega ^2}{2}}𝒓^2,`$ (17)
the particle-hole excitation energy $`ϵ`$ is taken to be $`1\mathrm{}\omega `$ for $`1^{}`$ and $`2\mathrm{}\omega `$ for $`0^+`$ and $`2^+`$. The corresponding collective states exhaust the energy-weighted sum rules, because the transition strengths vanish outside these p-h excitation spaces. Therefore, the transition densities of charge to the collective states of this fictitious nucleus are given by (13). When the ground-state charge density is approximated by a uniform charge density, the transition density of charge becomes identical with that of the collective model employed in NP calculations for hydrogenlike atoms. However, the gauge invariant electromagnetic interaction of this schematic microscopic model is given by the minimal substitution $`𝒑_i𝒑_ie_i𝑨`$ to the Hamiltonian $`H=H_{HO}+V_S`$. Hence the lowest-order contributions to NP with this model are given by the three Feynman diagrams in Fig. 1, where two photons are exchanged between a bound electron and a nucleus. The nuclear vertices are understood to have no diagonal matrix elements for the ladder and cross diagrams, and no nuclear intermediate states for the seagull diagram. It is well known that the NP results with this model is gauge invariant provided these three diagrams are taken into account. Although $`J_{LL+1}(r)`$ current density appears in this model, $`J_{LL1}(r)`$ dominates in the transition to the collective state.

Thus we can conclude that the gauge invariance of the collective model is also guaranteed with the charge-current density satisfying the continuity equation (8), provided the contributions from the three diagrams are taken into account. It should be noted that the seagull contribution is given in terms of the ground-state charge distribution and does not depend on the details of the model for nuclear excitations. In the actual NP calculations NE96 ; PL89 ; PL95 ; YA01  with the collective model, the assumption that each nuclear intermediate state saturates the sum-rule is not strictly obeyed, because the observed nuclear data are used for the low-lying states. However, since the gauge violation is serious only in the dipole giant resonance, this does not invalidate our arguments as is confirmed by the numerical results in the following.

The formulas to calculate the NP energy shifts due to the three diagrams of Fig. 1 were given in HHT02  for arbitrary nuclear models. In the present NP calculations of the $`1s_{1/2}`$ states in hydrogenlike $`{}_{82}{}^{}{}_{}{}^{208}`$Pb and $`{}_{92}{}^{}{}_{}{}^{238}`$U, the parameters of the collective model are the same as those given in Refs. NE96 ; YA01 . The same low-lying states and giant resonances are taken into account. In addition, the contributions from the $`4^{}`$ and $`5^{}`$ giant resonances are also calculated in order to see the effects of higher multipoles neglected previously. The $`B(EL)`$ values are adjusted to the observed values for low-lying states and the $`B(EL)`$ are estimated through the energy-weighted sum rule for giant resonances.

Tables I and II show the results for the $`1s_{1/2}`$ states in $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+ and $`{}_{92}{}^{}{}_{}{}^{238}`$U91+, where the sum of the contributions from the three diagrams of Fig.1 is given for each multipole. The second and the third columns are the results including the transverse effects in the Feynman and Coulomb gauges, respectively. The values in the parentheses are the contributions from the seagull diagram. The NP energy shifts due to the ladder and crossed diagrams only are obtained by subtraction of the seagull contributions given in the parentheses. The fourth column gives the results of the present Coulomb nuclear polarization (CNP). The last two columns are the results of the previous calculations.

The results with the collective model, as with the microscopic RPA model HHT02 , also lead to large violations of gauge invariance if ladder and crossed diagram contributions only are considered. The seagull corrections are considerable in the $`1^{}`$ contributions for both of $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+ and $`{}_{92}{}^{}{}_{}{}^{238}`$U91+. Note that, in the limit of point nucleus, which is not unrealistic even for heavy hydrogenlike ions, the seagull collection occurs only in the dipole mode which involves the current density $`J_{10}(r)`$.

In $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+, the contributions from low-lying states are about 10% of the total results and the NP energy shift is mainly determined by the giant resonance contributions. The most dominant contribution comes from the giant dipole resonance, where a large violation of gauge invariance occurs if the seagull contributions in the parentheses are neglected: $`22`$ meV becomes $`+20`$ meV and $`14`$ meV in Feynman and Coulomb gauges, respectively. The column 5 gives the previous results in the Feynman gauge without seagull contributions. The differences between the two results in the Feynman gauge without seagull contribution come from the accuracy of numerical integration over the continuum threshold region of electron intermediate states and from the differences of the electron wave functions: here we have used wave functions in a finite charge distribution, while YA01  employs point Coulomb solutions.

In $`{}_{92}{}^{}{}_{}{}^{238}`$U91+, the dominant contribution comes from the lowest excited states $`2^+`$ with a large $`B(E2)`$ value. Since the transition density of current in the present model given by (11) is proportional to the excitation energy, the transverse contribution of the lowest $`2^+`$ is negligible due to its exceptionally small excitation energy $`\mathrm{\Delta }E_2=44.9`$ keV. Apart from this large Coulomb contribution, the contributions from other states show similar tendencies as in $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+. Namely, the contributions from other low-lying states are small compared with the giant resonance contributions, and a large gauge violation occurs in the giant dipole resonance when the seagull contribution is omitted.

To summarize, the transverse effects with the collective model are estimated gauge invariantly by inclusion of the seagull contribution. The gauge invariance is satisfied to a few percent levels in both $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+ and $`{}_{92}{}^{}{}_{}{}^{238}`$U91+ for each of the multipoles separately. Without the seagull correction, the Feynman gauge in particular does not give reliable predictions of NP, although numerical calculation in this gauge is easier than in the Coulomb gauge. Hence the conclusion of YA01  on the transverse effects is no longer tenable. The NP energy shifts are $`35.0(35.4)`$ meV in $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+ and $`205(203)`$ meV in $`{}_{92}{}^{}{}_{}{}^{238}`$U91+ for Coulomb (Feynman) gauge. The net transverse effect is about $`1415\%`$ of the Coulomb energy shift of $`30.7`$ meV in $`{}_{82}{}^{}{}_{}{}^{208}`$Pb81+. This is similar to the conclusion of the microscopic model HHT02 , and should be compared with the transverse effect of the $`1s_{1/2}`$ state in muonic $`{}_{82}{}^{}{}_{}{}^{208}`$Pb, which is about 6 % of the Coulomb contribution HHT02a . The agreement between the two models provides stability of the prediction of the NP effects with respect to the choice of the nuclear models. The percentage of the transverse effect in the total shift in $`{}_{92}{}^{}{}_{}{}^{238}`$U91+ is reduced to about $`6`$% of the Coulomb effect due to the dominant Coulomb contribution from the lowest $`2^+`$ state.

The authors wish to acknowledge Prof. Y. Tanaka for generous support and useful discussions during our research on NP effects. They appreciate Drs. N. Yamanaka and A. Ichimura for collaboration on the NP effects with the collective model, which motivated the present work.