text:
# Nonlinear Beltrami equation and 𝜏-function for dispersionless hierarchies

## 1 Introduction

Dispersionless integrable hierarchies attracted a considerable interest during the last ten years (see e.g. -). Recently it became clear that they play an important role in various problems of hydrodynamics and complex analysis -.

Dispersionless integrable hierarchies can be described in different forms within different approaches. In the papers  it was shown that such hierarchies can be introduced starting with the nonlinear Beltrami equation (quasi-classical $`\overline{}`$-problem)

$`S_{\overline{z}}=W(z,\overline{z},S_z),`$ (1)
where $`z`$, bar means complex conjugation, $`S_z=\frac{S(z,\overline{z})}{z},`$ $`S_{\overline{z}}=\frac{S(z,\overline{z})}{\overline{z}},`$ and $`W`$ (quasi-classical $`\overline{}`$-data) is an analytic function of $`S_z`$. Applying the quasi-classical $`\overline{}`$-dressing method based on equation (1), one can get dispersionless integrable hierarchies and the corresponding addition formulae in a very regular and simple way. Such an approach reveals also the connection of dispersionless hierarchies with the quasi-conformal mappings on the plane.

In the present paper we demonstrate that the quasi-classical $`\overline{}`$-dressing method based on equation (1) leads to explicit formula for the $`\tau `$-function for dispersionless hierarchies, which is connected with the Lagrangian for equation (1) and the corresponding action, evaluated on the solution of the boundary problem for this equation. We will also derive determinant form of the generating equations for the function $`S_z`$, defining infinitesimal deformations of the $`\tau `$-function. In this paper we concentrate on the dispersionless Kadomtsev-Petviashvili (dKP) hierarchy, but we present also the basic formulae for the dispersionless 2DTL hierarchy.

## 2 $`\tau `$-function as an action for nonlinear Beltrami equation

It was shown in  that the dKP hierarchy is connected with Beltrami equation (1) with the $`\overline{}`$-data equal to zero outside the unit disc. This problem can be formulated as a boundary problem for equation (1) in the unit disc as follows. Let the function $`S_0(z)`$ analytic in the unit disc $`D`$ be given. The problem is to find the function $`S=S_0+\stackrel{~}{S}`$, satisfying (1), with $`\stackrel{~}{S}`$ analytic outside the unit disc and decreasing at infinity (this is in fact a boundary condition on the unit circle, which can be written down using standard projection operator). We suggest that the function $`W`$ is of the form

$`W(z,\overline{z},S_z)={\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}w_p(z,\overline{z})(S_z)^p,`$ (2)
where $`w_p(z,\overline{z})`$ are arbitrary smooth functions in the unit disc vanishing on the boundary.

Introducing parameterization of the function $`S_0(z)`$ in terms of times, $`S_0(z)=_{n=1}^{\mathrm{}}t_nz^n`$, and using the technique of quasi-classical $`\overline{}`$-dressing method, it is possible to demonstrate that $`S(z,𝐭)`$ is a solution of Hamilton-Jacobi equations for dKP hierarchy, and the first coefficient of expansion of $`\stackrel{~}{S}(z,𝐭)`$ as $`z\mathrm{}`$ satisfies equations of dKP hierarchy (see ).

Here we establish a relation between the action for the problem (1) and the $`\tau `$-function for dKP hierarchy. This relation illustrates a well-known observation that a transition from dispersionfull to dispersionless hierarchies resembles a transition from quantum mechanics to classical mechanics.

It was noted in  that equation (1) is a Lagrangian one. It can be obtained by variation of the action (for the boundary problem in the unit disc)

$`f={\displaystyle \frac{1}{2\pi \text{i}}}{\displaystyle _D}\left({\displaystyle \frac{1}{2}}\stackrel{~}{S}_{\overline{z}}\stackrel{~}{S}_zW_1(z,\overline{z},S_z)\right)𝑑zd\overline{z},`$ (3)
where

$`W_1(z,\overline{z},S_z)={\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}w_p(z,\overline{z}){\displaystyle \frac{(S_z)^{p+1}}{p+1}},_\eta W_1(z,\overline{z},\eta )=W(z,\overline{z},\eta ).`$
One should consider independent variations of $`\stackrel{~}{S}`$, possessing required analytic properties (analytic outside the unit circle, decreasing at infinity), keeping $`S_0`$ fixed.

###### Proposition 1

The function

$`F(𝐭)={\displaystyle \frac{1}{2\pi \text{i}}}{\displaystyle _D}\left({\displaystyle \frac{1}{2}}\stackrel{~}{S}_{\overline{z}}(𝐭)\stackrel{~}{S}_z(𝐭)W_1(z,\overline{z},S_z(𝐭))\right)𝑑zd\overline{z},`$ (4)
i.e., the action (3) evaluated on the solution of the problem (1), is a $`\tau `$-function of dKP hierarchy.

Proof In order to prove that $`F(𝐭)`$ is a $`\tau `$-function of dKP hierarchy, it is sufficient to demonstrate that (see, e.g., )

$$\stackrel{~}{S}(z,𝐭)=D(z)F(𝐭),$$

where $`D(z)`$ is the quasiclassical vertex operator, $`D(z)=_{n=1}^{\mathrm{}}\frac{1}{n}\frac{1}{z^n}\frac{}{t_n}`$, $`|z|>1`$. Applying the operator $`D(z)`$ to the r.h.s. of formula (4), one gets (we change the variable of integration to $`y`$)

$`D(z)W_1(y,\overline{y},S_y(𝐭))=\stackrel{~}{S}_{\overline{y}}D(z)(\stackrel{~}{S}_y+S_{0}^{}{}_{y}{}^{}).`$
Using the formula

$$D(z)S_{0}^{}{}_{y}{}^{}=\frac{1}{zy},$$

we get

$`{\displaystyle \frac{1}{2\pi \text{i}}}{\displaystyle _D}\stackrel{~}{S}_{\overline{y}}D(z)S_{0}^{}{}_{y}{}^{}𝑑yd\overline{y}={\displaystyle \frac{1}{2\pi \text{i}}}{\displaystyle \frac{1}{zy}\stackrel{~}{S}(y)𝑑y}=\stackrel{~}{S}(z).`$
Thus we have obtained the crucial term for our proof, and now we should demonstrate that the combination of all other terms, namely,

$`{\displaystyle \frac{1}{4\pi \text{i}}}{\displaystyle _D}(\stackrel{~}{S}_yD(z)\stackrel{~}{S}_{\overline{y}}\stackrel{~}{S}_{\overline{y}}D(z)\stackrel{~}{S}_y)𝑑yd\overline{y}`$
is equal to zero. Indeed, using Green’s formula and taking into account that $`\stackrel{~}{S}_{\overline{y}}`$ is equal to zero on the unit circle, and the function $`\stackrel{~}{S}(y)`$ is analytic outside the unit circle and decreases at infinity, we immediately come to the conclusion that this combination is equal to zero. QED
Remark 1. Note that integral formulae for the $`\tau `$-function of different type has been derived also in .
Remark 2. Let us consider also the integral

$`F_s(𝐭)={\displaystyle \frac{1}{2\pi \text{i}}}{\displaystyle _D}\left({\displaystyle \frac{1}{2}}S_{\overline{z}}(𝐭)S_z(𝐭)W_1(z,\overline{z},S_z(𝐭))\right)𝑑zd\overline{z}.`$ (5)
Since $`\overline{}S_0=0`$ for $`zD`$, using Green’s formula, one obtains

$$F_s=F+\frac{1}{4\pi \text{i}}𝑑z\stackrel{~}{S}(z)S_{0}^{}{}_{z}{}^{}.$$

Taking into account the relation $`\stackrel{~}{S}(z)=D(z)F`$, one gets

$`F_s=F{\displaystyle \frac{1}{2}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}t_n{\displaystyle \frac{F}{t_n}}.`$ (6)
The dKP hierarchy and addition formula for $`F`$ admits scale invariance

$$F(𝐭)F^{}(\lambda 𝐭)=\lambda ^2F(𝐭).$$

A full infinitesimal variation of $`F`$ under this transformation is

$$\delta _sF=\delta _{\text{form}}F+\delta _\text{t}F=2ϵF,$$

where $`\delta _{\text{form}}F`$ denotes a variation of the form of $`F`$, while $`\delta _\text{t}F`$ stands for the variation due to the infinitesimal variation of times $`𝐭`$. In virtue of (6), one has $`\delta _{\text{form}}F=2ϵF_s.`$ In the particular case $`W(z,\overline{z},S_z)=\mu (z,\overline{z})S_z`$ we have $`F_s=0`$ and thus $`F`$ is a homogeneous function of times of the second order. Such a class of $`\tau `$-functions has been considered within different approaches in . Starting with linear Beltrami equation, we obtain only a subclass of $`\tau `$-functions, which are quadraic with respect to times.

It is well known that dKP $`\tau `$-function $`F`$ obeys the dispersionless addition formula

$`(z_1z_2)e^{D(z_1)D(z_2)F}+(z_2z_3)e^{D(z_2)D(z_3)F}+(z_3z_1)e^{D(z_3)D(z_1)F}=0,`$ (7)
$`z_1,z_2,z_3D.`$
The formula (4) gives a solution to this equation in terms of solution of nonlinear Beltrami equation (1).

## 3 Variations of the $`\tau `$-function

The function $`W`$ is the $`\overline{}`$ data for the dKP hierarchy. Its variations provide us with a wide class of variations of the function $`F`$. For the functions $`W`$ of the form (2), varying $`w_n(z,\overline{z})`$, one has

$`\delta W={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\delta w_n(z,\overline{z})\left(S_z\right)^n,\delta W_1={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\delta w_n}{n+1}}\left(S_z\right)^{n+1},`$
and

$`\delta F={\displaystyle \frac{1}{2\pi \text{i}}}{\displaystyle _D}(\delta W_1)(z,\overline{z},S_z)𝑑zd\overline{z}.`$ (8)
Considering elementary variation $`\delta w_{n_0}=ϵ\alpha _{n_0}\delta (zz_0)`$, $`\delta w_n=0,nn_0`$, one gets

$`\delta F={\displaystyle \frac{ϵ}{2\pi \text{i}}}{\displaystyle \frac{\alpha _{n_0}}{(n_0+1)}}\left(S_z\right)^{n_0+1}|_{z=z_0},`$ (9)
and, respectively,

$`\delta \stackrel{~}{S}={\displaystyle \frac{ϵ}{2\pi \text{i}}}{\displaystyle \frac{\alpha _{n_0}}{(n_0+1)}}D(z)\left(S_z(z_0)\right)^{n_0+1}.`$ (10)
Taking superposition of elementary variations (9), we obtain a general variation of the form

$`\delta F=ϵf(S_z(z_0)),`$ (11)
where $`f`$ is an arbitrary analytic function (summation over different points and integration over $`z_0`$ are also possible).

The formulae (10), (9) can be also derived considering the deformations of nonlinear Beltrami equation (1).
Remark 3. Since a variation of the $`\overline{}`$-data $`W`$ transforms solution of the dKP hierarchy into another solution, then the formula (11) defines an infinitesimal symmetry transformation for the function $`F`$ satisfying equation (7). It is possible to prove this statement directly starting with the formula

$$p(z_0)p(z)+z\mathrm{exp}(D(z)S(z_0))=0,z_0,zD.$$
(12)
Derivation of this formula can be found in, e.g., .

### 3.1 Determinant form of equation for $`\varphi `$

Existence of symmetry transformation of the form (11) leads us directly to equation for the function $`\varphi =S_z(z_0)`$. Indeed, let us consider a special symmetry transformation (11) of the form

$$F^{}=F+ϵ\mathrm{exp}(\mathrm{\Theta }\varphi ),$$

where $`\mathrm{\Theta }`$ is an arbitrary parameter, and substitute it to (7). Then we get a system of linear equations

$`\{\begin{array}{c}x+y+z=0,\hfill \\ (D_2D_3\varphi )x+(D_1D_3\varphi )y+(D_1D_2\varphi )z=0,\hfill \\ (D_2\varphi )(D_3\varphi )x+(D_1\varphi )(D_3\varphi )y+(D_1\varphi )(D_2\varphi )z=0,\hfill \end{array}`$
where we use notations $`D_i=D(z_i)`$ and

$$x=(z_2z_3)e^{D_2D_3F},y=(z_3z_1)e^{D_3D_1F},z=(z_1z_2)e^{D_1D_2F}.$$

The condition that determinant of this system is equal to zero gives the equation for the function $`\varphi `$, ($`\varphi _i=D_i\varphi `$)

$`\text{det}\left(\begin{array}{ccc}1& 1& 1\\ \varphi _2\varphi _3& \varphi _1\varphi _3& \varphi _1\varphi _2\\ \varphi _{23}& \varphi _{13}& \varphi _{12}\end{array}\right)=0.`$ (17)
Expanding the l.h.s. of this equation into powers of parameters $`z_1^1`$, $`z_2^1`$, $`z_3^1`$, in the order $`z_1^1z_2^2z_3^3`$ one gets the equation

$`_x\left({\displaystyle \frac{\varphi _t}{\varphi _x}}{\displaystyle \frac{3}{8}}\left({\displaystyle \frac{\varphi _y}{\varphi _x}}\right)^2\right)={\displaystyle \frac{3}{4}}_y\left({\displaystyle \frac{\varphi _y}{\varphi _x}}\right),`$ (18)
which is the KP singular manifold equation in dispersionless limit.

It is interesting to note that equation (17) written in the form

$`D_1\mathrm{log}\left({\displaystyle \frac{D_2\varphi }{D_3\varphi }}\right)+D_2\mathrm{log}\left({\displaystyle \frac{D_3\varphi }{D_1\varphi }}\right)+D_3\mathrm{log}\left({\displaystyle \frac{D_2\varphi }{D_3\varphi }}\right)=0`$ (19)
has been derived in  as a naive continuous limit of the discrete Möbius-invariant KP equation  (see equation (33)) in connection with the Menelaus theorem. It is easy to check that equation (17) (or, equivalently, (19)) is invariant under conformal transformation of dependent variable $`\varphi f(\varphi )`$, where $`f`$ is an analytic function.

The determinant form (17) gives some hint for the geometric interpretation of this equation. Indeed, as it is known (see, e.g., ), the area of the plane triangle with the coordinates of vertices given by the pairs $`(x_1,y_1)`$, $`(x_2,y_2)`$, $`(x_3,y_3)`$, can be written as

$`A={\displaystyle \frac{1}{2}}\left|\text{det}𝐀\right|,𝐀=\left(\begin{array}{ccc}1& 1& 1\\ x_1& x_2& x_3\\ y_1& y_2& y_3\end{array}\right).`$
Thus equation (17) means that area (may be complex) of corresponding triangle vanish, that is, obviously, the only possibility for the geometry admitting arbitrary transformation $`\varphi f(\varphi )`$.

More general generating equation for the gauge-invariant function $`𝒮=S(z_1)S(z_0)`$

$`{\displaystyle ϵ_{ijk}D_j\mathrm{log}(\mathrm{exp}(D_i𝒮)1)}=0`$ (21)
derived in  (from which equation (19) can be easily obtained by the limit $`z_1z_0`$) can be also written in the determinant form,

$`\text{det}\left(\begin{array}{ccc}1& 1& 1\\ (e^{𝒮_2}1)(e^{𝒮_3}1)& (e^{𝒮_3}1)(e^{𝒮_1}1)& (e^{𝒮_1}1)(e^{𝒮_2}1)\\ 𝒮_{23}& 𝒮_{13}& 𝒮_{12}\end{array}\right)=0,`$ (25)
where the subscript $`i`$ refers to the derivative $`D_i`$. We believe it is possible to obtain this equation in a manner similar to the derivation of (17).

### 3.2 Determinant form of discrete SKP equation

The basic idea of derivation of the determinant form (17) is applicable to the dispersionfull case too (in fact, originally we used it in the dispersionfull case first). This idea gives another way to obtain discrete Swartzian KP equation, which was introduced in . Instead of equation (7), we start from the well-known addition formula for the KP $`\tau `$-function

$`c_1(T_1\tau )(T_2T_3\tau )+c_2(T_2\tau )(T_1T_3\tau )+c_3(T_3\tau )(T_1T_2\tau )=0,`$ (26)
where $`T_i`$ denotes a Sato shift, $`c_i`$ are certain coefficients.

###### Proposition 2

If for some function $`\mathrm{\Phi }`$ the function

$`\stackrel{~}{\tau }=\tau (1+\mathrm{\Theta }\mathrm{\Phi }),`$ (27)
satisfies equation (26) for arbitrary $`\mathrm{\Theta }`$ (i.e., formula (27) defines a Bäcklund transformation for the $`\tau `$-function), then the function $`\mathrm{\Phi }`$ satisfies the equation

$`det\left(\begin{array}{ccc}1& 1& 1\\ (T_1\mathrm{\Phi }+T_2T_3\mathrm{\Phi })& (T_2\mathrm{\Phi }+T_1T_3\mathrm{\Phi })& (T_3\mathrm{\Phi }+T_1T_2\mathrm{\Phi })\\ (T_1\mathrm{\Phi })(T_2T_3\mathrm{\Phi })& (T_2\mathrm{\Phi })(T_1T_3\mathrm{\Phi })& (T_3\mathrm{\Phi })(T_1T_2\mathrm{\Phi })\end{array}\right)=0.`$ (31)
Proof Substituting (27) into (26), we obtain a system of linear equations

$`\{\begin{array}{c}x+y+z=0,\hfill \\ (T_1\mathrm{\Phi }+T_2T_3\mathrm{\Phi })x+(T_2\mathrm{\Phi }+T_1T_3\mathrm{\Phi })y+(T_3\mathrm{\Phi }+T_1T_2\mathrm{\Phi })z=0,\hfill \\ (T_1\mathrm{\Phi })(T_2T_3\mathrm{\Phi })x+(T_2\mathrm{\Phi })(T_1T_3\mathrm{\Phi })y+(T_3\mathrm{\Phi })(T_1T_2\mathrm{\Phi })z=0,\hfill \end{array}`$
where

$$x=c_1(T_1\tau )(T_2T_3\tau ),y=c_2(T_2\tau )(T_1T_3\tau ),z=c_3(T_3\tau )(T_1T_2\tau ).$$

Then, from the requirement that determinant of this linear system should be equal to zero, we get equation (31), QED.

It is easy to check that equation (31) coincides with discrete SKP equation, introduced in  in multiplicative form

$`(T_2\mathrm{\Delta }_1\mathrm{\Phi })(T_3\mathrm{\Delta }_2\mathrm{\Phi })(T_1\mathrm{\Delta }_3\mathrm{\Phi })=(T_2\mathrm{\Delta }_3\mathrm{\Phi })(T_3\mathrm{\Delta }_1\mathrm{\Phi })(T_1\mathrm{\Delta }_2\mathrm{\Phi }),`$ (33)
where $`\mathrm{\Delta }_i=T_i1`$ is a difference operator.

## 4 Dispersionless 2DTL hierarchy

Most of the results presented above for dispesionless KP hierarchy are valid after some modification for dispersionless 2DTL hierarchy. First we will outline the basic notations, following the work , and then we will discuss how the main formulae are modified.

For the dispersionless 2DTL hierarchy the $`\overline{}`$-data are localized on the domain $`G`$ which is an annulus $`a<|z|<b`$, where $`a,b`$ ($`a,b`$, $`a,b>0`$; $`b>a`$) are arbitrary (instead of the unit disc in KP case). To set the quasi-classical $`\overline{}`$-problem (1) correctly, in general we do not need to require analyticity of the function $`S_0`$ in $`G`$, it is enough to have analyticity of its derivative $`S_{0}^{}{}_{z}{}^{}`$. A generic function $`S_0`$ with $`S_{0}^{}{}_{z}{}^{}`$ analytic in $`G`$ can be represented as

$$S_0(t,𝐱,𝐲)=t\mathrm{log}z+\underset{n=1}{\overset{\mathrm{}}{}}z^nx_n+\underset{n=1}{\overset{\mathrm{}}{}}z^ny_n,$$

where $`t,x_n,y_n`$ are free parameters . We assume that $`\stackrel{~}{S}(z)_{n=1}^{\mathrm{}}\frac{S_n}{z^n}`$ as $`z\mathrm{}`$ and denote $`\stackrel{~}{S}(0)=\varphi `$, $`G_+=\{z,|z|>b\}`$, $`G_{}=\{z,|z|<a\}`$. Relations, characterizing the $`\tau `$-function $`F`$, are $`\varphi =DF`$, $`\stackrel{~}{S}(z_1)=D_+(z_1)F`$ $`(z_1G_+)`$, $`\stackrel{~}{S}(z_2)=\varphi D_{}(z_2)F`$ $`(z_2G_{})`$, where $`D_+(z)=_{n=1}^{\mathrm{}}\frac{1}{n}\frac{1}{z^n}\frac{}{x_n}`$, $`D_{}(z)=_{n=1}^{\mathrm{}}\frac{1}{n}z^n\frac{}{y_n}`$, $`D=\frac{}{t}`$.

### 4.1 $`\tau `$-function for dispersionless 2DTL hierarchy

Boundary problem for nonlinear Beltrami equation (1) in this case is formulated on the boundary of the annulus $`G`$, and integration in the formula (4) goes over the annulus.

###### Proposition 3

The function

$`F(t,𝐱,𝐲)={\displaystyle \frac{1}{2\pi \text{i}}}{\displaystyle _G}\left({\displaystyle \frac{1}{2}}\stackrel{~}{S}_{\overline{z}}(t,𝐱,𝐲)\stackrel{~}{S}_z(t,𝐱,𝐲)W_1(z,\overline{z},S_z(t,𝐱,𝐲))\right)𝑑zd\overline{z}`$
is a $`\tau `$-function of dispersionless 2DTL hierarchy.

The proof is completely analogous to the KP case.

The formula for the $`\tau `$-function may be also considered in a more general context, for arbitrary domain $`G`$. In this case it defines $`F`$ as a functional on the space of functions $`S_0(z)`$, having the derivative $`S_{0}^{}{}_{z}{}^{}`$ analytic in $`G`$. The form of corresponding hierarchy depends on parameterization of this space in terms of ‘times’.

### 4.2 Symmetries and singular manifold equations

Variations of the function $`F`$, preserving the hierarchy, have the same form (11) as in KP case.

Dispersionless 2DTL equations can be obtained from KP case (equations (21), (19), (17) ) by the transformations $`D_1D_+`$, $`D_2D_{}+_t`$, $`D_3_t`$ or, equivalently, $`D_1D_{}`$, $`D_2D_+_t`$, $`D_3_t`$, connected by the the symmetry $`D_+D_{}`$, $`D_{}D_+`$, $`_t_t`$.

Equation (21) for $`𝒮`$ takes the form

$`\left(e^{D_+𝒮}(e^{_t𝒮}1)+e^{D_{}𝒮}(e^{_t𝒮}1)\right)D_+D_{}𝒮`$
$`(e^{D_+𝒮}1)(e^{D_{}𝒮}1)_t(D_{}D_++_t)𝒮=0.`$ (34)
In the order $`z_+(z_{})^1`$ of expansion of equation (34) one gets

$`𝒮_{xy}𝒮_x𝒮_y{\displaystyle \frac{𝒮_{tt}}{(e^{𝒮_t}1)(1e^{𝒮_t})}}=0,`$ (35)
that is the dispersionless 2DTL singular manifold equation.

Generating equation for conformally invariant case reads

$`(_t\varphi )(D_+D_{}\varphi )(_t+D_{}D_+)\varphi =(D_+\varphi )(D_{}\varphi )_t(_t+D_{}D_+)\varphi .`$ (36)
The order of expansion $`z_+(z_{})^1`$ gives the conformally invariant dispersionless 2DTL equation

$`{\displaystyle \frac{\varphi _{xy}}{\varphi _x\varphi _y}}={\displaystyle \frac{\varphi _{tt}}{\varphi _t\varphi _t}}.`$ (37)
### Acknowledgments

LVB was supported in part by RFBR grant 01-01-00929 and President of Russia grant 1716-2003; BGK was supported in part by the grant COFIN 2000 ‘Sintesi’.


format: text

text:
# Strings in the 𝜂 Carinae nebula: hypersonic radiative cosmic bullets

## 1 INTRODUCTION

Of the many intriguing questions surrounding the luminous blue variable $`\eta `$ Carinae, the origin of the long thin “strings” found in the outer nebular regions remains particularly vexing. First observed by Meaburn et al. (1996) and further studied by Weis, Duschl, & Chu (Weis et al.1999) the strings exhibit two remarkable properties: very high values of the length-to-width ratios and the presence of the Hubble-type flow, i.e. linear velocity increase from the base of the strings to their tip. The strings also show a surface brightness decrease toward the string head and the true tips possibly may be invisible in optical wavelengths. These features represent the key challenges that must be met by any model that hopes to offer a successful explanation of the strings’ origin and evolution.

Several models were suggested for the strings including jets (García-Segura et al., 1999) and ionization shadows (Soker, 2001), see Redman, Meaburn, & Holloway (Redman et al.2002) for a complete discussion of the suggested models. A particularly promising model has been proposed by (Redman et al.2002) who suggest that a single hypersonic bullet propagating in the ambient medium is generating the strings observed in the $`\eta `$ Carinae nebula. Since the Great Eruption that ejected the main nebula was know to be an impulsive event with most of the mass and momentum directed in the polar directions (Smith et al., 2003), and since the strings seem to be coeval with the overall nebula, it is plausible that fragments or bullets were produced in that same event.

In this work we do not address the issue of bullet origin focusing instead on the nature of their evolution and observable properties. In particular, we investigate the feasibility of the (Redman et al.2002) model and we address two questions: (1) can a radiatively cooled hypersonic bullet produce structures with large length-to-width ratios, and (2) do such systems exhibit the Hubble-type flow behaviour in their downstream wake.

In Section 2 we briefly discuss the numerical code used, the key dimensionless parameters characterizing the system, and the setup of the numerical simulation that was carried out. In Section 3 we discuss the results of our numerical study, and in Section 4 we present the conclusions including possible answers to the questions posed above as well as applicability of our results to other astrophysical systems.

## 2 NUMERICAL STUDY

The simulation presented in this paper was performed with the AstroBEAR code based on the BEARCLAW adaptive mesh refinement package (Mitran, 2001; Berger & LeVeque, 1998). The code allows for arbitrary levels of refinement within computer memory constraints, as well as adaptive numerics and multiphysics. The Euler hydrodynamic equations with cooling source terms were solved using explicit second-order–accurate scheme. The wave propagation method (LeVeque, 1997) was applied with the full nonlinear Riemann solver. An operator split set of ODE’s, describing the source term effects, was solved using the fourth-order–accurate fully implicit Kaps-Rentrop–type scheme (Kaps & Rentrop, 1979) (see also (Press et al., 1997)).

There are three dimensionless parameters that describe the evolution of the system studied. The first two are the density contrast between the bullet and the ambient medium $`\chi _b\rho _b/\rho _a`$ and the Mach number $`M_b`$ of the bullet velocity at time $`t=0.0`$. Those two parameters determine the overall hydrodynamic regime of the system evolution. The third one, the cooling parameter $`\psi _b`$, describes the effects of cooling and is defined as

$$\psi _b=\frac{t_{cool}}{t_{hydro}}.$$
(1)
Here the hydrodynamic timescale $`t_{hydro}`$ is defined as the bullet crushing time (or clump crushing time $`t_{cc}`$ as described by Poludnenko, Frank, & Blackman (Poludnenko et al.2002))

$$t_{hydro}=\left(\chi _b^{1/2}\left(F_{c1}F_{st}\right)^{1/2}\right)\frac{2r_b}{v_b},$$
(2)
where $`r_b`$ is bullet radius, $`v_b`$ is bullet velocity, $`F_{c1}1.3`$ and

$$F_{st}1+\frac{2.16}{1+6.55\chi ^{1/2}}.$$
(3)
The two factors $`F_{c1}`$ and $`F_{st}`$ relate the unperturbed upstream conditions with the internal bullet post-shock ones and are described in (, Poludnenko et al.2002). The interaction of the bullet with the ambient medium launches a shock propagating inside the bullet with the post-shock temperature $`T_{ps}`$ and $`\rho _{ps}`$. The cooling rate $`\mathrm{\Lambda }_{ps}`$ for such a temperature can be found based on the cooling curve given by Dalgarno & McCray (1972). Thus the cooling timescale is estimated as follows

$$t_{cool}=\frac{T_{ps}T_{min}}{(\gamma 1)\rho _{ps}\mathrm{\Lambda }_{ps}}km_H,$$
(4)
where $`T_{min}=100`$ K is the minimum temperature gas can cool down to. The internal bullet post-shock temperature $`T_{ps}`$ and density $`\rho _{ps}`$ can be found in the usual manner using the Rankine-Hugoniot relations. $`T_{ps}`$ and density $`\rho _{ps}`$ are the functions of the unshocked bullet temperature and density respectively as well as of the internal bullet shock Mach number $`M_{is}`$ described by the expression

$$M_{is}=M_b\left(F_{c1}F_{st}\right)^{1/2}.$$
(5)
We carried out simulations covering bullet Mach numbers in the range $`M_b=10200`$. In terms of the cooling parameter we explored the regimes ranging from the quasi-adiabatic one $`\psi _b>1`$ to the strongly cooling one $`\psi _b10^5`$. In this paper we will report only on the case that allows us to best address the issues related to $`\eta `$ Carinae strings. From our parameter space explorations this turns out to be the case with the strongest cooling $`\psi _b=2.810^5`$. We defer broader conclusions gleaned from the parameter space exploration to a subsequent paper.

Here we present a simulation of bullet propagation where the initial bullet radius is $`r_b=3.010^{15}`$ cm. The initial bullet density is $`\rho _b=10^5cm^3`$ in the inner 80% of the bullet radius and in the outer 20% of the radius density is smoothed out through the $`tanh`$type function (cf. expression (1) in (, Poludnenko et al.2002)). This translates into the initial bullet mass $`m_b0.510^5\mathrm{M}_{}`$. The initial bullet Mach number is $`M_b=20`$, or $`v_b235`$ km $`\mathrm{s}^1`$. The computational domain is of size $`1.810^{17}`$ cm $`\times 2.410^{16}`$ cm or expressed in terms of bullet radii $`60r_b\times 8r_b`$. The ambient density is $`\rho _a=1000cm^3`$ and the ambient temperature $`T_a=10^4K`$. It is assumed that initially the bullet is in pressure equilibrium with the ambient material.

The run was carried out in cylindrical symmetry with the x-axis being the symmetry axis. We use an adaptive grid with 4 levels of refinement. This provided a resolution of 128 cells per cloud radius or an equivalent resolution of $`1024\times 7680`$ cells. Although the above resolution is sufficient to place the simulation in the converged regime in the adiabatic case (Klein, McKee, & Colella, Klein et al.1994), in the cooling case one has to be more precise in one’s definition of convergence. Since, as it will be discussed below, the key process in the system is bullet fragmentation via instabilities at the upstream surface, our criterion of convergence was the constancy of the initial fragmentation spectrum. In this sense the above resolution ensured that the instability wavelength, and therefore the number of fragments into which the bullet breaks up, does not change with increasing resolution. On the other hand, it should be mentioned that such resolution might not be sufficient to track the contraction of each fragment due to radiative cooling to its smallest size defined by the equilibrium between cooling and external heating. Finally, outflow boundary conditions were prescribed on all boundaries.

For the simulation described above the cooling parameter was $`\psi _b=2.810^5`$. The hydrodynamic timescale, i.e. the bullet crushing time, was $`t_{hydro}=2.5410^9`$ s $`80.5`$ yrs., the cooling timescale was $`t_{cool}=70.910^3`$ s $`19.7`$ hours. The simulation total run time was $`9.4910^9`$ s $`251.1`$ yrs111Note that the simulation run time is larger than the age of the $`\eta `$ Carinae strings ($`150`$ yrs.). However, since it is the dimensionless numbers that determine the system evolution, the results can be appropriately scaled to $`\eta `$ Carinae conditions. We discuss that in further detail in Section 3.2.. Thus the evolution of the system was in a regime strongly dominated by cooling.

## 3 RESULTS

### 3.1 Overall Morphology and Evolution of the System

When the interaction of an inhomogeneity (bullet or clump) proceeds in the adiabatic regime the dominant process defining the evolution of the system is the lateral clump re-expansion (, Klein et al.1994, Poludnenko et al.2002). The internal shock compresses the clump and once such compression is completed the re-expansion begins into regions of lowest total pressure, i.e. in the lateral direction. Such re-expansion, acting in combination with instabilities at the upstream clump surface, is responsible for the destruction of the clump. When more than one clump is present, lateral re-expansion also drives interclump interactions via merging and formation of larger structures that subsequently alter the global flow (, Poludnenko et al.2002).

The fundamental distinction between the radiatively cooled inhomogeneous systems, including individual bullets, and the adiabatic ones is the minimal role of lateral re-expansion. Instead, the dominant process is the formation of instabilities at the upstream surface of the bullet with a wavelength significantly smaller than the one observed in the adiabatic case. As the bullet begins to drive through the ambient medium, hydrodynamic instabilities (Richtmeyer-Meshkov, Rayleigh-Taylor) produce the initial instability seed. The resulting density variations quickly trigger the onset of thermal instabilities which thereafter become the dominant process responsible for bullet fragmentation. In our simulation such instabilities started developing with a wavelength of about 8 - 10% of the initial bullet radius or $`2.410^{14}3.010^{14}`$ cm $`1620`$ a.u. This initial scale is crucial for subsequent evolution since it determines not only the continuing process of bullet fragmentation but also the structure of the downstream flow. As it was mentioned above, we did not observe any significant changes in the fragmentation spectrum with changing grid resolution, therefore we believe that we observed the true fragmentation spectrum corresponding to the given flow conditions. Of course this conclusion may be altered with fully 3-D simulations. The question of predicting the properties of the initial fragmentation spectrum for the given characteristics of the system is important for the understanding of the physics of radiative hypersonic inhomogeneous systems, however it is outside the scope of the current paper and it will be investigated in the subsequent work.

Since the wavelength of the initial fragmentation spectrum is approximately constant along the bullet radius, the fragments produced by the instability are of different mass with the most massive ones being closest to the symmetry axis and the outermost ones being the lightest. As a result the fragments are “peeled off” from the bullet one by one starting with the outermost ones in radius. Each fragmentation event results in the formation of a distinct “ring-like” feature in the bullet wake. The formation of rings is a consequence of the axisymmetry of these simulations. In 3-D it is likely that the rings would themselves fragment (Klein et al., 2003; Robey et al., 2002).

One noteworthy point is that the evolution of the bullet and its fragments suggests the presence of steady “mass loading” of the downstream flow via the hydrodynamic ablation. Mass loading has been claimed to be an important process in all clumpy hydrodynamic systems (Hartquist et al., 1986; Hartquist & Dyson, 1988). The extent and nature of such mass loading will be investigated in the subsequent work.

Figure 1 shows the computational domain at time $`t=251.1`$ yrs. In Figure 1a the synthetic Schlieren image (gradient of the density logarithm) is shown illustrating the shock and vortex sheet structure in the flow. Figure 1b shows the synthetic observation image of the computational domain. Since our simulation did not track the full ionization dynamics of the flow, the image represents the total radiative energy losses summed along each ray. The 2D distribution of the state vector obtained in the simulation was extended using cylindrical symmetry to a $`2048\times 2048\times 7680`$ cells uniform grid 3D data cube. Thus the synthetic observation image represents the 2D projected distribution of the logarithm of the emissivity $`I`$ integrated in the z-direction according to the formula

$$I_{ij}=\underset{k}{}n_{ijk}^2\mathrm{\Lambda }_{ijk}(T_{ijk}),$$
(6)
where $`i,j,`$ and $`k`$ are the cell indices in the x-, y-, and z-direction respectively, and the cooling rate $`\mathrm{\Lambda }(T)`$ here, as well as in the simulation, was determined based on the cooling curve described by Dalgarno & McCray (1972).

There are several ring-like structures in Figure 1 in the bullet wake. These result from the fragmentation events mentioned above. The most prominent ring occurs at the distance of $`1.3810^{17}`$ cm from the bullet head (ring 1) with the width of $`2.9410^{16}`$ cm $`1965`$ a.u. which is the widest part of the bullet wake. In Figure 1b several other bright rings, resulting from fragmentation events, are visible closer to the bullet head. They gradually increase in radius further downstream with the radius of the largest ring (ring 2), located at the distance of $`5.810^{16}`$ cm from the bullet head, being $`1.8210^{16}`$ cm $`1216`$ a.u. Note that the width of the bullet wake stays practically constant, increasing only slightly, for a distance of $`8.010^{16}`$ cm $`5350`$ a.u. from ring 2 to ring 1. Thus in our simulations the length-to-width ratio of the bullet wake is $`6.19.9`$ and the bullet wake width is $`46.5`$ times larger than the width of String 5 while they share the same length. Note also that the width of the bullet head is only $`1.0510^{15}`$ cm $`70`$ a.u. We followed the simulation beyond the time $`t=251.1`$ yrs, shown in Figures 1 and 2, up to the time $`t=300.8`$ yrs. and we did not observe any appreciable changes in the aforementioned values of the length-to-width ratio.

The second process that determines the structure of the flow in the wake is gas re-expansion into the cavity excavated by the bullet. The highest temperature reached by gas in the system is at the tip of the central bullet fragment and is about $`3.510^5`$ K. Gas passing through the bullet bowshock cools down very rapidly, weakening the bowshock and making it more oblique, which in its turn prevents further heating of the ambient gas. As it can be seen in Figure 1a, the bowshock practically disappears half-way downstream from the bullet head. Beyond the turbulent region immediately behind the bullet head and in between the stripped bullet fragments, the gas tends to re-expand essentially at the sound speed of the ambient gas and fill the cavity. Such re-expansion causes gas to rebound on the axis resulting in a reflected shock. The effect of this shock can be seen in Figure 1a around the symmetry axis as a complex shock region that is narrower than the bowshock.

Note, that in Figure 1b the surface brightness of the system increases toward the bullet head. This differs from the surface brightness behaviour observed in the $`\eta `$ Carinae strings. As was noted by (Weis et al.1999) the string surface brightness decreases toward the head so that the true string tip might not even be detected. Recall, however, that Figure 1b is not a true synthetic observation image in a certain line but rather a projected map of total radiative energy losses in the system. In addition illumination by stellar photons is not included in our calculation. Thus in the simulation the bullet head and its vicinity, where the density and temperature are highest so that $`n^2\mathrm{\Lambda }(T)`$ are large, appear brightest. In reality gas in the bullet head and near it may be too hot to be visible optically as was suggested by (Redman et al.2002). In our further work we plan to include tracking of full ionization dynamics that will allow us to produce true synthetic observation images in a particular emission line or a combination of lines.

### 3.2 Velocity Distribution and Hubble-type Flow

Another key property of the strings is the presence of Hubble-type flows in the bullet wake. The left panel of Figure 2 shows the distribution of the total velocity $`v_{tot}(x)`$ along the symmetry axis of the bullet/wake as a function of distance from the bullet head. The linear velocity decrease from the maximum value of $`210`$ km $`\mathrm{s}^1`$from head to base is clear aside from some minor fluctuations arising due to the unsteadiness of the downstream flow222It should be noted that we observed similar behaviour of the total velocity distribution in other systems with the cooling parameter up to $`\psi _b0.25`$. For systems approaching the adiabatic regime with $`\psi >1.0`$ the velocity distribution takes a more complicated form.. Note that very little deceleration of the bullet head has occurred. Given that the initial bullet velocity is 235 km $`\mathrm{s}^1`$we see that after 251.1 yrs., in which the bullet material traveled the distance of $`1.810^{17}`$ cm $`0.058`$ pc, the bullet material lost only $``$ 10% of its original velocity. That is due to the fact that the central fragment retains most of the original bullet mass while it has a rather small cross-section due to the cooling-induced contraction.

(Weis et al.1999) quote a velocity slope in the $`\eta `$ Carinae strings ranging from 2590 km $`\mathrm{s}^1`$pc-1 (String 5) to 3420 km $`\mathrm{s}^1`$pc-1 (String 2). In our simulation we find a velocity slope of 3600 km $`\mathrm{s}^1`$pc-1 which is about 5% - 40% higher. Note, that the length of String 5 is practically equal to the distance over which we allowed the bullet to propagate, however String 2 is about 25% shorter (, Weis et al.1999).

One should be cautious, however, in directly comparing the total velocity distribution shown in the left panel of Figure 2 to the observationally determined velocity distributions in the strings of $`\eta `$ Carinae (cf. Figure 6 in (, Weis et al.1999)). Quantities, more closely resembling the ones that are obtained observationally, should be employed in order to make the comparison between the numerical results and observations relevant. For example, one can create emissivity-weighed total velocity $`v_{emis}(x)`$ maps of the flow. Distribution of such a quantity along the x-axis is shown in the right panel of Figure 2. $`v_{emis}(x)`$ was obtained according to the formula

$$v_{tot}(x)=\frac{_jn_j^2\mathrm{\Lambda }_jv_j}{_jn_j^2\mathrm{\Lambda }_j}\frac{_{y_{min}}^{y_{max}}n^2(x,y)\mathrm{\Lambda }(x,y)v_{tot}(x,y)𝑑y}{_{y_{min}}^{y_{max}}n^2(x,y)\mathrm{\Lambda }(x,y)𝑑y},$$
(7)
where $`i`$ and $`j`$ are the cell indices in the x- and y-directions, $`v_j`$ is the total velocity in the cell $`j`$, and in the summation in the numerator we included only the cells with non-zero total velocity. The distribution of this quantity is significantly more noisy than the total velocity cross-cut. However, the local maxima of $`v_{emis}(x)`$ roughly tend to fall on the same line as that in the previous figure giving some indication of the presence of the Hubble-type flow.

The distance traveled by the bullet in our simulation is practically equal to the length of String 5 and is factors of $`1.53.0`$ shorter than the rest of the strings (except for String 2) (, Weis et al.1999). If we assume that the true length of strings is not their observed length but rather the distance of the tip from the star, than the $`\eta `$ Carinae strings are the factors of $`3.25.56`$ longer (, Weis et al.1999) than the computational domain in our simulation. However, the kinematic age of the bullet in our simulation is about 67% larger than that of the $`\eta `$ Carinae strings: 251.1 years vs. 150 years. This indicates that in our simulation the bullet velocity was approximately a factor $`2.55.0`$ (or a factor $`5.359.3`$ if we assume that the string base is located near the central star) smaller than in the case of $`\eta `$ Carinae strings, i.e. the bullet velocity should be $`5901175`$ km $`\mathrm{s}^1`$(or $`12602185`$ km $`\mathrm{s}^1`$) with respect to the central star. An increased bullet velocity would decrease the width of the bullet wake and increase the length-to-width ratio making a better match between our simulations and the strings of $`\eta `$ Carinae. However, a significant increase in velocity, and therefore Mach number, results in much higher temperatures in the bullet wake and the bow shock. That, in turn, increases the cooling parameter $`\psi _b`$ eventually sending the system into the adiabatic regime which no longer exhibits the high length-to-width ratio of the wake and the Hubble-type downstream flow333We observed that effect in a simulation of Mach 200 bullet propagation. The cooling parameter for that case was $`\psi _b=3.7`$. We do not discuss this simulation in this paper however its animation can be accessed at the web site mentioned in the acknowledgment section.. The easiest way to reconcile the need for the significantly higher bullet velocity, required in order to obtain the correct kinematic age, with the need for the relatively low post shock temperatures, required to decrease the bullet wake width, is to embed the bullet in the ambient flow that itself moves with significant speed relative to the central star. There is some observational evidence for the presence of such high velocity ambient flow (Weis, Duschl, & Bomans, Weis et al.2001). For example, string 1 may be embedded in the material moving with velocity $`500600`$ km $`\mathrm{s}^1`$with respect to the central star, whereas the maximum velocity of the string material is $`995`$ km $`\mathrm{s}^1`$(, Weis et al.1999).

## 4 CONCLUSIONS

In this paper we have presented a numerical study of a hypersonic radiative cosmic bullet in an attempt to model the strings observed in the $`\eta `$ Carinae Homunculus nebula. Our simulations follow from the discussion by (Redman et al.2002).

Our principal conclusions are as follows: (1) hypersonic radiative bullets are capable of producing structures with high length-to-width ratios (between 6 and 10 for our study); (2) the dominant process responsible for bullet destruction is instability formation at the bullet upstream interface leading to separate fragmentation episodes; these result in the formation of periodic “ring-like” structures in the bullet wake; (3) the simulations do show the presence of Hubble-type flows along the axis in the bullet wake in terms of a linear decrease of the total velocity downstream from the bullet head. How these flows appear observationally remains an open question.

Thus we conclude that the bullet model appears as a good candidate for the strings of $`\eta `$ Carinae, though further work is needed. If this model finds continued success then theorists will confront the issue of how such high velocity bullets are generated, i.e. at the star or at larger distances. Confirming the existence of such bullets may be helpful in determining the nature of the processes, occurring within the star, which were part of its various eruptions.

The astrophysical applicability of the model, discussed in this paper, might extend beyond the $`\eta `$ Carinae strings. One of the most spectacular examples of nebular systems with long thin structures is the protoplanetary nebula CRL 618 and, in particular, its shocked lobes. The most prominent features of the lobes are the periodic rings similar to the ones present in our simulation. As it was discussed, such rings arise naturally as a consequence of the bullet fragmentation events. There is also some evidence for the velocity increase in the lobes from the base to the tip (Sánchez Contreras, Sahai, & Gil de Paz, Sánchez Contreras et al.2002). Other examples include, but are not limited to, other planetary nebulae, e.g. strings in NGC 6543 (, Weis et al.1999), HH objects, etc.

Further developments of the model should include the study of the details of the formation of the initial fragmentation spectrum, investigation of the importance of mass-loading due to the ablation of the bullet head and bullet fragments, and inclusion of more realistic description of ionization dynamics and radiative cooling in the system. The latter will allow us to produce more realistic synthetic observation images and synthetic spectra for better comparison with observational data.

This work was supported in part by the NSF grant AST-9702484, NASA grant NAG5-8428 and the Laboratory for Laser Energetics under DOE sponsorship. The most recent results and animations of the numerical experiment, described above as well as the ones not mentioned in the current paper, can be found at www.pas.rochester.edu$/^{}$wma.


text:
# Reply to “Comment on ‘Origin of combination frequencies in quantum magnetic oscillations of two-dimensional multiband metals’ ” by A.S. Alexandrov and A.M. Bratkovsky [cond-mat/0207173]

## Abstract

In their comment on the paper (Phys. Rev. B 65, 153403 (2002); cond-mat/0110154), Alexandrov and Bratkovsky (cond-mat/0207173) argue that they correctly took into account the chemical potential oscillations in their analytical theory of combination frequencies in multiband low-dimensional metals by expanding the free energy in powers of the chemical potential oscillations. In this reply, we show that this claim contradicts their original paper (Phys. Rev. B 63, 033105 (2001)). We demonstrate that the condition given for the expansion is mathematically incorrect. The correct condition allows to understand the limits of validity of the analytical theory.

It is well-known that, due to the quantization of the electronic energy spectrum of metals into discrete Landau levels in the presence of a magnetic field, the chemical potential is expected to oscillate with the magnetic field when the number of electrons is kept constant. In two-dimensional (2D) multiband metals, it has been predicted initially in the framework of numerical studies  that such chemical potential oscillations are responsible for the appearance of additional oscillations of the magnetization whose frequencies are combinations of the independent band frequencies $`f_\alpha `$. The number of electrons is fixed experimentally, independently of the dimensionality of the energy spectrum. For this reason, combination frequencies are in principle possible as well in 3D metals as in 2D multiband metals. However, chemical potential oscillations effects are not observed in multiband 3D metals. So, the mechanism of appearance of combination frequencies has to clearly point out the difference between low-dimensional metals and 3D metals.

Alexandrov and Bratkovsky are the first authors to propose an analytical derivation for these combination frequencies . However, they did not mention clearly the difference between low-dimensional metals and 3D metals in the framework of their theory. More precisely, an important point noticed for one-band 2D metals was missing in their derivation : the analytical treatment of the chemical potential oscillations effects leads to a system of nonlinear equations . In 3D metals, the resolution of these equations is trivial because the oscillating part $`\stackrel{~}{\mu }`$ of the chemical potential is of the order of $`\mathrm{}\omega _c\sqrt{\mathrm{}\omega _c/\epsilon _F}`$ (where $`\omega _c`$ is the cyclotron pulsation and $`\epsilon _F`$ the Fermi energy). As a result, the magnetization oscillations are not sensitive to the chemical potential oscillations effects: the expression for the magnetization is (with high accuracy) the same for a fixed chemical potential $`\mu `$ as for a fixed number of electrons $`N`$.

In 2D metals $`\stackrel{~}{\mu }\mathrm{}\omega _c`$, that leads to the observable difference between low-dimensional metals and 3D metals. This fact has not been noticed in the article  and has motivated the analytical work  and its principal statement that the chemical potential oscillations appearing in the arguments of the Fourier components (of the grand canonical potential or of the magnetization oscillations) were not taken into account by Alexandrov and Bratkovsky.

In their comment , these authors have addressed principally two criticisms to the paper . First, (i) the chemical potential oscillations were correctly taken into account in their previous work  : they argue that they made in fact an expansion (as in the commented paper ), but did not judge to mention explicitly this technical point . Secondly, (ii) the analytical formula for the combination frequencies amplitude they derived is accurate even at zero temperature in a clean two-band 2D metal. We will show in this reply that the argument (i) is in contradiction with the original paper  and is thus not receivable. The point (ii) is far from being correct.

The main point of the Ref.  is to express the oscillating part $`\stackrel{~}{F}`$ of the free energy as

$$\stackrel{~}{F}=\stackrel{~}{\mathrm{\Omega }}(\mu )\frac{1}{2\rho }\left[\left(\frac{\stackrel{~}{\mathrm{\Omega }}}{\mu }\right)_H(\mu )\right]^2$$
(1)
where $`\rho `$ is the total density of states of a two-dimensional multiband metal, and $`\stackrel{~}{\mathrm{\Omega }}(\mu )`$ is the oscillating part of the grand canonical potential, which is an explicit function of the chemical potential $`\mu `$. At a constant number of electrons $`N`$, the chemical potential $`\mu `$ oscillates with the magnetic field and can be written as the sum of a constant part $`\mu _0`$ (independent of the magnetic field) and an oscillating part $`\stackrel{~}{\mu }`$ given by

$$\stackrel{~}{\mu }=\frac{1}{\rho }\left(\frac{\stackrel{~}{\mathrm{\Omega }}}{\mu }\right)_H(\mu )=\mu \mu _0$$
(2)
As stressed in the Ref. , this equation (2) is a self-consistent nonlinear equation to solve in order to determine the dependence of the chemical potential $`\mu `$ on the magnetic field $`H`$.

The grand canonical potential $`\stackrel{~}{\mathrm{\Omega }}(\mu )`$ depends on the chemical potential $`\mu `$ through the expression

$$\stackrel{~}{\mathrm{\Omega }}(\mu )=\underset{\alpha }{}\underset{r=1}{\overset{\mathrm{}}{}}(1)^lA_\alpha ^r\mathrm{cos}\left(2\pi r\frac{\mu _\alpha }{\mathrm{}\omega _{c\alpha }}\right)$$
(3)
where

$$2\pi \frac{\mu _\alpha }{\mathrm{}\omega _{c\alpha }}=2\pi \frac{\mu \mathrm{\Delta }_\alpha }{\mathrm{}\omega _{c\alpha }}=\frac{f_\alpha }{H}$$

are the arguments entering in the Fourier components. Here $`\omega _{c\alpha }`$ is the cyclotron pulsation with the effective mass $`m_\alpha `$, $`\mathrm{\Delta }_\alpha `$ is the $`\alpha `$ band-edge, and

$$A_\alpha ^r=\frac{m_\alpha \omega _{c\alpha }}{2\pi ^3}R_T(r)R_D(r)$$

is the amplitude for the harmonic $`r`$ in the band $`\alpha `$ with

$$R_T(r)=\frac{\lambda _r}{\mathrm{sinh}\lambda _r},R_D(r)=\mathrm{exp}\left(2\pi r\frac{\mathrm{\Gamma }_0}{\mathrm{}\omega _{c\alpha }}\right),$$

$`T`$ the temperature, $`\lambda _r=2\pi ^2rk_BT/\mathrm{}\omega _{c\alpha }`$, and $`\mathrm{\Gamma }_0`$ the relaxation rate at $`H=0`$. Since for a constant number of electrons the chemical potential oscillates with the magnetic field, so does the quantity $`f_\alpha `$. For this reason, the explicit expression (3) can not be seen as a Fourier series ($`f_\alpha `$ is not a frequency, which by definition has to be independent of the magnetic field).

This difficulty can be overcomed by expanding the free energy $`\stackrel{~}{F}(\mu )`$ (expression (1)) in powers of the oscillating part $`\stackrel{~}{\mu }`$ of the chemical potential. By developing separately the first term and the last term of the right-hand side of the equation (1), we straightforwardly obtain keeping terms up to the order of $`\rho \stackrel{~}{\mu }^2`$

$$\stackrel{~}{\mathrm{\Omega }}(\mu )\stackrel{~}{\mathrm{\Omega }}(\mu _0)+\left(\mu \mu _0\right)\left(\frac{\stackrel{~}{\mathrm{\Omega }}}{\mu }\right)_H(\mu _0)=\stackrel{~}{\mathrm{\Omega }}(\mu _0)+\frac{1}{\rho }\left[\left(\frac{\stackrel{~}{\mathrm{\Omega }}}{\mu }\right)_H(\mu _0)\right]^2$$
(4)
and

$$\frac{1}{2\rho }\left[\left(\frac{\stackrel{~}{\mathrm{\Omega }}}{\mu }\right)_H(\mu )\right]^2\frac{1}{2\rho }\left[\left(\frac{\stackrel{~}{\mathrm{\Omega }}}{\mu }\right)_H(\mu _0)\right]^2.$$
(5)
Therefore the free energy becomes

$$\stackrel{~}{F}=\stackrel{~}{\mathrm{\Omega }}(\mu _0)+\frac{1}{2\rho }\left[\left(\frac{\stackrel{~}{\mathrm{\Omega }}}{\mu }\right)_H(\mu _0)\right]^2.$$
(6)
Obviously after equations (4) and (5), Fourier harmonics with combination harmonics are produced by both terms of the right-hand side of the expression (1), the contribution of the first term being reduced to one half due to a partial cancellation with the second term.

Alexandrov and Bratkovsky argue  that the free energy $`\stackrel{~}{F}`$ is expanded in powers of $`\stackrel{~}{\mu }\mu _0`$ in their original work ; in fact, the authors asserted just after giving the expression (1) (their Eq. (12)) in the Ref.  : “It is the last term, which yields combination Fourier harmonics with the combination frequencies $`f=rf_\alpha \pm r^{}f_\alpha ^{}`$”. Thus, this sentence directly contradicts the fact that the free energy has been expanded in the paper : the chemical potential oscillations were explicitly not taken into account in the trigonometric arguments.

It is worth noting that a similar claim has been addressed independently by Kishigi and Hasegawa , who wrote concerning the authors of Ref. :“their result for the free energy \[…\] is formally correct but they did not take account of the magnetic-field dependence of the $`f_\alpha `$, which cannot be neglected in two-dimensional systems. As a result their analysis of the de Haas-van Alphen oscillation for the fixed $`N`$ system (canonical ensemble) is insufficient and their conclusions on Fourier-transform intensities are incorrect”.

The similitude of the formulae obtained in Ref.  and Ref.  has a simple explanation. We note by comparing the expressions (1) and (6) that the squared term changes sign after the expansion. It implies that if we ignore by hand the oscillations of the chemical potential in the quantity $`f_\alpha `$, the amplitude obtained for the combination harmonics is fortuitously the same (a sign apart) after and before the development.

The second point (ii) of discordance is in fact related to the validity of the analytical expansion of $`\stackrel{~}{F}`$ (or of $`\stackrel{~}{\mathrm{\Omega }}`$, or of the magnetization oscillations $`\stackrel{~}{M}`$) in powers of the oscillating part $`\stackrel{~}{\mu }`$ of the chemical potential $`\mu `$. The condition of validity for it, $`|\stackrel{~}{\mu }|\mu _0`$ (which is fulfilled in the whole regime of magnetic quantum oscillations) given in , is actually mathematically incorrect. The correct condition for performing the developments (4-6) is revealed when considering the explicit form (3) for $`\stackrel{~}{\mathrm{\Omega }}(\mu )`$. The basic point is that the oscillating part $`\stackrel{~}{\mu }`$ enters in the arguments of trigonometric functions for which the approximation

$$\mathrm{cos}(x_0+\stackrel{~}{x})\mathrm{cos}(x_0)$$

is valid under the condition $`|\stackrel{~}{x}|1`$ (and not $`|\stackrel{~}{x}|x_0`$, take e.g. $`\pi /2x_0`$). The expansion (6) is therefore valid provided that

$$2\pi r\left|\frac{\stackrel{~}{\mu }}{\mathrm{}\omega _{c\alpha }}\right|1$$
(7)
for all significant harmonics $`r`$. Obviously, this latter condition (7) is much more restrictive than the condition given by the authors of the Ref.  and depends on the degree $`r`$ of the harmonic (it is stronger for higher harmonics). It holds simultaneously for the expansion of all quantities which are functions of $`\mu /\omega _{c\alpha }`$ through trigonometric arguments such as $`\stackrel{~}{\mathrm{\Omega }}`$, the magnetization oscillations $`\stackrel{~}{M}`$ or even $`\stackrel{~}{\mu }`$ itself (see equation (2)).

It has already been noticed in Ref.  that $`\stackrel{~}{\mu }`$ is naturally reduced by the presence of multiple bands independently of the temperature or the impurity reduction factors. However, at very low temperatures and in clean 2D two-band metals (the most unfavorable case), many harmonics $`r`$ are significant, and the condition (7) is not fulfilled for all terms because in this case $`|\stackrel{~}{\mu }|`$ is of the order of a few tenths of $`\mathrm{}\omega _{c\alpha }`$. Then, higher powers of $`\stackrel{~}{\mu }`$ have to be considered and the analytical treatment of the chemical potential oscillations effects is not obvious. In this regime of so-called strong chemical potential oscillations the validity of the analysis of the oscillations in terms of Fourier series is questioned at small but finite temperatures (this is not a property established a priori; the use of the Fourier analysis has thus to be justified).

The importance of higher powers of $`\stackrel{~}{\mu }`$ has been demonstrated analytically and numerically in the case of one-band 2D metals . In these papers  the full nonlinear equation (2) is considered and solved at zero temperature. As a result, the drops of the magnetization oscillations occur at integer values of the ratio $`\mu _0/\mathrm{}\omega _c`$, while they occur at half-integer values when neglecting $`\stackrel{~}{\mu }`$ in the trigonometric arguments. Furthermore, harmonic amplitudes are found to differ strongly in the two situations, especially for high harmonics .

In the two-band 2D metals, the combination frequencies appear in the first order expansion in powers of $`\stackrel{~}{\mu }`$. Alexandrov and Bratkovsky argue  with the help of an analytical and a numerical estimates that higher powers of $`\stackrel{~}{\mu }`$ can be neglected to describe the magnetization oscillations with a relatively good accuracy even at zero temperature and in clean samples. In their analytical estimate, they only keep the first harmonic in the self-consistent equation (2) for the chemical potential oscillations, which is valid in fact at zero temperature only at a small Dingle factor $`R_D(1)`$. In this regime, the oscillations of $`\stackrel{~}{\mu }`$ are small and particularly smooth: for this reason the linear approximation is naturally expected to be quite good since the condition (7) is then fulfilled. On the contrary, in their numerical study they consider all the harmonics in Eq. (2), which is relevant in the regime of strong oscillations of $`\stackrel{~}{\mu }`$. We can however cast some doubts on their result shown in the Fig. 1 of Ref. : surprisingly, at small Dingle factors (i.e. when the expansion is possible), the accuracy of the linear approximation is increasing when the Dingle factor $`R_D(1)`$ decreases; this is completely opposite to that can be expected.

In conclusion, we have demonstrated that the chemical potential oscillations were not correctly taken into account in the Ref. . Careful numerical studies are still needed to analyze the regime of strong chemical potential oscillations where the expansion of the quantities in powers of $`\stackrel{~}{\mu }`$ is not convergent.

I thank V.P. Mineev for his advices and careful reading of the manuscript.


text:
# Quantum size effects in Pb islands on Cu(111): Electronic-structure calculations

## Abstract

The appearance of ”magic” heights of Pb islands grown on Cu(111) is studied by self-consistent electronic structure calculations. The Cu(111) substrate is modeled with a one-dimensional pseudopotential reproducing the essential features, i.e. the band gap and the work function, of the Cu band structure in the  direction. Pb islands are presented as stabilized jellium overlayers. The experimental eigenenergies of the quantum well states confined in the Pb overlayer are well reproduced. The total energy oscillates as a continuous function of the overlayer thickness reflecting the electronic shell structure. The energies for completed Pb monolayers show a modulated oscillatory pattern reminiscent of the super-shell structure of clusters and nanowires. The energy minima correlate remarkably well with the measured most probable heights of Pb islands. The proper modeling of the substrate is crucial to set the quantitative agreement.

The confinement of valence electron states in low-dimensional systems has a strong influence on the size distributions of nanostructures grown or produced in experiments. In clusters of alkali metal atoms the confinement occurs in all three dimensions and the structures corresponding to closed electronic shells are the most stable and therefore the most abundant ones. Knight et al. (1984) Alkali metal nanowires Yanson et al. (1999) and metallic overlayers on solid surfaces Hinch et al. (1989); Crommie et al. (1993); Otero et al. (2002) exhibit two- and one-dimensional (1D) confinement, respectively. In these systems the sinking of the bottom of a new subband below the Fermi level is accompanied with an increase in the total energy destabilizing the system. The detailed understanding of the mechanisms controlling the growth of nanoobjects is of vital importance when producing, and eventually manufacturing, highly-organized atomic-scale structures in nanodevices.

In this work we focus on the 1D confinement in Pb islands on the Cu(111) surface. The growing of Pb occurs in the  direction. Hinch et al. Hinch et al. (1989) studied the island height distribution of Pb using the He atom scattering (HAS). Recently, Otero et al. Otero et al. (2002) have determined the height distribution of Pb islands using the scanning tunneling microscope (STM) up to heights over 20 Pb monolayers (ML’s). The lateral dimension of the islands is large, of the order of 500 Å, justifying the 1D modeling of their energetics. The vertical confinement in the Pb overlayer is due to the potential barrier between Pb and the vacuum and the energy gap in the projection of the Cu bulk bands in the  direction. A closed shell occurs periodically when the thickness $`D`$ of the Pb layer satisfies $`D=n\lambda _F/2`$, where $`\lambda _F`$ is the Fermi wavelength of Pb and $`n`$ is an integer. The effect of this shell structure on the physical properties was recognized already in early calculations Schulte (1976) of unsupported jellium slabs and the stability as a function of the thickness was discussed on the basis of density-functional theory (DFT) calculations already two decades ago. Feibelman (1983) Later jellium and pseudopotential calculations Sarría et al. (2000); Kiejna et al. (1999); Saalfrank (1992); Materzanini et al. (2001); Wei and Chou (2002) have dealt with quantum size effects in unsupported metallic slabs, while the effects of the substrate have not been explicitly considered.

We model the whole system consisting of the Pb overlayer and the Cu(111) substrate by using the stabilized jellium (SJ) modelPerdew et al. (1990) for the overlayer and a 1D pseudopotential for the substrate and obtain results with a quantitative predictive power. We calculate electronic structure using local-density approximation (LDA). The delocalized electron character of Pb at the Fermi level justifies the use of the SJ model with $`r_s=2.3a_0`$ to describe the Pb overlayer. The SJ model gives the work function of 4.1 eV, which is close to the experimental value of 4.0 eV for Pb. Ashcroft and Mermin (1981) The relevant feature giving physical insight is that the SJ model allows us to simulate overlayers of any thickness. For the Cu(111) surface we have constructed a new 1D pseudopotential. We start from the 1D model potential by Chulkov et al., Chulkov et al. (1999) cosinelike bulk part of which reproduces correctly the experimental energy gap of Cu in the (111) direction. The position of the Fermi level with respect to the model potential determines the electron density. We subtract then the interactions between valence electrons from the effective potential within the LDA and obtain a local 1D pseudopotential. The pseudopotential is constructed to give the experimental Cu(111) work function of 4.94 eV.Chulkov et al. (1999) The self-consistent screening of this 1D pseudopotential within the LDA accounts correctly for the energy gap in the (111) projection, both in the width and in the position with respect to the Fermi level. The details of the construction are described in a forthcoming paper. Ogando et al. (2003)

In the numerical calculations a semi-infinite crystal is simulated by the slab geometry in which 25 layers of Cu are covered on both sides by Pb-SJ with the desired thickness. In Fig. 1 we show the effective potential in the case of 3 Pb ML. One ML corresponds to the thickness $`d=5.441a_0`$. In Cu, the potential maxima due to the 1D-pseudopotential mimics the orthogonalization of the valence states to the core states. At the Pb/Cu interface the potential is affected by the alignment of the Fermi levels of the substrate and the overlayer accompanied by the charge transfer from Pb to Cu. In the Pb layer the potential is a rather flat well and at the Pb-vacuum interface it forms a dipole-layer step. The dark gray region denotes the energy gap induced by the Cu potential and the light gray region gives the energy range between the potential in Pb and the bottom of the Cu bulk band. These two regions will accommodate the quantum-well states (QWS) localized perpendicular to the surface and mainly in the Pb overlayer. QWS’s at the Cu band gap determine the energy shell structure.

The QWS eigenenergies (bottoms of two-dimensional paraboloid bands) are shown in Fig. 2 as a function of the number of Pb ML’s completed. They are compared with the results measured by STM. Otero et al. (2000) The agreement is good for thick Pb coverages (thicker than 6 ML of Pb). In particular, both the theory and experiment give a QWS at $`0.65`$ eV for every even number of ML’s (no experimental data are available for 18, 20, and 22 ML). For coverages less than 7 ML the correspondence is worse. This disagreement can reflect the compression of the Pb layers in both the vertical (perpendicular to the interface) and in horizontal (parallel to the interface) directions omitted in our simple model. Materzanini et al. (2001); Sarría et al. (2000) The effect of the compression is to push up the eigenenergies the stronger the thinner the Pb overlayer is. Another important factor that can influence the QWS energies for thin coverages of Pb and which is not taken into account in our model is the interaction of 3d electrons of Cu with delocalized electrons of Pb. Due to the very localized character of 3d electrons this interaction takes place mostly at the Cu-Pb interface and disappears at distances of a few ML’s of Pb because of screening by Pb electrons. The overall agreement between the theoretical and measured eigenenergies for thick Pb overlayers gives confidence to our model when studying the behaviour of the total energy.

The total energy of the Pb/Cu(111) system can be written as the sum $`E_{Cu}+E_{Pb}+E_{CuPb}+E_{Pbvac}`$, where $`E_{Cu}`$ and $`E_{Pb}`$ are contributions calculated by using the bulk energies corresponding to the 1D-Cu(111) pseudopotential and the Pb-SJ, respectively. The sum of the energies due to the Cu-Pb and Pb-vacuum interfaces,$`E_{CuPb}+E_{Pbvac}`$, contains the oscillating part of the total energy. In order to see the oscillations more clearly also at high coverages Fig. 3 shows the oscillating part of the energy multiplied by the Pb thickness. As a continuous function of the Pb thickness the oscillations are regular and their wavelength is half the Fermi wavelength, $`\frac{\lambda _F}{2}=3.77a_0`$, of the Pb jellium.

The solid line in Fig. 3 links the points corresponding to completed Pb ML’s in the  direction. Because the interlayer spacing $`d3/2\frac{\lambda _F}{2}`$ an even-odd staggering of the energy as function of the ML’s is obtained. However, the above relation is not exact and therefore the staggering amplitude diminishes regularly and the phase of the staggering changes at the beats from that corresponding to minima at an even number of ML’s to minima at an odd number of ML’s or vice versa. The pattern resembles that of the supershell structures for atomic clusters de Heer (1993) and metallic nanowires. Yanson et al. (2000) Therefore we will speak in the present context also about the supershell structure although now the physical origin is different.

The arrows pointing upward in Fig. 3 denote the most abundant island heights measured by Otero et al. Otero et al. (2002) with the STM. Actually, the measurements give the percentage of area covered with islands of a given height and the highlighted heights correspond to maxima in this distribution. The correlation with the energy minima of our solid curve is remarkable. Even the even-odd phase changes (supershell structure) agree. This means that the minima of the total energy determines the most abundant island heights at the growing temperature which is around the room temperature. According to Fig. 3 the energy differences between adjacent systems with even and odd ML’s of Pb are of the order of 0.05 eV per surface atom, i.e. higher than the thermal energy. Menzel et al. Menzel et al. (2003) estimated the evolution of a seven-layer island from a five-layer island on Si(111) to occur through a barrier of 0.32 eV. They concluded that the thermodynamics determines the most abundant island heights but the kinetics is important in determining the height distribution.

However, there is some mismatch between the theory and experiment in Fig. 3. In experiments, islands of 13 ML’s of Pb are not especially abundant. The experimental situation in which the substrate is not flat but contains steps or terraces may affect the height distribution at this point Miranda and de Parga (2003). At the second beat our model predicts similar energies for the island heights of 17 and 18 ML. But according to experiments the 17 ML island should be more abundant. The disagreement may be caused by a small error in the phase of the energy oscillations in our model. In the beat regions the relative energies are very sensitive on this phase.

Hinch et al. Hinch et al. (1989) have measured the shell and supershell structure of the Pb island heights with the HAS technique. During the Pb deposition, a high intensity of the reflected He beam indicates smooth surfaces of islands of completed ML’s. The heights of these islands are obtained indirectly from the deposition rate and the exposition time. The experiment shows the even-odd staggering which has an opposite phase in the comparison with the STM results by Otero et al. Otero et al. (2002) and our calculations. The beat positions by Hinch et al. are located at coverages of 10-11 and 20-22 ML, which are slightly larger than those by Otero et al. and our predictions. Thus, it seems evident that the indirect measurement overestimates the Pb overlayer thickness by one ML.

We have calculated the total energy of free-standing Pb-SJ slabs. The results show oscillations with the same wavelength and a similar amplitude as the dashed curve in Fig. 3. But the oscillating pattern is now shifted about $`0.8a_01/7`$ Pb ML forward. This changes the shell and supershell structure so that energy minima are shifted and the beat positions occur at 6-7 and 14-15 ML’s, i.e., in a clear disagreement with the measurements. We have compared our unsupported Pb-SJ results with those obtained for unsupported Pb films from pseudopotential calculations. Materzanini et al. (2001); Wei and Chou (2002) These more sophisticated evaluations give similar wrong position of the first beat as our unsupported Pb-SJ slabs. The damping of the energy oscillations makes the second beat hardly recognizable in these calculations. In conclusion, the energy oscillations are very sensitive to the correct determination of the potential barriers and it is crucial to take the Cu(111) substrate into account in the electronic structure calculations.

The shell structure of simple metal overlayers is a result of the periodic sinking of new QWS’s below the Fermi level when the thickness of the overlayer increases. New QWS’s increase the density of states (DOS) at the Fermi level giving the oscillations in the energy with the wavelength of $`\lambda _F/2`$. Neglecting the damping we can write the oscillating part of the total energy as a function of the thickness $`D`$ of the slab as $`E_{osc}(z)=\mathrm{cos}(2\pi \nu D+\theta )`$, where $`\nu =2/\lambda _F`$ and $`\theta `$ is a phase that shifts the energy to the correct position depending on both the Pb-vacuum and the Cu-Pb interface properties. The energy values corresponding to $`N`$ completed ML’s are obtained from

$$E_{osc}(N)=\mathrm{cos}(2\pi \nu dN+\theta ).$$
(1)
For a fcc metal grown in the (111) direction the interlayer spacing is $`d=a/\sqrt{3}`$, where $`a`$ is the lattice parameter. The Fermi wavelength is $`\lambda _F=(\frac{2\pi }{Z_v})^{1/3}a`$, where $`Z_v`$ is the number of valence electrons. These relations give

$$\nu d=\frac{1}{\sqrt{3}}\left(\frac{12Z_v}{\pi }\right)^{1/3}0.903Z_v^{\frac{1}{3}},$$
(2)
and substituting this into Eq. 1 gives a general formula for $`E_{osc}`$ of a fcc metal grown in the (111) direction.

Using Eqs. (1) and (2) and setting $`\theta =0`$ we have plotted in Fig. 4 the oscillating part of the energy as a function of the number of ML’s for $`Z_v`$ = 2 and 4. The experimental counterparts of these curves are the abundance spectra of islands heights. For monovalent and divalent overlayer metals the model predicts strong peaks occurring regularly after a certain number of ML’s. However, this pattern may be suppressed in the real experiments, because the energy differences between neighboring completed ML’s are small and because the growing of a metastable island from the previous one requires several new ML’s. For tri- and tetravalent metals the energy oscillates more rapidly. For $`Z_v`$ = 4, $`\nu d`$ is close to a half integer and the supershell structure with beats is clear. Note that the increase of electron density as a function of $`Z_v`$ decreases the wavelength of $`E_{osc}(z)`$ as shown by the dashed curves in Fig. 4. This makes the determination of the stable island heights very sensitive to the variations of parameters such as the interplane distances determining the overlayer thickness or the penetration of the QWS wavefunctions into vacuum or/and into the Cu substrate.

In summary, we have analysed the electronic structures and stability of Pb islands grown on the Cu(111) surface as a function of the island height. The Cu(111) substrate is described by a new 1D pseudopotential and the Pb overlayer by the stabilized jellium model. As a function of Pb completed ML’s our model gives quantum well states in a good agreement with measurements. Otero et al. (2000) The total energy shows modulated odd-even oscillations resembling the supershell structure of simple-metal atomic clusters and nanowires. The pattern correlates well with the height abundance spectrum measured by Otero et al. Otero et al. (2002). We demonstrate that a proper modeling of the Cu(111) substrate plays a crucial role in predicting the beat positions of the abundance spectra.

We thank R. Miranda for useful discussions. This work was partially supported by the University of the Basque Country, Departamento de Educación del Gobierno Vasco, Spanish MCyT (MAT 2001-0946, PB98-0870-C02 and MAT 2002-04087-C02-01) and by the Academy of Finland through its Centre of Excellence Program (2000-2005).