text:
# Model study of adsorbed metallic quantum dots: Na on Cu(111)

## I Introduction

At certain faces of metals, such as the (111) face of noble metals, the surface electron states are confined to the vicinity of the top layer by the vacuum barrier on the vacuum side and the bandgap on the substrate side Luth . The electrons in these surface states form a two-dimensional nearly free electron gas Gartland ; Heimann ; Zangwill ; Ashcroft-Mermin . It has also been observed that when adsorbing one to several monolayers of alkali atoms on these surfaces a manifold of discrete standing wave states, so called quantum well states (QWS), perpendicular to the surface are formed. Lindgren-Wallden80 ; Lindgren-Wallden87 . These states can be detected, for instance, in photoemission spectroscopy (PES) Carlsson97 , inverse photoelectron spectroscopy (IPES) Dudde , two photon photoemission spectroscopy (2PPES) Fisher  and scanning tunneling microscopy (STM) Kliewer00 . A large amount of experimental data is available for the system Na on Cu(111) Lindgren-Wallden80 ; Lindgren-Wallden87 ; Carlsson97 ; Dudde ; Fisher ; Kliewer00 . The electronic structure and dynamics for this system have also been investigated by first-principles theoretical calculations Carlsson ; Hellsing2000 .

These localized surface states are of great interest since they play an important role in many physical processes like epitaxial growth Memmel , surface catalysis Bertel ; Lauritsen , molecular ordering Stranick  and adsorption Bertel97 . Experimental tools like STM and PES play an important role in the investigations, since they enable spatial and spectroscopic resolution of the electron states.

One important discovery is the confinement of surface state electrons in so called quantum corrals. These man-made nanoscale structures are formed by deliberately assembling adatoms to enclosed structures by STM Crommie . Due to the small size of the corrals, quantum effects are present, and both spatial and spectroscopic properties of the confined states can be studied experimentally. A natural way of forming low-dimensional structures on metal surfaces is by controlled growth of epitaxial layers. With an appropriate choice of deposition and annealing temperatures small islands, so called quantum dots (QD), with variable shapes and sizes may form Roder . The advantage of these structures, in comparison with the corrals, is that they are relatively stable at low temperatures. This enables the imaging and investigation of their properties without inducing structural damage. One quantum mechanical effect of the confinement is the increase of the surface state energy which in turn may lead to the depopulation of the surface state band and thereby changes in the surface properties.

In this paper we present calculated results for the electronic structure of Na on Cu(111), with the emphasis on describing the real-space resolved density of states nearby a sodium QD adsorbed on a sodium-covered Cu(111) surface. Previously, an all-electron density-functional theory (DFT) study of a free-standing Na layer in vacuum has been presented Wimmer . More recently, a DFT calculation using ultrasoft pseudopotentials for the free-standing Na layer as well as for the layers adsorbed on Cu(111) have been presentedCarlsson . A simple free electron model calculation for a free-standing Na QD has already been published by two of the present authorsvanja1 .

The calculations in this work are made in the context of the DFT lang83 ; dft . More specifically, the Rayleigh-quotient multigrid (RQMG) method in axial symmetrymgarticle1 ; mika  is used for the numerical solution of the ensuing Kohn-Sham equations. The electron-ion interaction is simplified using the jellium model lang83 , where the ions are replaced by a rigid positive background charge of constant density. This model has provided basic physical understanding of the electronic structures of simple metal surfaces lang83  thin metal films Schulte , vacancies and voids inside metals manninen75  and finite clusters of simple metal atoms heerbrack . Recently, also uniform cylindrical nanowires have been studied within the jellium model zabala1 ; zabala2 ; zabala3 .

The paper is organized as follows: Section II gives a short review over experimental results for the system Na on Cu(111), Section III describes the computational method used in the calculations, Section IV discusses the details of our jellium model and in Section V the results are presented and comparisons are made with experimental findings. Finally, Section VI gives the conclusions.

## II Na on Cu(111)

Alkali metals adsorbed on the closed-packed (111) surface of metals form hexagonal structures at saturated monolayer coverages, following approximately the underlying substrate structure Diehl . The first monolayer of Na on Cu(111) is observed to saturate at the coverage of $`\mathrm{\Theta }=4/90.44`$ Tang , corresponding to 4 Na atoms per 9 surface Cu atoms. The Na atoms thus form a hexagonal (3/2 $`\times `$ 3/2) structure and the Na atom spacing of $`7.43\mathrm{a}_0`$ is comparable to the atomic distance of $`6.92\mathrm{a}_0`$ in bulk Na.

The adsorption of Na atoms on the Cu(111) surface will induce a charge redistribution at the interface between the adlayer and the substrate. It has been seen from photoemission experiments Fisher ; Lindgren-Wallden88  that when the Na coverage is increased the Cu Shockley surface state decreases in energy. For coverages above $`\mathrm{\Theta }0.11`$ the surface state is shifted below the lower band edge of the local band gap of the Cu(111) surface, and is no longer visible in photoemission experiments. Two-photon photoemission experiments Fisher  indicate other unoccupied Na-induced states in the local bandgap at the Cu(111) surface, which will also decrease in energy with increasing Na coverage. For higher coverages, the lowest of these states will be down-shifted below the Fermi energy, and thus get occupied. At the saturated monolayer coverage, this state will be located about $`0.1\mathrm{eV}`$ Fisher ; Carlsson ; Lindgren-Wallden88  below the Fermi energy at the $`\overline{\mathrm{\Gamma }}`$-point of the surface Brillouin zone (SBZ). The corresponding next lowest state will be located $`2.1\mathrm{eV}`$ above the Fermi level Fisher  at the $`\overline{\mathrm{\Gamma }}`$-point. The lower of the two states has one node in the $`z`$-direction, while the higher has two.

If the Na atom deposition continues after the first monolayer (ML) is completed, a second layer will start to grow. Recent STM measurements Kliewer\_Berndt\_SS01 ; Kliewer\_Berndt\_PRB01  indicate that the second monolayer of Na grows via the formation of compact islands with hexagonal atomic arrangement. Normal photoemission experiments Carlsson97 ; Carlsson94  indicate that when the second monolayer grows the emission intensity due to one-monolayer states decreases gradually, and for coverages above 1.3 ML a new peak, approximately $`0.1\mathrm{eV}`$ above the Fermi energy appears. This peak is ascribed to the two-monolayer thick parts, and the energy is shifted to somewhat lower values as the coverage is increased.

Some theoretical attention has also been paid to Na on Cu(111), including the island growth. Free-electron model calculations have been performed for circular Na vanja1  as well as hexagonal Ag CrampinPRL ; CrampinSS  and Na Berndt  free-standing islands. All-electron calculations for an unsupported monolayer of Na Wimmer  and first-principles slab calculations for one atomic Na layer in (2 $`\times `$ 2) and (3/2 $`\times `$ 3/2) adsorbate structures on Cu(111) Carlsson  have been presented. In the present paper, we report jellium model calculations for an unsupported monolayer of Na and a cylinder shaped free-standing Na QD. We also present two-density-jellium calculations for the system Na on Cu(111), where we have modeled the underlying Cu(111) substrate by using a lower density slab to mimic the decay of the surface states into the substrate. Comparison is made with experiments and previous theoretical calculations.

## III Computational methods

In the Kohn-Sham scheme of DFT, one solves the electron density $`n(𝐫)`$ of the system selfconsistently from a set of equations. One of these equations is the single particle Schrödinger equation. The models used in this work are axially symmetric. Thus the Schrödinger equation is separable, and the wave functions can be written as products

$$\psi _{m𝐤n}(r,z,\varphi )=e^{im\varphi }U_{m𝐤n}(r,z).$$
(1)
Above, $`m`$ is the azimuthal quantum number implied by the axial symmetry while $`n`$ differentiates between orthogonal states with same $`m`$ and $`𝐤`$. In the calculations involving the infinite monolayer, two k-vectors are used as explained in Sec. IV. The external potential of the systems studied in this work is caused by the positive background charge $`n_+(𝐫)`$. The effective potential $`V_{\mathrm{eff}}`$ includes also the Hartree potential of the electron density and the exchange-correlation potential $`V_{\mathrm{XC}}`$, which we treat in the local density approximation (LDA) xc . The electron density $`n(𝐫)`$ is obtained by summing single-electron densities with the occupation numbers $`f_{m𝐤n}`$. The degeneracies of the states are taken into account by the factor $`(2\delta _{0m})`$ and the occupation numbers $`f_{m𝐤n}`$ obey the Fermi-Dirac statistics with a Fermi level ($`E_F`$) so that the system is neutral. A finite temperature of 1200 K is used to stabilize the solution of the set of equations. Thus, in the present axial symmetry ($`𝐫=(r,z)`$) the Kohn-Sham equations read as

$$\frac{1}{2}\left(\frac{1}{r}\frac{}{r}+\frac{^2}{r^2}\frac{m^2}{r^2}\frac{^2}{z^2}+2V_{eff}\right)U_{m𝐤n}(𝐫)=\epsilon _{m𝐤n}U_{m𝐤n}(𝐫)$$
(2)
$$n(𝐫)=2\underset{m𝐤n}{}(2\delta _{0m})f_{m𝐤n}|U_{m𝐤n}(𝐫)|^2,$$
(3)
$$V_{\mathrm{eff}}(𝐫)=\mathrm{\Phi }(𝐫)+V_{\mathrm{XC}}(𝐫),$$
(4)
$$\left(\frac{1}{r}\frac{}{r}+\frac{^2}{r^2}+\frac{^2}{z^2}\right)\mathrm{\Phi }=4\pi \left[n_{}(𝐫)n_+(𝐫)\right].$$
(5)
The Schrödinger equation (2) is solved using the Rayleigh-quotient multigrid (RQMG) method mgarticle1  which has been implemented in various geometries, including the axial symmetrymgarticle1 ; mika . In the RQMG method, the Rayleigh quotient $`\psi |H|\psi /\psi |\psi `$ on the finest level grid is directly minimized, the orthogonality constraint being taken into account by a penalty functional. The Poisson equation (5) is solved for the electrostatic potential $`\mathrm{\Phi }(𝐫)`$ using a standard multigrid method brandt1 .

To obtain selfconsistency, we use the simplest possible potential mixing scheme,

$$V_{in}^{i+1}=AV_{out}^i+(1A)V_{in}^i.$$
(6)
The largest system of this work contains 2550 electrons, the diameter of the supercell being $`170`$Å. Obtaining selfconsistency in such a system requires a very small $`A`$ value of 0.005. Otherwise the charge sloshing results in divergence. More sophisticated mixing strategiesWangPRB01  will be indispensable in the future calculations. However, because of the simplicity of our model systems, we can very accurately estimate an initial guess for the selfconsistent effective potential of large systems using the more easily convergent smaller systems as reference.

In our largest calculation a grid of $`319\times 95`$ points is used for the presentation of the wave functions, potential and density. Taking into account the unoccupied states needed in the modeling, up to 2400 different states have to be solved at every selfconsistency iteration. Luckily, it is straightforward to parallelize the calculation over the 65 different $`m`$ values and the two k-points (see below). Moreover, the RQMG methodmgarticle1  handles this part of the calculation with optimal efficiency.

The two-jellium model for the surface described in Sec. IV results in an asymmetric density distribution with a surface dipole. Thus the electrostatic potential on the substrate side is a constant different from that on the vacuum side. For the Poisson equation, we thus use the boundary condition of zero derivative on the substrate side, and that of zero value on the vacuum side. Solving the Poisson equation, the boundaries above and below the system are extended to a distance five times greater than in the case of the wave functions.

In this work, we calculate the local density of states (LDOS) above a surface at distances corresponding to those typical in STM measurements (of the order of $`20\mathrm{a}_0`$). At such distances, the amplitude of the wave function decreases by several orders of magnitude. This kind of modeling is thus a serious test for the RQMG-method. We have checked the accuracy of our method for the spherical harmonic oscillator and a model hydrogen atom potential, for which the wave-functions are known analytically. The evanescent tails of the wave functions solved with the RQMG method agree with the analytical ones, even when the amplitude of the wave functions has dropped by 20 orders of magnitude. This level of accuracy is beyond the reach of plane-wave methods, where periodic boundary conditions are necessary, and which provide a uniform accuracy across the calculation volume, resulting in spurious oscillations in the vacuum parts of the system.

## IV Modeling the system

We are interested in the system of a monolayer-thick Na QD on the complete Na monolayer on the Cu(111) surface. We know from experiments that these islands are approximately hexagonal in shape, following the underlying structure of the Na monolayer.

In order to interpret recent STM data for these types of systems, mapping the energy resolved real-space electron density near a QD is necessary. First, it is of interest to find what level of theoretical modeling is required. It has been shown previously CrampinPRL ; CrampinSS ; Berndt  that simple two-dimensional ’particle-in-a-box’ calculations give qualitatively good results, in the sense that the peak structure of LDOS resembles spectra obtained in the STM $`dI/dV`$ measurements. In this work we improve the theoretical description by performing selfconsistent three-dimensional DFT calculations, where the effects of the underlying monolayer and substrate are introduced. The hexagonal QD is modeled by a cylindrical jellium QD, and the underlying Na monolayer and Cu(111) substrate by the two-density jellium slab, as described below. Comparisons between calculations of free-standing QD’s and QD’s on a substrate show indeed that the underlying monolayer and substrate induce a new type of states that the simple ’particle-in-a-box’ calculations cannot account for. Since the $`z`$-dependence of wave-functions is included in our calculations, we can calculate the tunneling current at realistic STM-tip distances above the system, and estimate the energy dependence of the step height from the calculated constant current topographs.

The different model systems studied in this work are shown in Fig. 1. Our model is readily applicable to the case of a free-standing cylindrical quantum dot, where we use zero Dirichlet boundary conditions for the wave functions and for the Coulomb potential. The next step is to model a free-standing monolayer. A uniform planar system cannot be exactly reproduced in the axial symmetry. We adopt an approximation scheme analogous to the Wigner-Seitz methodWigner-Seitz . We imagine the plane being filled by hexagons, and then approximate these hexagons by area-covering circles. In order to sample the Brillouin zone of the lattice of circles we use two $`𝐤`$-points, $`𝐤=0`$ and $`𝐤`$ at the Brillouin zone boundary. The wave functions with $`𝐤=0`$ are required to have a vanishing radial derivative at the radius of the circle whereas the wave functions with $`𝐤`$ at the Brillouin zone boundary vanish there. According to our calculations the model gives a uniform ($`r`$-independent) charge distribution for the monolayer. It also minimizes the interactions between a QD inside a circle with its periodic images makov .

The next step is to place the Na monolayer on top of the Cu(111) substrate. The effect of the substrate is modeled using the two-jellium model which is illustrated in Fig. 1c and Fig. 2. We do not model the electrons of the bulk Cu. The density of electrons per unit area in the two-jellium model is kept the same as in the jellium model for a free-standing monolayer. We add a layer of lower density jellium, in order to mimic the different wave function decays into the substrate and into the vacuum. The thickness $`w_2`$ and density (via $`r_{s2}`$) give two free parameters of the lower density jellium, which we adjust in order to reproduce the relevant experimental values of the first and second surface band bottoms at the coverages of 1 ML and 2 ML, respectively. Here the first and second bands correspond to wave functions with one and two nodes in the vertical direction, respectively. The thickness $`w_1^{}`$ of the higher density jellium in the two-jellium model is given by

$$w_1^{}=w_1\left(\frac{r_{s1}}{r_{s2}}\right)^3w_2,$$
(7)
where $`w_1`$ is the thickness of the free-standing Na monolayer.

### IV.1 The underlying monolayer and substrate

To test our model, we first study the systems of a free-standing Na monolayer and that of a Na monolayer on Cu(111), and compare the results with other theoretical results and experimental findings. The Na jellium density is determined from the bulk nearest neighbor distance of $`6.92\mathrm{a}_0`$ and the experimental height of $`5.5\mathrm{a}_0`$ ($`2.9`$Å) of 1 ML of Na on Cu(111) Kliewer\_Berndt\_PRB01 ; Berndt . The resulting density parameter, $`r_s=3.79\mathrm{a}_0`$, gives a slightly higher density than its bulk value of $`3.93\mathrm{a}_0`$ for Na. The thickness and the density of the lower density slab have been chosen by fitting the bottom of the second band for the 1 ML Na coverage and that of the third band for the 2 ML Na coverage on Cu(111) to the experimental valuesFisher ; Lindgren-Wallden88 ; Carlsson94 . The values of $`r_{s2}=6.0\mathrm{a}_0`$ and $`w_2=6.3\mathrm{a}_0`$ give (using the unit cell of radius $`72.8\mathrm{a}_0`$ containing 400 electrons per monolayer) in the 1 ML case the bottom of the second band at $`75\mathrm{meV}`$ below the Fermi level and in the 2 ML case the bottom of the third band at $`50\mathrm{meV}`$ above the Fermi level. These values are reasonably close to the experimental values of about $`100\mathrm{meV}`$ below and above the Fermi level, respectively Lindgren-Wallden88 ; Carlsson94 ; Carlsson97 . The correct positions relative to the Fermi level are important, because we solve for the electronic structures selfconsistently, so that the occupancies of the single-electron states affect the potential and the character of the states themselves.

### IV.2 The quantum dot

We start by studying a free-standing monolayer-thick Na QD (Fig. 1a), since this system shows close resemblance to the simple particle-in-a-box system often used as a first approximation when describing the electronic structure of a QD on a surface. The number of atoms in the QD is chosen to 100, which corresponds to a QD radius of about $`R=36.40\mathrm{a}_0`$. The uppermost panel in Fig. 3 shows the corresponding energy spectrum relative to the Fermi energy $`E_F`$. Note that the the discrete energy eigenvalues are plotted as a function of the quantum number $`m`$ and not as a function of $`𝐤`$. There are three bands below the vacuum level, but only the first band (no horizontal nodal planes) is occupied. The emergence of the succeeding second and third bands can be seen as the condensation of the energy levels at around 1.1 eV and 2.5 eV, respectively.

The LDOS calculated at the cylinder axis at the jellium edge and at $`8\mathrm{a}_0`$ above the edge are shown in the lowest and the middle panel of Fig. 3, respectively. Only states with $`m=0`$ contribute, since they are the only ones with nonzero contributions at the axis. The discrete energy levels are broadened to Lorenzians with the width $`\mathrm{\Gamma }=8\mathrm{meV}`$. The LDOS at the jellium edge can easily be resolved in terms of the contributions from the different bands: The peaks corresponding to first, second, and third bands form series with quadratically increasing intervals and smoothly increasing peak amplitudes. At the distance of $`8\mathrm{a}_0`$ above the QD, the contribution due to the first band states is diminished and the contribution due to the third band states with high quantization in the z-direction (two horizontal nodal planes) is dominating the LDOS. Comparison with the results of the simple ’particle-in-a-box’ calculation by Lindberg and Hellsing vanja1  shows that these two calculations give qualitatively the same results.

We now compare these results for the free-standing QD to the system of the QD adsorbed on a Na monolayer on the Cu(111) surface (Fig. 1d). In our calculation, the substrate is a two-jellium cylindrical supercell containing 400 electrons. The energy spectrum is shown in the uppermost panel of Fig. 4 in which the different $`m`$ levels corresponding to the $`𝐤=0`$ -points are given with the $`𝐤`$ dispersion calculated using the $`𝐤`$ point at the Brillouin zone boundary. In comparison with the spectrum of the monolayer-thick QD in Fig. 3 the bands are shifted downwards, because of the larger jellium thickness at the QD. The lowest energy states (in the bulb of the level diagram) have no $`𝐤`$-dispersion and they are localized at the QD and the substrate slab below it (See Fig. 5). Their dispersion as a function of $`m`$ is similar to that in Fig. 3. Introducing the underlying substrate gives rise to a new type of states, which are not localized to the QD region. In fact, these states form the overwhelming majority. As a result, new bands are induced in the energy spectrum and they are less dispersive as a function of $`m`$. The reduced $`m`$ dispersion reflects the fact that the states are extended over the entire circular supercell. For the same reason, these bands have a larger dispersion in the $`𝐤`$-space than the localized QD bands. A simple particle-in-a-box or free-standing QD calculation cannot provide states of this kind, and it is therefore interesting to see to what extent they contribute to the local electronic structure above the QD.

The LDOS for the QD adsorbed on a Na monolayer on the Cu(111) surface is given in the lowest and the middle panel of Fig. 4 at the jellium edge and at $`14\mathrm{a}_0`$ above the edge, respectively. In order to avoid complications due to the interactions between the supercells we calculate the LDOS using only the $`𝐤=0`$ states and the LDOS is then calculated as in the case of the free-standing monolayer-thick Na QD. To enable a thorough comparison with the free-standing QD results we have to study first the wave functions in more details.

Fig. 5 gives all the wave functions in the interesting energy region for the $`m=0`$ states of the QD adsorbed on the Na monolayer. The first three states, $`n`$ = 1, 2, and 3, correspond to states within the first band. They are localized to the QD and the substrate slab below it and they have no nodes in the z-direction. The $`n`$ = 4 state shows another character with a density no longer localized to the QD region but spread also to the slab region around. This is the first state belonging to the new type of bands, induced by the slab. State $`n`$ = 8 is the beginning of the next band consisting of states with one node in the z-direction (second band in the QD). These states are resonance states, the amplitude of which is strongly enhanced in the QD region, but due to the hybridization with the delocalized slab states they are actually delocalized to the whole system. The $`n`$ = 14 state starts the next band consisting of delocalized states with one horizontal node in the slab region, i.e. it is a second band in the slab. Finally, state $`n`$ = 16 represents the first resonance state with two nodes in the z-direction (third band in the QD). In Fig. 5 one notes that the states $`n`$ = 16 and 17 and also the states $`n`$ = 19 and 20 form pairs. The state lower in energy in the pair is a bonding combination of a QD state and a surrounding slab state whereas the state higher in energy is an antibonding combination.

The electronic structure of the QD adsorbed on the Na monolayer, discussed above in terms of the $`m=0`$ wave functions, is reflected in the LDOS in Fig. 4. First, in the LDOS at the jellium edge we notice that, in comparison with the free-standing QD model, the underlying slab both introduces new type of bands and squeezes the other bands more tightly together to fit more states below the Fermi level (the lowest panels of Figs. 3 and 4). Therefore the present LDOS looks qualitatively different from that for the monolayer-thick unsupported QD in Fig. 3. Moreover, the hybridization of the QD and surrounding slab states to bonding-antibonding pairs causes the splitting of peaks seen clearly for the third band states at the distance of $`14\mathrm{a}_0`$ above the jellium edge (the middle panel of Fig. 4). If we had a continuous spectrum of slab states we would have a single resonance peak with a finite energy width.

Having the modeling of the STM results in mind, the interesting question arising is whether or not the localized states calculated by the free-standing QD model give the same LDOS far above the QD as the states in the model including the substrate slab. Studying the LDOS plots, matching each peak with the corresponding wave function, we notice that the resonance states with strongly enhanced amplitude in the QD region are dominant at large distances above the QD. The contribution of the more delocalized slab states is small. Therefore the free-standing QD model is expected to preserve validity in predicting LDOS at large distances above the QD. The too broad energy spectrum in the free-standing QD model can be corrected for by increasing the dot height with a monolayer of Na jellium or with a two-jellium layer.

## V Comparison with experiment

One of the most useful instruments in surface science is the Scanning Tunneling Microscope (STM) binnig ; tersoff . It can be used for measuring the real-space electronic distribution with atomic resolution, as well as the local energy distribution of electrons and the lifetimes of excited electron states. The real-space distribution is achieved when scanning the surface either in the constant current mode, where the tunneling current is kept constant by changing the tip-surface distance using a feedback mechanism, or in the constant height mode, where the tunneling current is measured when scanning the surface at a constant tip height. The resulting image then displays the topography of the surface. Information about the local electronic structure is obtained by measuring the current variation with the applied voltage. This quantity, the differential conductance $`dI/dV`$, is proportional to the product of the local density of states (LDOS) and the transmission coefficient $`T`$selloni ; LiPRB97 . However, if the applied voltage is small, the bias dependence of $`T`$ is small, and with Eq. 3 we approximately have

$$\frac{dI}{dV}\underset{m\mathrm{𝐤𝐧}}{}(2\delta _{0m})|U_{m\mathrm{𝐤𝐧}}(𝐫)|^2\delta (ϵ_{m\mathrm{𝐤𝐧}}eV).$$
(8)
In the STM study by Kliewer and Berndt Berndt , constant current topographs and $`dI/dV`$ measurements are presented for a Na island on Na monolayer on Cu(111). The size of the island is $`230\times 170\mathrm{a}_{0}^{}{}_{}{}^{2}`$ ($`120\times 90`$Å<sup>2</sup>). We have studied a cylindrical jellium dot with similar dimensions, i.e. having the radius of $`85\mathrm{a}_0`$ and containing thus about 550 electrons. The Na/Cu substrate is described in our calculations by a cylindrical two-density-jellium supercell with the radius of $`160\mathrm{a}_0`$ and containing 2000 electrons. The radius of $`85\mathrm{a}_0`$ is actually fixed to reproduce the peak structure of $`dI/dV`$ spectra by Kliewer and Berndt Berndt  as well as possible (See Fig. 6 and discussion below). For comparison, Kliewer and Berndt Berndt  used in their modeling two-dimensional hard-wall hexagons with the radius of $`83\mathrm{a}_0`$.

We show in Fig. 6 the LDOS at $`18\mathrm{a}_0`$ above the jellium edge both at the $`z`$-axis of the QD and at $`r=20\mathrm{a}_0`$ away from the axis. The height corresponds to a typical tip-sample distance in the STM experiments. The LDOS is calculated as for the smaller systems in Sec. IV.2. The peaks in the figure correspond to states with two horizontal nodes in the 2 ML part of the system (third band in the QD). At the axis, only the $`m=0`$ states contribute, while away from the axis also peaks with $`m0`$ occur. The LDOS peaks can be labeled with the ”quantum number” $`N`$ by counting for the number of radial nodes of the corresponding wave functions in the 2 ML part (See Fig. 6). For $`m`$=0, the $`N=1`$ state has no radial nodes in the dot region whereas $`N=2`$ has one radial node and so on for larger $`N`$. As in the case of the smaller dot discussed above in Sec. IV.2, the states strongly peaked in the QD are resonance states due to the hybridization with the states of the surrounding monolayer and span the whole system. Besides the delocalization of the states, the resonance character causes the fact that in the LDOS (Fig. 6) several peaks may correspond to the same resonance state as discussed earlier for the smaller system. We have identified the LDOS peaks by examining the wave functions. The horizontal lines below the quantum numbers $`m`$ and $`N`$ connect the peaks belonging to the resonance in question.

The relative positions of peaks appearing in the experimental $`dI/dV`$ spectra by Kliewer and Berndt Berndt  are shown in Fig. 6 as arrows pointing downwards. The experimental data is shifted so that the lowest experimental peak coincides with the lowest calculated peak. The experimental spectrum is recorded slightly off from the center of the hexagonal QD which should be taken into account when comparing with the calculated results. Previous two-dimensional free-electron calculations for hexagonal potential boxes have reproduced well the experimental peak positions Berndt ; CrampinSS . In the present modeling, the experimental peak positions agree with the calculated $`m=0`$ resonance positions with the exception of the third and fifth experimental peak. The calculated $`m=0`$ resonances obey the pattern $`\epsilon _n=E_0+AN^2`$, as would be expected for a free particle in a hard-wall cylinder. We have fitted $`E_022\mathrm{meV}`$, $`A15\mathrm{meV}`$. It is gratifying to note that in the LDOS recorded off the cylinder axis strong $`m=1`$ resonance peaks appear so that the third and fifth experimental peaks can be explained. Thus, our model can reproduce quantitatively the experimental peak positions. According to our calculations the resonance width increases toward to higher energies. The increase is maybe slightly stronger than in experiment, indicating a somewhat too weak confinement of the resonance states in our model.

We have also calculated the isosurfaces of the LDOS at the energies corresponding to the dominant peaks in Fig. 6. The results are shown along with the total electron density in Fig. 7. In order to see clearly the nodal structures of the different states, the LDOS is calculated using a smaller lorenzian width of 0.8 meV. The density is smooth in the interior of the QD and shows minor oscillations at the perimeter of the QD. The development of the nodal structure is clear and compares qualitatively well with that found in the experimental $`dI/dV`$ maps Kliewer\_Berndt\_PRB01 . It can also be seen that the isosurfaces corresponding to the two $`(m,N)=(0,4)`$ peaks differ from each other mainly near the perimeter of the dot. The higher peak shown at the energy corresponding to the highest (0,3) state, on the other hand, does not show equally clear (0,3)-character. The explanation is a state with quantum numbers (2,2) at almost exactly the same energy.

From Eq. 8 we obtain a simple formula for the tunneling current,

$$I(U,r,z)_{E_F}^{E_F+U}\mathrm{LDOS}(E,r,z)𝑑E.$$
(9)
From these numerical constant current topographs, we can estimate the apparent step height at the perimeter of the QD. The inset of Fig. 8 shows the constant LDOS height as a function of the distance from the cylinder axis for the bias voltages of -400, 0, and +400 meV. The value of the LDOS is $`10^{11}`$ arb. units on the scale of Fig. 6. The absolute height of the isosurface from the jellium edge depends naturally on the LDOS value chosen, but, according to our calculations, the relative changes are insensitive to the LDOS value over a wide range of values. In order to construct the apparent step height we first obtain the numerical constant current topographs (Eq. 9), then average the profiles over the oscillations above the 2 ML and 1 ML parts of the system and take the difference. The results are shown in Fig. 8 as a function of the bias voltage. The trends seen can be explained by studying the LDOS-isosurfaces (see the inset of Fig. 8), and then noting that the constant current surfaces are obtained by simple integation (Eq. 9). At -100 meV, the second band starts to contribute in the 1 ML part of the system rising the height there and thereby lowering the step height. Then the onset of the third band in the 2 ML part rises the step again. This rise is similar to that seen in the experiment by Kliewer and BerndtKliewer\_Berndt\_PRB01  as well as the decline at higher voltages. However, the comparison with the experimental result shows differences: our step height is too low by a factor of 2, and the raising of the step at negative bias voltages is not seen in the experiment. There may be several reasons for the differences in the step heights. One is that the experimental step height of $`5.5\mathrm{a}_0`$ ($`2.9`$Å), which is determined at a voltage just before the rise in the step height, is directly used as the thickness of the jellium describing the second monolayer of Na. A more consistent procedure might be to take the voltage dependence into account. Moreover, the apparent step height depends on the relative vacuum decay rates of the second and third band states of the 1 ML and 2 ML systems, respectively. Their correct description may be too demanding for our simple model.

## VI Conclusions

In this paper we have presented a model for the electronic structures of alkali metal islands or quantum dots adsorbed on metal surfaces. In particular, we have focused on the system of Na on the Cu(111) surface, where approximately hexagonal Na quantum dots have been observed to form during the epitaxial growth of the second Na monolayer.

We have modeled the quantum dots as small cylindrical jellium islands, and the underlying Na monolayer and Cu substrate as a two-density jellium slab. The parameters of the model have been chosen to fit experimental spectroscopic data and calculated first-principles band structures for one and two completed monolayers of Na on the Cu(111) surface. The calculations were performed in the context of the density-functional theory, using a real-space electronic structure calculation method.

The calculated results are compared with experimental findings from scanning tunneling microscope and photoemission experiments. The model gives local densities of states which are in a quantitative agreement with constant current topographs and $`dI/dV`$ spectra and maps. Thereby the idea of surface states which are localized as resonances at the quantum dots is supported. The future applications of the model will include studies of the adsorption and dissociation of molecules in the vicinity of alkali metal quantum dots.

###### Acknowledgements.

We would like to thank Lars Wallden and S.Å. Lindgren for sharing their knowledge on the system based on PES experiments. We would like to thank R. M. Nieminen for many useful discussions. We acknowledge the generous computer resources from the Center for Scientific Computing, Espoo, Finland. One of the authors (T.T.) acknowledges financial support by the Vilho, Yrjö and Kalle Väisälä foundation. This research has been supported by the Academy of Finland through its Centers of Excellence Program (2000-2005).


format: text

text:
# Vortex line representation for flows of ideal and viscous fluids []

## Abstract

It is shown that the Euler hydrodynamics for vortical flows of an ideal fluid coincides with the equations of motion of a charged compressible fluid moving due to a self-consistent electromagnetic field. Transition to the Lagrangian description in a new hydrodynamics is equivalent for the original Euler equations to the mixed Lagrangian-Eulerian description - the vortex line representation (VLR) . Due to compressibility of a ”new” fluid the collapse of vortex lines can happen as the result of breaking (or overturning) of vortex lines. It is found that the Navier-Stokes equation in the vortex line representation can be reduced to the equation of the diffusive type for the Cauchy invariant with the diffusion tensor given by the metric of the VLR.

PACS: 47.15.Ki, 47.32.Cc

1. Collapse as a process of a singularity formation in a finite time from the initially smooth distribution plays the very important role being considered as one of the most effective mechanisms of the energy dissipation. For hydrodynamics of incompressible fluids collapse must play also a very essential role. It is well known that appearance of singularity in gasodynamics, i.e., in compressible hydrodynamics, is connected with the phenomenon of breaking that is the physical mechanism leading to emergence of shocks. From the point of view of the classical catastrophe theory  this process is nothing more than the formation of folds. It is completely characterized by the mapping corresponding to transition from the Eulerian description to the Lagrangian one. Vanishing the Jacobian $`J`$ of this mapping means emergence of a singularity for spatial derivatives of velocity and density of a gas. In the incompressible case breaking as intersection of trajectories of Lagrangian particles is absent because the Jacobian of the corresponding mapping is fixed, in the simplest case equal to unity. By this reason, it would seem that there were no any reasons for existence of such phenomenon at all. In spite of this fact, as it was shown in , breaking, however, is possible in this case also. It can happen with vortex lines. Unlike the breaking in gasodynamics, the breaking of vortex lines means that one vortex line reaches another vortex line. For smooth initial conditions breaking happens first time while touching vortex lines at a single point. In the touching point the vorticity becomes infinite. And this is possible in spite of incompressibility of both divergence-free fields, i.e., vorticity and velocity. To describe the breaking of vortex lines in the papers  it was suggested the vortex line representation – a mixed Lagrangian-Eulerian description when each vortex line is labeled by a two-dimensional marker and another parameter defines the vortex line itself.

This paper is devoted to development of this method to apply to both ideal and viscous fluids. We clarify the role of the Clebsch variables in the vortex line representation: these variables can be used as Lagrangian markers of vortex lines. However, as well known, these variables can be introduced always only locally and, generally speaking, can not be extended for the whole space. In the general situation we demonstrate in this paper that transition to the vortex line representation is equivalent to consideration of a new compressible hydrodynamics of a charged fluid flowing under action of a self-consistent electromagnetic field. In this case the electric and magnetic fields satisfy the Maxwell equations. The most essential property of a new hydrodynamics is a compressibility of a new fluid that for its Lagrangian description means compressibility of the corresponding mapping and, respectively, a possibility of a breaking. In terms of the Eulerian characteristics this results in the breaking of vortex lines when the vorticity $`𝛀=\text{curl}𝐯`$ takes infinite value. In the framework of the new hydrodynamics of a charged fluid the role of density plays the quantity inverse to $`J`$ which is naturally called as a density of vortex lines. This quantity appears from the Cauchy formula for the vorticity $`𝛀`$. Evolution of the vortex line density in time and space is defined by the velocity component normal to the vorticity. As it is shown in this paper the Cauchy formula can be obtained from a “new” Kelvin theorem as well as from the analog of the Weber transformation. As the result, the Euler equations turn out to be resolved with respect to the Cauchy invariants, i.e., relative to the infinite number of integrals of motion. In this case one can consider the Euler equations as the partially integrated equations. This circumstance is very important for numerical solution of the Euler equation.

The vortex line representation can be applied not only to ideal hydrodynamics but also to flow description of viscous incompressible fluids in the framework of the Navier-Stokes equation. In the paper we obtain the equation of the diffusion type describing dynamics of the Cauchy invariant in the viscous case with the “diffusion tensor” determined by the VLR metric. In its form this equation coincides with the equation derived in . In this case the equations of motion of vortex lines in its original (for ideal fluids) form are understood as the equations given the transformation to a new curvilinear system of coordinates. The obtained exact equations for description of viscous flows can be considered as the result of exact separation of two different temporal scales: the inertial (in fact, nonlinear) scale and the viscous one.

2. As well known ( see, for instance, , ) the Euler equations for an ideal incompressible fluid,

$$\frac{𝐯}{t}+(𝐯)𝐯=p,\text{div}𝐯=0,$$
(1)
in both two-dimensional and three-dimensional cases possess the infinite (continuous) number of integrals of motion. These are the so called Cauchy invariants. The most simple way to derive the Cauchy invariants is one to use the Kelvin theorem about conservation of the velocity circulation,

$$\mathrm{\Gamma }=(𝐯d𝐥),$$
(2)
where the integration contour $`C[𝐫(t)]`$ moves together with a fluid. If in this expression one makes a transform from the Eulerian coordinate $`𝐫`$ to the Lagrangian ones $`𝐚`$ then Eq. (2) can be rewritten as follows:

$$\mathrm{\Gamma }=\dot{x}_i\frac{x_i}{a_k}𝑑a_k,$$

where a new contour $`C[𝐚]`$ is already immovable. Hence, due to arbitrariness of the contour $`C[𝐚]`$ and using the Stokes formula one can conclude that the quantity

$$𝐈=\text{rot}_a\left(\dot{x}_i\frac{x_i}{𝐚}\right)$$
(3)
conserves in time at each point $`𝐚`$. This is just the Cauchy invariant. If the Lagrangian coordinates $`𝐚`$ in (3) coincide with the initial positions of fluid particles the invariant $`𝐈`$ is equal to the initial vorticity $`𝛀_0(𝐚)`$.

Conservation of these invariants, as it was shown first by Salmon , is consequence of the special (infinite) symmetry - the so-called relabeling symmetry. The Cauchy invariants characterize the frozenness of the vorticity into fluid. This is a very important property according to which fluid (Lagrangian) particles can not leave its own vortex line where they were initially. Thus, the Lagrangian particles have one independent degree of freedom – motion along vortex line. From another side, such a motion as it follows from the equation for the vorticity

$$\frac{𝛀}{t}=\text{rot}[𝐯\times 𝛀],$$
(4)
does not change its value. From this point of view a vortex line represents the invariant object and therefore it is natural to seek for such a transformation when this invariance is seen from the very beginning. Such type of description - the vortex line representation - was introduced in the papers  by Ruban and the author of this paper.

3. Consider the vortical flow ($`𝛀0`$) of an ideal fluid given by the Clebsch variables $`\lambda `$ and $`\mu `$:

$$𝛀=[\lambda \times \mu ].$$
(5)
The geometrical meaning of these variables is well known: intersection of two surfaces $`\lambda =\text{const}`$ and $`\mu =\text{const}`$ yields the vortex line. It is known also that in the incompressible case the Clebsch variables are Lagrangian invariants, being unchanged along trajectories of fluid particles:

$$\frac{\lambda }{t}+(𝐯)\lambda =0;\frac{\mu }{t}+(𝐯)\mu =0.$$
(6)
Therefore these variables can be taken as markers for vortex lines. It is easily to establish that transition in (5) to new variables

$$\lambda =\lambda (x,y,z),\mu =\mu (x,y,z),s=s(x,y,z),$$
(7)
where $`s`$ is the parameter given the vortex line, leads to the expression

$$𝛀(𝐫,t)=\frac{1}{J}\frac{𝐑}{s}$$
(8)
where

$$J=\frac{(x,y,z)}{(\lambda ,\mu ,s)}$$
(9)
is the Jacobian of the mapping

$$𝐫=𝐑(\lambda ,\mu ,s).$$
(10)
The transform (10) inverse to (7) defines the corresponding transition to the curvilinear, connected with vortex lines, system of coordinates.

The equations of motion of vortex lines - the equations for $`𝐑(\lambda ,\mu ,s,t)`$ – can be obtained directly from the equation of motion for the vorticity (4). However, the most simple way to derive them is to use the combination of the equations (6):

$$\mu \left[\frac{\lambda }{t}+(𝐯)\lambda \right]\lambda \left[\frac{\mu }{t}+(𝐯)\mu \right]=0,$$
(11)
which is identical to (6) due to a linear independence of the vectors $`\lambda `$ and $`\mu `$.

Performing in (11) the transformations (7), we arrive at the equation of motion for vortex lines :

$$\left[\frac{𝐑}{s}\times \left(\frac{𝐑}{t}𝐯(𝐑,t)\right)\right]=0.$$
(12)
This equation has one important property: any motion along a vortex line does not change the line itself. It is easily to check that Eq. (12) is equivalent to the equation

$$\frac{𝐑}{t}=𝐯_n(𝐑,t),$$
(13)
where $`𝐯_n`$ is the velocity component normal to the vorticity vector.

In accordance with the Darboux theorem, the Clebsch variables can be introduced locally always but not globally. It is well known also that the flows parameterized by the Clebsch variables has a zero helicity integral $`(𝐯\text{rot}𝐯)𝑑𝐫`$ – the topological invariant which characterizes a degree of knottiness of vortex lines. Therefore to introduce the vortex line representation for flows with nontrivial topology it is necessary to come back to the original equations of motion (1) and (4) for velocity and vorticity.

4. According to the equation (4) the tangent to the vector $`𝛀`$ velocity component $`𝐯_\tau `$ does not effect (directly) on the vorticity dynamics, i.e., in (4) we can put, instead of $`𝐯`$, its transverse component $`𝐯_n`$.

The equation of motion for the transverse velocity $`𝐯_n`$ follows directly from the equation (1). It has the form of the equation of motion of charged particle moving in an electromagnetic field:

$$\frac{𝐯_n}{t}+(𝐯_n)𝐯_n=𝐄+[𝐯_n\times 𝐇],$$
(14)
where the effective electric and magnetic fields are given by the expressions:

$$𝐄=\left(p+\frac{v_\tau ^2}{2}\right)\frac{𝐯_\tau }{t},$$
(15)
$$𝐇=\text{rot}𝐯_\tau .$$
(16)
Interesting to note that the electric and magnetic fields introduced above are expressed through the scalar $`\phi `$ and vector $`𝐀`$ potentials by the standard way:

$$\phi =p+\frac{𝐯_\tau ^2}{2},𝐀=𝐯_\tau ,$$
(17)
so that two Maxwell equations

$$\text{div}𝐇=0,\frac{𝐇}{t}=\text{rot}𝐄$$

satisfy automatically. In this case the vector potential $`𝐀`$ has the gauge

$$\text{div}𝐀=\text{div}𝐯_n,$$

which is equivalent to the condition $`\text{div}𝐯=0`$.

Two other Maxwell equations can be written also but they can be considered as definition of the charge density $`\rho `$ and the current $`𝐣`$ which follow from the relations (15) and (16). The basic equation in the new hydrodynamics is the equation of motion (14) for the normal component of the velocity which represents the equation of motion for nonrelativistic particle with a charge and a mass equal to unity, the light velocity in this units is equal to $`1`$.

The equation of motion (14) is written in the Eulerian representation. To transfer to its Lagrangian formulation one needs to consider the equations for ”trajectories” given by the velocity $`𝐯_n`$:

$$\frac{d𝐑}{dt}=𝐯_n(𝐑,t)$$
(18)
with initial conditions

$$𝐑|_{t=0}=𝐚.$$

Solution of the equation (18) yields the mapping

$$𝐫=𝐑(𝐚,t),$$
(19)
which defines transition from the Eulerian description to a new Lagrangian one.

The equations of motion in new variables are the Hamilton equations:

$$\dot{𝐏}=\frac{h}{𝐑},\dot{𝐑}=\frac{h}{𝐏},$$
(20)
where dot means differentiation with respect to time for fixed $`𝐚`$, $`𝐏=𝐯_n+𝐀𝐯`$ is the generalized momentum, and the Hamiltonian of a particle $`h`$ being a function of momentum $`𝐏`$ and coordinate $`𝐑`$ is given by the standard expression:

$$h=\frac{1}{2}(𝐏𝐀)^2+\phi p+\frac{𝐯^2}{2},$$

i.e., coincides with the Bernoulli ”invariant”.

The first equation of the system (20) is the equation of motion (14), written in terms of $`𝐚`$ and $`t`$, and the second equation coincides with (18).

For new hydrodynamics (14) or for its Hamilton version (20) it is possible to formulate a ”new” Kelvin theorem (it is also the Liouville theorem):

$$\mathrm{\Gamma }=(𝐏d𝐑),$$
(21)
where integration is taken along a loop moving together with the ”fluid”. Hence, analogously as it was made before while derivation of (3) we get the expression for a new Cauchy invariant:

$$𝐈=\text{rot}_a\left(P_i\frac{x_i}{𝐚}\right).$$
(22)
Its difference from the original Cauchy invariant (3) consists in that in the equation of motion (18) instead of the velocity $`𝐯`$ stands its normal component $`𝐯_n`$. As consequence, the ”new” hydrodynamics becomes compressible: $`\text{div}𝐯_n0`$. Therefore on the Jacobian $`J`$ of the mapping (19) there are imposed no restrictions. The Jacobian $`J`$ can take arbitrary values.

From the formula (22) it is easily to get the expression for the vorticity $`𝛀`$ in the given point $`𝐫`$ at the instant $`t`$ (compare with ):

$$𝛀(𝐫,t)=\frac{(𝛀(𝐚)_a)𝐑(a,t)}{J},$$
(23)
where $`J`$ is the Jacobian of the mapping (19) equal to

$$J=\frac{(x_1,x_2,x_3)}{(a_1,a_2,a_3)}.$$

Here we took into account that the generalized momentum $`𝐏`$ coincides with the velocity $`𝐯`$, including the moment of time $`t=0`$: $`𝐏_0(𝐚)𝐯_0(𝐚)`$. $`𝛀_0(𝐚)`$ in this relation is the ”new” Cauchy invariant with zero divergence: $`\text{ }\mathrm{div}_a𝛀_0(a)=0`$.

The representation (23) generalizes the relation (5) to an arbitrary topology of vortex lines. The variables $`𝐚`$ in this expression can be considered locally as a set of $`\lambda `$, $`\mu `$ and $`s`$.

As known (see, for instance, ), expression for the Cauchy invariant can be obtained from the Weber transformation. This is the representation of velocity in terms of the initial data which can be obtained by integrating the Cauchy invariant (23).

Consider the following one-form $`\omega =(𝐏d𝐑)`$ and calculate its time derivative. By means of the equations of motion (20) we get:

$$\dot{\omega }=d[h+(𝐏\dot{𝐑})].$$

Hence it follows that the vector function

$$u_k=\frac{x_i}{a_k}P_i,$$

dependent on $`t`$ and $`𝐚`$, will obey the following equation of motion:

$$\dot{u}_k=\frac{}{a_k}\left(p+\frac{v_n^2}{2}\frac{v_\tau ^2}{2}\right).$$

Integration of this equation in time gives the Weber-type transformation:

$$u_k(𝐚,t)=u_{k0}(𝐚)+\frac{\mathrm{\Phi }}{a_k},$$
(24)
where the potential $`\mathrm{\Phi }`$ satisfies the nonstationary Bernoulli equation:

$$\dot{\mathrm{\Phi }}=p+\frac{v_n^2}{2}\frac{v_\tau ^2}{2}.$$

If $`\mathrm{\Phi }|_{t=0}=0`$ then the time independent vector $`𝐮_0(𝐚)`$ coincides with the initial velocity $`𝐯_0(𝐚)`$. By applying the operator curl to the relation (24) we arrive again at the Cauchy invariant (22).

Thus, in the general situation the equation of motion of vortex lines has the form (18) which is completed by the relation (23) and the equation

$$𝛀(𝐫,t)=\text{rot}_r𝐯(𝐫,t)$$
(25)
with additional constraint $`\text{div}_r𝐯(𝐫,t)=0`$.

The equations of motion (18), (25) together with the relation (23) can be considered as the result of partial integration of the Euler equation (1). These new equations are resolved with respect to the Cauchy invariants – an infinite number of integrals of motion, that is a very important issue for numerical integration (see ). For the partially integrated system the Cauchy invariants conserve automatically that, however, for direct numerical integration of the Euler equation one needs to test in which extent these invariants remain constant. Probably, this is one of the main restrictions defining accuracy of discrete algorithms of direct integration of the Euler equations.

Another very important property of the vortex line representation is absence of any restrictions on the value of the Jacobian $`J`$ which do exist, for instance, for transition from the Eulerian description to the Lagrangian one in the original Euler equation (1) when Jacobian in the simplest situation is equal to unity. The value $`1/J`$ for the system (18), (25), (23) has a meaning of a density $`n`$ of vortex lines. This quantity as a function of $`𝐫`$ and $`t`$ , according to (18, obeys the discontinuity equation:

$$\frac{n}{t}+\text{div}_r(n𝐯_n)=0.$$
(26)
In this equation $`\text{div}_r𝐯_n0`$ because only the total velocity has zero divergence.

5. Consider now the question about application of the VLR to flows of viscous fluids. Write down the Navier-Stokes equation for vorticity $`𝛀`$:

$$\frac{𝛀}{t}=\text{rot}[𝐯\times 𝛀]\nu \text{rot}\text{rot}𝛀,$$
(27)
and perform in this equation the transform to new variables $`𝐚`$ and $`t`$ by means of changes defined by the equation (18) together with the Cauchy relation (23) where $`𝛀_0`$ is assumed a function of not only $`𝐚`$ but also time $`t`$: $`𝛀_0=𝛀_0(𝐚,t)`$.

Then after substitution (23) into (27) the first term in the right hand side is cancelled because of (18). At the result, the equation (27) is written in the form:

$$\frac{1}{J}\left(\frac{𝛀_0}{t}_a\right)𝐑=\nu \text{rot}\text{rot}\left\{\frac{1}{J}(𝛀_0_a)𝐑\right\}.$$
(28)
Next, change differentiation relative to $`𝐫`$ in the r.h.s. of (28) to differentiation against $`𝐚`$. After simple, but cumbersome calculations the equation (28) transforms into the equation for $`𝛀_0(𝐚,t)`$:

$$\frac{𝛀_0}{t}=\nu \text{rot}_a\left(\frac{\widehat{g}}{J}\text{rot}_a\left(\frac{\widehat{g}}{J}𝛀_\mathrm{𝟎}\right)\right).$$
(29)
Formally it is a linear equation for $`𝛀_0`$, here $`\widehat{g}`$ is the metric tensor equal to

$$g_{\alpha \beta }=\frac{x_i}{a_\alpha }\frac{x_i}{a_\beta }.$$

The equation (29) for the Cauchy invariant formaly coincides with that obtained by Zenkovich and Yakubovich for incompressible hydrodynamics  in which the variables $`𝐚`$ are assumed to be Lagrangian markers of fluid particles. In the Zenkovich-Yakubovich equation the Jacobian $`J`$ is proposed to be independent on time, in the simplest case equal to 1. Just this is a principle difference between the Zenkovich-Yakubovich equation and the equation (29). $`J`$ in (29) is a function of time $`t`$ and coordinates $`𝐚`$.

Remarkable peculiarity of the obtained system is the exact separation of two different temporal scales, responsible for the inertial (in fact, nonlinear) processes and for the viscous processes. The former ones are described by the equation (18), and the latter by the equation of diffusive type (29) in which the diffusion ”coefficient”, proportional to viscosity $`\nu `$, is defined by the metric of the mapping $`𝐫=𝐑(𝐚,t)`$.

The author is grateful to E.I. Yakubovich for possibility of acquaintance with the paper  before its publication. This work was supported by the RFBR (grant no. 00-01-00929), by the Program of Support of the Leading Scientific Schools of Russia (grant no. 00-15-96007) and by INTAS (grant no. 00-00797).


text:
# Topological flux sectors in extended U(1) gauge theory on 𝑇⁴

## 1 The flux in abelian gauge theories

Consider an abelian gauge theory defined on a hypercube of size $`L`$ and periodic boundary conditions. In the continuum the definition of flux through any $`(\mu ,\nu )`$ plane is the following:

$$\mathrm{\Phi }_{\mu \nu }=F_{\mu \nu }𝑑\sigma $$
(1)
This quantity, because of periodic boundary conditions, is $`2\pi k`$ valued ($`kZ`$), and configurations with different $`k`$ values are topologically disconnected (so we talk of topological superselection $`sectors`$).

Consider now the same system, but with a discretized space-time. The definition of flux is now

$$\mathrm{\Phi }_{\mu \nu }=\frac{1}{L^2}\underset{\mu \nu planes}{}\underset{P}{}[\theta _{P_{\mu \nu }}]_{\pi ,\pi }$$
(2)
where $`[\theta _P]_{\pi ,\pi }`$ is the plaquette angle reduced to the interval $`(\pi ,\pi )`$, its so called $`physical`$ part. A double sum is present: the internal $`_P`$ is the sum over the plaquettes in a single $`(\mu ,\nu )`$ plane; this quantity is a multiple of $`2\pi `$, as in the continuum case. But when we consider the external average over all parallel planes, we observe that the flux can change from plane to plane, due to the presence of magnetic monopoles, specific to the lattice, so that the allowed values for $`\mathrm{\Phi }_{\mu \nu }`$ are multiples of $`(2\pi )/L^2`$.

Having in mind the continuum limit, configurations whose flux is $`2\pi k`$ play a special role: in fact $`tunneling`$ between different ‘sectors’ is now possible, but higher and higher barriers are expected between them as $`L\mathrm{}`$.

## 2 The model

We work with the following action

$$S=\beta \underset{P}{}\mathrm{cos}\theta _P+\gamma \underset{P}{}\mathrm{cos}2\theta _P$$
(3)
$`\beta `$ is the Wilson coupling. This model was introduced by Bhanot  and has aroused interest also recently . In Fig. 1 we indicate three unknown aspects of the phase diagram with a question mark:
i) the nature of the phase indicated with $`(b)`$;
ii) the nature of the phase $`(d)`$;
iii) the fate of the two phase boundaries in the bottom right corner of the phase diagram (the same holds on the left side, due to the symmetry $`\beta \beta `$).

## 3 Characterization of the phase transition

Let us consider the behavior of flux sectors across a phase boundary; here we work on the Wilson axis, but similar results hold everywhere.
The numerical strategy is very simple: using a $`multicanonical`$ approach, we measure the values of the flux through one $`(\mu ,\nu )`$ orientation along a Monte Carlo simulation and plot the (inverse) histogram of these measurements. The results are the following:
i) In the Coulomb phase flux sectors are well defined (Fig. 2), and as the thermodynamic limit is approached we observe higher and higher barriers between them.
ii) In the confined phase (Fig. 3) the high density of monopoles hides the sectors and gives a Gaussian flux distribution. If we check the thermodynamic limit, the distribution does $`not`$ change: in fact more and more planes are present (this would give a narrower distribution), but also more and more monopoles (which give a broader distribution) and these two effects compensate each other.

## 4 Characterization of the phases

We extend to the abelian context the ideas of ’t Hooft  about non-abelian gauge theories: we probe the response of the system to a variation of the flux. In the context of $`SU(N)`$ pure gauge theories, we know that it is possible to modify the flux (by twisting the boundary conditions) according to the discrete group $`Z_N`$. We can implement a similar modification also in the abelian context, but a continuous variation of the flux is now possible.

Our strategy is the following: we consider a $`stack`$ of plaquettes, one in each plane of a given orientation, and on these we change $`\theta _P\theta _P+\varphi `$, where $`\varphi [0,2\pi ]`$. This corresponds to imposing an extra flux $`\varphi `$ through the chosen orientation. The new partition function is

$$Z(\varphi )=D\theta e^{\beta (_{stack}\mathrm{cos}(\theta _p+\varphi )+_{\overline{stack}}\mathrm{cos}\theta _p)}$$
(4)
and analogously away from the Wilson axis; $`Z(\varphi )`$ is clearly $`2\pi `$ periodic. We measure the free energy of the flux, $`F(\varphi )=\mathrm{log}\frac{Z(\varphi )}{Z(0)}`$.

We make the following observations:
i) In the Coulomb phase we recover the $`2\pi `$ sectors (Fig. 4, upper curve).
ii) In the confined phase (Fig. 4, lower curve), within our statistical error, the free energy is independent of $`\varphi `$. This is due to the decoupling of flux values through different planes, so that the effects of the variation of the flux are screened.
iii) In the unknown phases $`(b),(d)`$ we observe the appearance of an extra $`\pi `$ periodicity (Fig. 5). We present an argument to interpret this result: consider the following quantity

$$P_q(x)P_q^{}(x+L\widehat{e}_\mu )_{\stackrel{~}{Z}}$$
(5)
where $`P_q`$ is the Polyakov loop of charge $`q`$, $`\stackrel{~}{Z}=_0^{2\pi }𝑑\varphi Z(\varphi )`$, and $`\widehat{e}_\mu `$ is the unit vector in direction $`\widehat{\mu }`$. This correlator is non trivial due to ’twisted’ b.c.. Gauss law ($`P_q(x+L\widehat{e}_\mu )=e^{iq\varphi }P_q(x)`$) gives

$$P_q(x)P_q^{}(x+L\widehat{e}_\mu )_{\stackrel{~}{Z}}=\frac{_0^{2\pi }𝑑\varphi e^{iq\varphi }Z(\varphi )}{_0^{2\pi }𝑑\varphi Z(\varphi )}$$
(6)
If $`Z(\varphi )`$ has a $`\pi `$ periodicity, it follows that when $`q=1`$, or more generally odd, the correlator is zero, while if $`q`$ is even, it can be different from zero. We then remember that $`P_qP_q^{}e^{V_{eff}(L)T}`$, so we deduce that odd charges are confined (in even combinations). We therefore claim that these two unknown phases $`(b),(d)`$ are just Coulomb phases for even charges.

## 5 Fate of the phase boundaries

The last question we address regards the fate of the two phase boundaries at the lower corners of Fig. 1. It is possible that they meet at some point or that they just become asymptotically close to each other. What we can exclude is that they end somewhere, letting the two different Coulomb phases communicate.

Suppose that they communicate somewhere, then we expect the twisted free energy $`F(\varphi )`$ to be the same in both: this leads to a contradiction. In fact $`F(\varphi )`$ would have $`least`$ periodicity $`\pi `$ and $`2\pi `$ at the same time. This is possible if the distribution is flat; but this again implies confinement.

## 6 Conclusions and perspectives

We observe that in the literature about $`U(1)`$ the issue of ergodicity through flux sectors is usually not considered, despite the early observation of ; in the region close to criticality this phenomenon could be relevant.

Possible perspectives of this work are the following: the signal related to the different periodicities of the partition function is very strong; it is clearly visible on small lattices also at large $`\beta `$. We can use it as a tool to study the position of the phase boundaries, and so to extend the quantitative picture of the phase diagram.
Secondly, we want to try an FSS analysis of the quantity $`F(\pi )`$ (the free energy with twist-flux $`\mathrm{\Phi }=\pi `$) as we cross a phase boundary: indications about the nature of the phase transition ($`1^{st}`$ or $`2^{nd}`$ order) could be obtained.

We gratefully acknowledge Oliver Jahn for useful discussions.


text:
# 1 Test of Eq. () for a series of “twin experiments” in the 4d Gaussian Ising spin glass strictly following protocols and analysis of Ref. []. 56 shift experiments are presented with 𝑇_𝑖/𝑇_𝑐∈[0.4,0.9], Δ𝑇/𝑇_𝑐=0.05, 0.1, 0.2, ⋯, 0.5, 𝑡_𝑤∈[80,57797]. System size is 𝐿=25≫ℓ_𝑇(𝑡). ℓ_eff was defined from the maximum of ∂{𝐶(𝑡+𝑡_𝑤,𝑡_𝑤)/∂𝑡}, where 𝐶 is the spin autocorrelation function (instead of TRM) and growth laws ℓ_𝑇(𝑡) taken from Ref. []. The exponent 𝜁 of the fictitious overlap length ℓ_{Δ𝑇}^𝐹 is chosen to collapse all our data.

\[

\]

Comment on ‘Symmetrical Temperature-Chaos effect with Positive and Negative Temperature Shifts in a Spin Glass’.

In a very interesting paper, Jönsson, Yoshino and Nordblad  have shown that the effect of small temperature shifts on the aging behaviour of an Heisenberg-like spin glass could be interpreted quantitatively using the ideas of temperature chaos and overlap length. In this comment, we show that the same analysis can be performed in a case where temperature chaos is known to be irrelevant, weakening the main conclusion of Ref. .

As now well established, aging at low temperatures $`T`$ in spin glasses is associated with the slow growth with time $`t`$ of a coherence length, $`\mathrm{}_T(t)`$. This length can be measured in simulations, but can also be inferred from experimental data, using plausible assumptions, leading to a rather consistent determination of $`\mathrm{}_T(t)`$. “Cumulative aging”  means that the same coherence length grows at different temperatures, although at different rates. In this case, the value of $`\mathrm{}_{T_i}(t_w)`$ after staying a certain time $`t_w`$ at a first temperature $`T_i`$ serves as the ‘initial condition’ for the growth of $`\mathrm{}`$ after a temperature shift.

Rejuvenation effects, on the other hand, demonstrate that cumulative aging cannot be the only story in spin glasses. The “temperature chaos” scenario postulates that typical equilibrium configurations in a spin glass at two temperatures differing by $`\mathrm{\Delta }T`$ are strongly correlated only up to the “overlap length”, $`\mathrm{}_{\mathrm{\Delta }T}`$, beyond which these correlations rapidly decay to zero. From scaling arguments, one expects $`\mathrm{}_{\mathrm{\Delta }T}|\mathrm{\Delta }T|^{1/\zeta }`$, with $`\zeta 1`$ for the Ising spin glass. If this scenario holds, one expects that the initial condition for the growth of $`\mathrm{}`$ after a temperature shift, as encoded by an effective length $`\mathrm{}_{\mathrm{eff}}`$, will obey the following scaling form:

$$\frac{\mathrm{}_{\mathrm{eff}}}{\mathrm{}_{\mathrm{\Delta }T}}=F\left(\frac{\mathrm{}_{T_i}(t_w)}{\mathrm{}_{\mathrm{\Delta }T}}\right),$$
(1)
with $`F(x1)=x`$ (cumulative aging) and $`F(x1)=1`$ (temperature chaos). Rejuvenation, i.e. deviations from cumulative aging, can thus be accounted for by temperature chaos. The verification of Eq. (1) using experimentally determined values of $`\mathrm{}_{T_i}`$ and $`\mathrm{}_{\mathrm{eff}}`$ for different $`\mathrm{\Delta }T`$ and $`t_w`$ is the central result of Ref. , thereby providing strong support for the temperature chaos scenario and a numerical value for the chaos exponent for the AgMn spin-glass, $`1/\zeta 2.6`$. In the above argument, one assumes from the start that rejuvenation is induced by temperature chaos, and finds self-consistent results.

However, other scenarii have been proposed in the literature to explain rejuvenation effects, such as the progressive freezing of smaller and smaller length scale modes . In this respect, we recently demonstrated in a numerical simulation of the 4d Ising spin glass that strong rejuvenation effects can be observed in conditions where the overlap length is independently observed to be much larger than all relevant length scales .

To understand these contradictory results , we have reproduced the protocol of Ref.  in an extensive new series of simulations, and followed the very same steps to determine the length scales $`\mathrm{}_{T_i}(t_w)`$ and $`\mathrm{}_{\mathrm{eff}}`$. Although we do know that all dynamic length scales in our simulations are $`5`$ and that the overlap length is much larger (probably larger than $`20`$), we tried to rescale all our results using Eq. (1). Surprisingly, as shown in Fig. 1, this works very well with the expected value $`\zeta 1`$. Rejuvenation results then in our case in the appearance of a fictitious overlap length, $`\mathrm{}_{\mathrm{\Delta }T}^F`$. This shows that although suggestive, the analysis of Ref.  cannot be viewed as definitive evidence for temperature chaos.

Simulations were performed on OSWELL at the Oxford Supercomputing Center, Oxford University, UK.

Ludovic Berthier1 and Jean-Philippe Bouchaud2

1Theo. Phys., 1 Keble Road, Oxford, OX1 3NP, UK.

2SPEC, Orme des Merisiers, 91191 Gif/Yvette, France.