text:
# Effective Coupling for Open Billiards

## I Introduction

During the last years quantum chaotic scattering was a field of intense research. A great deal of the results obtained was based on the projection operator formalism due to Feshbach, Weidenmüller and others MW69 ; Feshbach ; Dit00 . In this approach the scattering system is decomposed into a closed subsystem described by the internal Hamiltonian $`H_{\mathrm{in}}`$ with discrete bound states $`n=1\mathrm{}N`$ and a continuum of external scattering states labeled by the energy $`E`$ and an index $`\lambda =1\mathrm{}\mathrm{\Lambda }(E)`$ corresponding to different open scattering channels. The coupling between the internal and external subsystems is then incorporated by an operator with matrix elements $`W_{n,\lambda }(E)`$. The S-matrix of the complete system can be expressed in terms of these matrix elements and the Hamiltonian $`H_{\mathrm{in}}`$. This relation can be cast into the form

$$S=\frac{I\mathrm{i}K}{I+\mathrm{i}K}$$
(1)
$$K=\pi W^{}\frac{I}{EH_{\mathrm{in}}}W.$$
(2)
Here, $`S`$ and $`K`$ are energy-dependent square matrices of dimension $`\mathrm{\Lambda }\times \mathrm{\Lambda }`$, and $`W`$ has dimension $`N\times \mathrm{\Lambda }`$. While this setting is very general, the tools developed for the subsequent analysis of the properties of the S-matrix require additional assumptions. In particular, the energy thresholds for the opening of new scattering channels are usually neglected. As a consequence the energy dependence of the coupling matrix $`W(E)`$ can be considered weak and then the S-matrix (1) can be rewritten in terms of an $`N\times N`$ effective non hermitian Hamiltonian $`H_{\mathrm{eff}}`$

$$S=I2\pi \mathrm{i}W^{}\frac{I}{EH_{\mathrm{eff}}}W$$
(3)
$$H_{\mathrm{eff}}=H_{\mathrm{in}}\mathrm{i}\pi WW^{}.$$
(4)
This canonical formalism MW69 ; Feshbach ; Dit00  for expressing the S-matrix is sometimes denoted as the Heidelberg description of scattering, and we use this name to distinguish it from an S-Matrix obtained directly, i. e. without reference to any auxiliary closed system.

When the internal Hamiltonian in (4) describes a chaotic system, it is justified to replace it by a random matrix BGS84 , and by performing an average over the appropriate ensemble a statistical theory for the S-matrix is obtained which allows to calculate quantities of interest such as correlation functions or the distribution of Wigner delay times and resonance poles MPS85 ; LW91 ; FS97 ; GMW98 . It was found that the results of such an approach are to a large extent independent of the detailed structure of the matrix $`W`$, but they do depend on the dimensionless mean coupling strength

$$g=\pi ^2\frac{|W_{n,\lambda }|^2_{n,\lambda }}{D}.$$
(5)
Here, $`D`$ is the mean energy level spacing of the internal subsystem, and the average $`\mathrm{}_{n,\lambda }`$ is taken over the internal states $`n`$ and all open scattering channels $`\lambda `$.

For example, when the coupling constant (5) exceeds the critical value $`g=1`$ and the number of scattering channels is small compared to the total number of states, a counterintuitive shrinking of the widths of most resonances with increasing coupling is observed SZ89 ; SZ92 ; Rot91 . For each attached scattering channel only one of the resonance widths grows further with the coupling $`g`$. The resulting redistribution of S-matrix poles was coined resonance trapping.

Chaotic billiards with attached waveguides are considered as paradigm for chaotic scattering DS92a . They are relevant as theoretical models for understanding the transport properties of mesoscopic semiconductor structures M+92  or experimental results on microwave scattering in flat resonators SS90 ; P+00 . Also signatures of resonance trapping were recently observed in billiards both, numerically P+98 ; R+00 ; S+00  and in microwave resonator experiments P+00 .

However, to our knowledge there is no theory which maps a given billiard to an effective Hamiltonian with overcritical coupling, thus really establishing a connection between the numerically and experimentally observed phenomena and the results on resonance trapping obtained within the formalism (1)-(5).

Motivated by this situation, it is the purpose of the present paper to discuss the application of the Heidelberg approach to open billiards in some detail and to answer questions such as: How can (1), (3) be derived for a billiard, what kind of approximations are involved and what is the resulting expression for the coupling constants $`W_{n,\lambda }(E)`$? What is the influence of the choice of the internal subsystem which is not unique for a given scattering system?

Using the expression for the coupling constants $`W_{n,\lambda }`$ to be derived in section II we will then address the effective coupling constant for a typical chaotic billiard with an attached waveguide. We show in section III that in the semiclassical regime, and when no tunneling barriers obstruct the waveguides, the coupling strength is fixed at the critical value $`g=1`$, independently of the size or the precise geometry of the billiard and of its openings. A numerical verification of our results is contained in section IV, followed by a short discussion of the implications of our findings.

## II Coupling Constants for Individual Levels

We consider a situation as shown in Fig. 1. A scattering system is formed in two dimensions by an infinite waveguide of width $`b`$ and and an arbitrary cavity. Inside the system the potential is identically 0. We set $`\mathrm{}=2m=1`$ and $`E=k^2`$ such that the stationary Schrödinger equation reduces to the Helmholtz equation

$$(\mathrm{\Delta }+k^2)\mathrm{\Psi }(x,y)=0,$$
(6)
with $`\mathrm{\Delta }=^2/_x^2+^2/_y^2`$. On the boundary of the scattering system (bold solid line in Fig. 1) we require Dirichlet b. c. $`\mathrm{\Psi }=0`$. This boundary condition and also the precise geometry of the system are by no means essential, the following generalizes e. g. immediately to a cavity with more than one attached lead or Neumann b. c. $`/_n\mathrm{\Psi }=0`$. In Fig. 1 we have shown several possibilities to define a closed billiard which corresponds to the scattering system in question (solid and dotted thin lines). We will restrict the discussion to the case shown with a solid line: We require that the boundary of the internal system is located inside the attached waveguide and that it consists of a transversal straight line on which either Dirichlet or Neumann boundary conditions are imposed to close the system. Clearly, even this restriction makes the correspondence between the scattering system and the auxiliary internal system not unique because the exact position of the closure along the waveguide is variable. We use coordinates where this closure is at $`x=0`$ while the matching between waveguide and cavity is at $`x=L`$ $`(L>0)`$.

In the region of the attached waveguide ($`xL`$) we can decompose any function into transversal modes

$$\varphi _\lambda (y)=\sqrt{\frac{2}{b}}\mathrm{sin}\left(\lambda \frac{\pi y}{b}\right)(\lambda =1,2,\mathrm{}),$$
(7)
because these functions form a complete and orthonormal basis on the interval $`(0,b)`$ according to

$$\underset{\lambda =1}{\overset{\mathrm{}}{}}\varphi _\lambda (y)\varphi _\lambda (y^{})=\delta (yy^{})(0<y,y^{}<b)$$
(8)
and

$$_0^bdy\varphi _\lambda (y)\varphi _\lambda ^{}(y)=\delta _{\lambda \lambda ^{}}.$$
(9)
The most general solution of the Helmholtz equation is a superposition of the scattering states $`\mathrm{\Psi }_\lambda (x,y)`$ which consist of a single incoming wave in transversal mode $`\lambda `$ and the corresponding outgoing modes given by the S-matrix of the system

$$\mathrm{\Psi }_\lambda (x,y)=\varphi _\lambda (y)\frac{\mathrm{e}^{\mathrm{i}k_\lambda x}}{\sqrt{k_\lambda }}+\underset{\lambda ^{}}{}S_{\lambda ^{}\lambda }\varphi _\lambda ^{}(y)\frac{\mathrm{e}^{+\mathrm{i}k_\lambda ^{}x}}{\sqrt{k_\lambda ^{}}}.$$
(10)
The longitudinal wave number $`k_\lambda =\sqrt{k^2(\lambda \pi /b)^2}`$ is real for $`\lambda \mathrm{\Lambda }=[kb/\pi ]`$, where $`[\mathrm{}]`$ denotes the integer part. These modes are called open or travelling, and the $`\mathrm{\Lambda }\times \mathrm{\Lambda }`$ matrix $`S_{\lambda ^{}\lambda }`$ corresponding to the open modes is the unitary S-matrix we are interested in. For $`\lambda >\mathrm{\Lambda }`$ the momentum along the waveguide is imaginary. These modes are called closed or evanescent. In the scattering state $`\mathrm{\Psi }_\lambda `$ with $`\lambda \mathrm{\Lambda }`$ the evanescent outgoing modes describe exponentially decaying contributions which modify the wave function in the vicinity of the mouth of the waveguide. Evanescent incoming modes are exponentially increasing into the waveguide and thus unphysical for the scattering system.

When the evanescent modes are included, the S-matrix becomes an infinite-dimensional operator which is no longer unitary. It is possible to construct an eigenstate of the closed billiard by a superposition of the scattering states (10) including evanescent modes: Suppose $`S(E)`$ has an eigenvalue unity at some energy $`E_n=k_n^2`$ and let $`a_{n,\lambda }`$ be the components of the corresponding eigenvector. Then the linear combination of scattering states

$`\mathrm{\Psi }_n^{(\mathrm{N})}(x,y)`$ $`=`$ $`{\displaystyle \underset{\lambda }{}}{\displaystyle \frac{a_{n,\lambda }^{(\mathrm{N})}}{2}}\mathrm{\Psi }_\lambda (x,y)`$ (11)
$`=`$ $`{\displaystyle \underset{\lambda }{}}{\displaystyle \frac{a_{n,\lambda }^{(\mathrm{N})}}{\sqrt{k_{n,\lambda }}}}\varphi _\lambda (y)\mathrm{cos}(k_{n,\lambda }x)`$
$`(Lx0)`$ satisfies Neumann b. c. at $`x=0`$ and it is thus indeed an eigenfunction of the billiard because it satisfies (6) and the remaining boundary conditions by construction. The normalization of the S-matrix eigenvector $`a_{n,\lambda }`$ in (11) is not unity but rather determined by the normalization of the billiard eigenfunction $`\mathrm{\Psi }_n^{(\mathrm{N})}(x,y)`$. If Dirichlet boundary conditions are required at $`x=0`$ the same argument can be repeated for an eigenvalue $`1`$ of the S-matrix and we have

$`\mathrm{\Psi }_n^{(\mathrm{D})}(x,y)`$ $`=`$ $`{\displaystyle \underset{\lambda }{}}{\displaystyle \frac{a_{n,\lambda }^{(\mathrm{D})}}{2\mathrm{i}}}\mathrm{\Psi }_\lambda (x,y)`$ (12)
$`=`$ $`{\displaystyle \underset{\lambda }{}}{\displaystyle \frac{a_{n,\lambda }^{(\mathrm{D})}}{\sqrt{k_{n,\lambda }}}}\varphi _\lambda (y)\mathrm{sin}(k_{n,\lambda }x)`$
$`(Lx0)`$. Consequently the spectrum of the billiard closed with Neumann or Dirichlet b. c. can be found from the secular equation

$$det(IS(E))=0,$$
(13)
which was first derived by Doron and Smilansky DS92b . In a sense we will in the following invert this so-called scattering approach to the quantization of billiards: we will express the S-matrix in terms of the eigenvalues and eigenfunctions of the closed system.

For this purpose consider the Green function of the closed billiard which is defined as the resolvent of $`\mathrm{\Delta }`$ in the space of functions which satisfy the boundary conditions of the billiard. In position representation this definition can be expressed by the inhomogeneous Helmholtz equation

$$(\mathrm{\Delta }+k^2)G(𝐫;𝐫^{},k)=\delta (𝐫𝐫^{}).$$
(14)
In the eigenbasis of the billiard the Green function reads

$$G(𝐫;𝐫^{},k)=\underset{n=1}{\overset{\mathrm{}}{}}\frac{\mathrm{\Psi }_n^{}(𝐫^{})\mathrm{\Psi }_n(𝐫)}{k^2k_n^2},$$
(15)
which can be verified using $`\mathrm{\Delta }\mathrm{\Psi }_n(𝐫)=k_n^2\mathrm{\Psi }_n(𝐫)`$ and the completeness of the functions $`\mathrm{\Psi }_n`$ inside the billiard. For $`𝐫`$ and $`𝐫^{}`$ inside the waveguide we can expand the Green function with respect to the transversal modes $`\varphi _\lambda (y)`$ and find as the general form of a solution of (14)

$`G(𝐫,𝐫^{},k)`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{i}}}{\displaystyle \underset{\lambda \lambda ^{}}{}}{\displaystyle \frac{\varphi _\lambda (y)}{\sqrt{k_\lambda }}}{\displaystyle \frac{\varphi _\lambda ^{}(y^{})}{\sqrt{k_\lambda ^{}}}}`$
$`\times \left[\delta _{\lambda \lambda ^{}}\mathrm{e}^{\mathrm{i}k_\lambda \left|xx^{}\right|}+{\displaystyle \underset{s,s^{}=\pm }{}}G_{\lambda \lambda ^{}}^{ss^{}}(k)\mathrm{e}^{\mathrm{i}sk_\lambda x+\mathrm{i}s^{}k_\lambda ^{}x^{}}\right].`$
Indeed, the first term inside the brackets gives rise to one particular solution of the inhomogeneous equation (14), while the second term with the unknown matrices $`G^{++}`$, $`G^+`$, $`G^+`$ and $`G^{}`$ represents the most general solution of the corresponding homogeneous Helmholtz equation (6). The unknown coefficients must be determined such that the Green function satisfies also the boundary conditions inside the cavity and on the transversal closure of the waveguide. For this purpose assume first $`xx^{}`$ and consider $`𝐫`$ as a fixed parameter. Then $`G(𝐫;𝐫^{},k)`$ as a function of $`𝐫^{}`$ should satisfy the homogeneous Helmholtz equation with the boundary conditions of the scattering system, i. e. it can be written as a superposition of the scattering states $`\mathrm{\Psi }_\lambda (𝐫^{})`$ defined in Eq. (10). On the other hand, when $`𝐫^{}`$ is fixed, the Green function as a function of $`𝐫`$ satisfies the boundary conditions (Neumann or Dirichlet) at $`x=0`$ where the closed billiard is separated from the waveguide by the additional straight wall. Thus, it must be a superposition of the functions

$$\mathrm{\Psi }_\lambda ^{\mathrm{N}/\mathrm{D}}(𝐫)=\frac{\varphi _\lambda (y)}{\sqrt{k_\lambda }}(\mathrm{e}^{\mathrm{i}k_\lambda x}\pm \mathrm{e}^{\mathrm{i}k_\lambda x}),$$
(17)
which are in fact the scattering states for a semi-infinite waveguide with Neumann or Dirichlet b.c. at one end. Consequently, the Green function has the form

$$G(𝐫;𝐫^{},k)=\frac{1}{2\mathrm{i}}\underset{\lambda ,\lambda ^{}}{}\mathrm{\Psi }_\lambda ^{\mathrm{N}/\mathrm{D}}(𝐫)g_{\lambda \lambda ^{}}(k)\mathrm{\Psi }_\lambda ^{}(𝐫^{}),$$
(18)
with another set of undetermined coefficients $`g_{\lambda \lambda ^{}}`$. Expanding (18) into transversal modes and comparing to (II) we obtain

$$\begin{array}{cc}G^{++}=gS\hfill & G^+=gI\hfill \\ G^+=\pm gS\hfill & G^{}=\pm g\hfill \end{array}.$$
(19)
We can now repeat this argumentation under the opposite assumption $`x<x^{}`$ and find again the relations (19) but with $`G^+`$ and $`G^+`$ exchanged. This can be regarded as a consequence of the symmetry of the Green function with respect to its two arguments which in turn follows from time-reversal symmetry. We conclude $`G^+=G^+=gI=\pm gS`$ and hence

$$g(k)=(IS(k))^1.$$
(20)
Note that $`g(k)`$ and thus the Green function diverges as expected at the solutions of the secular equation (13), i. e., when $`k`$ corresponds to an eigenvalue of the closed billiard.

Using Eq. (20) we can now directly relate the S-matrix to the transversal expansion coefficients of the Green function at the closure of the billiard. We define for $`x,x^{}`$ inside the waveguide

$$G_{\lambda ,\lambda ^{}}(x,x^{})=_0^bdydy^{}\varphi _\lambda (y)G(x,y;x^{},y^{};k)\varphi _\lambda ^{}(y^{})$$
(21)
and

$$K_{\lambda \lambda ^{}}^{(\mathrm{N})}=\sqrt{k_\lambda k_\lambda ^{}}G_{\lambda ,\lambda ^{}}(0,0)$$
(22)
$$K_{\lambda \lambda ^{}}^{(\mathrm{D})}=\frac{1}{\sqrt{k_\lambda k_\lambda ^{}}}\frac{^2}{_x_x^{}}G_{\lambda ,\lambda ^{}}(x,x^{})|_{x=x^{}=0}$$
(23)
and derive from (18) using (20)

$$\mathrm{i}K^{(\mathrm{N})}=\frac{I+S}{IS}(\text{Neumann b. c. at }x=0)$$
(24)
$$\mathrm{i}K^{(\mathrm{D})}=\frac{IS}{I+S}(\text{Dirichlet b. c. at }x=0).$$
(25)
Obviously, $`K^{(\mathrm{N})}=0`$ for Dirichlet b. c. and $`K^{(\mathrm{D})}=0`$ for Neumann b. c. at $`x=0`$. Eqs. (24), (25) can be inverted and yield

$$S=\frac{I\mathrm{i}K^{(\mathrm{N})}}{I+\mathrm{i}K^{(\mathrm{N})}}=+\frac{I\mathrm{i}K^{(\mathrm{D})}}{I+\mathrm{i}K^{(\mathrm{D})}},$$
(26)
which is now in the form of (1) (for Neumann b. c. up to an irrelevant constant phase). We can now proceed to determine the corresponding coupling constants $`W_{n,\lambda }^{(\mathrm{N}/\mathrm{D})}`$ by representing the K-matrix in the eigenbasis of the billiard. From (2) we have

$$K_{\lambda \lambda ^{}}=\pi \underset{n=1}{\overset{\mathrm{}}{}}\frac{W_{n,\lambda }^{}W_{n,\lambda ^{}}}{k^2k_n^2}$$
(27)
and from (15), (22), (23) we find

$$K_{\lambda \lambda ^{}}^{(\mathrm{N})}=\sqrt{k_\lambda k_\lambda ^{}}\underset{n=1}{\overset{\mathrm{}}{}}\frac{\mathrm{\Psi }_{n,\lambda }^{(\mathrm{N})}(0)\mathrm{\Psi }_{n,\lambda ^{}}^{(\mathrm{N})}(0)}{k^2k_n^2}$$
(28)
$$K_{\lambda \lambda ^{}}^{(\mathrm{D})}=\frac{1}{\sqrt{k_\lambda k_\lambda ^{}}}\underset{n=1}{\overset{\mathrm{}}{}}\frac{\frac{}{x}\mathrm{\Psi }_{n,\lambda }^{(\mathrm{D})}(0)\frac{}{x}\mathrm{\Psi }_{n,\lambda ^{}}^{(\mathrm{D})}(0)}{k^2k_n^2}$$
(29)
where we have introduced the projections

$$\mathrm{\Psi }_{n,\lambda }^{(\mathrm{N}/\mathrm{D})}(x)=_0^bdy\varphi _\lambda (y)\mathrm{\Psi }_n^{(\mathrm{N}/\mathrm{D})}(x,y)$$
(30)
of the eigenfunctions of the closed billiard onto the transversal modes of the waveguide. The values of the coupling constants follow from comparing (27) to (28) and (29)

$$W_{n,\lambda }^{(\mathrm{N})}=\sqrt{\frac{k_\lambda }{\pi }}\mathrm{\Psi }_{n,\lambda }^{(\mathrm{N})}(0),$$
(31)
$$W_{n,\lambda }^{(\mathrm{D})}=\frac{1}{\sqrt{k_\lambda \pi }}\frac{}{x}\mathrm{\Psi }_{n,\lambda }^{(\mathrm{D})}(0).$$
(32)
This form of the dependence of the coupling constants on the internal wave functions is not surprising: also within perturbation theory the coupling depends on the value of the wave function at the point where the system is opened or on its normal derivative for Neumann and Dirichlet boundary conditions, respectively. However, in the situation we consider perturbation theory is not applicable and, in particular, the precise value of the prefactor in the coupling constants (31), (32) could only be obtained from the derivation given in this section.

The representation (26) of the S-matrix in terms of the K-matrices (28), (29) is exact, when all transversal modes are included. However, usually one is interested only in the $`\mathrm{\Lambda }\times \mathrm{\Lambda }`$ unitary part of the S-matrix, and this is only approximately given by (26) when the K-matrix is restricted to open modes. In DS92b ; SS95  the effect of this so-called semiquantal approximation for the accuracy of eigenvalues within the scattering approach to quantization was investigated numerically. It becomes negligible when the energy $`E`$ is sufficiently far from the threshold for the opening of a new channel. Under this restriction we can consider $`W_{n,\lambda }^{(\mathrm{N}/\mathrm{D})}`$ as the coupling constants corresponding to the unitary part of the S-matrix. When $`E`$ approaches a threshold, Eq. (3) breaks down, since the energy-dependence of the coupling constants can no longer be neglected. We will not consider this case here.

## III The Mean Coupling Strength

Given the explicit values (31), (32) for the coupling constants between individual states and individual scattering channels we are now going to derive an estimate for the dimensionless coupling strength (5) in the semiclassical limit and neglecting evanescent modes. Since the concept of a mean coupling strength is not well-defined for infinitely many internal states, the internal Hamiltonian $`H_{\mathrm{in}}`$ entering (2), (4) should for this purpose be cut to some finite matrix including only states which are close enough in energy $`k_nk`$. In particular this means that we can replace the momenta along the waveguide $`k_{n,\lambda }`$ in (11), (12) by their on-shell values $`k_\lambda `$. The resulting approximate expansions of the billiard eigenfunctions are projected onto the transversal modes according to (30) and inserted into (31), (32) which simplify to

$$W_{n,\lambda }^{(\mathrm{N}/\mathrm{D})}=\frac{a_{n,\lambda }^{(\mathrm{N}/\mathrm{D})}}{\sqrt{\pi }}.$$
(33)
At this point it is necessary to determine the average magnitude of the coefficients $`a_{n,\lambda }^{(\mathrm{N}/\mathrm{D})}`$. We assume that the classical dynamics of the billiard is chaotic. In the semiclassical limit this means that the quantum ergodicity theorem applies to the eigenstates of the billiard, i. e. in particular the probability density integrated over an arbitrary region of the billiard tends to the relative area of that region. Applied to the part of the billiard inside the waveguide we find

$`{\displaystyle \frac{bL}{A}}`$ $``$ $`{\displaystyle _0^b}dy{\displaystyle _L^0}dx|\mathrm{\Psi }_n^{(\mathrm{N})}(x,y)|^2`$ (34)
$`=`$ $`{\displaystyle _L^0}dx{\displaystyle \underset{\lambda }{}}{\displaystyle \frac{|a_{n,\lambda }|^2}{k_\lambda }}\mathrm{cos}^2(k_\lambda x)`$
$``$ $`|a|^2{\displaystyle \frac{L}{2}}{\displaystyle _0^{kb/\pi }}d\lambda \sqrt{{\displaystyle \frac{1}{k^2(\lambda \pi /b)^2}}}`$
$`=`$ $`|a|^2{\displaystyle \frac{bL}{4}},`$
where $`A`$ denotes the total area of billiard. In the second line we have inserted the normal mode decomposition into $`|\mathrm{\Psi }_n^{(\mathrm{N})}(x,y)|^2`$. The orthonormalization of the transversal modes was then used to restrict the resulting double sum over modes to diagonal terms. In the third line $`\mathrm{cos}^2=1/2`$ was used ($`\mathrm{sin}^2=1/2`$ in the completely analogous calculation for Dirichlet b.c.), and the sum over modes was approximated by a continuous integral. This is justified when the number of modes is large, i. e. in the semiclassical limit. Using the resulting constraint on the normalization of the coefficients of the billiard eigenfunctions in the transversal basis implied by (34), $`|a_{n,\lambda }^{(\mathrm{N}/\mathrm{D})}|^2=4/A`$, we find

$$|W_{n,\lambda }^{(\mathrm{N}/\mathrm{D})}|^2=\frac{4}{A\pi }$$
(35)
According to (5) the average coupling between the internal states and the continuum must be normalized by the mean level spacing $`D`$ of the billiard which is the only independent energy scale of the system. To leading semiclassical order we have Weyl’s law $`D=4\pi /A`$ BH76  which finally results in

$$g=1.$$
(36)
## IV Numerical Results

To check the validity of the Heidelberg approach for a quantum billiard we have performed direct numerical calculations for a Sinai billiard connected to a single waveguide (see Fig. 2). First, we evaluated numerically the $`1300`$ lowest eigenvalues and eigenvectors of the closed system with Dirichlet and Neumann b.c. at the boundary segment to which the waveguide was attached. Using the expressions (30)-(32) we calculated the elements of the coupling matrix $`W_{n,\lambda }^{(\mathrm{N}/\mathrm{D})}`$. The first observation is that the distribution of the numerically obtained values of $`W_{n,\lambda }^{(\mathrm{N}/\mathrm{D})}`$ is fairly close to a Poisson distribution, that is expected for a fully random coupling. Hence the elements of the coupling matrix $`W`$ evaluated with help of the formulae (30)-(32) resemble the random coupling standardly used in the Heidelberg approach.

Knowing $`W`$ and using Eqs. (1), (2) we evaluated—as the next step—the $`S`$-matrix and compared it to the $`S`$-matrix obtained by a direct method based on the numerical solution of the underlying Schrödinger equation (see Ref. A91  for details). In order to visualize the differences between these two $`S`$-matrices we compared first the corresponding Wigner-Smith time delays

$$\tau (E)=\frac{\mathrm{i}}{M}\mathrm{Tr}\left(\frac{S^{}(E)}{E}S(E)\right),$$
(37)
where $`M`$ is the number of the open channels inside the waveguide. The results are plotted in Fig. 3.

The figure shows an outstanding agreement of the two approaches for low energies. The Heidelberg approach describes the positions as well as the widths of the narrow resonances with high accuracy. The agreement is good even close to the threshold energies of the individual channels. For higher energies, however, the difference between the two time delay functions increases due to the limited number of internal states included into the evaluation of the Heidelberg $`S`$-matrix. Similar results (not displayed) were obtained also with the waveguide attached to the boundaries No. $`2`$, $`3`$ and $`4`$ respectively. In all these cases the Neumann matching procedure was used.

On the other hand, for Dirichlet matching (Eq. 32) the results change drastically. In this case the agreement is not good even for low energies. This may seem surprising because the derivations of the previous sections were entirely parallel for Neumann and Dirichlet b. c. However, an important difference is hidden in the convergence properties of the spectral decompositions of the K-matrix (28), (29) as we shall explain now. Projecting (11), (12) onto transversal mode $`\lambda `$ we find

$$\mathrm{\Psi }_{n,\lambda }^{(\mathrm{N})}(0)=\frac{a_{n,\lambda }^{(\mathrm{N})}}{\sqrt{k_{n,\lambda }}}\frac{}{x}\mathrm{\Psi }_{n,\lambda }^{(\mathrm{D})}(0)=\sqrt{k_{n,\lambda }}a_{n,\lambda }^{(\mathrm{D})}.$$
(38)
Since the coefficients $`a_{n,\lambda }^{(\mathrm{N}/\mathrm{D})}`$ are according to (34) of order $`2/\sqrt{A}`$ which is independent of $`n`$, we have from (31), (32)

$$W_{n,\lambda }^{(\mathrm{N})}\sqrt{\frac{k_\lambda }{k_{n,\lambda }}}W_{n,\lambda }^{(\mathrm{D})}\sqrt{\frac{k_{n,\lambda }}{k_\lambda }}.$$
(39)
For a 2D billiard $`k_n,k_{n,\lambda }=O(\sqrt{n})`$ $`(n\mathrm{})`$ such that the terms in the infinite spectral sum (27) decay asymptotically as $`n^{3/2}`$ for Neumann and as $`n^{1/2}`$ for Dirichlet b.c. Hence the convergence is absolute for Neumann b.c. while (29) converges only conditionally. As a consequence, the numerically necessary cut-off in the summation over the internal states $`n`$ introduces large errors for Dirichlet b.c.

For the following considerations we will concentrate on Neumann b.c. As already mentioned, the choice of the internal and external part of the system is not unique, since an arbitrary part of the ideal waveguide can be considered part of the internal system. Increasing the length of the waveguide included, we decrease in fact the influence of the evanescent modes since they are exponentially vanishing inside waveguide. We have checked this relation and evaluated the time delay functions also for various waveguide parts included into the internal system. The results remain practically unchanged regardless on the length of the included part. This demonstrates the small influence of the evanescent modes on the resulting Heidelberg $`S`$-matrix.

Knowing the coupling matrix $`W`$ and using the relation Eq. (5) we evaluated numerically the value of the coupling constant $`g`$. The obtained result is in an excellent agreement with the estimated value (36) for all considered types of the waveguide attachment leading to $`g0.98`$. We have evaluated the coupling constant $`g`$ also for a different shape of the billiard P+98  obtaining similar values for $`g`$. It has to be stressed that the estimate (36) was obtained using semiclassical arguments. Our calculation shows, however, that it remains valid even in the deep quantum region.

The eigenvalues of $`H_{\mathrm{eff}}`$ (Eq. 4) are usually interpreted as the resonance poles and are used for the study of the statistical properties of resonances in open quantum chaotic systems FS97 . To check the validity of this approach we have evaluated the resonance poles of the system independently using the complex scaling method S+00  that provides a direct access to the positions of the poles of the analytically continued $`S`$-matrix. The obtained results were compared with the eigenvalues of the effective Hamiltonian (4). However, a direct comparison is obscured by the fact that the coupling matrix $`W`$ is in fact energy dependent. In the standard random matrix approach the coupling matrix is treated as being energy independent—a simplification that is well justified inside a small energy interval. To mimic this situation and to minimize the influence of the energy dependence of the coupling matrix $`W`$ we have compared the eigenvalues of $`H_{\mathrm{eff}}`$ with the directly evaluated resonance poles always within a small energy interval the center of which was equal to the energy used to evaluate the coupling matrix $`W`$.

The energy dependence of $`W`$ can be taken into account more precisely using the relation (1) and evaluating the Heidelberg resonance poles as zeros of the function $`I+\mathrm{i}K^{(\mathrm{N})}(E)`$. We have evaluated the Heidelberg resonance poles using both of the above described approaches and compared the results to the resonances obtained by complex scaling. The results are shown in Fig. 4. From this figure we see that for narrow resonances the eigenvalues of effective Hamiltonian $`H_{\mathrm{eff}}`$ represent a good approximation to the resonance poles of the system and the energy dependence of the coupling matrix $`W`$ can be omitted. For broad resonances the situation changes and the complex eigenvalues of the effective Hamiltonian have nothing in common with the directly evaluated resonance poles. This discrepancy can be explained as follows: A resonance localized at $`E_R=E+i\mathrm{\Gamma }`$ represents in fact a collective mode of all bound states $`E_n`$ of the internal Hamiltonian $`H_{in}`$ that are located inside the energy interval $`(E\mathrm{\Gamma },E+\mathrm{\Gamma })`$. For a broad resonance with $`\mathrm{\Gamma }`$ significantly larger then the mean spacing between the bound states $`E_n`$ the number of the internal states to be included into $`H_{\mathrm{eff}}`$ must be very high.

## V Conclusions

To summarize, we have shown that the Heidelberg approach to scattering, which is the basis of many important random-matrix results on quantum chaotic scattering, leads to reasonably good agreement when compared with the results of a direct calculation of the S-matrix in the case of billiards with Neumann boundary conditions. We have explained the somewhat unexpected finding that the accuracy is much worse for Dirichlet boundary conditions, while it does not depend very much on other possibilities of varying the auxiliary closed system used, such as the position of the attached waveguide.

Even for Neumann boundary conditions the effective Hamiltonian $`H_{\mathrm{eff}}`$ based on the Heidelberg approach seems to be not very well suited for the computation of broad resonances of the system. This result does not contradict the fact, that also for systems with time-reversal symmetry the statistical properties of billiard resonances follow the predictions of random-matrix theory based upon the effective Hamiltonian approach quite well Ish00 , because our test goes way beyond a purely statistical analysis.

Moreover, we have shown numerically and with semiclassical arguments that the mean dimensionless coupling for a chaotic billiard is the critical value $`g=1`$—irrespective of the precise form of the billiard, the size of the attached waveguide and other details of the model. Interestingly, in our model an effective coupling near $`g=1`$ is observed already deep in the quantum regime. Nevertheless, fluctuations around the mean value $`g=1`$ should in general be largest for small energies and can possibly result locally in overcritical coupling. This might be an explanation for the observed resonance trapping in billiards P+98 ; R+00 ; S+00 ; P+00 .

Our results concerning the value of the effective coupling for chaotic systems are not restricted to billiards and apply, e. g., to quantum graphs as well. In these systems, a systematic way to achieve overcritical coupling for many states is to modulate the density of states, e. g. by considering systems with band spectra Dit00 ; SD01 .

###### Acknowledgements.

This work has been partially supported by the Czech grand GAAV A1048804 and by the ”Foundation for Theoretical Physics” in Slemeno, Czech Republic.


format: text

text:
# Ab initio calculations for bromine adlayers on the Ag(100) and Au(100) surfaces: the 𝑐(2×2) structure

## I Introduction

Anion adsorption on metals can strongly modify surface morphology and electronic structure and chemical reactivity. It is therefore of great scientific and technological importance. In particular, the halide-adsorbed noble-metal systems play a significant role in electrochemistry. From the fundamental point of view, halide-adsorbed noble-metal surfaces are important model systems for adsorption on metal surfaces with formation of ordered two-dimensional adsorbate structures. It is thus not surprising that the adsorption of halides on noble metals has been extensively investigated.

The systems we selected to study are the Br-chemisorbed Ag(100) and Au(100) surfaces. The adsorption of bromine on the Ag(100) and Au(100) surfaces both in vacuum and in solution have been widely studied by experiments  and by classical simulations.  Experimentally, bromine has been found to form different bonding structures on the Ag(100) and Au(100) surfaces. While bromine chemisorbed on the Ag(100) surface occupies the four-fold hollow site (hereafter referred to as $`H_4`$), the most stable chemisorption structure on Au(100) is the configuration with bromine at the two-fold bridge site (hereafter referred to as $`B_2`$). These different chemisorption structures have been verified by various experimental measurements. However, theoretical studies have not yet reproduced these different adsorption behaviors.

Kleinherbers et al. performed angle-resolved photoemission, low-energy electron diffraction (LEED), and X-ray photoemission measurements for the interaction of halides with Ag surfaces.  They found that the adsorption of Cl, Br, and I on the Ag(100) surface in vacuum all resulted in the formation of a $`c(2\times 2)`$ overlayer with the adsorbates in the $`H_4`$ sites. Using in situ surface X-ray scattering, Ocko et al. studied the adsorption of bromide on an Ag(100) electrode. They observed a disordered phase at lower coverages and an ordered $`c(2\times 2)`$ phase at a coverage of half a monolayer.  The Br was determined to bond at the $`H_4`$ site in both the $`c(2\times 2)`$ and disordered phases. The atomic geometry of the Br/Ag(100)-$`c(2\times 2)`$ surface was further investigated by Endo et al. with the in situ X-ray absorption fine structure (XAFS) method.  The $`H_4`$ site was confirmed to be the bonding site of bromine. The disordered phase of the Br/Ag(100) surface at lower bromine coverages was also recently further investigated experimentally.  In that study, it was suggested that while most of the bromide ions occupy the $`H_4`$ sites, there are additional bromide ions adsorbed slightly off the $`H_4`$ sites.

The LEED data reported by Bertel et al.  have shown that the chemisorption of Br on the Au(100) surface in vacuum results in the rearrangement of the top-layer Au atoms of the original clean reconstructed Au(100)-$`(5\times 20)`$ surface and the formation of an unreconstructed $`(1\times 1)`$ substrate structure. Several ordered structures of the Br adlayer, including $`c(2\times 2)`$, $`(\sqrt{2}\times 4\sqrt{2})R45^{}`$, and $`c(4\times 2)R45^{}`$, were obtained after bromine exposure on Au(100) surfaces, with the former two structures being metastable. It was concluded from the experimental data that Br adsorbed at the $`B_2`$ site on the Au(100) surface in all the observed phases. This is in contrast to the case of the Br-adsorbed Ag(100) surface. Under electrochemical in situ conditions, surface X-ray scattering and scanning tunneling microscopy (STM) experiments showed that bromide adsorbed on the unreconstructed Au(100)-$`(1\times 1)`$ surface forms a commensurate $`c(\sqrt{2}\times 2\sqrt{2})R45^{}`$ structure and an incommensurate $`c(\sqrt{2}\times 2p)R45^{}`$ ($`p2\sqrt{2}`$, depending on the applied potential).  In this case, too, the bromide ions were determined to reside at the $`B_2`$ sites.

In contrast to the considerable progress of the experimental measurements, theoretical studies employing ab initio methods to these systems are still at an early stage. Several groups have performed ab initio Hartree-Fock (HF) and density-functional-theory (DFT) calculations for the Br/Ag(100) and Br/Au(100) interfaces using cluster models.  While these investigations have provided useful information about the interaction between Br and the surfaces, as we discuss below, many of the results are not yet sufficiently accurate. For example, the preferred bonding site of Br on the Au(100) surface was incorrectly predicted by ab initio DFT cluster calculations, which showed that Br would prefer to bond at the $`H_4`$ site on both the Ag(100) and Au(100) surfaces. Similarly, ab initio HF studies with small clusters predicted the bridge site as the preferred adsorption site for Br on Ag(100).  Given that these calculations with small clusters cannot reproduce such a fundamental property as the binding site, all other results (e.g., energy barriers) obtained from such calculations for both the Br/Ag(100) and Br/Au(100) systems are questionable.

Here we present results of total-energy DFT calculations in which we used supercell models for the Ag(100) and Au(100) surfaces. The detailed atomic structures and electronic properties of the chemisorbed surfaces and the preferred bonding site of the adsorbate have been determined. Our theoretical approach has reproduced the different behavior of Br on the Ag(100) and Au(100) surfaces. The most stable adsorption sites for Br chemisorbed on the Ag(100) and Au(100) surfaces are determined by our calculations to be the $`H_4`$ and the $`B_2`$ sites, respectively. Our results are in excellent agreement with the experimental data. The obtained results for the electronic properties also enable us to analyze the nature of the bonding between Br and the substrates and understand the different adsorption behavior of Br on the Ag(100) and Au(100) surfaces.

The remainder of this paper is organized as follows. In Sec. II we outline in detail the computational method and the supercell models that we used. In Sec. III we present and discuss the results for bulk, clean surfaces (Sec. III A), and adsorbed surfaces (Sec. III B). The adsorption geometries and atomic relaxations are discussed in Sec. III B1, and the electronic properties and bonding character in Sec. III B2. We also give comparisons of our results with previous calculations and experimental data. Finally, in Sec. IV, we summarize the main results of our calculations.

## II Method and model

On the Ag(100) and Au(100) surfaces, there are three different symmetric adsorption sites, known as $`H_4`$, $`B_2`$, and $`T_1`$ (on-top) sites. These three sites are shown in Fig. 1. We have studied a $`c(2\times 2)`$ structure in which Br

forms an adlayer with a coverage of $`\frac{1}{2}`$ monolayer on the surface for each of the three bonding configurations.

The metal surface is modeled by repeated slabs with five, seven, and nine metal layers separated by a vacuum region equivalent to five or seven metal layers. Each metal layer in the supercell contains two metal atoms. Br is adsorbed symmetrically on both sides of the slab. All the metal atoms were initially located at their bulk positions, with the equilibrium lattice constant of the bulk determined by our calculations.

The calculations were performed within density-functional theory, using the pseudopotential (PP) method and a plane-wave basis set. The results reported in this paper were obtained using the Vienna ab-initio simulation package (VASP).  The exchange-correlation effects were treated with the generalized gradient-corrected exchange-correlation functionals (GGA) given by Perdew and Wang.  We adopted the scalar-relativistic Vanderbilt ultrasoft pseudopotentials supplied by Kresse and Hafner . A plane-wave energy cutoff of 20 Ry and 56 special k points in the irreducible part of the two-dimensional Brillouin zone of the $`c(2\times 2)`$ surface were used for calculating both the Br/Ag(100) and Br/Au(100) surfaces. Optimization of the atomic structure was performed for each supercell via a conjugate-gradient technique using the total energy and the Hellmann-Feynman forces on the atoms.  All the structures were fully relaxed until the change in total energy was smaller than 1 meV between two ionic steps. The convergence of the total energies was checked with different values of the plane-wave cutoff and different numbers of special k points. A series of test calculations with different slab thicknesses (from five to nine metal layers) and vacuum-gap widths (equivalent to five and seven metal layers) were also carried out to check convergence. The calculations on which this paper is based represent approximately 300 CPU hours on an IBM SP2 computer.

## III Results and discussion

### A Bulk and clean surface

We first present the calculated properties for bulk silver and bulk gold, and the relaxed but unreconstructed clean Ag(100) and Au(100) surfaces.

Calculations for bulk Ag and Au were conducted with 408 special k points and cutoff energies ranging from 20 Ry to 40 Ry. The total energy convergence with respect to the cutoff energy was shown to be within a few tenths of 1 meV. We obtained lattice constants of 4.17 Å and 4.18 Å for bulk Ag and Au, about 2.0% and 2.5% larger than the corresponding experimental values  at room temperature, respectively. Previous total-energy DFT calculations at the GGA level found lattice constants between 4.13 Å and 4.19 Å for bulk Ag,  and between 4.19 Å \[Ref. \] and 4.20 Å \[Ref. \] for bulk Au. Our results are in good agreement with these calculations.

The properties of the clean Ag(100) and Au(100) surfaces were calculated using supercells containing a 7-layer metal slab and a vacuum gap equivalent to 7 bulk-metal layers. A $`1\times 1`$ surface cell was used with 66 special k points in the surface Brillouin zone. The kinetic energy cutoff for the calculations was 20 Ry. All the layers except for the central one were relaxed. Surface reconstruction was not considered.

The surface energies for the Ag(100) and Au(100) surfaces obtained from our calculations are 0.43 eV/atom and 0.47 eV/atom, respectively. Both results are in good agreement with recent pseudopotential GGA calculations by Yu and Scheffler, who reported the corresponding values of 0.48 and 0.45 eV/atom. A recent calculation with linear-muffin-tin-orbital (LMTO)-GGA methods, however, obtained much larger values of 0.65 and 0.90 eV/atom for the unrelaxed Ag(100) and Au(100) surfaces, respectively.  The reason for the large discrepancy from the other GGA results is not clear. Pseudopotential local-density approximation (LDA),  LMTO-LDA,  and linearized augmented-plane-wave LDA  calculations have provided values of 0.59–0.7 eV/atom, and 0.69–0.72 eV/atom for the surface energies of the Ag(100) and Au(100) surfaces, respectively. These DFT-LDA values for the surface energy are generally larger than those calculated with the DFT-GGA calculations reported in this paper and Refs.  and . This is consistent with the previous observation  that LDA surface energies are normally larger than the corresponding GGA values due to the different treatment of the exchange-correlation functional. The calculated values for the surface energies are thus seen to be quite sensitive to the computational method and the form of the exchange-correlation functional.

Table I shows the results of the surface relaxation. While no significant structural relaxation is found for either surface, the Ag(100) surface shows a slightly more relaxed geometry than the Au(100) surface. Both surfaces show an inward relaxation of the top layer and slight outward relaxation of the second and third layers. LEED measurements  showed insignificant relaxation of the Ag(100) surface with $`\mathrm{\Delta }d_{12}/d_0=0\pm 1.5\%`$ and $`\mathrm{\Delta }d_{23}/d_0=0\pm 1.5\%`$, where $`\mathrm{\Delta }d_{12}`$ and $`\mathrm{\Delta }d_{23}`$ are the changes in spacing between the top and the second layer and between the second and the third layer, and $`d_0`$ is the bulk interlayer distance. Our results are thus in good agreement with the experimental data, and basically consistent with other ab initio calculations, the results of which are also listed in Table I for comparison.

### B Br-adsorbed Ag(100) and Au(100) surfaces

#### 1 Relaxations and energetics

The results of our calculations for the Br/Ag(100)-$`c(2\times 2)`$ and the Br/Au(100)-$`c(2\times 2)`$ surfaces are shown in Tables II and III. If not otherwise indicated, the results reported in this section (and in Tables II and III) were calculated with a cutoff energy of 20 Ry, 56 special k points, and supercells containing a 9-layer metal slab with a vacuum region equivalent to seven metal layers.

In-plane relaxations of the top-layer metal atoms were found to result in changes of the distance between Br and its nearest-neighbor metal atoms and of the total-energy difference between two different configurations within only 0.01 Å and a few meV, respectively. The effects of in-plane relaxations are thus negligible, and such relaxations were not considered in the calculations for the Br-adsorbed surfaces.

Table II shows total-energy differences between different bonding configurations of both surfaces. Each structure was optimized. We found that while the total energy of the $`H_4`$ configuration is lower by 213 meV than the $`B_2`$ configuration for the Br/Ag(100) surface, it is higher by 58 meV for the Br/Au(100) surface. The $`T_1`$ configuration for both surfaces is found to be higher in total energy than both the corresponding $`H_4`$ and $`B_2`$ configurations. Thus, we conclude that while Br adsorbed on the Ag(100) surface prefers the $`H_4`$ site, it is adsorbed at the $`B_2`$ site on the Au(100) surface. This conclusion is in agreement with experimental observations.

It is interesting to note that the magnitude of the total-energy difference between the $`H_4`$ and $`B_2`$ structures for Br/Ag(100) is significantly larger than the corresponding value for the Br/Au(100) surface. This suggests that diffusion of Br on the Au(100) surface may occur much more easily than on the Ag(100) surface since the total-energy difference between the most stable configuration (the global minimum) and the less favorable configuration (probably a saddle point) is directly relevant to adsorbate diffusion.

Previous theoretical studies employing cluster models also determined the preferred bonding sites of Br on the Ag(100) and Au(100) surfaces, as mentioned in Sec. I. Ab initio HF calculations showed that Br would prefer to bond at the $`B_2`$ site on the Ag(100) surface (by 370 meV/adatom over the $`H_4`$ site, and by 570 meV/adatom over the $`T_1`$ site).  This is inconsistent with both our DFT-supercell calculations and the experimental data.  DFT cluster calculations predicted that the binding energy of Br at the $`H_4`$ site on both the Ag(100) and Au(100) surfaces was larger than at the $`B_2`$ and $`T_1`$ sites by 120 meV for Ag(100) \[89 meV for Au(100)\] and

202 meV for Ag(100) \[202 prefers to bond at the $`H_4`$ site on the Au(100) surface is in disagreement with experimental measurements,  as well as with our results. We believe that the main problem is that these previous calculations were limited to small clusters, containing only up to 13 metal atoms. It is well known that a small metal cluster has a very different electronic structure than an extended metal surface, yielding very significant differences in adsorbate binding energies and reaction pathways. Large clusters or extended surface models (e.g., supercell models) are therefore needed to simulate real metal surfaces accurately.

The structural parameters of the optimized geometries for the $`H_4`$, $`B_2`$, and $`T_1`$ configurations of the Br/Ag(100) and Br/Au(100) surfaces are presented in Table III. The vertical distances ($`d_z`$) between the Br centers and the plane of the centers of the top-layer atoms were calculated to be 1.91 Å on the Ag(100) surface and 2.01 Å on the Au(100) surface for the $`H_4`$ structure, 2.16 Å on Ag(100) and 2.18 Å on Au(100) for the $`B_2`$ configuration, and 2.48 Å on Ag(100) and 2.46 Å on Au(100) for the $`T_1`$ structure. While the values of $`d_z`$ for the $`B_2`$ and $`T_1`$ configurations of the Br/Ag(100) surface are very close to the corresponding values for the Br/Au(100) surface, the distance between Br and the surface in the $`H_4`$ configuration is observed to be significantly longer (by 0.1 Å) on Br/Au(100) than on Br/Ag(100). Accurate in-situ XAFS measurements for Br/Ag(100) in NaBr solution by Endo et al. showed that the bond length between Br and its four nearest-neighbor Ag atoms in the $`H_4`$ configuration is $`2.82\pm 0.05`$ Å, and the distance between the Br and the surface is $`1.94\pm 0.07`$ Å.  Our results (2.82 Å and 1.91 Å, respectively) are thus in excellent agreement with the experimental data, provided that the solution has only a minor influence on the bond lengths between the adsorbate and the surface.

The bond lengths between Br and Ag and Au clusters of varying size have been obtained with both HF and DFT calculations. Illas et al., using the HF method with a cluster of 5 Ag atoms simulating the $`H_4`$ configuration of the Ag(100) surface, obtained a value of 3.43 Å for the length of the Br-Ag bond. Paccioni, also using the HF method with slightly larger clusters, found that the bond lengths between a Br ion and the surface Ag atom were 3.24 Å, 2.97 Å, and 2.94 Å in clusters of Br--Ag13 (modeling the $`H_4`$ structure), Br--Ag8 (simulating the $`B_2`$ geometry), and Br--Ag13 (representing the $`T_1`$ configuration), respectively. Ignaczak and Gomes performed DFT calculations with clusters containing a Br ion and 12 metal atoms and determined the bond lengths to be 3.2 Å, 3.0 Å, and 2.9 Å for the $`H_4`$, $`B_2`$, and $`T_1`$ configurations of Br--Ag12 and Br--Au12 clusters, respectively. All of these values are much larger that those obtained from our supercell calculations and the XAFS measurements, suggesting that small clusters do not represent the metal surfaces properly.

The $`B_2`$ configuration of the Br/Ag(100) surface shows a very similar relaxation of the surface metal layers as that of the same configuration for the Br/Au(100) surface. Both undergo an inward relaxation of the top layer and slight outward relaxations of the second and third layers. Similar relaxed structures are also found for the clean Ag(100) and Au(100) surfaces (see Table I).

Our calculations show that the second metal layer undergoes a small buckling with the adsorption of Br in the $`H_4`$ configurations. The atoms in the second layer that are immediately below the $`H_4`$ sites are observed to shift slightly up towards the surface, while the other atoms in that layer shift up by only on the order of 0.001 Å and hence essentially keep their bulk positions. The spacing between these two sub-layers is found to be 0.02 Å and 0.04 Å for the Br/Ag(100) and Br/Au(100) surfaces. The distance between Br and the second-layer metal atom just below it is still far larger than the bond length between Br and its nearest-neighbor metal atoms in the top layer. Thus a pseudo-five-fold coordination, which has been observed in $`c(2\times 2)`$ overlayer structures on bcc metal surfaces,  does not exist for the Br/Ag(100) and Br/Au(100) surfaces. This buckling may give rise to an effective Br-Br interaction, mediated through the surface strain field. The top metal layer in the $`H_4`$ configurations still shows a slight inward relaxation, similar to the cases of the clean surfaces and the $`B_2`$ configurations.

The top metal layer of the $`T_1`$ configuration of the Br/Au(100) surface also shows a small buckling. The Au atoms in the top layer that are bonded to Br are observed to undergo a larger inward relaxation than the other half of the Au atoms in the top metal layer. The corresponding buckling is, however, very large for the $`T_1`$ configuration of the Br/Ag(100) surface. The distance between the two sublayers formed from the top Ag layer is 0.75 Å, indicating a zigzag surface reconstruction.

Finally, in Table IV we show results of convergence checks for the total-energy differences. Such checks are particularly important for the Br/Au(100) surface due to the small value of the total-energy difference between the $`H_4`$ and $`B_2`$ configurations. Calculations with a higher cutoff energy (30 Ry) obtained total-energy differences within 1 meV of those from calculations with a cutoff energy of 20 Ry. Increasing the number of special k points from 36 to 56, increasing the slab thickness from 7 to 9 metal layers, and increasing the vacuum region in the supercell from 5 to 7 layers, all changed the results by only a few meV. Supercells with 5 metal layers are seen to cause errors in the total-energy differences of $``$30 meV for Br/Ag(100) and $``$10 meV for Br/Au(100). The use of 20 special k points also causes an error of $``$10 meV. Therefore, it is necessary to employ supercells with at least 7 metal layers and 36 special k points for obtaining the total-energy differences with errors smaller than 10 meV. The distances between Br and its nearest-neighbor metal atoms were also checked. We found that the changes of these distances were smaller than 0.01 Å over the ranges of cut-off energies between 20 and 30 Ry, numbers of k points between 20 and 56, and numbers of metal layers between 5 and 9 in the supercells, indicating that the bond lengths are not very sensitive to the choice of computational parameters.

#### 2 Electronic properties and bonding character

In order to better understand the differences between the bonding of bromine on the Ag(100) and Au(100) surfaces, we calculated the total electronic density of states (DOS), the DOS projected onto individual atoms and specific atomic states, and the charge transfer between bromine and the substrate.

Figure 2 shows the total DOS for the Br-adsorbed Ag(100) and Au(100) surfaces. For comparison, the total DOS for the clean relaxed Ag(100) and Au(100) surfaces are also shown. The peaks of the DOS curves for the clean surfaces represent the main features of the $`s`$ and $`d`$ states of the substrates and remain essentially at the same positions when bromine atoms are adsorbed. New states are, however, found to be located at between $``$15 eV and $``$13 eV relative to the Fermi level in the DOS curves of the Br-adsorbed Ag(100) and Au(100) surfaces. These states are predominantly the bromine $`3s`$ states with small contributions from the $`s`$ and $`d`$ states of the substrate, as seen from the curves for the DOS projected onto the specific atomic states of the adsorbate and the substrate (shown in Fig. 3). Significant changes of the total DOS in the higher energy range close to the Fermi level (above $``$2.5 eV and $``$1.5 eV for the Br/Ag(100) and the Br/Au(100) surfaces, respectively) are also observed

when the DOS curves for the clean surfaces are compared with those for the $`H_4`$, $`B_2`$, and $`T_1`$ configurations of the Br-adsorbed surfaces (see Fig. 2). The electronic states in the higher energy range are composed mostly of the bromine $`3p`$ states and the Ag 4$`d`$ (or Au 5$`d`$) states, with some contributions coming from the Ag 5$`s`$ (or Au 6$`s`$) states (see Fig. 3).

In Fig. 4, we show the DOS projected onto bromine for the systems before and after adsorption. The results for the systems before adsorption were calculated by employing the supercell of a 7-layer slab and a 7-layer vacuum region. A bromine layer with a c($`2\times 2`$) periodicity was kept fixed in the middle of the vacuum region (located at $``$ 7 Å above the surface) so that there was essentially no interaction between Br and the substrate. The peaks located in the lower and higher energy ranges in the projected DOS before adsorption are due to the Br $`3s`$ and $`3p`$ states, respectively. The slight broadening of the $`3p`$ states reflects the weak $`3p`$-$`3p`$ interaction between neighboring bromine atoms. When bromine is adsorbed, both the $`3s`$ and $`3p`$ states shift down in energy due to the bonding between bromine and the substrate. A broadening of the Br $`3p`$ states is also observed and can be attributed to the hybridization of the bromine $`3p`$ states with the $`s`$\- and $`d`$-bands of the substrates (see Fig. 3). The Br $`3s`$ states are also seen to mix slightly with the $`s`$ and $`d`$ states of the substrate (see also Fig. 3). The hybridization of the Br $`3s`$ and $`3p`$ states with the electronic states of the substrate suggests covalent bonding between bromine and the Ag(100) and Au(100) surfaces.

The bonding of Br with the Ag(100) and Au(100) surfaces is also found to be associated with a charge transfer from the substrate to Br. To obtain a rough estimate of the charge transfer, we calculated the change of the charge for a bromine atom upon adsorption by integrating the difference of the corresponding charge densities over a sphere with a radius of 1.28 Å around the atom.  We found that 0.15 and 0.14 electrons were transfered from the Ag(100) and Au(100) surfaces, respectively, to the bromine atom. The amount of the charge transfer was found to be basically the same for the $`H_4`$, $`B_2`$, and $`T_1`$ configurations. These results are consistent with the data for the DOS projected onto the Br atom. By integrating the $`3p`$ contributions up to the Fermi level, we observe that more $`3p`$ states are occupied in the Br-adsorbed surfaces than in the systems before adsorption (see Fig. 4). A recent periodic GGA calculation with a local basis set for the adsorption of chlorine on the Ag(111) surface also found that a slight charge ($``$ 0.2 electrons) was transferred from the Ag(111) substrate to the chlorine atom.  Experimental measurements of the electrosorption valency of Br adsorbed on Ag(100) report values of approximately $``$0.70 to $``$0.75,  corresponding to a residual charge of 0.25 to 0.30 electrons on the adsorbed Br. These values are considerably larger than our calculated charge of 0.15 electrons. The discrepancy may be due to the fact that our calculations were performed for systems in vacuum. In an electrochemical environment, the net charge associated with the adsorbate might be very different from that in vacuum, due to solvation. However, the discrepancy might also be attributed to inaccuracies in the theoretical and experimental methods used to estimate the charge.

The difference in the bonding strength of the bromine with the substrate between the different configurations directly affects their relative stability. Based on our DOS data, we provide a qualitative explanation of the difference in bonding strength between the $`H_4`$ and $`B_2`$ configurations. The Br $`3s`$ and $`3p`$ states in the Br/Ag(100)-$`H_4`$ configuration are significantly lower in energy than in the Br/Ag(100)-$`B_2`$ configuration (see Fig. 4). In addition, the intensity of the lower part of the $`3p`$ states (below approximately $``$3 eV) is larger for Br/Ag(100)-$`H_4`$ than for Br/Ag(100)-$`B_2`$. Both facts suggest a stronger bonding for the $`H_4`$ configuration on Ag(100). This is expected since there are more direct bonding neighbors for the bromine atom at the $`H_4`$ site. On the other hand, the Br

$`3s`$ states in the Br/Au(100)-$`H_4`$ configuration are only slightly lower in energy than in the Br/Au(100)-$`B_2`$ configuration. While the Br $`3p`$ states extend over almost the same range in energy for the Br/Au(100)-$`H_4`$ and -$`B_2`$ configurations, they have slightly larger intensity in the lower part (below $``$ $``$2 eV) and smaller intensity in the higher part (above $``$ $``$2 eV) for the Br/Au(100)-$`H_4`$ configuration than for the Br/Au(100)-$`B_2`$ configuration. While the $`H_4`$ configuration thus has a stronger covalent bonding for both the Br-adsorbed Ag(100) and Au(100) surfaces, the difference in bonding strength between the $`H_4`$ and the corresponding $`B_2`$ configurations is smaller for the Au(100) surface than for the Ag(100) surface. This is probably due to the fact that the Au $`6s`$ and $`5d`$ electrons are more delocalized than the Ag $`5s`$ and $`4d`$ electrons, and the bonding strength is expected to be less sensitive to the bonding sites for substrates with more delocalized electrons. In addition to the stronger covalent bonding, the Coulomb attraction resulting from the charge transfer in the $`H_4`$ configuration is also stronger than the corresponding $`B_2`$ configuration for both the Br/Ag(100) and Br/Au(100) systems, due to the shorter distance between bromine and the surface in the $`H_4`$ configuration (see Table III).

It is clear that there is a delicate competition between the attractive and repulsive interactions in each configuration. In particular, the core-core repulsion between bromine and the substrate, which is irrelevant to the electronic DOS but makes a contribution to the total energy of the system, in the $`H_4`$ configuration is stronger than in the $`B_2`$ configuration. The core-core repulsive energy is calculated as an Ewald sum  (see the third column of Table V). The total energy, as determined in our DFT calculations, contains as separate parts the electronic and the core-core Coulomb contributions. Differences of the electronic and the core-core contributions to the total energy between the $`H_4`$ and the $`B_2`$ configuration are presented in Table V. We note that consideration of the electronic contributions alone does not properly address the opposite order of the total-energy difference (ie, the binding-energy difference) between the $`H_4`$ and $`B_2`$ configurations for the the Br/Ag(100) and Br/Au(100) surfaces. The core-core interactions need to be included. For both the Br/Ag(100) and Br/Au(100) surfaces, the electronic contribution favors the $`H_4`$ configuration, while the core-core contribution favors the $`B_2`$ configuration. In the Ag(100) case, the core-core energy, which is higher for $`H_4`$ than for $`B_2`$, is more than compensated by the lower electronic energy for the $`H_4`$ configuration, resulting in $`H_4`$ being the preferred bonding site for Br on Ag(100). For the Br/Au(100) system, however, the lower electronic energy for the $`H_4`$ configuration only partially compensates the higher core-core energy for Br/Au(100)-$`H_4`$. As a result, the $`B_2`$ configuration is lower in total energy than the $`H_4`$ configuration for the Br-adsorbed Au(100) surface. The small magnitudes of the total-energy differences, compared to the individual electronic and core-core contributions, strongly emphasize the need for very accurate energy calculations and careful convergence checks.

## IV Conclusions

The theoretical approach of supercell models combined with first-principles total-energy DFT pseudopotential methods has reproduced experimental measurements of preferred adsorption sites for Br-chemisorbed Ag(100) and Au(100) surfaces.

We have shown that while the hollow-site configuration is more stable on the Br/Ag(100) surface (by 210 meV/adatom over the bridge-site structure), the bridge-site configuration is more stable than the corresponding hollow-site structure by 60 meV/adatom on the Br/Au(100) surface. The calculations also predict that the one-fold on-top configuration is the least stable structure on both surfaces (560 meV and 300 meV higher than the corresponding most stable structure for the Br/Ag(100) and Br/Au(100) surfaces, respectively). Other aspects of the geometries of the Br/Ag(100) and Br/Au(100) systems have also been determined and are shown to be in excellent agreement with the available experimental data.

The bond between Br and the substrate is found to be covalent with a slight polarization due to a small charge transfer from the substrate to the bromine. The chemical bonding between Br and the substrate is shown to be stronger in the $`H_4`$ configuration than in the $`B_2`$ configuration. Compared with the Br/Ag(100) surface, however, the Br/Au(100) surface exhibits a reduced difference in the bonding strength between the $`H_4`$ and $`B_2`$ configurations. The core-core Coulomb interaction is found to be higher for the $`H_4`$ configuration than for the $`B_2`$ configuration. The detailed balance between the electronic and the core-core contributions to the total energy determines $`H_4`$ and $`B_2`$ as the the preferred bonding site on the Ag(100) and Au(100) surfaces, respectively.

Our work demonstrates that the use of extended surface models and careful convergence checks are critical for obtaining reliable information on the Br/Ag(100) and Br/Au(100) systems from ab initio calculations.

###### Acknowledgements.

We thank S. J. Mitchell and L. G. Wang for helpful discussions. We also thank S. P. Lewis, M. A. Novotny and G. Brown for comments on the manuscript. This work was supported by the National Science Foundation under grant No. DMR-9981815 and by Florida State University through the Center for Materials Research and Technology and the School of Computational Science and Technology.


text:
# Search for the Lepton Number Violating Process 𝜈̄_𝜇𝑒⁻→𝜇⁻𝜈̄_𝑒

## Abstract

The NuTeV experiment at Fermilab has used a sign-selected neutrino beam to perform a search for the lepton number violating process $`\overline{\nu }_\mu e^{}\mu ^{}\overline{\nu }_e`$, and to measure the cross-section of the Standard Model inverse muon decay process $`\nu _\mu e^{}\mu ^{}\nu _e`$. NuTeV measures the inverse muon decay asymptotic cross-section $`\sigma /E`$ to be $`(13.8\pm 1.2\pm 1.4)\times 10^{42}`$ cm2/GeV. The experiment also observes no evidence for lepton number violation and places one of the most restrictive limits on the cross-section ratio $`\sigma (\overline{\nu }_\mu e^{}\mu ^{}\overline{\nu }_e)/\sigma (\nu _\mu e^{}\mu ^{}\nu _e)`$ $`1.7\%`$ at 90% C.L. for V$``$A couplings and $`0.6\%`$ for scalar couplings.

Neutrino-lepton interactions provide an excellent tool to study the properties of the weak interaction. Such purely leptonic processes experience no interference from strong coupling terms, and thus provide a direct channel to investigate the nature of the weak force. The inverse muon decay (IMD) process:

$$\nu _\mu +e^{}\mu ^{}+\nu _e$$
(1)
allows one to make an accurate determination of the vector/axial-vector (V$``$A) nature of the weak interaction . This process is also sensitive to scalar couplings and right-handed currents.

An experiment with separate neutrino and anti-neutrino beams can search for the process:

$$\overline{\nu }_\mu +e^{}\mu ^{}+\overline{\nu }_e$$
(2)
Such an interaction is forbidden by the Standard Model, since it violates lepton family number conservation ($`\mathrm{\Delta }L_e=\mathrm{\Delta }L_\mu =2`$). Theories which incorporate multiplicative lepton number conservation , left-right symmetry , or the existence of bileptons  allow for such processes to occur.

The NuTeV neutrino experiment at Fermilab has investigated these processes in its high-energy, sign-selected neutrino beamline. Although the NuTeV inverse muon decay measurement is dominated by systematic uncertainties, the search for lepton number violation (LNV) processes is very sensitive because the relevant backgrounds are highly suppressed.

The experiment collected data during the 1996-1997 fixed target run, receiving a total of $`2.9\times 10^{18}`$ 800 GeV protons striking a BeO target. Pions and kaons produced in the interaction were focused using the Sign-Selected Quadrupole Train (SSQT) and aimed toward the NuTeV detector at a 7.8 mrad angle relative to the primary proton beam direction. The SSQT enabled the detector to be exposed to either pure neutrino or pure anti-neutrino beams. NuTeV received $`1.3\times 10^{18}`$ and $`1.6\times 10^{18}`$ protons on target for neutrino and anti-neutrino running modes, respectively. The fractional contamination from wrong-sign meson decays was below $`5\times 10^3`$ . Pions and kaons decay to neutrinos as they travel through a 440 m vacuum pipe; undecayed hadrons are filtered out in a beam dump at the end of the pipe. The neutrinos pass through about 900 m of earth berm shielding before reaching the NuTeV neutrino detector.

The NuTeV detector , located 1.4 km downstream of the primary target, consists of a segmented iron-scintillator sampling calorimeter, followed by a toroid spectrometer (see Fig. 1). The calorimeter is composed of 42 segments, each segment consisting of four 2 inch thick steel plates, two liquid scintillator counters, and one drift chamber. The calorimeter serves as a neutrino target with a fiducial mass of 350 tons. The scintillation counters measure the deposited hadronic energy and the drift chamber determine the position and direction of the outgoing muon. The toroid spectrometer uses a 15 kG toroid magnetic field to measure the charge and energy of muons exiting from the calorimeter. The toroid magnetic field is configured so as to always focus muons coming from the selected neutrino beam ($`\mu ^{}`$ for neutrinos, $`\mu ^+`$ for anti-neutrinos). The energy resolution and response of the detector is measured directly using a separate beam of hadrons, muons, and electrons at varying energies. The hadronic energy resolution of the calorimeter is $`\sigma /E=(0.024\pm 0.001)(0.874\pm 0.003)/\sqrt{E}`$, and the electromagnetic energy resolution is $`\sigma /E=(0.04\pm 0.001)(0.52\pm 0.01)/\sqrt{E}`$. The resolution of the muon energy as determined by the toroid spectrometer is $`\mathrm{\Delta }p/p=11\%`$, limited predominantly by multiple scattering.

The selection criteria for the inverse muon decay measurement and the lepton number violation search were similar, since the characteristic signatures of the processes are nearly identical. Candidate events were selected based on the following criteria: the event occurred during the beam gate, had its interaction vertex within the fiducial volume, and had a single $`\mu ^{}`$ reaching the toroid spectrometer. The muon was required to be well contained within the toroid and to have an energy between 15 and 600 GeV. The muon angle was also required to be less than 150 mrad with respect to the beam axis. To reduce the number of cosmic ray muons entering the selection sample, events which contained significant activity upstream of the reconstructed vertex were removed. The hadronic energy of the interaction was required to be less than 3 GeV. Finally, the neutrino beam running mode determined the sample into which the events were placed. For IMD, we required a $`\mu ^{}`$ in neutrino mode; for LNV candidates, a $`\mu ^{}`$ in anti-neutrino mode. For the LNV sample, we also placed an additional requirement on the $`\chi ^2`$ of the muon track within the toroid, in order to minimize events where the charge of the muon was misidentified.

Because both IMD and LNV events involve neutrino scattering from an electron, there exists a kinematic limit on the transverse momentum of the muon: $`p_t^22m_eE_\nu `$, where $`m_e`$ is the mass of the electron and $`E_\nu `$ is the neutrino energy. We therefore apply an energy-dependent requirement on the transverse momentum of the event in order to further isolate signal events. The cut requires $`p_t^2p_t^{2max}`$, where $`p_t^{2max}(0.059+E_\mu /671)`$ GeV2. This cut, which is based on Monte Carlo signal studies, was designed to retain 90% of the signal. The efficiency after all cuts for IMD events was 79.6%. The final efficiency for LNV events was 89.8% for V$``$A couplings and 79.9% for scalar couplings. Figures 2 and 3 show the $`p_t^2`$ distributions for right-sign and wrong-sign events respectively. Right-sign events are neutrino (anti-neutrino) events with an outgoing $`\mu ^{}`$ $`(\mu ^+)`$; wrong-sign events are the opposite: neutrino (anti-neutrino) events with a outgoing $`\mu ^+`$ $`(\mu ^{})`$.

For right-sign events, the primary backgrounds that enter the IMD sample come from low hadronic energy neutrino charged-current interactions in the detector. These include quasi-elastic events, resonance events, and some small fraction of deep-inelastic scattering (DIS) events with a very small momentum transfer . Understanding the background levels is essential for the IMD measurement, since the signal to background ratio for these events is one to eight. To accomplish this, we perform a full Monte Carlo simulation of low hadronic energy neutrino processes. To simulate neutrino resonance production, we used a low-$`Q^2`$ higher-twist approximation. We found this method more accurate in averaging over all low-multiplicity states than the single-pion production model from Rein and Sehgal. Nuclear effects such as Fermi motion and Pauli suppression were also applied to the Monte Carlo simulations. The Monte Carlo was absolutely normalized to data DIS events with hadronic energies above 30 GeV for each running mode. The normalization sample contained 0.83(0.25) million neutrino (anti-neutrino) interactions with a mean energy of 140(120) GeV.

The dominant systematic uncertainties for the right-sign events are related to the modeling of these low hadronic energy processes. Systematic errors include effects from muon energy and angular resolution, background cross-section uncertainties, Pauli suppression, and MC normalization. In addition, we take into account radiative correction errors which affect the IMD cross-section. A complete list of systematic errors is shown in Table I.

The validity of the background modeling was checked directly against the data by looking at the right-sign, anti-neutrino process $`\overline{\nu }_\mu +N\mu ^++N^{}`$. This particular configuration selects only background events, and thus is an ideal platform to test the data to Monte Carlo agreement and systematics. A fit to the anti-neutrino $`p_t^2`$ distribution (Fig. 2) is performed where the backgrounds are allowed to vary within the uncertainties shown in the first column of Table I. The fit gives an excellent $`\chi ^2`$/d.o.f. of 44.9/50, indicating that the background estimate agrees well with the anti-neutrino data within the systematic uncertainties.

Having verified the size and spectrum of the background, a fit to the neutrino data is performed to extract the IMD signal. The fit includes the previously mentioned backgrounds plus an IMD signal contribution with the proper $`p_t^2`$ distribution. As before, the backgrounds are allowed to vary within the uncertainties shown in the first column of Table I. As shown in Fig. 2, the data are well described by the combination of an IMD signal at low $`p_t^2`$ plus the background. From the fit, we extract a total of $`1050\pm 139`$ IMD events, where 1311 events were expected based on Standard Model predictions, taking into account radiative corrections  (see Table II).

The differential cross-section for IMD can be written as:

$$\frac{d\sigma }{dy}=\sigma _0E_\nu (1r)$$
(3)
where $`y=E_\mu /E_\nu `$, $`r=m_\mu ^2/s`$, $`\sigma _0=\frac{2m_eG_F^2}{\pi }`$, and $`s`$ is the center-of-mass energy of the interaction. For inverse muon decay, the NuTeV measurement for the IMD asymptotic cross-section ($`E_\nu m_\mu `$) is:

$$\sigma _0=(13.8\pm 1.2\pm 1.4)\times 10^{42}\mathrm{cm}^2/\mathrm{GeV}$$
(4)
where the first error is statistical, and the second is systematic. The average neutrino energy for the IMD events sampled in the NuTeV experiment is 130 GeV. This measurement is in agreement with the theoretically predicted value of $`17.2\times 10^{42}`$ cm2/GeV and is also consistent with the CHARM II measurement of $`(16.5\pm 0.9)\times 10^{42}`$ cm2/GeV .

By requiring that the muon charge not match the neutrino running mode (wrong-sign events), the analysis immediately becomes sensitive to lepton number violation. The dominant backgrounds in this case arise from beam impurities and muon charge mis-identification. Beam impurities come mainly from charmed meson decays and decays of wrong-sign hadrons produced in secondary interactions. Beam impurities constitute about 72% of the total LNV background. Charge mis-identification backgrounds are often associated with $`\delta `$-ray production or multiple scattering of the muon in the toroid spectrometer. These backgrounds can be greatly reduced by imposing quality cuts on the muon track in the toroid spectrometer. The total fraction of charge mis-identification is 0.06% for anti-neutrino running mode. This source contributes 14% of the LNV background. Finally, there exists an irreducible background from $`\overline{\nu }_ee^{}\mu ^{}\overline{\nu }_\mu `$, which contributes about 14% of the LNV background.

The generic expression for the differential LNV cross-section is given by:

$$\frac{d\sigma }{dy}=\lambda \frac{G_F^2s}{\pi }[A_Vy(yr)+A_S(1r)]$$
(5)
where $`\lambda `$ represents the strength of the interaction and $`A_V`$ and $`A_S`$ determine whether the reaction is V$``$A or scalar. Integrating over all allowed values of $`y`$, and normalizing to the Standard Model IMD cross-section, allows the LNV cross-section to be written as:

$$\frac{\sigma (\overline{\nu }_\mu e^{}\mu ^{}\overline{\nu }_e)}{\sigma (\nu _\mu e^{}\mu ^{}\nu _e)}=\lambda [A_V(\frac{1+r/2}{3})+A_S]$$
(6)
We can make a consistency check on the background estimation by looking at the neutrino process $`\nu _\mu +N\mu ^++N^{}`$. Momentum distributions of these events are shown in Fig. 3. A total of 28 data events were seen in this sample where $`23.5\pm 3.7`$ (stat. + sys.) were expected, consistent with the background estimate.

Looking in the LNV signal channel $`\overline{\nu }_\mu +e^{}\mu ^{}+\overline{\nu }_e`$ yields a total of 24 data events. A fit of the $`p_t^2`$ distribution to only background sources yields an acceptable $`\chi ^2`$/d.o.f. of 2.5/5, showing no indication of a LNV signal. Including a possible LNV signal in the fit yields a total LNV contribution of $`0.6\pm 3.1\pm 1.1`$ events for a V$``$A coupling and $`0.6\pm 3.1\pm 1.1`$ events for a scalar coupling. As shown in Table I, the LNV analysis is dominated by statistical uncertainty. These fit results can be recast in the form of 90% C.L. limits on the LNV cross-section as a function of $`(A_VA_S)/(A_V+A_S)`$, as shown in Figure 4. If we assume a pure V$``$A coupling, this yields $`\lambda 1.7\%`$ while a scalar coupling yields a limit of $`\lambda 0.6\%`$.

This limit is currently the most stringent limit obtained directly from neutrino-electron scattering. Previous results limited the pure V$``$A coupling to below 5%. The LAMPF experiment has set an earlier limit from muon decay rates at $`1.2\%`$ for pure V$``$A couplings.

In conclusion, NuTeV has performed a measurement of the inverse muon decay cross-section and a direct search for lepton number violation. The IMD asymptotic cross-section is measured to be $`(13.8\pm 1.2\pm 1.4)\times 10^{42}`$ cm2/GeV. The LNV search limits the strength of the interaction to be $`1.7\%`$ for V$``$A and $`0.6\%`$ for scalar couplings.

This research was supported by the U.S. Department of Energy and the National Science Foundation. We thank the staff of FNAL for their contributions to the construction and support of this experiment during the 1996-1997 fixed target run.


text:
# Criticality in a model of banking crises.

## I Introduction

Systemic failure in banking arises when one bank’s failure triggers off further failures. The history of modern banking is full of examples of systemic failure at both moderate and large scales, with the 1997 East Asian crisis being the most recent example of large scale bank failure.

Several channels through which systemic failure may arise have been explored . In this paper, we focus on one channel, namely interbank lending. Interbank lending allows banks which face temporary shortfalls in funds to borrow from banks which have surpluses. Such lending creates a network of credit and debt relationships within the banking system. If one debtor bank fails, it adversely affects the balance sheet of its creditors and can trigger off their subsequent failures.

At the same time, interbank lending allows for the pooling of risk arising from the stochastic pattern of deposits and withdrawals by each bank’s customers. Without interbank lending, banks might become even more vulnerable to failure, although such failure would be idiosyncratic and free of the symptoms of systemic collapse. Judgment on the effects of interbank lending for a banking system’s stability should therefore take into account both the ex ante risk-sharing effect and the ex post systemic failure effect. In this paper, we study these effects by developing and simulating a dynamic model of a banking system, linked together through an interbank credit market.

Our approach abstracts from the rationality of individual behavior and resembles the approach taken in the statistical mechanics of disordered systems. Episodes of systemic failure in banking appear analogous to the periodic collapses or ‘avalanches’ which arise in natural systems in the form of earthquakes, microfracturing, epidemics, magnetization of ferro-magnetic systems, etc. The capacity of a natural system to generate avalanches of various sizes seems to depend on certain general features, such as the heterogeneity of the elementary constituents of a system and the nature of interaction between them. By analogy in a banking system, we focus on how the individual characteristics of banks and their mutual interaction on an interbank market can affect the stability of the system.

## II  The Model.

Time is discrete. At the initial time, $`t=0`$, the system starts with $`N_0`$ banks. At any subsequent time, there are $`N_t`$ banks operating in the system. The number of banks may go down as a result of failures, but failing banks are not replaced by new entrants.

Credit linkages between banks are defined by a connectivity matrix, $`J_{ij}`$. $`J_{ij}`$ is either one or zero; a value of one indicates that a credit linkage exists between banks $`i`$ and $`j`$ and zero indicates no relationship. $`J_{ij}`$ are randomly chosen at the beginning of the simulation. $`c`$ denotes the probability that $`J_{ij}`$ is one for any two banks. At one extreme, $`c=0`$ represents the case of no interbank lending, while $`c=1`$ represents a situation in which all banks can potentially borrow and lend from each other.

The primary purpose of a bank is to channel funds received from depositors towards productive investment. Deposits, $`A_t^k`$, and investment opportunities, $`\omega _t^k`$, are randomly and independently drawn for each bank $`k`$ at each time $`t`$ from Gaussian distributions with means $`\overline{A}^k`$ and $`\overline{\omega }^k`$ respectively, and variances $`\overline{A}^k\sigma _A`$ and $`\overline{\omega }^k\sigma _\omega `$ respectively:

$$A_t^k=|\overline{A}^k(1+\sigma _Aϵ_t^1)|$$

$$\omega _t^k=|\overline{\omega }^k(1+\sigma _\omega ϵ_t^2)|$$

$$ϵ^{1,2}N(0,1)$$

Once tied up in investment, resources become relatively illiquid; any investment made at time $`t`$ fully matures only at some time $`t+\tau `$. A return, $`\rho `$, is realized at each point, $`t+1`$, $`t+2`$, $`\mathrm{}`$ $`t+\tau `$. $`\rho `$ is exogenous and risk-free. However, since withdrawals are unpredictable, a bank may find it itself unable to repay its depositors due to the illiquidity of its investment.

Each bank $`k`$ starts a period, $`t`$, with cash holdings inherited from the past, $`M_{t1}^k`$. These are constituted by

$$M_{t1}^k=A_{t1}^k+B_{t1}^k+V_{t1}^k\underset{s=1}{\overset{\tau }{}}I_{ts}^k,$$

where $`A^k`$ are its deposits, $`V^k`$ represents its equity, $`I^k`$ represents past investments and $`B^k`$ represents borrowings. $`B^k`$ can be negative or positive and satisfies: $`_{k=1}^NB^k=0`$. Borrowing consists of one-period loans, which require repayment in full in the period after which they are undertaken.

Each bank receives income, $`\rho _{s=1}^{\tau 1}I_{ts}^k`$ and $`(1+\rho )I_{t\tau }^k`$ from investments made over the last $`\tau `$ periods. It then receives a new level of $`A_t^k`$, as defined above. Its cash holdings then adjust to reflect $`(A_t^kA_{t1}^k)`$, the net change in its deposits and $`\rho _{s=1}^\tau I_{ts}^k+I_{t\tau }^k`$, the receipt of incomes from past investment. If the adjusted cash holdings become negative, the bank can issue negotiable debt certificates to cover any excess of payments over its cash reserves. However, the certificates have to be redeemed at the end of the period through borrowing from other banks. If this is not done, the bank fails and its debt certificates become worthless.

If a bank has borrowed in the past period, its priority is to repay its creditors. Such payments are made in cash. Hence, banks with cash holdings in excess of their debt obligations, which equal $`(1+r_b)B_{t1}^j`$ for bank $`j`$, make payments to their respective creditors. Here, $`r_b`$ is the exogenous interest rate on interbank borrowing.

At this point, two types of banks can be distinguished, those with positive cash and those with negative cash. Accordingly, they get classified as potential lenders and potential borrowers respectively. Borrowing banks issue demands for loans equaling their debt obligations (interest plus principal) minus their current cash (which if negative adds to their demand).

Lending banks give priority to dividend payments (to shareholders) and investment. Dividend payments are made for maintaining a target equity:deposit ratio, $`\chi `$, and for preventing excessive levels of capitalization. Only banks whose adjusted cash holdings plus illiquid assets exceed a certain fraction of their total deposits make dividend payments. After dividends have been paid, the bank is assumed to undertake investment on the basis of its available liquid resources on the one hand, and its stochastic investment opportunity on the other. Available liquid resources comprise the bank’s current cash minus any statutory reserve requirement imposed on it by a regulator. If the available investment opportunity exceeds this limit, the bank cannot exploit the full opportunity. It then invests upto the extent of its available cash resources.

After investment, any excess left over is made available to those borrowing banks with whom a credit linkage exists. Each borrowing bank contacts each lending bank with which it is linked in a random order. The two banks exchange an amount of credit equal to the minimum of the two banks’ respective demand and supply. If the borrowing bank is left with an unfulfilled trade it contacts another lender with which it is linked.

A borrowing bank does not receive actual funds until it has lined up enough credit to ensure that it will not fail during the current period. Once a bank has obtained sufficient credit, funds are transferred and the cash positions of all banks involved is updated. This continues until either all loanable funds are exhausted or all demands for credit are satisfied.

At this point, the process reiterates itself through the following steps: banks which had not repaid creditors in the first round but now have borrowed enough cash to pay off past debt entirely, do so; these payments go to their creditors and money holdings get updated accordingly; potential borrowers and lenders are determined for the next round; if lenders had underutilized investment opportunities from the previous round, they make further investment. Finally, a fresh round of borrowing and lending takes place. The process repeats itself until reiteration produces no further exchange of credit.

All banks which are left with negative cash holdings or cash holdings which fall short of their remaining debt obligations are deemed to be in default. These banks are removed from the system. If, at the time of closure, a failing bank has illiquid assets in the form of investments made previously, these are liquidated at a fraction, $`\gamma `$ (which is exogenous), of their true value. The proceeds are distributed, first to the depositors, then, if there is still some value left, this goes to creditors from the previous period and finally to the shareholders.

## III Simulations and Results.

In all the simulations to be presented, the return on investment, $`\rho `$, was risk-free and equal to $`1.0\%`$ , the interest rate on interbank borrowing, $`r_b`$, was $`0.5`$, $`\chi `$, the equity:deposit ratio was $`30\%`$, the period of maturity $`\tau `$ was set at 3 and the recovery rate $`\gamma `$ was set at 0. Parameters which vary across the reported simulations are identified separately.

At $`t=0`$ each bank holds deposits, $`A_0`$, worth 1000 units and initial equity, $`V_0`$, equal to $`0.3`$ times the initial deposit.

Initially, $`\overline{A}_k`$ and $`\overline{\omega }_k`$ were chosen as identical for each bank. This led to homogeneity in average size of deposits and investment opportunities across banks.

Figure 1 displays results on the main question of the paper. It compares bank failures with different degrees of linkage, $`c`$, in the interbank market. Increasing linkage adds stability, in the sense that in any period, there are more surviving banks the greater the degree of linkage. This pattern was very robust to changing the parameters of the simulation.

A policy issue raised by the results of Figure 1 concerns the role of statutory reserve requirements. By preventing banks from investing more than a certain fraction of the deposits placed by customers in the form of illiquid and/or risky assets, reserve requirements reduce the exposure of individual banks to the risk of failure. But they also inhibit interbank lending activity. This can reduce the risk sharing provided by interbank credit and destabilise the system. We conducted experiments along the lines of those reported in Figure 1 and found that, indeed, while without any interbank credit higher reserve requirements always led to fewer bank failures, with interbank credit linkages, similar increases in reserve requirements could increase the incidence of bank failures (see ).

In the homogeneous case depicted in Figure 1, episodes of bank failure affect only a small number of banks at a time and there appears to be no correlation across the episodes. This is because there is relatively little borrowing and lending taking place at any time. Even this small amount appears to save banks which face short run liquidity problems, but is not enough to generate contagion effects.

To allow for systematically larger volumes of interbank activity, we try to differentiate banks in a way that some banks tend to be lenders and others to be borrowers. We therefore make banks differ according to investment opportunities, by choosing: $`\overline{\omega }^k=\overline{\omega }|z^k|`$ with $`z^kN(0,\sigma _{\overline{\omega }})`$. Given that banks remained identical in terms of their deposit fluctuations, a large $`\overline{\omega }_k`$ meant that a bank would invest large amounts and consequently face a greater exposure to liquidity risk. Such a bank would be more likely to act as a large borrower on the credit market, while for analogous reasons, a bank with small $`\overline{\omega }_k`$ would be more likely to act as a large lender.

Figures 2 and 3 show some of the results from this experiment. Figure 2 shows the effect of greater heterogeneity on the incidence of bank failures when banks are 100 percent connected to each other. As the variance of types increases from 5 to 40, the incidence of failures goes up. Note that now the time path shows ‘avalanches’, i.e. many banks collapse together over brief periods of time. Since the aggregate environment is on average constant, avalanches suggest that knock-on effects are taking place.

Figure 3 shows that, fixing the variance of types at 10, increasing connectivity from 5 to 50 percent leads to less failures during the interval studied, but increasing connectivity to 80 percent leads to more failures than at 50 percent. Hence, in contrast to the case of homogeneous banks, increasing the extent of the interbank market need not increase overall stability.

We repeated the experiments discussed in Figures 2 and 3 by differentiating banks both in terms of size and in terms of investment opportunity, as follows: $`\overline{A}^k=\overline{A}|z^k|`$, and $`\overline{\omega }^k=\overline{\omega }|z^k|`$, with $`z^kN(0,\sigma _{\overline{\omega }})`$. With this formulation, all banks face the same liquidity risk, but heterogeneity exists in the volume of borrowing and lending by different banks. A big borrowing bank might need to contact many small lending banks and a big lending bank might be a source of liquidity to many small banks. In both cases, the collapse of a big bank can create ripples through the system. The simulation results verified that this indeed does happen . These results are also consistent with empirical observations made in the literature (see, e.g., ).

The results suggest that when banks are heterogeneous in the volume of their activity on the interbank market, contagion effects can arise which can undermine the overall risk sharing influence of interbank credit. It is interesting to observe that heterogeneity plays a similar role in a well known physical model, the Random Field Ising Model (RFIM). In particular, the RFIM shows a ’disorder’ induced phase transition between a regime displaying only finite size avalanches and a regime corresponding to avalanches which cover the whole system. At a critical level of disorder, avalanches of all sizes, up to and including the system as a whole, may occur. The statistical distribution of avalanches takes the form of a power-law at the critical value .

A similar phase transition is found in our model. For example, Figure 4 shows that when linkages are fixed at 100 percent and setting $`\sigma _{\overline{\omega }}=40`$, the system reaches a critical state where the size of avalanches becomes distributed according to a power-law. In order to better identify this power-law behaviour, we increased the number of banks to 1600 and repeated the simulations several times, resampling the investment opportunities each time. A preliminary investigation indicates that for each value of linkage across banks, it is possible to find a critical value of heterogeneity in the system which leads to power-law distributions of avalanches sizes. The variance of sizes needed to attain the power-law distribution increases as the connectivity decreases.

## IV Conclusions:

Our simulations have identified certain characteristics which can lead interbank lending to be associated with complex forms of instability in a banking system. Whether or not the crises which periodically grip actual banking systems represent such complex behaviour is difficult to establish empirically. The data needs for such measurement outstrip the recorded evidence on the subject. It is nonetheless useful for policy makers to be able to distinguish circumstances under which interbank credit can destabilise a system from those under which it stabilises. In this respect, our simulations suggest that when bank lending is extended across banks which are heterogeneous in size and risk exposure, it can exarcebate the potential for avalanche effects and destabilise the system.