format: text

text:
# Abstract

## Abstract

In X-rays, the Galactic Centre (GC) Region appears as a complex of diffuse thermal and non-thermal X-ray emission intermixed with a population of luminous discrete sources. Here we present some new findings from the XMM-Newton GC survey, focusing particularly on Sgr A, Sgr A East, the Radio Arc region, and a newly discovered X-ray non-thermal filament.

## 1 Introduction

Our Galactic Centre Region provides, arguably, the best available laboratory in which to study the central activity of galaxies. Recent measurements of the proper motion of infrared stars around the dynamical centre of the Galaxy coupled with the enhanced radio brightness of the region, strongly support the presence of a super massive black hole of mass (2–3)$`\times 10^6`$ M (e.g., Genzel et al. 2000) in Sgr A. A number of unusual phenomena are also observed in the Galactic Centre Region: e.g., magnetic radio filaments, high density, fast-moving molecular clouds, massive star clusters, and high-temperature plasma extending over scales in excess of a few hundred parsecs, etc.

A wide-angle XMM-Newton survey of the Galactic Centre Region was carried out in the period 2000–2001 utilising 11 pointings in total, each with 10–25 ksec exposures (see Warwick 2002). This XMM-Newton Galactic Centre survey along with a complementary Chandra survey (Wang et al. 2002) has revealed the complex X-ray structure of the region. Here we report some of the recent discoveries which have stemmed from the XMM Galactic Centre Survey, focusing particularly on Sgr A, Sgr A East, the Radio Arc region, and an X-ray non-thermal filament designated XMM J174540$``$2904.5.

## 2 Global image

Fig. 2 shows a false colour X-ray image of the Galactic Centre Region taken with XMM-Newton in which Galactic north is up. Most of red dots are presumably foreground stars. Two persistent (1E 1740.7$``$2942 and 1E 1743.1$``$2843) and one transient (SAX J1747.0$``$2853) X-ray binaries are apparent as very bright sources in these observations. There is extended high surface brightness emission coincident with the Sgr A and Radio Arc (or G0.1$``$0.1) regions. The Arches star cluster and SNR G0.9$`+`$0.1 are also prominent in the hard energy band. In addition to these identifiable sources, low-surface brightness emission is evident over the whole region. The configuration is, however, highly asymmetric with the region to the Galactic west of Sgr A being brighter and having a more complex morphology than that to the east. We also point out that there is an extended and clumpy structure extending from Sgr A towards the (Galactic) south-west direction on a scale of several tens of parsecs.

## 3 Sgr A

With the moderate angular resolution of XMM-Newton, Sgr A cannot be resolved. However, we have detected an X-ray flare from Sgr A. Since the flare was detected in the last tens of minutes of an observation and was increasing in brightness right up to the end of the observation, the peak luminosity is unknown; the maximum luminosity observed was 4$`\times 10^{34}`$erg s-1. The flare nature is basically consistent with those detected by Chandra (Baganoff et al. 2001). Full details are given in Goldwurm et al. (2002).

## 4 Sgr A East

### 4.1 Images

The Sgr A region which is characterised by extremely bright radio structure, consists of Sgr A East and West. Sgr A East is an oval shell-like radio morphology, and in projection surrounds Sgr A West, which includes Sgr A. From its shell-like structure, Sgr A East is supposed to be a supernova remnant (SNR), SNR 000.0+00.0, although alternative interpretations have also been proposed.

Fig. 2 shows the XMM-Newton MOS1+2 image of the Sgr A region. The X-ray emitting region is contained within the radio shell of Sgr A East, as is also readily apparent from the Chandra observation (Maeda et al. 2002). We also made a 6.7-keV line narrow-band image corresponding to helium-like iron, since the X-ray spectrum of Sgr A East apparently shows that line (see Sec.4.2). The underlying continuum is subtracted using an adjacent bandpass and assuming an averaged spectral shape (Fig. 2 right panel). The 6.7-keV line is found to be clearly more concentrated in the core of Sgr A East than the continuum. This implies that the core of Sgr A East is more abundant in iron or higher in temperature, or perhaps a combination of these factors (see Sec.4.2 for detail).

### 4.2 Spectra

We accumulated the source spectrum of Sgr A East from a region of radius 100′′ (see Fig. 2), but excluding both Sgr A West and the region around a bright soft point source. In the spectrum, several strong emission lines at energies corresponding to K$`\alpha `$ lines from highly ionized ions, can be seen (Sakano et al. 2002b).

From the line ratios of K$`\alpha `$ lines from helium(He)-like and hydrogen(H)-like atoms, we found that the spectrum requires a multi-temperature plasma. When we approximate the spectrum with a two-temperature thin thermal plasma model by fitting all the spectra (MOS1, 2, and pn) simultaneously, we obtained the best-fitting temperatures to be $``$1 keV and $``$4 keV. The best-fit abundances for silicon, sulfur, argon, calcium, and iron vary slightly from element to element, but are mostly in the range 1–2 solar units.

We further examined the possible spectral variation within Sgr A East using spatially resolved spectra, and found that only the iron abundance varies significantly from place to place: $``$0.5 solar in the outer region ($`r>60^{\prime \prime }`$) to $``$3.5 solar in the central region ($`r<28^{\prime \prime }`$). Thus, the iron line concentration in Fig. 2 is due to an abundance gradient.

### 4.3 Discussion

The concentration of the iron at the centre suggests that the plasma originated in ejecta. From the derived spectral parameters and spatial extent, we estimated the total energy and mass to be $`1.5\times 10^{49}\eta ^{1/2}`$ erg and 1.9$`\eta ^{1/2}`$ M, respectively, assuming pressure equilibrium between the low- and high-temperature components, where $`\eta `$ is a filling factor of the total plasma. These results, therefore, favour an origin for Sgr A East in a type-II SN event. The most intriguing problem is the extremely high temperature of 4 keV. The special environment in the Galactic Centre Region may be responsible for it. The detailed analysis and discussion are found in Sakano et al. (2002b).

## 5 The Radio Arc region

This region is unique due to the bright structure both in the radio and X-ray bands. The brightest X-ray emission is coincident with the position of the molecular cloud G 0.1$``$0.1. On the other hand, the radio brightest structure, the Radio Arc, is not particularly luminous in the X-ray band. Another unusual characteristic is the extremely high surface brightness of the 6.4-keV line, a fluorescent K$`\alpha `$ line from neutral iron. Whereas the X-ray reflection nebula, Sgr B2, is dim in the continuum but bright only in the 6.4-keV line (Murakami et al. 2001), this region is bright both in the continuum and 6.4-keV line. The origin of this X-ray structure remains an open question.

## 6 An X-ray non-thermal filament XMM J174540$``$2904.5

A hard extended source, XMM J174540$``$2904.5, newly discovered with XMM, is marked in Fig. 2 This source has a corresponding radio non-thermal structure, the Sgr A-E ‘wisp’ (Ho et al. 1985) (= 1LC 359.888$``$0.086 = G359.88$``$0.07). Both the X-ray and radio images show similar elongated structure, although the X-ray emission is more localised. We found that the X-ray emission from this source is clearly non-thermal with an energy index of 1.0$`{}_{0.9}{}^{}{}_{}{}^{+1.1}`$ and heavily absorbed ($`N_\mathrm{H}=38_{11}^{+7}\times 10^{22}`$H cm$`{}_{}{}^{}2`$). So far, this is the first X-ray/radio filament in the Galactic Centre Region which has been shown, unequivocally, to have a non-thermal X-ray spectrum.

The spectral energy distribution shows that the X-ray spectrum is steeper than that observed in the radio band but is, most likely, smoothly connected to it. Thus, the emission is presumably due to synchrotron radiation. The lifetime of the high energy particles ($``$7 yr) emitting X-rays is roughly comparable with the spatial extent of the structure. Hence, the X-ray source is probably the site of particle acceleration. The position of the X-ray source coincides with the peak of the molecular cloud and it is quite possible that dense environment is one factor in determining this as an efficient acceleration site. This represents one of only a few X-ray/radio non-thermal filament known in the Galactic Centre Region. The details are discussed in Sakano et al. (2002a).

M.S. are financially supported by JSPS.


text:
# References

ABOUT SOME RELATIVISTIC EFFECTS IN THE $`A(\stackrel{}{d},p)X`$ REACTION WITH THE PROTON EMISSION IN A FORWARD DIRECTION

L.S. Azhgirey1, N.P. Yudin2,

1 JINR, 141980 Dubna, Moscow Region, Russia
2 Moscow State University, 119899 Moscow, Russia

## Abstract

In a nonrelativistic quantum theory the longitudinal and transversal components of the internal momentum in the deuteron wave function (DWF) are bound together as they have to form the well known superposition of the S- and D-waves. In a relativistic case — at least in light cone dynamics — the dependence of the DWF on the longitudinal and transversal momenta can by considerably changed, because these components turn out in large part to be bound in another way. Previously such a possibility was pointed out by Blankenbecler et al. Later Karmanov et al. have developed a quantitative approach that can be used to describe this situation. Here we present some results of calculations of the tensor analyzing power of $`(d,p)`$ reaction at relativistic deuteron momenta using Karmanov’s DWF. The momentum behaviour of the tensor analyzing power obtained is compared with that calculated with standard DWFs.

In recent years attempts to gain some insight into a manifestation of the structure of the relativistic deuteron in nuclear interactions are concerned mainly with an investigation of spin phenomena in reactions initiated by relativistic polarized deuterons:

$`\stackrel{}{d}+p`$ $``$ $`p+d,`$ (1)
$`\stackrel{}{d}+p`$ $``$ $`p+p+n,`$ (2)
$`\stackrel{}{d}+p(A)`$ $``$ $`p+X,`$ (3)
where detected protons are emitted in forward directions in the lab. syst. The simplest mechanism of these processes can be represented schematically by the diagram shown in Fig. 1. When protons are detected at non-zero angles, they may have a rather high transversal component of the momentum. At usually considered mechanisms these reactions are concerned with the deuteron wave function.

Already upon the first measurements of the tensor analyzing power $`T_{20}`$ of reaction (3) at 9 GeV/$`c`$ with the emission of protons at 0 \[1-5\] a significant discrepancy between the values calculated in the relativistic impulse approximation and experimental ones has come to light. A comparison of currently available experimental data on $`A_{yy}(0^{})`$ (recall, that at $`0^{}`$ $`A_{yy}=T_{20}/\sqrt{2}`$) with calculations in the frame of the impulse approximation with the standard DWFs  is made in fig. 2. More recently, on the measurement of the tensor analyzing power $`A_{yy}`$ of reaction (3) at 9 GeV/$`c`$ with the emission of protons with large transverse momenta , this discrepancy has been compounded.

> Fig. 1. The simplest diagram to describe reactions initiated by relativistic deuterons with proton emission in a forward direction. Along with the one-nucleon exchange ($`N`$) there also the exchange of baryon resonances ($`R`$) can take place .

The discrepancy between the theory and experiment in principle can be overcome by the addition of a P-wave into the deuteron wave function (DWF). True enough, the P-wave is conditioned by the different mechanisms in a number of approaches. In the paper  the data on $`T_{20}`$ at 0 have been successfully described taking account of the P-wave arising because of the production of six-quark configuration that gives odd parity resonances at their fragmentation to the baryon channel. At the same time the calculations on the basis of the Bethe-Salpeter equation give the P-wave admixture that is inadequate to eliminate the discrepancy with experiment . In relation to the reaction (1) it is notable that taking account of an admixture of the odd parity baryon resonances to the DWF allows the $`T_{20}`$ data at 0 to be better described than within the usual one-nucleon exchange approximation .

In this report we would like to point out one more factor related to the wave function of a fast moving deuteron. These are the rotational properties of the relativistic DWF in the light cone (LC) dynamics. . Let us consider shortly the calculations in the frame of this approach.

A general expression for the parameter $`T_{20}`$ has the form

$$T_{20}=\frac{_{M,M^{}}Sp\{\psi _Mt_{20}\psi _M^{}^{}\}}{(1/3)Sp\{\psi _M\psi _M^{}^{}\}},$$
(4)
where the operator $`t_{20}`$ is defined by expression

$$<m^{}|t_{20}|m>=(1)^{jm^{}}<1m1m^{}|\mathrm{2\; 0}>,$$
(5)
$`\psi _M`$ is DWF with spin $`J=1`$ and its projections on the quantization axis $`M=0,\pm 1`$, $`\psi ^{}`$ is the Hermitian conjugate function, and $`<1m1m^{}|\mathrm{2\; 0}>`$ is a Clebsh-Gordan coefficient. In order to calculate $`T_{20}`$ we should know $`\psi _M`$ with different projections $`M`$. However $`\psi _M`$ in LC dynamics has several unusual properties. The first one consists in ”strange” relations of $`\psi _M`$ with different $`M`$. Let us suppose that we have $`\psi _M`$ with $`M=0`$. Then $`\psi _M`$ with $`M=+1`$ should be obtained from $`\psi _M`$ with $`M=0`$ by the action of the standard combination of operators $`J_x,J_y`$:

$$\psi _{M=1}(J_x+iJ_y)\psi _{M=0}.$$
(6)
But in LC dynamics these operators of angular momentum are dependent on the interaction, and to calculate the sum in the expression for $`T_{20}`$ the Schrödinger equation for all components of $`\psi _M`$ has to be solved independently. Of course it is an objectionable problem. The dependence of angular momentum on interaction involves difficulties in forming states with definite angular momentum.

The second complication centres around the dependence of of the wave function on variables. Indeed, the wave function of two particles is in general dependent on two variables,

$$\psi =\psi (k_l,k_T),$$
(7)
where $`k_l`$ and $`k_T`$ are longitudinal and transversal components of momentum, respectively. In the nonrelativistic physics $`k_l`$ and $`k_T`$ are rigidly bound together by the necessity to form the usual superposition of $`S`$\- and $`D`$-waves in DWF, with the result that the nonrelativistic DWF is dependent on the only variable $`k`$. In LC dynamics this relation between $`k_l`$ and $`k_T`$ may be much weaker, and in principle DWF can be dependent on two variables. Attention to such a situation was first called by Schmidt and Blankenbecler .

These difficulties of the LC-dynamics were to a great extent overcame by Karmanov et al. . At first sight the approach of Karmanov complicates the situation since there is introduced an additional dynamical variable $`\omega ^\mu `$ — a four-vector of the normal to the chosen light cone surface. But this new variable permits the necessity dealing with operators of the angular momentum dependent on the interaction to be excluded: it turns out that the rotation of the light cone surface is equivalent to the interaction dependent rotation of the system. Therefore states with the definite angular momentum become formally independent on the interaction.

The dependence on the interaction does not completely disappear, it arises in another context: there are two many states with the same angular momentum and the role of the interaction now is to choose the right combination among these states. But this new problem is much more easier than the problem of solving the Schrödinger equation. The difficulties existing in Karmanov’s approach were overcame in , where the parameters of the relativistic DWF were given. This function was used in our calculations of the tensor analysing power $`A_{yy}`$ of the deuteron breakup reaction (3).

The DWF, now the function of nucleon spins, an internal deuteron momentum, and an orientation of the quantization plane, has the form:

$$\mathrm{\Psi }_{\sigma _2\sigma _1}^M=w_{\sigma _2}^{}\psi ^M(\stackrel{}{k},\stackrel{}{n})\sigma _yw_{\sigma _1},$$
(8)
where

$`\stackrel{}{\psi }(\stackrel{}{k},\stackrel{}{n})=f_1{\displaystyle \frac{1}{\sqrt{2}}}\stackrel{}{\sigma }+f_2{\displaystyle \frac{1}{2}}[{\displaystyle \frac{3}{\stackrel{}{k}^2}}\stackrel{}{k}(\stackrel{}{k}\stackrel{}{n}\stackrel{}{\sigma }]+f_3{\displaystyle \frac{1}{2}}[3\stackrel{}{n}(\stackrel{}{n}\stackrel{}{\sigma })\stackrel{}{\sigma }]+`$ (9)
$`f_4{\displaystyle \frac{1}{2k}}(3\stackrel{}{k}(\stackrel{}{n}\stackrel{}{\sigma })+3\stackrel{}{n}(\stackrel{}{k}\stackrel{}{\sigma })2(\stackrel{}{k}\stackrel{}{n})\stackrel{}{\sigma })+f_5\sqrt{{\displaystyle \frac{3}{2}}}{\displaystyle \frac{i}{k}}(\stackrel{}{k}\times \stackrel{}{n})+f_6{\displaystyle \frac{\sqrt{3}}{2k}}((\stackrel{}{k}\times \stackrel{}{n})\times \stackrel{}{\sigma }).`$
Here $`\stackrel{}{n}`$ is the normal to the light front surface, $`\stackrel{}{k}`$ is the momentum of nucleons in deuteron in their rest frame, $`\stackrel{}{\sigma }`$ are the Pauli matrices, $`w_{\sigma _1(\sigma _2)}`$ are the spin functions of nonrelativistic nucleons, and $`f_{1,\mathrm{},6}`$ are the invariant about rotations functions of the kinematical variables, that define the deuteron state.

> Fig. 2. Tensor analyzing power $`A_{yy}`$ of the deuteron breakup with the emission of protons at 0 vs internal momentum $`k`$. Experimental data: triangles — , squares — , diamonds — , empty circles —  ($`H`$ target), full circles —  ($`C`$ target). Calculations are made in relativistic impulse approximation with Paris  (broken curve), Bonn  (dotted curve) and Karmanov’s  DWFs.

Previously we described $`dppX`$ and $`dppd`$ reactions  in the coordinate system where the $`z`$ axis was aligned with the beam direction. For the wave function (8), too, let us assume that $`\stackrel{}{n}`$ is directed along $`z`$ axis, i.e. $`\stackrel{}{n}=(0,0,1)`$.

The expression (4) may be written in the form:

$$T_{20}=\sqrt{\frac{3}{2}}\frac{[3Sp(\psi _z\psi _z^{})+Sp(\stackrel{}{\psi }\stackrel{}{\psi }^{})]}{Sp(\stackrel{}{\psi }\stackrel{}{\psi }^{})},$$
(10)
where we went from spherical notations to the vector ones. For the traces we have:

$`Sp(\stackrel{}{\psi }\stackrel{}{\psi }^{})`$ $`=`$ $`3[f_1^2+f_2^2+(3z^21)f_2f_3+f_3^2+4zf_4(f_2+f_3)+`$ (11)
$`f_4^2(z^2+3)+(1z^2)f_5^2+(1z^2)f_6^2],`$
$`Sp(\psi _z\psi _z^{})`$ $`=`$ $`f_1^2+\sqrt{2}[3(kz)^21]f_1f_2+2\sqrt{2}f_1f_3+4\sqrt{2}kzf_1f_4+`$
$`{\displaystyle \frac{1}{2}}[3(kz)^2+1]f_2^2+2[3(kz)^21]f_2f_3+kz[3(kz)^2+5]f_2f_4+`$
$`2f_3^2+8kzf_3f_4+3sqrt3kz[1(kz)^2]f_2f_6+3sqrt3[1(kz)^2]f_4f_6+`$
$`\{[{\displaystyle \frac{9}{2}}[(kz)^2+1](kz)^2\}f_4^2+{\displaystyle \frac{3}{2}}[1(kz)^2]f_6^2,`$
where $`z=cos\theta `$.

The results of calculations of $`A_{yy}(0^{})`$ with the wave function (9) are shown in fig. 2 by the solid curve. It is seen that as opposed to the calculations with the standard nonrelativistic DWFs  the solid curve does not cross the horizontal axes and is in reasonable good agreement with experimental data in the region of $`k`$ from 0.4 to 0.8 GeV/$`c`$. This result in our opinion is due to the fact that Karmanov’s model establishes a new link between $`k_l`$ and $`k_T`$ that is different from those of the S- and D-wave superposition in nonrelativistic DWF. Similar effect was discussed in  on the example of Wick-Cutcosky model, where it was shown in the clear form that a $`S`$-wave two-particle system becomes dependent on the angle in the LC dynamics. In other words it is a manifestation of the intimate connection between the internal motion and the motion of the system as a whole.

This work was supported in part by the Russian Foundation for Fundamental Research (grant No. 01-02-17299).


text:
# Exact scaling functions for one-dimensional stationary KPZ growth

## 1 Introduction

In their well-known work  Kardar, Parisi, and Zhang argue that surface growth through random ballistic deposition can be modeled by a stochastic continuum equation, which in the case of a one-dimensional substrate reads

$$_th=\frac{1}{2}\lambda (_xh)^2+\nu _x^2h+\eta .$$
(1.1)
Here $`h(x,t)`$ is the height at time $`t`$ at location $`x`$ relative to a suitable reference line. $`\eta (x,t)`$ is space-time white noise of strength $`D`$, $`\eta (x,t)\eta (x^{},t^{})=D\delta (xx^{})\delta (tt^{})`$, and models the randomness in deposition. $`\nu _x^2h`$ is a not further detailed smoothening mechanism. The important insight of  is to observe that the growth velocity is nonlinear, in general, and is relevant for the large scale properties of the solution to (1.1). To simplify, the growth velocity is expanded in the slope. The first two terms can be absorbed through a suitable choice of coordinate frame. The quadratic nonlinearity in (1.1) is relevant and higher orders can be ignored, unless $`\lambda =0`$.

The one-dimensional KPZ equation (1.1) is regarded as exactly solved in the usual terminology. In fact, what can be obtained is the dynamic scaling exponent $`z=3/2`$ . No other universal quantity has been computed exactly so far. In our contribution we will improve the situation and explain how to extract the scaling function for the stationary two-point function. A few other universal quantities can be computed as well. But they have been discussed already elsewhere .

In  a mode-coupling equation for the two-point function is written down, in essence following the scheme from critical dynamics and kinetic theory. At the time only $`z=3/2`$ and a few qualitative properties could be extracted from the mode-coupling equation. In  this equation is solved numerically. Such computations are repeated in  with greatly improved precision and using a more convenient set of coordinates. Thus for the 1D KPZ equation we are in the unique position of an exact solution and an accurate numerical solution to the mode-coupling equation with no adjustable parameters. As will be explained below, given the uncontrolled approximation, mode-coupling does surprisingly well.

To attack (1.1) directly does not seem to be feasible, a situation which is rather similar to the one for two-dimensional models in equilibrium statistical mechanics. For example, the Ginzburg-Landau $`\varphi ^4`$-theory is given through the (formal) functional measure

$$Z^1\underset{x^2}{}d\varphi (x)\mathrm{exp}\left[d^2x\left((\varphi )^2+g\varphi ^2+\varphi ^4\right)\right]$$
(1.2)
for the scalar field $`\varphi `$. (1.2) is not the proper starting point for computing the exact two-point scaling function at the critical coupling $`g_c`$. Rather one discretizes through the lattice $`^2`$ and replaces the $`\varphi `$-field by Ising spins $`\pm 1`$. Then, following e.g. , the scaling function at and close to criticality can be obtained. By universality this scaling function is the one of (1.2). (While certainly true, to establish universality is difficult and carried out in a few cases only .) In the same spirit we replace (1.1) by a discrete model, where the most convenient choice seems to be the polynuclear growth (PNG) model.

Before explaining the PNG model let us review the standard scaling theory for (1.1). If the initial conditions $`h(x,0)`$ of the KPZ equation are distributed according to two-sided Brownian motion, then formally the distribution of $`h(x,t)h(0,t)`$ is again two-sided Brownian motion. Therefore it is natural to define the stationary time correlation

$$C(x,t)=\left(h(x,t)h(0,0)t_th\right)^2,$$
(1.3)
where from the height difference the average displacement is subtracted. By assumption

$$C(x,0)=A|x|$$
(1.4)
with roughness amplitude $`A=D/\nu `$ to ensure stationarity in time. If $`z=3/2`$, then $`C(x,t)`$ scales as

$$C(x,t)t^{2/3}g(constx/t^{2/3}),\text{as }x,t\mathrm{}$$
(1.5)
with a universal scaling function $`g(y)`$ having the asymptotics $`g(y)c_0>0`$ for $`y0`$ and $`g(y)c_{\mathrm{}}|y|`$ for $`|y|\mathrm{}`$. In order to define $`g`$ as a dimensionless function we fix the proportionality constants in (1.5) as appropriate combinations of $`\lambda `$ and $`A`$,

$$g(y)=\underset{t\mathrm{}}{lim}\frac{C((2\lambda ^2At^2)^{1/3}y,t)}{(\frac{1}{2}\lambda A^2t)^{2/3}},$$
(1.6)
where the particular choice of numerical prefactors is in principle arbitrary. The factor $`2^{2/3}`$ in the denominator is chosen in order to conform with the convention for the GUE Tracy-Widom distribution . The factor $`2^{1/3}`$ in the argument of the numerator differs from the convention used by Baik and Rains  by a factor $`2`$ but conforms with the definition of the closely related Airy process  and has the further advantage to absorb a lot of prefactors in the equations defining $`g(y)`$. Note however, that the exponents for the parameters $`\lambda \left[\frac{x^2}{th}\right]`$, $`\nu \left[\frac{x^2}{t}\right]`$, and $`D\left[\frac{h^2x}{t}\right]`$ are fixed uniquely by dimensional reasoning.

We remark that the slope $`_xh(x,t)`$ is space-time stationary in the usual sense. For fixed $`t`$, $`x_xh(x,t)`$ is white noise with strength $`A`$. Since $`_xh=0`$, the standard 2-point function is

$$_xh(0,0)_xh(x,t)=\frac{1}{2}_x^2C(x,t).$$
(1.7)
This relation and the asymptotic behavior of $`g`$, $`g(y)/|y|2`$ as $`y\mathrm{}`$, motivates the definition of a second scaling function,

$$f(y)=\frac{1}{4}g^{\prime \prime }(y),$$
(1.8)
which by definition has integral normalized to one and which will be shown to be positive in the next section.

In the sequel we will analyze the distribution function for the height differences in the stationary PNG model. As shown in , they can be represented in terms of certain orthogonal polynomials, which lead to recursion relations connected to the Painlevé II differential equation . The asymptotic analysis is carried out in . Our own contribution is twofold: (i) We observe that the stationary PNG model maps to a last passage percolation with boundaries . (ii) The expressions in  are given in terms of certain differential equations and the extraction of the scaling function $`g`$ requires a careful numerical integration. This is one central point of our article. We will provide then plots of the structure function and give a comparison with the mode-coupling theory.

## 2 The polynuclear growth model

The polynuclear growth (PNG) model is a model for layer-by-layer growth through deposition from the ambient atmosphere. The surface is parameterized by a time dependent integer-valued height function $`h(x,t)`$, $`t`$, above a one-dimensional substrate, $`x`$. Thus the height function consists of terraces bordered by steps of unit height. The up-steps move to the left and the down-steps to the right with speed $`1`$. Steps disappear upon collision. In addition to this deterministic dynamical rule new islands of unit height are nucleated randomly with space-time density $`2`$ on top of already existing terraces. The corresponding stochastic process $`h(x,t)`$ is well defined even in infinite volume (cf.  for the closely related Hammersley particle process).

Of interest to us here is the stationary growth process, which means that the slope $`_xh(x,t)=\rho (x,t)`$ is stationary in space-time. One can think of $`\rho (x,t)`$ as the density of a particle/antiparticle process. The particles are located at the up-steps and thus move with velocity $`1`$, the antiparticles are located at the down-steps and move with velocity $`1`$. Upon collision particle/antiparticle pairs annihilate. In addition, with space-time density $`2`$, a particle/antiparticle pair is created with the particle moving to the left, the antiparticle to the right. To make $`\rho `$ stationary, one prescribes at $`t=0`$ up-steps Poisson distributed with density $`\rho _+`$ and down-steps independently Poisson distributed with intensity $`\rho _{}`$ such that

$$\rho _+\rho _{}=1.$$
(2.1)
This measure for steps is stationary under the PNG dynamics. The mean slope is given by

$$u=\rho _+\rho _{}=_xh(x,t),$$
(2.2)
which is the only remaining free parameter. For fixed $`t`$, $`xh(x,t)h(0,t)`$ is a (two-sided) random walk with rate $`\rho _\pm `$ for a jump from $`n`$ to $`n\pm 1`$. It has average $`u`$ and variance $`\rho _++\rho _{}`$, which implies for the roughness amplitude

$$A(u)=\sqrt{4+u^2}.$$
(2.3)
For the growth velocity one obtains

$$v(u)=_th=\rho _++\rho _{}=\sqrt{4+u^2}.$$
(2.4)
Given $`\rho (x,t)`$ the height $`h(x,t)`$ is determined only up to a constant which we fix as $`h(0,0)=0`$. To emphasize that only height differences count, $`h(0,0)`$ is sometimes kept in the formulas.

The stationary process with slope $`u`$ transforms to the stationary process with slope $`0`$ through the Lorentz transformation

$$x^{}=(1c^2)^{1/2}(xct),t^{}=(1c^2)^{1/2}(tcx),$$
(2.5)
with the speed of “light” equal to $`1`$ and the velocity parameter $`c=u/\sqrt{4+u^2}`$. Thus it suffices to restrict ourselves to $`u=0`$ which we do from now on. In particular $`\rho _+=1=\rho _{}`$. $``$ and $`𝔼`$ refer to the stationary density field at slope $`u=0`$.

The central objects are the height-height correlation

$$C(x,t)=(h(x,t)h(0,0)2t)^2$$
(2.6)
and the closely related two-point function for the density,

$$S(x,t)=\rho (x,t)\rho (0,0).$$
(2.7)
They are related as

$`\frac{1}{2}_x^2C(x,t)`$ $`=`$ $`\frac{1}{2}_x^2𝔼\left(\left(h(x,t)h(0,0)2t\right)^2\right)`$ (2.8)
$`=`$ $`_x𝔼\left(\rho (x,t)\left(h(x,t)h(0,0)2t\right)\right)`$
$`=`$ $`_x𝔼\left(\rho (0,t)\left(h(0,t)h(x,0)2t\right)\right)`$
$`=`$ $`𝔼\left(\rho (0,t)\rho (x,0)\right)`$
$`=`$ $`S(x,t).`$
The height correlation is convex, equivalently

$$S(x,t)0.$$
(2.9)
To prove this property we show that the structure function $`S(x,t)`$ can be regarded as the transition probability for a second class particle starting at the origin. Its initial velocity is $`\pm 1`$ with probability $`\frac{1}{2}`$, as for the “first-class” up/down-steps. In contrast to an ordinary step the second class particle is never destroyed upon colliding with another step. Rather it eats up the step encountered and, by reversing its own direction of motion, continues along the trajectory of the absorbed step, cf. Figure 1.

Let $`\rho (x,t)`$ be a given realization of the PNG process. The second class particle is added as

$$\rho ^{(\sigma )}(x,0)=\rho (x,0)+\sigma \delta (x),\sigma =\pm 1.$$
(2.10)
$`\rho ^{(\sigma )}(x,0)`$ evolves to $`\rho ^{(\sigma )}(x,t)`$ with nucleation events identical to the one for $`\rho (x,t)`$. By construction, if $`X_t`$ denotes the position of the second class particle at time $`t`$,

$$\rho ^{(\sigma )}(x,t)\rho (x,t)=\sigma \delta (xX_t).$$
(2.11)
Noting that by the Poisson property $`\rho ^{(\sigma )}(x,0)`$ is given by $`\rho (x,0)`$ conditioned on the presence of either an up-step ($`\sigma =+1`$) or down-step ($`\sigma =1`$) at the origin, we obtain

$`0p_t(x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma =\pm 1}{}}𝔼\left(\sigma \left(\rho ^{(\sigma )}(x,t)\rho (x,t)\right)\right)={\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma =\pm 1}{}}𝔼\left(\sigma \rho ^{(\sigma )}(x,t)\right)`$ (2.12)
$`=`$ $`\underset{\delta 0}{lim}{\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma =\pm 1}{}}\sigma 𝔼\left(\rho (x,t)|{\displaystyle _\delta ^\delta }\rho (y,0)𝑑y=\sigma \right)`$
$`=`$ $`\underset{\delta 0}{lim}{\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma ,\sigma ^{}=\pm 1}{}}\sigma \sigma ^{}{\displaystyle \frac{\left\{_\delta ^\delta \rho (x+y,t)𝑑y=\sigma ^{},_\delta ^\delta \rho (y,0)𝑑y=\sigma \right\}}{2\delta \left\{_\delta ^\delta \rho (y,0)𝑑y=\sigma \right\}}}`$
$`=`$ $`\frac{1}{2}𝔼\left(\rho (x,t)\rho (0,0)\right)=\frac{1}{2}S(x,t).`$
For arbitrary slope $`u`$ the normalization of $`S(x,t)`$ would be given by $`v(u)=\sqrt{4+u^2}`$ and the mean of $`p_t(x)`$ evolves along the characteristics of the macroscopic evolution equation $`_tu=_xv(u)`$. Thus

$$S(x,t)𝑑x=\sqrt{4+u^2},\text{and}xS(x,t)𝑑x=tu.$$
(2.13)
## 3 The distribution functions for the height differences

For the PNG model the distribution function for the height difference $`h(x,t)h(0,0)`$ satisfies certain recursion relations, which are the tool for analyzing the scaling limit when $`t\mathrm{}`$ and $`x=yt^{2/3}`$ with $`y=𝒪(1)`$. The second moments yield $`C(x,t)`$ and therefore by (2.8) also $`S(x,t)`$.

Since the nucleation events are Poisson, $`h(x,t)h(0,0)`$ depends only on the events in the backward light cone $`\{(x^{},t^{})^2,\mathrm{\hspace{0.17em}0}t^{}t,|xx^{}|tt^{}\}`$ and the initial conditions at $`t=0`$. Along the line $`\{x^{}=t^{}\}`$ the down-steps are Poisson distributed with line density $`\sqrt{2}`$ and correspondingly for the up-steps along the line $`\{x^{}=t^{}\}`$. This property can be deduced from the uniqueness of the stationary state at given slope and the Lorentz invariance (2.5) in the limit $`c\pm 1`$. Thus $`h(x,t)h(0,0)`$ is determined by the nucleation events in the rectangle $`R_{x,t}=\{(x^{},t^{})^2,|x^{}||t^{}|,|xx^{}|tt^{}\}`$ together with the said boundary conditions. $`h(x,t)h(0,0)`$ can be reexpressed as a directed last passage percolation according to the following rules: Inside $`R_{x,t}`$ there are Poisson points with density $`2`$. Along the two lower edges of $`R_{x,t}`$ there are independently Poisson points with line density $`\sqrt{2}`$. A directed passage from $`(0,0)`$ to $`(x,t)`$ is given through a directed path (polymer). It is a piecewise linear path in the plane, starts at $`(0,0)`$ and ends at $`(x,t)`$, alters its direction only at Poisson points, and is time-like in the sense that for any two points $`(x^{},t^{})`$, $`(x^{\prime \prime },t^{\prime \prime })`$ on the path one has $`|x^{}x^{\prime \prime }||t^{}t^{\prime \prime }|`$. Note that, once the directed path leaves one of the lower edges to move into the bulk, it can never return. By definition the length of a directed path equals the number of Poisson points traversed. With these conventions

$$h(x,t)h(0,0)=\text{maximal length of a directed path from }(0,0)\text{ to }(x,t).$$
(3.14)
We remark that in general there are several maximizing paths, their number presumably growing exponentially with $`t`$.

Under the Lorentz transformation (2.5) the distribution for the height differences (3.14) does not change. Therefore, we might as well transform $`R_{(x,t)}`$ to a square. By an additional overall scaling by $`\sqrt{2}`$ one arrives at a $`v\times v`$ square, $`v=\sqrt{t^2x^2}`$, with bulk density $`1`$ and the line densities $`\alpha _{}=\alpha =\sqrt{(tx)(t+x)}`$ for the lower left, resp. $`\alpha _+=1/\alpha `$ for the lower right edge. In this way we have recovered precisely the setting in  with $`t`$ replaced by $`v`$. Baik and Rains derive an explicit expression for the height distribution in terms of Toeplitz determinants, which can be further simplified by means of corresponding orthogonal polynomials.

Let us state the result for the distribution function of $`h(x,t)h(0,0)`$,

$`F_{x,t}(n)`$ $`=`$ $`\{h(x,t)h(0,0)n\}`$ (3.15)
$`=`$ $`G_n(\alpha )F(n)G_{n1}(\alpha )F(n1).`$
For fixed $`v`$ the functions $`g`$ and $`F`$ are given in terms of the monic polynomials $`\pi _n(z)=z^n+𝒪(z^{n1})`$, which are pairwise orthogonal on the unit circle $`|z|=1`$ with respect to the weight $`e^{v(z+z^1)}`$. Their norm $`N_n`$ is given by

$$\pi _n,\pi _m=\delta _{n,m}N_n$$
(3.16)
with $`p,q=p(z)q(z^1)e^{v(z+z^1)}(2\pi iz)^1𝑑z`$. One has

$$F(n)=e^{v^2}\underset{k=0}{\overset{n1}{}}N_k,$$
(3.17)
where $`F(n)`$ itself is the distribution function of the maximal length of a directed path in the case $`\alpha _+=0=\alpha _{}`$. Thus $`lim_n\mathrm{}F(n)=1`$ and

$`G_n(\alpha )`$ $`=`$ $`e^{v(\alpha +\alpha ^1)}N_n{\displaystyle \underset{k=0}{\overset{n}{}}}N_k^1\pi _k(\alpha )\pi _k(\alpha ^1)`$ (3.18)
$`=`$ $`e^{v(\alpha +\alpha ^1)}((1n)\pi _n(\alpha )\pi _n(\alpha ^1)`$
$`\alpha \pi _n^{}(\alpha )\pi _n(\alpha ^1)\alpha ^1\pi _n(\alpha )\pi _n^{}(\alpha ^1)).`$
Defining the dual polynomials $`\pi _n^{}(z)=z^n\pi _n(z^1)`$, the second equality in (3.18) is an easy consequence of the Christoffel-Darboux formula ,

$`N_n{\displaystyle \underset{k=0}{\overset{n1}{}}}{\displaystyle \frac{\pi _k(a)\pi _k(b)}{N_k}}`$ $`=`$ $`{\displaystyle \frac{\pi _n^{}(a)\pi _n^{}(b)\pi _n(a)\pi _n(b)}{1ab}},`$ (3.19)
valid for $`a,b`$, $`ab1`$ and extended by l’Hospital’s rule to $`ab=1`$, and the trivial relation

$$\pi _n^{}(z)z^1\pi _{n}^{}{}_{}{}^{}(z^1)+z\pi _n^{}(z)\pi _n(z^1)=n\pi _n(z)\pi _n(z^1)=n\pi _n^{}(z)\pi _n^{}(z^1).$$
(3.20)
Taking only the leading order of $`a`$ in (3.19) one obtains the well-known relations

$`\pi _{n+1}(z)`$ $`=`$ $`z\pi _n(z)+p_{n+1}\pi _n^{}(z),`$
$`\pi _{n+1}^{}(z)`$ $`=`$ $`zp_{n+1}\pi _n(z)+\pi _n^{}(z),`$ (3.21)
$`N_{n+1}`$ $`=`$ $`N_n\left(1p_{n+1}^2\right)`$
which are closed given $`p_n=\pi _n(0)`$ for $`n0`$. For the particular weight function $`e^{v(z+z^1)}`$ one can derive a nonlinear recursion relation for the $`p_n`$’s,

$$p_n=\frac{v}{n}(p_{n+1}+p_{n1})(1p_n^2),$$
(3.22)
with initial values $`p_0=1`$, $`p_1=\frac{I_1(2v)}{I_0(2v)}`$. $`I_k(2v)=(2\pi )^1_0^{2\pi }e^{ik\theta }e^{2v\mathrm{cos}\theta }𝑑\theta `$ is the modified Bessel function of order $`k`$ and thus $`N_0=I_0(2v)`$. Eq. (3.22) is the discrete Painlevé II equation. It has been derived in the context of orthogonal polynomials for the first time in , and later on more or less independently in . The differential equations for $`\pi _n`$, $`\pi _n^{}`$,

$`\pi _{n}^{}{}_{}{}^{}(z)`$ $`=`$ $`(n/z+v/z^2p_{n+1}p_nv/z)\pi _n(z)+(p_{n+1}v/zp_nv/z^2)\pi _n^{}(z)`$
$`\pi _{n}^{}{}_{}{}^{}(z)`$ $`=`$ $`(p_{n+1}v/z+p_nv)\pi _n(z)+(v+p_{n+1}p_nv/z)\pi _n^{}(z),`$ (3.23)
can be shown to hold by a tedious but straightforward induction, using (3) and (3.22). They are implicitly derived in , from which we learned their actual form. In  an integral expression is obtained for the derivative of orthogonal polynomials on the circle with respect to (up to some technical conditions) an arbitrary weight function. Specializing to the weight $`e^{v(z+z^1)}`$ results in a differential-difference equation equivalent to (3).

Of course, the mean of the probability distribution $`F_{x,t}(n)F_{x,t}(n1)`$ is $`2t`$ and its variance, the correlation function (2.6), is given by

$$C(x,t)=\underset{n0}{}\left(2(n2t)1\right)F_{x,t}(n).$$
(3.24)
Thus to establish (1.5), one has to understand the scaling properties of the distribution function $`F_{x,t}(n)`$. Let us introduce the new variables $`s`$, $`y`$ defined by

$`n`$ $`=`$ $`2v+v^{1/3}s,`$ (3.25)
$`x`$ $`=`$ $`v^{2/3}y,`$ (3.26)
where $`v=\sqrt{t^2x^2}`$ is regarded as fixed when varying $`n`$ and $`x`$. In  the different scaling variable $`w=\frac{1}{2}y`$ is used, which leads to a string of factors of $`2`$, avoided by our convention. Setting

$$R_n=(1)^np_n,$$
(3.27)
we rewrite (3.22) as

$$R_{n+1}2R_n+R_{n1}=\frac{(\frac{n}{v}2)R_n+2R_n^3}{1R_n^2}.$$
(3.28)
Under the scaling (3.26), $`R_n=v^{1/3}u\left(v^{1/3}(n2v)\right)+𝒪(v^1)`$, it becomes the Painlevé II equation

$$u^{\prime \prime }(s)=2u(s)^3+su(s),$$
(3.29)
in the limit $`v\mathrm{}`$. The starting value $`R_0=1`$ is consistent with the left asymptotics of $`u(s)`$ only if

$$u(s)\sqrt{s/2}\text{as }s\mathrm{},$$
(3.30)
which singles out the Hastings-McLeod solution to (3.29) . This particular solution will be denoted by $`u(s)`$ and we conclude that

$$u(s)=\underset{v\mathrm{}}{lim}v^{1/3}R_{[2v+v^{1/3}s]},$$
(3.31)
provided the limit exists (A complete proof is the main content of ). $`u(s)<0`$ and $`u`$ has the right asymptotics

$$u(s)\text{Ai}(s)\text{as }s\mathrm{}.$$
(3.32)
At this point we can derive the scaling limit for $`F(n)`$. It is the GUE Tracy-Widom distribution function

$$F_{\text{GUE}}(s)=e^{V(s)},V(s)=_s^{\mathrm{}}v(x)𝑑x,v(s)=\left(u(s)^2+s\right)u(s)^2u^{}(s)^2,$$
(3.33)
which appears already as the limiting height distribution for nonstationary curved selfsimilar growth . Since $`v^{}(s)=u(s)^2`$, one has $`v^{1/3}\mathrm{log}N_{[2v+v^{1/3}s]}v(s)`$ and $`F([2v+v^{1/3}s])F_{\text{GUE}}(s)`$ as $`v\mathrm{}`$.

Next we turn to the scaling limit for the orthogonal polynomials. For it to be nontrivial we set

$`P_n(\alpha )`$ $`=`$ $`e^{v\alpha }\pi _n^{}(\alpha ),`$
$`Q_n(\alpha )`$ $`=`$ $`e^{v\alpha }(1)^n\pi _n(\alpha ).`$ (3.34)
(3.26) implies $`\alpha =1v^{1/3}y+𝒪(v^{2/3})`$. Setting $`n`$ as in (3.25), we claim that

$$a(s,y)=\underset{v\mathrm{}}{lim}P_n(\alpha ),b(s,y)=\underset{v\mathrm{}}{lim}Q_n(\alpha ).$$
(3.35)
If so, the limit functions $`a`$, $`b`$ satisfy the differential equations

$`_sa`$ $`=`$ $`ub,`$
$`_sb`$ $`=`$ $`uayb,`$ (3.36)
as a consequence of (3), and

$`_ya`$ $`=`$ $`u^2a(u^{}+yu)b,`$
$`_yb`$ $`=`$ $`(u^{}yu)a+(y^2su^2)b,`$ (3.37)
as a consequence of (3). From (3) one immediately obtains $`\pi _n^{}(1)=(1)^n\pi _n(1)=_{k=1}^n(1R_k)`$. One has the limit $`_{k=1}^{\mathrm{}}(1R_k)=e^v`$ since $`e^{v^2+v}N_0_{k=1}^nN_k(1R_k)^1`$ has an interpretation as a probability distribution function . Therefore the initial conditions to (3) are

$$a(s,0)=b(s,0)=e^{U(s)},U(s)=_s^{\mathrm{}}u(x)𝑑x.$$
(3.38)
The scaling limit of $`G_n(\alpha )`$ as defined in Eq. (3.18) is the function

$`G(s,y)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^s}a(s^{},y)a(s^{},y)𝑑s^{}`$ (3.39)
$`=`$ $`a(s,y)_ya(s,y)b(s,y)_yb(s,y),`$
where the second equality can be verified by differentiation with respect to $`s`$ and using the identity

$$a(s,y)=b(s,y)e^{\frac{1}{3}y^3sy},$$
(3.40)
itself being a direct consequence of (3) and (3.38). Putting these pieces together we obtain as scaling limit for the distribution functions $`F_{x,t}(n)`$,

$$F_y(s)=\frac{d}{ds}\left(G(s+y^2,y)F_{\text{GUE}}(s+y^2)\right).$$
(3.41)
The shift in (3.41) by $`y^2`$ comes from the fact that $`F_y(s)`$ is evaluated for constant $`t=v+\frac{1}{2}v^{1/3}y+𝒪(v^{1/3})`$.

In conclusion we arrive at the scaling function $`g(y)`$ as defined in the Introduction. From (1.6), with $`\lambda =\frac{1}{2}`$ and $`A=2`$ for the PNG model, and (3.24) we obtain

$$g(y)=s^2𝑑F_y(s).$$
(3.42)
As already mentioned, except for (3), all our relations are derived in , and the existence of limits is proven with Riemann-Hilbert techniques. For completeness let us collect some more properties of $`a`$ shown in :

$`a(s,y)`$ $``$ $`1,\text{as }s+\mathrm{},`$
$`a(s,y)`$ $``$ $`0,\text{as }s\mathrm{},`$
$`a((2y)^{1/2}x+y^2,y)`$ $``$ $`1,\text{as }y+\mathrm{},`$
$`a((2y)^{1/2}x+y^2,y)`$ $``$ $`{\displaystyle \frac{1}{(2\pi )^{1/2}}}{\displaystyle _{\mathrm{}}^x}e^{\frac{1}{2}\xi ^2}𝑑\xi ,\text{as }y\mathrm{}.`$ (3.43)
Therefore $`F_y(s)`$ is asymptotically Gaussian and we recover $`g(y)2|y|`$ for large $`y`$.

## 4 Numerical determination of the scaling function

The key object in determining the scaling functions $`g(y)`$, $`f(y)`$ is the Hastings-McLeod solution  to Painlevé II, $`u(s)`$, which is the unique solution to

$$u^{\prime \prime }=2u^3+su$$
(4.1)
with asymptotic boundary conditions (3.30) and (3.32). Tracy and Widom  integrate (4.1) numerically with conventional differential equation solvers using the known asymptotics at $`s=\pm \mathrm{}`$. The precision achieved with this technique does not suffice for our purposes, since we need $`u(s)`$ as starting values (3.38) for the differential equations (3). We develop here a different method to obtain $`u(s)`$, in principle with arbitrary precision. Next the functions $`a(s,y)`$ and $`b(s,y)`$ have to be determined, which directly leads to values for the distribution functions $`F_y(s)`$. They have to be further integrated with respect to $`s`$ in order to obtain their variance, which is the desired scaling function $`g(y)`$. The Taylor expansion method to be explained intrinsically produces not only function values at a point but also higher derivatives. Therefore we obtain $`f(y)`$ not by numerically differentiating $`g(y)`$ but rather by direct calculation via the knowledge of $`_y^2F_y(s)`$.

In a first step, to obtain reliable approximations to the Hastings-McLeod solution, we need to guess its initial data at some finite $`s_0`$ by using asymptotic expansions around $`\pm \mathrm{}`$. Any initial data $`u(s_0)=u_0`$, $`u^{}(s_0)=u_1`$ give rise to a maximal solution, $`\stackrel{~}{u}(s)`$ of (4.1), which admits analytic continuation to a meromorphic function on $``$. The only essential singularity for Painlevé II solutions lies at $`\mathrm{}`$. If $`\stackrel{~}{u}(s)`$ is close to $`u(s)`$ we can estimate the difference $`\mathrm{\Delta }(s)=\stackrel{~}{u}(s)u(s)`$ by linearizing (4.1) around $`u(s)`$. One obtains that $`\mathrm{\Delta }(s)/\mathrm{\Delta }(0)`$ is of order $`\mathrm{exp}(D(s))`$ if $`(\mathrm{\Delta }(s),\mathrm{\Delta }^{}(s))`$ is in the unstable subspace and of order $`\mathrm{exp}(D(s))`$ for the stable subspace, where $`D(s)\frac{1}{3}(2s)^{3/2}`$ for $`s0`$ and $`\frac{2}{3}s^{3/2}`$ for $`s0`$. Note that on the exponential scale we are looking at, derivation with respect to $`s`$ leaves invariant the order. Therefore the exponential orders of $`\mathrm{\Delta }(s)`$ and $`\mathrm{\Delta }^{}(s)`$ are the same. Since generically initial values at $`s_0`$ have a component in both subspaces one obtains that $`\mathrm{\Delta }(s)/\mathrm{\Delta }(s_0)`$ is of order $`\mathrm{exp}(|D(s)D(s_0)|)`$ in the range of validity of the linear approximation, $`|\mathrm{\Delta }(s)||u(s)|`$.

It turns out that the left asymptotic expansion in $`(s)^{1/2}`$, optimally truncated at large negative $`s_0`$, gives rise to initial values with $`\mathrm{\Delta }(s_0)`$ only of order $`\mathrm{exp}(\frac{1}{3}(2s_0)^{3/2})`$. Thus control of the approximation always breaks down near $`s=0`$. Approximations of the right asymptotics on the other hand allow a, in principle, arbitrary precision on any given finite interval.

For $`s\mathrm{}`$ the deviations of $`u(s)`$ from the Airy function can be expanded in an alternating asymptotic power series with exponentially small prefactor,

$$u_{\text{right},n}(s)=\text{Ai}(s)\frac{e^{2s^{3/2}}}{32\pi ^{3/2}s^{7/4}}\underset{k=0}{\overset{n}{}}\frac{(1)^ka_k}{(\frac{2}{3}s^{3/2})^k}.$$
(4.2)
The coefficients are $`a_0=1`$, $`a_1=\frac{23}{24}`$, $`a_2=\frac{1493}{1152}`$, $`\mathrm{}`$, and can be obtained via the recursion relation

$$a_n=\text{Ai}_n^{(3)}+\frac{3}{4}na_{n1}\frac{1}{8}(n\frac{1}{6})(n\frac{5}{6})a_{n2}\text{for }n0$$
(4.3)
with initial conditions $`a_1=a_2=0`$.

$$\text{Ai}_n^{(3)}=\underset{0kln}{}\text{Ai}_{nl}\text{Ai}_{lk}\text{Ai}_k.$$
(4.4)
are the coefficients in the asymptotic expansion of $`\text{Ai}(x)^3`$ and

$$\text{Ai}_n=\frac{(6n1)(6n5)}{72n}\text{Ai}_{n1},\text{Ai}_0=1$$
(4.5)
are the coefficients of the asymptotic expansion of the Airy function itself ,

$$\text{Ai}(s)\frac{e^{\frac{2}{3}s^{3/2}}}{2\sqrt{\pi }s^{1/4}}\underset{n0}{}\frac{(1)^n}{(\frac{2}{3}s^{3/2})^n}\text{Ai}_n.$$
(4.6)
Empirically we observe that for $`s_00`$ the optimal truncation in (4.2) is $`n\frac{4}{3}s_0^{3/2}`$ leading to an exponentially improved (relative) precision

$$\left|\frac{u_{\text{right},n}(s_0)u(s_0)}{u(s_0)}\right|\mathrm{exp}(\frac{8}{3}s_0^{3/2}).$$
(4.7)
The linear perturbation argument for $`\stackrel{~}{u}(s)`$ with initial values $`\stackrel{~}{u}(s_0)=u_{\text{right},n}(s_0)`$, $`\stackrel{~}{u}^{}(s_0)=u_{\text{right},n}^{}(s_0)`$, $`n=[\frac{4}{3}s_0^{3/2}]`$, is now valid for a slightly smaller interval than $`[2s_0,3^{3/2}s_0]`$. For example, by choosing the interval $`[2^{1/3}s_0,s_0]`$ and gluing $`\stackrel{~}{u}(s)`$ at the boundaries smoothly to the optimally truncated asymptotic expansions, the maximal relative error of $`\stackrel{~}{u}(s)`$ with respect to the Hastings-McLeod solution is of order $`e^{4/3s_0^{3/2}}`$. For our purpose it turns out that we do not need values with $`s<20`$. On the other hand, to access large values of $`y`$ in (3.41), because of the shift in (3.41), we need $`u(s)`$ for large $`s`$ with high precision. $`u_{\text{right},n}(s)`$ is numerically costly to evaluate for large $`s`$, so we finally choose $`s_0=100`$ and integrate in the interval $`[20,200]`$. This requires a maximal working precision of $`1500`$ digits and, given the integration of (4.1) is precise enough, $`\stackrel{~}{u}(s)`$ and $`u(s)`$ coincide in the first $`1000`$ digits for $`s[20,115]`$ and still up to $`50`$ digits at $`s=200`$, where $`u(s)10^{820}`$. In the sequel we drop the distinction between $`u(s)`$ and its numerical approximation. The arithmetic computing is done partially with Mathematica$`^\text{®}`$ and for the computationally intensive tasks with the C++-based multiprecision package MPFUN++ .

To solve initial value problems for ordinary differential equations highly sophisticated iteration schemes are available, like Runge-Kutta, Adams-Bashford and multi-step methods. For arbitrary high (but fixed) precision results, all these methods become ineffective, since the step size is a decreasing function of the required precision goal for the solution and tends to become ineffectively small. The only remaining choice is to Taylor expand the solution at a given point. The step size is limited by the radius of convergence only and the precision is controlled by the error made in truncating the Taylor series at some order .

$`u(s)`$ is expanded at $`s_0`$ as

$$u(s)=\underset{n0}{}u_n(ss_0)^n.$$
(4.8)
For the Painlevé II equation the expansion coefficients $`u_n`$ at $`s_0`$ are determined by $`u_0=u(s_0)`$, $`u_1=u^{}(s_0)`$ and

$$u_{n+2}=\frac{2u_n^{(3)}+s_0u_n+u_{n1}}{(n+2)(n+1)},$$
(4.9)
where $`u_n^{(k)}=_{j=0}^nu_{nj}u_j^{(k1)}`$ are the expansion coefficients of $`u(s)^k`$ at $`s_0`$, $`u_n^{(1)}=u_n`$. We include the factorial into the expansion coefficients instead of taking the bare Taylor coefficients, in order to reduce the workload from multiplications by binomials when multiplying two expansions numerically.

Numerically we find that the Hastings-McLeod solution does not have any pole in a strip $`|\text{Im}(s)|<2.9`$. To have a safety margin we choose a step size one for the extrapolation of the expansion (4.8).

We take the starting values $`u(s_0)`$, $`u^{}(s_0)`$ from (4.2) at $`s_0=100`$. The coefficients of the functions $`U(s)`$, $`V(s)`$, see (3.33) and (3.38), when expanded around $`s_0`$, as in (4.8), are given by

$$U_{n+1}=\frac{u_n}{n+1},V_{n+2}=\frac{u_n^{(2)}}{(n+2)(n+1)},n0,$$
(4.10)
and $`V_1=u_0^4u_1^2+s_0u_0^2`$, leaving unspecified the yet unknown integration constants $`U_0=U(s_0)`$ and $`V_0=V(s_0)`$. The higher expansion coefficients of $`u`$, $`u^{}`$, $`U`$, and $`V`$ are independent of the values for $`U_0`$ and $`V_0`$ and are calculated with the recursion relations (4.9) and (4.10). The values of $`u`$, $`u^{}`$, $`U`$, and $`V`$ at $`s=s_0\pm 1`$ are extrapolated, the expansion coefficients at $`s_0\pm 1`$ iterated. Then values are calculated at $`s=s_0\pm 2`$ by extrapolation, and so on. A posteriori we assign to $`U(s_0)`$ and $`V(s_0)`$ values, such that $`U(200)=0=V(200)`$. The numerical errors from iterating (4.9) and from truncating (4.8) can be neglected compared to the uncertainty originating from the initial conditions. At the end of this first step we have at our disposal the expansion coefficients for $`u,u^{},U,V`$ at the integers in the interval $`[20,200]`$. For the convenience of the interested reader let us just state the values at $`s=0`$ up to $`50`$ digits,

$`u(0)`$ $`=`$ $`\text{-0.367061551548078427747792113175610961512192053613139},`$
$`u^{}(0)`$ $`=`$ $`\text{0.295372105447550054557007047310237988227233798735629},`$
$`U(0)`$ $`=`$ $`\text{0.336960697930551393597884426960964843885993886628226},`$
$`V(0)`$ $`=`$ $`\text{0.0311059853063123536659591008775670005642241689547838},`$
which might be used as starting values for a quick conventional integration of Painlevé II to reproduce parts of our results with much less effort but also less precision. Tables can be found at .

The next step is to determine $`a(s,y)`$, $`b(s,y)`$ at $`s_0\{20,\mathrm{},200\}`$ in the interval $`y[9,9]`$ employing (3) and (3.38). Setting

$`a(s,y)={\displaystyle \underset{m,n0}{}}a_{m,n}(ss_0)^m(yy_0)^n,`$
$`b(s,y)={\displaystyle \underset{m,n0}{}}b_{m,n}(ss_0)^m(yy_0)^n,`$ (4.11)
(3) becomes a recursion relation for the expansion coefficients,

$`a_{m,n+1}`$ $`=`$ $`{\displaystyle \frac{1}{n+1}}{\displaystyle \underset{k=0}{\overset{m}{}}}\left(u_k^{(2)}a_{mk,n}(k+1)u_{k+1}b_{mk,n}u_kb_{mk,n1}\right),`$
$`b_{m,n+1}`$ $`=`$ $`{\displaystyle \frac{1}{n+1}}(b_{m,n2}b_{m1,n}`$
$`+{\displaystyle \underset{k=0}{\overset{m}{}}}(u_k^{(2)}b_{mk,n}+(k+1)u_{k+1}a_{mk,n}u_ka_{mk,n1})),`$
$`n0`$, allowing one to determine $`a_{0,n}`$, $`b_{0,n}`$, $`n0`$ upon the knowledge of $`a_{0,0}`$, $`b_{0,0}`$. We integrate along $`\pm y`$ with an extrapolation step size of $`\frac{1}{8}`$. From (3) one obtains the recursions

$`a_{m+1,n}`$ $`=`$ $`{\displaystyle \frac{1}{m+1}}{\displaystyle \underset{k=0}{\overset{m}{}}}u_kb_{mk,n}`$
$`b_{m+1,n}`$ $`=`$ $`{\displaystyle \frac{1}{m+1}}\left(b_{m,n1}+{\displaystyle \underset{k=0}{\overset{m}{}}}u_ka_{mk,n}\right).`$ (4.13)
The expansion coefficients $`G_{m,n}`$ of $`G(s,y)`$ at $`(s_0,y_0)`$, are determined from (3.39) as

$$G_{m,n}=(n+1)\left(a_{m,n}^{}a_{m,n+1}b_{m,n}^{}b_{m,n+1}\right)$$
(4.14)
where $`a_{m,n}^{}`$, $`b_{m,n}^{}`$ are the corresponding expansion coefficients of $`a`$ and $`b`$ at $`(s_0,y_0)`$.

To finally determine $`g(y)`$ and its derivatives we write

$`g^{(n)}(y_0)`$ $`=`$ $`{\displaystyle \frac{d^n}{dy_0^n}}{\displaystyle \underset{s_0}{}}{\displaystyle _{s_0}^{s_0+1}}(sy_0^2)^2{\displaystyle \frac{d^2}{ds^2}}\left(G(s,y_0)F_{\text{GUE}}(s)\right)𝑑s`$ (4.15)
$`=`$ $`{\displaystyle \underset{s_0}{}}n!{\displaystyle \underset{m1}{}}c_{m,n}.`$
$`c_{m,n}`$ are the expansion coefficients of $`(s,y)_{s_0}^s(ry^2)^2\frac{d^2}{dr^2}\left(G(r,y)F_{\text{GUE}}(r)\right)𝑑r`$ at $`(s_0,y_0)`$,

$$c_{m,n}=(m2)(GF)_{m1,n}2m(GF)_{m,n2}+\frac{m(m+1)}{m1}(GF)_{m+1,n4}.$$
(4.16)
Here

$$(GF)_{m,n}=\underset{k=0}{\overset{m}{}}F_kG_{mk,n}$$
(4.17)
are the expansion coefficients of $`G(s,y)F_{\text{GUE}}(s)`$ and $`F_n=_{k=1}^n\frac{k}{n}V_kF_{nk}`$ are the expansion coefficients of $`F_{\text{GUE}}`$. Numerically the sum over $`s_0`$ in (4.15) is truncated to values inside $`[15,200]`$, since outside contributions turn out to be negligible at the chosen precision goal. After accomplishing this program we keep values for $`g(y)`$ at $`y\frac{1}{128}[9,9]`$ and for $`g^{(n)}(y)`$, $`n=0,\mathrm{},4`$, at $`y\frac{1}{8}[9,9]`$ with an accuracy of about $`100`$ digits (a table in ASCII format is available online at ). For interpolating these values we used the Interpolation-function of the Mathematica$`^\text{®}`$ package yielding best results due to the high precision data with an interpolation order of $`57`$.

## 5 Discussion of the scaling function

There have been numerous attempts to approximately determine $`g()`$ . For historical reasons a different scaling function, $`F()`$, is analyzed in some of these works. The relation to our scaling function $`g()`$ is

$`F(\xi )=(\xi /2)^{2/3}g\left((2\xi ^2)^{1/3}\right),\text{resp.}g(y)=2yF\left(1/(2^{1/2}y^{3/2})\right).`$ (5.18)
Note that by (1.6) the large $`y`$ behavior of $`g`$ is fixed by definition as $`g(y)2|y|`$. The special value $`g(0)=1.1503944782594709729961`$ is the Baik-Rains constant . In the literature the universal amplitude ratio $`R_G=2^{2/3}g(0)=0.7247031092`$ and the universal coupling constant $`g^{}=g(0)^{3/2}=0.810456700`$ have been investigated. Approximate values have been determined by means of Monte-Carlo simulations for the single step model , numerically within a mode-coupling approximation , and even experimentally for slowly combusting paper  yielding estimates for $`g(0)`$ within reasonable ranges around the (numerically) exact value indicated.

In Figure 2 the scaling function $`f(y)=\frac{1}{4}g^{\prime \prime }(y)`$ is

shown as determined by the multiprecision expansion method explained in the previous section. We estimate its large $`y`$ asymptotics as

$$\mathrm{log}f(y)c|y|^3+o(|y|)\text{ for }y\mathrm{}.$$
(5.19)
The cubic behavior is very robust and numerical fits yield about $`2.996`$$`2.998`$ quite independently of the assumed nature of the finite size corrections. The prefactor $`c=0.295(5)`$ has a relatively high uncertainty because of the unknown subleading corrections. Even though unaccessible in nature we estimate the error term, as indicated in (5.19), to be sublinear or even only logarithmic from the numerical data. Possibly, the exact asymptotic behavior could be extracted from a refined asymptotic analysis of the Riemann-Hilbert problem.

Colaiori and Moore  tackled the same scaling function by completely different means. Starting from the continuum version of the KPZ equation they numerically solved the corresponding mode-coupling equation , which contains an uncontrolled approximation, since diagrams which would renormalize the three-point vertex coupling are neglected. Nevertheless a qualitative comparison of their result with the exact scaling function $`f(y)`$ shows reasonable similarity, cf. Figure 3. Both functions are normalized to integral $`1`$ by definition. The mode

coupling solution oscillates around $`0`$ for $`|y|>3`$, whereas $`f(y)>0`$ for the exact solution. We do not know whether this is a numerical artifact or an inherent property of the mode-coupling equation. On the other hand, the second moments are reasonably close together, $`0.510523`$ for $`f(y)`$, and $`0.4638`$ for the mode-coupling approximation. So is the value of the Baik-Rains constant $`g(0)=2|y|f(y)𝑑y`$ for which mode-coupling predicts the value $`1.1137`$.

From the solution to the mode-coupling equations one does not directly obtain $`f(y)`$, but rather its Fourier transform. The function $`G(\tau )`$ from  is defined through

$$G(k^{3/2}/2^{7/2})=\widehat{f}(k)=2_0^{\mathrm{}}\mathrm{cos}(ky)f(y)𝑑y.$$
(5.20)
Moore and Colaiori predict a stretched exponential decay of $`G(\tau )`$ as $`\mathrm{exp}(c|\tau |^{2/3})`$  and numerically find a superimposed oscillatory behavior on the scale $`|\tau |^{2/3}`$ . In Figure 4

$`\widehat{f}(k)`$ is plotted as obtained by a numerical Fourier transform of $`f(y)`$. Indeed it exhibits an oscillatory behavior as can be seen in Figure 5 where the modulus of $`\widehat{f}(k)`$ is shown on a semilogarithmic scale. The dotted line in the plot is the modulus of the function

$$10.9k^{9/4}\mathrm{sin}(\frac{1}{2}k^{3/2}1.937)e^{\frac{1}{2}k^{3/2}},$$
(5.21)
shifted by a factor of $`1000`$ for visibility, which fits $`\widehat{f}(k)`$ very well in phase and amplitude for $`k15`$. This behavior is not in accordance with the results of Colaiori and Moore, since the oscillations and the exponential decay of $`G(\tau )`$ for the exact solution are apparently on the scale $`\tau `$ and not $`\tau ^{2/3}`$.

Note that $`\widehat{f}(k)`$ is the scaling function for the intermediate structure function

$$S(k,t)=𝑑xe^{ikx}S(x,t)2\widehat{f}(t^{2/3}k).$$
(5.22)
By Fourier transforming with respect to $`t`$ we determine the dynamical structure function,

$$S(k,\omega )=𝑑x𝑑te^{i(kx+\omega t)}S(x,t)2k^{3/2}\stackrel{}{f}(\omega /k^{3/2}),$$
(5.23)
where

$`\stackrel{}{f}(\tau )={\displaystyle 𝑑se^{i\tau s}\widehat{f}(s^{2/3})}=2{\displaystyle _0^{\mathrm{}}}𝑑y\tau ^1L^{}(y/\tau ^{2/3})f(y)`$ (5.24)
and $`L`$ has the convenient representation

$$L(\kappa )=23^{2/3}\text{Ai}(3^{4/3}\kappa ^2)\mathrm{sin}(2\kappa ^3/27).$$
(5.25)
The correlation function (2.6) in Fourier space is given by

$$C(k,\omega )=2k^2S(k,\omega )C^{\text{KPZ}}(k,\omega )\stackrel{\text{def}}{=}4k^{7/2}\stackrel{}{f}(\omega /k^{3/2}),$$
(5.26)
describing the asymptotic behavior at $`k,\omega =0`$. Note that $`C(k,\omega )>0`$ by definition, since $`h_{k,\omega }h_{k^{},\omega ^{}}=\delta _{k,k^{}}\delta _{\omega ,\omega ^{}}C(k,\omega )`$ for $`(k,\omega )(0,0)`$. The anomalous scaling behavior in real space is reflected by the exponents for the divergence of $`C^{\text{KPZ}}(k,\omega )`$ at $`k=\omega =0`$. In the linear case, the Edwards-Wilkinson equation $`\lambda =0`$ in (1.1), one easily obtains

$$C^{\text{EW}}(k,\omega )=\frac{D}{\omega ^2+\nu ^2k^4}.$$
(5.27)
A 3d-plot of $`C^{\text{KPZ}}(k,\omega )`$ is shown in Figure 6.

Its striking features are the smooth behavior away from $`k,\omega =0`$, especially on the lines where $`k=0`$ and $`\omega =0`$ and the two symmetric maxima of $`kC^{\text{KPZ}}(k,\omega )`$ for constant $`\omega `$. Our numerical data yield for the singular behavior at $`k=0`$, $`\omega =0`$,

$`C^{\text{KPZ}}(k,\omega )`$ $`=`$ $`\omega ^{7/3}\left(2.10565(1)+0.85(1)k^2w^{4/3}+𝒪(k^4\omega ^{8/3})\right),`$ (5.28)
$`=`$ $`k^{7/2}\left(19.4443(1)52.5281(1)\omega ^2k^3+𝒪(\omega ^4k^6)\right).`$
## 6 Conclusions and Outlook

For systems close to equilibrium many properties valid in generality rely on detailed balance, amongst them in particular the link between correlation and response functions. The KPZ equation does not satisfy detailed balance, since the growth is directed. However, it has been speculated that in $`1+1`$ dimensions detailed balance is recovered in the scaling regime. With our exact scaling function at hand, such a claim can be tested.

Detailed balance implies that the eigenvalues of the generator in the master equation lie on the negative real axis. Thus autocorrelations in the form $`X(t)X(0)`$ can be written as the Laplace transform of a positive measure. The structure function $`S(k,t)`$ at fixed $`k`$ is such an autocorrelation. Using the scaling form (5.22) detailed balance would imply

$$S(k,t)=_{\mathrm{}}^0\nu \left(|k|^{3/2}d\lambda \right)e^{\lambda |t|}$$
(6.1)
with $`\nu (d\lambda )0`$. In particular $`S(k,t)0`$. From (5.21) we know that $`S(k,t)`$ oscillates around zero. Definitely, at $`|k|5`$ there is a negative dip, cf. Fig. 4. Thus (6.1) cannot be correct.

The Bethe ansatz  indicates that, for large system size, the density of states is concentrated on an arc touching $`0`$. If so, the integration in (6.1) would have to be replaced by a corresponding line integral. It is not clear to us how to extract from the numerical knowledge of $`S(k,t)`$ such a representation.

Our main result is the exact scaling function $`f`$, see Figure 3, for the two-point function of the stationary KPZ equation in $`1+1`$ dimensions. “Exact” must be qualified in two respects. Firstly $`f`$ is given indirectly through the solution of certain differential equations, which can be solved numerically only with considerable effort. The errors are well controlled, however. Secondly, we rely on universality, in the sense that the scaling function is derived for the PNG model, which is one rather particular model within the KPZ universality class. Of course, it would be most welcome to establish the scaling limit also for other models in this class.

The KPZ equation (1.1) is a two-dimensional field theory and, in spirit, belongs to the same family as two-dimensional models of equilibrium statistical mechanics, one-dimensional quantum spin chains, and other $`(1+1)`$-dimensional quantum field theories at zero temperature. While in the latter cases, there are a number of models for which the two-point function can be computed, in the dynamical context such solutions are scarce. In addition, the KPZ equation does not satisfy the condition of detailed balance. Such nonreversible models are known to be difficult and we believe that the PNG model is the first one in the list of exact solutions, disregarding noninteracting field theories.

For the nonstationary KPZ equation with a macroscopic profile of nonzero curvature the analogue of $`F_0`$ is the Tracy-Widom distribution function. In that case the full statistics of $`xh(x,t)`$ for large but fixed $`t`$ is available . It is conceivable that an extension of the techniques used there also admits a more detailed study of, say, the joint distribution of $`h(x,t)h(0,0)`$, $`h(x^{},t)h(0,0)`$. On the other hand the joint distribution of $`h(0,t)h(0,0)`$, $`h(0,t^{})h(0,0)`$ does not seem to be accessible. In the representation through the directed polymers it means that space-like points, even several of them, can be handled, whereas time-like points remain a challenge.
Acknowledgments. We greatly enjoyed the collaboration with Jinho Baik at the early stage of this project and are most grateful for his important input. We also thank Mike Moore for providing us with the numerical solution of the mode-coupling equation.