text:
# Metastability and the Casimir Effect in Micromechanical Systems

\[

## Abstract

Electrostatic and Casimir interactions limit the range of positional stability of electrostatically-actuated or capacitively-coupled mechanical devices. We investigate this range experimentally for a generic system consisting of a doubly-clamped Au suspended beam, capacitively-coupled to an adjacent stationary electrode. The mechanical properties of the beam, both in the linear and nonlinear regimes, are monitored as the attractive forces are increased to the point of instability. There ”pull-in” occurs, resulting in permanent adhesion between the electrodes. We investigate, experimentally and theoretically, the position-dependent lifetimes of the free state (existing prior to pull-in). We find that the data cannot be accounted for by simple theory; the discrepancy may be reflective of internal structural instabilities within the metal electrodes.

\] The technology of micro electro mechanical systems (MEMS) now routinely allows fabrication of miniature movable structures on-chip. This new branch of microelectronics enables the realization of a variety of miniature, fully integrated sensors and actuators, with a rapidly growing range of applications. One common building block within MEMS is the capacitively-tunable resonator (CTR), comprising a movable electrode that serves as a mechanical resonator, and an adjacent stationary electrode (see  and references therein). The capacitive coupling between the electrodes allows electrostatic control of both the mean position of the resonator and its resonance frequency. In addition to the electrostatic interaction, such electrodes are coupled via the Casimir effect , . This quantum electrodynamical effect originates from the dependence of the ground state energy of the electromagnetic field upon boundary conditions and leads to an observable forces between macroscopic bodies. Due to its relatively short range, this force has only a small effect on the dynamics of macroscopic mechanical systems. However, the Casimir force can play a major role in modern MEMS where typical distances between neighboring surfaces can be on the sub-micron length scale . These attractive electrostatic and Casimir forces give rise to a mechanical metastability that may cause stiction (for a review see , ). In the present case it is manifested as collapse (pull-in) of the movable resonator onto the nearby electrode, resulting in their permanent adhesion. This phenomena can be a principal cause of malfunctioning in MEMS. Moreover, it limits the range of stable operation of CTRs.

In the present work we study experimentally and theoretically the range of tunability and mechanical metastability in a CTR system made of a doubly clamped Au beam separated from an adjacent counterelectrode by a small vacuum gap. Once the beam is brought to contact with the electrode, permanent adhesion occurs. This indicates that the free state of this system is merely metastable, and that the state of contact (after pull-in) has lower energy due to the strongly attractive Casimir interaction. The potential barrier separating these two states determines the lifetime of the free state. The barrier can be reduced by introducing electrostatic attraction, namely by applying a DC voltage, $`V_{dc}`$ between the beam and the electrode. In our measurements we study the mechanical properties of the beam as $`V_{dc}`$ is gradually increased to the point of pull-in. Near this critical point all of the measured properties of the beam show strong dependence on $`V_{dc}`$. Using a simple model we calculate the shape of the potential barrier confining the metastable mechanical state, and then estimate the escape rate via both thermal excitation and quantum tunneling.

Fig. 1(a) is a micrograph displaying a side view of our device. The details of the fabrication process are given elsewhere . The structure is designed to allow full characterization of the beam’s properties to yield straightforward and unambiguous interpretation of our results. We use bulk (rather than surface) micromachining, which allows the substrate to be completely removed beneath the sample. This greatly simplifies the boundary conditions of the electromagnetic field in the vicinity of the sample. Moreover, we avoid using multilayered structures, since their internal stresses are generally important, and are difficult to model theoretically. The beam has length $`l=200`$ $`\mu `$m, width $`a=0.28\mu `$m and thickness $`t=0.25\mu `$m (measured using an atomic force microscope). At its center, the stationary electrode has a 20$`\mu `$m long rectangular protrusion. This is separated from the beam by a vacuum gap of $`g=1.28`$ $`\mu `$m.

All measurements are done at room temperature, in-situ within a commercial scanning electron microscope (SEM). A voltage $`V=V_{dc}+V_{ac}\mathrm{cos}\left(2\pi \nu t\right)`$ is applied between the beam and the electrode (see Fig. 1(b)). The static response of the beam is measured by operating the SEM in scanning mode, imaging the device, and using digital image processing to extract experimental data. The dynamical response is measured with a stationary electron beam focused at a point near the edge of the Au beam. To detect mechanical displacement we use a spectrum analyzer to monitor the output signal from the photomultiplier serving as a secondary electron detector. Note that this detection scheme is sensitive almost exclusively to motion in the plane of the sample.

To characterize the mechanical properties of the beam, namely stress and stiffness, we measure its resonance frequencies. For small stiffness they are given by:

$$\nu _n=n\nu _0\left[1+2\zeta +\left(4+n^2\pi ^2/2\right)\zeta ^2\right],$$
(1)
where $`\zeta ^2=Ea^3t/12Tl^2,`$ $`\nu _0=\sqrt{T/\rho A}/2l`$, with $`E`$ being Young’s modulus, $`T`$ is the stress, and $`\rho `$ is the mass density ,. The dimensionless parameter $`\zeta `$ indicates the relative effect of stiffness compared to stress on the dynamics of the beam. Note that the terms that make the spectrum unequally spaced are of order $`O\left(\zeta ^2\right).`$

We excite the beam by applying AC voltage ($`V_{dc}=0)`$ and find three modes with frequencies $`\nu _1=185.53kHz`$, $`\nu _2=372.4kHz`$, and $`\nu _3=563.8kHz`$. Using Eq. (1) we find $`\zeta `$=0.017$`\pm 0.002`$, thus stiffness has only a small effect on the dynamics of the beam. In what follows, we ignore altogether corrections due to stiffness. The stress is found from the measured resonance frequencies to be $`T=7.4\times 10^6N.`$

Next we gradually increase $`V_{dc}`$ and measure the mechanical properties of the beam until pull-in occurs. Figure 2(a) shows the normalized static deflection of the beam’s center, $`x_1/g`$, as a function of $`V_{dc}`$. After completing the measurements at $`V_{dc}=29`$ V, pull-in occurred with this particular beam. Simultaneously applying the small AC excitation with the DC bias allows us to measure the frequency of the fundamental mode $`\nu _1=\omega _1/2\pi `$. Figure. 2(b) shows the ratio $`\mathrm{\Delta }\omega /\omega _0=\left(\omega _1\omega _0\right)/\omega _0`$, where $`\omega _0`$ is the angular frequency at $`V_{dc}=0`$. Both $`x_1`$ and $`\omega _1`$ show strong dependence on $`V_{dc}`$ near pull-in.

Further insight can be gained by studying the dynamics of the beam in the nonlinear regime. For oscillations with a large amplitude the harmonic approximation breaks down and effects due to higher order terms of the potential are observable. Figure 3 shows the measured response of a beam as a function of frequency around the fundamental resonance for a discrete range of ac drive ($`V_{ac})`$ levels, at fixed DC bias, $`V_{dc}=20`$ V. As $`V_{ac}`$ is increased, the response becomes asymmetric, with a sharp drop occuring on the high frequency side. Moreover, a strong hysteresis is found for high $`V_{ac}`$ when the frequency is scanned back and forth across the resonance (see inset of Fig. 3). Theoretically, the response function in the nonlinear regime is given to lowest order by :

$$a^2\left[\left(\omega \omega _0\left(1+\kappa a^2\right)\right)^2+\left(\omega _0/2Q\right)^2\right]=\left(\frac{f}{2m\omega _0}\right)^2,$$
(3)
where $`f`$ is the amplitude of the external force, $`a`$ is the amplitude of oscillations, $`m`$ is the mass, and $`\kappa `$ is a parameter describing the amplitude-dependent frequency shift due to nonlinearity.

From the shape of the resonance peaks in the nonlinear regime we extract the nonlinear frequency-shift parameter, $`\kappa `$, as a function of $`V_{dc}`$ (Fig. 2(c)). For relatively small $`V_{dc}`$ we find $`\kappa >0`$. For higher $`V_{dc}`$ $`\kappa `$ becomes negative, and close to the point of pull-in $`\kappa `$ drops quite rapidly as a function of $`V_{dc}`$, indicating the strong nonlinearity of the potential in this regime.

We describe our system using a simple one dimensional model as proposed in Ref.  (see Fig. 4(a)). This model employs rather simplistic approximations and thus cannot be expected to provide accurate quantitative predictions – it is extremely useful, however, in providing physical insight and a qualitative description of the dynamics. A mass $`m`$ coupled to a string with a constant $`k=m\omega _0^2`$ represents the doubly-clamped beam. An adjacent surface at distance $`g`$ interacts with the mass electrostatically, as well as by the Casimir effect. For our device geometry it is appropriate to model these forces, $`F`$, as applied locally to the beam’s center. In this case the displacement of the beam’s center is $`Fl/4T`$, and the effective spring constant is thus $`k=4T/l`$.

The total potential energy is given by:

$$U_t\left(x\right)=U_E\left(x\right)+U_V\left(gx\right)+U_C\left(gx\right).$$
(4)
The first term $`U_E\left(s\right)=\left(1/2\right)ks^2`$ is the beam’s elastic potential. The second term is the electrostatic potential introduced by applying a DC voltage, $`V`$, between the beam and its counterelectrode. We estimate the capacitance $`C`$ between the beam and the electrode using a simple parallel plate approximation, thus $`U_V\left(s\right)=ϵ_0AV^2/2s`$, where $`ϵ_0`$ is the vacuum permittivity and $`A`$ is the area. The third term is the Casimir potential. As in the electrostatic case, we employ a parallel plate approximation leading to $`U_C\left(s\right)=\left(\pi ^2/720\right)\left(\mathrm{}cA/s^3\right)`$, where $`\mathrm{}`$ is the Planck’s constant and $`c`$ is the velocity of light . Using these expressions we rewrite $`U_t`$ in a dimensionless form:

$$U_t=U_E\left(g\right)\left[\left(\frac{x}{g}\right)^2\frac{\alpha }{1x/g}\frac{\beta }{1\left(x/g\right)^3}\right],$$
(5)
where $`\alpha =U_V\left(g\right)/U_E\left(g\right)=ϵ_0AV^2/kg^3`$ represents the relative importance of the electrostatic interaction (compared with the elastic term), and similarly $`\beta =U_C\left(g\right)/U_E\left(g\right)=\left(\pi ^2/360\right)\left(\mathrm{}cA/kg^5\right)`$ for the ratio of Casimir interaction to the elastic term.

For relatively small values of $`\alpha `$ and $`\beta `$ (see figure 4(b)) the potential energy has a local minimum at $`x=x_1`$ associated with a metastable state. Figure 4(c) shows a typical example of $`U_t`$ in the metastable regime. The angular frequency of small oscillations around the local minimum is $`\omega _1`$. The barrier that separates the local minimum at $`x=x_1`$ from the global minimum at $`x=g`$ has height $`U_b`$. Below, we calculate the parameters $`\omega _1`$ and $`U_b`$ for two special cases relevant to the present experiment.

Consider first the case of $`\alpha 1`$ and $`\beta 1`$. In this case we find $`x_1/g=\left(\alpha +3\beta \right)/2`$ and $`\omega _1/\omega _0=\left(1\alpha /23\beta \right)`$. Since the escape rate for this case is very small, we do not estimate the associated $`U_b`$.

In the next case we consider the instability occurring with a strong electrostatic and a weak Casimir interaction. Note first that if $`\beta =0`$ then a local minimum and a potential barrier exist only for $`\alpha <8/27`$. At the critical point, $`\alpha =8/27`$, the local minimum and the local maximum points coincide at $`x/g=1/3`$. We therefore expand the potential near that point, assuming $`ϵ8/27\alpha 1`$ ($`ϵ>0`$) and also $`\beta 1`$. To lowest order we find

$`x_1/g`$ $`=`$ $`1/3\left(1/8\right)\sqrt{32ϵ216\beta }`$ (7)
$`\omega _1/\omega _0`$ $`=`$ $`\left(3/4\right)\left(32ϵ216\beta \right)^{1/4}`$ (8)
$`U_b/kg^2`$ $`=`$ $`\left(3/512\right)\left(32ϵ216\beta \right)^{3/2}`$ (9)
The critical point is thus shifted to $`\alpha =8/2727\beta /4`$ due to the Casimir effect. The nonlinear frequency-shift parameter $`\kappa `$ is calculated from the coefficients of the Taylor expansion of the potential around the point $`x=x_1\text{[9]}`$. To lowest order we find $`\kappa g^2\left(32ϵ216\beta \right)^{3/4}`$.

To compare the above results with theory we substitute the area, $`A,`$ with an effective area, $`A_{eff}=2.3A`$, chosen to ensure that the measured and calculated critical points coincide. This substitution serves to account, in part, for the crude approximations employed in the model. The solid lines in Fig. 2(a) and 2(b) show the theoretical results of Eq. (7). Agreement with experiment for both $`x_1`$ and $`\omega _1`$ is moderately good despite the simplicity of the model. Note that our measurements within the nonlinear regime are strongly affected by nonlinearity in the detector’s response. Given that this is not properly accounted for in our model, we do not attempt comparison with theory for the parameter $`\kappa `$.

We now turn our attention to the lifetime of the metastable mechanical state before pull-in. At finite temperatures, thermal excitation above the energy barrier can allow the system to ”escape” from the freely-moving to the pulled-in state. The rate for this generic type of process was first found by Kramers, via the Fokker-Planck equation. The thermal escape rate for the case of small damping is

$$\mathrm{\Gamma }_t=\frac{\omega _1}{Q}\frac{U_b}{k_BT}\mathrm{exp}\left(\frac{U_b}{k_BT}\right).$$
(11)
This result is valid when $`U_b/k_BT1`$ and $`U_b/k_BT1/Q`$. Alternatively, at low temperatures the system escapes via quantum tunneling . Caldeira and Leggett (CL)   have explored what is termed macroscopic quantum tunneling (MQT) in systems with many degrees of freedom. In the case of small damping they find that the MQT rate is

$$\mathrm{\Gamma }_Q=\omega _1\xi _1\sqrt{\frac{U_b}{\mathrm{}\omega _1}}\mathrm{exp}\left[\frac{U_b}{\mathrm{}\omega _1}\left(\xi _2+\xi _3/Q\right)\right].$$
(12)
Here $`\xi _1`$, $`\xi _2`$ and $`\xi _3`$ are positive constants of order unity. The rate for the case of no dissipation ($`1/Q=0`$) coincides with the conventional WKB result for tunneling within a noninteracting, single-particle picture. The CL model yields a description of the effect of dissipation upon MQT involving only a single parameter, $`Q`$, that subsumes all of the precise details of the interaction between the beam and the environment. Equation (12) explicitly demonstrates that the effect of dissipation is to suppress the tunneling rate.

In the present experiment the temperature is relatively high ($`k_BT\mathrm{}\omega _1`$) therefore the rate of thermal escape is much larger compared to the rate of quantum tunneling. However, the opposite situation, namely $`\mathrm{\Gamma }_Q\mathrm{\Gamma }_t`$, may be realized with present-day technology. For high frequency mechanical oscillators (see for example ) the quantum limit, $`k_BT\mathrm{}\omega _1`$, is reachable at cryogenic temperatures. The model employed in the present paper may be useful in designing an experiment to study phenomena in this regime.

Comparison between the simple theory and experiment indicates significant disagreement. We extract the parameter $`32ϵ216\beta `$ from Eq. (5b) and the value $`\omega _1/\omega _0=0.81`$ found experimentally when the pull-in occurs. This parameter is substituted in Eq. (5c) to find $`U_b`$. This yields $`U_b/k_BT10^7`$, a value which, when substitited into Eq. (11), yields a lifetime that is essentially infinite for all practical purposes. Experimentally, however, we observed pull-in after only a few minutes, in striking contradiction with this theoretical estimate. Similar devices, having $`g=`$1.29, 0.71 and 0.65 $`\mu `$m, exhibited pull-in when $`\omega _1/\omega _0=`$ 0.84, 0.987, and 0.995 respectively. In these cases, again, the simple theory is far from explaining what is experimentally observed. We note, in general, that any theory predicting $`U_b=Akg^2\left(\omega _1/\omega _0\right)^B`$ near the critical point, where both $`A`$ and $`B`$ are of order unity, will exhibit similar discrepancy. This is due to the exponential dependence in Eq. (11), in concert with the fact that $`kg^2/k_BT10^71`$.

Our results do not appear to be anomalous. We have analyzed data from previous experiments involving CTR devices made from stress-free silicon, in which pull-in was also observed upon application of a DC voltage between electrodes. Calculating the expected lifetime for the experimental conditions of Ref.  using Eq. (11), we find lifetimes that are in similarly gross contradiction with experimental observations. In principle, slow processes of stress relaxation in both cases may be the underlying mechanism inducing early pull-in . Further investigations are warranted. In our work we have taken care to rule out spurious effects from environmental electrical noise sources that could, in principle, drive the system out of the metastable state. Measured electrical noise was far smaller than the steady bias voltages applied. Also, ”early” pull-in occurred even when the electron beam was blanked, ruling out the possibility that some form of parasitic local charging may be operative.

The present work shows that MEMS can provide an ideal approach to studies of mechanical metastability. Further experimental and theoretical work will elucidate the underlying mechanisms governing the instability of these systems. Such work is important not only from a theoretical point of view, but also from the standpoint of practical device engineering. Moreover, further enhancement of experimental techniques will enable new studies of macroscopic quantum tunneling in mechanical systems.

The authors are grateful to K. Schwab for his assistance in sample fabrication. This research was supported by DARPA MTO/MEMS under grant DABT63-98-1-0012. E.B. gratefully acknowledges suppport from a Rothschild Fellowship and the R.A. Millikan Postdoctoral Fellowship at Caltech.


format: text

text:
# A lattice gas model of II–VI(001) semiconductor surfaces

## Abstract

We introduce an anisotropic two–dimensional lattice gas model of metal terminated II–VI(001) semiconductor surfaces. Important properties of this class of materials are represented by effective NN and NNN interactions, which result in the competition of two vacancy structures on the surface. We demonstrate that the experimentally observed $`c(2\times 2)`$$`(2\times 1)`$ transition of the CdTe(001) surface can be understood as a phase transition in thermal equilibrium. The model is studied by means of transfer–matrix and Monte Carlo techniques. The analysis shows that the small energy difference of the competing reconstructions determines to a large extent the nature of the different phases. Possible implications for further experimental research are discussed.

Two–dimensional lattice gases have served as models of atoms adsorbed to a singular crystal surface, or the terminating layer of such a surface itself, respectively. The interplay of attractive and repulsive short range interactions can result in highly non–trivial features, see e.g.  and references therein. For instance, square lattice systems with infinite NN–repulsion (hard squares) and NNN–attraction display tricritical behavior. At low temperatures a dense, $`c(2\times 2)`$ ordered phase coexists with a disordered phase of low coverage. Here we will investigate a particular model with highly anisotropic attractive and repulsive interactions, which result in a $`c(2\times 2)`$ groundstate, as well. However, this ordering competes with a $`(2\times 1)`$ structure which can prevail locally in the disordered regime.

The model parameters are chosen as to represent certain properties of metal terminated II–VI(001) semiconductor surfaces. This class of materials has attracted considerable attention due to their potential technological relevance in the development of optoelectronic devices, for a recent overview see . Frequently, (001) surfaces serve as substrates for the growth of II–VI crystals  by means of Molecular Beam Epitaxy or Atomic Layer Epitaxy, for instance. Surface reconstructions play an important role in this context and have been the target of experimental studies . In contrast to most III–V materials, II–VI(001) surfaces exhibit a fairly small number of possible reconstructions, which are less complex than their III–V counterparts, in general.

In the following we will mainly address the CdTe(001) surface, see  for a detailed discussion. Apparently, only Cd–terminated (001) surfaces are observed in vacuum . The underlying, complete Te half–layer provides potential Cd–sites which form a simple square lattice. Electron counting rules  and similar considerations  show that the simultaneous occupation of NN–sites in the \[$`1\overline{1}0`$\]–direction (termed the $`y`$–direction in the following) is excluded in the terminating Cd–layer, whereas NN–neighbors along the $`[110]`$–direction (or $`x`$–axis, for short) are possible. Therefore, unless excess Cd is deposited, the surface is characterized by a vacancy structure with a maximum Cd–coverage of $`\theta =1/2`$.

Figure 1 (a) illustrates the structure of the two relevant configurations which satisfy this constraint at $`\theta =1/2`$. The $`c(2\times 2)`$ reconstruction is characterized by a staggered (checkered) occupation of the square lattice sites. In the $`(2\times 1)`$ structure, Cd–atoms arrange in rows along the $`x`$–direction which alternate with rows of vacancies. In principle, the configurations can be transformed into one another by shifting every other column of Cd–atoms by one lattice site.

Density functional (DF) calculations have shown that the surface energies of the two competing structures at $`\theta =1/2`$ and $`T=0`$ differ only by a small amount $`\mathrm{\Delta }E`$, with the $`c(2\times 2)`$ reconstruction having the slightly lower energy. Qualitatively this preference can be unterstood in terms of electron Coloumb interactions, as the distances of neighboring metal atoms are smaller in the $`(2\times 1)`$ arrangement . For ZnSe, a value of $`\mathrm{\Delta }E0.03eV`$ per potential Zn–site is given in . According to , the energy difference is even smaller ($`\mathrm{\Delta }E0.016eV`$) for the CdTe(001) surface.

This factor should play a crucial role in a phase transition which has been studied for CdTe : in vacuum at temperatures below a critical value of about $`T_c=270^oC\pm 10^oC`$, the surface displays a mixed $`c(2\times 2)`$$`(2\times 1)`$ structure with a clear prevalence of the checkered configuration close to (but below) $`T_c`$. Above $`T_c`$, the $`(2\times 1)`$ arrangement of Cd–atoms dominates the surface. The observed coverage is in the vicinity of $`\theta 0.4`$ in both regimes . The situation is complicated by the fact that the material begins to sublimate at about the same temperature $`T_c`$. However, it has been argued that sublimation through step flow would not hinder the surface to achieve an effective equilibrium configuration on terraces .

The aim of our theoretical investigation is to clarify, whether the nature of the above discussed transition can be explained within a thermodynamic equilibrium framework at all, or if non–equilibrium effects should play a crucial role.

The modeling of reconstructions which are characterized by displacement of atoms from their regular lattice positons, usually requires continuous two– or three-dimensional degrees of freedom. A prominent example is the description of W(100) surfaces by XY–models, see e.g.  and references therein. Here, however, reconstruction occurs via the rearrangement of atoms in vacancy structures and a description in terms of occupation variables is appropriate.

We present here a lattice gas model which takes into account important features of the above discussed II–VI(001) surfaces. We will loosely speak of Cd–atoms in the following, without claiming to reproduce particular properties of CdTe faithfully. In fact, the basic structure of the model would be the same for other II–VI(001) surfaces. In our simplifying picture we consider only the terminating Cd–layer, represented by a square lattice of sites $`(x,y)`$ which can be either occupied $`(n_{x,y}=1)`$ or empty $`(n_{x,y}=0)`$. The influence of the underlying crystal structure is accounted for by effective pairwise interactions of atoms. In the $`y`$–direction, an infinite repulsion excludes the simultaneous occupation of NN–sites, i.e. $`n_{x,y}=1`$ always implies $`n_{x,y\pm 1}=0`$. In the $`x`$–direction, an attractive interaction favors the occupation of NN–pairs, the strength of which is denoted by $`J_x<0`$. A competing attractive interaction of diagonal neighbors (NNN) $`J_d<0`$ tends to stabilize the $`c(2\times 2)`$ arrangement of atoms. The total energy of the system is given by

$$H=\underset{x,y}{}n_{x,y}(J_d[n_{x+1,y+1}+n_{x+1,y1}]+J_xn_{x+1,y}\mu ),$$
(1)
where the sum is over all lattice sites and the (effective) chemical potential $`\mu `$ controls the mean coverage $`\theta =n_{x,y}1/2`$. Without loss of generality we can choose $`J_d=1`$ and thus fix the energy scale. Then $`J_x`$ controls the energy difference $`\mathrm{\Delta }E`$ (in units of $`\left|J_d\right|`$) between a perfectly ordered $`c(2\times 2)`$ and a perfect $`(2\times 1)`$ arrangement at $`\theta =1/2`$: $`\mathrm{\Delta }E=\left|\mathrm{\hspace{0.17em}2}+J_x\right|/2`$ (per lattice site). The groundstate of the system is a $`c(2\times 2)`$ ordered configuration with $`\theta =1/2`$, whenever $`J_x>2`$ (and $`\mu >2`$).

The free energy of the system is obtained from the partition function $`Z=_{\left\{n_{x,y}\right\}}e^{\beta H}`$, where the temperature $`T=1/\beta `$ is also measured in units of $`\left|J_d\right|=1`$. The sum is restricted to configurations $`\left\{n_{x,y}\right\}`$ which obey the NN–exclusion in $`y`$–direction. We have applied standard transfer matrix (TM) techniques  to evaluate the logarithm of $`Z_L`$, the partition sum of a system with $`M=L\times N`$ lattice sites in the limit $`N\mathrm{}`$. Strips of width $`L`$ with periodic boundary conditions were aligned with the $`x`$-axis. Hence, only even $`L`$ allow for the perfect $`c(2\times 2)`$ ordering of the groundstate. Note that the TM is of dimension $`2^L\times 2^L`$, but with a much smaller number $`3^L`$ of non–zero elements due to the anisotropic repulsion.

As a first example we consider the model with $`J_x=1.96`$. Figure 1 (b) shows results for strip width $`L=10`$ at different temperatures and constant chemical potential $`\mu =1.96`$. We have evaluated the coverage $`\theta =n_{x,y}=_{x,y}n_{x,y}/M`$ as well as the correlations

$$c_d=\frac{1}{2}n_{x,y}\left(n_{x+1,y+1}+n_{x+1,y1}\right)\text{ and }c_x=n_{x,y}n_{x+1,y}.$$
(2)
These measure the probabilities of finding an occupied NN–pair ($`c_x`$) or NNN–pair ($`c_d`$) of Cd–atoms, i.e. the contribution of $`(2\times 1)`$– or $`c(2\times 2)`$– dominated regions in the system. Coverage and correlations can be obtained from proper derivatives of $`\mathrm{ln}Z_L`$, or, as in the case of $`\theta `$ and $`c_x`$, directly from the relevant eigenvector of the TM .

In addition, Figure 1 (b) displays results of Monte Carlo simulations of a system with $`M=64\times 64`$ sites. In order to achieve reasonably fast equilibration we have applied a rejection-free algorithm , the results are in good agreement with the TM–calculation. In addition to the correlations (2) we determine order parameters which are associated with a perfect $`c(2\times 2)`$ or $`(2\times 1)`$ structure on one of the sublattices:

$$m_{2\times 1}=\frac{1}{M}\underset{x,y}{\overset{x\text{even}}{}}n_{x,y}\text{ and }m_{2\times 2}=\frac{1}{M}\underset{x,y}{\overset{(x+y)\text{even}}{}}n_{x,y}$$
(3)
Large values ($`\theta `$) of these quantities indicate long range order, whereas a homogeneously disordered occupation of the lattice would yield $`m_{2\times 2}=m_{2\times 1}=\theta /2`$. For the sake of breaking the sublattice symmetry, we have initialized the system with $`m_{2\times 2}=\theta `$ for the equilibration dynamics. We have refrained from determining the order parameters within the TM–approach, which would require the introduction of additional staggered fields to the energy function (1). The TM–formalism offers a more suitable method to localize the phase transition .

In the considered example, one observes a sudden drop of the coverage at $`T0.3`$ when $`\mu =1.96`$ is held constant. Simultaneously, the system looses its long range order as indicated by values $`m_{2\times 2}=m_{2\times 1}=\theta /2`$ in the simulations. This is also signaled in the properties of the relevant eigenvector in the TM-analysis . The behavior is consistent with a first order transition, as it was investigated for similar models with isotropic or anisotropic interactions, see e.g.  and references therein.

Here, however, also the NNN–correlation $`c_d`$ decreases rapidly at the coverage drop, while $`c_x`$ displays a sudden increase and $`c_x>c_d`$ in the high temperature regime. This indicates that the phase transition also affects the short range correlations in the system: atoms order in rows of the $`(2\times 1)`$–type without long range order. At $`\theta =1/2`$ the $`c(2\times 2)`$ ordering is always preferred energetically. For significantly smaller coverages, however, the local rearrangement of atoms is possible and can be favorable if $`J_x2J_d`$. Indeed, the degree of the prevalence of $`c_x`$ over $`c_d`$ depends strongly on the actual coverage as will be discussed below.

We have followed the prescription outlined by Bartelt et al.  for estimating the coverage discontinuity and phase boundaries for $`L\mathrm{}`$ from three different strip widths. The results as obtained from $`L=6,8,10`$ are shown in Figure 2 for the models with $`J_x=1.90`$ and $`J_x=1.60`$, i.e. $`\mathrm{\Delta }E=0.05`$ and $`0.2`$, respectively. At low temperatures (III), an ordered phase with $`\theta 1/2`$ coexists with a disordered phase of low coverage. At higher temperatures, the system becomes homogeneously disordered (II) or ordered (I) depending on the coverage. For $`T\mathrm{}`$, we expect the phase boundary (I/II) to approach the $`\theta =1/2`$ axis. In this limit the infinite repulsion should be the only relevant interaction, columns of lattice sites decouple and the system is always disordered. This is in contrast to hard square models with isotropic NN–repulsion, where an extended regime (I) persists for arbitrary temperature.

As an additional characteristics of the system we have determined the line $`T(\theta )`$ where $`c_x=c_d`$ and extrapolated for $`L\mathrm{}`$. Right of the dashed lines in Figure 2, the $`c(2\times 2)`$–structure is prevalent and vice versa. For small coverage, this characteristic line coincides with the boundary (II/III) of the coexistence region. Hence, for a range of coverages, the transition into disorder is accompanied by a simultaneous and discontinuous change of local ordering from $`c(2\times 2)`$ to $`(2\times 1)`$ arrangement of Cd–atoms.

We obtain also a rough estimate of the phase diagram from additional Monte Carlo simulations at constant coverage. For this purpose, we apply a non–local algorithm which exchanges empty with occupied sites according to Kawasaki–like rates . The system is again initialized in an ordered $`c(2\times 2)`$-configuration for equilibration, and a rapid decrease of $`m_{2\times 2}`$ with increasing $`T`$ marks the transition into the homogeneously disordered phase. Figure 2 shows in both diagrams the results for $`M=128\times 128`$, which are in good agreement with the TM–prediction. Within error bars, we obtain the same results by searching for a pronounced maximum in the fluctuations of order parameters, correlations, or energy. Note that this method is not suitable for detecting the transition into the homogeneously ordered region (I): simulations slow down considerably at almost maximal coverage and, furthermore, (I) and (III) become virtually indistinguishable in small systems.

Figure 2 demonstrates the crucial role that the energy difference $`\mathrm{\Delta }E`$ plays for the nature of the phase transition. With increasing $`\mathrm{\Delta }E`$, the tricritical point shifts to smaller coverage and higher temperature. Even more so does the line which separates $`c(2\times 2)`$ from $`(2\times 1)`$ prevalence. This feature might offer a qualitative explanation for the remarkable fact that the $`c(2\times 2)`$$`(2\times 1)`$ transition, which was investigated for CdTe in great detail, has not been found in ZnSe, so far. There, $`\mathrm{\Delta }E`$ is expected to be significantly larger than for CdTe and the region of noticeable $`(2\times 1)`$–dominance should indeed be smaller. Note that in the experimental investigation, integrated HRLEED–peak intensities provide information about local correlations, similar to $`c_x`$ and $`c_d`$, rather than about long range ordering .

In summary, our model offers an interpretation of the $`c(2\times 2)`$$`(2\times 1)`$–transition in CdTe(001) as an equilibrium phase transition. At medium coverage the transition is, with increasing $`T`$, from a coexistence regime into a homogeneously disordered phase. For small enough energy difference $`\mathrm{\Delta }E`$, this phase transition is accompanied inevitably by a rearrangement of the vacancy structure from $`c(2\times 2)`$– to local $`(2\times 1)`$–ordering.

Of course, some of the detailed experimental findings cannot be accounted for in our simple model, see for instance  for particular phenomena related to the relaxation of surface strain. For a more quantitative comparison with experiments, additional information is needed. A precise measurement of $`\theta `$ as a function of the temperature is difficult, but would reveal the path on which the system enters the $`(2\times 1)`$–dominated region in the phase diagram.

In a naive attempt to interpret our results quantitatively one would identify the dimensionless critical temperature (in units of $`\left|J_d\right|=1`$) with $`T_c270^oC`$, thus setting the scale for expressing the energy difference $`\left|\mathrm{\hspace{0.17em}2}+J_x\right|/2`$ in physical units. For example, the model with $`J_x=1.94`$ exhibits the desired transition with $`\theta 0.4`$ at a temperature $`T0.3`$. This would translate into $`\mathrm{\Delta }E0.005eV`$ which is significantly smaller than the value $`(0.03eV)`$ given in . DF–calculations yield $`\mathrm{\Delta }E`$ at $`T=0`$ and the precise effect of higher temperatures on the relation of (free) energies is unknown. Furthermore, recent calculations have shown that the DF–results are very sensitive (up to a factor of about $`2`$) to the number of atomic layers considered in the calculation . Hence, a serious quantitative matching is not feasible unless more reliable estimates of $`\mathrm{\Delta }E`$ become available.

Another open question is, if and how our results for small values of $`\theta `$ can be interpreted in the experimental context. Terminating layers of metal atoms with very low coverage are unstable in vacuum and the next (metal) layer is uncovered, see e.g. . However, the presence of excess group VI atoms might stabilize an effective equilibrium situation with small metal coverage. As a test for this hypothesis we suggest to search for the structural transition of the ZnSe(001) surface under mildly Se–rich conditions.

Our model also opens the possibility to study the shapes and sizes of domains, e.g. the regions of local $`(2\times 1)`$–dominance in the disordered phase. Experimental data is available for the pronounced anisotropy of such domains . Furthermore, we will study the equilibrium shape of isolated islands of atoms and its dependence on the temperature. This should allow for further comparison with experimental results as reported in , for instance.

Acknowledgment: We would like to thank S. Gundel for fruitful discussions. M. Ahr was supported by the Deutsche Forschungsgemeinschaft.


text:
# 1 Introduction

## 1 Introduction

Analytic torsion  arises in a quantum field theoretic context as (the square of) the partition function of Schwarz’s topological field theory  (see  for a detailed review). This has turned out to be an important result in topological quantum field theory; for example it is used to evaluate the semiclassical approximation for the Chern–Simons partition function , which gives a QFT-predicted formula for an asymptotic limit of the Witten–Reshetikhin–Turaev 3-manifold invariant  since this invariant arises as the partition function of the Chern–Simons gauge theory on the 3-manifold . See also  for a review of Schwarz’s topological field theory in a general context, and  for some explicit results in the case of hyperbolic 3-manifolds.

The partition function, $`Z`$, of Schwarz’s topological field theory is a priori a formal, mathematically ill-defined quantity and its evaluation  is by formal manipulations which in the end lead to a mathematically meaningful result: $`Z=\tau ^{1/2}`$ where $`\tau `$ is the analytic torsion of the background manifold. In this paper we show (§2) that there is an alternative formal evaluation of the partition function which results in the trivial answer $`Z=1`$. This result amounts to a novel perspective on analytic torsion: we find that it can be formally written as a certain ratio of volumes of spaces of differential forms which is formally equal to 1 by Hodge duality.

Reidemeister combinatorial torsion (R-torsion)  arises as the partition function of a discrete version of Schwarz’s topological field theory . This is of potential interest if one is to attempt to capture the invariants of topological QFT in a discrete, i.e. combinatorial, setting. In §3 an analogue of the above-mentioned result is derived for combinatorial torsion.

## 2 Schwarz’s topological field theory and analytic torsion

We begin by recalling the evaluation of the partition function

$`Z={\displaystyle \frac{1}{V}}{\displaystyle 𝒟\omega e^{S(\omega )}}`$ (2.1)
of Schwarz’s topological field theory . Here $`V`$ is a normalisation factor to be specified below. The background manifold (“spacetime”) $`M`$ is closed, oriented, riemannian, and has odd dimension $`n=2m+1`$. For simplicity we assume $`m`$ is odd; then the following variant of Schwarz’s topological field theory can be considered : The field $`\omega \mathrm{\Omega }^m(M,E)`$ is an $`m`$-form on $`M`$ with values in some flat $`O(N)`$ vectorbundle $`E`$ over $`M`$. The action functional is

$`S(\omega )={\displaystyle _M}\omega d_m\omega .`$ (2.2)
Here $`d_p:\mathrm{\Omega }^p\mathrm{\Omega }^{p+1}`$ ($`\mathrm{\Omega }^p\mathrm{\Omega }^p(M,E)`$) is the exterior derivative twisted by a flat connection on $`E`$ (which we suppress in the notation), and a sum over vector indices is implied in (2.2) 111Note that (2.2) vanishes if $`m`$ is even.. A choice of metric on $`M`$ determines an inner product in each $`\mathrm{\Omega }^p`$, given in terms of the Hodge operator $``$ by

$`\omega ,\omega ^{}={\displaystyle _M}\omega \omega ^{}`$ (2.3)
Using this, the action (2.2) can be written as $`S(\omega )=\omega ,d_m\omega `$. Let $`\text{ker}(S)`$ denote the radical of the quadratic functional $`S`$ and $`\text{ker}(d_p)`$ the nullspace of $`d_p`$. Then $`\text{ker}(S)=\text{ker}(d_m)`$, and after decomposing the integration space in (2.1) as $`\mathrm{\Omega }^m=\text{ker}(S)\text{ker}(S)^{}`$ the partition function can be formally evaluated to get

$`Z={\displaystyle \frac{V(\text{ker}(S))}{V}}det{}_{}{}^{}((d_m)^2)_{}^{1/4}={\displaystyle \frac{V(\text{ker}(S))}{V}}det{}_{}{}^{}(d_m^{}d_m)_{}^{1/4}`$ (2.4)
(we are ignoring certain phase and scaling factors; see  for these). Here $`V(\text{ker}(S))`$ denotes the formal volume of $`\text{ker}(S)`$. The obvious normalisation choice $`V=V(\text{ker}(S))`$ does not preserve a certain symmetry property which the partition function has when $`S`$ is non-degenerate ; therefore we do not use this but instead proceed, following Schwarz, by introducing a resolvent for $`S`$. For simplicity we assume that the cohomology of $`d`$ vanishes, i.e. $`\text{Im}(d_p)=\text{ker}(d_{p+1})`$ for all $`p`$ ($`\text{Im}(d_p)`$ is the image of $`d_p`$). Then $`S`$ has the resolvent

$`0\mathrm{\Omega }^0\stackrel{d_0}{}\mathrm{\Omega }^1\stackrel{d_1}{}\mathrm{}\mathrm{\Omega }^{m1}\stackrel{d_{m1}}{}\text{ker}(S)0`$ (2.5)
which we use in the following to formally rewrite $`V(\text{ker}(S))`$. The orthogonal decompositions

$`\mathrm{\Omega }^p=\text{ker}(d_p)\text{ker}(d_p)^{}`$ (2.6)
give the formal relations

$`V(\mathrm{\Omega }^p)=V(\text{ker}(d_p))V(\text{ker}(d_p)^{}).`$ (2.7)
The maps $`d_p`$ restrict to isomorphisms $`d_p:\text{ker}(d_p)^{}\stackrel{}{}\text{ker}(d_{p+1})`$, giving the formal relations

$`V(\text{ker}(d_{p+1}))=|det{}_{}{}^{}(d_p)|V(\text{ker}(d_p)^{}).`$ (2.8)
Combining (2.7)–(2.8) we get

$`V(\text{ker}(d_{p+1}))=det{}_{}{}^{}(d_p^{}d_p)_{}^{1/2}V(\mathrm{\Omega }^p)V(\text{ker}(d_p))^1.`$ (2.9)
Now a simple induction argument based on (2.9) and starting with $`V(\text{ker}(S))=V(\text{ker}(d_m))`$ gives the formal relation

$`V(\text{ker}(S))={\displaystyle \underset{p=0}{\overset{m1}{}}}(det{}_{}{}^{}(d_p^{}d_p)_{}^{1/2}V(\mathrm{\Omega }^p))^{(1)^p}.`$ (2.10)
A natural choice of normalisation is now 222This choice can be motivated by the fact that, in an analogous finite-dimensional setting, the partition function then continues to exhibit a certain symmetry property which it has when $`S`$ is non-degenerate .

$`V={\displaystyle \underset{p=0}{\overset{m1}{}}}V(\mathrm{\Omega }^p)^{(1)^p}.`$ (2.11)
Substituting (2.10)–(2.11) in (2.4) gives

$`Z=[{\displaystyle \underset{p=0}{\overset{m1}{}}}det{}_{}{}^{}(d_p^{}d_p)_{}^{\frac{1}{2}(1)^p}]det{}_{}{}^{}(d_m^{}d_m)_{}^{1/4}.`$ (2.12)
These determinants can be given well-defined meaning via zeta-regularisation , resulting in a mathematically meaningful expression for the partition function. As a simple consequence of Hodge duality we have $`det^{}(d_p^{}d_p)=det{}_{}{}^{}(d_{np1}^{}d_{np1})`$, which allows to rewrite (2.12) as

$`Z=\tau (M,d)^{1/2}`$ (2.13)
where

$`\tau (M,d)={\displaystyle \underset{p=0}{\overset{n1}{}}}det{}_{}{}^{}(d_p^{}d_p)_{}^{\frac{1}{2}(1)^p}.`$ (2.14)
This is the Ray–Singer analytic torsion ; it is independent of the metric, depending only on $`M`$ and $`d`$. This variant of Schwarz’s result is taken from ; it has the advantage that the resolvent (2.5) is relatively simple. The cases where $`m`$ need not be odd, and the cohomology of $`d`$ need not vanish, are covered in  (see also  for the latter case). Everything we do in the following has a straightforward extension to these more general settings, but for the sake of simplicity and brevity we have omitted this.

We now proceed to derive a different answer for $`Z`$ to the one above. Our starting point is (2.13)–(2.14) which we consider as a formal expression for $`Z`$, i.e. we do not carry out the zeta regularisation of the determinants. Instead, we use (2.8) to formally write

$`det{}_{}{}^{}(d_p^{}d_p)_{}^{1/2}={\displaystyle \frac{V(\text{ker}(d_{p+1}))}{V(\text{ker}(d_p)^{})}}`$ (2.15)
Substituting this in (2.14) and using (2.7) we find 333This relation is obtained without any restriction on $`m`$, i.e. for arbitrary odd $`n`$.

$`\tau (M,d)={\displaystyle \frac{V(\mathrm{\Omega }^1)V(\mathrm{\Omega }^3)\mathrm{}V(\mathrm{\Omega }^n)}{V(\mathrm{\Omega }^0)V(\mathrm{\Omega }^2)\mathrm{}V(\mathrm{\Omega }^{n1})}}`$ (2.16)
Formally, the ratio of volumes on the r.h.s. equals 1 due to

$`V(\mathrm{\Omega }^p)=V(\mathrm{\Omega }^{np}),`$ (2.17)
which is a formal consequence of the Hodge star operator being an orthogonal isomorphism from $`\mathrm{\Omega }^p`$ to $`\mathrm{\Omega }^{np}`$. (Recall $`\omega ,\omega ^{}=\omega ,\omega ^{}`$ for all $`\omega ,\omega ^{}\mathrm{\Omega }^p`$.) This implies $`Z=1`$ due to (2.13).

The formal relation (2.16) shows that analytic torsion can be considered as a “volume ratio anomaly”: The ratio of the volumes on the r.h.s. of (2.16) is formally equal to 1, but when $`\tau (M,d)`$ is given well-defined meaning via zeta regularisation of (2.14) a non-trivial value results in general.

It is also interesting to consider the case where $`n`$ is even: In this case, using (2.7)–(2.8) we get in place of (2.16) the formal relation

$`{\displaystyle \frac{V(\mathrm{\Omega }^0)V(\mathrm{\Omega }^2)\mathrm{}V(\mathrm{\Omega }^n)}{V(\mathrm{\Omega }^1)V(\mathrm{\Omega }^3)\mathrm{}V(\mathrm{\Omega }^{n1})}}={\displaystyle \underset{p=0}{\overset{n1}{}}}det{}_{}{}^{}(d_p^{}d_p)_{}^{\frac{1}{2}(1)^p}=1`$ (2.18)
The last equality is an easy consequence of Hodge duality and continues to hold after the determinants are given well-defined meaning via zeta regularisation . On the other hand, the ratio of volumes on the l.h.s. is no longer formally equal to 1 by Hodge duality.

## 3 The discrete analogue

Given a simplicial complex $`K`$ triangulating $`M`$ a discrete version of Schwarz’s topological field theory can be constructed which captures the topological quantities of the continuum theory . The discrete theory uses $`\widehat{K}`$, the cell decomposition dual to $`K`$, as well as $`K`$ itself. This necessitates a field doubling in the continuum theory prior to discretisation: An additional field $`\omega ^{}`$ is introduced and the original action $`S(\omega )=\omega ,d_m\omega `$ is replaced by the doubled action,

$`\stackrel{~}{S}(\omega ,\omega ^{})=\left({\displaystyle \genfrac{}{}{0pt}{}{\omega }{\omega ^{}}}\right),\left({\displaystyle \genfrac{}{}{0pt}{}{0}{d_m}}{\displaystyle \genfrac{}{}{0pt}{}{d_m}{0}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{\omega }{\omega ^{}}}\right)=2{\displaystyle _M}\omega ^{}d_m\omega .`$ (3.1)
This theory (known as the abelian BF theory ) has the same topological content as the original one; in particular its partition function, $`\stackrel{~}{Z}`$, can be evaluated in an analogous way to get $`\stackrel{~}{Z}=Z^2=\tau (M,d)`$. The discretisation prescription is now :

$`(\omega ,\omega ^{})`$ $``$ $`(\alpha ,\alpha ^{})C^m(K)\times C^m(\widehat{K})`$ (3.2)
$`\stackrel{~}{S}(\omega ,\omega ^{})`$ $``$ $`\stackrel{~}{S}(\alpha ,\alpha ^{})=\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha }{\alpha ^{}}}\right),\left({\displaystyle \genfrac{}{}{0pt}{}{0}{^Kd_m^K}}{\displaystyle \genfrac{}{}{0pt}{}{^{\widehat{K}}d_m^{\widehat{K}}}{0}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha }{\alpha ^{}}}\right)`$ (3.3)
Here $`C^p(K)=C^p(K,E)`$ is the space of $`p`$-cochains on $`K`$ with values in the flat $`O(N)`$ vectorbundle $`E`$ and $`d_p^K:C^p(K)C^{p+1}(K)`$ is the coboundary operator twisted by a flat connection on $`E`$, with $`C^q(\widehat{K})`$ and $`d_q^{\widehat{K}}`$ being the corresponding $`\widehat{K}`$ objects; $`^K:C^p(K)C^{np}(\widehat{K})`$ and $`^{\widehat{K}}:C^q(\widehat{K})C^{nq}(K)`$ are the duality operators induced by the duality between $`p`$-cells of $`K`$ and $`(np)`$-cells of $`\widehat{K}`$. The $`p`$-cells of $`K`$ and $`\widehat{K}`$ determine canonical inner products in $`C^p(K)`$ and $`C^p(\widehat{K})`$ for each $`p`$, and with respect to these $`^K`$ and $`^{\widehat{K}}`$ are orthogonal maps. (The definitions and background can be found in ; see also  and .) As in §2 we are assuming that $`m`$ is odd and that the cohomology of the flat connection on $`E`$ vanishes: $`H^{}(M,E)=0`$. Then the partition function of the discrete theory, denoted $`\stackrel{~}{Z}_K`$, can be evaluated by formal manipulations analogous to those in §2 (see ) and the resulting expression can be written as either

$`\stackrel{~}{Z}_K=\tau (K,d^K)\text{or}\stackrel{~}{Z}_K=\tau (\widehat{K},d^{\widehat{K}})`$ (3.4)
where

$`\tau (K,d^K)={\displaystyle \underset{p=0}{\overset{n1}{}}}det{}_{}{}^{}(_{p+1}^Kd_p^K)_{}^{\frac{1}{2}(1)^p}`$ (3.5)
and $`\tau (\widehat{K},d^{\widehat{K}})`$ is defined analogously. Here $`_{p+1}^K`$ denotes the adjoint of $`d_p^K`$ (it can be identified with the boundary operator on the $`(p+1)`$-chains of $`K`$). The quantities $`\tau (K,d^K)`$ and $`\tau (\widehat{K},d^{\widehat{K}})`$ coincide; in fact (3.5) is the Reidemeister combinatorial torsion of $`M`$ determined by the given flat connection on $`E`$, and is the same for all cell decompositions $`K`$ of $`M`$ . (This is analogous to the metric-independence of analytic torsion.) Moreover, the analytic and combinatorial torsions coincide , so the discrete partition function in fact reproduces the continuum one:

$`\stackrel{~}{Z}_K=\stackrel{~}{Z}.`$ (3.6)
We now present an analogue of the formal argument which led to $`Z=1`$ in §2. Consider

$`\tau (K,d^K)\tau (\widehat{K},d^{\widehat{K}})={\displaystyle \underset{p=0}{\overset{n1}{}}}det{}_{}{}^{}(_{p+1}^Kd_p^K)_{}^{\frac{1}{2}(1)^p}det{}_{}{}^{}(_{p+1}^{\widehat{K}}d_p^{\widehat{K}})_{}^{\frac{1}{2}(1)^p}.`$ (3.7)
Using the analogues of (2.15) and (2.7) in the present setting,

$`det{}_{}{}^{}(_{p+1}^Kd_p^K)_{}^{1/2}={\displaystyle \frac{V(\text{ker}(d_{p+1}^K))}{V(\text{ker}(d_p^K)^{})}}`$ (3.8)
and

$`V(C^p(K))=V(\text{ker}(d_p^K))V(\text{ker}(d_p^K)^{}),`$ (3.9)
and the corresponding $`\widehat{K}`$ relations, we find an analogue of the formal relation (2.16):

$`\tau (K,d^K)\tau (\widehat{K},d^{\widehat{K}})`$
$`=`$ $`{\displaystyle \frac{V(C^1(K))V(C^3(K))\mathrm{}V(C^n(K))}{V(C^0(K))V(C^2(K))\mathrm{}V(C^{n1}(K))}}{\displaystyle \frac{V(C^1(\widehat{K}))V(C^3(\widehat{K}))\mathrm{}V(C^n(\widehat{K}))}{V(C^0(\widehat{K}))V(C^2(\widehat{K}))\mathrm{}V(C^{n1}(\widehat{K}))}}`$
Formally, the r.h.s. equals 1 due to

$`V(C^p(K))=V(C^{np}(\widehat{K})),`$ (3.11)
which is a formal consequence of the duality operator being an orthogonal isomorphism from $`C^p(K)`$ to $`C^{np}(\widehat{K})`$ (i.e. $`^K\alpha ,^K\alpha ^{}=\alpha ,\alpha ^{}`$ for all $`\alpha ,\alpha ^{}C^p(K)`$). This implies that, formally,

$`\stackrel{~}{Z}_K=\left(\tau (K,d^K)\tau (\widehat{K},d^{\widehat{K}})\right)^{1/2}=1.`$ (3.12)
Thus we see that combinatorial torsion can also be considered as a “volume ratio anomaly” in an analogous way to analytic torsion.

Finally, in the $`n`$ even case it is straightforward to find a combinatorial analogue of the formal relation (2.18) –we leave this to the reader.

Acknowledgements. D.A. is supported by an ARC fellowship. E.P. is supported by funding from FORBAIRT and Trinity College Dublin.


text:
# VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at 𝑧=3.572

## 1. Introduction

VLBI111Very Long Baseline Interferometry observations of high redshift quasars at centimeter wavelengths allow to address two connected topics:

(1) Frequency-dependent properties of the milliarcsecond structures in quasars. An image obtained at a received frequency $`\nu _\mathrm{r}`$ represents the source structure at a higher emitted frequency $`\nu _{em}=\nu _{obs}(1+z)`$ thus facilitating a direct comparison between radio structures at a given frequency in high redshift quasars and their low redshift counterparts observed at a higher frequency $`\nu _{obs}`$.

(2) Cosmological tests using the milliarcsecond structure of QSOs as a “standard” object. As shown in several recent studies, the “apparent angular size – redshift” ($`\theta z`$, Kellermann (1993), Gurvits (1994), Gurvits et al. (1999)) and “proper motion – redshift” ($`\mu z`$, Vermeulen & Cohen (1994), Vermeulen (1996), Kellermann et al. (1999)) relation on milliarcsecond scales have a potential to provide meaningful cosmological information. In this respect, sources at extremely high redshifts are of special interest.

In this paper we present a VSOP (VLBI Space Observatory Programme, Hirabayashi et al. (1998)) observation of the high-redshift quasar PKS 2215+020. The observations have been conducted as part of a VSOP study of extremely high redshift quasars. The complete results of this project will be presented elsewhere (Gurvits et al. (2000)).

The VSOP is a Space VLBI (SVLBI) mission utilizing an orbiting 8–meter antenna deployed on the Japanese satellite HALCA222The Highly Advanced Laboratory for Communication and Astronomy (Hirabayashi et al. (1998)) and the worldwide array of radio telescopes. The satellite has an elliptical orbit, with the apogee and perigee heights at $``$ 21 000 km and $``$ 560 km, respectively, and orbiting period of roughly 6 hours. In each observation, data from the satellite are recorded by a network of up to 5 tracking stations, and subsequently correlated with the data from participating ground radio telescopes. Regular VSOP observations started in September 1997 at 1.6 and 5 GHz. A 1.6 GHz VSOP observation of PKS 2215+020 described here was made on September 14-15, 1997.

PKS 2215+020 is an optically faint ($`m_\mathrm{B}=21.97`$), radio-loud quasar ($`S_{5\mathrm{GHz}}=0.5`$-0.6 Jy; Griffith & White (1993), Gregory et al. (1996)) included in the Parkes half-jansky flat-spectrum sample (Drinkwater et al. (1997)). It has an emission redshift, $`z=3.572`$ (Drinkwater et al. (1997)); a corresponding linear scale of 3.38 $`h^1`$ pc mas-1 (for the Hubble constant $`H_0=100h`$ km s-1Mpc-1 and deceleration parameter $`q_0=0.5`$). The radio spectral index of 2215+020 is $`\alpha _{2.7\mathrm{GHz}}^{5\mathrm{GHz}}=0.15`$ ($`S_\nu \nu ^\alpha `$). Siebert & Brinkmann (1998) have observed X-ray emission from 2215+020 with the ROSAT High Resolution Imager (HRI), and found only marginal evidence for possible elongation along the $`\mathrm{P}.\mathrm{A}.=60`$-$`70^{}`$.

## 2. VSOP observation and data reduction

We observed the quasar PKS 2215+020 at 1.6 GHz, using the HALCA satellite and a global VLBI array of 15 ground-based radio telescopes (Table 1). The observation lasted for 10 hours, with HALCA data recorded for 7 hours. The data were recorded in the VLBA format (Rogers (1995)), with a total observing bandwidth of 32 MHz divided in two intermediate frequency (IF) bands, each having 256 spectral channels. The tracking stations in Green Bank and Goldstone (USA) were used for the HALCA data acquisition. The data were correlated at the VLBA correlator in Socorro (Benson (1995)), with a pre–averaging time of 1.966 and 0.524 sec for the ground–ground and space–ground baselines respectively. Fringe visibilities were detected in the HALCA data recorded at both tracking stations. The resulting sampling function (uv-coverage) of the final correlated dataset is shown in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$, indicating an improvement of resolution by a factor of $``$2.5 compared to ground–based VLBI observations at the same frequency, owing to the ground–space baselines.

We used AIPS333The NRAO Astronomical Image Processing System and the DIFMAP program (Shepherd et al. (1994)) for post-processing of the correlated dataset. The a priori amplitude calibration was determined from the antenna gain and system temperature measurements, including that of the HALCA radio telescope (Moellenbrock et al. (2000)). After inspecting the IF bandpasses, the last 46 channels in each IF were deleted, owing to a cutoff in the filter bandpass. This has reduced the effective observing bandwidth to 26.2 MHz. We corrected the residual delays and rates, using single-band (SB) and multi-band (MB) fringe fitting (Cotton (1995)) with solution intervals of 5 minutes, and accepted all solutions with $`\mathrm{SNR}7`$. To improve the SNR on the space–ground baseline solutions, we fringe fitted the HALCA data, using the source model obtained by imaging the ground–ground VLBI data. After the fringe fitting, the residual phase variations on the ground–ground baselines were within 3, and within $`15^{}`$, on the space–ground baselines. We then averaged over all frequency channels, and calibrated the phases with a point source model (to enable time averaging). Finally, the data were exported into DIFMAP and further time–averaged into 60–second bins. The amplitude and phase errors were re-calculated in DIFMAP, from the scatter in the unaveraged data. The estimated RMS noise on the HALCA baselines was about 4 times higher than on the ground–ground baselines, consistent with the smaller antenna size and higher system temperature of the HALCA radio telescope. The amplitude visibility distribution obtained after the amplitude calibration and fringe–fitting is shown in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$.

## 3. Ground VLBI array and VSOP images of 2215+020

From the final fringe–fitted dataset, we produced two images: a) an image with all data included (VSOP image, hereafter), and b) an image using only the ground–ground baseline data included (GVLBI image, hereafter). In both cases, uniform weighting is used, and gridding weights are scaled inversely with the amplitude errors. The resulting images are shown in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$, and their basic characteristics are presented in Table 2.

The images in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$ show a clear core–jet type morphology. We note that the general milliarcsecond-scale morphology of PKS 2215+020 described here is consistent with the VLBA 15 GHz image (Kellermann et al. (2000)). In the GVLBI image, we identify the core, C, and 7 enhanced emission regions (jet components J1–J7, hereafter) which we will discuss below. The jet extends up to almost 80 mas from the core, and shows a particularly wide, bright section between 45 and 60 mas (150–200 $`h^1`$pc projected distance), suggesting a possible working surface of a young jet propagating in a dense ambient medium.

The basic structures appear to be similar, in both images shown in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$, taking into account the larger noise and scatter of the data on the space–ground baselines, which results in the increased RMS noise and stronger deconvolution effects in the VSOP image. These effects can be countered by using the natural weighting for the VSOP dataset. We find that the resolution and sensitivity of an image of 2215+020 obtained from the naturally weighted VSOP data are practically identical to those of the uniformly weighted ground VLBI data from our observation (which is commonly expected for VSOP observations; Murphy et al. (1999)). Since our goal is to study the jet at highest resolution, we will discuss exclusively the uniformly weighted data represented by the images in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$.

The most interesting difference between the two images shown in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$ is the absence of the component J4 in the VSOP image. The adjacent features J3 and J5 can be identified in both images, even though, in the GVLBI image, they appear to be weaker than J4. It is not immediately clear whether this effect results from unidentified imaging problems, or reflects in fact the true nature of J4 (which may be completely resolved in the VSOP image). Deconvolution effects are visible elsewhere in the two images (for instance, in the VSOP image, J6 appears being split into two features). On the other hand, deconvolution alone is not likely to completely eliminate a feature. In order to provide quantitative criteria for deciding upon this issue, we turn now to modeling the source structure using the interferometric visibility data (Pearson (1995)).

### 3.1. Model fitting the jet structure in 2215+020

To model the jet structure in 2215+020, we first use DIFMAP to fit the ground baseline data by 8 circular Gaussian components (marked in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$). Hereafter, we shall call this model “the GVLBI model”. The GVLBI model is then used as a starting guess for modeling the visibilities in the entire VSOP dataset. The resulting model (“the VSOP model,” hereafter) contains only 7 of the original 8 components. The component corresponding to feature J4 is indeed completely resolved in the VSOP data. In the VSOP model, the formal size of J4 exceeds 20 mas, with the flux of the feature being similar to that in the VLBI model. We take this as evidence in favor of our suggestion that J4 is completely resolved in the VSOP data.

The derived component parameters are listed in Table 3, for both the GVLBI and VSOP models. To estimate the uncertainties of the fitted component parameters, we measure, for each feature, the peak flux density, $`S_{\mathrm{peak}}`$ and its RMS, in the respective image. We then calculate the uncertainties using approximations given by Fomalont (1989).

In several features, both in the GVLBI and VSOP data, the fitted sizes appear to be smaller than the resolution limit. In order to determine which of the fitted values should be considered only as upper limit, we derive, for each feature, its smallest detectable size of a circular Gaussian component in an image with RMS noise $`\sigma _{\mathrm{RMS}}`$:

$$d_{\mathrm{lim}}=\frac{\pi }{4}\left[\frac{b_{\mathrm{maj}}b_{\mathrm{min}}}{\mathrm{ln}(2)\mathrm{ln}(S/\sigma _{\mathrm{RMS}})}\right]^{1/2},$$
(1)
where $`S`$ is the flux density of the feature, and $`b_{\mathrm{maj},\mathrm{min}}`$ are the major and minor axes of the restoring beam, respectively. For those features in which $`d<d_{\mathrm{lim}}`$, we take $`d_{\mathrm{lim}}`$ as the upper limit estimate of the size.

Using the calculated flux densities and sizes, we calculate brightness temperatures, $`T_\mathrm{b}`$, for all jet features.444We note that our choice of Gaussian components for fitting the visibility data (as opposed to, e.g., optically thick spheroids) does not introduce qualitative changes to the results obtained. Models with non–Gaussian brightness distribution would require only a small, constant correction factor to be applied to the values of brightness temperature presented in this paper. For most of the calculated $`T_\mathrm{b}`$, formal errors are of the order of 100% (perhaps indicating that the measured RMS errors of $`S_{\mathrm{peak}}`$ are too conservative). We therefore regard the derived values of $`T_\mathrm{b}`$ as first order estimates. The component sizes and calculated brightness temperatures are plotted in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$.

## 4. Discussion

The VSOP observation shows several very interesting properties of the jet in 2215+020. The observed extent of the jet ($`>80\mathrm{mas},250h^1`$ pc) is remarkable. The structures observed in 2215+020 are almost ten times larger than in any other quasars at $`z>3`$ observed with VLBI (c.f. Paragi et al. (1999) and references therein). The morphology of the source suggests that the component J1 may be a working surface of the jet. From the jet expansion presented in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$, we estimate that the jet opening angle should vary between 2and 7(with the smaller value obtained not taking into account the two most peculiar features: J4, and J1). These values must also be corrected by a deprojection factor, $`\mathrm{sin}\theta _\mathrm{j}^1`$. However, considering the unusually large extent of the jet, the effect of deprojection is not likely to be strong.

### 4.1. Brightness temperature

The observed brightness temperatures plotted in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$ can be used to probe the physics of the jet. Following Marscher (1990), we assume that each of the jet components is an independent plane shock in which the radio emission is dominated by adiabatic energy losses. The jet plasma has a power law energy distribution, $`N(E)dEE^sdE`$, and the magnetic field is assumed to vary as $`Bd^a`$. The Doppler factor is assumed to be constant or varying weakly throughout the jet. Under these assumptions, one can relate the brightness temperature, $`T_{\mathrm{b},\mathrm{J}}`$, of each jet component to the brightness temperature of the core, $`T_{\mathrm{b},\mathrm{C}}`$:

$$T_{\mathrm{b},\mathrm{J}}=T_{\mathrm{b},\mathrm{C}}\left(d_\mathrm{J}/d_\mathrm{C}\right)^ϵ,$$
(2)
where $`d`$ denotes the respective measured sizes of the core and jet feature, and $`ϵ=[2(2s+1)+3a(s+1)/6`$. We take $`s=2.0`$ (corresponding to synchrotron emission with a spectral index $`\alpha =0.5`$) and $`a=1`$ (which corresponds to the transverse orientation of magnetic field in the jet). The measured value, $`T_{\mathrm{b},\mathrm{C}}=9\times 10^{11}`$, is used for the core. We use the sizes and upper limits measured in the VSOP image, and add the GVLBI measurement for J4. The resulting brightness temperatures are indicated in Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$ by thick dashed line. The measured and model values of brightness temperatures agree rather well, suggesting that the jet components may indeed be a collection of relativistic shocks. The discrepancy between the model and observed brightness temperatures is significant only in the outermost jet feature, J1, for which we find a ratio $`\xi T_\mathrm{b}^{\mathrm{obs}}/T_\mathrm{b}^{\mathrm{model}}20`$. This discrepancy suggests that J1 is a region of the jet where physical conditions may have become substantially different. If J1 interacts with the ambient medium, its material may be going through additional compression and reacceleration, which could lead to a localized increase of the brightness temperature. The observed discrepancy between our model and the measured brightness temperature of J1 may also be explained by variations of the jet Doppler factor, resulting from changes of the speed or orientation of the jet. For this case, one can derive for the Doppler factor of J1, $`\delta _{\mathrm{J1}}=\delta _\mathrm{j}\xi ^{1/(\alpha 3)}`$ (here $`\delta _\mathrm{j}`$ is the Doppler factor in the rest of jet). This condition implies a jet bulk Lorentz factor $`\mathrm{\Gamma }_\mathrm{j}\xi ^{1/(3\alpha )}=2.5`$ and viewing angle $`\theta _{\mathrm{jet}}\mathrm{arccos}[(\mathrm{\Gamma }_\mathrm{j}\xi ^{1/(3\alpha )})(\mathrm{\Gamma }_\mathrm{j}^21)^{1/2}]`$. We find then that the Lorentz factor in J1 must be 40–50% lower than in the rest of the jet ($`\mathrm{\Gamma }_{\mathrm{J1}}`$ is lower than $`\mathrm{\Gamma }_\mathrm{j}`$ because $`\delta _\mathrm{j}<1`$, so the jet has to be Doppler–deboosted, and its viewing angle must be $`40^{^{}}`$, even for $`\mathrm{\Gamma }_\mathrm{j}`$ as high as 15). Alternatively, if the bulk Lorentz factor is constant throughout the entire jet, the direction of the flow entering the region of J1 must change by 20–30$`^{^{}}`$ with respect to the line of sight.

### 4.2. The transverse size of the jet

For the feature J4, which is resolved out in the VSOP image, the linear size derived from the GVLBI data is $`3.4`$ mas ($`12h^1`$ pc). If this size is related to the physical transverse dimension of the jet, the mass of the central object can be estimated, assuming that the jet is collimated by the ambient magnetic field of the host galaxy (e.g. Appl & Camenzind 1992, 1993, Beskin & Malyshkin (2000)). On parsec scales, the width of a jet can often be underestimated, owing to limited resolution and dynamic range of VLBI observations. In 2215+020, the jet feature J4 is arguably the best measure of the width of the jet, as this feature is the only one resolved out in the VSOP image, and its size (6.2 mas) measured in the GVLBI image is significantly larger than the respective limiting resolution (2.3 mas). The measured sizes of other jet features (except J1) are most likely affected by the limited resolution and dynamic range of the data, and therefore are likely to be smaller than the true jet width. Finally, J1 is unsuitable for estimating the central black hole mass because it is most likely a strongly perturbed region, actively interacting with the external medium.

For a jet collimated by the ambient magnetic field, $`B_{\mathrm{ext}}`$, of the host galaxy, the mass of the central object, $`M_{\mathrm{BH}}`$, can be related to the width of the jet, $`r_\mathrm{j}`$, so that (Beskin (1997)):

$`M_{\mathrm{BH}}`$ $`0.5r_\mathrm{j}c^2G^1(B_{\mathrm{ext}}/B_\mathrm{g})^{1/2}`$ (3)
$`r_{\mathrm{j},\mathrm{pc}}h^1(B_{\mathrm{ext}}/B_\mathrm{g})^{1/2}10^{13}\mathrm{M}_{_{\mathrm{}}}.`$
where $`G`$ is the gravitational constant, $`R_\mathrm{g}`$ is the Schwarzschild radius of the the central black hole, and $`B_\mathrm{g}`$ is the magnetic field measured at $`R_\mathrm{g}`$. Equation (3) refers to the transverse dimension of the jet measured at distances comparable to the collimation scale (typically expected to be located at 100–1000$`R_\mathrm{g}`$).

A typical galactic magnetic field is $`B_{\mathrm{ext}}10^5`$ G (Beck (2000)), and one can expect to have $`B_\mathrm{g}10^4`$ G (Field & Rogers (1993)). With these estimates, we obtain $`M_{\mathrm{BH}}4\times 10^9h^1\mathrm{M}_{_{\mathrm{}}}`$. Since J4 is located at about 80 pc distance from the jet core, its size is most likely affected by the expansion of the jet. This effect may result in overestimating the mass of the central black hole by $`<30`$% (assuming a jet opening angle of $`4^{^{}}`$ ).555Similar calculations yield $`6\times 10^8\mathrm{M}_{_{\mathrm{}}}`$ for the mass of the central black hole in M87, based on images from ground VLBI observations made at 1.6 GHz (Reid et al. (1989)). This estimate is comparable with the value of $`3\times 10^9\mathrm{M}_{_{\mathrm{}}}`$ deduced from HST observations of M87 (Harms et al. (1994)).

### 4.3. Large scale radio and X-ray structures in 2215+020

The jet direction observed in our VLBI images is consistent with the extension seen in a VLA observation (Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$) on a scale of 7 arcseconds (23$`h^1`$ kpc). Most interestingly, a similar extension is suggested by the ROSAT HRI image (ROSAT archive directory 701900h666ftp://ftp.xray.mpe.mpg.de/rosat/archive/700000/701900h/). Siebert & Brinkmann (1998) argued that the extension of the emission in the ROSAT image may result from the asymmetrical point–spread function of the HRI. In Figure VSOP observation of the quasar PKS 2215+020: a new laboratory for core-jet physics at $`z=3.572`$, we overlay the ROSAT and VLA images of 2215+020. Given the remarkable morphological similarity of the X–ray and radio emission, we tend to conclude that the extension seen in the ROSAT may be real. This would imply that 2215+020 contains an X–ray emitting gas extending out to about 50$`h^1`$ kpc. The X-ray emission may then arise from interaction between the relativistic outflow and a dense circumgalactic gas. It is however not possible to say at present whether the emission on this scales is directly related to the outflow observed on parsec scales. Recent observations with Chandra X–Ray satellite (Marshall et al. 2000) have indicated similar coincidence of radio and X–ray emission distributions in several extragalactic radio sources, most notably in NGC 1275, allowing to speculate that radio lobes may be pushing the X–ray emitting gas outwards. A similar scenario may be at work in 2215+020.

## 5. Summary

The main results from the observations reported in this paper can be summarized as follows:

1. The radio source PKS 2215+020 has a rich core–jet morphology and unusually large size of the jet which can be traced up to 300$`h^1`$ parsecs. This is by far the longest jet observed at redshifts above 3.

2. The high resolving power of the VSOP allows us to measure the transverse width of the jet, which provides a rough estimate for the mass of the central black hole in the object $`M_{\mathrm{BH}}4\times 10^9h^1\mathrm{M}_{_{\mathrm{}}}`$.

3. The sizes and brightness temperatures derived from the model of the source structure are consistent with the adiabatic energy losses in relativistic shocks embedded in the jet.

4. The arcsecond–scale structure detected in PKS 2215+020 with the VLA coincides with the extension observed in the ROSAT X–ray image, implying a possible presence of an X–ray emitting gas extending out to about 50$`h^1`$ kiloparsecs.

The findings summarized above allow us to present the quasar PKS 2215+020 as a promising laboratory for future studies of core–jet physics at very high redshifts.

We gratefully acknowledge the VSOP Project, which is led by the Institute of Space and Astronautical Science (Japan) in cooperation with many agencies, institutes, and observatories around the world. We thank W. Brinkmann and M. Gliozzi for their helpful remarks on the interpretation of the ROSAT data. We would like to thank the anonymous referee for comments and suggestions for improving the manuscript. SF acknowledges support from the Hungarian Space Office, NWO and OTKA. LIG acknowledges partial support from the European Commission TMR Programme, Access to Large–Scale Facilities under contract ERBFMGEXT950012, the European Commission, TMR Programme, Research Network Contract ERBFMRXCT 96–0034 “CERES”. The National Radio Astronomy Observatory is operated by Associated Universities, Inc. under a Cooperative Agreement with the National Science Foundation.


text:
# Higher twist corrections and maxima for DIS on a proton in the high density QCD region

## 1 Introduction

In this letter we continue our study of phenomena associated with gluon saturation in deep inelastic scattering (DIS) following our original study in Ref. . Our goal is to find the kinematic region where the evolution of parton densities can no longer be correctly described by linear evolution equations. To solve this problem analytically one has to find a solution of the non-linear evolution equation which was derived in Refs. 111Strictly speaking this equation is derived for a nucleus target, however there are reasons to assume that it holds for the nucleon as well. and in Ref.  in different approaches. In spite of the fact that the asymptotic solutions are known, a detailed description of the boundary of the saturation region is still lacking. We would like to determine how gluon saturation can be detected in current and future experimental data.

The idea suggested in Ref.  was that there is a breakdown of the Operator Product Expansion as one gets sufficiently close to the saturation region. This enabled us to find the scale $`\widehat{Q}^2`$ where the twist expansion breaks down. Since, at sufficiently large $`Q^2`$, the saturation scale is the only scale of the problem, it should be expressed unambiguously in terms of $`\widehat{Q}^2`$. We refer the reader to our previous publication for more details.

In Ref.  we used the eikonal model to calculate the various structure functions relevant in DIS on nuclei, and then to expand them into the twist series. In this letter we will address the subject of DIS on a nucleon and apply the same approach developed for DIS on nuclei.

The paper is organized as follows: We start with a brief review of our approach without going into details of the calculation (Sec. 2). In Sec. 3, we discuss numerical estimates regarding the nucleon target. Finally, in Sec. 4 we conclude by proposing possible ways to measure the different twists at HERA and LHC.

## 2 The Model

As was pointed out in the introduction we use the eikonal (or Mueller-Glauber) model to calculate shadowing corrections. A complete description of the eikonal approach to shadowing corrections (SC) was given in Ref. . It was shown in Ref.  that the eikonal approximation leads to a reasonable theoretical as well as experimental approximation in the region of not too small $`x_B`$. There are weighty arguments in favor of using the eikonal model to calculate contributions of different twists in various cross sections. This simple eikonal model is a reasonable starting point for understanding how the shadowing corrections occur. The use of the Glauber-Mueller approach was justified for a nuclear target, and it should hold for a nucleon as well.

Following Ref.  and assuming that all dipoles in the virtual photon wave function have the size $`r_t^2=4/Q^2`$ and neglecting the quark masses we obtain for the longitudinal and transverse photon cross sections

$`\sigma _L(\gamma ^{}p)`$ $`=`$ $`2\alpha _{em}\pi R^2{\displaystyle \frac{Q_s^2}{4Q^2}}G_{\mathrm{3\hspace{0.17em}4}}^{\mathrm{3\hspace{0.17em}2}}\left({\displaystyle \frac{Q^2}{Q_s^2}}|{\displaystyle \genfrac{}{}{0pt}{}{0,1;\frac{5}{2}}{1,1,2;\frac{1}{2}}}\right)`$ (1)
$`\sigma _T(\gamma ^{}p)`$ $`=`$ $`4\alpha _{em}\pi R_A^2{\displaystyle \frac{Q_s^2}{4Q^2}}[G_{\mathrm{3\hspace{0.17em}4}}^{\mathrm{3\hspace{0.17em}2}}\left({\displaystyle \frac{Q^2}{Q_s^2}}\right|{\displaystyle \genfrac{}{}{0pt}{}{1,1;\frac{5}{2}}{1,1,1;\frac{1}{2}}}){\displaystyle \frac{1}{2}}G_{\mathrm{3\hspace{0.17em}4}}^{\mathrm{3\hspace{0.17em}2}}\left({\displaystyle \frac{Q^2}{Q_s^2}}\right|{\displaystyle \genfrac{}{}{0pt}{}{0,1;\frac{5}{2}}{1,1,1;\frac{1}{2}}})`$ (2)
$`+{\displaystyle \frac{Q_s^2}{Q^2}}\{{\displaystyle \frac{1}{4}}G_{\mathrm{3\hspace{0.17em}4}}^{\mathrm{3\hspace{0.17em}2}}\left({\displaystyle \frac{Q^2}{Q_s^2}}\right|{\displaystyle \genfrac{}{}{0pt}{}{1,1;\frac{7}{2}}{2,2,2;\frac{1}{2}}}){\displaystyle \frac{1}{2}}G_{\mathrm{3\hspace{0.17em}4}}^{\mathrm{3\hspace{0.17em}2}}\left({\displaystyle \frac{Q^2}{Q_s^2}}\right|{\displaystyle \genfrac{}{}{0pt}{}{1,2;\frac{7}{2}}{2,2,2;\frac{3}{2}}})\}]`$
where $`{}_{3}{}^{}F_{1}^{}`$ is the generalized hypergeometric function, $`G_{\mathrm{3\hspace{0.17em}4}}^{\mathrm{3\hspace{0.17em}2}}`$ is the Meijer function (see Ref.), and we have defined222In our previous work  we used different notation $`4C=Q_s^2`$.

$$Q_s^2(x_B,Q^2)\frac{2\pi \alpha _Sx_BG(x_B,Q^2)}{3R^2},$$
(3)
where $`G(x_B,Q^2)`$ denotes the gluon distribution in the nucleon, $`R`$ the radius of the nucleon, and $`Q_s^2`$ is used as a shorthand for $`Q_s^2(x_B,Q^2)`$ defined in Eq. (3). $`Q_s^2`$ is a saturation scale in the framework of our model.

Diffractive cross sections can be calculated using the simple relations following from the unitarity constraint:

$$\sigma (\gamma ^{}p)_{T,L}^D(Q_s^2/Q^2)=\sigma (\gamma ^{}p)_{T,L}(Q_s^2/Q^2)\frac{1}{2}\sigma (\gamma ^{}p)_{T,L}(2Q_s^2/Q^2)$$
(4)
Formulae Eq. (1), Eq. (2) and Eq. (4) make it possible to estimate the contribution of any order twist to the particular cross section by closing the contour of integration in a Meijer function over the corresponding pole.

We argued in the Ref. , that studying ratios of the cross sections can yield very useful information about the saturation region. While the cross sections receive large corrections from the next-to-leading order of perturbative expansion, their ratios do not. In this paper we will consider two ratios $`\sigma _L/\sigma _T`$ and $`\sigma _L^D/\sigma _T^D`$. Both of these ratios display a remarkable property as a function of photon virtuality $`Q^2`$. At small $`Q^2`$ the ratios vanish since the real photon is transverse. At large $`Q^2`$ they vanish as well, as implied by pQCD. This leads us to suggest that both ratios have maximum $`Q_{max}^2`$. We would expect that this maximum is a certain function of the saturation scale, since at small $`x_B`$ and large $`Q^2`$ this is the only scale important in the process of dipole scattering. So, studying the behaviour of the $`Q_{max}^2`$ with energy we can determine the $`Q_s^2(x_B)`$ as well. Two problems occur. First, the confinement scale affects the position of the maximum of ratios. We can only hope that at large energies (small $`x_B`$) the maximum will occur at sufficiently large $`Q^2`$ where this influence is small, though as we will see below, it cannot be completely excluded.

Second, we do not know the precise form of the function $`Q_s^2=f(Q_{max}^2)`$ which has to be determined using the full twist expansion. We can, however, roughly estimate the form of this function using a simple approach. Namely, we illustrate the most important features (for finding the $`Q_{max}^2`$) of Eq. (1) and Eq. (2) by the following formulae which have the correct asymptotic behaviour for $`Q^2m_f^2<Q_s^2`$ and $`Q^2Q_s^2>m_f^2`$, and which are valid for any model:

$`\sigma _T`$ $``$ $`R^2{\displaystyle \frac{Q_s^2}{Q^2+Q_s^2}},`$
$`\sigma _L`$ $``$ $`R^2{\displaystyle \frac{Q^2}{Q^2+m_f^2}}{\displaystyle \frac{Q_s^2}{Q^2+Q_s^2}},`$
$`\sigma _T^D`$ $``$ $`R^2{\displaystyle \frac{Q_s^2}{Q^2+Q_s^2}},`$
$`\sigma _L^D`$ $``$ $`R^2{\displaystyle \frac{Q^2}{Q^2+m_f^2}}{\displaystyle \frac{Q_s^4}{(Q^2+Q_s^2)^2}},`$
where the logarithmic contributions are neglected. Note, that we explicitly assumed here that at large energies the only relevant scale of the process is $`Q_s^2`$.

The ratios of the cross sections then read

$`R`$ $`=`$ $`{\displaystyle \frac{\sigma _L}{\sigma _T}}={\displaystyle \frac{Q^2}{Q^2+m_f^2}},`$
$`R^D`$ $`=`$ $`{\displaystyle \frac{Q^2}{Q^2+m_f^2}}{\displaystyle \frac{Q_s^2}{Q^2+Q_s^2}}.`$
Taking the derivative of $`R`$ and $`R^D`$ with respect to $`Q^2`$ and equating it to zero, we find that the ratio $`R(Q^2)`$ has no maximum, while the maximum of the ratio $`R^D(Q^2)`$ is at $`Q_{max}^2=\sqrt{m_f^2Q_s^2}`$. There is a maximum of $`R(Q^2)`$ which appears when we take the logarithmic contributions into account ($`R`$ decreases logarithmically with $`Q^2`$) and turns out to be proportional to $`Q_s^2`$. Our simple arguments show that

1. The maximum of $`R`$ is much shallower than that of $`R^D`$ .
2. The maximum of $`R^D`$ is proportional to $`\sqrt{Q_s^2}`$ .
3. Existence of the non-perturbative scale $`m_f^2`$ is crucial for appearance of maxima.

The maxima that we have just discussed in general occur due to the well-known limiting behaviour of the longitudinal and transverse cross sections at small and large $`Q^2`$. What we showed is that at small enough $`x_B`$ those maxima are situated at a scale, which can be expressed unambiguously in terms of the saturation scale.

To check these conclusions, we plot the $`Q_{max}^2`$ and $`Q_s^2`$ versus $`\mathrm{log}1/x_B`$ in Fig. 3(a). Recall, that we obtained formulae Eq. (1) and Eq. (2) by assuming that all dipoles in the virtual photon wave function are of the size $`4/Q^2`$. Thus, to get the correct (up to logarithmic contributions) numerical value of the saturation scale at fixed $`x_B`$ we evaluated the gluon distribution at sufficiently large value of $`Q^2`$ (in fact we took $`Q^2=100`$ GeV2) where to a good approximation it is a constant.

Note, that in making the simple estimations above we did not rely on a particular model. So, we also check our conclusions using the Golec-Biernat–Wusthoff model which is similar to ours, but for which we know exactly what the saturation scale is. Using their $`\widehat{\sigma }`$ we calculated the ratios $`R`$ and $`R^D`$ and found the maxima. The result is shown in Fig. 3(b). Surprisingly, both results agree quite well with our rough estimations333Recently it was shown in Ref.  that the experimental data at small $`x_B`$ indeed exhibit scaling, the total cross section depends on the unique scaling variable $`\tau =Q^2/Q_s^2`$. (estimations of Fig. 3(a) are valid at not too large $`\mathrm{log}Q^2`$). Moreover, the maximum of $`R`$ is proportional to $`Q_s^2`$ in contrast to the maximum of $`R^D`$. This result can also be obtained from our simple estimation by including logarithmic corrections in $`\sigma _T`$.

## 3 Numerical estimations

We now present the results of the numerical calculations. We performed the calculation using GRV’94 parameterization for the gluon structure function $`x_BG(x_B,Q^2)`$ for two reasons. First, we hope that the non-linear corrections to the DGLAP evolution equation are not included in it. Second, it enables us to perform calculations at values of $`Q^2`$ as small as $`0.4`$ GeV2.

We noted in Ref.  that the serious drawback of GRV’94 is the existence of a domain in $`Q^2`$ where it gives the anomalous dimension of $`x_BG(x_B,Q^2)`$ larger than 1. To overcome this we use the following amended gluon structure function

$`x_BG(x_B,Q^2)`$ $`=`$ $`x_BG^{GRV}(x_B,Q^2)\theta (Q^2\stackrel{~}{Q}^2(x_B))`$ (5)
$`+{\displaystyle \frac{Q^2}{\stackrel{~}{Q}^2(x_B)}}x_BG^{GRV}(x_B,\stackrel{~}{Q}^2(x_B))\theta (\stackrel{~}{Q}^2(x_B)Q^2)`$
where $`\theta `$ is a step function, and $`\stackrel{~}{Q}^2(x_B)`$ is a scale at which the anomalous dimension equals 1 (it is a saturation scale within our model accuracy). Eq. (5) can be used as a simplified parameterization of the behavior of the gluon structure function in the saturation region. Our results turn out not to be sensitive to the precise form of the behaviour of $`xG(x_B,Q^2)`$ in the vicinity of the saturation scale.

In order to compare our calculations with the experimental data one has to multiply the theoretical formulae for the cross sections (structure functions) by factor $`0.5`$. This factor stems from the next-to-leading order corrections which can be taken into account by modeling the anomalous dimension as described in Ref. .

To estimate contributions of different twists to various structure functions defined by

$$F(x_B,Q^2)=\frac{Q^2}{4\alpha _{em}\pi ^2}\sigma (\gamma ^{}p),$$

we performed calculations using Eq. (1), Eq. (2) and Eq. (4). Note, that Eq. (1) and Eq. (2) are valid only in the massless limit $`m_f^2=0`$. As was explained above, this is a reasonable assumption as long as we are not concerned with the behavior at small $`Q^2`$. All calculations are made for the THERA kinematical region. The results of calculations for different twist contributions are shown in Fig. 1. Similarly to the nuclear target case the following important effects are seen, when we compare contributions of the first two non-vanishing twists:

1. In all structure functions except $`F_L`$ there is a scale $`\widehat{Q}^2`$ for the particular structure function where these two twist corrections are equal. $`\widehat{Q}^2`$ is largest for $`F_L^D`$;
2. $`\widehat{Q}^2`$ grows as $`x_B`$ decreases;
3. Twist-4 corrections in $`F_L`$ cancels numerically with twist-4 in444This effect was also noted in Ref.  where another, though similar, model for SC was used. $`F_T`$ in the region where their separate contributions are quite large. This is a reason to measure the polarized structure functions instead of the total one.

In analogy with the scale $`Q_{max}^2`$ the scale at which OPE breaks down is expected to be a certain function of the saturation scale, this function differs for the various structure functions. The question of what the analytic expression of this function is warrants further study. We, however, can assert that the scale at which the OPE breaks down is not smaller than $`\widehat{Q}^2`$. The latter can easily be estimated by equating the analytic expressions for twist-2 and twist-4 corrections obtained in Ref. . The result is $`\widehat{Q}^2`$ $`=fQ_s^2`$ up to logarithmic corrections, where the numerical factor $`f`$ is of order unity (e.g. for $`\sigma _T`$, $`f=0.79`$).

In Figs. 2 we plot the ratios $`F_L/F_T(x_B,Q^2)`$ and $`F_L^D/F_T^D(x_B,Q^2)`$ versus $`Q^2`$. The calculations were performed without neglecting the quark mass $`m_f`$. As is expected on theoretical grounds the ratio has maximum which increases as $`x_B`$ decreases. In the kinematical region of THERA it is expected, that $`Q_s^21/x_B`$. The maxima in Fig. 2 occur at $`Q_{max}^2m_f^2`$ which implies that we are sufficiently far from the confinement region. The maximum of the ratio of the diffractive structure functions $`R^D`$ is more obvious, as we saw in sec. 2.

We show in Figs. 3(a) the behavior of $`Q_{max}^2(x_B)`$ as a function of $`x_B`$. The increase is still smaller than predicted theoretically. Our conclusion is that the kinematical region of THERA is intermediate between the linear and saturation regimes.

## 4 Conclusions

The aim of our work is to understand to what extent the high parton density regime of QCD could be detected in DIS experiments at HERA, there the density of partons could be sufficiently large to make SC significant. Using the eikonal model we observed two manifestations of the saturation scale:

1. The position of the maxima $`Q_{max}^2`$ in ratios $`F_L/F_T`$ and $`F_L^D/F_T^D`$ is intimately related to the saturation scale $`Q_s^2(x_B)`$. $`Q_{max}^2`$ shows $`x_B`$ dependence typical for the saturation scale (see Fig. 3). On the contrary, the ratios of structure functions satisfying DGLAP equation have no maxima, as is readily seen from the leading twist formulae given in Ref. .
2. The maximum of the ratio $`F_L^D/F_T^D`$ is much more apparent than that of $`F_L/F_T`$, since the former occurs due to leading twist behaviour of the structure functions, whereas the later is due to logarithmic corrections to this behaviour.
3. The results of our calculations show that there exists a scale at which the twist expansion breaks down since all twist corrections become of the same order. We can use the DGLAP evolution equations only for $`Q^2`$ larger than this scale. The energy ($`x_B`$) dependence of this scale suggests that the experimental data for deep inelastic scattering on a proton will help us separate leading and higher twist contributions. In the case of nucleon deep inelastic scattering such a separation appears to be a rather difficult task and has not yet been performed.

Acknowledgments

We thank Jochen Bartels, Krzysztof Golec-Biernat, Dima Kharzeev and Yura Kovchegov for very fruitful discussions concerning saturation phenomena.

This paper in part was supported by BSF grant # 9800276. The work of L.McL. was supported by the US Department of Energy (Contract # DE-AC02-98CH10886). The research of E.L. and K.T. was supported in part by the Israel Science Foundation, founded by the Israeli Academy of Science and Humanities.

E.G. and E.L. are very grateful to DESY Theory Group for hospitality extended to them during their stay at DESY where this paper was completed.