text:
# Multifragmentation - what the data tell us about the different models

## 1 Models for describing multifragmentation

About 15 years ago it has been found that in nucleus-nucleus collisions at intermediate energies up to 15 intermediate mass fragments (IMF’s) with $`Z3`$ are created. Since then many dedicated experiments have been performed to ”nail down” the reaction mechanism but despite of all experimental and theoretical efforts we have not found an answer yet. The different contributions in these proceedings present a clear evidence for this fact. Over the years two different conjectures have been launched and developed, improved and further improved without giving a conclusive answer. The reason is, as we will see, that the two most advanced but also almost opposite approaches predict the same behavior for several key observables. Before we discuss this in detail we will introduce shortly the two principal approaches.

### 1.1 Multifragmentation as a statistical process

One may assume that multifragmentation is a statistical process. This means that during the reaction at least a subsystem comes to statistical equilibrium and maintains this equilibrium during the expansion until it has reached the freeze out volume. Then it disintegrates INSTANTANEOUSLY into neutrons, light charged particles (LCP’s) and IMF’s. Because it is assumed that at freeze out the system is in statistical equilibrium the disintegration pattern is determined by phase space. This means that one searches for all microstates (consisting of fragments and nucleons with a given kinetic and excitation energy) which are compatible with a total energy, a total proton and neutron number and a given volume. Each microstate has the same weight and hence one samples all microstates in order to obtain $`\frac{d^4\sigma }{d^3pdA}`$. There are three different models ( \- ) which perform this task. Sequential decay (i.e. the formation of a compound nucleus which decays by subsequent fragment emission) has been ruled out because it gives too few fragments.

Unfortunately this task is not as easy as it seems to be, both conceptually and numerically. Therefore the results of these models differ considerable despite of the fact that one is tempted to believe that counting the number of microstates is a well defined task.

Conceptually there are two problems:
a) there is no theory which tells us how to treat the unstable states. Usually the number of excited states of fragments is calculated via a level density parameter but some of these excited states are unstable against particle (usually neutron) emission. Therefore it is not evident, whether they should be treated as an excited fragment A or as a fragment A-1 plus one neutron. Both descriptions give rise to different microstates. The different answers to this problem yield a quite different neutron yield as has been recently found out by Toeke and Schroeder .
b) the freeze out volume can be defined quite differently. Whereas in ref  a constant freeze out volume is assumed in ref.  the freeze out volume depends on the microstate. There is no convincing argument why the one assumption should be better than the other.

Also numerically (especially concerning the Coulomb energy) the approaches differ considerably. It is the merit of Gross and Sneppen  to have discussed all the differences in detail and to have demonstrated how it comes that two seemingly identical approaches produce quite different results.

### 1.2 Multifragmentation as a dynamical process

On the other hand there are dynamical models. Originally developed to describe proton spectra, particle production and the influence of the nuclear equation of state on observables like the flow, they have later been further developed to describe multifragmentation as well. In 1985 the Quantum Molecular Dynamics model  has been introduced which describes the time evolution of the n-body Wignerdensity and hence allows for the investigation of fragments which are n-body correlations in this context. In this approach the nucleons are presented as Gaussian wave functions

$$\varphi _\alpha (p_1,t)=(\frac{L}{2\pi })^{3/4}e^{(p_1p_\alpha )^2L/4}e^{ix_\alpha p_1}e^{2ip_1p_\alpha t/2m+ip_\alpha ^2t/2m}$$
(1)
where $`\alpha `$ characterizes the nucleon. The time evolution of the individual nucleons is determined by a variational principle. At the very end a minimum spanning tree is applied to identify the clusters. In these models the whole reaction, starting from the initially separated projectile and target nuclei until the finally observable fragments and single particles is simulated. Besides the Hamiltonian which is, however, well determined by the requirement to have stable nuclei with the right binding energy, the only parameters to change are the impact parameter and the beam energy. The different final states, created in each individual simulation, are caused by the different possibilities to realize the same single particle distribution and the same binding energy by choosing different positions $`x_\alpha `$ and momenta $`p_\alpha `$ for the individual nucleons. For details of this model we refer to the references .

There are two conceptual criticisms of this model. First, that the wave function of the nucleons is not antisymmetrized. Therefore one argued, making reference to a Fermi gas, that the specific heat will be different as compared to an antisymmetrized system. Hence, if the system comes to equilibrium one expected a quite different fragmentation pattern. Although this is true of course for a Fermi gas, the Aladin collaboration has found  that this is not true for an interacting system containing nucleons and fragments. The above mentioned incertainties of the statistical models have a larger influence on the fragmentation pattern than the different specific heats. Second, being a semiclassical approach, the fragments have the properties of the Weizsäcker mass formula but are not real quantum states. In addition effective charges are used. Hence isotopic distributions cannot be calculated.

Nevertheless with two parameters which have to be adjusted, the width of the wave function and the fragmentation radius of the minimum spanning tree, the QMD model has described many low energy fragmentation data . At these energies the whole nucleus disintegrates into fragments. At higher beam energies, where only the weakly excited spectator matter disintegrates into fragments the model in its original version failed . Only recently the reason for this failure has been identified and the model describes now spectator fragmentation as well .

Careful investigations have shown that in these simulation programs the system passes never a state of equilibrium . In almost all cases where multifragmentation has been observed the systems is, on the contrary, far away from it. This can easily be seen from the videos which have been produced with help if this program .

## 2 Multifragmentation - a theory invariant phenomenon?

In view of the quite different dynamics assumed or obtained in these both approaches to multifragmentation it seems to be astonishing that experiments have not yet decided which of the both conjectures is correct. This fact becomes more understandable if one considers in similar but somewhat simpler models, what one expects to see. For this purpose we assume that the system can be described by canonical thermodynamics, although this is not the case, as we will see later. For the question discussed here it serves the purpose.

First we will suppose a fast fragmentation process as it is the case for a dynamical process. If the disintegration of the nucleus is instantaneous each nucleon keeps its momentum due to the Fermi motion and one expects an average fragment kinetic energy of $`3/5E_F`$ independent of the fragment size , where $`E_F`$ is the Fermi energy. The same independence one expects if the fragments are formed very late, after the system has been expanded while maintaining thermal equilibrium as assumed in the statistical model. This requires the opposite, that the disintegration is very slow. Here the average kinetic energy is 3/2 T, where T is the temperature at freeze out. Hence it is only the absolute value but not the mass dependence of the mean kinetic energy which is different.

The same ambiguity one finds for the mass yield. Fischer  has pointed out that finite systems close to a phase transition show a fragment distribution $`\sigma (A)A^\tau `$. But this power law dependence appears in many physical processes which are not thermal at all . Figure 1 shows schematically this behavior.

This ambiguity remains valid if one employs the more sophisticated models discussed above. For our studies we use data and simulations for central collisions of the reaction 50 AMeV Xe + Sn, measured by the INDRA collaboration . This is the experiment where the most precise data are available at the moment. The SMM calculation we are discussion about are adjusted to reproduce the subset of events fulfilling the condition of completeness (more the 80the total charge and the total momentum are detected) and they show a flow angle $`\theta _flow>60^o`$ . When comparing with a dynamical calculations, we will use the subset of events where the total transverse energy of light charged particles $`E_{trans}>450MeV`$. Both models, SMM and QMD, have been extensively used to interpret this reaction -. We would like to mention that an analysis of the QMD calculations shows that in this approach the system never passes through a state of thermal equilibrium  whereas the application of SMM is only justified if the fragments are formed from an equilibrized subsystem. Thus the reaction scenarios in these models are orthogonal. For the comparison with the data both calculations have been filtered with the experimental acceptance and a centrality cut is employed . The increase of the average transverse energy of the fragments as a function of their mass, observed in the experiment, is larger than predicted in SMM  calculations. Therefore the statistical model calculations have been modified by adding a fourth system parameter, a collective energy, which is parameterized as $`E_{coll}=cA`$, where c is a parameter which remains to be determined and A is the fragment mass. The best agreement the subset of experimental data fulfilling the centrality condition $`\theta _flow>60^o`$ and SMM calculations is obtained with the following set of input parameter:

| freeze out density: | $`1/3\rho _0`$ |
| --- | --- |
| source size: | ZS=78 AS = 186 |
| excitation energy: | Ethermal=7A MeV Ecoll=2.2A MeV |

Even central collisions at intermediate energies have a less or more binary character and consequently emission close to beam or target velocity spoils the spectra of particle emission from a possible thermal source at rest in the center of mass. Therefore a meaningful comparison between statistical model calculations and experimental data is only possible around $`\theta _{CM}=90^o`$. We subdivide therefore the experimental data and the simulations into two equal size $`2\pi `$ intervals:

$`B_{obs}`$ : $`60^o\theta _{CM}120^o`$

$`B_{unobs}`$ : $`\theta _{CM}<60^o,\theta _{CM}>120^o`$.

In $`B_{obs}`$ we observe a flat angular distribution and a constant average energy of IMF’s and LCP’s  as a function of the emission angle, both being prerequisites for a statistical equilibrium. Applying the statistical models, one assumes that in $`B_{unobs}`$ a preequilibrium component is superimposed to the thermal component. Therefore the statistical model cannot be compared to data there. In the dynamical model there is of course nothing like a distinction between an equilibrium and a preequilibrium part therefore we can compare the results with the data. As fig. 2 shows, QMD describes there the fragment spectra quite well .

Now we come to the comparison in $`B_{obs}`$. In fig. 3 we display on top the average kinetic energy of the fragments. On the bottom we show the charge yield distribution. We display the results for QMD and SMM calculations in comparison with the INDRA data. As one can see, these are well reproduced in both theories, underlining the above mentioned observation that these observables are not sensitive to the reaction mechanism.

Hence we are confronted with the fact that completely different underlying reaction dynamics produces the same key observables for the fragments seen in $`60^o\theta _{CM}120^o`$. Therefore more complicated observables have to be employed to distinguish between the different possible reaction mechanisms. The next more complicated is $`\frac{d^2\sigma }{dEdA}`$.

### 2.1 Fragment energy spectra

In fig.4, left hand side, we compare the experimental fragment kinetic energy spectra for emission at $`60^o\theta _{CM}120^o`$ with the SMM calculation for three different fragments. On the right hand side a surprisal analysis is presented to show the structure of the deviation between the experimental and theoretical energy distribution.

We see that the increase of the average fragment kinetic energy as a function of the fragment mass - seen is fig. 3 \- has a quite different reason in the SMM calculations as compared to the experiment. In SMM the maximum of the energy distribution increases due to the flow but the high energy slope remains almost constant. This is expected because the disintegration due to phase space produces the same slope for all particles. In the INDRA data the increase of the average kinetic energy is due to an increase of the high energy slope of the spectra. Because the emission of fragments at midrapidity is a rare process which is logarithmically suppressed for large charges we do not have the statistics in QMD to compare in detail the kinetic energy spectra of fragments with a charge larger than 12. For higher charges the fluctuations render the analysis meaningless, for the slopes as well as for the spectra. The spectra for selected charges are presented in fig. 4. We see a rough agreement of the high energy slope but there are too many fragments at low energy. This effect is understood  but there is presently no remedy to it. It is caused by the nuclear interaction range. It is too large as compared to the standard values and screens therefore a part of the repulsive Coulomb interaction.

Hence one has to conclude that even the most complete experiment, available up to now, is not sufficiently precise to distinguish in a convincing way between the different theories. The more detailed the comparison is made the more the shortcomings of the approaches become apparent. As seen, precise high energy fragment data may discard the statistical models because the high energy part of the spectra has to have there a mass independent slope. To measured them is, however, a very cumbersome task. Therefore it is better to look for a more intelligent solution. A possible way to overcome this deadlock we will discuss in the next section.

In general the statistical model calculations provide excellent agreement with data as far as multiplicity distributions, correlations between the fragment masses and other quantities, which have nothing to do with the kinetic energy of the fragments, are concerned. They fail usually if they are compared with kinetic energy distributions, which show an apparent slope temperature of about 15 MeV. Such an high apparent temperature can only be created if the excitation energy/N is higher than the binding energy/N and hence if the fragments are not stable anymore. Hence at these excitation energies the statistical codes produce configurations consisting LCP’s only.

## 3 Isospin tracing- a possible method to overcome the present deadlock

The above discussion has shown that it is all but easy to use the fragment distributions, obtained in a single reaction, to determine the reaction mechanism. Therefore it is certainly worthwhile to think about a different approach. The idea to use a combination of different projectile/target combinations to address this questions, has been advanced quite a while ago by the Texas group . Changing the isospin of projectile and/or target, one may investigate the multiplicity of the different fragments and can compare the results with the prediction from a statistical model calculation which assumes that the system comes to equilibrium. The problem is, however, that the variation of the isospin is rather limited and therefore the effects are tiny. Therefore it is necessary to find a convincing analysis of the results. Recently the FOPI collaboration came up with an analysis which allows for a direct answer to the question of equilibration . Following earlier work on smaller systems of the Texas group  they used an equal mass projectile/target combination (Ru and Zr) and performed all 4 possible reactions (Zr+Zr, Zr+Ru, Ru+Zr, Ru+Ru). Because the center of mass is the same many detector acceptance problems disappear, and in addition a comparison of Zr+Ru and Ru+Zr allows to address directly the question of equilibration as a function of the rapidity. If the system comes to equilibrium, in both combinations the same rapidity distribution is expected. If, on the contrary, there is a difference in the rapidity distribution, the assumption that the particles are emitted from a statistically equilibrated source is excluded. For the first results of this experiment performed at 400 AMeV we refer to the contribution of Rami.

Here we only repeat the main result and discuss the interpretation in the framework of the models discussed above. In fig.6 we see $`R_Z=\frac{2ZZ^{Zr}Z^{Ru}}{Z^{Zr}Z^{Ru}}`$ as a function of the rapidity for the 4 different projectile/target combinations. Z is the number of protons observed at the given rapidity. $`R_Z`$ gives 1 for Zr+Zr, -1 for Ru+Ru and 0 if equilibrium is reached in the mixed combinations. We see clearly that for $`y^{(o)}=\frac{y}{y_{beam}}>.25`$ the protons do not show a sign of equilibrium. The same has been observed for the ratio of $`He^3`$ and triton . Fig. 6,right, displays the result of the QMD calculations which have been performed by the FOPI-collaboration. For testing the sensitivity of the results they have employed different nucleon nucleons cross sections by rescaling the fitted experimental values by a given factor. We see that the proton rapidity distribution depends quite sensitively on the cross section. Already an increase by a factor of 2 brings the system close to equilibrium. The experimental data are best described by assuming that the cross section in the nuclear medium is close to the free one.

Although this beautiful experiment has shown that at a beam energy of 400 AMeV the protons and the light fragments are not emitted by an equilibrated source we have no definite answer about the mechanism of multifragmentation yet. Therefore we have to continue along this proposition. We have to study the rapidity distribution of fragments not only at 400 MeV but at lower energies as well. But this idea to use the isospin degree of freedom as a tag for equilibration, may free us from the present deadlock and may give a final answer concerning the mechanism of multifragmentation.

## 4 Canonical models-applicable to nuclear physics?

In recent times it has been conjectured that the system formed in heavy ion reactions can be described by canonical thermodynamics as well. This conjecture is the basis for the caloric curve, which displays the temperature as a function of the excitation energy of the system. In case that the system suffers a first order phase transition, this function becomes flat, i.e. an increase in the excitation energy does not corresponds anymore to an increase in temperature because the energy is used as latent heat. The Aladin collaboration has claimed to have seen such a sign for a first order phase transition. The prerequisite for the meaningfulness of such a plot is the possibility to measure the temperature and the excitation energy at the same time. Since the temperature is extracted from the fragment yield it is the temperature of the IMF’s which counts. Here we do not discuss the question how experimentally these quantities can be measured. We are concerned with the question whether for a system as small as that created in nuclear reactions the temperature fluctuations for a given excitation energy are sufficiently small to make a caloric curve possible. In a large system, the fluctuation of the excitation energy for a given temperature is $`\frac{1}{\sqrt{N}}`$ where N is the number of particles in the system. In a large system, with large N these fluctuations are therefore negligible and for a given excitation energy we find a precise value for the temperature of the system.

In order to discuss this question we would like to start with an explanation of the differences between the microcanonical description of a system, as done in the statistical multifragmentation models, and canonical thermodynamics. The basis of all statistical mechanics is the counting of microstates. To make life easier we take a concrete example. Consider an harmonic oscillator whose energy levels i ($`i0`$) can be occupied by particles which have energies $`ϵ_i=\mathrm{}\omega i`$ We have suppressed the zero point motion here. The total energy of the particles is given by $`E=ϵ_in_i`$ where $`n_i`$ is the number of particles which are on the level i. Now we can calculate the number of different repartition which exists for a given number of particles $`N=n_i`$ and a given energy $`E=ϵ_in_i`$. For example, if we have N=2 and $`E=2\mathrm{}\omega `$ we can either have both particles on the level i=1 or one particle on the level i=0 and the other on i=2. Each of these repartitions is called a microstate. $`g(N,E)`$ is the number of microstates for a given E,N and $`\mathrm{\Omega }(N,E)=klng(N,E)`$ is the entropy of the system (k is the Boltzmann constant).

Now, assume that we have two oscillators which can transfer energy and particles and have the total energy $`E=E_1+E_2`$ and the total particle number $`N=N_1+N_2`$. The number of microstates of the combined system is given by

$$g(N,E)=\underset{E_1,N_1}{}g_2(N_1,E_1)g_1(NN_1,EE_1)$$
(2)
where $`g_1(E_1,N_1)`$ and $`g_2(N_2,E_2)`$ are the number of microstates of the oscillators 1 and 2 for given energies and particle numbers $`E_1,N_1,E_2,N_2`$, respectively.

The temperature is defined as $`\frac{klng}{E}`$. For our combined oscillators

$$g(N,E,E_1)=\underset{N_1}{}g_2(N_1,E_1)g_1(NN_1,EE_1)$$
(3)
one finds

$`{\displaystyle \frac{klng(N,E,E_1)}{E_1}}`$ $`=`$ $`{\displaystyle \underset{N_1}{}}{\displaystyle \frac{klng_1(N_1,E_1)}{E_1}}+{\displaystyle \frac{klng_2(NN_1,EE_1)}{E_1}}`$ (4)
$`=`$ $`{\displaystyle \underset{N_1}{}}{\displaystyle \frac{klng_1(N_1,E_1)}{E_1}}{\displaystyle \frac{klng_2(NN_1,E_2)}{E_2}}`$ (5)
$`=`$ $`{\displaystyle \underset{N_1}{}}{\displaystyle \frac{1}{T_1}}{\displaystyle \frac{1}{T_2}}`$ (6)
What can one conclude? Only if $`g(N,E,E_1)`$ is maximal (at $`E_1=E_1^{max}`$) the temperatures in the two subsystems are identical. Otherwise the temperatures are different. Hence for small systems it is the generic case that both subsystems have different temperatures. This is completely counterintuitive because our experience is based on the large systems we encounter in real life where - in equilibrium - the two subsystems have the same temperature. There the sum in eq. 2 can be very well approximated by the largest term

$$\underset{N_1}{}g_2(N_1,E_1)g_1(NN_1,EE_1)g_2(N_1,E_1^{max})g_1(NN_1,EE_1^{max})$$
(7)
Only if this approximation is justified it makes sense to describe the system as a canonical system, i.e. in the variables T and $`\mu `$. Then the system obeys all the laws of equilibrium thermodynamics, that the entropy becomes additive and that subsystems in equilibrium have the same temperature and the same chemical potential.

Is this approximation justified in heavy ion collisions where the systems are still rather small containing at most around 400 particles? We use the above mentioned microcanonical SMM calculations to find an answer. We identify the two subsystems with the fragments with $`Z3`$ and the LCP’s, respectively, and calculate the temperature difference between the subsystems.

In fig. 7, top, we display the number of events (left) and the average number of IMF’s (right) as a function of the total number of nucleons entrained in the fragments and the total fragment energy. We see a rather broad distribution. In the middle we present the difference of the temperature and the chemical potential in the two subsystems (LPC’s and IMF’s)as a function of the total mass and the total energy of all IMF’s. For calculating the temperature difference we have assumed that $`\frac{1}{T_1}\frac{1}{T_2}\frac{\mathrm{\Delta }T}{T^2+\mathrm{\Delta }T^2/4}`$ where T is assumed to be 7 MeV. We see a rather broad distribution of the chemical potentials and of the temperatures. The probability distribution of fragments to be emitted in a microstate which shows a $`\mathrm{\Delta }T`$ and a chemical potential difference of $`\mathrm{\Delta }\mu `$, respectively, between the two subsystems is plotted in the bottom row. The chemical potential $`\mathrm{\Delta }\mu `$ in units of the temperature fluctuates by about to 40% around the mean value 0 whereas the variation of the distribution of the temperature difference between the subsystems is about 1.7 MeV. More precisely the temperature in the subsystem of the fragments cannot be determined. Consequently, the temperature values in the caloric curve have to get an error bar of this size which renders any conclusion about a phase transition useless because one cannot distinguish whether T increases slowly as a function of the excitation energy (cross over) or whether it stays constant (first order phase transition).

This variance of the temperature in the subsystems renders as well the application of canonical thermodynamics to heavy ion reaction useless. There is nothing like a well defined temperature of the matter created in these reactions and therefore the hope to measure the system temperature is in vain.


format: text

text:
# Observation of power-law scaling for phase transitions in linear trapped ion crystals

## Abstract

We report an experimental confirmation of the power-law relationship between the critical anisotropy parameter and ion number for the linear-to-zigzag phase transition in an ionic crystal. Our experiment uses laser cooled calcium ions confined in a linear radio-frequency trap. Measurements for up to 10 ions are in good agreement with theoretical and numeric predictions. Implications on an upper limit to the size of data registers in ion trap quantum computers are discussed.

Ions confined in linear radio-frequency traps, and cooled by laser radiation, will condense into a crystalline state. Such crystals are the most rarefied form of condensed matter known . Besides being of inherent scientific interest for this reason, cold trapped ions have a growing number of applications, notably spectroscopy , frequency standards , and quantum computing . The existence of different kinds of phase transitions of these crystals has been known for some time  and has been the subject of various theoretical and numeric studies . Previous experimental work identified different crystal phases/configurations in a quadrupole ring trap . Here we explicitly investigate the $`\mathrm{𝑡𝑟𝑎𝑛𝑠𝑖𝑡𝑖𝑜𝑛}`$ between two of these phases: the linear and the zigzag configurations. We report the first experimental confirmation of one of the key theoretical/numeric predictions for the linear-to-zigzag transition, namely the existence of a power-law relating the critical anisotropy parameter to the number of ions in the crystal. Further, we discuss the usefulness of this power-law expression in determining the ultimate size of a quantum logic register realizable using a single ion trap.

The potential energy of a crystal of $`N`$ identical ions of mass $`M`$ and charge $`e`$ confined in an effective three-dimensional harmonic potential is

$`U(𝐫_1,𝐫_2,\mathrm{}𝐫_N)`$ $`=`$ $`{\displaystyle \frac{M(2\pi )^2}{2}}{\displaystyle \underset{n=1}{\overset{N}{}}}\left(\nu _x^2x_n^2+\nu _y^2y_n^2+\nu _z^2z_n^2\right)`$ (2)
$`+{\displaystyle \frac{e^2}{8\pi ϵ_0}}{\displaystyle \underset{\stackrel{n,m=1}{mn}}{\overset{N}{}}}{\displaystyle \frac{1}{|𝐫_n𝐫_m|}},`$
where $`𝐫_n=(x_n,y_n,z_n)`$ is the postion vector of the $`nth`$ ion, and $`\nu _x`$, $`\nu _y`$ and $`\nu _z`$ are the trapping potential frequencies in the three directions. Assume the trapping potentials are approximately equal in the two transverse directions ($`x`$ and $`y`$), so that $`\nu _x\nu _y\nu _r`$, and the trapping potential in the axial ($`z`$) direction is different than in the other two directions. This anisotropy is characterized by the parameter $`\alpha =(\nu _z/\nu _r)^2`$. Work by Schiffer  and Dubin  predicts that the ions undergo a phase transition from a linear to a zigzag configuration at a critical anisotropy value $`\alpha _{crit}`$. The predicted power-law scaling for $`\alpha _{crit}`$ versus the number of trapped ions $`N`$ is:

$$\alpha _{crit}=cN^\beta ,$$
(3)
where the constants found empirically by Schiffer are $`c=2.53`$ and $`\beta =1.73`$ . Later in this paper we provide a simple alternative approach which produces qualitative agreement with Schiffer’s results but is based on a stability analysis of the transverse oscillatory modes, rather than on numerical simulations of the equilibrium crystalline configurations.

We studied the zigzag phase transition using strings of 40Ca+ ions confined in a linear radio-frequency quadrupole trap. Trap electrodes (see Fig. 1 or Ref. ) consist of four rods and two conical endcaps. An rf drive at $`\mathrm{\Omega }/2\pi =6`$ MHz is applied to two diagonally opposite rods to produce a radial confining pseudo-potential  with frequency $`\nu _r`$ proportional to the rf trapping voltage $`V_0`$ (for example, ). Axial confinement is provided by a static potential applied to the two endcap electrodes, resulting in a harmonic potential with frequency $`\nu _z`$ proportional to the square root of the endcap voltage. Anharmonic contributions to the axial potential energy are estimated to be less than $`3\times 10^7`$ the harmonic contribution, and are thus ignored.

Calcium ions are introduced into the trap by intersecting an atomic beam with an electron beam. Laser cooling on the S1/2 to P1/2 transition at 397 nm (see Fig. 1) cools the ions so that a crystallized string is formed along the trap axis. The resulting 397 nm ion fluorescence is detected in the horizontal plane along a direction perpendicular to the trap axis by both a CCD camera and a photomultiplier tube. An 866 nm laser returns ions falling into the D3/2 metastable state to the cooling cycle. The 397 nm light is produced by doubling 794 nm diode-laser light which is intensified by a tapered amplifier; the 866 nm beam is also produced by a diode laser.

For each measurement of $`\alpha _{crit}`$, the rf trapping potential is lowered, keeping the endcap voltage constant, until the transition to a zigzag pattern is observed. The rf potential is then repeatedly raised and lowered to determine the reproduceability of the transition point and to rule out hysteresis. At the identified critical rf voltage, we measure the radial and axial resonant frequencies. The measurements are repeated for a range of endcap voltages and then for different numbers of ions $`N`$.

To determine the rf trapping voltage at which the linear crystal becomes zigzag, the voltage is lowered until the equilibrium position of at least one of the ions moves visibly (0.5-1.5 pixels $`0.31.0`$ $`\mu `$m) on the CCD camera image. Axial shifts are detected more easily than radial shifts; therefore, movement in the axial direction is used as the diagnostic to identify the phase transition. Moreover, observations show the spacing between the two innermost ions to be most sensitive to this axial reorganization. The critical rf voltage is determined to within $`\pm `$0.1 V out of 3-8 V on the synthesizer generating the rf trapping voltage. This 0.1 V resolution, corresponding to 8-16 kHz uncertainty in radial frequency, is limited by the voltage change necessary to move an ion visibly on the CCD camera rather than by the synthesizer itself.

The radial and axial frequencies are individually measured by applying an external drive and observing melting of the crystal into a diffuse but still stable cloud when the drive frequency is resonant. Note that although the ions are driven at a center-of-mass frequency, coupling to higher order modes heats the string until it melts. The external drive is the output of a function generator, attenuated and capacitively coupled through 1000 pF to one endcap electrode. Resonant frequencies are measured reproduceably to within $`\pm `$1-2 kHz by stepping repeatedly back and forth through the typical 1-4 kHz range in which melting is observed. Since the azimuthal symmetry of the trap is not perfect, we measure two separate radial frequencies, differing from each other by 1-2$`\%`$. The data presented includes only the smaller radial frequency (corresponding to the weaker axis of the radial well) which is responsible for the zigzag instability onset.

Figure 2 shows the radial versus axial frequencies at which the linear-to-zigzag phase transition was observed for string lengths up to $`N=10`$. Figure 3 shows the average $`\alpha _{crit}`$ as a function of $`N`$, and demonstrates the predicted scaling behavior for these critical parameters. In both figures, the solid lines are results from our theoretical analysis (presented below) with no free parameters. Although the measurements of $`\alpha _{crit}`$ lie slightly above theory, the overall agreement within the 5-6$`\%`$ error bars is quite remarkable for a second order phase transition.

In Fig. 2, vertical error bars are dominated by the 8-16 kHz resolution in determining the onset of the zigzag mode as described above, but also include (in quadrature) the 1-2 kHz measurement resolution of the radial frequency. Horizontal errors representing the axial frequency measurement resolution are smaller than the point size and not shown, but are included when calculating the error for a measured $`\alpha _{crit}`$. In Fig. 3, each plotted value is a weighted average of measured $`\alpha _{crit}`$’s. Corresponding error bars are taken to be weighted average errors rather than quadrature sums of the individual $`\alpha _{crit}`$ errors, because they represent detection resolution rather than statistical errors and performing many measurements would not be expected to reduce their size.

An unaccounted source of error arises from the threshold axial shift of 0.5-1.5 pixels required to detect a phase transition. Any smaller axial shifts due to onset of the zigzag mode at higher rf voltages go undetected. Investigations tracking the axial position of an ion near the center of a linear crystal as a function of rf trapping voltage were performed to estimate the size of this detection threshold error. Examples such as those shown in Fig. 4 reveal that once a phase transition is detected, the axial equilibrium moves roughly linearly with rf voltage at a rate of 0.3-1.1 pixels per 0.1 V step of the drive synthesizer. Not surprisingly, the weaker trap settings – low endcap voltage and small ion number – produce the larger amount of movement per voltage step but also require the larger threshold shift. Underestimates in the critical trapping voltage due to this effect, therefore, end up of the same order as the 0.1 V resolution error already assigned. This detection threshold error accounts in part for the slightly high trend exhibited by the measured $`\alpha _{crit}`$’s. Alternatively, ignoring the detection threshold error, our measurement can instead be interpreted as an upper bound on $`\alpha _{crit}`$.

Trapping and cooling parameters were varied outside their normal ranges to rule out systematic effects. Detuning the 397 nm laser as far as possible to longer or shorter wavelengths, while still maintaining crystallization, altered neither the measured resonant frequencies nor the measured critical rf voltage. Similar tests involving the 866 nm laser wavelength showed no effects. Finally, for each new endcap voltage, care was taken to move the ion string to the radial center of the rf quadrupole potential. Ions not centered in this potential experience micromotion at the 6 MHz trap frequency. This driven motion, detected as a modulation on the ion fluorescence intensity, was used as a diagnostic to find the optimal position along one radial direction . The diagnostic of micromotion for the orthogonal radial direction was the observation of a radial shift in ion position upon weakening or strengthening the pseudo-potential well, because an ion’s position in the pseudo-potential will depend on the potential well strength unless the ion is centered . Ions were moved radially by applying voltages to the compensation electrodes, which were typically tuned to $`\pm `$1 V; tests detuning them by 5-10 V revealed no change in measured resonant frequencies or in critical rf voltages.

We now present a simple theoretical analysis of the onset of zigzag instability. The equilibrium positions of the ions are determined by the condition that the potential be a minimum, i.e. $`_nU=0(n=1,\mathrm{}N)`$. For the case of strong anisotropy ($`\alpha 1`$), the ions are configured along a line in the $`z`$-direction. The solutions to the equilibrium equations in this case have been investigated by various authors (see, for example ); we denote the equilibrium position of the $`nth`$ ion by $`\overline{𝐫}_n=(0,0,\overline{z}_n)`$.

Small oscillations of the ions about their equilibrium positions are described by the Lagrangian

$`L`$ $`=`$ $`{\displaystyle \frac{M}{2}}\left[{\displaystyle \underset{n=1}{\overset{N}{}}}\dot{\zeta }_n^2\omega _z^2{\displaystyle \underset{n,m=1}{\overset{N}{}}}A_{n,m}\zeta _n\zeta _m\right]`$ (4)
$`+`$ $`{\displaystyle \frac{M}{2}}\left[{\displaystyle \underset{n=1}{\overset{N}{}}}\dot{\xi }_n^2\omega _x^2{\displaystyle \underset{n,m=1}{\overset{N}{}}}B_{n,m}\xi _n\xi _m\right]`$ (5)
$`+`$ $`{\displaystyle \frac{M}{2}}\left[{\displaystyle \underset{n=1}{\overset{N}{}}}\dot{\eta }_n^2\omega _y^2{\displaystyle \underset{n,m=1}{\overset{N}{}}}B_{n,m}\eta _n\eta _m\right],`$ (6)
where $`(\xi _n,\eta _n,\zeta _n)`$ are the displacements of the $`nth`$ ion from its equilibrium position in the $`(x,y,z)`$ directions, respectively. The coupling matrices $`A_{n,m}`$ and $`B_{n,m}`$ are given by the formulae (following Ref. )

$`A_{n,m}`$ $`=`$ $`\{\begin{array}{cc}1+2{\displaystyle \underset{\stackrel{p=1}{pm}}{\overset{N}{}}}{\displaystyle \frac{\mathrm{}^3}{\left|\overline{z}_m\overline{z}_p\right|^3}}\hfill & \mathrm{if}n=m,\hfill \\ \text{}{\displaystyle \frac{2\mathrm{}^3}{\left|\overline{z}_m\overline{z}_n\right|^3}}\hfill & \mathrm{if}nm.\hfill \end{array}`$ (9)
$`B_{n,m}`$ $`=`$ $`\left({\displaystyle \frac{1}{\alpha }}+{\displaystyle \frac{1}{2}}\right)\delta _{n,m}{\displaystyle \frac{1}{2}}A_{n,m},`$ (10)
where $`\mathrm{}=\left[e^2/4\pi ϵ_0M(2\pi \nu _z)^2\right]^{1/3}`$ is a length scale and $`\delta _{n,m}`$ is the Kronecker delta.

The eigenvectors and eigenvalues of the real, symmetric, positive-definite matrix $`A_{n,m}`$ define the normal modes of oscillation of the ions along the $`z`$-direction. Because of their importance to quantum computing, they have been studied in some detail . The eigenvectors are defined by the formula $`_{m=1}^NA_{n,m}b_m^{(p)}=\mu _pb_n^{(p)}`$, where $`\mu _p>0`$ is the eigenvalue, $`b_n^{(p)}`$ the normalized eigenvector, and $`p`$ $`(=1,\mathrm{}N)`$ the mode index (the modes being enumerated in order of increasing eigenvalue).

From the definition of the coupling matrix for the radial oscillations, $`B_{n,m}`$, we see that it has identical eigenvectors to $`A_{n,m}`$, but has different eigenvalues:

$$\underset{m=1}{\overset{N}{}}B_{n,m}b_m^{(p)}=\left(\frac{1}{\alpha }+\frac{1}{2}\frac{\mu _p}{2}\right)b_n^{(p)}.$$
(11)
For these oscillations, the eigenvalues are no longer always positive. Indeed, as $`\alpha `$ is increased, a critical value occurs for which one of the eigenvalues of the radial oscillations becomes zero. Beyond this point the radial oscillation is unstable and so is the linear configuration: this marks the onset of the zigzag mode. The critical value of $`\alpha `$ for which this occurs is given by the exact result $`\alpha _{crit}(N)=2/(\mu _N1)`$, where $`\mu _N`$ is the largest eigenvalue of the matrix $`A_{n,m}`$ for a given number of ions $`N`$. The value of $`\mu _N`$ in general must be determined by numerically diagonalizing the matrix $`A_{n,m}`$ for each value of $`N`$, but the expression for $`\alpha _{crit}(N)`$ is exact.

We can approximate $`\alpha _{crit}(N)`$ as the power-law in Eq. 3 in order to compare it to our experimental measurements. Table I compares the constants $`c`$ and $`\beta `$ (see Eq. 3) derived from the measurements, from a fit to our theory (over $`N=310`$), and from Schiffer’s fit. We emphasize that no matter how well any model calculates $`\alpha _{crit}`$ for a specific $`N`$, the power-law deduced for $`\alpha _{crit}`$ versus $`N`$ is still approximate and must be obtained from a fit of many calculated $`\alpha _{crit}`$’s. In light of this, the three sets of coefficients are close, particularly for $`\beta `$. Small discrepancies may be explained by the fact that the exact expression for $`\alpha _{crit}(N)`$ is $`\mathrm{𝑛𝑜𝑡}`$ a power-law, and that Schiffer’s predicted values were deduced from a fit of 10 $`\alpha _{crit}`$ values over the range $`N=2500`$ whereas we fit over just the experimental range $`N=310`$. Fitting, instead, over all values in the range $`N=2100`$, yields the constants $`c=2.88\pm 0.03`$ and $`\beta =1.773\pm 0.003`$, which are more appropriate (as are Schiffer’s) for large $`N`$ applications of the power-law discussed below.

Using the power-law expression, we illustrate how to estimate the maximum number of ions it is reasonable to confine along the trap axis. This gives one possible upper limit to the size of data registers for quantum computation with on-axis ions in a linear trap. Rewriting Eq. 3 gives an upper limit $`N_{crit}`$ to the number of ions one can trap in the linear configuration: $`N_{crit}=(\nu _z/\sqrt{c}\nu _r)^{2/\beta }`$. Experimental limitations will dictate how large the ratio $`\nu _z/\nu _r`$ can be made, but this expression should prove useful in optimizing future designs of ion traps.

In summary, a simple theoretical analysis has produced the critical anisotropy parameter $`\alpha _{crit}`$ for the linear-to-zigzag phase transition. Experimental measurements have provided a confirmation of the predicted power-law scaling for $`\alpha _{crit}`$ versus $`N`$ which had until now been untested. From this scaling we obtain an expression for an upper limit to the size of linear data registers for quantum computation in these traps.

We are grateful to D. J. Berkeland and A. G. Petschek for helpful conversations. This work was supported by the National Security Agency.


text:
# 1 Introduction

## 1  Introduction

Symmetry principles have generally origin from experimental evidence and therefore their validity has to be verified to an ever increasing degree of precision, or falsified. Lorentz Invariance (LI) is no exception, being based on experimental facts (e.g. the impossibility of verifying the motion of the Earth from laboratory experiments) and having innumerable experimentally testable consequences (e.g. constancy of speed of light, equivalence of physics in different reference frames ….).

In a recent paper (Aloisio et al. 2000) we discussed the possibility of using Cosmic Ray (CR) experiments to put very stringent constraints on the validity of LI. That this is the case can be intuitively motivated in the following way. Consider the process giving rise to the Greisen, Kuzmin, Zatsepin cut-off (Greisen 1966, Zatsepin and Kuzmin 1966), i.e. pion photoproduction. In the terrestrial laboratory the reaction is $`\gamma p\pi N`$ and has a threshold at $`E_{th}100`$ MeV for a (proton) target at rest. In the same frame the reaction initiated by UHE CRs on Cosmic Microwave Background (CMB) photons as a target ($`p\gamma _{CMB}\pi N`$) has a threshold of $`510^{19}`$ eV. In this frame the two reactions appear very different: in particular a $`10^{20}`$ eV proton needs only an extremely tiny fraction of its energy ($`10^{23}`$) to make a transition in the final state containing a pion. Obviously there is nothing misterious in this, LI just implying that the two reactions are exactly the same taking place in two reference frames in motion with a relative Lorentz factor $`\gamma 10^{11}`$. However it is clear that even very small deviations from strict relativistic invariance are likely to profoundly modify the value of the thresholds and of the associated absorption cut-offs, in a way in principle experimentally verifiable. A similar process, giving rise to an absorption threshold for VHE-UHE Cosmic $`\gamma `$ rays, is $`\gamma \gamma _{BCKG}e^+e^{}`$ where again $`\gamma _{BCKG}`$ is a low energy background photon: IR (Stecker 1999), microwave (Nikishov 1962, Goldreich and Morrison 1964, Gould and Schreder 1966) or radio (Protheroe and Biermann, 1996, 1997).

On the other hand, it has been been conjectured since the 50’s (Wheeler 1957) that strict Lorentz invariance is likely to be profoundly modidified at the scales at which quantum grvitational effects begin to be relevant: at Planck distances (times) the topology of space-time may become highly non trivial, making the definition of distance essentially imposible thus profoundly modifing our picture of the physical world.

In this talk I present arguments showing, on very general grounds and without referring to specific models, that in the processes which give rise to absorption thresholds for the propagation of cosmic rays in the Universe (due the onset of particle production on universal background radiation) it is possible to test the validity of LI to a very high degree of precision; in particular, if deviations from LI are ascribed to QG effects, they can be studied down to length scales orders of magnitude smaller than Planck length, even using beam particles with energy far less than the Planck mass!

I want however to remark that we do not advocate violation of LI to explain the features of the spectrum of UHE Cosmic Rays (e.g. the apparent absence of the GZK cut-off). Our knowledge of the sources of the highest energy Cosmic and $`\gamma `$ rays is rather poor, so that we prefer to stress the capability of CR experiments to test such a fundamental symmetry. However the situation is going to change drastically with the new generation of extremely large collecting area Cosmic Ray Experiments (Cronin 1992), or is already changing, for instance due to the high statistics studies of the VHE spectrum of extragalactic $`\gamma `$-sources like Mkn421 and 501 (see for instance Aharonian et al. 1999, Guy et al. 2000).

## 2  Absorption thresholds in non-LI world

As intuitively motivated above, the reactions that can lead to a verification of LI are those in which VHE and UHE Cosmic and $`\gamma `$ rays interact with a universal diffuse (photon) background, that can be the very well known Cosmic Microwave Background Radiation (energy of maximum radiance $`ϵ10^3`$ eV), the less known Far Infrared Radiation ($`ϵ10^2`$ eV) or the hypotetical Radio Background ($`ϵ10^7`$ eV). On this diffused radiation UHE protons and VHE $`\gamma `$’s can interact and produce secondary particles, mostly $`\pi `$ in the case of primary protons, $`e^+e^{}`$ in the case of $`\gamma `$’s.

In all the cases the processes can be seen as Lorentz boosted from terrestrial laboratory with boost factors ranging from $`10^7`$ to $`10^{14}`$. The cross sections for these processes are large ($`10^{25}`$ $`cm^2`$) so that the energetic particles are rapidly degraded in energy with an absorption length of the order of tens Kpc for $`\gamma `$ absorption on CMBR to tens of Mpc for protons absorption in CMBR and $`\gamma `$ absorption on Radio background. Both these facts are important quality factors for using CR’ s to test LI. The first gives the range of parameters in which SR can be tested, while the second suggests that possible modifications of these reactions are expected to produce a (in principle easily) detectable signal.

To make the test quantitative one needs a parametrization of possible Lorentz violations 222for earlier discussions of LI breaking see: Kirzhnits and Chechin 1971, Gonzalez-Mestres 1997, Coleman and Glashow 1997, Amelino-Camelia et al. 1997. A very similar approach is in Amelino-Camelia and Piran 2000.. The constancy of speed of light ($`c(=1)`$ in the following) implies the existence of an invariant interval $`ds^2=dt^2dx^2`$ (=0 for light signals), and that the norm of any four-vector is invariant: in particular for the Energy-Momentum four-vector of a particle of mass $`m`$ one has the dispersion relation:

$$P_\mu P^\mu =m^2=E^2|\stackrel{}{p}^2|$$
(1)
This is the relation that we modify in order to parametrize violations of LI. We follow a phenomenological approach and we do not refer to any specific model (see however next section). Our guiding principles are:

1) Violations are universal, $`i.e.`$ do not depend on particle type;

2) Preserve rotational invariance;

3) Violations are an high energy phenomenon, vanishing at low momenta.
Finally we only consider the $`p<<M_P`$ range (relevant for the experiments we cosider). We therefore modify the dispersion relation as ($`p=|\stackrel{}{p}|`$):

$$E^2p^2=m^2+p^2f(\frac{p}{M})+.\mathrm{}.$$
(2)
The dots stand for terms that are subleading in the regime we are considering; also in this regime the difference between $`E`$ and $`p`$ in the RHS of eq. 2 is a higher order correction. $`M`$ is a mass parametrizing the violation of LI, which is expected to be of the order of the Planck mass $`M_P`$ if violations originate from quantum gravity effects, and $`f(0)=0`$. Under these conditions we can expand $`f(p/M)`$ and, reabsorbing numerical coefficients in M, we can write, for the leading correction:

$$E^2p^2m^2\pm p^2(\frac{p}{M})I_\pm $$
(3)
$$E^2p^2m^2\pm p^2(\frac{p^2}{M^2})II_\pm $$
(4)
being of first (second) order in the (small) parameter $`p/M`$.

It is worth noticing that we could have modified the Lorentz transformations of four-momenta

$$E^{}=\gamma MF(\frac{E+\beta p}{M});p^{}=\gamma MG(\frac{\beta E+p}{M})(F(0)=G(0)=0)$$
(5)
which, with the appropriate choice of F and G con lead to (3) or (4). In this procedure there is however a large arbitrariety (eq. (5) is not even the most general) and we do not pursue it here 333Notice that if $`F(\mathrm{}),G(\mathrm{})\mathrm{}`$ eq. (5) implies that there is a maximum value for energy (momentum) (Garay 1999)..

We also $`assume`$ Energy-Momentum conservation. This is not granted, if the violation of LI is associated to a violation of traslational invariance, and has to be checked in specific models. However this assumption is legitimate if we want to put bounds on the violations of LI from the $`observation`$ of the absorption thresholds.

We have now all the ingredients to compute the value of the threshold momenta in this new framework. It is important to notice that the computation must be performed in a specific frame, which we take the one in wich the diffuse radiation is isotropic (neglecting the motion of the Earth), taking into account energy-momentum conservation and the modifications of the dispersion relations. This leads to algebraic equations (Aloisio et al. 2000)

$$\pm \alpha _Ix^3+x1=0$$
(6)
$$\pm \alpha _{II}x^4+x1=0$$
(7)
where $`x=p_{th}/p_0`$ and $`p_0`$ are the thresholds as $`M\mathrm{}`$ 444$`p_0(m_\pi m_p)/(2\omega )`$ for pion production by protons, $`p_0m_e^2/(2\omega )`$ for pair production by $`\gamma `$. $`\omega `$ is the energy of background photons., and (in parentesys the values for pion production)

$`\alpha _I`$ $`=`$ $`{\displaystyle \frac{p_0^3}{8m_e^2M}}({\displaystyle \frac{p_0^3}{m_p^2M}})`$ (8)
$`\alpha _{II}`$ $`=`$ $`{\displaystyle \frac{3p_0^4}{16m_e^2M^2}}({\displaystyle \frac{3p_0^4}{m_p^2M^2}})`$
Notice that if we require that $`p_{th}p_0`$ then we must have $`\alpha _I`$ ($`\alpha _{II}`$) $`<<1`$, which, taking as example the pion production by protons, implies $`M>10^{14}M_P`$ ($`M>10^3M_P`$) respectively. Namely, a verification of the LI threshold momentum value would put a limit on the violation parameter (much) larger than the Planck mass.

The qualitative behaviour of the solution is presented in Fig. 1. In case of positive modification the threshold moves towards lower values as $`M`$ moves away from $`\mathrm{}`$, while it increases, and becomes rapidly unphysical ($`i.e.`$ negative or complex) for negative modifications; in this case the reaction becomes kinematically forbidden.

In a more quantitative way, in the following table we present the values of $`p_{th}/p_0`$ for pair production, assuming $`M=M_P`$:

|  | Infrared | Microwave | Radio |
| --- | --- | --- | --- |
| $`I_+`$ | $`1`$ | $`0.06`$ | $`510^7`$ |
| $`I_{}`$ | no solution | no solution | no solution |
| $`II_+`$ | $`1`$ | $`1`$ | $`210^3`$ |
| $`II_{}`$ | $`1`$ | $`1`$ | no solution |

Notice however that, for positive modifications the process $`\gamma e^+e^{}`$ becomes allowed, giving rise to absorption even in absence of target (Coleman and Glashow 1997).

If we leave M as free parameter, assuming experimental verification of thresholds (with a $`100\%`$ uncertainty in momentum determination) we have

|  | Infrared | Microwave | Radio |
| --- | --- | --- | --- |
| $`I_+`$ | $`\frac{M}{M_P}0.2`$ | $`\frac{M}{M_P}800`$ | $`\frac{M}{M_P}10^{18}`$ |
| $`I_{}`$ | $`\frac{M}{M_P}6`$ | $`\frac{M}{M_P}310^4`$ | $`\frac{M}{M_P}810^{19}`$ |
| $`II_+`$ | $`(\frac{M}{M_P}310^8)`$ | $`(\frac{M}{M_P}710^6)`$ | $`\frac{M}{M_P}10^5`$ |
| $`II_{}`$ | $`(\frac{M}{M_P}310^7)`$ | $`(\frac{M}{M_P}10^4M_P)`$ | $`M10^6`$ |

Analogously for pion production (the process giving rise to the GZK cut-off), for $`M=M_P`$

|  | GZK |
| --- | --- |
| $`I_+`$ | $`210^5`$ |
| $`I_{}`$ | no solution |
| $`II_+`$ | $`0.02`$ |
| $`II_{}`$ | no solution |

Also in this case for positive modifications the process $`p\pi N`$ becomes kinematically allowed. If we leave M as free parameter, assuming experimental verification of thresholds (again within a factor 2 in energy)

|  | GZK |
| --- | --- |
| $`I_+`$ | $`M310^{13}M_P`$ |
| $`I_{}`$ | $`M10^{15}M_P`$ |
| $`II_+`$ | $`M500M_P`$ |
| $`II_{}`$ | $`M610^3M_P`$ |

It is important to notice that these solutions confirm the intuitive expectations described in the introduction: in fact, the values of threshold momenta are (in most cases) profoundly modified by the modifications of dispersion relations produced by violations of Lorentz Invariance.

## 3  Theoretical Motivation

The above analisys is quite general and does not refer to specific models. It is however important to notice that there are, in the literature, models which lead violations of LI in the form discussed above (see for instance Amelino-Camelia et al. 1997, Amelino-Camelia et al. 1998 and for a more complete list Aloisio et al. 2000).

More generally, if violations of LI are ascribed to Quantum Gravitational effects we can classify the violations in a general way. In fact QG effects imply that the metric of space-time is non trivial when examined near the Planck scale

$$g_{\mu \nu }=\eta _{\mu \nu }+h_{\mu \nu }$$
(9)
where $`\eta _{\mu \nu }=`$ diag(1,-1,-1,-1) is the flat metrics and $`h_{\mu \nu }`$ is a term fluctuating in the vicinity of the Planck length 555Notice that we cosider (8) as a phenomenological description. In QG the concept of a background (flat) metric might be ill-posed..

Consider $`g_{\mu \nu }P^\mu P^\nu `$ in such a metric so that

$$g_{\mu \nu }P^\mu P^\nu =\eta _{\mu \nu }P^\mu P^\nu +h_{\mu \nu }P^\mu P^\nu $$
(10)
A particle traveling with energy $`E`$ averages the fluctuations over a scale $`\lambda /l_P`$ where $`\lambda `$ is its de-Broglie wavelength ($`\lambda /l_P10^8`$ for a $`10^{20}`$ eV proton) so that:

$$g_{\mu \nu }P^\mu P^\nu _{\frac{\lambda }{l_P}}=\eta _{\mu \nu }P^\mu P^\nu +h_{\mu \nu }P^\mu P^\nu _{\frac{\lambda }{l_P}}=m^2E^2p^2+(\frac{l_p}{\lambda })^n\overline{h}_{\mu \nu }P^\mu P^\nu +\mathrm{}$$
(11)
where the two parameters $`n`$, $`\overline{h}_{\mu \nu }`$ describe the (possible) violations of LI due to QG effects; clearly in the spirit of the discussion of Sect. 2, if $`n>2`$ then the effects are in any case likely to be negligible when $`p<<M_P`$. On the other hand $`\overline{h}`$ may be:

1. $`\overline{h}_{\mu \nu }=0`$: the dispersion relation is not modified;
2. $`\overline{h}_{\mu \nu }\eta _{\mu \nu }E^2p^2f(m^2),f(0)=0`$: this gives the mildest violation, that does not affect photon propagation.
3. $`\overline{h}_{\mu \nu }`$ non diagonal: $`E^2p^2m^2`$. For instance: $`\overline{h}_{ij}=\delta _{ij},\overline{h}_{00}=\overline{h}_{0i}=0`$
is a possible (among many others) choice that leads to $`E^2p^2=m^2\pm p^2(p/M)^n+..`$ and for $`n=1(2)`$ generates the violations of type $`I_\pm (II_\pm )`$.

Even in case 1, in general, one can have (Ford and Yu 1999)

$$\overline{h}_{\mu \nu }\overline{h}_{\mu ^{}\nu ^{}}0(E^2p^2m^2)^2_{\frac{\lambda }{l_P}}0$$

with possibly observable effects.

The vacuum of Quantum Gravity is often described as filled by virtual, Planck mass Black Holes (Hawking 1995). Even if the full vacuum metric might be difficult to manage, it is interesting to notice that in the field of a single BH the dipersion relation becomes: $`E^2p^2=m^2\pm \frac{l_P}{L}(E^2+p^2)`$ which holds when $`(\frac{l_P}{L}<<1)`$, and we expect $`L\lambda `$. The average over the ensemble of fluctuating BHs is non trivial, since it depends on the QG dynamics; it seems however natural to expect that $`(E^2p^2m^2)^2^{\frac{1}{2}}𝒪(\frac{p^3}{M_P})`$

## 4  Conclusions

Cosmic Ray experiments already in operation have the capability to investigate the validity of Lorentz Invariance to a high degree of sensitivity. There are already a number of events above the GZK cut-off and the situation will improve dramatically in a few years with the beginning of data taking of the Auger experiment. And $`\gamma `$-telescopes are already in the position of studying the spectrum of a few (up to now) extargalactic sources up to tens of TeV.

These experiments can test LI down to length scales in principle much smaller than the Planck length. This is extremely important, since it means that we do not need Planck energies to study the Planck physics if we appropriately chose the processes to study.

The experimental situation is still rather unclear: no sign of the GZK cut-off has been seen up in the proton spectrum up to a few in $`10^{20}`$ eV (Takeda et al 1999, Abu-Zayyad et al. 1999) while the $`\gamma `$ spectrum of Mkn501 although showing a bend, might be inconsistent with expectations based on a new estimate of IR background (Protheroe and Meyer 2000).

It is important however to use much caution in interpreting the experimental data: our knowledge of the possible sources of highest energy CRs is rather poor, and in general it is certainly premature to invoke violations of LI to explain experimental data. This might be true also in the case of VHE $`\gamma `$ astronomy, although, at least in the case of Mkn501, also due to (almost) simultaneous multi-wavelength observation (Guy et al. 2000), the experimental data are constraining the source spectrum. Again, the statistics is at present low, and there is uncertainty on the IR flux, but the situation is likely to improve. And, $`\gamma `$-experiments performed in the PeV range (Catanese etal., Ghia et al.), where the knowledge of the background is extremely better, are possibly the best arena for testing Lorentz Invariance.

Due to the large cross-sections involved, this kind of experiments might not need to be only of cosmic nature. In fact a terrestrial photon target, containing $`10^{21}`$ infrared ($`ϵ0.01`$ eV) photons/cm3 and 1 cm long, would have an efficiency of $`10^4`$ to convert TeV photons into $`e^+`$ $`e^{}`$ pairs. This target does not seem unfeasible: a 1 W monochromatic source would produce as many IR photons in one second, and if TeV photons could be produced (at LHC?) with sufficient intensity, this device could test models of LI violation of type I up to the Planck scale.

## 5 Acknowledgements

While the content of section 3 is responsibility of the present authors, this talk uses much of the material of Aloisio et al. (2000). We thank Piera L. Ghia and Pasquale Blasi for continuos collaboration. We thank also Angelo Galante for collaboration and Giuseppe Di Carlo for discussions.

## 6  References

Abu-Zayyad T. et al., 1999 Proc. 26th ICRC (Salt Lake City, USA), 3, 264.

Aharonian F. A. et al., 1999, Astron. & Astroph. 349, 11.

Aloisio, R., Blasi, P., Ghia, P.L. and Grillo, A.F., 2000, Phys. Rev. D62, 53010.

Amelino Camelia G., Ellis J., Mavromatos N.E. and Nanopoulos D.V., 1997, Int. J. Mod. Phys. A12, 2029.

Amelino Camelia G., Ellis J., Mavromatos N.E., Nanopoulos D.V. and Sarkar S., 1998, Nature 393, 763.

Amelino Camelia G. and Piran T., 2000, preprint hep-ph/0006210.

Catanese M. et al., 1996, Astrophys. J. 469, 572.

Coleman S. and Glashow S.L., 1997, Phys. Lett. B405, 249.

Cronin J. W., 1992, Nucl. Phys. B (Proc. Suppl.) 28B, 213.

Hawking S.W., 1995 Phys. Rev. D53, 3099.

Ford L.H. and Yu H., 1999, preprint gr-qc/9907037.

Garay J.L., 1999, Int. J. Mod. Phys. A14, 4079.

Ghia P.L. et al., 1999, Nucl. Phys. B (Proc. Suppl.) 70, 506.

Gonzalez-Mestres L., 1997, preprint physics/9702026.

Goldreich P. and Morrison P., 1964, Sov. Phys. - JETP 18, 239.

Gould R.J. and Schreder G.P., 1966, Phys. Rev. Lett. 16, 252.

Guy J., Renault C., Aharonian F.A., Rivoal M. and Tavernet J.P., 2000, preprint astro-ph 0004355

Greisen K., 1966, Phys. Rev. Lett. 16, 748.

Kirzhnits D.A. and Chechin V.A., 1971, Sov. Jour. Nucl. Phys. 15, 585.

Nikishov A.I., 1962, Sov. Phys. - JETP 14, 393.

Protheroe R.J. and Biermann P.L., 1996, Astropart. Phys. 6, 45 ; 1997, Erratum-ibid 7, 181.

Protheroe R.J. and Meyer H., 2000, preprint astro-ph/0005349.

Stecker F.W., 1999, preprint astro-ph/9904416.

Takeda M. et al., 1999, Astrophys. J. 522, 225.

Wheeler J.A., 1957, Ann. Phys. (N.Y.) 2, 604.

Zatsepin G.T. and Kuzmin V.A., 1966, Pis’ma Zh. Ekps. Teor. Fiz. 4, 114 \[JETP Lett. 4, (1966) 78\].


text:
# Beyond the Wigner Distribution: Schrödinger Equations and Terrace Width Distributions

## I Introduction

Vicinal surfaces consist of terraces divided by steps which interact with each other via a variety of mechanisms, including elastic, dipolar, and indirect electronic interactions. Since these interactions directly determine the distribution of terrace widths, the use of terrace width distributions (TWDs) to determine the step-step interactions is an important goal. Theoretical attempts towards this goal have produced differing results depending upon specific approximations Gruber67 ; Bartelt90 ; PM98 ; IMP98 ; Masson94 ; Barbier96 ; LeGoff99 .

In many cases, the potential energy $`V`$ due to the interaction between neighboring steps can be written in terms of the distance $`L`$ between the steps as

$$V(L)=\frac{A}{L^2}.$$
(1)
This model of interacting steps can be mapped directly onto the Calogero-Sutherland model Calogero69 ; Sutherland71  of interacting spinless fermions (or hard-core bosons) in one spatial dimension (Fig. 1); in this picture, the direction along the steps (the $`y`$-direction in “Maryland notation”) is interpreted as time, and the steps themselves are interpreted as the world-lines of the spinless fermions. The static properties of this system (such as the TWD) depend on $`A`$ only through the dimensionless interaction strength

$$\stackrel{~}{A}\frac{\stackrel{~}{\beta }A}{(k_\mathrm{B}T)^2},$$
(2)
where $`\stackrel{~}{\beta }`$ is the step stiffness, $`k_\mathrm{B}`$ is Boltzmann’s constant, and $`T`$ is the absolute temperature. The TWD can be calculated exactly Calogero69 ; Sutherland71  when $`\stackrel{~}{A}=1/4`$, $`0`$, or $`2`$.

An argument QuantumTransport  from random matrix theory originally due to Wigner suggests that these three exact TWDs can be approximated by

$$P(s)=a_\varrho s^\varrho \mathrm{exp}(b_\varrho s^2),$$
(3)
where $`sL/L`$, the constraint $`s=1`$ gives $`b_\varrho `$ \[see Eq. (32)\], $`a_\varrho `$ is a normalization constant, and the relationship between $`\stackrel{~}{A}`$ and $`\varrho `$ is Calogero69 ; Sutherland71

$$\stackrel{~}{A}=\frac{\varrho }{2}\left(\frac{\varrho }{2}1\right).$$
(4)
These “Wigner distributions” have proven quite successful and have been widely used MehtaRanMat ; Haake ; Guhr98  in nuclear physics (in which case $`s`$ is the spacing between energy levels), random matrix theory, quantum transport, and universal conductance fluctuations. There have been a number of attempts Brody73 ; Berry84 ; Hasegawa88 ; Izrailev89 ; Casati91 ; AbulMagd96 ; Kota99  to interpolate between these three special cases.

The so-called generalized Wigner distribution EP99  is defined by the simple assumption that Eq. (3) and Eq. (4) remain valid for all values of $`\stackrel{~}{A}`$. (This is only reasonable for positive or weakly negative values of $`\stackrel{~}{A}`$ L96 .) The generalized Wigner distribution appears to be in better agreement with computer simulations of vicinal surfaces Einstein01  than its competitors, and is in reasonable agreement with many experimental TWDs Giesen00 . Furthermore, it has been proved that for rational values of $`\varrho `$, $`P(s)s^\varrho `$ as $`s0`$ Forrester92 ; Haldane94 ; ZNCHa95 . However, for arbitrary values of $`\stackrel{~}{A}`$ there is no clear relation between random matrix theory and the Calogero-Sutherland model. This paper provides a formal and physically intuitive justification for Eq. (3) for arbitrary values of $`\stackrel{~}{A}`$; the resulting formalism allows us to consider more general step-step interactions than Eq. (1). Some of these results have been mentioned briefly elsewhere Einstein01  but are explained here in detail.

The organization of this paper is as follows. In Sec. II, we show that the generalized Wigner distribution can be derived from a phenomenological model, which reduces, under certain circumstances, to the two-particle Calogero model Calogero69 . This leads, by separation of variables, to a one-variable Schrödinger equation that determines the TWD. We extend this treatment to non-trivial examples in Sec. III. The significant process of extracting step-step interactions from experimental TWDs is discussed in Sec. IV. Finally, we conclude in Sec. V with some possible extensions of this work.

## II Physical Explanation of the Generalized Wigner Distribution

### II.1 Deriving a Schrödinger equation from the generalized Wigner Distribution

Given a system of steps with interactions of the form specified by Eq. (1) and a TWD given by the corresponding generalized Wigner distribution \[Eq. (3)\], we can define a wave function such that $`\psi _0^2(s)P(s)`$:

$$\psi _0(s)=a_\varrho ^{1/2}s^{\varrho /2}\mathrm{exp}(b_\varrho s^2/2).$$
(5)
Differentiating twice, we find

$`{\displaystyle \frac{\mathrm{d}^2}{\mathrm{d}s^2}}\psi _0(s)`$ $`=`$ $`\left[{\displaystyle \frac{\varrho }{2}}\left({\displaystyle \frac{\varrho }{2}}1\right)s^2(\varrho +1)b_\varrho +b_\varrho ^2s^2\right]\psi _0(s)`$ (6)
$`=`$ $`\left[\stackrel{~}{A}s^2(\varrho +1)b_\varrho +b_\varrho ^2s^2\right]\psi _0(s),`$
where we have used Eq. (4). Eq. (6) allows us, in somewhat the same spirit as the Gruber-Mullins approximation Gruber67 , to propose the following dimensionless Schrödinger equation:

$$\frac{\mathrm{d}^2}{\mathrm{d}s^2}\psi _n(s)=[\stackrel{~}{U}(s)+\stackrel{~}{V}(s)\stackrel{~}{E}_n]\psi (s),$$
(7)
where

$$\stackrel{~}{V}(s)=\stackrel{~}{A}s^2$$
(8)
is the explicit step-step interaction potential in dimensionless form and

$$\stackrel{~}{U}(s)=b_\varrho ^2s^2$$
(9)
is a dimensionless potential due to other steps not explicitly considered. (This idea is explained in the next subsection.) By inspection we see that $`\psi _0(s)`$ is the ground state wave function (it has only one antinode), with an eigenvalue given by

$$\stackrel{~}{E}_0=(\varrho +1)b_\varrho .$$
(10)
In comparing different experimental TWDs with a range of values of $`L`$ but with presumably the same step-step interaction $`V(L)`$, it is often useful to rewrite Eq. (7) in its dimensional form

$$\frac{(k_\mathrm{B}T)^2}{\stackrel{~}{\beta }}\frac{\mathrm{d}^2}{\mathrm{d}L^2}\mathrm{\Psi }_n(L)=[U(L)+V(L)E_n]\mathrm{\Psi }_n(L),$$
(11)
where $`V(L)`$ is given by Eq. (1),

$$U(L)\frac{(k_\mathrm{B}T)^2}{\stackrel{~}{\beta }L^2}\stackrel{~}{U}\left(\frac{L}{L}\right)$$
(12)
and

$$E_n\frac{(k_\mathrm{B}T)^2}{\stackrel{~}{\beta }L^2}\stackrel{~}{E}_n.$$
(13)
Throughout this paper we will alternate freely and without comment between dimensional and dimensionless representations.

Eq. (7) can be solved for all the eigenfunctions and eigenvectors QuantumProblems , so one can conveniently consider perturbations from the purely inverse-square interaction given by Eq. (1). The eigenfunctions are

$$\psi _n(s)=c_ns^{\varrho /2}\mathrm{exp}\left(b_\varrho s^2\right){}_{1}{}^{}F_{1}^{}(n,\frac{\varrho +1}{2},b_\varrho s^2),$$
(14)
where $`n`$ is any nonnegative integer, $`c_n`$ is a normalization constant, and $`{}_{1}{}^{}F_{1}^{}(a,b,x)`$ is Kummer’s confluent hypergeometric function Kummer . The corresponding eigenvalues are

$$\stackrel{~}{E}_n=(\varrho +1+4n)b_\varrho .$$
(15)
### II.2 Deriving the generalized Wigner Distribution from a Schrödinger equation

An undesirable element of the previous derivation was that it implicitly required that one of the steps be held fixed, which introduced an artificial (and unnecessary) difference between steps. Furthermore, the quadratic nature of $`\stackrel{~}{U}(s)`$ was not well explained. In this subsection we explore these issues.

The Calogero-Sutherland Calogero69 ; Sutherland71  model with an infinite number of interacting spinless fermions would be an ideal model of the vicinal surface, but as mentioned above, it is only integrable for three values of $`\stackrel{~}{A}`$. A useful alternative is to consider only two adjacent fermions explicitly and model the effects of the other fermions phenomenologically through the pressure they exert, yielding the Hamiltonian

$$=\frac{1}{2}\left(\frac{^2}{x_1^2}+\frac{^2}{x_2^2}\right)+\stackrel{~}{V}(x_2x_1)x_1𝒫_1(x_1)+x_2𝒫_2(x_2).$$
(16)
Here $`𝒫_i`$ is the pressure exerted on fermion $`i`$, which we allow to be position-dependent because the limited correlation length in the $`x`$-direction leads to an effectively finite system size comment . Accordingly, we place the $`𝒱`$ fermions immediately to the left of fermion 1 in a “box” with a fixed left wall at $`x=𝒱`$; the right wall is at $`x_1`$, the position of fermion 1. The volume (length) of this box is $`𝒱_1=𝒱+x_1`$. The pressure $`𝒫_1`$ can be expanded as follows:

$`𝒫_1(x_1)`$ $`=`$ $`𝒫+x_1\left({\displaystyle \frac{𝒫_1}{x_1}}\right)|_{x_1=0}+𝒪(x_1^2)`$ (17)
$`=`$ $`𝒫+x_1\left({\displaystyle \frac{𝒫_1}{𝒱_1}}\right)|_{𝒱_1=𝒱}+𝒪(x_1^2)`$ (18)
$`=`$ $`𝒫x_1(𝒱\kappa )^1+𝒪(x_1^2),`$ (19)
where $`𝒫𝒫_1(0)`$ and the compressibility $`\kappa `$ is given by

$$\kappa \frac{1}{𝒱}\left(\frac{𝒫_1}{𝒱_1}\right)^1|_{𝒱_1=𝒱}.$$
(20)
Likewise, we place the $`𝒱`$ fermions immediately to the right of fermion 2 in a “box” of volume (length) $`𝒱_2=𝒱x_2`$, yielding

$$𝒫_2(x_2)=𝒫+x_2(𝒱\kappa )^1+𝒪(x_22).$$
(21)
Combining Eqs. (16), (19), and (21), we obtain

$``$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{^2}{x_1^2}}+{\displaystyle \frac{^2}{x_2^2}}\right)+\stackrel{~}{V}(x_2x_1)+(x_2x_1)𝒫`$ (22)
$`+(x_1^2+x_2^2)(𝒱\kappa )^1+𝒪(x_2^3x_1^3).`$
Finally, we can make the following change of variables:

$`x_{\mathrm{cm}}`$ $`=`$ $`{\displaystyle \frac{x_1+x_2}{2}}`$ (23)
$`s`$ $`=`$ $`x_2x_10.`$ (24)
This allows us to rewrite the Hamiltonian as follows:

$``$ $`=`$ $`\left({\displaystyle \frac{1}{4}}{\displaystyle \frac{^2}{x_{\mathrm{cm}}^2}}+{\displaystyle \frac{^2}{s^2}}\right)+\stackrel{~}{V}(s)`$ (25)
$`+[𝒫s+(𝒱\kappa )^1s^2]+4(𝒱\kappa )^1x_{\mathrm{cm}}^2+𝒪(sx).`$
If terms of quadratic and higher order in Eqs. (19) and (21) can be neglected, the Hamiltonian can be separated into a part that depends only on $`s`$ and a part that depends only on $`x_{\mathrm{cm}}`$:

$$=_{\mathrm{cm}}+_s,$$
(26)
where

$$_{\mathrm{cm}}=\frac{1}{4}\frac{^2}{x_{\mathrm{cm}}^2}+4(𝒱\kappa )^1x_{\mathrm{cm}}^2$$
(27)
and

$$_s=\frac{^2}{s^2}+\stackrel{~}{V}(s)+\stackrel{~}{U}(s)$$
(28)
with

$$\stackrel{~}{U}(s)=𝒫s+(𝒱\kappa )^1s^2.$$
(29)
None of the argument so far requires the assumption that $`\stackrel{~}{V}(s)`$ is given by Eq. (8). In the next section we will examine how well this phenomenological formalism works when applied to other forms of $`\stackrel{~}{V}(s)`$.

If $`(𝒱\kappa )^1𝒫`$, we can neglect the term $`𝒫s`$ in $`\stackrel{~}{U}(s)`$. The resulting Hamiltonian $``$ can easily be seen to be the Hamiltonian of the two-fermion Calogero model Calogero69 ,

$$=\left(\frac{^2}{x_1^2}+\frac{^2}{x_2^2}\right)+\stackrel{~}{V}(x_2x_1)+(x_1^2+x_2^2)\omega ^2,$$
(30)
where we identify from Eq. (9)

$$\omega ^2=(𝒱\kappa )^1.$$
(31)
More importantly, with the identification

$$\omega =b_\varrho \left[\frac{\mathrm{\Gamma }\left(\frac{\varrho +2}{2}\right)}{\mathrm{\Gamma }\left(\frac{\varrho +1}{2}\right)}\right]^2,$$
(32)
which again comes from the constraint $`s=1`$, we see that

$$(_s\stackrel{~}{E}_n)\psi _n(s)=0$$
(33)
is just Eq. (7) and again leads to the generalized Wigner distribution.

Finally, it should be pointed out that the one particle “Calogero model,” with a Hamiltonian given by

$$=\frac{\mathrm{d}^2}{\mathrm{d}x^2}+x^2\omega ^2,$$
(34)
is really what in essence is used in the many Gaussian approximations Gruber67 ; Bartelt90 ; PM98 ; IMP98 ; Masson94 ; Barbier96 ; LeGoff99 . These approximations, however, produce conflicting functional relationships between $`\omega `$ and $`\stackrel{~}{A}`$. This is hardly surprising, since $`\stackrel{~}{V}`$ does not appear explicitly in Eq. (34). By using a two-particle Calogero model, we are able to state unambiguously the relationship between $`\omega `$ and $`\stackrel{~}{A}`$.

## III Nontrivial Extensions

Here and in the remainder of the text, for our numerical work we use models that are discrete in both the $`y`$-direction and the $`x`$-direction Bartelt90 ; Bartelt92 ; Masson94 ; LeGoff99 , in contrast to the continuum step mode, which corresponds to the description given above. Specifically, in Monte Carlo simulations we use the terrace-step-kink (TSK) model, in which the kinks can be any integral number of lattice units long. Our transfer-matrix calculations use the restricted TSK model, in which the kinks can only be of one lattice unit. In both models there are exact expressions for the stiffness of an isolated step as a function of temperature Bartelt92 . The interaction between steps is given by a specified function $`V(L)`$, as with the continuum step model. Further details of our numerical work can be found in Refs. Einstein01, ; Hailu04, .

Eq. (7) can be extended to step-step interactions of forms other than Eq. (8). In this section we will present a few arbitrary but physically motivated step-step interactions and show that the phenomenological method of the preceding section still yield TWDs in excellent agreement with numerical simulations of the TSK model.

As a first example, we consider interactions of the form

$$\stackrel{~}{V}(s)=\stackrel{~}{A}_3s^3,$$
(35)
which would be a plausible subdominant term in a Taylor expansion of the step-step interaction in $`s^1`$. Figure 2 shows TWDs obtained by Monte Carlo simulations of steps interacting with this potential for three different values of $`\stackrel{~}{A}_3`$. Since the Wigner distribution worked so well with nearest-neighbor interactions of the form given by Eq. (8), we again take $`\stackrel{~}{U}(s)s^2`$. The solutions of the Schrödinger equation are in excellent agreement with the Monte Carlo data.

Interestingly, when all steps interact (via the potential given by Eq. (8), the variance of a TWD with a fixed value of $`\stackrel{~}{A}`$ will be smaller than the variance of a TWD with the same value of $`\stackrel{~}{A}`$ and just nearest-neighbor steps interacting. Numerical evidence supports this intuition, and in fact it appears that the Wigner distribution describes TWDs from nearest-neighbor interactions better than TWDs from all steps interacting (see Fig. 3). On the other hand, the phenomenological approach which leads to the generalized Wigner distribution depends only on the interaction between nearest neighbor steps, which is identical in both cases.

The difference appears to come from terms in Eq. (25) that are neglected in Eq. (30). Since the variance is decreased by interactions with more steps, it is clear that increasing the number of interacting steps will not increase the linear term $`𝒫`$ in Eq. (25) — that would increase the variance. Instead, it appears that increasing the number of interacting steps increases the quadratic terms in Eqs. (19) and (21). This would decrease the variance and correspond to the intuitive notion that the steps “interact more strongly” while still observing the constraint $`s=1`$. Unfortunately, it also destroys the separability of the Schrödinger equation.

A variety of circumstances, such as the presence of adsorbates Akutsu  or electronic surface states Pai94 ; EinsteinPSSS , can give rise to short-ranged attractive forces between steps, oscillating step-step interaction potentials, and other complicated interactions. Under such circumstances, one typically finds terraces of two more-or-less well-defined widths; interactions with steps beyond nearest neighbors can segregate these widths, leading to “step bunching” Shenoy . This situation is somewhat analogous to liquid-gas coexistence, and as with liquid-gas coexistence the compressibility diverges. As a result, $`\stackrel{~}{U}(s)`$ is linear in $`s`$ \[Eq. (29)\].

Vicinal surfaces with both elastic and surface-state mediated electronic interactions may be characterized by potentials of the form Pai94 ; EinsteinPSSS

$$V(L)=AL^2+BL^{3/2}\mathrm{cos}(2k_FL+\varphi ),$$
(36)
where $`A`$ is determined by the elastic interactions, $`B`$ is determined by the coupling to the surface state, $`k_F`$ is the Fermi energy, and $`\varphi `$ is a phase shift. The oscillations of this potential can give rise to a coexistence as described above, together with a diverging compressibility and a linear $`U(L)`$. As a dramatic — albeit rather unphysical — example, in Fig. 4 we compare transfer-matrix TWDs derived from the corrugated interaction potential

$$\stackrel{~}{V}(s)=4\mathrm{cos}\left(2\pi s\right)$$
(37)
\[a special case of Eq. (36)\] with corresponding solutions of Eq. (7). Since the TWDs now have multiple peaks, clearly neither the generalized Wigner distribution nor a Gaussian distribution is appropriate. Again, considering that there are no free parameters, the phenomenological TWD is in excellent agreement with the transfer-matrix TWD. The slight disagreement is probably due to finite-size effects in the transfer-matrix calculation (which involved only five steps with $`L=5`$).

## IV Solving for the step-step interaction

The chief practical interest in terrace width distributions comes from a desire to use experimental data to determine the interaction potential $`V(L)`$, or equivalently $`\stackrel{~}{V}(s)`$. In this section we discuss three progressively more sophisticated methods. The first two are appealing in their simplicity, but they suffer from serious flaws, and therefore should be avoided. By building on them, we introduce the third method, which is more computationally intensive than the preceding two, but which has the flexibility necessary to be appplied to real experimental data. For more details see Ref. CuTWDs, .

### IV.1 Naive direct numerical approach

In Sec. II, we derived a Schrödinger equation that included $`\stackrel{~}{V}(s)`$ in an obvious way, and in Sec. III we have seen that similar Schrödinger equations yield TWDs for more complicated interactions that are in excellent agreement with numerical simulations of the TSK model. It is tempting to follow an analogous path when dealing with experimental data, using the following procedure:

1. Connect the experimental data points with a smooth interpolating function $`P(s)`$, such as a cubic spline.
2. Calculate $`\psi (s)\sqrt{P(s)}`$.
3. Find $`\mathrm{d}^2\psi /\mathrm{d}s^2`$.
4. The total potential is $`\stackrel{~}{V}(s)+\stackrel{~}{U(s)}E[\psi (s)]^1[\mathrm{d}^2\psi (s)/\mathrm{d}s^2]`$.
5. At large $`s`$, $`\stackrel{~}{V}(s)`$ is negligible, so the total is given by $`\stackrel{~}{U}(s)E`$, which should be constant plus either a quadratic function of $`s`$ or a linear function of $`s`$.
6. Subtract $`\stackrel{~}{U}(s)E`$ from the potential to recover the step-step interaction $`\stackrel{~}{V}(s)`$.

This method has two serious problems. Experimental TWDs will contain measurements indicating $`P(s)=0`$, which lead to division by zero in determining the total potential. Furthermore, typical statistical fluctuations make the numerical estimates of $`\mathrm{d}^2\psi /\mathrm{d}s^2`$ from experimental data extremely unreliable.

### IV.2 Fitting the TWD to a preconceived form of $`P(s)`$

One may eliminate both of these problems by first performing a least-squares fit of the TWD to some positive definite and twice-differentiable function. For example, one might use

$$P(s)=c_0s^\varrho (1+c_1s+c_2s^2)\mathrm{exp}(bs^2),$$
(38)
where $`c_0`$ is a normalization constant and $`c_1`$, $`c_2`$, $`\varrho `$, and $`b`$ all are parameters to be fitted. The fitted $`P(s)`$ can then be analyzed as above.

The shortcomings of this second approach are more subtle. Eq. (38), for example, will always yield a potential that diverges as $`s^2`$ as $`s`$ approaches zero. In fact, we have fitted Eq. (38) to the Monte Carlo data shown in Fig. 2. Instead of correctly reproducing the potential $`\stackrel{~}{V}(s)=\stackrel{~}{A}_3s^3`$, the analysis of the fitted $`P(s)`$ produced an interaction that is well approximated by $`\stackrel{~}{V}(s)\stackrel{~}{A}_2s^2`$. The intersection of this derived potential with the known, true potential is for a value of $`L`$ near $`L`$, which is not surprising — the TWD contains more information about the interaction where $`P(s)`$ is large.

### IV.3 Recommended: Fitting the interaction $`V(L)`$ to a preconceived form

The method we do recommend is more computationally demanding, but it suffers from none of the above defects:

1. Parameterize $`V(L)`$, and make a crude initial estimate of the values of the parameters. At least some information on the form of the interaction is often available — whether it is influenced by a surface electronic state, for example. In principle, $`U(L)`$ contains a term linear in $`L`$ and a term quadratic in $`L`$, as in Eq. (29); in practice, one can often use either a purely quadratic form of $`U(L)`$ (as was done to derive the generalized Wigner distribution) or a purely linear form of $`U(L)`$ (if the compressibility diverges).
2. Solve Eq. (11) numerically for the specified parameters to find $`\mathrm{\Psi }_0(L)`$.
3. Find

$$\chi ^2=\underset{L}{}\left[\mathrm{\Psi }_0^2(L)P_{\mathrm{exp}}(L)\right]^2$$
(39)
for this set of parameters, where $`P_{\mathrm{exp}}(L)`$ is the experimental TWD.
4. Iterate procedural steps 2 and 3 in a minimization routine to find parameters for which $`\mathrm{\Psi }_0^2(L)`$ best fits the experimental TWD.

In Fig. 6 we apply this approach to the Monte Carlo TWD shown in Fig. 5. We parameterize $`\stackrel{~}{V}(s)`$ by

$$\stackrel{~}{V}(s)=\stackrel{~}{A}_1s^1+\stackrel{~}{A}_2s^2+\stackrel{~}{A}_3s^3+\stackrel{~}{A}_4s^4;$$
(40)
this form amounts to a generalization of the two-parameter fit discussed in Ref. statsteps, , where the only allowed term in $`\stackrel{~}{V}(s)`$ was $`\stackrel{~}{A}_2s^2`$. The values of the fitted parameters are given in Table 1; they clearly are quite different from those of the true potential. The resulting $`\stackrel{~}{V}(s)`$ is shown in Fig. 6, where it is compared with the actual potential used to generate the Monte Carlo data. What is important, however, is that the fitted and true potentials are close over the range in $`s`$ corresponding to significantly nonzero $`P(s)`$ — in this case, $`0.5s1.5`$. The phenomenological TWD (Fig. 5) derived from these fitted parameters is in good agreement with the Monte Carlo TWD and is very difficult to distinguish from the TWD derived from the true potential via Eq. (33). This result underscores the sensitivity of the process of extracting $`\stackrel{~}{V}(s)`$ and the necessity of very good statistics in experimental TWDs statsteps  and a good idea of the functional form of $`\stackrel{~}{V}(s)`$.

For this reason, we strongly recommend making a simultaneous fit of $`V(L)`$ to experimental TWDs from several different misorientations, and preferably from several different temperatures. We are currently applying this approach to Cu surfaces vicinal to the (001) plane CuTWDs . We have selected 10 experimental TWDs to be fitted, with mean step separations ranging from $`5.44a`$ to $`9.53a`$, where the lattice constant $`a=0.255\text{ nm}`$, and with temperatures ranging from 285 K to 360 K. Our preliminary studies Yancey  show that even with experimental data of good quality, the extraction of real potentials is not a trivial matter.

Using all 10 experimental TWDs, we were able to obtain good agreement between fitted TWDs and experimental TWDs; two examples are shown in Fig. 7. The corresponding interaction potential is shown in Fig. 8; clearly, it deviates markedly from Eq. (1). The irregularities in $`V(L)`$ make it plausible that a coexistence exists between two step widths, which would result in a diverging compressibility. This is supported by the fact that the experimental TWDs are fitted with $`U(L)L`$ than with $`U(L)L^2`$; furthermore, the linear form of $`U(L)`$ would lead to “fatter” tails, as were observed in Ref. Giesen00, .

Although these results show internal consistency, the fitted interaction potential is markedly different from potentials derived from simulations of surface relaxation using an empirical many-body potential RaouafiSS ; Barreteau ; Desjonqueres ; RaouafiPRB , which are in reasonable agreement with Eq. (1). We believe that our interaction potentials are “over-fitted” to the experimental data, with the irregularities of the $`V(L)`$ being too strongly influenced by statistical noise in our experimental data. This is not really surprising, because $`V(L)`$ depends quite sensitively on $`P(L)`$.

A similar problem Robertson93 ; Mishin99  occurs in determining embedding functions for the Embedded Atom Model from experimental data; if too many parameters are fitted, the resulting embedding functions of course can fit the supplied experimental data, but they give worse agreement with experimental data not included in the fitting process. The recommended remedy Robertson93 ; Mishin99  is to divide the available data into two groups, one of which is used to perform the fits, and the other of which is used to test the robustness of the fits.

Presumably, the same remedy will work in our case. However, a more complete treatment of the problem, including estimates of the uncertainties of fitted functions, is deferred to a later paper CuTWDs .

## V Conclusions

We have seen that the generalized Wigner distribution can be derived from a Schrödinger equation somewhat in the spirit of the Gruber-Mullins approximation, and that straightforward extensions of this method work for general step-step interactions.

In Sec. II we saw that the generalized Wigner distribution can be derived exactly from the two-fermion Calogero model Calogero69 . This, in turn, can be justified from a phenomenological model in which the force on two adjacent fermions is derived from the position-dependent pressure exerted by other fermions confined in a box, the size of which is presumably related to the correlation length in the $`x`$-direction. It is worth noting that the “entropic repulsion” is handled implicitly by the uncertainty principle in the quantum mapping, so we need explicitly consider only the energetic interactions.

Nothing in this phenomenological picture requires the step-step interaction to be given by Eq. (8), and in Sec. III we demonstrate numerally that for very general step-step interactions the phenomenological picture yields TWDs in excellent agreement with numerical simulations of the TSK model. This success is particularly impressive when a coexistence between two well-defined step widths occurs, as can happen when steps bunch Shenoy . Under this circumstance the compressibility diverges, which causes the tail of the TWD to decay exponentially with $`L`$ rather than exponentially with $`L^2`$, as it does in the case of the generalized Wigner distribution. This phenomenon could explain the curious slowly decaying tails of TWDs mentioned in Ref. Giesen00, .

Since exact solutions are available for both equations in the important case in which $`\stackrel{~}{V}(s)`$ is given by Eq. (8), it is possible to find not only the TWD but also, using the methods of Ref. Bartelt92, , an improved estimate of step wandering. This will be undertaken in a separate paper.

Finally, although we have addressed this work primarily to its surface science applications, it may be of interest to research in random matrix theory as well. The Calogero-Sutherland model with $`\stackrel{~}{A}`$ = $`1/4`$, $`0`$, or $`2`$ corresponds to random matrices with specific symmetries, and attempts to interpolate between them have simply varied the fraction of matrices belonging to each symmetry in the ensemble of random matrices. It is far from clear, however, what relation the Calogero-Sutherland model with more general values of $`\stackrel{~}{A}`$ or our treatment of it has to such mixed ensembles.

## Acknowledgment

Work was supported by the NSF-MRSEC at U. of Maryland Grant No. NSF-DMR 00-80008 and benefited from interactions with J. D. Weeks, M. Uwaha, and E. D. Williams. The authors also thank M. Haftel for useful discussions on numerical methods, and M. Giesen and R. van Gastel for experimental TWDs for copper surfaces and helpful comments.