text:
# Magnetic Anisotropy in the Molecular Complex V15

## I Introduction

Magnetic molecules have attracted significant attention in recent years as systems where macroscopic quantum phenomena are displayed, such as relaxation by tunneling of magnetization through a potential well -. These molecules form crystals with very weak intermolecular interactions, thus experiments reflect directly the properties of individual molecules, which can be described by magnetic Hamiltonians. The most studied molecules are a cluster of manganese ions known as Mn12 , and clusters with iron ions, Fen .

Mn12 and Fe8 possess strong easy axis anisotropy and their states can be characterized to a very good approximation by $`<S^z>`$, the expectation value of the total spin along the easy axis, $`z`$. Their ground state has a total spin $`S=10`$, and to a good approximation they are described by a single-spin model, where only the $`S=10`$ multiplet states are considered . Restriction to the $`S=10`$ multiplet reduces drastically the Hilbert space of the problem. Iron atoms also form ferric wheel clusters with predominantly antiferromagnetic interactions. In these the total spin of the ground state is $`S=0`$-. Spin dynamics of Fe6 and Fe8 have recently been calculated by Honecker et al.  using exact diagonalization.

Another molecule of the same class is a cluster with $`15`$ Vanadium ions, known as V15 , each of which has $`S=\frac{1}{2}`$. This is similar to some of the Fen systems in that the principal interactions are antiferromagnetic. This system has also been investigated using effective spin models in which sets of spins are replaced by one effective spin. However, effective three-spin models which contain only isotropic terms - can not explain the broadening of the spin transitions as a function of magnetic field in equilibrium magnetization measurements, as was pointed out by Chiorescu et al. . Miyashita and Nagaosa  considered a three-spin model for V15 which included a term, $`_{ij}\alpha _{ij}S_i^zS_j^x`$. This term is responsible for the quantum tunneling of magnetization states and for the smoothing out of transitions between magnetization states in applied magnetic fields. However, a more microscopic model is needed if a quantitative comparison with data is to be made. In particular, the tunnel splittings, which are important for magnetic relaxation, are very sensitive to Hilbert space truncation and the full Hilbert space of the molecule has to be taken into account . Here we solve for the lowest energy states using degenerate perturbation theory in the full Hilbert space, which has $`2^{15}=32,768`$ states.

Apart from the opportunity to study the magnetic properties of finite systems which these systems provide, Leuenberger and Loss  have proposed the use of Mn12 as a qubit, a quantum bit, in a quantum computer. In this type of application it is important that prepared spin states do not decohere on the time scale of the operations on the system necessary to read, write and carry out the manipulations necessary for a computation. Consequently it is important to understand the mechanisms for magnetic relaxation in these molecules. These mechanisms can be probed by the application of time dependent magnetic fields .

The magnetic hysteresis studies for V15 show so-called “butterfly” curves which are different from the ones of Mn12 and Fe8 . This has been explained in terms of spin-phonon interaction using the Landau-Zener model for transitions between two states . However, the low number of phonons at very low temperatures points to a different mechanism for magnetization relaxation, intrinsic to the molecule. Here we examine the role of magnetic anisotropy described by the Dzyaloshinskii-Moriya (DM) interactions as a relaxation mechanism. This interaction arises from spin-orbit coupling and was introduced by Dzyaloshinskii phenomenologically by using symmetry considerations . Moriya derived it microscopically by extending the Anderson theory of superexchange to include spin-orbit coupling . The DM term arises in low symmetry crystals. The fact that finite systems have surfaces where symmetry is reduced ensures that the DM term is an important part of any microscopic description of magnetic molecules . There are other anisotropy sources such as dipolar and hyperfine fields which have very small magnitudes, about $`1`$ mT and $`40`$ mT respectively . Magnetic anisotropy has also been modeled in Ising Hamiltonians by $`S^+`$ and $`S^{}`$ terms to reproduce the experimental data . Single site anisotropy terms do not contribute in V15 since $`\sigma _i^2=1`$, $`i=x,y,z`$ for $`S=\frac{1}{2}`$. Therefore the DM interaction is the leading order anisotropy term for this system.

The model Hamiltonian, which includes a Heisenberg exchange interaction and a DM anisotropy term, is solved here with degenerate perturbation theory . It is found that inclusion of the DM terms quantitatively explains the widening of the spin transition in the equilibrium magnetization curve. The DM terms are determined quantitatively by comparing with experimental data of Chiorescu et al. , and the resulting splitting of the low energy groups of states is in agreement with experimental estimates, $`10^2`$ K. Furthermore, it is shown with a model calculation that the non-equilibrium magnetization curve is intrinsically hysteretic due to the presence of the DM interaction. This interaction leads to rapid oscillations wth frequencies given by its magnitude. The form of the hysteresis curve depends on when during these oscillations the magnetic field is reversed.

The paper is organized as follows. In section II we discuss the magnetic description of the V15 molecule and in section III we introduce the model including the DM interaction determined by symmetry. In section IV we present the properties of the model and compare with experimental data. Finally in section V we give our conclusions.

## II V15

The chemical structure of V15 is given by the formula K6 \[V$`{}_{}{}^{IV}{}_{15}{}^{}`$As6O42(H2O)\] $``$ 8H2O . There are fifteen VIV ions with spin $`S=\frac{1}{2}`$ placed in a quasi-spherical layered structure (figure 1). These spins sit on two hexagons with a triangle sandwiched between them and all interactions are antiferromagnetic in the simplest model . The spins in each hexagon are paired via a strong coupling $`J800`$ K . These pairs are connected with a nearest neighbor interaction $`J^{}150`$ K. There is also a diagonal interaction $`J^{\prime \prime }300`$ K. The triangle spins do not interact directly with each other. They interact with spins from the top and bottom hexagons, and the interactions are $`J_1J^{}`$ and $`J_2J^{\prime \prime }`$. The molecule has trigonal symmetry and the space group is $`R_{\overline{3}c}`$. The unit cell includes two molecules but intermolecular interactions are negligible.

Since all interactions are antiferromagnetic the total spin of the ground state is $`S=\frac{1}{2}`$. This is in contrast with Mn12 and Fe8, where $`S=10`$. In the absence of anisotropy the ground state is four-fold degenerate, and it is made up of two doublets. The lowest lying excited states form a quartet which is separated from the ground state doublets by $`3.7`$ K. The splitting between the lowest energy multiplets from the rest of the spectrum has been found from EPR spectra to be $`500`$ K . Magnetic anisotropy leads to tunnel splittings between the two doublets, which are $`10^2`$ K , in contrast with the very small values for the other two molecules ($`10^710^{11}`$ K). There is no macroscopic quantum tunneling for V15 according to the Kramers theorem , since the system has an odd number of $`S=\frac{1}{2}`$ spins and the total spin is half-integer , giving a doubly degenerate ground state.

From equilibrium magnetization measurements, it is clear that the ground state has $`S\frac{1}{2}`$ up to a critical value of the field, where a transition to $`S\frac{3}{2}`$ is detected . This is true for temperatures smaller than $`0.9`$ K. The broadening of the spin transitions can not be attributed to temperature alone, and various factors like dipolar or nuclear hyperfine field distributions are too small, about $`1`$ mT and $`40`$ mT respectively . This broadening can be explained by anisotropic DM interactions, which were also found to be important for the description of neutron scattering data and EPR measurements in Mn12 . More recently, analysis of neutron scattering data on V15 by Chaboussant et al.  has shown the need to include DM interactions. Garanin and Chudnovsky  have suggested that the width of EPR lines could also be explained by dislocations in the crystal lattice. Here we assume that there are no defects in the crystals so that data reflects the properties of the V15 clusters.

## III Model

A microscopic treatment of the problem of the origin of magnetic anisotropy is considered using the full Hilbert space with $`2^{15}=32,768`$ states. This approach based on degenerate perturbation theory makes it possible to include the effect of spin-orbit coupling through the DM interaction. Initially perturbation theory is applied in the absence of magnetic anisotropy. The unperturbed Hamiltonian, $`H_0`$, takes into account the spin singlet-triplet correlations present due to the strength of the bond between alternating pairs of spins in the hexagons being more than twice as large as any other coupling. The ground state of $`H_0`$ is eight-fold degenerate due to uncoupled spins in the triangle (figure 1). Then the perturbation includes the remaining terms for bonds on which exchange constants are smaller. An effective Hamiltonian, $`H^{eff}`$, is generated for the degenerate subspace. Since singlet-triplet correlations of the prevalent $`J`$ bond are built into $`H_0`$, the series expansions generated by perturbation theory are absolutely convergent. Diagonalizing $`H^{eff}`$ gives the eight lowest energy states of the system which are separated from the next set of states roughly by the energy splitting between singlet and triplet states in $`H_0`$. The Hamiltonian is of the form

$$H=H_0+\lambda H_1$$
(1)
where $`\lambda `$ is the perturbation parameter with

$$H_0=J(\stackrel{}{S}_1\stackrel{}{S}_2+\stackrel{}{S}_3\stackrel{}{S}_4+\stackrel{}{S}_5\stackrel{}{S}_6+\stackrel{}{S}_{10}\stackrel{}{S}_{11}+\stackrel{}{S}_{12}\stackrel{}{S}_{13}+\stackrel{}{S}_{14}\stackrel{}{S}_{15})$$
(2)
and

$`H_1=`$ $`J^{}(\stackrel{}{S}_2\stackrel{}{S}_3+\stackrel{}{S}_4\stackrel{}{S}_5+\stackrel{}{S}_1\stackrel{}{S}_6+\stackrel{}{S}_{11}\stackrel{}{S}_{12}+\stackrel{}{S}_{13}\stackrel{}{S}_{14}+\stackrel{}{S}_{10}\stackrel{}{S}_{15})+`$ (3)
$`J^{\prime \prime }(\stackrel{}{S}_1\stackrel{}{S}_3+\stackrel{}{S}_3\stackrel{}{S}_5+\stackrel{}{S}_1\stackrel{}{S}_5+\stackrel{}{S}_{10}\stackrel{}{S}_{12}+\stackrel{}{S}_{12}\stackrel{}{S}_{14}+\stackrel{}{S}_{10}\stackrel{}{S}_{14})+`$ (4)
$`J_1(\stackrel{}{S}_2\stackrel{}{S}_7+\stackrel{}{S}_7\stackrel{}{S}_{11}+\stackrel{}{S}_4\stackrel{}{S}_8+\stackrel{}{S}_8\stackrel{}{S}_{15}+\stackrel{}{S}_6\stackrel{}{S}_9+\stackrel{}{S}_9\stackrel{}{S}_{13})+`$ (5)
$`J_2(\stackrel{}{S}_1\stackrel{}{S}_7+\stackrel{}{S}_7\stackrel{}{S}_{10}+\stackrel{}{S}_3\stackrel{}{S}_8+\stackrel{}{S}_8\stackrel{}{S}_{14}+\stackrel{}{S}_5\stackrel{}{S}_9+\stackrel{}{S}_9\stackrel{}{S}_{12})`$ (6)
The ground state of the unperturbed Hamiltonian is composed of a singlet at each of the six hexagon pairs. When $`\lambda 0`$ the perturbation mixes in states where one or more of the pairs are excited and generate the elements of $`H^{eff}`$. Since $`J`$ is a very strong bond, this lowest multiplet of $`8`$ states will be well-separated in energy from the rest of the spectrum, and the low temperature behavior is well described by the eigenfunctions of $`H^{eff}`$.

In equation (1) the exchange constants given in  are varied when $`\lambda =1`$ so that the gap from the doublets to the quartet agrees with the experimental value of $`3.7`$ K. The following parameters are kept:

$$J=1\text{ , }J_1=J^{}=\frac{225}{800}J\text{ , }J_2=J^{\prime \prime }=\frac{350}{800}J$$
(7)
The unit for $`J`$ is $`800`$ K and the energy gap to the excited state with these values is equal to $`\mathrm{\Delta }_1=3.61`$ K. Thus, the low energy diagram of V15 is as shown in figure 2. $`\mathrm{\Delta }_010^2`$ K is the tunnel splitting between the two doublets. This splitting is not described by the present Hamiltonian. A test for the success of perturbation theory is the calculation of the expectation values of the square of the total spin $`\stackrel{}{S}^2`$ and its projection on the $`z`$ axis, $`S^z`$, for the calculated wavefunctions. These expectation values are half-integer for $`S^z`$ and integer for $`\stackrel{}{S}^2`$ at $`\lambda =1`$ with a precision equal to the one used to generate the perturbation expansions (double precision-15 digits). This confirms the success of perturbation theory.

### A Magnetic Anisotropy

The anisotropic interaction between two spins $`\stackrel{}{S}_i`$ and $`\stackrel{}{S}_j`$ due to spin-orbit coupling is :

$$\stackrel{}{D}_{i,j}\stackrel{}{S}_i\times \stackrel{}{S}_j+\stackrel{}{S}_i𝚪_{𝐢,𝐣}\stackrel{}{S}_j$$
(8)
The first term is the DM interaction and it is first order in spin-orbit coupling and antisymmetric. We refer to $`\stackrel{}{D}_{i,j}`$ as the DM vector. The second term is second order in spin-orbit coupling and symmetric, with $`𝚪_{𝐢,𝐣}`$ a second-rank tensor. Thus the DM vector provides the lowest order anisotropy term.

As a minimal model, the DM interaction will be considered on the $`J`$ bonds only, the strongest exchange constants in the system. The DM interaction is:

$$\stackrel{}{D}_{i,j}\stackrel{}{S}_i\times \stackrel{}{S}_j=D_{i,j}^x(S_i^yS_j^zS_i^zS_j^y)+D_{i,j}^y(S_i^zS_j^xS_i^xS_j^z)+D_{i,j}^z(S_i^xS_j^yS_i^yS_j^x)$$
(9)
Symmetry operations of the V15 molecule, rotations of $`\frac{2\pi }{3}`$ and $`\frac{4\pi }{3}`$ around the axis that passes through the centers of the hexagons ($`z`$), have to leave the form of the Hamiltonian invariant. This constraints the $`\stackrel{}{D}_{i,j}`$ to be related to one another. If

$$\stackrel{}{D}_{1,2}=D_{1,2}^x\text{ }\widehat{x}+D_{1,2}^y\text{ }\widehat{y}+D_{1,2}^z\text{ }\widehat{z}$$
(10)
then these three parameters determine the DM interaction in the other two $`J`$ bonds of the bottom hexagon as:

$`\stackrel{}{D}_{3,4}={\displaystyle \frac{1}{2}}\text{ }(D_{1,2}^x+\sqrt{3}\text{ }D_{1,2}^y)\text{ }\widehat{x}{\displaystyle \frac{1}{2}}(\sqrt{3}\text{ }D_{1,2}^x+D_{1,2}^y)\text{ }\widehat{y}+D_{1,2}^z\text{ }\widehat{z}`$ (11)
$`\stackrel{}{D}_{5,6}={\displaystyle \frac{1}{2}}\text{ }(D_{1,2}^x+\sqrt{3}\text{ }D_{1,2}^y)\text{ }\widehat{x}+{\displaystyle \frac{1}{2}}(\sqrt{3}\text{ }D_{1,2}^xD_{1,2}^y)\text{ }\widehat{y}+D_{1,2}^z\text{ }\widehat{z}`$ (12)
The sites $`1,\mathrm{},6`$ are those in the bottom hexagon. Due to absence of reflection symmetry with respect to the triangle plane the DM terms in the upper hexagon will differ from the ones in the lower one. However, to minimize the number of free parameters in the model, $`\stackrel{}{D}_{10,11}=\stackrel{}{D}_{1,2}`$. Then from symmetry $`\stackrel{}{D}_{12,13}=\stackrel{}{D}_{5,6}`$ and $`\stackrel{}{D}_{14,15}=\stackrel{}{D}_{3,4}`$.

In the same manner symmetry determines the elements of $`𝚪_{𝐢,𝐣}`$. Considering the symmetric part of the spin-orbit coupling again only on the strong bonds ($`J`$), and if

$$𝚪_{\mathrm{𝟏},\mathrm{𝟐}}=\left(\begin{array}{ccc}\mathrm{\Gamma }_{1,2}^{xx}& \mathrm{\Gamma }_{1,2}^{xy}& \mathrm{\Gamma }_{1,2}^{xz}\\ \mathrm{\Gamma }_{1,2}^{xy}& \mathrm{\Gamma }_{1,2}^{yy}& \mathrm{\Gamma }_{1,2}^{yz}\\ \mathrm{\Gamma }_{1,2}^{xz}& \mathrm{\Gamma }_{1,2}^{yz}& \mathrm{\Gamma }_{1,2}^{zz}\end{array}\right)$$
(13)
then symmetry determines $`𝚪_{𝐢,𝐣}`$ by:

$$𝚪_{𝐢,𝐣}=U\text{ }𝚪_{\mathrm{𝟏},\mathrm{𝟐}}\text{ }U^T$$
(14)
where $`U`$ is a rotation by $`\varphi `$ around the $`z`$ axis:

$$U=\left(\begin{array}{ccc}cos\varphi & sin\varphi & 0\\ sin\varphi & cos\varphi & 0\\ 0& 0& 1\end{array}\right)$$
(15)
$`𝚪_{\mathrm{𝟑},\mathrm{𝟒}}`$ and $`𝚪_{\mathrm{𝟓},\mathrm{𝟔}}`$ are determined in this manner and, as in the case of the DM terms, to minimize the number of free parameters in the model $`𝚪_{\mathrm{𝟏𝟎},\mathrm{𝟏𝟏}}=𝚪_{\mathrm{𝟏},\mathrm{𝟐}}`$. Then from symmetry $`𝚪_{\mathrm{𝟏𝟐},\mathrm{𝟏𝟑}}=𝚪_{\mathrm{𝟓},\mathrm{𝟔}}`$ and $`𝚪_{\mathrm{𝟏𝟒},\mathrm{𝟏𝟓}}=𝚪_{\mathrm{𝟑},\mathrm{𝟒}}`$.

The DM term added to the $`H_1`$ term in equation (1) is of the form:

$`H_{DM}=`$ $`\stackrel{}{D}_{1,2}\stackrel{}{S}_1\times \stackrel{}{S}_2+\stackrel{}{D}_{3,4}\stackrel{}{S}_3\times \stackrel{}{S}_4+\stackrel{}{D}_{5,6}\stackrel{}{S}_5\times \stackrel{}{S}_6+`$ (16)
$`\stackrel{}{D}_{10,11}\stackrel{}{S}_{10}\times \stackrel{}{S}_{11}+\stackrel{}{D}_{12,13}\stackrel{}{S}_{12}\times \stackrel{}{S}_{13}+\stackrel{}{D}_{14,15}\stackrel{}{S}_{14}\times \stackrel{}{S}_{15}`$ (17)
This term is considered as part of the perturbation, H1, and is scaled with the perturbation strength $`\lambda `$. From now on $`D^xD_{1,2}^x`$, $`D^yD_{1,2}^y`$ and $`D^zD_{1,2}^z`$. This minimum model is investigated for magnetic anisotropy and the parameters, $`D^x`$, $`D^y`$ and $`D^z`$, are determined by comparing with experiment.

We have found that the symmetric $`𝚪_{𝐢,𝐣}`$ terms, which are second order in spin-orbit coupling, have little effect in the magnetization curves. Once we have determined the parameters of the minimal model discussed above, we will return to this point. In the case where the moments lie on an open linear array it is possible to introduce local axes so that the DM term, $`\stackrel{}{D}_{ij}\stackrel{}{S}_i\times \stackrel{}{S}_j`$, is absorbed into an effective symmetric term . However, the connectivity of the bonds between the spins on the V sites ensures that this term can not be transformed away in this manner.

## IV Results

The effect of the DM terms in the energy spectra is investigated by varying the strength of $`D^x`$, $`D^y`$ and $`D^z`$. The average energies of the two doublets and the quartet are plotted for different values of the $`\stackrel{}{D}`$ vector in figure 3. The changes in the energy become more pronounced as $`|\stackrel{}{D}|`$ increases. However, one sees that even DM terms $`100`$ K lead to energy shifts $`20`$ K. The deviations from the average energies of the two doublets are plotted in figure 4, while for the quartet in figure 5. It can be seen that even when $`|\stackrel{}{D}|100`$ K the sizes of the energy splittings in the two doublets and the quartet are modest. This suggests that the energy level differences estimated from the experiment are comparatively insensitive to $`\stackrel{}{D}`$ so that it can not alone determine $`\stackrel{}{D}`$.

These spectra arise from the form of the effective Hamiltonian for the triange spins, $`H^{eff}`$, calculated with perturbation theory. The effective Hamiltonian for the triangle spins is of the form $`H^{eff}=H^{diag}+H^{offdiag}`$, where $`H^{diag}`$ is a constant diagonal matrix in the basis of total $`S^z`$ states which is three orders of magnitude larger than the terms in $`H^{offdiag}`$ and depends on the magnitude of the anisotropy term, as can be seen in figure 3. All the total $`S^z`$ sectors of the triangle spins are coupled in $`H^{offdiag}`$ even when the anisotropy term is $`\stackrel{}{D}^z_{ij}\stackrel{}{S}_i\times \stackrel{}{S}_j`$ on the $`800`$ K hexagonal bonds, case (I) in figure 4 and figure 5.

By contrast, when one considers a three-spin model with a Heisenberg coupling, $`J_0_{ij}\stackrel{}{S}_i\stackrel{}{S}_j`$, and an anisotropy term, $`\stackrel{}{D}^z_{ij}\stackrel{}{S}_i\times \stackrel{}{S}_j`$, between the triangle spins, the quartet remains degenerate while the doublets are split into time-reversed pairs separated by $`\sqrt{3}D^z`$. Unlike the result of this three-spin model, the energies of all quartet levels are changed by the anisotropy term, as can be seen from figures 4 and 5, and the splitting of the higher energy multiplet is larger than the low energy one. This points to the difficulty in understanding the effect of the anisotropy terms on $`H^{eff}`$ in terms of a simple model for the triangle spins. More generally it supports the conclusions of De Raedt et al.  on the dangers of Hilbert space truncation.

### A Magnetization

The magnetization of V15 was calculated for various magnetic fields and temperatures as a function of the DM terms. In the absence of DM terms and at zero temperature, $`S^z`$ is a good quantum number and the magnetization stays constant except at critical fields where there are level crossings. In the first level crossing the ground state selects the state of the zero field ground state doublet with the higher expectation value of the spin, while in the second it switches from this doublet to one of the quartet states. When the temperature differs from zero, the jumps are smeared out and at high enough temperatures there are no sudden changes in the magnetization.

When the DM terms are non-zero at zero temperature $`\stackrel{}{S}^2`$ does not commute with the Hamiltonian, and the same is true for $`S^z`$ if either $`D^x`$ or $`D^y`$ is non-zero. The magnetization along the field was calculated for $`D^x=50`$ K and $`D^x=100`$ K at a very low temperature, $`T=0.001`$ K, in the case of a field along the $`x`$ axis and a field along the $`z`$ axis and is plotted in figure 6. The gyromagnetic ratio used is $`g=2.0023`$ . It observed that the magnetization depends on the degree of anisotropy, $`D^x`$. In addition, it is seen that the jump is more smeared out as the magnitude of $`D^x`$ increases from $`0`$ to $`100`$ K. Magnetic anisotropy is demonstrated by the dependence of the magnetization along the direction of the field on the direction of the field itself.

The broadening of the spin transition can also be observed in figure 7. There the magnetization $`<S^z>`$ is plotted as a function of the magnetic field for four different temperatures and a DM term of magnitude $`\sqrt{3}`$ $`100`$ K. The magnetization curve is almost identical for $`T=0.001`$ K and $`T=0.1`$ K. This shows that the broadening of the transition is not a temperature effect, but it is associated with the DM terms. For high enough temperatures the susceptibility is almost constant with the field.

In figure 8 the calculated magnetization is plotted for three different values of $`D^x`$ with $`D^y=D^z=0`$. It is seen that the width of the transition gets bigger with the increase of the DM terms. The origin of the width can be seen in the energy spectrum plotted in figure 9. The presence of matrix elements between the states in the vicinity of the transition field ($`3`$ T) leads to level anticrossings. Pairs of degenerate states from the lower multiplet in zero field split up with one hybridizing with a member of the quartet multiple with $`S^z\frac{3}{2}`$, while the other has no matrix element with this state.

The width of the magnetization transition is an indication of the effect of the DM terms, as can be seen from figures 7 and 8. The transition is characterized by its midpoint $`B_{transition}`$ and width-difference $`\mathrm{\Delta }B`$ in magnetic field values where it has changed by 10 $`\%`$ from its low and high field average spin values on either side of the transition. The results for various DM terms are shown in table I, where it can be seen that $`B_{transition}`$ depends sensitively on the form of the DM terms. This dependence on the DM terms is illustrated in figure 10 where we plot $`B_{transition}`$ versus $`\mathrm{\Delta }B`$ for two different forms of the DM vector, one in which $`D^x0`$ K, $`D^y=D^z=0`$ and the other in which $`D^z0`$ K, $`D^x=D^y=0`$. The values of $`B_{transition}`$ and the corresponding values of $`\mathrm{\Delta }B`$ shown are for different magnitudes of $`D^x`$ and $`D^z`$ at two different temperatures. Although the magnetization curves for these two temperatures appear to lie on top of one another in figure 7, there is a difference which should be taken into account in quantitative comparison with data. The effect of increasing temperature is seen to push $`B_{transition}`$ to higher fields, while $`\mathrm{\Delta }B`$ is seen to be non-monotonic. For small magnitudes of the DM term $`\mathrm{\Delta }B`$ is larger at $`T=0.1`$ K than at $`T=10^3`$ K but this is reversed for a sufficiently large magnitude depending on the two forms considered.

These temperature dependences can be understood as follows. $`B_{transition}`$ moves to higher applied fields with increasing temperature because the higher energy low magnetization state contributes significantly to the thermodynamic average for larger energy differences between the ground and excited state, that is for larger applied fields beyond the field at which the two states are degenerate. The temperature dependence of $`\mathrm{\Delta }B`$ is more complicated. Both the temperature and the $`\stackrel{}{D}_{ij}`$ vectors contribute to $`\mathrm{\Delta }B`$. Because the temperatures considered here are so low compared to the variation in the magnitudes of the $`\stackrel{}{D}_{ij}`$ vectors, the difference in the smearing effects of the different temperatures, $`0.1`$ K, can be neglected compared to that due to the hybridization of the levels near $`B_{transition}`$. The increase in $`\mathrm{\Delta }B`$ for $`T=10^3`$ K compared to $`T=0.1`$ K arises because the transition takes place at lower values of the applied field at $`T=10^3`$ K, $`2.60`$ T, instead of $`2.68`$ T at $`T=0.1`$ K, and because the effect of the $`\stackrel{}{D}_{ij}`$ vectors is larger at lower values of the applied field once the energy levels are close. The matrix elements associated with the $`\stackrel{}{D}_{ij}`$ vectors are larger compared to the Zeeman energies in the effective Hamiltonian. As a result increasing the magnitude of $`\stackrel{}{D}_{ij}`$ vectors leads to a larger effect on $`\mathrm{\Delta }B`$ at lower temperatures. It was also observed that the results do not depend on the sign of the DM terms.

The above results were calculated with a gyromagnetic ratio $`g=2.0023`$. A similar analysis for experimentally determined temperature dependence of the magnetization at $`T=0.1`$ K gives $`\mathrm{\Delta }B=0.77\pm 0.10`$ T and $`B_{transition}=2.84\pm 0.05`$ T . In order to compare with this data, an average of the values found from EPR spectra for the gyromagnetic ratio is used . They are $`g_a=g_b=1.95`$ and $`g_c=1.98`$, where $`a`$, $`b`$ and $`c`$ are the crystallographic axes of the molecule. Using $`g=1.97`$ and $`D^x=70`$ K, $`D^y=D^z=0`$, it is found that $`\mathrm{\Delta }B=0.75`$ T and $`B_{transition}=2.81`$ T, which is within the uncertainty of the experimental value. Therefore, the inclusion of DM terms accounts for the rounding of the magnetization at very low temperatures and provides a description of the magnetic anisotropy of V15. In addition, the gap between the doublets and quartet, $`3.7`$ K, is consistent with the estimates given by Chiorescu et al.  and Chaboussant et al.  for the effective coupling of the spins in the triangle, induced by interactions between these and the spins in the hexagons. Both groups found $`J_02.4`$ K which gives a gap $`\frac{3}{2}J_03.7`$ K. We also found that the tunnel splitting between the two doublets is $`\mathrm{\Delta }_0=13`$ mK, in agreement with the $`10^2`$ K estimate . This is an order of magnitude smaller than the estimates given by Chiorescu et al.  ($`0.05`$ K) or Chaboussant et al. ($`0.1`$ K). We note that the calculated splittings of the quartets, shown in figure 5, are an order of magnitude larger than those of the doublets, shown in figure 4, which may account for the discrepancy between experimental estimates and our results.

In arriving at the form and magnitude of the $`\stackrel{}{D}_{ij}`$ we have relied on a simple model in which the $`\stackrel{}{D}_{ij}`$ are non-zero only on the strongest bonds and it is not obvious how model dependent our results are in the absence of more detailed data. The original estimates by Moriya  relied on a superexchange model where magnetic order is driven by strong short-range correlations due to strongly suppressed hopping between sites. Estimates of the $`\stackrel{}{D}_{ij}`$’s then come from considering the mixing of orbitals by angular momentum matrix elements. In V15 interactions between ions on hexagon sites are mediated by pairs of oxygen atoms. Depending on whether one of these oxygen atoms is linked to an arsenic atom or another vanadium atom, the Heisenberg coupling is either $`800`$ K or $`150`$ K. In these circumstances it is not clear how reliable estimates are based on the simple tight-binding approach. The magnitude of the $`\stackrel{}{D}_{ij}`$ vectors found here is comparable to the Heisenberg couplings. This is larger than the estimate based on Moriya’s original approach which gives $`\frac{|\stackrel{}{D}_{ij}|}{J}\frac{\mathrm{\Delta }g}{g}`$, where $`\mathrm{\Delta }g`$ is the shift in the gyromagnetic ratio from that of a free electron. However, Katsnelson et al.  also found strong Dzyaloshinskii-Moriya interactions with values comparable to the isotropic constants in their analysis, which was based on an effective eight spin model of neutron scattering and EPR data on Mn12.

The determination of $`\stackrel{}{D}_{ij}`$ is carried out in the absence of the symmetric anisotropy term, $`𝚪_{𝐢,𝐣}`$. We investigated the effect of including $`𝚪_{𝐢,𝐣}`$’s with a range of magnitudes for its elements with different $`\stackrel{}{D}_{ij}`$’s. Since the $`\stackrel{}{D}_{ij}`$’s are first order in the spin-orbit interaction while the $`𝚪_{𝐢,𝐣}`$’s are second order, we have assumed the magnitudes of the symmetric anisotropy terms to be smaller than the antisymmetric ones. For illustration we show in fig. 11 the dependence of $`B_{transition}`$ and the corresponding values of $`\mathrm{\Delta }B`$ for different forms of $`𝚪_{𝐢,𝐣}`$ given in table II. In this case $`D^x=70`$ K, $`D^y=D^z=0`$ with $`g=2.0023`$. One sees that the values of $`B_{transition}`$ and $`\mathrm{\Delta }B`$ are not significantly changed from the values given by the DM terms alone. It was observed here that the results depend on the sign of the $`𝚪_{𝐢,𝐣}`$ terms.

### B Magnetic Relaxation

Apart from equilibrium magnetization measurements, the magnetization of V15 in a varying external field has been measured. The data shows a so-called “butterfly” hysteresis, where the magnetization initially increases with the field, then reaches a plateau, and eventually approaches saturation . This has been attributed to non-equilibrium distribution of phonons, and the spin-phonon mechanism for magnetic relaxation is called “hole-burning”. Here the role of DM terms in magnetic relaxation is investigated.

In order to calculate the magnetization in a swept field as in the experiments it is necessary to calculate the density matrix for the eight low energy states. However, this is a prohibitively expensive numerical calculation if one is to reach applied fields of $`7`$ T with the sweeping rate of $`0.14`$ T/s. This is because the characteristic time, $`\tau _0=\frac{\mathrm{}}{k_BK}=7.64\times 10^{12}`$ s, associated with anisotropy terms is so short compared to the $`50`$ s required to reach $`7`$ T, that the number of time steps needed is very large. The presence of rapidly varying phases in the elements of the density matrix requires a very fine time step $`\tau _0`$. The computational effort to get accurate results is enormous. Therefore, rather than consider a time dependent magnetic field which changes at a constant rate, we consider the case in which a $`10`$ T field along the $`z`$ direction is suddenly switched on, off, or reversed. In the absence of magnetic anisotropy there is no change in $`<S^z>`$ with time. However, this is not the case in the presence of DM terms. As can be seen in figure 9, there is hybridization between the zero field states as a function of applied field, and matrix elements responsible for this lead to a time dependent $`<S^z>`$ when the applied field is suddenly changed.

We illustrate this point by calculating the magnetization for the case where the magnetic field has an initial strength of $`10`$ T, and $`D^x=80`$ K, $`D^y=30`$ K, $`D^z=0`$ and $`g=2.0023`$. The ground state magnetization is $`<S^z>=1.49`$. Then the magnetic field is suddenly switched off, and the evolution of the magnetization is calculated in the basis of the eight lowest states of the same problem in zero magnetic field. Any excited state higher in energy is separated by a gap of order $`J`$ from the lowest energy manifold, and the probability of the state of the system belonging to this manifold is $`0.9997`$, so this approximation is well justified.

The magnetization varies harmonically with time, as seen in figure 12, and the frequencies of its variation are determined by the energy differences of the system at zero magnetic field. After a time equal to $`100\tau _010^{10}`$ s, the field is turned back on to $`10`$ T and $`<S^z>`$ now varies with different frequencies and amplitude. The presence of more than one frequency in the time dependence points to the inadequacy of a two-level model and suggests that there are significant corrections to any Landau-Zener treatment of the relaxation. The average value of the magnetization is different in the two cases, and in particular the average value after the field is turned on again depends on the wavefunction at the exact moment the field is switched. The origin of this effect is the magnetic anisotropy of the DM terms, a source of magnetic decoherence for the molecule. The calculation can also be done starting from the zero field ground states in a non-zero magnetic field, and after $`10^{10}`$ s switching off the field suddenly. The time dependence is similar to the previous case, and once again the average value of $`<S^z>`$ after the field has been switched off depends on the state of the system at the moment it was switched off. The DM terms therefore provide a source of spin decoherence with a decoherence time $`10^{12}`$ s.

In both cases the calculation was done assuming instantaneously switched on and off fields and demonstrates magnetic decoherence and hysteresis in the presence of DM terms. If a time varying field is introduced in the calculation, the presence of rapidly varying phases requires a very fine time step for the calculation of the time evolution of the magnetization. This time step has to be a fraction of $`10^{12}`$ s, and a calculation of the evolution of the magnetization for realistic times would require an enormous amount of computational effort. However, the small time scale calculation demonstrates the intrinsic magnetic decoherence in the presence of magnetic anisotropy. Hysteresis is seen in experiment and attributed to the spin-phonon interaction . The present mechanism is independent of temperature and can be expected to provide the dominant relaxation process at very low temperatures.

## V Conclusion

Degenerate perturbation theory was used to study magnetic anisotropy in the molecule V15. The DM term introduced in the model quantitatively accounts for the rounding of transitions in the magnetization at very low temperatures, their width at midpoint, through the avoided level crossing, as well as observed zero field splitting of the lowest four states.

We have shown that the DM interaction is also responsible for spin decoherence and magnetic hysteresis. This is present even at very low temperatures where thermal mechanisms are absent. This source of decoherence may pose a problem for the application of magnetic molecules in quantum computing, since it is necessary to control the time evolution of states with external magnetic fields and electromagnetic pulses . In order for a molecule with a DM interaction to be suitable for quantum computing it would be necessary to carry out operations on a time scale much shorter than $`\frac{\mathrm{}}{\mathrm{\Delta }E}`$, where $`\mathrm{\Delta }E`$ is comparable to the magnitude of the DM interaction. In the present example this is $`10^{12}`$ s.

The use of perturbation theory to calculate the low energy and temperature properties of V15 is an alternative to exact diagonalization which has been discussed by Honecker et al.  for Fe6 and Fe8. Perturbation theory works particularly well in the case of V15 because of the strength of the hexagonal bonds. However, we have also applied perturbation theory to a system of $`20`$ $`S=\frac{1}{2}`$ spins by calculating the effect of quantum fluctuations around the classical ground state. The dimension of the Hilbert space in that case is almost exactly the same as Fe8, $`10^6`$. It is possible that this approach will also succeed in systems with larger, more classical moments and larger Hilbert spaces. This would allow a more microscopic approach to the treatment of Mn12, which may resolve some of the questions raised by Zhong et al. .

## VI Acknowledgements

This work was carried out at the Center for Computational Research at SUNY Buffalo. N.P.K. is supported by a Marie Curie Fellowship of the European Community program Development Host Fellowship under contract number HPMD-CT-2000-00048.

Current address: Department of Physics, Buffalo State College, 1300 Elmwood Avenue, Buffalo, NY 14222.


format: text

text:
# Tuning the dipolar interaction in quantum gases

## Abstract

We have studied the tunability of the interaction between permanent dipoles in Bose-Einstein condensates. Based on time-dependent control of the anisotropy of the dipolar interaction, we show that even the very weak magnetic dipole coupling in alkali gases can be used to excite collective modes. Furthermore, we discuss how the effective dipolar coupling in a Bose-Einstein condensate can be tuned from positive to negative values and even switched off completely by fast rotation of the orientation of the dipoles.

In the Bose-Einstein condensates (BECs) created thus far  atoms interact essentially only at very short distance. At low temperatures, the interaction between these neutral atoms can be typically characterized by a single parameter, the s-wave scattering length. The magnitude and sign of the s-wave scattering length can be changed by tuning the external magnetic fields in the vicinity of a so-called Feshbach resonance . This technique has opened new possibilities for the study and manipulation of BECs. In particular, the value of the s-wave scattering length can be varied in real time, allowing, for example, the observation of collapsing and exploding BECs .

In recent years an increasing interest is devoted to dipole-dipole interactions in ultracold gases . Dipolar interactions would largely enrich the variety of phenomena to be observed due to their long-range and anisotropic character. If atomic dipole moments are sufficiently large, the resulting dipole-dipole forces may influence, or even completely change the properties of Bose gases , the conditions for BCS transition in Fermi gases , or the phase diagram for quantum phase transitions in ultracold dipolar gases confined in optical lattices . The interplay of short-range scattering and long-range interaction may give rise to phenomena like ferromagnetic order and spin waves . Moreover the dipolar particles are considered to be promising candidates for the implementation of fast and robust quantum-computing schemes .

In this Letter, we consider the case of permanent (magnetic) dipoles. For alkali atoms with a magnetic moment $`m`$ of $`\mu _B`$, where $`\mu _B`$ is a Bohr magneton, the long range part of the magnetic dipole-dipole interaction is generally neglected since it is very small compared to the s-wave pseudo-potential. Here, we show that even in alkali BECs this dipolar interaction can be made visible if the dipole-dipole interaction is modulated by rotating the atomic dipoles. For appropriate rotation frequencies, a coupling to elementary excitations can be achieved leading to a dynamic growth of the amplitude of the excitation. In particular, we will show how the dipolar coupling in an alkali condensate can be used to excite the quadrupole mode and how it is possible to extract quantitative information about the ratio of the dipolar coupling and the usual s-wave coupling.

In ensembles of atoms with a larger magnetic moment, such as chromium where $`m=6\mu _B`$, the magnetic dipole-dipole interaction energy may become comparable in size with the s-wave contribution to the mean-field energy. Therefore, even the static properties of a condensate such as the aspect ratio are significantly altered. Here, we show how the effective strength of the interaction and thus its effect on the static properties of a condensate can be varied by fast rotation of the atomic dipoles. In the case of a dipolar interaction which is stronger than the s-wave interaction only certain combinations of trap geometry and orientation of the dipoles lead to a stable (or metastable) condensate . By tuning the dipolar interaction unstable condensates can thus be transformed into stable ones and viceversa. In this paper, we study the interaction between permanent magnetic dipoles. However, the methods we propose can also be adapted for particles with electric dipole moment such as heteronuclear molecules. In this case the dipolar interaction is expected to dominate over the s-wave interaction and thus the tuning of the interaction might be necessary to reach stability.

*Manipulating the magnetic dipolar coupling in a BEC*—The long-range part of the interaction between two magnetic dipole moments $`𝐦_1`$ and $`𝐦_2`$ of atoms located at $`𝐫_1`$ and $`𝐫_2`$ takes the familiar dependence

$`U_{\mathrm{dd}}(𝐫)={\displaystyle \frac{\mu _0}{4\pi }}{\displaystyle \frac{3(𝐦_1\widehat{𝐫})(𝐦_2\widehat{𝐫})(𝐦_1𝐦_2)}{r^3}}`$ (1)
as a function of the interatomic distance $`𝐫=𝐫_2𝐫_1`$, where $`\widehat{𝐫}=𝐫/r`$ and $`\mu _0`$ is the magnetic permeability of the vacuum. We consider the case of a general *time-dependent* homogeneous magnetic field

$`𝐁(t)`$ $`=`$ $`B𝐞(t),`$ (2)
slowly-varying with respect to the Larmor frequency $`\omega _{Larmor}=mB/\mathrm{}`$, where $`m=|𝐦_1|=|𝐦_2|`$ for identical atoms. A spin-polarized atomic ensembles will adiabatically follow the external magnetic field $`B(t)`$ and the resulting interatomic energy becomes

$`U_{\mathrm{dd}}(𝐫,t)={\displaystyle \frac{\mu _0m^2}{4\pi }}{\displaystyle \frac{3(𝐞(t)\widehat{𝐫})^21}{r^3}}`$ (3)
and is thus *time-dependent*.

Our analysis is based on the mean-field approach in the Thomas-Fermi limit already used in the context of dipolar Bose gases . Such a description can be accomplished through the hydrodynamic equations for superfluids :

$`{\displaystyle \frac{n}{t}}`$ $`=`$ $`(n𝐯),`$ (4)
$`{\displaystyle \frac{𝐯}{t}}`$ $`=`$ $`\left({\displaystyle \frac{v^2}{2}}+{\displaystyle \frac{\delta \mu }{M}}\right),`$ (5)
where $`n`$ is the density, $`𝐯`$ the velocity field and $`\delta \mu `$ the difference between the mean-field plus external fields and the chemical potential $`\mu `$ given by

$`\delta \mu `$ $`=`$ $`gn+\mathrm{\Phi }_{\mathrm{dd}}+{\displaystyle \frac{M\omega _0^2r^2}{2}}\mu .`$ (6)
Here $`M`$ is the atomic mass, $`M\omega _0^2r^2/2`$ is an isotropic harmonic potential which could be created using optical dipole forces, $`g=4\pi \mathrm{}^2a/M`$ and $`a`$ is the s-wave scattering length which we assume to be positive. $`\mathrm{\Phi }_{\mathrm{dd}}(𝐫,t)`$ is the mean-field potential generated by the (time dependent) dipolar interaction

$`\mathrm{\Phi }_{\mathrm{dd}}(𝐫,t)={\displaystyle 𝑑𝐫^{}U_{\mathrm{dd}}(𝐫𝐫^{},t)n(𝐫^{})}.`$ (7)
In the Thomas-Fermi limit it is useful to define a dimensionless quantity

$`\epsilon _{\mathrm{dd}}={\displaystyle \frac{\mu _0m^2M}{12\pi \mathrm{}^2a}},`$ (8)
which is a measure of the strength of the dipole-dipole interaction relative to the s-wave scattering energy. Modifications in the ground state condensate density can only be expected if $`\epsilon _{\mathrm{dd}}`$ is appreciably different from zero.

For chromium the s-wave scattering length $`a_{\mathrm{Cr}}`$ is not known. If we assume that $`a_{\mathrm{Cr}}`$ is equal to the sodium scattering length ($`a=2.8`$ nm) the dipolar strength parameter becomes $`\epsilon _{\mathrm{dd}}^{\mathrm{Cr}}=0.29`$ ($`m=6\mu _B`$). Therefore, the dipolar force is expected to lead to visible modification in the density of a chromium BEC. Differently, in both rubidium and sodium the magnetic dipolar energy is rather small compared to the mean-field s-wave scattering energy. The values of $`\epsilon _{\mathrm{dd}}`$ are $`\epsilon _{\mathrm{dd}}^{\mathrm{Rb}^{87}}=0.0064`$ and $`\epsilon _{\mathrm{dd}}^{\mathrm{Na}}=0.0035`$, in the doubly spin-polarized ground-state, where $`m=\mu _B`$. For heteronuclear molecules with a permanent electric dipole moment of one Debye the corresponding dipolar coupling will be increased by approximatively the inverse square of the fine structure constant. This implies that $`\epsilon _{\mathrm{dd}}`$  is on the order of $`10^2`$.

*Probing the quadrupole modes in sodium and rubidium*—The time modulation of the dipolar coupling (3) can be used to resonantly excite a quadrupole mode even for sodium and rubidium condensates. For simplicity we restrict ourselves to the case of a condensate in the Thomas-Fermi regime confined in an isotropic harmonic trap with frequency $`\omega _0`$ (for example $`\omega _0/2\pi =100`$ Hz). Because $`\epsilon _{\mathrm{dd}}`$ is small the dipole-dipole mean-field potential $`\mathrm{\Phi }_{\mathrm{dd}}`$ can be calculated using the symmetric equilibrium density distribution $`n_0=M\omega _0^2(R^2r^2)/2g`$, where $`R`$ is its radius

$`\mathrm{\Phi }_{\mathrm{dd}}(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{\epsilon _{\mathrm{dd}}M\omega _0^2}{5}}\left[3(𝐞(t)\widehat{𝐫})^21\right]r^2r<R,`$ (9)
$`=`$ $`{\displaystyle \frac{\epsilon _{\mathrm{dd}}M\omega _0^2}{5}}\left[3(𝐞(t)\widehat{𝐫})^21\right]{\displaystyle \frac{R^5}{r^3}}r>R.`$ (10)
The potential outside the condensate region ($`r>R`$) corresponds exactly to the field generated by N dipoles located in the center of the condensate. Differently, inside the condensate ($`r<R`$) the mean-field potential has an anisotropic but harmonic dependence.

If the interaction is driven on a resonance frequency the evolution of the condensate can be modified even for very small dipolar coupling. Consider a magnetic field rotating in the plane $`xy`$ as $`𝐁(t)=B𝐞(t)`$, where

$`𝐞(t)=\mathrm{cos}(\mathrm{\Omega }t)\widehat{𝐱}+\mathrm{sin}(\mathrm{\Omega }t)\widehat{𝐲}.`$ (11)
The resulting dipolar interatomic energy (3) and its corresponding contribution to the mean-field potential (9-10), will have a ’rotating’ part . For small dipolar coupling within alkali gases the non oscillating part of (9) can be neglected and the time-dependent part will resonantly couple to the ’transversal’ quadrupole mode when $`\mathrm{\Omega }\omega _0/\sqrt{2}`$ (see Ref. ). After linearization of the equations of motion (4) and (5) we obtain

$`{\displaystyle \frac{^2\delta n}{t^2}}=\left({\displaystyle \frac{gn_0}{M}}\left(\delta n\right)\right)+\left({\displaystyle \frac{n_0}{M}}\left(\mathrm{\Phi }_{\mathrm{dd}}\right)\right),`$ (12)
where $`\delta n=nn_0`$ is the density displacement around the equilibrium density $`n_0`$. The time-dependent part of $`\mathrm{\Phi }_{\mathrm{dd}}`$ corresponds in Eq. (12) to the rotation of an anisotropic harmonic trap with deformation

$`\epsilon ={\displaystyle \frac{\omega _X^2\omega _Y^2}{\omega _X^2+\omega _Y^2}}={\displaystyle \frac{3}{10}}\epsilon _{\mathrm{dd}},`$ (13)
where $`\omega _X`$ and $`\omega _Y`$ are the maximum and minimum values of the rotating trap frequencies. Neglecting the time-independent part of (9) the density will evolve as

$`\delta n={\displaystyle \frac{3\epsilon _{\mathrm{dd}}M\omega _0^2}{10g}}r^2\mathrm{sin}^2\theta \mathrm{sin}(2\varphi 2\mathrm{\Omega }t)\mathrm{\Omega }t,`$ (14)
where $`\varphi `$ is the azimuthal angle ($`r<R`$). It is noteworthy that the rotating part  of (9) couples *only* the first quadrupole mode in the Thomas-Fermi limit . From the study of the evolution of the condensate deformation

$`\alpha ={\displaystyle \frac{x^2y^2}{x^2+y^2}}={\displaystyle \frac{3}{5}}\epsilon _{\mathrm{dd}}\mathrm{\Omega }t\mathrm{sin}(2\mathrm{\Omega }t)`$ (15)
we can extract information about the ratio of the dipolar coupling versus the standard mean-field energy. If the quality factor $`Q=\sqrt{2}\omega _0\tau `$ of such a resonance, where $`\tau `$ is the corresponding decay time, is assumed to be 50  the amplitude of the oscillation of the condensate anisotropy $`\alpha `$ is expected to be $`|\alpha |=3Q\epsilon _{\mathrm{dd}}/109.5\times 10^2`$ for rubidium (see Fig. 1). We have performed our calculation in the linear regime and for an isotropic trap. The non-linearity of the mean-field is expected to give a significant correction for larger $`Q`$. The extension of the calculation to an axially symmetric cigar-shaped trap is expected to yield similar results.

*Tunability of the magnetic dipolar interaction:* — We consider a condensate in an axially symmetric harmonic trap, with radial and axial trap frequencies $`\omega _r`$ and $`\omega _z`$. The magnetic field

$`𝐁(t)=B\left[\mathrm{cos}(\phi )\widehat{𝐳}+\mathrm{sin}(\phi )\left(\mathrm{cos}(\mathrm{\Omega }t)\widehat{𝐱}+\mathrm{sin}(\mathrm{\Omega }t)\widehat{𝐲}\right)\right]`$ (16)
is a combination of a static magnetic field $`B_z`$ directed along the $`z`$ direction and a fast rotating field $`B_\rho `$ in the radial plane. The frequency is chosen such that the atoms are not significantly moving during the time $`\mathrm{\Omega }^1`$; while the magnetic moments will follow adiabatically the external field $`𝐁(𝐭)`$; this corresponds to $`\omega _{Larmor}\mathrm{\Omega }\omega _{trap}`$, where $`\omega _{trap}`$ is any of the external trap frequencies. In this limit we can consider the average of the interaction (3) in the period $`2\pi /\mathrm{\Omega }`$, resulting in the cylindrically symmetric interatomic potential

$`U_{\mathrm{dd}}(\stackrel{}{r})`$ $`=`$ $`{\displaystyle \frac{\mu _0m^2}{4\pi }}\left({\displaystyle \frac{3\mathrm{cos}^2\phi 1}{2}}\right)\left(3\mathrm{cos}^2\theta 1\right).`$ (17)
The averaged interaction energy (17) equals $`(3\mathrm{cos}^2\phi 1)/2`$ times the interatomic energy of dipoles aligned along the $`z`$-axes. Note that this factor can be changed continuously from $`1`$ to $`1/2`$, providing the possibility to change the dipolar interaction from attractive to repulsive. At a particular angle $`\phi _M=54.7^o`$ the dipolar interaction averages to zero. This angle, the so-called magic-angle, is well known in solid-state nuclear magnetic resonance (NMR) technology . The tunability of the dipolar interaction will allow to observe changes in the condensate eccentricity which should be a $`10\%`$ effect in Cr assuming $`a_{\mathrm{Cr}}=a_{\mathrm{Na}}`$ or in the collective excitation frequencies .

*Stabilizing strongly interacting dipolar gases:*—Let us now apply this tunability to the stability diagram of dipolar BECs. We consider the case in which the dipolar coupling is larger than the s-wave pseudopotential ($`\epsilon _{\mathrm{dd}}1`$). Here the tunability can stabilize strongly interacting dipolar gases like heteronuclear molecules which otherwise exhibit instability over a wide range of parameters. This regime, considered also in Ref. , may be achieved in a cloud of magnetic dipolar atoms, by making use of a Feshbach resonance to reduce the s-wave scattering length.

In a spin polarized cloud the stability diagram as well as the condensate anisotropy can be studied in the Thomas-Fermi limit as function of $`ϵ_{\mathrm{dd}}`$ (see Ref. ). When $`ϵ_{\mathrm{dd}}1`$ the cloud may be unstable against collapse, i.e. can be always possible to find a density configuration of arbitrary large negative energy. The stability as $`ϵ_{\mathrm{dd}}`$ and the trap anisotropy $`l=\omega _z/\omega _r`$ are varied can be investigated in analogy to Ref. .

Within the Thomas-Fermi approximation, the stability analysis can be made using a scaling variational approach . Assuming a density for the condensate of the form

$`n(x,y,z)={\displaystyle \frac{1}{\lambda _r^2\lambda _z}}n_0({\displaystyle \frac{x}{\lambda _r}},{\displaystyle \frac{y}{\lambda _r}},{\displaystyle \frac{z}{\lambda _z}}),`$ (18)
where $`\lambda _r`$ and $`\lambda _z`$ are scaling variational parameters and $`n_0`$ is a spherically symmetric density distribution (for instance the usual Thomas-Fermi inverted parabola). The condensate anisotropy $`\lambda _r/\lambda _z`$ can be found by minimizing the expectation value of the total energy $`H_{\mathrm{tot}}=H_{\mathrm{ho}}+H_{\mathrm{int}}`$ where

$`H_{\mathrm{ho}}=c_1Nm\left(2\omega _r^2\lambda _r^2+\omega _z^2\lambda _z^2\right)`$ (19)
is the harmonic potential energy and

$`H_{int}`$ $`=`$ $`c_2{\displaystyle \frac{N^2g}{\lambda _r^2\lambda _z}}\left[1\epsilon _{\mathrm{dd}}\left({\displaystyle \frac{3\mathrm{cos}^2\phi 1}{2}}\right)f\left({\displaystyle \frac{\lambda _r}{\lambda _z}}\right)\right]`$ (20)
the mean-field interaction energy, where $`c_1`$ and $`c_2`$ are numerical constants  and

$`f(\kappa )`$ $`=`$ $`{\displaystyle \frac{1+2\kappa ^2}{1\kappa ^2}}{\displaystyle \frac{3\kappa ^2\mathrm{tanh}^1\sqrt{1\kappa ^2}}{\left(1\kappa ^2\right)^{3/2}}}`$ (21)
derives from an angular integration.

The stability can be varied for a fixed value of $`\epsilon _{\mathrm{dd}}`$ and the trap anisotropy $`l`$ just varying one easy experimental parameter, namely the angle $`\phi `$. In Fig. 3 we show one example of how the stability phase diagram for $`\epsilon _{\mathrm{dd}}=4`$ can be varied as function of the angle $`\phi `$. For very large $`\epsilon _{\mathrm{dd}}`$ the stable region in the $`\phi `$-$`l`$ plane becomes a thin region around the magic angle line.

*Conclusion*—We have shown that the anisotropy of the dipole-dipole interaction of a polarized gas can be used to tune the strength of the dipolar coupling and to excite collective excitations. Tuning the dipolar coupling will become an important tool for designing atomic quantum gases with novel properties and dynamical aspects. This approach in optically trapped Bose condensed gases should allow to make even very weak dipolar coupling visible with experimentally feasible parameters (like $`B`$ fields in the Gauss region and spinning frequencies up to several $`10`$ kHz). The tuning of the interaction might become a relevant technique for stabilizing strongly interacting dipolar systems like heteronuclear molecular gases.

*Acknowledgement:*—Valuable discussion with G. Denninger are acknowledged; funding was provided by the RTN network ”cold quantum gases” under the contract number HPRN-CT-2000-00125, and the Deutsche Forschungsgemeinschaft.


text:
# 1 Introduction

## 1 Introduction

Glassy systems are very common in nature, yet not quite well understood. Lattice models may serve as good candidates to describe some properties of these systems. In recent years some interest has been raised by remarkably simple Ising spin systems that originated from the consideration of a model of random surfaces in the context of string theory , . The version of this model in a discretized space defines the so-called gonihedric spin model which consists, in general, of an Ising model with finely tuned nearest, next-to-nearest neighbour and plaquette interactions. The relation among the couplings of the hamiltonian depends also on the dimensionality of the system. The geometric origins of the model show up in a remarkable simple way of writing the energy of a given configuration: the surfaces corresponding to the interfaces between up and down spins are weighed by $`E=n_2+\kappa n_4`$, where $`n_2`$ is the number of edges of such an interface and $`n_4`$ is the number of four plaquettes that share a common link. The parameter $`\kappa `$ can thus be interpreted as an indicator of the self-avoidance of the model. Notice that there is no microscopic surface tension.

Up to now there has been a considerable amount of numerical work on the three–dimensional case, to which we shall refer in the following and also some peliminiary results in four–dimensional case . In two dimensions, the model with $`\kappa =0`$ is actually trivial (no phase transition ) but the solution with $`\kappa 0`$ is unknown.

For the three dimensional case which is of our concern in this paper the Hamiltonian of the model takes the form

$$H(\sigma )=2\kappa \underset{\stackrel{}{r},\stackrel{}{\alpha }}{}\sigma _\stackrel{}{r}\sigma _{\stackrel{}{r}+\stackrel{}{\alpha }}+\frac{\kappa }{2}\underset{\stackrel{}{r},\stackrel{}{\alpha },\stackrel{}{\beta }}{}\sigma _\stackrel{}{r}\sigma _{\stackrel{}{r}+\stackrel{}{\alpha }+\stackrel{}{\beta }}\frac{1\kappa }{2}\underset{\stackrel{}{r},\stackrel{}{\alpha },\stackrel{}{\beta }}{}\sigma _\stackrel{}{r}\sigma _{\stackrel{}{r}+\stackrel{}{\alpha }}\sigma _{\stackrel{}{r}+\stackrel{}{\alpha }+\stackrel{}{\beta }}\sigma _{\stackrel{}{r}+\stackrel{}{\beta }},$$
(1)
where $`\stackrel{}{\alpha }`$ and $`\stackrel{}{\beta }`$ are lattice unit vectors. The model is defined on an cubic lattice. The system exhibits a very high degree of symmetry due to the particular ratio of the couplings. This symmetry implies that there is no cost by the flipping of any plane of spins. This results in a highly degenerate ground state. Note in passing that this last feature is common in every glassy system.

Some numerical evidence regarding dynamical and equilibrium properties of the hamiltonian (1) has been accumulated in the past, particularly for $`\kappa =0`$, but also for $`\kappa 0`$. These studies revealed many interesting features and provided evidence that a glassy phase is present in the phase diagram For the interesting case of the four–spin model where randomly distributed couplings are considered, see ().. First of all, it is known that the model has a first order phase transition at a temperature $`T_c`$ where the solid to liquid transition is present . There is also good evidence that a dynamical transition exists at $`T_g<T_c`$ which seems to mark the onset of the glassy behaviour . The existence of $`T_g`$ and the study of the dynamical properties of the system above and below that value have been considered in -. Preliminary results for $`\kappa 0`$ were given in  in what respects the equilibrium properties, and in  in what respects the dynamical properties of the system .

In section 2 we provide a qualitative view of the slow dynamics behaviour of the cold phase by considering the relaxation properties of the model and we estimate with good accuracy $`T_g`$ by measuring the spin-spin autocorrelation functions. We show that there is a dramatic change in the behaviour of this function above and below $`T_g`$. In section 3 we study a different case of the spin gonihedric action by taking $`\kappa =1`$ in (1), apparently a much simpler system (the plaquette term in the hamiltonian is absent for this value). For that value of the parameter $`\kappa `$ the system is described by nearest and next-to-nearest interactions. We confirm, by studying the energy susceptibility, that the system has a second order phase transition Note that from some preliminiary results referred in () the first order transition present for $`\kappa =0`$ gets weaker and possibly becomes second order at $`\kappa 0.5`$. . Furthermore we find slow dynamics behaviour anywhere in the cold phase (below the critical temperature). By studying spin–spin autocorrelation function and the overlap spin distribution function, we provide evidence that the system exhibits type II aging , which is a feature of the glassy systems.

The physical interest of this model is twofold. On the formal side, the model is of such simplicity that a theoretical understanding of the mechanisms underlying slow dynamics and glassy behaviour appears possible. On a more practical side, it would be extremely interesting to be able to understand and produce magnetic materials and coatings with such finely tuned (or approximately so) couplings. The extremely long relaxation times would make them very robust agains thermal noise and fluctuations, yet encoding information there would be as simple as in a normal magnetic material. This possibility has been suggested in . While the plaquette term seems hard to imitate in real materials, the fact that many of the interesting features persist for $`\kappa =1`$ makes perhaps such possibility less remote.

## 2 Four–spin interaction ($`\kappa =0`$)

In this section we will study the case $`\kappa =0`$, i.e, a spin model with only four–spin (plaquette) interaction. In this case the hamiltonian (1) takes the form:

$$H=\frac{1}{2}\underset{\stackrel{}{r},\stackrel{}{\alpha },\stackrel{}{\beta }}{}\sigma _\stackrel{}{r}\sigma _{\stackrel{}{r}+\stackrel{}{\alpha }}\sigma _{\stackrel{}{r}+\stackrel{}{\alpha }+\stackrel{}{\beta }}\sigma _{\stackrel{}{r}+\stackrel{}{\beta }}$$
(2)
This form of the interaction leads to a highly degenerated ground state. Flipping every spin in any plane of the cubic lattice implies invariant ground state energy. Taking this symmetry into account, the degenerancy of the ground state is equal to $`2^{3L}`$ due to the $`3L`$ diferent planes in a cubic lattice. This degeneracy survives even at $`T0`$.

It is well known that this model exhibits a first order phase transition at $`T_c=1.95`$ along with a dynamical transition at $`T_g1.7`$, which is the temperature where the glassy behaviour shows up , , . Even though our main interest is the study of the glassy characteristics by looking at the relaxation as well as the autocorrelation of the order parameters (to be defined below) in the glassy phase, the region $`T_g<T<T_c`$ is interesting as well. In this region, numerical simulations clearly indicate the presence of metastability. This is exemplified by the result presented in Fig. 1a where the time evolution of an $`8^3`$ volume at a temperature value just below $`T_c`$ is shown \**\**\**This simulation refers to the system expressed in dual variables  and it has been carried out using a cluster algorithm.. This figure also gives clear evidence that there is a first order phase transition. Fig. 1b consists of two different curves. One of them corresponds to a heating process starting from an initial ordered congiguration. The other one describes the result of quenching a random initial configuration for each one of the temperature values shown. In both cases the simulation has been carried out using the Metropolis algorithm on a $`24^3`$ volume by performing $`10^4`$ measurements at each temperature value. As we will see in a while, the approach to equilibrium in the region $`T_g<T<T_c`$ is non-standard and it is well described by a stretched exponential, instead of a simple exponential. For $`T<T_g`$ one immediately sees that the results from the quench of the random configuration after $`10^4`$ MC steps differ from those obtained starting from an ordered configuration (any of the $`2^{3L}`$ vacua) . The difference appears to be constant all over this region (for a fixed number of thermalization steps). This clearly hints to the coexistence of two different dynamics. Initially fast dynamics quickly brings an initial configuration which is badly out of equilibrium to some sort of approximate equilibrium. At that point slow dynamics takes over and the evolution of the system is considerably freezed.

Before getting into the more quantitative aspects of these results, it is perhaps interesting to turn to one of our motivations, namely to test whether the appearence of slow dynamics makes the transition between two approximate ground states so slow as to make a given configuration virtually indestructible by thermal fluctuations, thus providing a convenient way of storing information.

To this end, we simulate the system on a cubic lattice and we use a Glauber algorithm , which is assumed to provide a good approximation to the thermal mechanism of fluctuations. Metropolis or heat bath give very similar results. We look at the decay of an artificial initial configuration consisted of an inner volume with a chessboard-like arrangement of the spins, which is one of the ground states, while the outer volume spins are fixed at $`+1`$ i.e, they form a ferromagnetic ground state. We use fixed boundary conditions in order to make the system decay to a ferromagnetic ground state. Our measured quantity is the number of the minus spins, $`N_{}(t)`$, in terms of the computing time divided by the number of spins at time $`t=0`$, $`N_{}(0)`$. The quantity $`N_{}(t)`$ is clearly related with the magnetization, $`M`$, since: $`MN_++N_{}=N_{tot}2N_{}`$. Two examples of our results are given in Fig. 2a and Fig. 2b which have been produced from simulations over two lattice volumes namely $`16^3`$ and $`30^3`$ which at time $`t=0`$ have enclosed $`8^3`$ and $`20^3`$ lattice volumes respectively with a chessboard–like arrangement for the spins.

The two curves in each subfigure have been generated after averaging over a sample of 50–100 copies where each one evolves starting from the same initial configuration, in order to reduce the noise. Two cases are shown. One is for temperature $`T=1.83`$, in the supercooled phase, while the other one, $`T=1.0`$, lies in the glassy phase . The difference on the relaxation time between the two temperatures is obvious. For the $`T=1.83`$ case the system reaches a stable value relatively fast for both lattice volumes, but this is not the case for $`T=1.0`$. Note that the slope of the curve keeps taking non-zero value even for remarkably large times. This is perhaps clearer in Fig. 3 where long runs are depicted for $`T=1.0`$, confirming that we are in presence of slow dynamics.

This behaviour can be seen in a more apparent way in Fig. 4 where the results of Fig. 3 are presented in a logarithmic time scale. In Fig. 4 the logarithmic decay is present clearly enough.

The very distinct dynamics between the glassy phase and the evolution of the supercooled phase in the metastability region $`T_g<T<T_c`$ can be seen in a more definite way in Fig. 5. In this plot the results for four temperatures are presented for the bigger volume used i.e. $`30^3`$. All of them correspond to a random starting configuration. For $`T=1.83`$, which lies in the metastability region the evolution seems very fast all the way to the equilibrium value. On the contrary for the other three temperatures, a fast evolution is initially observed followed by a very slow one which persists up to very long times.

These results are very suggestive and indeed show that the system finds very difficult to overcome dynamical energy barriers that are created along the evolution and this is undoubtedly the reason for the slow dynamics. Recall that the energy of the model is concentrated on the edges; the system has vanishing microscopic surface tension. To reduce the volume of the excitation with local moves, the total edge length must temporarily increase by a substantial amount. That makes excitations such as the one we have been analysing virtually stable.

In spite of the unambiguity of the previous results, in order to determine some properties associated to the very slow dynamics observed at low temperatures we shall proceed to analysing several dynamical correlators.

We shall first study the spin–spin autocorrelation function for temperatures lying in the metastability region after a random start, thus forcing the system to be in a supercooled phase. Its definition is given by

$$C_{spin}(t,t_w)=\frac{1}{N}\underset{\stackrel{}{r}}{}\sigma _\stackrel{}{r}(t_w)\sigma _\stackrel{}{r}(t+t_w)$$
(3)
The brackets mean that we take the average value over copies starting from a random configuration (200-400 in our case). The waiting time $`t_w`$ is the time for which the system is being thermalized before taking the measurements at subsequent times denoted by $`t`$. The waiting time $`t_w`$ is taken to be about 300 by noticing that from that value on the resultant values for $`C_{spin}`$ are quite identical as long as the temperature lies in the metastability region. In Fig. 6 we present some of our results for $`C_{spin}`$ in a $`40^3`$ lattice volume and for four temperature values. The fittings are stretched exponentials of the form $`ae^{(t/\tau )^b}`$. We denote by $`\tau `$ the relaxation time. In all cases we found $`0.60<b<0.80`$ with an error smaller than 0.004. Since the fits seem quite good we go on and make a plot of the resulting values of $`\tau `$ versus the temperature $`T`$. The result is shown in Fig. 7 where the corresponding fit to a function having the form, $`\text{const}/(TT_g)^c`$, is quite good and leads to the prediction $`T_g=1.698(1)`$ with $`c=0.41`$. At $`T_g`$ the autocorrelation time is expected to diverge because of the onset of the slow dynamics which turns the stretched exponential behaviour into a power law (with a small exponent) or a logarithm.

Although that result is in good agreement with previous simulations ,  the method relying on the stretched exponential fit may prove to be too risky for an exact prediction of $`T_g`$ due to ambiguities in the fitting process and perhaps is not trustable to that accuracy. As an alternative and a cross check, for two different temperature values, namely $`T=1.720`$ and $`T=1.695`$ we show the behaviour of the $`C_{spin}`$ for two very different values of the waiting time, $`t_w=300,2000`$. One sees that for the higher temperature the behaviour of the $`C_{spin}`$ is identical but for the smaller one there is a strong dependence on $`t_w`$ showing that for larger waiting times the system exhibits different, and in particular, slow dynamics behaviour which is a strong signal that the system has passed to the glassy phase. Hence we are able to estimate the value of $`T_g`$ to lie in between the interval $`1.695<T_g<1.720`$ which indeed agrees also with the previously obtained value.

## 3 Next and next-to-nearest neighbour interaction ($`\kappa =1`$)

As we have indicated the gonihedric spin model is actually a family of models which are parametrized by a real quantity $`\kappa `$. All the members of this family share the common feature of having their spin interfaces weighed with the total edge length. The value of $`\kappa `$ simply indicates the degree of self-avoidance of the surfaces. The plaquette model whose dynamical properties we just discussed corresponds to $`\kappa =0`$. An interesting member of this family is given by $`\kappa =1`$. The degeneracy of the ground state in this case is $`32^L`$.

In this case the Hamiltonian is given by

$$H=2\underset{\stackrel{}{r},\stackrel{}{\alpha }}{}\sigma _\stackrel{}{r}\sigma _{\stackrel{}{r}+\stackrel{}{\alpha }}+\frac{1}{2}\underset{\stackrel{}{r},\stackrel{}{\alpha },\stackrel{}{\beta }}{}\sigma _\stackrel{}{r}\sigma _{\stackrel{}{r}+\stackrel{}{\alpha }+\stackrel{}{\beta }}$$
(4)
As we see, the spin plaquette term has disappeared. From a practical point of view this model may be particularly interesting as the plaquette term is obviously hard to get in real materials. From the standpoint of a spin system this is just a model with nearest and next-to-nearest neighbour interactions (though a finely tuned one).

We shall proceed to studying the dynamical properties of the system by using Monte–Carlo methods. In performing the simulations we used the Metropolis algorithm for several lattice volumes, namely $`10^3`$, $`16^3`$, $`20^3`$, $`24^3`$, $`30^3`$, $`40^3`$ and $`46^3`$. We imposed either fixed or periodic boundary conditions depending on the kind of the measurement that was performed. During the presentation of our results it will be explicitly mentioned when we used one or the other.

We begin by giving the behaviour for the susceptibility of the energy with the lattice volume which is given by

$$S(E)=V(E^2E^2)$$

We denote by $`V`$ the lattice volume and by the symbol $``$, the average value over sweeps. In Fig. 9 the susceptibility for the system energy as a function of the temperature for four lattice volumes is depicted. For every point in the figure we performed $`10^5`$ thermalization sweeps followed by more than $`10^5`$ measurements, using periodic boundary conditions and starting from the ordered configuration. The peak for each volume clearly increases with it although with an exponent less than one, which is a signal for a second order phase transition. Also, the positions of the peaks for the bigger volumes are seen to concetrate around the value $`T_c=2.329`$ which is the pseudocritical value for the up to the volumes we used. It should be mentioned that while the second order character of the phase transition was already mentioned in ,  however our prediction for $`T_c`$ gives a somewhat smaller value.

Having estimated the critical temperature value for the system we want now to study the dynamic behaviour at low temperatures, i.e. below $`T_c`$. To this end we use first the same method as for the $`\kappa =0`$ case, described in the previous section, which consists of choosing an initial chessboard–like configuration defined in a cubic subvolume of the system. We impose fixed boundary conditions and we study the behaviour of the dynamical quantity

$$\frac{E(t)E_{eq}}{E(0)E_{eq}}$$
(5)
with time for various temperatures in the cold phase and for different sizes of the inner chessboard–like volume. In Eq.(5) $`E(t)`$ and $`E(0)`$ denote the energy of the system at time $`t`$ and just after the first sweep performed, respectively and $`E_{eq}`$ is the system energy when an ordered configuration is taken as the initial one. In Fig. 10a we give an example for lattice volume $`16^3`$. We fix the temperature at the value $`T=1.1`$ and we are interested in the behaviour as the dimension of the inner chessboard–like volume increases from $`L=2`$ to $`L=8`$. Each curve has been produced after averaging several hundreds of repetitions starting from the same initial configuration in order for the noise to be reduced. The results show a dramatic increase of the relaxation time with increasing system size. For the maximum value of $`L`$ used, the observable (5) shows very slow variation with the time and it follows a power law behaviour with estimated exponent equal to 0.18, so a logarithmic behaviour cannot be really excluded.

We can reach similar conclusions by studying the above relaxation quantity for decreasing values of temperature at fixed $`L=3`$. In Fig. 10b the corresponding results are depicted for a $`16^3`$ total volume. In particular for the value $`T=0.2`$ the dynamics is so slow that it seems really hard for the system to get to the ground state for accessible computing times.

Next we shall move to a more quantitative study following the same lines of the previous section. We look for aging features by considering the autocorrelation spin–spin function defined by (3).

Brackets indicate that we take the average value over 100-200 trials starting from a random configuration. $`t_w`$ is the time for which the system is being thermalized before taking the measurements at subsequent times denoted by $`t`$. We follow the same lines of analysis as in  which concerns the $`\kappa =0`$ case. In Fig. 11 we show an example of the behaviour of the autocorrelation function (3) for the temperature value T=2.2 and for a $`24^3`$ lattice volume. Six different curves for the $`C_{spin}`$ are depicted each corresponding to a different value of $`t_w`$. The very slow evolution of the autocorrelation function combined with the clear dependence on the time $`t_w`$ give strong evidence in favour of aging. In order to determine the aging (whether I or II) we consider the overlap function :

$$Q(t+t_w,t+t_w)=\frac{1}{N}\underset{\stackrel{}{r}}{}\sigma _\stackrel{}{r}^{(1)}(t+t_w)\sigma _\stackrel{}{r}^{(2)}(t+t_w)$$
(6)
The measuring process for the above consists of the relaxation of the system for time $`t_w`$ starting from the random configuration. At that moment we take two copies $`\sigma ^{(1)}`$ and $`\sigma ^{(2)}`$ each evolving independently for time $`t`$. In Fig. 12 the behaviour of $`Q(t+t_w,t+t_w)`$ in terms of $`C_{spin}(t,t_w)`$ for the same values for the volume and temperatute as in Fig. 11 and for four values of the waiting time $`t_w`$ is shown. It clearly follows that as $`C_{spin}(t,t_w)0`$, $`Q(t+t_w,t+t_w)0`$. The fact that the two system copies are moving independently since they show zero overlapping indicates strong evidence for type II aging.

## 4 Conclusions

At this point we would like to summarize our main conclusions.

We have carried out a rather complete analysis of the dynamical properties of the gonihedric spin model. The model naturally exhibits large potential barriers that are dynamically generated along the evolution process when trying to adjust to the environment. This is due to the fact that the area term of the interfaces is completely subleading and the dynamics is driven by the total edge length of the excitations. This creates basins of stability from where it is virtually impossible to kick out the system by local thermal fluctuations.

Several algorithms have been used: Glauber dynamics, heat bath and standard Metropolis. Modulo an overall rescaling of Monte–Carlo time, the results are fully equivalent.

In this paper we have studied numerically two cases of the gonihedric spin model. For the $`\kappa =0`$ case we give evidence of a glassy transition and we estimate the temperature value $`T_g`$ at which the glassy phase arises being in the interval $`[1.695,1.720]`$. We show that the dynamics is most likely logarithmic for $`T<T_g`$. In the mestatability region, above $`T_g`$ but below the true thermodynamical transition $`T_c`$ (where the system has a first order transition) the approach to equilibrium is well described by a stretched exponential, whose characteristic time scale diverges as $`TT_g^+`$.

We confirm that for the case $`\kappa =1`$ the system exhibits a second order phase transition. The critical properties of this transition have not been elucidated yet; a tentative determination of the exponents using finite size scaling carried out in  gives a rather non-standard set of values. Work on this is in progress, but the problem is agravated by the presence of slow dynamics. Indeed we have detected in the cold phase of this system very slow dynamical behaviour, quite similar to what takes place with $`\kappa =0`$ for $`T<T_g`$. This dynamical behaviour has all the features of type II aging and is again compatible with logarithmic evolution.

Possible applications of these simple, but interesting, spin models are discussed.

Acknowledgements

D.E. thanks D. Johnston and A. Lipowski for discussions and the hospitality of the Department of Mathematics of Heriot–Watt University. P.D. is grateful to G. Koutsoumbas for useful discussions and acknowledges support from ”EUROGRID– Discrete random geometries: from solid state physics to quantum gravity” (HPRN–CT–1999–00161). The support from grant MCyT FPA 2001–3598 and CIRIT 200156R–00065 is also acknowledged. Some of this work was carried out at CESCA.


text:
# 1 The Situation - Open Questions after HERA2

## 1 The Situation - Open Questions after HERA2

HERA is currently scheduled to run until the end of 2006, the aim being to collect polarised $`\stackrel{}{e}p`$ scattering data corresponding to an integrated luminosity of 1 fb-1 at the highest possible energy and to obtain some data at lower energy. This programme will allow exploration of proton structure at the smallest possible dimensions and will also produce large amounts of precision data at low Bjorken $`x`$. These data may result in the discovery of new effects at low $`x`$ or at high momentum transfers. In the latter case, the comparison of precision measurements with predictions from the QCD evolution equations, well understood in this kinematic regime, leads to sensitivity to the effects of large extra dimensions or contact interactions in the multi-TeV range. The production of new heavy particles, such as the squarks of $`R`$-parity violating supersymmetric theories, leptoquarks or excited fermions is also possible up to masses of about $`320`$ GeV and the effects of these particles can be detected for masses significantly above this. The observation of novel effects with the envisaged five-fold increase in $`e^+p`$ and thirty-fold increase in $`e^{}p`$ integrated luminosity would of course make further investigations beyond the year 2006 mandatory. Even if no new effects are discovered, the current HERA programme will leave a number of fundamental questions unanswered. Issues of prime importance are:

* What is the partonic structure of the neutron at low $`x`$, and at large $`Q^2`$ and $`x`$?
* What is the origin of confinement in long distance (small $`x`$) and short distance (large $`x`$) processes?
* How does the total deep inelastic scattering cross-section for lepton-proton and lepton-nucleus interactions change in the high density regime where cross sections may saturate and how are these changes reflected in the hadronic final state?
* How do the partons conspire to ensure that the spin of the nucleon is 1/2?

Understanding these problems will require data from lepton-deuteron and lepton-heavy ion scattering, as well as a programme of high energy and high luminosity colliding beam and fixed target experiments involving polarised leptons, polarised protons and polarised deuterons. In addition to addressing the above questions, this programme will be vital to our understanding of the astrophysical significance of extremely high-energy neutrino scattering, as well as being necessary to ensure full exploitation of the physics potential of the LHC and of heavy ion collisions.

The Durham meeting brought together about 40 physicists, including the Spokespersons, Physics and Technical Coordinators of the HERA experiments H1, ZEUS, and HERMES, of COMPASS at CERN, leaders of the Electron-Ion Collider (EIC) Community and leading machine and theory experts. The programme of the workshop (available at http://hep.ph.liv.ac.uk/$``$green/HERAfuture) comprised sessions devoted to the current HERA programme (HERA2), QCD, polarised $`\mathrm{}N`$ physics, detectors, deuterons and heavier nuclei and machine developments at HERA and the EIC. It ended with an extended discussion about the plans, options and prospects for future deep-inelastic scattering experiments prior to the operation of TESLA.

## 2 Physics Subjects - Nuclei and Polarisation

Much of the physics of nuclei and polarisation has been studied at previous workshops, and frequent reference was made to the HERA $`eA`$ workshops , to the HERA spin workshops , to the EIC white book , to the THERA book , to the TESLA-N  and to the ELFE  proposals.

As mentioned above, after completion of the currently approved HERA programme, there will be a continued interest in precision measurements in the low $`x`$ region. This concerns in particular the measurement of jets very close to the proton beam direction, crucial to the understanding of the emission of gluons at low $`x`$, and the kinematic region in which $`Q^21`$ GeV2, where the energy dependence of the $`\gamma ^{}p`$ cross-section becomes hadron-hadron like. Precision inclusive and exclusive cross-section measurements in this kinematic region may yield insight into saturation physics and the confinement problem. This programme will require that the instrumentation in the forward and backward regions very close to the beam pipe at HERA be upgraded, that efficient proton tagging is possible and that the beam divergence is small.

Electron-deuteron scattering appears to be the natural next step at HERA following the $`ep`$ programme. Deuterons are a source of quasi-free neutrons, measurements of which are perturbed slightly in the region $`x<0.1`$ by nuclear shadowing effects. These latter are related to the diffractive parton densities and hence the effects of shadowing on lepton-deuteron cross-sections can be determined to an accuracy of better than 1-2%. The quark distribution asymmetry, $`(u+\overline{u})(d+\overline{d})`$, can thus be accurately measured from the difference $`F_2^pF_2^n`$. Deuteron data are also essential for determining the individual flavour decomposition of the nucleon parton distributions at large $`x`$, allowing measurement of quantities such as $`sc`$ or the $`d/u`$ ratio, as well as for determining charged current structure functions and for precision tests of $`Q^2`$ evolution in perturbative QCD.

The rise of the proton structure function $`F_2`$ towards low $`x`$ in the deep-inelastic scattering region is due to the high sea quark density in the proton. This is related to the gluon distribution $`xgF_2/\mathrm{ln}Q^2`$, which also rises as $`x^\lambda `$ towards low $`x`$, with $`\lambda 0.1\mathrm{}0.4`$ for $`Q^21\mathrm{}100`$ GeV2. Since $`xg`$ is expected to increase with atomic number $`A`$ as $`A^{1/3}`$ (modulo shadowing effects), electron-nucleus scattering allows an equivalent Bjorken $`x=x_N/(A^{1/3})^{1/\lambda }`$ to be accessed: such collisions can thus be used to investigate a kinematic regime which in $`ep`$ scattering would require a considerable increase in the centre-of-mass energy. Diffractive processes may constitute up to 50% of the inclusive lepton-nucleus cross section, the maximum value allowed by unitarity considerations. Such an observation would represent an unambiguous signal for a new regime in deep-inelastic scattering. Since the contribution of small-size configurations to diffraction goes as the square of $`xg`$, non-linear effects which damp the growth of the total cross-section may first be seen in diffraction. The nucleus is widely considered to be an ideal laboratory for deconfinement and high density QCD analyses.

It is envisaged that spin physics with fixed polarised targets will be continued at DESY beyond 2006 with a high luminosity measurement programme. With modest improvements to the HERMES apparatus, high resolution experiments with increased luminosity will become possible. Such measurements will be devoted to a detailed study of exclusive reactions, in particular deeply virtual Compton scattering, and will lead to the first precise experimental information on generalised parton distributions.

Collisions of beams of polarised electrons and polarised nucleons will make possible the first investigations of spin phenomena at high $`Q^2`$ and at low $`x`$. These will be enhanced by the comprehensive reconstruction of the final state possible in colliding beam experiments. The behaviour of the spin structure function $`g_1(x,Q^2)`$ at low $`x`$ is unknown, but it is expected to change even more dramatically than $`F_2`$. The $`Q^2`$ dependence of $`g_1`$ determines the gluon spin distribution $`\mathrm{\Delta }G`$, one of the components of the spin of the proton. Semi-inclusive measurements can be used to explore the flavour structure of spin: dijets give access to $`\mathrm{\Delta }G`$ and deeply virtual Compton scattering to the generalised parton distributions. Photoproduction processes give insight into the polarised gluonic structure of the photon and diffractive polarised scattering into pomeron exchange. Beyond the fixed target $`\stackrel{}{e}\stackrel{}{p}`$ domain, there is a vast unexplored kinematic region and polarised colliding beam experiments may radically change our view on the spin structure of the nucleon.

## 3 Detector Aspects

The status and upgrade plans for H1, ZEUS and HERMES at DESY, and COMPASS at CERN were summarised by the technical coordinators of these experiments. It was pointed out that in 2006 the HERA collider detectors will be 15 years old. Continued operation will require some detector and electronics upgrades. For example, whereas the H1 liquid argon calorimeter has proved to be extremely stable, in all probability an upgrade to the central tracking facilities will be necessary, including replacement and possibly extension of the silicon detectors. More accurate judgements on the long-term prospects for operation of the H1 and ZEUS detectors will be possible following the first years of high luminosity running at HERA2.

Currently, the HERA interaction regions are optimised for high luminosity running, with focussing magnets placed close to the interaction region. This sets a limit on the minimum measurable four-momentum transfer squared of about 4 GeV2. These high luminosities are essential for the measurement of small spin-induced asymmetries and also for high statistics $`eD`$ measurements. At low $`Q^2`$, event rates will be high and substantial upgrades of the front-end electronics may be required.

For accessing the region around $`Q^2=1`$ GeV2 and the smallest $`x`$ values in $`eA`$ scattering, it will be necessary to return to an arrangement which leaves room for new or upgraded detectors close to the beam pipe. An essential requirement for diffractive measurements and for deuteron spectator tagging is large acceptance forward nucleon tagging.

It is expected that the HERMES detector, perhaps with some upgrades to the data acquisition system and silicon detectors, will be able to run beyond 2006. The COMPASS experiment uses a variety of modern detector technologies which may be of use for future HERA3 detectors.

## 4 Accelerators - HERA and the EIC

Following the recent upgrade, HERA is scheduled to run with a luminosity of $`7\times 10^{31}`$ cm-2s-1 until the end of 2006. The maximum attainable luminosity is estimated to be $`13\times 10^{31}`$ cm-2s-1. This corresponds to an annual integrated luminosity of about 500 pb-1, making the measurement of small asymmetries at low $`x`$ in polarised $`\stackrel{}{e}\stackrel{}{p}`$ scattering statistically feasible if intense polarised proton sources can be constructed. Precise measurements require a high degree of polarisation which must be transferred through the accelerator chain, necessitating the use of Siberian Snakes. Polarimeters allowing accurate polarisation measurements are also essential. The small anomalous magnetic moment of the deuteron may make realisation of polarised $`\stackrel{}{e}\stackrel{}{D}`$ collisions easier than the $`\stackrel{}{e}\stackrel{}{p}`$ case using a scheme which relies on resonant driving of the spin using horizontal RF fields.

Pilot $`eA`$ scattering studies can be performed using deuterium, oxygen and calcium. Light ions can be accelerated in HERA with moderate modifications to the accelerator chain. Electron cooling becomes necessary if heavier nuclei are to be accelerated to counter intra-beam scattering (IBS) effects. The luminosity is then expected to scale as $`L_AL_p/A`$. For deuterons, the IBS time exceeds two hours and a luminosity of $`3.5\times 10^{31}`$ cm-2s-1 may be achieved without cooling. Thus, high luminosity $`eD`$ running can be realised immediately after completion of the HERA2 programme and requires no major modifications to HERA beyond the necessary changes to the source.

The EIC, formal proposals for which are expected to be presented in about 2005, will be an intense polarised electron-ion collider using electron beam cooled ions in RHIC (of up to $`E_p=250`$GeV proton and $`E_{Au}=100`$GeV/A gold energy) colliding with electrons from a ring or linear accelerator of about 10 GeV maximum energy. Thus the EIC has 10 times less centre-of-mass energy than HERA, but due to its high luminosity it greatly extends the range of current polarised fixed target experiments, and with heavy nuclei allows access to high parton densities. The ring-ring accelerator has an estimated luminosity of $`25(0.7)\times 10^{31}`$ cm-2s-1 for protons (gold), which is comparable to HERA after run 2. With an energy recovery linac, it is expected that a still higher luminosity of $`100(1.0)\times 10^{31}`$ cm-2s-1 for protons (gold) may be obtained. The linac, while accelerating $`e^{}`$ only, has the further advantage of providing high polarisation and avoiding synchrotron radiation background.

## 5 Conclusions and Outlook

If answers are to be found to the questions posed in the introduction, three experimental programmes must be pursued:

* Electron-deuteron scattering experiments, with both high luminosity (more than 100 pb-1) and detectors that have dedicated low $`x`$ and forward nucleon tagging capability.
* Extended measurements at low $`x`$: studies of the problems of confinement require precision diffractive and non-diffractive data at $`Q^2`$ near to 1 GeV2; high density gluon effects leading to saturation may be studied in electron-nucleus scattering.
* High luminosity polarised $`\stackrel{}{e}\stackrel{}{p}`$ and $`\stackrel{}{e}\stackrel{}{D}`$ scattering: spin physics requires further operation of fixed target experiments and the investigation of the unknown world of low $`x`$ and high $`Q^2`$ deep-inelastic spin phenomena using collider detectors.

Work in the coming years will be directed towards producing proposals for the operation of HERA beyond 2006 and for the EIC which may start operation in 2012, after completion of the polarised proton phase of RHIC. A programme will be devised which produces maximum physics return from HERA and the EIC, with their complementary reach in energy, luminosity and nuclear mass. An extension by another order of magnitude in energy in $`ep`$ scattering can be achieved with THERA, the future $`ep`$ collider operating at TeV energies.

## 6 Acknowledgements

The authors would like to thank the DESY Directorate and the Institute for Particle Physics Phenomenology for their support and Linda Wilkinson and the other members of the organising committee for helping to arrange what proved to be a most interesting and productive meeting, namely Jochen Bartels, Alan Martin, Martin McDermott, Klaus Rith, Jim Whitmore, Ferdinand Willeke and Giulia Zanderighi. This summary has benefited from the suggestions of Allen Caldwell, Eckhard Elsen, Robert Klanner, Alan Martin, Martin McDermott, Dirk Ryckbosch, Peter Schleper, Mark Strikman and Ferdinand Willeke. We would also like to thank the representatives of the theory community and the COMPASS, HERMES, H1 and ZEUS Collaborations who contributed to the workshop.