text:
# Untitled Document

STATISTICAL MECHANICS AND CAPACITY-APPROACHING

ERROR-CORRECTING CODES

Nicolas Sourlas

Laboratoire de Physique Théorique de l’ Ecole Normale Supérieure **UMR 8549, Unité Mixte de Recherche du Centre National de la Recherche Scientifique et de l’ Ecole Normale Supérieure.

24 rue Lhomond, 75231 Paris CEDEX 05, France.

e-mail: sourlas@lpt.ens.fr

Abstract: I will show that there is a deep relation between error-correction codes and certain mathematical models of spin glasses. In particular minimum error probability decoding is equivalent to finding the ground state of the corresponding spin system. The most probable value of a symbol is related to the magnetization at a different temperature. Convolutional codes correspond to one-dimensional spin systems and Viterbi’s decoding algorithm to the transfer matrix algorithm of Statistical Mechanics.

I will also show how the recently discovered (or rediscovered) capacity approaching codes (turbo codes and low density parity check codes) can be analysed using statistical mechanics. It is possible to show, using statistical mechanics, that these codes allow error-free communication for signal to noise ratio above a certain threshold. This threshold depends on the particular code, and can be computed analytically in many cases.

LPTENS 01/ 31

It has been known$`[14]`$ that error-correcting codes are mathematically equivalent to some theoretical spin-glass models. As it is explained in Forney’s paper in this volume, there have been recently very interesting new developments in error-correcting codes. It is now possible to approach practically very close to Shannon’s channel capacity. First came the discovery of turbo codes by Berrou and Glavieux$`[5]`$ and later the rediscovery of low density parity check codes$`[6]`$, first discovered by Gallager$`[7,8]`$, in his thesis, in 1962. Both turbo codes and low density parity check (LPDC) codes are based on random structures. It turns out, as I will explain later, that it is possible to use their equivalence with spin glasses, to analyse them using the methods of statistical mechanics.

Let me start by fixing the notations. Each information message consists of a sequence of $`K`$ bits $`\stackrel{}{u}=\{u_1,\mathrm{},u_K\},u_i=0`$ or $`1`$. The binary vector $`\stackrel{}{u}`$ is called the source-word. Encoding introduces redundancy into the message. One maps $`\stackrel{}{u}\stackrel{}{x}`$ by encoding. $`\stackrel{}{u}\stackrel{}{x}`$ has to be a one to one map for the code to be meaningful. The binary vector $`\stackrel{}{x}`$ has $`N>K`$ components. It is called a code-word. The ratio $`R=K/N`$ which specifies the redundancy of the code, is called the rate of the code. One particularly important family of codes are the so-called linear codes. Linear codes are defined by

$$\stackrel{}{x}=G\stackrel{}{u}$$

$`G`$ is a binary (i.e; its elements are zero or one) $`(N\times K)`$ matrix and the multiplication is modulo two. $`G`$ is called the generating matrix of the code. Obviously by construction all the components $`x_i`$ of a code-word $`x`$ are not independent. Of all the $`2^N`$ binary vectors only $`2^K=2^{NR}`$, those corresponding to a vector $`\stackrel{}{u}`$, are code-words. Codewords satisfy the linear constraints (called parity check constraints) $`H\stackrel{}{x}=0`$ (modulo two), where $`H`$ is a $`(K\times N)`$ binary matrix, called the parity check matrix. The connection with spin variables is straightforward. $`u_i\sigma _i=(1)^{u_i}`$, $`x_iJ_i=(1)^{x_i}`$. It follows that $`u_i+u_j=\sigma _i\sigma _j`$ and

$$J_i=(1)^{_jG_{ij}u_j}=C_{k_1\mathrm{}k_i}^i\sigma _{k_1}\mathrm{}\sigma _{k_i}$$
$`(1)`$
The previous equation defines the “connectivity matrix ” $`C`$ in terms of the generating matrix of the code $`G`$. Similarly one can write the parity check constraints in the form:

$$(1)^{_jH_{lj}x_j}=1M_{k_1\mathrm{}k_l}^lJ_{k_1}\mathrm{}J_{k_l}=1$$
$`(2)`$
This defines the “parity constraint matrix ” $`M`$ in terms of the parity check matrix $`H`$ of the code.

Codewords are sent through a noisy transmission channel and they get corrupted because of the channel noise. If a $`J_i=\pm 1`$ is sent, the output will be different, in general a real number $`J_i^{out}`$. Let us call $`Q(\stackrel{}{J}^{out}|\stackrel{}{J})d\stackrel{}{J}^{out}`$ the probability for the transmission channel’s output to be between $`\stackrel{}{J}^{out}`$ and $`\stackrel{}{J}+d\stackrel{}{J}^{out}`$, when the input was $`\stackrel{}{J}`$. The channel “transition matrix” $`Q(\stackrel{}{J}^{out}|\stackrel{}{J})`$ is supposed to be known. We will assume that the noise is independent for any pair of bits (“memoryless channel”), i.e.

$$Q(\stackrel{}{J}^{out}|\stackrel{}{J})=\underset{i}{}q(J_i^{out}|J_i)$$
$`(3)`$
Communication is a statistical inference problem. Knowing the noise probability i.e. $`q(J_i^{out}|J_i)`$, the code (i.e. in the present case of linear codes knowing the generating matrix $`G`$ or the parity check matrix $`H`$) and the channel output $`\stackrel{}{J}^{out}`$, one has to infer the message that was sent. The quality of inference depends on the choice of the code.

We will now show that there exists a close mathematical relationship between error-correcting codes and theoretical models of disordered systems. To every possible information message (source word) $`\stackrel{}{\tau }`$ we can assign a probability $`P^{source}(\stackrel{}{\tau }|\stackrel{}{J}^{out})`$, conditional on the channel output $`\stackrel{}{J}^{out}`$. Or, equivalently, to any code-word $`\stackrel{}{J}`$ we can assign a probability $`P^{code}(\stackrel{}{J}|\stackrel{}{J}^{out})`$.

Because of Bayes theorem, the probability for any code-word symbol (“letter”) $`J_i=\pm 1`$, $`p(J_i|J_i^{out})`$, conditional on the channel output $`J_i^{out}`$, is given by

$$\mathrm{ln}p(J_i|J_i^{out})=c1+\mathrm{ln}q(J_i^{out}|J_i)=c2+h_iJ_i$$
$`(4)`$
where $`c1`$ and $`c2`$ are constants (non depending on $`J_i`$) and

$$h_i=\frac{1}{2}\mathrm{ln}\frac{q(J_i^{out}|+1)}{q(J_i^{out}|1)}$$
$`(5)`$
It follows that

$$P^{code}(\stackrel{}{J}|\stackrel{}{J}^{out})=c\underset{l}{}\delta (M_{k_1\mathrm{}k_l}^lJ_{k_1}\mathrm{}J_{k_l};1)\mathrm{exp}(\underset{i}{}h_iJ_i)$$
$`(6)`$
where $`c`$ is a normalising constant. The Kronecker $`\delta `$’s enforce the constraint that $`\stackrel{}{J}`$ obeys the parity check equations (Equ. (2) ), i.e. that it is a code-word. The $`\delta `$’s can be replaced by a soft constraint,

$$P^{code}(\stackrel{}{J}|\stackrel{}{J}^{out})=const\mathrm{exp}[u\underset{l}{}M_{k_1\mathrm{}k_l}^lJ_{k_1}\mathrm{}J_{k_l}+\underset{i}{}h_iJ_i]$$
$`(7)`$
where $`u\mathrm{}`$. We now define the corresponding spin Hamiltonian by:

$$H^{code}(\stackrel{}{J})=\mathrm{ln}P^{code}(\stackrel{}{J}|\stackrel{}{J}^{out})=u\underset{l}{}M_{k_1\mathrm{}k_l}^lJ_{k_1}\mathrm{}J_{k_l}+\underset{i}{}h_iJ_i$$
$`(8)`$
This is a spin system with multispin interactions and an infinite ferromagnetic coupling and a random external magnetic field.

Alternatively, one may proceed by solving the parity check constraints $`J_i=C_{k_1\mathrm{}k_i}^i\sigma _{k_1}\mathrm{}\sigma _{k_i}`$ (i.e. by expressing the code-words in terms of the source-words).

$$P^{source}(\stackrel{}{\sigma }|\stackrel{}{J}^{out})=const.\mathrm{exp}(\underset{i}{}h_iC_{k_1\mathrm{}k_i}^i\sigma _{k_1}\mathrm{}\sigma _{k_i})$$
$`(9)`$
where the $`h_i`$’s are given as before. The logarithm of $`P^{source}(\stackrel{}{\sigma }|\stackrel{}{J}^{out})`$,

$$H^{source}(\stackrel{}{\sigma })=\mathrm{ln}P^{source}(\stackrel{}{\sigma }|\stackrel{}{J}^{out})=\underset{i}{}h_iC_{k_1\mathrm{}k_i}^i\sigma _{k_1}\mathrm{}\sigma _{k_i}$$
$`(10)`$
has obviously the form of a spin glass Hamiltonian.

We have given two different statistical mechanics formulations of error correcting codes. One in terms of the souceword probability $`P^{source}`$ and the other in terms of the code-word probability $`P^{code}`$. Because of the one to one correspondence between code-words and source-words, the two formulations are equivalent. In practice however it may make a difference. It may be more convenient to work with $`P^{source}`$ or $`P^{code}`$, depending on the case. For the case of turbo codes (see later) it will be more convenient to define another probability, the “register word” probability.

It follows that the most probable sequence (“word MAP decoding”) is given by the ground state of this Hamiltonian ($`H^{code}`$ or $`H^{source}`$, depending on the case). Instead of considering the most probable instance, one may only be interested in the most probable value $`\tau _i^p`$ of the i’th “bit” $`\tau _i`$$`[9,10,11]`$ (“symbol MAP decoding”). Because $`\tau _i=\pm 1`$, the probability $`p_i`$ for $`\tau _i=1`$ is simply related to $`m_i`$, the average of $`\tau _i`$, $`p_i=(1+m_i)/2`$.

$$m_i=\frac{1}{Z}\underset{\{\tau _1\mathrm{}\tau _N\}}{}\tau _i\mathrm{exp}H(\stackrel{}{\tau })Z=\underset{\{\tau _1\mathrm{}\tau _N\}}{}\mathrm{exp}H(\stackrel{}{\tau })\tau _i^p=\mathrm{sign}(m_i)$$
$`(11)`$
In the previous equation $`m_i`$ is obviously the thermal average at temperature $`T=1`$. It is amusing to notice that $`T=1`$ corresponds to Nishimori’s temperature$`[12]`$.

When all messages are equally probable and the transmission channel is memoryless and symmetric, i.e. when $`q(J_i^{out}|J_i)=q(J_i^{out}|J_i)`$, the error probability is the same for all input sequences. It is enough to compute it in the case where all input bits are equal to one, i.e. when the transmitted code-word is the all zero’s code-word. In this case, the error probability per bit $`P_e`$ is $`P_e=\frac{1m^{(d)}}{2}`$, where $`m^{(d)}=\frac{1}{N}_{i=1}^N\tau _i^{(d)}`$ and $`\tau _i^{(d)}`$ is the symbol sequence produced by the decoding procedure.

This means that it is possible to compute the bit error probability, if one is able to compute the magnetization in the corresponding spin system.

Let me give a simple example of an $`R=1/2`$ “convolutional” code. From the $`N`$ source symbols (letters) $`u_i`$’s we construct the $`2N`$ code-word letters $`x_k^1`$, $`x_k^2`$, $`k=1,\mathrm{},N`$.

$$x_i^1=u_i+u_{i1}+u_{i2},x_i^2=u_i+u_{i2}$$
$`(12)`$
It follows that

$$J_k^1=\sigma _k\sigma _{k1}\sigma _{k2},J_k^2=\sigma _k\sigma _{k2}$$
$`(13)`$
$$C_{i_{k_1}i_{k_2}i_{k_3}}^{(1,k)}=\delta _{k,i_{k_1}}\delta _{k,i_{k_2}+1}\delta _{k,i_{k_3}+2},C_{i_{k_1}i_{k_3}}^{(2,k)}=\delta _{k,i_{k_1}}\delta _{k,i_{k_3}+2}$$
$`(14)`$
The corresponding spin Hamiltonian is

$$H=\frac{1}{w^2}\underset{k}{}J_k^{1,out}\tau _k\tau _{k1}\tau _{k2}+J_k^{2,out}\tau _k\tau _{k2}$$
$`(15)`$
Here I assumed a Gaussian noise. In that case, Equ. (5) reduces to $`h_k=J_k^{out}/w^2`$, where $`w^2`$ is the variance of the noise. This is a one dimensional spin glass Hamiltonian. In fact it is easy to see that convolutional codes correspond to one dimensional spin systems. Their ground state can be found using the $`T=0`$ transfer matrix algorithm. This corresponds to the Viterbi algorithm in coding theory. For symbol MAP (maximum a posteriori probability) decoding, one can use the $`T=1`$ transfer matrix algorithm. This in turn is the BCJR algorithm in coding theory$`[13]`$.

As it is explained in Forney’s paper in this volume, the newly discovered (or rediscovered) capacity approaching codes are based on random constructions. Using the equivalence explained above it has been possible to analyse them using the methods of statistical mechanics.

Gallager’s low density parity check ($`k,l`$) codes are defined by choosing at random a sparse parity check matrix $`H`$ as follows. $`H`$ has $`N`$ columns (we consider the case of code-words of length $`N`$). Each column of $`H`$ has $`k`$ elements equal to one and all other elements equal to zero. Each row has $`l`$ non zero elements. It follows that $`H`$ has $`Nk/l`$ rows and that the rate of the code is $`R=1k/l`$. It follows from equation (8) that Gallager’s $`k,l`$ codes correspond to diluted spin models with $`l`$-spin infinite strength ferromagnetic interactions in an external random field. It turns out that the belief propagation algorithm, used to decode LPDC codes, amounts to an iterative solution of the Thouless Anderson Palmer$`[14]`$ (TAP) equations for spin glasses. A detailed analysis of these codes is presented in Urbanke’s paper in this volume. Low density parity check codes have been analysed using Statistical Mechanics methods by Kabashima Kanter and Saad$`[15,16]`$ in the replica symmetric approximation. More recently Montanari$`[17]`$ was able to establish the entire phase diagram of LDPC codes. For $`k,l\mathrm{}`$ with $`k/l`$ fixed, he showed that $`k,l`$ codes correspond to a random energy model which can be solved without replicas. There is a phase transition in this model, which occurs at a critical value of the noise $`n_c`$. $`n_c`$ separates a zero error phase from a high error phase. $`n_c`$ in this case equals the value provided by Shannon’s channel capacity. For finite $`k`$ and $`l`$ he found an exact one step replica symmetry breaking solution. The location of the phase transition determines $`n_c`$. In this way he computed also for finite values of $`k`$ and $`l`$ the critical value of the noise below which error free communication is possible. A different value of $`n_c`$, $`n_c^{bp}`$ had already being computed by Richardson and Urbanke$`[18]`$ (see Urbanke’s paper in this volume). Richardson and Urbanke compute $`n_c^{bp}`$ by analysing the behaviour of the decoding algorithm, belief propagation in this case. Statistical mechanics provides a threshold $`n_c`$ which in principle is different from $`n_c^{bp}`$. $`n_c`$ is reached by the optimum (but unknown) decoder.

Turbo Codes also have been analysed using statistical mechanics$`[19,20]`$. Turbo Codes are based on recursive convolutional codes. An example of non recursive convolutional code was given in Equ. (12). The corresponding recursive code is given, most conveniently, in terms of the auxiliary bits $`b_i`$, defined below. The $`b_i`$’s are stored in the encoder’s memory registers, that’s why I call $`\stackrel{}{b}`$ the the “register word”.

$$x_i^1=u_i,x_i^2=b_i+b_{i2},b_i=u_i+b_{i1}+b_{i2}$$
$`(16)`$
It follows that the source letters $`u_i`$ are given in terms of the auxiliary “register letters” $`b_i`$

$$u_i=b_i+b_{i1}+b_{i2}$$
$`(17)`$
All additions are modulo two.

To construct a turbo code, one artificially considers a second source word $`\stackrel{}{v}`$, by performing a permutation, chosen at random, on the original code-word $`\stackrel{}{u}`$. So one considers $`v_i=u_{P(i)}`$ where $`j=P(i)`$ is a (random) permutation of the $`K`$ indices $`i`$ and a second “register word” $`c_i`$, $`c_i=v_i+c_{i1}+c_{i2}`$. Obviously

$$v_i=c_i+c_{i1}+c_{i2}=u_j=b_j+b_{j1}+b_{j2},j=P(i)$$
$`(18)`$
Equ. (18) can be viewed as a constraint on the two register words $`\stackrel{}{b}`$ and $`\stackrel{}{c}`$. Finally in the present example, a rate $`R=1/3`$ turbo code, one transmits the $`3K`$ letter code-word $`x_i^1=u_i`$, $`x_i^2=b_i+b_{i2}`$, $`x_i^3=c_i+c_{i2}`$, $`i=1,\mathrm{},K`$. Let’s call, as before,

$$J_i^\alpha =(1)^{x_i^\alpha },\alpha =1,2,3$$

the channel inputs and $`J_i^{out,\alpha }`$ the channel outputs. In the previous, for reasons of convenience, we formulated convolutional codes using the source-word probability $`P^{source}`$ and LDPC codes using the code-word probability $`P^{code}`$. The statistical mechanics of turbo codes is most conveniently formulated in terms of the “register words” probability $`P^{reg}(\stackrel{}{\sigma },\stackrel{}{\tau }|\stackrel{}{J}^{out})`$ conditional on the channel outputs $`\stackrel{}{J}^{out}`$, where $`\tau _i=(1)^{b_i}`$ and $`\sigma _i=(1)^{c_i}`$. The logarithm of this probability provides the spin Hamiltonian

$$H=\frac{1}{w^2}\underset{k}{}J_k^{out,1}\tau _k\tau _{k1}\tau _{k2}+J_k^{out,2}\tau _k\tau _{k2}+J_k^{out,3}\sigma _k\sigma _{k2}$$
$`(19)`$
Because of Equ. (18), the two spin chains $`\stackrel{}{\tau }`$ and $`\stackrel{}{\sigma }`$ obey the constraints

$$\sigma _i\sigma _{i1}\sigma _{i2}=\tau _j\tau _{j1}\tau _{j2},j=P(i)$$
$`(20)`$
(As previously, we have considered the case of a Gaussian noise of variance $`w^2`$.) This is an unusual spin Hamiltonian. Two short range one dimensional chains are coupled through the infinite range, non local constraint, Equ. (20). This constraint is non local because neighboring $`i`$’s are not mapped to neighboring $`j`$’s under the random permutation. It turns out that this Hamiltonian can be solved by the replica method. One finds a phase transition at a critical value of the noise $`n_{crit}`$. For noises less than $`n_{crit}`$, the magnetization equals one, i.e. it is possible to communicate error free. In this respect, turbo codes are similar to Gallager’s LDPC codes. The statistical mechanical models however, are completely different. Let me also mention that, under some reasonable assumptions, the iterative decoding algorithm for turbo codes (turbodecoding algorithm), which I am not explaining here, can be viewed$`[20]`$ as a time discretisation of the Kolmogorov, Petrovsky and Piscounov equation$`[21]`$. This KPP equation has traveling wave solutions. The velocity of the traveling wave, which is computable analytically, corresponds to the convergence rate of the turbodecoding algorithm. The agreement with numerical simulations is excellent.

So the equivalence between linear codes and theoretical models of spin glasses is quite general and we have established the following dictionary of correspondence.

$$\begin{array}{cc}\hfill Errorcorrectingcode& SpinHamiltonian\hfill \\ \hfill Signaltonoise& J_0^2/\mathrm{\Delta }J^2\hfill \\ \hfill MaximumlikelihoodDecoding& Findagroundstate\hfill \\ \hfill Errorprobabilityperbit& Groundstatemagnetization\hfill \\ \hfill Sequenceofmostprobablesymbols& magnetizationattemperatureT=1\hfill \\ \hfill ConvolutionalCodes& Onedimentionalspinglasses\hfill \\ \hfill Viterbidecoding& T=0Transfermatrixalgorithm\hfill \\ \hfill BCJRdecoding& T=1Transfermatrixalgorithm\hfill \\ \hfill GallagerLDPCcodes& Dilutedpspinferromagnetsinarandomfield\hfill \\ \hfill TurboCodes& CoupledspinchainsPC\hfill \\ \hfill Zeroerrorthreshold& Phasetransitionpoint\hfill \\ \hfill Beliefpropagationalgorithm& IterativesolutionofTAPequations\hfill \end{array}$$

I would like to conclude by pointing out some open questions.

What is the order of the phase transition? This question is particularly relevant for turbo codes and has important implications for decoding.

What are the finite size effects? This question is particularly relevant near the zero error noise threshold (i.e. near the phase transition). The answer will depend on the order of the transition.

How does the decoding complexity behave as one approaches the zero error noise threshold? Is there a critical slowing down? As it was said before, the decoding algorithms both for LDPC codes and turbo codes are heuristic and there are not known results as one approaches the phase transition.

Is there a glassy phase in decoding? In other terms, do the heuristic decoding algorithms reach the threshold of optimum decoding, computed by statistical mechanics, or is there a (lower) noise “dynamical” threshold where decoding stops reaching optimal performance?

I hope that at least some of the above questions will be answered in the near future.

References

1) Sourlas, N., Nature 339, 693 (1989)

2) Sourlas, N., in Statistical Mechanics of Neural Networks, Lecture Notes in Physics 368, ed. L. Garrido, Springer Verlag (1990)

3) Sourlas, N., Ecole Normale Supérieure preprint (April 1993)

4) Sourlas, N., in From Statistical Physics to Statistical Inference and Back, ed. P. Grassberger and J.-P. Nadal, Kluwer Academic (1994) p. 195.

5) C. Berrou, A. Glavieux, and P.Thitimajshima. Proc.1993 Int.Conf.Comm. 1064-1070

6) MacKay, D. J. C. Neal, R. M. Elect. Lett. 33, 457 (1997).

7) Gallager, R. G. IRE Trans. Inform. Theory , IT-8, 21 (1962).

8) Gallager, R. G. Low-Density Parity-Check Codes , MIT Press, Cambridge MA (1963).

9) Ruján, P., Phys. Rev. Lett. 70, 2968 (1993)

10) Nishimori, H., J. Phys. Soc. Jpn. , 62, 2973 (1993)

11) Sourlas, N., Europhys. Lett. 25, 169 (1994)

12) Nishimori, H., Progr. Theor. Phys. 66, 1169 (1981)

13) L. Bahl, J. Cocke, F. Jelinek, and J. Raviv. IEEE Trans.Inf.Theory IT-20(1974) 248-287

14) Thouless, D. J. Anderson, P. W. Palmer, R. G. Phil. Mag.  35, 593 (1977)

15) Kanter, I. and Saad, D. Phys. Rev. Lett. 83, 2660 (1999)

16) Kabashima, Y. Murayama T. and Saad, D. Phys. Rev. Lett. 84, 1355 (2000)

17) Montanari, A. cond-mat/0104079

18) Richardson, T. J. Urbanke, R. L. IEEE Trans. Inform. Theory  47, 638 (2001).

19) Montanari, A. Sourlas, N. Eur. Phys. J.  B 18, 107 (2000)

20) Montanari, A. Eur. Phys. J.  B 18, 121 (2000)

21) Kolmogorov, A. Petrovsky, I and Piscounov, N. Moscou Univ. Math. Bull.  1, 1 (1937).


format: text

text:
# Edge waves along a sloping beach

## 1 Introduction

Standing on a gently sloping straight beach, it is a matter of observation that various waveforms propagate on the surface of the sea. Among these we find edge waves - water waves that progress along the shoreline. These waves, often difficult to visualize (this has, probably, prevented the regarding of this waveform as important for a long time), are coastal trapped, i.e. their amplitude is maximal at the shoreline and decays rapidly offshore. They produce on the beach beautiful run-up patterns (highest points reached by a wave on the beach). Although propagation is along the straight shoreline and the waveform is sinusoidal in the longshore, these waves are not one-dimensional (Lighthill 1978).

While they were originally considered to be a curiosity (Lamb 1932), edge waves are now recognized to play a significant role in nearshore hydrodynamics. For shallow beaches empirical evidence shows that incident storm waves loose most of their energy trough wave breaking by the time they reach the shore. After breaking offshore, as the waves progress to shallower water, their height decreases reaching its least value at the shoreline. Since storms often result in pronounced shoreline erosion, the surf zone water processes with the onset of a storm are dominated by wave conditions other than the incident waves - role attributed by oceanographers to the edge waves ; see Komar (1998). There are other instances when edge waves are of significance. For example, processed data from the water waves created by an earthquake occuring in April 1992 in the ocean floor near the Californian coast show that two distinct wave packets (both directly generated in the nearshore by the vertical motion of the ocean bottom) were recorded at a coastal station about $`150km`$ far from the epicentre. At first, less than an hour after the occurence of the earthquake, relatively fast-moving swell with amplitude around $`15cm`$ striked from offshore. About two hours after the swell has subsided, relatively slow-moving edge waves with amplitudes around $`50cm`$ (Boss et al. 1995) arrived. Measurements performed on this occasion confirm the rapid decay of the amplitude of the edge waves: at an offshore distance of $`12km`$ the amplitude is reduced to $`10\%`$ of its maximal value (attained at the shoreline). Let us also mention that it has been observed (see Evans 1988) that hurricanes travelling approximately parallel to a nearby coastline sometimes give rise to edge waves. Interestingly, edge waves can in fact be generated directly in a laboratory wave tank (Yeh 1985).

The edge wave phenomenon has been extensively studied and discussed in the mathematical literature within the framework of linear theory. Due to the small displacements associated with these waves, the governing equations for water waves or the shallow-water equations are linearized (Minzoni & Whitham 1977 showed that both approximations are equally consistent) and this simplification permits a thorough analysis. We refer to Ehrenmark (1998) for an excellent up-to-date survey.

Despite the fact that the linearizing approximation lacks rigorous mathematical justification, it has been used with considerable success as a large variety of theoretical studies are confirmed in experimental contexts. The investigation of nonlinear edge waves can be seen as a natural extension to the linear theory. Whitham (1976) showed the existence of irrotational weakly nonlinear edge waves that propagate parallel to the shore using a formal Fourier series expansion for the full water-wave theory. A study of properties of nonlinear progressive edge waves based on the fact that the evolution is described by the nonlinear Schrödinger equation was performed by Yeh (1985). This paper describes an alternative approach; the main impetus for the results reported here comes from the belief that the need for a more rigorous theory remains thoroughly justified. A quest for an explicit edge wave solution for the governing equations appears to be of interest since the structure of the edge waves obscures their visual observation; moreover, a solution in closed mathematical form provides a background against which certain features which have been observed (and predicted) can be checked. It turns out that the deep water wave solution discovered by Gerstner (1809), can be adapted to construct edge waves propagating along a plane-sloping beach. This possibility was pointed out by Yih (1966) but the treatement therein, in essence followed also by Mollo-Christensen (1982), provides only an implicit form for the free water surface. We present a procedure by which exact edge wave solutions to the full water-wave equations are obtained. The closed form of the solution in Lagragian (material) coordinates permits us to provide clear illustrations of the structure of these edge waves. From an examination of the solution we also obtain the run-up pattern, an attractive feature being the occurence of cusps - these shapes are confirmed in both field and laboratory evidence (see Komar 1998). Thus, we establish with rigour the existence of rotational nonlinear edge waves, unraveling the detail structure of the wave pattern.

## 2 The edge wave

We take a plane beach and adopt a coordinate system as shown below, with the shoreline the $`x`$-axis, and the still sea in the region

$$R=\{(x,y,z):x,yb_0,\mathrm{\hspace{0.17em}0}z(b_0y)\mathrm{tan}\alpha \}$$

for some $`b_00`$; here $`\alpha (0,{\displaystyle \frac{\pi }{2}})`$ defines the uniform slope.

Let $`u=(u_1,u_2,u_3)`$ be the velocity field and let us recall (e.g. Crapper 1984) the governing equations for the propagation of gravity water waves when ignoring viscous effects. Homogeneity (constant density $`\rho `$) is a good approximation for water (see the numerical data in Lighthill 1978) so that we have the equation of mass conservation in the form

(1)
$$\frac{u_1}{x}+\frac{u_2}{y}+\frac{u_3}{z}=0.$$

The equation of motion is Euler’s equation,

(2)
$$\{\begin{array}{c}\frac{Du_1}{Dt}=\frac{1}{\rho }\frac{P}{x},\hfill \\ \frac{Du_2}{Dt}=\frac{1}{\rho }\frac{P}{y}g\mathrm{sin}\alpha ,\hfill \\ \frac{Du_3}{Dt}=\frac{1}{\rho }\frac{P}{z}g\mathrm{cos}\alpha .\hfill \end{array}$$

where $`P(t,x,y,z)`$ denotes the pressure, $`g`$ is the gravitational acceleration constant and $`D/Dt`$ is the material time derivative, $`{\displaystyle \frac{Df}{Dt}}={\displaystyle \frac{f}{t}}+u_1{\displaystyle \frac{f}{x}}+u_2{\displaystyle \frac{f}{y}}+u_3{\displaystyle \frac{f}{z}}`$, expressing the rate of change of the quantity $`f`$ associated with the same fluid particle as it moves about. The boundary conditions which select the water-wave problem from all other possible solutions of the equations (1)-(2) are (see Johnson 1997):

(i) the dynamic boundary condition $`P=P_0`$ at the free surface, where $`P_0`$ is the constant atmospheric pressure, decouples the motion of the air from that of the water;

(ii) the kinematic boundary condition at the free surface expresses the fact that the same particles always form the free water surface;

(iii) the kinematic boundary condition at the bottom, requiring the normal velocity component at the bed to be zero so that it is impossible for water to penetrate.

The general description of the propagation of a water wave is encompassed by the equations (1)-(2) and the three boundary conditions (i)-(iii), a distinctive feature being that the free surface is not known and must be determined as part of the solution.

We adopt the Lagrangian point of view by following the evolution of individual water particles. We suppose that the position of a particle at time $`t`$ is given by

(3)
$$\{\begin{array}{c}x=a\frac{1}{k}e^{k(bc)}\mathrm{sin}\left(ka+\sqrt{gk\mathrm{sin}\alpha }t\right),\hfill \\ y=bc+\frac{1}{k}e^{k(bc)}\mathrm{cos}\left(ka+\sqrt{gk\mathrm{sin}\alpha }t\right),\hfill \\ z=c+c\mathrm{tan}\alpha +\frac{\mathrm{tan}\alpha }{2k}e^{2kb_0}\left(1e^{2kc(1+\mathrm{cot}\alpha )}\right),\hfill \end{array}$$

where $`k>0`$ is fixed. It should be pointed out that the quantities $`a,b,c`$ do not stand for the initial coordinates of a particle, they are simply labeling variables serving to identify a particle. We may think of them as parameters which fix the position of a particular particle before the passage of the wave (in still water), despite the fact that the wave is not developing from the still state - otherwise, the flow would be irrotational in view of Helmholtz’s theorem (see Johnson 1997) but its vorticity is nonzero (see the last section). Let us explain the origin of (3). Gerstner (1809; see Constantin 2001) gave the only known nontrivial explicit solution to the full water-wave equations, showing that the two-dimensional particle motion ($`a,bb_00,k>0`$)

$$t(a+\frac{1}{k}e^{kb}\mathrm{sin}(ka+\sqrt{gk}t),b\frac{1}{k}e^{kb}\mathrm{cos}(ka+\sqrt{gk}t))$$

represents waves of finite amplitude in water of infinite depth. This suggests that it might be possible to construct an edge wave using an approach similar to the one for a Gerstner wave field. While the theoretical correctness of this conclusion was established by Yih (1966) and Mollo-Christensen (1982), the outcome in both treatements was an implicit form the water’s free surface which makes the obtained waveform graphically and computationally inaccessible. The closed form (3) provides the full details of the edge wave motion without considerable labour.

Our aim is to prove that the motion (3) is dynamically possible and that we can associate to it an expression for the hydrodynamical pressure $`P`$ such that the governing equations and boundary conditions are all satisfied. The resulting free surface of the water will be the edge wave we are looking for.

The map (3) is a diffeomorphism from the still water region $`R`$ to the water region bounded below by the rigid bed $`\{z=0\}`$ and above by the free water surface parametrized by

(4)
$$\{\begin{array}{c}x=a\frac{1}{k}e^{kb(1+\mathrm{tan}\alpha )kb_0\mathrm{tan}\alpha }\mathrm{sin}\left(ka+\sqrt{gk\mathrm{sin}\alpha }t\right),\hfill \\ y=b(1+\mathrm{tan}\alpha )b_0\mathrm{tan}\alpha +\frac{1}{k}e^{kb(1+\mathrm{tan}\alpha )kb_0\mathrm{tan}\alpha }\mathrm{cos}\left(ka+\sqrt{gk\mathrm{sin}\alpha }t\right),\hfill \\ z=(b_0b)(1+\mathrm{tan}\alpha )\mathrm{tan}\alpha +\frac{\mathrm{tan}\alpha }{2k}e^{2kb_0}\left(1e^{2k(bb_0)(1+\mathrm{tan}\alpha )}\right),\hfill \end{array}$$

with $`a`$, $`bb_0`$, and $`t0`$. Indeed, observing that $`(a,b,c)(a,bc,c)`$ defines a diffeomorphism of $`^3`$, it is enough to show that the map

$$\left(\begin{array}{c}a\\ b^{}\\ c\end{array}\right)\left(\begin{array}{c}a\frac{1}{k}e^{kb^{}}\mathrm{sin}\left(ka+\sqrt{gk\mathrm{sin}\alpha }t\right)\\ b^{}+\frac{1}{k}e^{kb^{}}\mathrm{cos}\left(ka+\sqrt{gk\mathrm{sin}\alpha }t\right)\\ c+c\mathrm{tan}\alpha +\frac{\mathrm{tan}\alpha }{2k}e^{2kb_0}\left(1e^{2kc(1+\mathrm{cot}\alpha )}\right)\end{array}\right)$$

is a diffeomorphism on $`\times _{}\times _+`$. To see this, observe that the third coordinate depends only on $`c`$, being an increasing function $`f(c)`$ with $`f(0)=0`$ and $`lim_c\mathrm{}f(c)=\mathrm{}`$. Therefore, we may slice $`\times _{}\times _+`$ by planes parallel to the plane $`c=0`$ and in each such plane the particle motion is precisely that of a Gerstner wave field (see Constantin 2001). This proves that (3) is a diffeomorphism and it is easy to identify the boundary of the image of the region $`R`$ under it.

The Lagrangian form of the equation of continuity (that is, the volume-preserving property of the flow) is fulfilled since the value of the Jacobian of the map (3) is independent of time. This, together with the previously proved fact that (3) defines at any fixed time a diffeomorphism, shows that the motion described by (3) is dynamically possible. To complete the proof that (3) describes the water motion induced by a gravity wave, we have to check Euler’s equation (2) and the boundary conditions (i)-(iii) for a suitably defined value of the hydrodynamical pressure.

The acceleration of a particular water particle is

$$\frac{Du}{Dt}=(g\mathrm{sin}\alpha e^{k(bc)}\mathrm{sin}(ka+\sqrt{gk\mathrm{sin}\alpha }t),g\mathrm{sin}\alpha e^{k(bc)}\mathrm{cos}(ka+\sqrt{gk\mathrm{sin}\alpha }t),\mathrm{\hspace{0.17em}0})$$

so that the equation of motion (2) is

$$\{\begin{array}{c}\frac{P}{x}=\rho g\mathrm{sin}\alpha e^{k(bc)}\mathrm{sin}(ka+\sqrt{gk\mathrm{sin}\alpha }t),\hfill \\ \frac{P}{y}=\rho g\mathrm{sin}\alpha e^{k(bc)}\mathrm{cos}(ka+\sqrt{gk\mathrm{sin}\alpha }t)\rho g\mathrm{sin}\alpha ,\hfill \\ \frac{P}{z}=\rho g\mathrm{cos}\alpha .\hfill \end{array}$$

Passing to Lagrangian coordinates, we obtain the system

$$\{\begin{array}{c}\frac{P}{a}=0,\hfill \\ \frac{P}{b}=\rho g\mathrm{sin}\alpha e^{2k(bc)}\rho g\mathrm{sin}\alpha ,\hfill \\ \frac{P}{c}=\rho g\mathrm{sin}\alpha e^{2k(bc)}\rho g\mathrm{cos}\alpha +\rho g\mathrm{cos}\alpha (1+\mathrm{tan}\alpha )e^{2kb_0}e^{2kc(1+\mathrm{cot}\alpha )},\hfill \end{array}$$

with the solution

$$P=P_0+\frac{\rho g\mathrm{sin}\alpha }{2k}e^{2k(bc)}\rho g\left(c\mathrm{cos}\alpha +(bb_0)\mathrm{sin}\alpha \right)\frac{\rho g\mathrm{sin}\alpha }{2k}e^{2kc(1+\mathrm{cot}\alpha )}e^{2kb_0}.$$

The obtained hydrodynamical pressure has the same value for any given particle as it moves about. At the free surface $`c=(b_0b)\mathrm{tan}\alpha `$ we have $`P=P_0`$ so that the dynamic boundary condition (i) is satisfied. The kinematic boundary condition at the free surface, (ii), holds as at any instant the free surface (4) is the image under (3) of the still water surface $`\{c=(b_0b)\mathrm{tan}\alpha :bb_0\}`$. That there is no velocity normal to the sloping shore - this takes care of the boundary condition (ii) - is obvious, because at $`z=0`$ we have $`c=0`$ and the motion (3) is planar, without any velocity component in the direction of $`z`$. The proof of the fact that (3) is an explicit solution to the governing equations for water waves on a plane-sloping beach is complete.

## 3 Discussion

We have presented an exact edge wave solution to the full water-wave problem, the graphical depiction of which is a fairly easy exercise. Let us now emphasize some of its significant properties. We present some simple observations which will provide a comprehensive description of this nonlinear wave and the particle motion it induces below the water surface.

The wavelength in the longshore direction $`\lambda =2\pi /k`$ is related to the wave frequency $`\omega `$ by

$$\omega ^2=gk\mathrm{sin}\alpha ,$$

while the wave period is

$$T=\frac{2\pi }{\sqrt{gk\mathrm{sin}\alpha }}.$$

We easily infer that

$$\lambda =\frac{gT^2}{2\pi }\mathrm{sin}\alpha ,$$

so that the length of the edge wave is strongly dependent on its period and to a smaller degree on the beach slope. The phase velocity $`U`$ of the edge wave (4) is given by

(5)
$$U=\sqrt{\frac{g\mathrm{sin}\alpha }{k}},$$

a fact consistent with the observation that if the bottom is flat ($`\alpha =0`$), then $`U=0`$ and no edge wave exists. The dispersion relation (5) for edge waves is obtained (see Johnson 1997) within the confines of the formal linear approximation to the governing equations, but in our case the relation is derived rigorously as a byproduct of (3).

From (3) it is clear that any water particle describes circles as the edge wave passes - all these circles lie in planes parallel to the sloped bottom. The radius $`{\displaystyle \frac{1}{k}}e^{k(bc)}`$ of the circle described counterclockwise by a particle is maximal for the particles at the shoreline (that is, for $`b=b_0,c=0`$).

As pointed out in the previous section, the motion of the water body induced by the passage of the edge wave (4) is rotational. The vorticity of the water flow defined by (3) is

$$\text{curl}u=(\frac{u_2}{z},\frac{u_1}{z},\frac{u_2}{x}\frac{u_1}{y}),$$

by the vanishing of $`u_3`$. Computing the inverse of the Jacobian matrix of the diffeomorphism (3) as

$$\left(\begin{array}{ccc}\frac{1+\mathrm{exp}[k(bc)]\mathrm{cos}k(a+Ut)}{1\mathrm{exp}[2k(bc)]}& \frac{\mathrm{exp}[k(bc)]\mathrm{sin}k(a+Ut)}{1\mathrm{exp}[2k(bc)]}& 0\\ \frac{\mathrm{exp}[k(bc)]\mathrm{sin}k(a+Ut)}{1\mathrm{exp}[2k(bc)]}& \frac{1\mathrm{exp}[k(bc)]\mathrm{cos}k(a+Ut)}{1\mathrm{exp}[2k(bc)]}& 0\\ 0& \frac{1}{(1+\mathrm{tan}\alpha )(1+\mathrm{exp}[2k(b_0cc\mathrm{cot}\alpha )])}& \frac{1}{(1+\mathrm{tan}\alpha )(1+\mathrm{exp}[2k(b_0cc\mathrm{cot}\alpha )])}\end{array}\right),$$

a straightforward calculation yields the expression of the vorticity

$$\text{curl}u=(0,\mathrm{\hspace{0.17em}0},\frac{2kU}{1e^{2k(bc)}}e^{2k(bc)})$$

at a particle whose parameters are $`(a,b,c)`$. Therefore the vorticity is in the opposite sense to the revolution of the particles in their circular orbits, decreasing rapidly with distance from the shoreline/bed. Observe that, despite the fact that the flow (3) is not two-dimensional, the vorticity of each individual water particle is conserved as the particle moves about.

The run-up pattern is obtained by setting $`z=0`$ in (4); this forces $`b=b_0`$ so that we have

(6)
$$\{\begin{array}{c}x=a\frac{1}{k}e^{kb_0}\mathrm{sin}\left(ka+\sqrt{gk\mathrm{sin}\alpha }t\right),\hfill \\ y=b_0+\frac{1}{k}e^{kb_0}\mathrm{cos}\left(ka+\sqrt{gk\mathrm{sin}\alpha }t\right),\hfill \\ z=0,\hfill \end{array}$$

with $`a`$. The above formula represents the parametrization of a smooth trochoid (if $`b_0<0`$) or of a cycloid with upward cusps (if $`b_0=0`$); it also explains why we imposed the condition $`b_00`$ as otherwise we would obtain a self-intersecting curve.

The cusp formation by some edge waves was clearly demonstrated in wave-tank experiments and found by field measurements on the ocean beach (see Komar 1998).

Another aspect of interest is the amplitude of the edge wave. To determine the elevation with respect to the reference plane

$$\{z=\frac{\mathrm{sin}\alpha }{2k}e^{2kb_0}+(b_0y)\mathrm{tan}\alpha :yb_0\},$$

we compute the distance of a point $`(x,y,z)`$ lying on the free surface (4) to this plane,

$$d=z\mathrm{cos}\alpha +(yb_0)\mathrm{sin}\alpha \frac{\mathrm{sin}\alpha }{2k}e^{2kb_0},$$

with the understanding that positive/negative values on the right-hand side mean that the point lies above/below the plane. Since (with $`bb_0`$)

(7)
$$d=\frac{\mathrm{sin}\alpha }{2k}\left(2e^{kb(1+\mathrm{tan}\alpha )kb_0\mathrm{tan}\alpha }\mathrm{cos}(ka+\sqrt{gk\mathrm{sin}\alpha }t)e^{2kb(1+\mathrm{tan}\alpha )2kb_0\mathrm{tan}\alpha }\right),$$

we see that the amplitude of the edge wave decays exponentially away from the shoreline (as $`b\mathrm{}`$). The same conclusion is reached by a formal linear approximation (see Whitham 1979) and explains why edge waves are called “trapped waves”.

As expected - and ensured by (7), the amplitude of the edge wave varies with the parameters $`a,bb_0`$. From (7) we also infer that, at a fixed $`bb_0`$, the crests and troughs correspond to the maximal/minimal values of $`\mathrm{cos}(ka+\sqrt{gk\mathrm{sin}\alpha }t)`$. At a fixed time $`t0`$, we obtain the crest curves (with $`m`$ fixed and $`bb_0`$ playing the role of a running parameter)

$$\{\begin{array}{c}x=\frac{2m\pi }{k}\frac{1}{k}\sqrt{gk\mathrm{sin}\alpha }t,\hfill \\ y=b(1+\mathrm{tan}\alpha )b_0\mathrm{tan}\alpha +\frac{1}{k}e^{kb(1+\mathrm{tan}\alpha )kb_0\mathrm{tan}\alpha },\hfill \\ z=(b_0b)(1+\mathrm{tan}\alpha )+\frac{\mathrm{tan}\alpha }{2k}e^{2kb_0}(1e^{2k(bb_0)(1+\mathrm{tan}\alpha )}),\hfill \end{array}$$

and the trough curves

$$\{\begin{array}{c}x=\frac{(2m+1)\pi }{k}\frac{1}{k}\sqrt{gk\mathrm{sin}\alpha }t,\hfill \\ y=b(1+\mathrm{tan}\alpha )b_0\mathrm{tan}\alpha \frac{1}{k}e^{kb(1+\mathrm{tan}\alpha )kb_0\mathrm{tan}\alpha },\hfill \\ z=(b_0b)(1+\mathrm{tan}\alpha )\mathrm{tan}\alpha +\frac{\mathrm{tan}\alpha }{2k}e^{2kb_0}(1e^{2k(bb_0)(1+\mathrm{tan}\alpha )}).\hfill \end{array}$$

Note that for both the crest and trough curves the value of $`x`$ (at a given time) is fixed: standing at that location and looking towards the sea, these curves, orthogonal to the shoreline, are fully visible at certain instants. Indeed, taking into account the fact that on a crest/trough curve the deviation from the reference plane is

$${}_{}{}^{+}\frac{\mathrm{sin}\alpha }{2k}\left(2e^{kb(1+\mathrm{tan}\alpha )kb_0\mathrm{tan}\alpha }e^{2kb(1+\mathrm{tan}\alpha )2kb_0\mathrm{tan}\alpha }\right),bb_0,$$

in view of (7), the monotonicity of the right-hand side shows that the deviation becomes smaller (in absolute value) with the distance from the shore. This feature can be easily recognized in the graphical representation of the edge wave given below.

The edge wave viewed from offshore. The sinusoidal longshore structure and the

exponential offshore decay in amplitude are clearly visible.

REFERENCES

Boss, E., Gonzalez, F., Sakate, K. & Mofjeld, H. 1995 Edge wave and non-

trapped modes of the 25 April 1992 Cape Mendocino tsunami. Pure Appl.

Geophys. 144, 409-426.

Constantin, A. 2001 On the deep water wave motion. J. Phys. A 34, 1405-1417.

Crapper, G. 1984 Introduction to Water Waves. Ellis Horwood, Chichester.

Ehrenmark, U. 1998 Oblique wave incidence on a plane beach: the classical problem

revisited. J. Fluid Mech. 368, 291-319.

Evans, D. 1988 Mechanisms for the generation of edge waves over a sloping beach.

J. Fluid Mech. 186, 379-391.

Gerstner, F. 1809 Theorie der Wellen samt einer daraus abgeleiteten Theorie der

Deichprofile. Ann. Physik 2, 412-445.

Johnson, R. 1997 A Modern Introduction to the Mathematical Theory of Water

Waves. Cambridge University Press.

Komar, P. 1998 Beach Processes and Sedimentation. Prentice-Hall, Inc.

Lamb, H. 1932 Hydrodynamics. Cambridge University Press.

Lighthill, J. 1978 Waves in Fluids. Cambridge University Press.

Minzoni, A. & Whitham, G. B. 1977 On the excitation of edge waves on beaches.

J. Fluid Mech. 79, 273-287.

Mollo-Christensen, E. 1982 Allowable discontinuities in a Gerstner wave field.

Phys. Fluids 25, 586-587.

Whitham, G. B. 1976 Nonlinear effects in edge waves. J. Fluid Mech. 74, 353-368.

Whitham, G. B. 1979 Lectures on Wave Propagation. Springer Verlag, Berlin.

Yeh, H. 1985 Nonlinear progressive edge waves: their instability and evolution. J.

Fluid Mech. 152, 479-499.

Yih, C. 1966 Note on edge waves in a stratified fluid. J. Fluid Mech. 24, 765-767.


text:
# New constraints on ultrashort-ranged Yukawa interactions from atomic force microscopy*footnote **footnote *To appear in Physical Review D.

## I Introduction

Motivated in part by the hierarchy problem, a number of theoretical models have been developed recently which raise the possibility that the energy scale at which gravity is unified with the other fundamental forces could be as low as $`1`$ TeV . In these models, large extra spatial dimensions (which for compact dimensions could be as large as $`1`$ mm) arise which might be accessible experimentally . While the effects of these new dimensions are generally model dependent, constraints on extra dimensional physics which are relatively independent of theory can be obtained from experiments searching for deviations from Newton’s law of gravity at sub-millimeter distances. Such experiments can be classified according to the dominant background force acting between the test bodies used. For separations $`10^4`$m gravity produces the dominant background and a number of groups are conducting experiments testing Newtonian gravity at $`0.1`$ mm scales. For smaller separations, Casimir forces provide the dominant background force, and the experiments setting the best limits in this regime are those testing the Casimir force law. In this paper, we will briefly review the theoretical motivation and phenomenology used in sub-millimeter force experiments. We then present new limits extracted from a recent experiment by the Riverside group , which uses an atomic force microscope to test the Casimir force. After comparing these results with existing limits, we conclude by describing future possibilities for using AFM experiments to set limits on new forces and extra dimensional physics.

## II Phenomenology

Since gravity is a theory of spacetime, it is inevitable that the effects of extra dimensions will modify Newton’s law of gravity which is only valid in a 4-dimensional spacetime. Although Newtonian gravity has been well-tested over long distance scales, our understanding of gravity at sub-millimeter scales is quite limited . Phenomenologically, the deviations from Newton’s law arise naturally in two different ways . In models with large compact dimensions such as originally proposed by Arkani-Hamed, Dimopoulos, and Dvali (ADD) , one finds that the potential energy between two point particles with masses $`m_1`$ and $`m_2`$ separated by a distance $`rR_{}`$ (the characteristic size of the compact dimensions) is given by the usual Newtonian gravitational potential energy with a Yukawa correction :

$$V(r)=\frac{Gm_1m_2}{r}\left(1+\alpha e^{r/\lambda }\right);(rR_{}),$$
(1)
where the dimensionless constant $`\alpha `$ depends on the nature of the extra dimensions, and where $`\lambda R_{}`$. For example, for a toroidal compactification with all $`n`$ extra dimensions having equal size, $`\alpha =2n`$ . When the separation between the masses decreases to the point where $`rR_{}`$, the usual inverse square law of gravity changes to a new power-law:

$$V(r)=\frac{G_{4+n}m_1m_2}{r^{n+1}}(rR_{}),$$
(2)
where $`G_{4+n}GR_{}^n`$ is the fundamental gravitational constant in the full $`4+n`$ dimensional spacetime. The size $`R_{}`$ of the extra dimensions, according to the ADD models, is then related to the energy scale $`M_{}(1/G_{4+n})^{1/(2+n)}`$ at which the unification occurs by

$$R_{}\frac{1}{M_{}}\left(\frac{M_{Pl}}{M_{}}\right)^{2/n}10^{32/n17}\text{cm},$$
(3)
where the Planck mass is $`M_{Pl}=1/G^{1/2}10^{19}`$ GeV, and we have assumed $`M_{}1`$ TeV . For $`n=1`$, $`R_{}10^{15}`$ cm is obviously excluded by solar system tests of Newtonian gravity . If $`n=2`$, $`R_{}1`$ mm, which is the scale currently being probed by a number of sub-millimeter gravity experiments . For $`n=3`$, Eq. (3) gives $`R_{}5`$ nm which is roughly the range accessible to experiments using atomic force microscopy like those discussed below. It is important to recognize, however, that while $`M_{}1`$ TeV is a natural scale, it is little more than a reasoned guess. Therefore, the scales suggested by Eq. (3) should be taken simply as a heuristic guide for designing experiments. It should also be noted that the Yukawa correction given in Eq. (1) arises naturally in ways unrelated to extra dimensional physics. For example, Yukawa potentials arise from new forces generated by the exchange of bosons of mass $`\mu =1/\lambda `$ (e.g., scalar, graviphoton, dilaton) .

A second class of extra dimensional models, such as that originally proposed by Randall and Sundrum , allow the extra dimensions to be non-compact, but warped. In these models, the leading order correction to the Newtonian potential energy between two point masses takes on a power-law form :

$$V(r)=\frac{Gm_1m_2}{r}\left(1+\frac{2}{3k^2r^2}\right),(r1/k),$$
(4)
where $`1/k`$ is the warping scale. The introduction of a new length scale $`r_0`$ such as $`1/k`$ or $`R_{}`$ into a theory allows one to generalize Eq. (4) to arbitrary powers $`N`$:

$$V(r)=\frac{Gm_1m_2}{r}\left[1+\alpha _N\left(\frac{r_0}{r}\right)^{N1}\right],$$
(5)
where $`\alpha _N`$ is a dimensionless constant. Power-law corrections to Newtonian gravity also arise from new interactions involving the exchange of 2 massless quanta. For example, $`N=2`$ may arise from the simultaneous exchange of two photons or two massless scalars , and $`N=3`$ characterizes the exchange of two massless pseudoscalars . Potentials with $`N=5`$ arise from massless axion exchange , and also from the exchange of a (massless) neutrino-antineutrino ($`\nu \overline{\nu }`$) pair .

All of the extra dimensional models discussed above are speculative at this stage, and much work remains to be done to test their consistency and viability. However, they are important in demonstrating that if extra spatial dimensions exist, new fundamental energy scales must also exist, and that these scales will modify Newtonian gravity at some level. Furthermore, from the models that have been developed we see that these modifications generally take the form of Yukawa or power-law potentials. Therefore, while there remain many unanswered questions regarding models of extra-dimensional physics, the phenomenology supporting short distance gravity experiments rests upon a much firmer footing.

As noted earlier, experiments testing Newtonian gravity over separations $`10^3`$m have set significant limits on the Yukawa coupling $`\alpha `$ for ranges $`\lambda 10^3`$ m . In addition, precise gravity experiments have set the best limits on power-law interactions for $`N<5`$ (see Table I). Motivated in part by the new extra dimensional models discussed above, several experimental groups  are now probing Newtonian gravity at sub-millimeter scales and have begun to set significant constraints on new physics at ranges $`0.1`$ mm. It is expected that within the next few years such experiments will extend our understanding of gravity down to separations where Casimir forces form the significant background.

## III New Forces and Tests of the Casimir Force

Coincidentally, the past few years have also witnessed the first precise tests of the Casimir force . Lamoreaux  used a torsion balance to measure the Casimir force between a disk and spherical lens, while Ederth  measured the force between two crossed cylinders. The practical implications of the Casimir force have been considered recently with a MEMS (microelectromechanical) device by Chan, et al. at Bell Labs . Finally, a series of force measurements specifically devised to test various aspects of the Casimir force has been conducted by the Riverside group using an atomic force microscope . This entire experimental effort has motivated theorists to more carefully examine the effects of surface roughness, temperature, and finite conductivity which produce significant deviations from the idealized case of the Casimir force between two perfectly conducting, smooth bodies . The inclusion of these corrections has been shown to be crucial in obtaining agreement between theory and experiment .

While the goal of all these experiments is to improve our understanding of the origin and behavior of Casimir forces, they have also been used to set limits on new interactions . These limits are in addition to those obtained in Refs.  from earlier, much less precise Casimir and van der Waals force measurements . It can be shown that Casimir force measurements are relatively ineffective compared to longer ranged gravity experiments in setting limits on power-law interactions of the form given by Eq. (5). The reason is that for $`N<5`$, such power-law forces between macroscopic bodies are relatively insensitive to the separation of the bodies when these are placed in close proximity . Thus, there is no advantage to performing very short distance experiments, which are in any case generally less precise than larger scale gravity experiments, to search for power-law corrections to Newtonian gravity.

On the other hand, experiments searching for Yukawa interactions of the form given by Eq. (1) are only sensitive when the separation $`r`$ of the test bodies is of order $`\lambda `$ . Hence, to constrain the Yukawa coupling $`\alpha `$ for a given value of $`\lambda `$, one must devise experiments such that $`r`$ is of order $`\lambda `$. Therefore, if one wishes to constrain extra dimensions of size $`R_{}10^5`$ m using force measurements, one must inevitably confront large Casimir forces which grow rapidly for $`r10^5`$ m . This is why Casimir force experiments yield the best constraints on $`\alpha `$ for these separations.

## IV Limits from Recent Mohideen AFM Experiment

We turn next to the focus of this paper, which is the recent measurement of the Casimir force performed by the Riverside group  who used an atomic force microscope. (We will henceforth refer to the work of the Riverside group as “Mohideen, et al.”) Like the previous AFM measurements carried out by this group , this experiment measured the force between a metallized polystyrene sphere mounted on an AFM tip and a substrate composed of a flat sapphire disk. However, this experiment differed from the previous efforts in several important respects: (1) The sphere and disk were each coated with a layer of gold of thickness $`\mathrm{\Delta }=86.6`$ nm. Previously, aluminum coatings had been used instead, but this also required an additional surface layer of Au/Pd on top of the Al to prevent oxidation. The use of only a single gold layer in the new experiment significantly simplifies the calculations of the forces involved, since one can treat the gold layer as infinitely thick as far as the Casimir force is concerned. Also, because Au is much more dense than the previously used layers, stronger limits on new forces can be obtained. (2) The gold layers were significantly smoother than previous coatings. The root-mean-square amplitude of the gold surfaces was measured to be only $`1.0\pm 0.1`$ nm which means that corrections arising from the surface roughness could be neglected in all force calculations. (3) Electrostatic forces which plague nearly all short distance force measurements were reduced to $`1`$% of the Casimir force at the shortest separation. This meant that no subtractions were required to separate Casimir forces from spurious electrostatic background forces. The electrostatic force was used to arrive at an independent measurement of the surface separation including separation on contact of the two surfaces. (4) The measurements were performed over smaller separations $`a`$, $`62\text{nm}a350`$ nm which means this experiment can search for Yukawa forces with smaller ranges $`\lambda `$. However, the absolute error of the force measurements was somewhat larger than in the previous experiments, $`\mathrm{\Delta }F=3.5\times 10^{12}`$N. This is due to the thinner gold coating used in  which led to poor thermal conductivity of the cantilever. At the smallest measured separations, this error was still less than 1% of the measured Casimir force. Thus, this experiment can be used to significantly tighten constraints on new Yukawa interactions as we will now show.

We begin by calculating the force arising between the AFM tip and the substrate using Eq. (1). First, it is easy to show that the Newtonian contribution is negligible. Since the diameter of the sphere (including the gold layer) was $`2R=191.3`$ $`\mu `$m, while the disk diameter was $`L=1`$ cm, it follows that $`RL`$ which allows us to consider each atom of the sphere as if it were placed directly above the center of the disk. In this case, the vertical component of the Newtonian gravitational force acting between the disk and a sphere atom of mass $`m_1`$ located at a distance $`lL`$ above the disk is

$`f_{N,z}(l)`$ $`=`$ $`{\displaystyle \frac{}{l}}\left[Gm_1\rho _{\mathrm{disk}}2\pi {\displaystyle \underset{0}{\overset{L}{}}}r𝑑r{\displaystyle \underset{l}{\overset{l+D}{}}}{\displaystyle \frac{dz}{\sqrt{r^2+z^2}}}\right]`$ (6)
$``$ $`2\pi Gm_1\rho _{\mathrm{disk}}D\left[1{\displaystyle \frac{D+2l}{2L}}\right],`$ (7)
where $`\rho _{\mathrm{disk}}`$ is the disk density, $`D=1`$ mm is the thickness of the disk, and only the first order terms in $`D/L`$ and $`l/L`$ have been retained. The Newtonian gravitational force acting between the disk and the sphere is obtained from Eq. (7) by integration over the sphere

$$F_{N,z}\frac{8}{3}\pi ^2G\rho _{\mathrm{disk}}\rho _{\mathrm{sphere}}DR^3\left(1\frac{D}{2L}\frac{R}{L}\right),$$
(8)
where $`\rho _{\mathrm{sphere}}`$ is the density of the sphere. Even if the sphere and disk were composed of solid vacuo-distilled gold with $`\rho _{\mathrm{disk}}=\rho _{\mathrm{sphere}}=18.88\times 10^3`$kg/m3, one finds from Eq. (8) that $`F_{N,z}6\times 10^{16}\text{N}\mathrm{\Delta }F`$, so the usual Newtonian gravitational force can be neglected.

The force arising from the Yukawa correction to Newtonian gravity given by the second term of Eq. (1) should be calculated taking into account the actual compositions of the test bodies. For a polystyrene sphere $`\rho _{\mathrm{sphere}}=1.06\times 10^3`$kg/m3, for a sapphire disk $`\rho _{\mathrm{disk}}=4.0\times 10^3`$kg/m3, and for the vacuo-distilled gold covering layers $`\rho _{\mathrm{Au}}=18.88\times 10^3`$kg/m3. The Yukawa force $`F_Y`$ arising from Eq. (1) can be easily obtained using the same procedure used to calculate the Newtonian gravitational force. The result is

$`F_Y(a)=4\pi ^2G\alpha \lambda ^3e^{a/\lambda }R`$ (9)
$`\times \left[\rho _{\mathrm{Au}}(\rho _{\mathrm{Au}}\rho _{\mathrm{disk}})e^{\mathrm{\Delta }/\lambda }\right]\left[\rho _{\mathrm{Au}}(\rho _{\mathrm{Au}}\rho _{\mathrm{sphere}})e^{\mathrm{\Delta }/\lambda }\right],`$ (10)
where $`a`$ is the sphere-disk separation distance, and $`\mathrm{\Delta }=86.6`$ nm is the thickness of the gold coatings.

The most stringent limits on new forces are obtained for the smallest possible separation. According to Ref. , the rms deviation of the measured force ($`F_{\mathrm{experiment}}`$) from the theoretical value ($`F_{\mathrm{theory}}`$) of the sum of all the known forces (Casimir and electrostatic) is given by

$$\sigma _F=\sqrt{\frac{\left(F_{\mathrm{experiment}}F_{\mathrm{theory}}\right)^2}{N_{\mathrm{trials}}}}=\text{3.8 pN},$$
(11)
where $`N_{\mathrm{trials}}=2583`$ is the total number of force measurements. This is slightly larger than 3.5 pN, the experimental uncertainty of the force measurement at $`a=62`$ nm. Therefore, to set constraints on the Yukawa coupling constant $`\alpha `$, we will assume $`|F_Y(a)|\sigma _F=`$ 3.8 pN. The resulting limit is represented by the bold curve labeled “Mohideen 2000” in Fig. 1.

Before comparing the constraint obtained here with limits obtained from other experiments, let us first discuss various sources of error which could affect our result. First, the uncertainty in the separation $`a`$ arising from the surface roughness is approximately $`\delta a1`$ nm which gives $`\delta a/a2\%`$. This will produce a significant effect when $`\lambda \delta a1`$ nm, but as we will see below, other experiments produce better limits in this range. For $`\lambda a`$, where the most stringent limits are obtained, the resulting errors are only $`2\%`$ which can be neglected. Second, it is important that all known forces be included in $`F_{\mathrm{theory}}`$ in using Eq. (11). In short distance force experiments, electrostatic effects can be significant. However, in the Mohideen experiment, the measured residual potential difference of 3 mV between the gold sphere and disk leads to a force only 0.1% of the Casimir force at $`a=62`$ nm. Although this effect is small, it was included in $`F_{\mathrm{theory}}`$.

Finally, the most significant difficulties arising in extracting limits from Casimir force experiments involve calculating the Casimir force between the interacting bodies. Due to practical difficulties, none of the most recent tests of the Casimir force use parallel plates for the interacting bodies, for which accurate calculations can be made. (The experiment by Carugno  uses parallel plates, but has not yet reached sufficient sensitivity to detect the Casimir force.) Rather, curved surfaces are used which requires the use of the proximity force theorem (PFT) to calculate the Casimir force . Since the PFT produces reasonable results when applied to Casimir forces, it is probably correct. Although there has not been a rigorous proof of the PFT in this context  (see Ref.  for an interesting semi-classical argument), we will assume that its use in the analysis of the Mohideen experiment is valid.

Other difficulties in Casimir force calculations arise from the use of tabulated values of the complex dielectric constant as functions of frequency for the materials comprising the test bodies. It has been pointed out that if different interpolation schemes are applied to these discrete data, the resulting calculated Casimir force can differ by as much as 4% . In addition, the tabulated dielectric properties were obtained using bulk samples and not the coatings which are actually used in the experiment. It would then be best to directly measure the dielectric properties of the actual test bodies to accurately determine $`F_{\mathrm{Casimir}}`$ . Since this cannot be done in the present context, and since there is no obvious discrepancy or inconsistency arising in the Mohideen experiment, we will not consider this problem here and leave it to future experiments to address this issue.

## V Comparison with Other Limits

As can be seen from Fig. 1, the Casimir force measurement between the gold surfaces carried out by Mohideen, et al. strengthens the previously known constraints obtained from their earlier experiment  using aluminum surfaces (labeled by “Mohideen 1999”) by as much as a factor of 19 within the range $`4.3\times 10^9\text{m}\lambda 1.5\times 10^7`$m, with the most significant improvement occurring at $`\lambda =`$(5–10) nm. These constraints are up to 4500 times more stringent than those obtained from older Casimir and van der Waals force measurements between dielectrics (curves labeled by “Casimir” and “vdW” respectively). In the same figure, the curve labeled by “Lamoreaux” exhibits the constraint obtained in Ref.  from the Casimir force measurement which used a torsion pendulum . The curve denoted by “Ederth” gives the new constraints obtained in Ref.  from the recent experiment which measured the Casimir force between two crossed cylinders . We see that curve “Mohideen 2000,” which exhibits the limits obtained in this paper, represents the most stringent bound in the interaction range $`1.1\times 10^8\text{m}\lambda 1.5\times 10^7`$m.

It is evident from Fig. 1 that for separations $`10^5`$ m much work is needed to improve the sensitivity of force experiments to Yukawa interactions of the strength predicted by extra dimensional models ($`\alpha 1`$–10). As was shown in Ref. , the constraints following from the experiment in Ref.  can be improved by up to four orders of magnitude in the range $`\lambda 10^4`$ m which is what is required to constrain physically interesting values of $`\alpha `$. Still, experiments using atomic force microscopy  remain almost fifteen orders of magnitude below the sensitivity needed to reach the value $`\alpha `$1–10 in the interaction range $`\lambda 10^7`$m. Furthermore, significant (and sometimes even more stringent) constraints on the extra dimensional mass scale $`M_{}`$ have been obtained from accelerator experiments , astrophysics , and cosmology . However, one must remember that gravity and Casimir force experiments are also sensitive to other non-extra dimensional effects such as dilaton and moduli exchange which can lead to Yukawa forces with $`\alpha 10`$ . Therefore, all experiments which can reduce the allowed region in the $`\alpha `$$`\lambda `$ plane are meaningful.

## VI Future AFM Experiments

We believe that in the future the best ultrashort range limits will come from a new class of experiments designed specifically to search for new forces rather than to test the Casimir force. In such experiments the Casimir effect is an unwanted background that needs to be suppressed. One suppression technique proposed recently  relies on the observation that the Casimir force depends on the electronic properties of interacting samples, whereas any gravitational interactions (either conventional or new), and virtually all other proposed interactions, involve couplings to nuclei (as well as to electrons). Hence by comparing the interactions of two test masses which have very similar electronic properties, but different nuclei, with a common attracting source mass, we may be able to subtract out the common Casimir background. This would leave a residual nuclear-nuclear interaction between each of the samples and the source, and these should be different for the two samples. We refer to this technique as the iso-electronic effect, and it can be implemented in two ways. One is to choose as the test masses elements such as Cu and Au whose nuclei are quite different, but which are known to have very similar Casimir properties . Alternatively, the test samples could be chosen to be isotopes of the same element whose nuclei would contain different numbers of neutrons, but whose electronic properties would be extremely similar. In either case, any residual differences in the electronic properties of the test masses would still have to be calculated, but this should in principle present no significant problems .

A possible design for an AFM experiment  which utilizes the iso-electronic effect, is shown schematically in Fig. 2. In this experiment, a gold-coated sphere would be attached to an AFM tip and oscillated with frequency $`\omega _{\mathrm{tip}}`$. The substrate would be composed of alternating strips of different isotopes of the same element, and would be oscillated horizontally with frequency $`\omega _{\mathrm{sub}}`$. The signal for a new force would then be a force on the tip that depends on $`\omega _{\mathrm{sub}}`$ which should not be present if surface roughness and electrostatic effects are the same for the strips. Since the Yukawa force of the AFM tip on the substrate couples only to the mass within a range $`\lambda `$ of the surfaces, the net sphere-substrate force will be proportional to the difference in the densities of the strips:

$$F_{\mathrm{Yukawa}}(\rho _{\mathrm{sub}}\rho _{\mathrm{sub}}^{}).$$
(12)
Therefore, it is important to select elements having isotopes with the largest possible mass density differences. Possible candidates include ruthenium, osmium, nickel, and palladium, all of which can have isotope mass density differences exceeding $`900`$ kg/m3. Other possible designs for ultrashort distance experiments can be found in Refs. .

## VII Conclusions

It is clear that experiments which search for new forces in the Casimir regime will be faced with significant challenges. However, there is strong theoretical motivation to conduct such searches, and new techniques (including, possibly, those suggested here) will certainly be developed to overcome these obstacles. The AFM is a natural tool to be used in this effort as demonstrated by the experiments of Mohideen, et al. . Furthermore, this instrument can be easily adapted to experiments whose intent is to specifically search for new forces. Such laboratory force experiments have provided, and will continue to provide, important, and relatively model-independent constraints on new physics which complement those obtained from high-energy experiments and astrophysical observations.

###### Acknowledgements.

The authors thank A. Lambrecht, S. Howell, G. L. Klimchitskaya, U. Mohideen, R. Reifenberger, and S. Reynaud for helpful discussions. We also thank S. Karunatillake and M. West for their calculations of the isotopic mass densities. V. M. M. is grateful to the Centro Brasileiro de Pesquisas Físicas for their kind hospitality. This work was supported in part by the U. S. Department of Energy under contract No. DE-AC02-76ER071428 (E. F.), the Wabash College Byron K. Trippet Research Fund (D. E. K.), FAPERJ and CNPq (V. M. M. and M. N.).


text:
# Localization properties of two interacting particles in a quasi-periodic potential with a metal-insulator transition

## 1 Introduction

The metal-insulator transition (MIT) in disordered electronic systems has been the subject of intense research activities over the last two decades and still continues to attract much attention. For free electrons in disordered systems KraM93  the scaling hypothesis of localization AbrALR79  can successfully predict many of the universal features of the MIT. However, the influence of many-particle interactions on the MIT is not equally well understood BelK94  and recent investigations of an apparent MIT in two-dimensional (2D) systems even question the main assumptions of the scaling hypothesis KraDS96 ; BelK98 ; SimKSP97 ; SimKSp98 ; KraSSK97 ; PudBPB97 . A simple theoretical approach to the interplay of interactions and disorder is based on the two-interacting-particles (TIP) problem in 1D random She94 ; She96a ; RomS97a  or quasiperiodic potentials She96b ; EilGRS99 . Furthermore, numerical results for spinless fermions at finite particle density have given additional insight SchSSS98 ; SchJWP98 ; WSchJP01 . In general, these investigations have shown that changes in the wave function interferences due to many-particle interactions RomSV99 ; RomSV00  can lead to a rather large enhancement of the localization lengths in 1D and 2D SchJWP98 ; LeaRS99 ; RomLS99b .

The standard approach for computing localization lengths in disordered, non-interacting systems is the transfer-matrix method MacK83 . It has been used for investigations of the enhancement of the TIP localization length in a 1D random potential RomS97a ; FraMPW95  where there is no MIT as all wave functions are always localized. Other numerical approaches to the TIP problem have been based on the time evolution of wave packets She94 ; BriGKT98 , exact diagonalization WeiMPF95  or Green function approaches LeaRS99 ; OppWM96 ; SonK97 .

In the single-particle case, the 1D quasiperiodic Aubry-André model is known rigorously to exhibit an MIT for all states in the spectrum as a function of the quasiperiodic potential strength $`\mu `$ AubA80 . The ground state wave function is extended for $`\mu <1`$ and localized for $`\mu >1`$. The system at $`\mu _\mathrm{c}=1`$ is critical: there the wave functions decrease algebraically, not exponentially as in the localized case. Recently, we examined this model for TIP by means of the transfer-matrix method together with a careful finite-size-scaling analysis EilGRS99  following earlier analytical work of Refs. BarBJS96 ; BarBJS97 . We showed that the model for TIP exhibits an MIT as a function of $`\mu `$ at $`\mu _\mathrm{c}=1`$ as in the single-particle case. Our finite-size-scaling results for onsite (Hubbard) interaction suggest that the critical behavior, i.e., the value for the critical exponent $`\nu `$ of the correlation length, is also not affected by the interaction EilGRS99 . However, it has been demonstrated RomS97a ; LeaRS99  that a transfer-matrix-method approach applied to the TIP problem without finite-size scaling leads to unreliable localization lengths, i.e., it systematically overestimates the TIP localization length $`\lambda _2`$ in finite-sized samples in the case of vanishing interaction ($`U=0`$). In addition, simple extrapolations to infinite sample size RomS97a ; FraMPW95  may lead to an underestimation of $`\lambda _2`$ SonO98 . An alternative approach, which does not suffer from the above problem, is based on the decimation method and has also been applied recently to TIP in a 1D random potential LeaRS99 . This encouraged us to reexamine the localization lengths for TIP in 1D quasiperiodic potentials with Hubbard interaction with the decimation method. As we shall show in the following, we find that the general conclusions of Ref. EilGRS99  remain valid, i.e., the MIT is not affected by the interaction. The critical properties of the single-particle transition at $`\mu _\mathrm{c}=1`$ are not altered within the accuracy of our calculation. One-parameter scaling is obeyed for onsite interaction strengths up to $`U=10`$.

As an independent extension of these low-density results, Chaves and Satija ChaS97  have studied a model of nearest-neighbor interacting spinless fermions YanY66a  at finite particle density in the same quasiperiodic potential by means of Lanczos diagonalization for small systems up to chain size $`M=13`$. They have found evidence for a critical region. In order to reach much larger system sizes for interacting systems, one can employ the numerical density-matrix renormalization group (DMRG) Whi93 . With the DMRG the ground state properties in 1D can be obtained very accurately PaiPR97 ; Sch99 . In a recent paper SchRS01 , we studied the quasiperiodic model of Ref. ChaS97  at various densities and interaction strengths $`V`$ by DMRG. We compare the results with the present TIP data at the end of the paper.

The paper is organized as follows. In section 2 we describe the Hamiltonian of our TIP system and explain how to obtain the TIP localization lengths via decimation method. In section 3, we comment on the particular finite-size-scaling method employed and present the estimated critical parameters. We summarize and conclude in section 4.

## 2 The TIP system and the numerical approach

The Hamiltonian for TIP in the 1D quasiperiodic potential of the Aubry-André model is given as

$`𝐇`$ $`=`$ $`{\displaystyle \underset{n,m}{}}|n,mn+1,m|+|n,mn,m+1|+\mathrm{h}.\mathrm{c}.`$ (1)
$`+|n,m\left[\mu _n+\mu _m+U(n,m)\right]n,m|.`$
Here $`\mu _m2\mu \mathrm{cos}(\alpha m+\beta )`$ is the quasiperiodic potential of strength $`\mu `$ with $`\alpha /2\pi `$ being an irrational number. $`\beta `$ is an arbitrary phase shift and we choose $`\alpha /2\pi =(\sqrt{5}1)/2`$, i.e., the inverse of the golden mean. This value of $`\alpha /2\pi `$ may be approximated by the ratio of successive Fibonacci numbers — $`F_n=F_{n2}+F_{n1}=0`$, $`1`$, $`2`$, $`3`$, $`5`$, $`8`$, $`13`$, $`\mathrm{}`$ — as is customary in the context of quasiperiodic systems Gri99 . The Hubbard onsite interaction matrix $`U(n,m)`$ is diagonal, i.e., $`U(n,m)=U\delta _{nm}`$. The indices $`n`$ and $`m`$ correspond to the positions of each particle on a chain of length $`M`$. Now we use the decimation method LeaRS99 ; LamW80  to construct an effective Hamiltonian for the diagonal of the $`M\times M`$ lattice along which the cigar-shaped TIP wave function has its largest extent WeiMPF95 ; SonO98 . The quantity of interest is the TIP localization length $`\lambda _2`$ defined by the TIP Green function $`𝐆_\mathrm{𝟐}(E)`$ OppWM96 :

$$\frac{1}{\lambda _2}=\frac{1}{|M1|}\mathrm{ln}|1,1|𝐆_\mathrm{𝟐}|M,M|.$$
(2)
For TIP in 1D and 2D random potentials, this approach has led to high precision results supporting the proposed increase of the TIP localization lengths due to the repulsive interaction LeaRS99 ; RomLS99b . We remark that similar data have also been obtained for nearest-neighbor OppWM96  and long-ranged interactions BriGKT98 .

The correlation length $`\xi _{\mathrm{}}`$ for the infinite system may be obtained from the localization lengths $`\lambda (M)`$ for finite system sizes by the using one-parameter scaling hypothesis $`\mathrm{\Lambda }_M=f(M/\xi _{\mathrm{}})`$ Tho74  for the reduced localization lengths $`\mathrm{\Lambda }_M=\lambda (M)/M`$. The MIT is characterized by a divergent correlation length $`\xi _{\mathrm{}}(\mu )|\mu \mu _\mathrm{c}|^\nu `$ KraM93 . In order to reliably extract the critical parameters from the calculated values of $`\lambda _2(M)`$ one may apply a finite-size-scaling procedure MacK83  that numerically minimizes deviations of the data from the common scaling curve $`f`$. The critical exponent $`\nu `$ can then be extracted by fitting the $`\xi _{\mathrm{}}`$ obtained from finite-size scaling cm0106005 ; cm0106006 . This method was used previously EilGRS99  for finding the critical parameters of the present model.

Higher accuracy can be achieved by a method applied recently cm0106005 ; cm0106006 ; SleO99a ; MilRSU00  to the MIT in the Anderson model of localization. We construct a family of fit functions which include corrections to scaling such as (i) nonlinearities of the dependence of the scaling variable on the quasiperiodic potential strength and (ii) an irrelevant scaling variable which accounts for a shift of the crossing point of the $`\mathrm{\Lambda }_M(\mu )`$ curves as a function of $`\mu `$, i.e.,

$$\mathrm{\Lambda }_M=\stackrel{~}{f}(\chi _\mathrm{r}M^{1/\nu },\chi _\mathrm{i}M^y).$$
(3)
where $`\chi _\mathrm{r}`$ and $`\chi _\mathrm{i}`$ are the relevant and irrelevant scaling variables, respectively. $`\stackrel{~}{f}`$ is then Taylor expanded up to order $`n_\mathrm{i}`$ in terms of the second argument

$$\mathrm{\Lambda }_M=\underset{n=0}{\overset{n_\mathrm{i}}{}}\chi _\mathrm{i}^nM^{ny}\stackrel{~}{f}_n(\chi _\mathrm{r}M^{1/\nu }),$$
(4)
and each $`\stackrel{~}{f}_n`$ is Taylor expanded up to order $`n_\mathrm{r}`$:

$$\stackrel{~}{f}_n=\underset{i=0}{\overset{n_\mathrm{r}}{}}a_{ni}\chi _\mathrm{r}^iM^{i/\nu }.$$
(5)
Nonlinearities are taken into account by expanding $`\chi _\mathrm{r}`$ and $`\chi _\mathrm{i}`$ in terms of $`u=(\mu _\mathrm{c}\mu )/\mu _\mathrm{c}`$ up to order $`m_\mathrm{r}`$ and $`m_\mathrm{i}`$, respectively,

$$\chi _\mathrm{r}(u)=\underset{n=1}{\overset{m_\mathrm{r}}{}}b_nu^n,\chi _\mathrm{i}(u)=\underset{n=0}{\overset{m_\mathrm{i}}{}}c_nu^n,$$
(6)
with $`b_1=c_0=1`$. The fit function is being adjusted to the data by choosing the orders $`n_\mathrm{i},n_\mathrm{r},m_\mathrm{r},m_\mathrm{i}`$ up to which the expansions are carried out. Of course, the orders have to be taken not too large to keep the number of fit parameters $`a_{ni}`$, $`b_n`$, and $`c_n`$ reasonably small.

## 3 Numerical results for TIP

We calculate $`\lambda _2`$ at energy $`E=0`$ for $`20`$ values of the Hubbard interaction, i.e., $`U=0`$ (the non-interacting single-particle case), $`0.1`$, $`\mathrm{}`$, $`0.9`$, $`1`$, $`2`$, $`\mathrm{}`$, $`10`$ for $`6`$ system sizes $`M=13`$, $`21`$, $`34`$, $`55`$, $`89`$, $`144`$. For $`U=0`$ and $`1`$, we also have data for $`M=233`$ and $`377`$. The quasiperiodic potential strengths $`\mu `$ is chosen close to the single-particle transition at $`\mu _\mathrm{c}1`$ and ranges typically from $`0.95`$ to $`1.05`$. As in Ref. EilGRS99  we average the results over different randomly chosen phase shifts $`\beta `$ in order to reduce the fluctuations. The number of $`\beta `$ values used in this averaging ranges from $`5000`$ for $`M=13`$ to $`1000`$ for $`M=377`$. In order to perform the non-linear fit necessary for the finite-size-scaling procedure as outlined in section 2, we use the Levenberg-Marquardt method MilRSU00 ; PreFTV92 . As the decimation-method data — like the transfer-matrix-method results EilGRS99  — are still rather noisy we have to suitably limit the ranges of the quasiperiodic potential strength $`\mu `$ and/or the system sizes $`M`$ used for fitting the data.

For $`U=0`$ and $`1`$, which were examined by the transfer-matrix method in detail EilGRS99 , the best fit is obtained for $`n_r=3`$, $`n_i=2`$, $`m_r=2`$ and $`m_i=1`$. For $`U=0`$ we use the data for $`\mu `$ ranging from $`0.96`$ to $`1.01`$ and $`M=55,89,144,233`$, and $`377`$; for $`U=1`$ we use all system sizes $`M=13`$, $`\mathrm{}`$, $`377`$ and $`0.97\mu 1.05`$. Figures 1 and 2 show the resulting TIP localization lengths for $`U=0`$ and $`1`$. Also shown are the fits of the finite-size-scaling curves to the data as given by Eq. (3) for $`U=0`$ and $`1`$, respectively. We find that for both $`U`$ values, there is an apparent transition close to $`\mu _\mathrm{c}=1`$. For the case $`U=0`$, we also observe a systematic shift of the crossing point with increasing system sizes necessitating the inclusion of an irrelevant scaling variable as discussed in section 2. The transition point is not so clearly distinguished for $`U=1`$, albeit the different behavior for $`\mu 1`$ and $`\mu 1`$, namely the increase and decrease, respectively, of $`\mathrm{\Lambda }_M`$ with increasing $`M`$, is clearly seen.

The corresponding plots of the scaling curves are displayed in the Figs. 3 and 4. The scaling curves are much better than reported previously EilGRS99  for the transfer-matrix-method data. The critical parameters can consequently be estimated to be $`\mu _\mathrm{c}=0.989\pm 0.001`$, $`\nu =1.00\pm 0.15`$ for $`U=0`$ and $`\mu _\mathrm{c}=0.997\pm 0.001`$, $`\nu =1.19\pm 0.16`$ for $`U=1`$. The irrelevant scaling exponents are close to $`y=1.8\pm 0.2`$ and $`y=0.15\pm 0.1`$ for $`U=0`$ and $`1`$, respectively. Note that the quoted errors correspond to the standard deviations estimated from the non-linear fit procedure. In this way the accuracy is significantly overestimated. Since it is apriori not clear, which values $`n_\mathrm{i},n_\mathrm{r},m_\mathrm{r},m_\mathrm{i}`$ to use, we estimate the true errors from a comparison of various fits with different $`n_\mathrm{i},n_\mathrm{r},m_\mathrm{r},m_\mathrm{i}`$. Even in the case of extremely high precision data close to the MIT in the Anderson model of localization, this has been shown MilRSU00  to increase the error by one order of magnitude. Therefore we conclude that the interaction strength $`U`$ for TIP does not influence the MIT in the quasiperiodic potential within the accuracy of the present calculation.

Further results for larger $`U`$ values are collected in Table 1. The expansion orders $`n_\mathrm{i},n_\mathrm{r},m_\mathrm{r},m_\mathrm{i}`$, the system sizes and ranges of the quasiperiodic potential strength have been chosen in order to minimize the $`\chi ^2`$ statistics and to get the most convincing scaling fit. Furthermore, one has to check that various initial parameters ($`a_n,b_n,c_n`$) converge to the same values of the critical quasiperiodic potential strength $`\mu _\mathrm{c}`$ and the critical exponent $`\nu `$. Figures 5 and 6 show the values obtained in this way. For almost all cases the critical quasiperiodic potential strength $`\mu _\mathrm{c}`$ remains close to $`1`$, the only exceptions are $`U=0`$ and $`0.1`$, when $`\mu _\mathrm{c}=0.99`$ and $`0.98`$, respectively. However, since we know that the transition in the single-particle case is exactly at $`\mu _\mathrm{c}=1`$ AubA80 , this observation can be used to estimate the true error of the estimate for $`\mu _\mathrm{c}`$. Thus comparing with the $`\mu _\mathrm{c}`$ estimates for $`U0`$, we find that the errors calculated within the non-linear fitting procedure are significantly underestimated as discussed above. We therefore conclude that within the accuracy of our calculation there is no change of the critical quasiperiodic potential strength $`\mu _\mathrm{c}`$ for the Hubbard interaction in the range $`0U10`$. The same argument leads to the conclusion that within the error bars the critical exponent $`\nu `$ does not change with the Hubbard interaction strength and is close to $`1`$. This is an agreement with the previous results obtained by the transfer-matrix method and finite-size scaling EilGRS99 . We stress that the critical exponents can only be obtained with much less accuracy than the transition point $`\mu _\mathrm{c}`$ as shown in Table 1.

## 4 Conclusions

In this work, we have studied the interplay of disorder and interactions for a quantum system at very low density (TIP). We calculated the pair localization lengths for a quasiperiodic potential and Hubbard interaction by means of the decimation method and extracted the critical parameters from the fit using the one-parameter scaling hypothesis. For both non-interacting particles as well as onsite interaction we obtain the value of the critical quasiperiodic potential strength $`\mu _\mathrm{c}=1`$ and the critical exponent $`\nu 1`$ in agreement with the previous results of transfer-matrix-method calculations and finite-size scaling EilGRS99 . The results for $`U>1`$ show that this conclusion remains valid also for much stronger interactions.

Let us briefly compare these results to the finite density situation. For $`N`$ interacting spinless fermions on a 1D ring of circumference $`M`$ with Aubry-André onsite potential $`\mu `$ and nearest-neighbour interaction $`V`$ it is possible to treat system lengths up to about $`M100200`$ using the DMRG. We applied SchRS01  the finite lattice algorithm for non-reflectionsymmetric models as described in Sch96 . For a system of free fermions at finite density like $`\rho =1/2`$ (incommensurate compared to the wave vector of the quasiperiodic potential — an irrational multiple of $`\pi `$), we reproduced SchRS01  the expected transition at $`\mu _\mathrm{c}=1`$ in agreement with Refs. EilGRS99 ; ChaS97 . For attractive and repulsive interactions at $`\rho =1/2`$ the numerial results are available for only two system sizes ($`M=34`$ and $`144`$), therefore conclusions about these regimes appear rather speculative. At commensurate densities $`\rho _ilim_n\mathrm{}F_{ni}/F_n0.618`$, $`0.382`$, $`0.236`$, and $`0.146`$ — corresponding to $`i=1`$, $`\mathrm{}`$, $`4`$ — and in the repulsive regime (nearest-neighbour interaction $`V>0`$), the ground state is localized for $`\mu >0`$ SchRS01  in agreement with previous studies for disordered and periodically disturbed systems SchSSS98 ; Sch99 . The above increase of the localization lengths as predicted by the arguments for TIP She94  is most likely too small SchJWP98  to be detected by the present accuracy. For attractive interactions $`V`$, all densities $`\rho _i`$ and $`\mu 0`$, the system shows a Peierls-like transition from insulating to metallic phase at $`V1.4`$ SchRS01  in agreement with the weak-coupling renormalization group treatment VidMG99  of spinless fermions on a Fibonacci lattice.

In conclusion, we have studied the influence of interactions on an MIT in a quasiperiodic model in 1D. Our results suggest that the delocalization found for low density TIP in the localized phase cannot simply be extrapolated to the finite-density situation. At finite densities, other effects such as a Peierls-like commensurability become important and dominate the transport properties.

###### Acknowledgements.

We thank M. Leadbeater for help with the decimation method and C. Schuster for stimulating discussions. We gratefully acknowledge the support of the SMWK and the Deutsche Forschungsgemeinschaft within Sonderforschungsbereich 393.


text:
# The peculiar X-ray properties of the low-redshift quasar MR 2251-178 and its environment: some surprises

## 1. Introduction

MR 2251-178 at redshift $`z`$=0.064 was the first quasar initially discovered by X-ray observations (Ricker et al. 1978), and the first one found to host an ionized absorber (Halpern 1984). The quasar turned out to be outstanding in many respects. It has a high ratio of $`L_\mathrm{x}/L_{\mathrm{opt}}`$, is surrounded by the largest quasar emission line nebula known, and is located in the outskirts of a cluster of galaxies (e.g., Bergeron et al. 1983). Our analysis of the ROSAT PSPC observations of this source and its environment revealed a number of surprises, summarized below.

## 2. X-ray results: some surprises

Results presented here are based on the analysis of two ROSAT PSPC observations of MR 2251-178, carried out in Nov. 1990 (all-sky survey) and Nov. 1993 (18.3 ksec duration). Details of the data analysis and the discussion of results is presented by Komossa (2001). The following results were obtained:

### X-ray variability.

There is evidence for an X-ray flaring event with variability by a factor $``$2 within 10 ksec during the ROSAT all-sky survey; rarely observed in such a luminous source, but similar to PDS 456 (Reeves et al. 2000) and PKS 0558-504 (Wang et al. 2001).

The source was a factor $``$3 brighter during the later pointed ROSAT observation. The mean observed soft X-ray luminosity at that epoch, $`L_{(0.12.4)\mathrm{keV}}10^{45}`$ erg/s, places the quasar among the most X-ray luminous AGN in the local universe.

### X-ray spectrum.

Remarkably, we do not detect any excess X-ray cold absorption expected to originate from the giant \[OIII\] gas nebula surrounding MR 2251-178. This finding excludes the presence of a huge HI extent of the HII emission line gas (along the line-of-sight), and constrains some formation scenarios (e.g., Shopbell et al. 2000) of the gas nebula.

On the other hand, we do find that the X-ray spectrum is modified by lots of amounts of gas along the line of sight - but this material is highly ionized. As indicated by some, but not all, earlier X-ray observations of MR2251-178, a single powerlaw does not provide an acceptable spectral fit ($`\chi _{\mathrm{red}}^2=5.4`$). The presence of an ionized absorber markedly improves the fit. We find an ionization parameter $`\mathrm{log}U=0.5`$ and a column density $`\mathrm{log}N_\mathrm{w}=22.6`$ of the ionized material (see Komossa 1999 for a general review on warm absorbers). The presence of highly ionized material along the line of sight is consistent with the detection of absorption lines in the UV spectrum of MR 2251-178 (Monier et al. 2001). Our best-fit warm-absorber model still leaves some residuals at the low-energy part of the spectrum, suggesting the presence of a very soft excess.

### Surrounding sources.

None of the other optically bright member galaxies of the cluster to which MR 2251-178 belongs, are detected in X-rays. However, east of the quasar there is a significant excess of X-ray sources (a factor of 4 compared to the log N – log S distribution of Hasinger et al. 1994), several of them without optical counterparts on the UK Schmidt plates.

### Intra-cluster gas emission:

The X-ray emission from the intra-cluster medium is weak or absent. We derive an upper limit on the X-ray luminosity of $`L_\mathrm{x}\mathrm{2\hspace{0.17em}10}^{42}`$ erg/s from the direction of the optical center of the cluster, weaker than other clusters of comparable richness. This may indicate that what appears as one single cluster could be just a chance projection of several poorer clusters or groups. Alternatively, the cluster X-ray emission might be off-set from the optical center.

## References

Bergeron J., et al. 1983, MNRAS 202, 125

Halpern J. 1984, ApJ 281, 90

Hasinger G., et al. 1994, A&A 275, 1

Komossa S. 1999, in ASCA-ROSAT workshop on AGN, ISAS Report, T. Takahashi & H. Inoue (eds), p. 149 \[also available at astro-ph/0001263\]

Komossa S. 2001, A&A 367, 801

Monier E.M., et al. 2001, ApJ, submitted

Reeves J.N., et al. 2000, MNRAS 312, L17

Ricker G.R., et al. 1978, Nature 271, 35

Shopbell P.L., Veilleux S., Bland-Hawthorn J. 1999, ApJ 524, L83

Wang T., et al. 2001, ApJ, in press