text:
# Search for a scalar bottom quark with mass 3.5–4.5 GeV/𝑐²

## Abstract

We report on a search for a supersymmetric $`\stackrel{~}{B}`$ meson with mass between 3.5 and 4.5 GeV/$`c^2`$ using 4.52 $`\mathrm{fb}^1`$ of integrated luminosity produced at $`\sqrt{s}=10.52`$ GeV, just below the $`e^+e^{}B\overline{B}`$ threshold, and collected with the CLEO detector. We find no evidence for a light scalar bottom quark.

preprint:   CLNS 00/1696  CLEO 00-20

V. Savinov,1 T. E. Coan,2 V. Fadeyev,2 Y. S. Gao,2 Y. Maravin,2 I. Narsky,2 R. Stroynowski,2 J. Ye,2 T. Wlodek,2 M. Artuso,3 R. Ayad,3 C. Boulahouache,3 K. Bukin,3 E. Dambasuren,3 S. Karamov,3 G. Majumder,3 G. C. Moneti,3 R. Mountain,3 S. Schuh,3 T. Skwarnicki,3 S. Stone,3 J.C. Wang,3 A. Wolf,3 J. Wu,3 S. Kopp,4 M. Kostin,4 A. H. Mahmood,5 S. E. Csorna,6 I. Danko,6 K. W. McLean,6 Z. Xu,6 R. Godang,7 G. Bonvicini,8 D. Cinabro,8 M. Dubrovin,8 S. McGee,8 G. J. Zhou,8 E. Lipeles,9 S. P. Pappas,9 M. Schmidtler,9 A. Shapiro,9 W. M. Sun,9 A. J. Weinstein,9 F. Würthwein,9,***Permanent address: Massachusetts Institute of Technology, Cambridge, MA 02139. D. E. Jaffe,10 G. Masek,10 H. P. Paar,10 E. M. Potter,10 S. Prell,10 D. M. Asner,11 A. Eppich,11 T. S. Hill,11 R. J. Morrison,11 R. A. Briere,12 G. P. Chen,12 A. Gritsan,13 J. P. Alexander,14 R. Baker,14 C. Bebek,14 B. E. Berger,14 K. Berkelman,14 F. Blanc,14 V. Boisvert,14 D. G. Cassel,14 P. S. Drell,14 J. E. Duboscq,14 K. M. Ecklund,14 R. Ehrlich,14 A. D. Foland,14 P. Gaidarev,14 R. S. Galik,14 L. Gibbons,14 B. Gittelman,14 S. W. Gray,14 D. L. Hartill,14 B. K. Heltsley,14 P. I. Hopman,14 L. Hsu,14 C. D. Jones,14 J. Kandaswamy,14 D. L. Kreinick,14 M. Lohner,14 A. Magerkurth,14 T. O. Meyer,14 N. B. Mistry,14 E. Nordberg,14 M. Palmer,14 J. R. Patterson,14 D. Peterson,14 D. Riley,14 A. Romano,14 J. G. Thayer,14 D. Urner,14 B. Valant-Spaight,14 G. Viehhauser,14 A. Warburton,14 P. Avery,15 C. Prescott,15 A. I. Rubiera,15 H. Stoeck,15 J. Yelton,15 G. Brandenburg,16 A. Ershov,16 D. Y.-J. Kim,16 R. Wilson,16 T. Bergfeld,17 B. I. Eisenstein,17 J. Ernst,17 G. E. Gladding,17 G. D. Gollin,17 R. M. Hans,17 E. Johnson,17 I. Karliner,17 M. A. Marsh,17 C. Plager,17 C. Sedlack,17 M. Selen,17 J. J. Thaler,17 J. Williams,17 K. W. Edwards,18 R. Janicek,19 P. M. Patel,19 A. J. Sadoff,20 R. Ammar,21 A. Bean,21 D. Besson,21 X. Zhao,21 S. Anderson,22 V. V. Frolov,22 Y. Kubota,22 S. J. Lee,22 R. Mahapatra,22 J. J. O’Neill,22 R. Poling,22 T. Riehle,22 A. Smith,22 C. J. Stepaniak,22 J. Urheim,22 S. Ahmed,23 M. S. Alam,23 S. B. Athar,23 L. Jian,23 L. Ling,23 M. Saleem,23 S. Timm,23 F. Wappler,23 A. Anastassov,24 E. Eckhart,24 K. K. Gan,24 C. Gwon,24 T. Hart,24 K. Honscheid,24 D. Hufnagel,24 R. Kass,24 T. K. Pedlar,24 H. Schwarthoff,24 J. B. Thayer,24 E. von Toerne,24 M. M. Zoeller,24 S. J. Richichi,25 H. Severini,25 P. Skubic,25 A. Undrus,25 S. Chen,26 J. Fast,26 J. W. Hinson,26 J. Lee,26 D. H. Miller,26 E. I. Shibata,26 I. P. J. Shipsey,26 V. Pavlunin,26 D. Cronin-Hennessy,27 A.L. Lyon,27 and E. H. Thorndike27

1Stanford Linear Accelerator Center, Stanford University, Stanford, California 94309

2Southern Methodist University, Dallas, Texas 75275

3Syracuse University, Syracuse, New York 13244

4University of Texas, Austin, TX 78712

5University of Texas - Pan American, Edinburg, TX 78539

6Vanderbilt University, Nashville, Tennessee 37235

7Virginia Polytechnic Institute and State University, Blacksburg, Virginia 24061

8Wayne State University, Detroit, Michigan 48202

9California Institute of Technology, Pasadena, California 91125

10University of California, San Diego, La Jolla, California 92093

11University of California, Santa Barbara, California 93106

12Carnegie Mellon University, Pittsburgh, Pennsylvania 15213

13University of Colorado, Boulder, Colorado 80309-0390

14Cornell University, Ithaca, New York 14853

15University of Florida, Gainesville, Florida 32611

16Harvard University, Cambridge, Massachusetts 02138

17University of Illinois, Urbana-Champaign, Illinois 61801

18Carleton University, Ottawa, Ontario, Canada K1S 5B6

and the Institute of Particle Physics, Canada

19McGill University, Montréal, Québec, Canada H3A 2T8

and the Institute of Particle Physics, Canada

20Ithaca College, Ithaca, New York 14850

21University of Kansas, Lawrence, Kansas 66045

22University of Minnesota, Minneapolis, Minnesota 55455

23State University of New York at Albany, Albany, New York 12222

24Ohio State University, Columbus, Ohio 43210

25University of Oklahoma, Norman, Oklahoma 73019

26Purdue University, West Lafayette, Indiana 47907

27University of Rochester, Rochester, New York 14627

There has recently been renewed interest in the possibility of a light scalar bottom quark , which ALEPH has searched for . It has been noted that such a squark can exist in certain regions of parameter space . In addition, the production cross section of scalar quarks in $`e^+e^{}`$ annihilations, well above threshold, is $`\frac{1}{4}`$ that of spin-$`\frac{1}{2}`$ quarks of the same charge; thus $`\stackrel{~}{b}\overline{\stackrel{~}{b}}`$ production would contribute $`\frac{1}{12}`$ unit to $`R`$, the ratio of hadronic cross section to $`\mu ^+\mu ^{}`$ cross section, and such an increase cannot be ruled out by existing measurements . Indirect evidence from the measured $`B`$ semileptonic branching fraction disfavors the existence of a light $`\stackrel{~}{b}`$ .

If $`\stackrel{~}{b}`$ is the lightest supersymmetric particle (LSP), and if $`R`$-parity is conserved, then $`\stackrel{~}{b}`$ would be stable. If, instead, a scalar neutrino $`\stackrel{~}{\nu }`$ is the LSP, then $`\stackrel{~}{b}`$ will decay $`\stackrel{~}{b}c\mathrm{}^{}\stackrel{~}{\nu }`$ and/or $`\stackrel{~}{b}u\mathrm{}^{}\stackrel{~}{\nu }`$. If $`R`$-parity is violated, $`\stackrel{~}{b}`$ will decay $`\stackrel{~}{b}c\mathrm{}^{}`$ and/or $`\stackrel{~}{b}u\mathrm{}^{}`$. We have searched for a light $`\stackrel{~}{b}`$ that decays $`\stackrel{~}{b}c\mathrm{}^{}\stackrel{~}{\nu }`$ and/or $`\stackrel{~}{b}c\mathrm{}^{}`$. Such a particle would dress itself as a supersymmetric $`\stackrel{~}{B}`$ meson. The dressed decays would be $`\stackrel{~}{B}DX\mathrm{}^{}\stackrel{~}{\nu }`$ and $`\stackrel{~}{B}DX\mathrm{}^{}`$, where $`X`$ represents possible additional conventional hadrons.

The decays we search for, characterized by leptons and charmed mesons, have much in common with conventional $`B`$ decays. We perform a direct search, avoiding the $`B`$ background by using a data sample collected below the $`B\overline{B}`$ threshold, at $`\sqrt{s}=10.52`$ GeV. Our search covers the $`\stackrel{~}{B}`$ mass range 3.5–4.5 GeV/$`c^2`$. Because our search is near $`\stackrel{~}{B}\overline{\stackrel{~}{B}}`$ threshold, the production cross section cannot be predicted with great precision. We include the $`\beta ^3`$ threshold factor and ignore strong interaction effects to obtain $`\sigma (e^+e^{}\stackrel{~}{b}\overline{\stackrel{~}{b}})=N_c\pi \alpha ^2Q_f^2\beta ^3/3s`$ (derived from Eq. 35.12 of Ref. ), where $`Q_f`$, the charge of the $`\stackrel{~}{b}`$, is $`\frac{1}{3}`$.

Below $`B\overline{B}`$ threshold, the major source of leptons and charmed mesons is $`e^+e^{}c\overline{c}`$, which predominantly produces two charged leptons when both charmed mesons decay semileptonically. To suppress this $`c\overline{c}`$ background, we search for events with two oppositely charged leptons as well as a fully reconstructed hadronic $`D`$ or $`D^{}`$ meson decay, where $`D^{()}`$ denotes either a $`D^{()0}`$ or $`D^{()+}`$ meson. Other sources of leptons include kaons that decay in flight, photon conversions, and $`\pi ^0`$ Dalitz decays. A search for wrong-sign $`D^{()}`$-lepton combinations was also conducted, but the resultant upper limits were significantly weaker than for the $`D^{()}`$-dilepton signature, and they are not discussed further in this paper.

The data sample used in this analysis was produced in symmetric $`e^+e^{}`$ collisions at the Cornell Electron Storage Ring (CESR) and collected with the CLEO detector in two configurations, known as CLEO II  and CLEO II.V . It comprises 4.52 $`\mathrm{fb}^1`$ of integrated luminosity produced at an $`e^+e^{}`$ center-of-mass energy of $`\sqrt{s}=10.52`$ GeV, 35 MeV below the $`B\overline{B}`$ threshold. Analysis of an additional 2.24 $`\mathrm{fb}^1`$ collected on the $`\mathrm{{\rm Y}}(4S)`$ resonance provides a sample of $`B`$ mesons, which can mimic the behavior of $`\stackrel{~}{B}`$ mesons, thus allowing us to verify our experimental technique. The response of the experimental apparatus is studied with a GEANT-based  simulation of the CLEO detector, where the simulated events are processed in a fashion similar to data.

In CLEO II, the momenta of charged particles are measured with a tracking system consisting of a six-layer straw tube chamber, a ten-layer precision drift chamber, and a 51-layer main drift chamber, all operating inside a 1.5 T superconducting solenoid. The main drift chamber also provides a measurement of specific ionization energy loss, which is used for particle identification. For CLEO II.V, the six-layer straw tube chamber was replaced by a three-layer double-sided silicon vertex detector, and the gas in the main drift chamber was changed from an argon-ethane to a helium-propane mixture. Photons are detected with a 7800-crystal CsI electromagnetic calorimeter, which is also inside the solenoid. Proportional chambers placed at various depths within the steel return yoke of the magnet identify muons.

Charged tracks are required to be well-measured and to satisfy track quality requirements based on the average hit residual and the impact parameters in both the $`r`$-$`\varphi `$ and $`r`$-$`z`$ planes. Muon candidates must penetrate the steel absorber to a depth of at least five nuclear interaction lengths, which effectively places a lower bound on the muon momentum of 1.2 GeV/$`c`$. Electrons are identified by a likelihood that includes the fraction of the particle’s energy deposited in the calorimeter and the spatial distribution of the deposited energy. To reduce contamination from low-momentum hadrons, the electron momentum is required to be greater than 1.0 GeV/$`c`$. For electrons that are combined with another lepton satisfying the above criteria, the momentum requirement is lowered to 600 MeV/$`c`$. Electron pairs from photon conversions are rejected by requiring the dielectron invariant mass to be greater than 200 MeV/$`c^2`$. We rely on the detector simulation to estimate the remaining contribution of hadrons to the lepton sample. $`\pi ^0`$ candidates are formed from pairs of photons with an invariant mass within 2.5 standard deviations of the known $`\pi ^0`$ mass that are kinematically fitted with the mass constrained to the known $`\pi ^0`$ mass.

$`D`$ mesons are reconstructed in the modes $`D^0K^{}\pi ^+`$ and $`D^+K^{}\pi ^+\pi ^+`$, where a sum over charge conjugate modes is implied. The daughters of the $`D`$ candidates must undergo ionization energy loss consistent with the particle hypothesis at the level of three standard deviations. To maximize detection efficiency, no requirements are placed on the momentum of the $`D`$ candidate. We reconstruct $`D^{}`$ mesons in the modes $`D^0D^0\pi ^0`$, $`D^+D^0\pi ^+`$, and $`D^+D^+\pi ^0`$. When used in the $`D^{}`$ modes, $`D`$ candidates must have invariant masses less than 15 MeV/$`c^2`$ (about 2.5 standard deviations) from the known mass.

Figure 1 shows the $`K^{}\pi ^+`$ and $`K^{}\pi ^+\pi ^+`$ invariant mass distributions of $`D^0`$ and $`D^+`$ candidates, respectively, associated with two oppositely charged leptons. Also shown are the normalizing distributions for $`D`$ candidates paired with a single lepton of either charge. The analogous distributions for the $`D^{}`$ modes, showing the $`D^{}`$-$`D`$ mass difference, are given in Figure 2. The $`D^{()}`$-dilepton distributions reveal a striking absence of signal. The $`D^{()}`$ yields are extracted from a fit of each histogram to a Gaussian distribution over a linear background for the $`D`$ modes and a quadratic background for the $`D^{}`$ modes. The means and widths of the Gaussians for the $`D^{()}`$-dilepton fit are fixed to values determined from the normalizing distributions, where the Gaussian widths are typically 6.5–7.0 MeV/$`c^2`$ for $`M(D)`$ and 0.5–1.0 MeV/$`c^2`$ for $`M(D^{})M(D)`$. Table I lists the yields observed in data as well as the expected yields determined from a simulation of $`e^+e^{}q\overline{q}`$ events, where $`q\{u,d,s,c\}`$, hereafter referred to as “continuum” events. No significant excess is observed in any mode.

In the normalizing distributions there is good agreement between the data and simulated continuum events. In addition, by analyzing data taken on the $`\mathrm{{\rm Y}}(4S)`$ resonance, we observe the closely related semileptonic $`B`$ decay with yields that agree well with predictions from the detector simulation.

We determine the efficiency for detecting a supersymmetric $`\stackrel{~}{B}`$ meson using Monte Carlo simulation. A range of values for the $`\stackrel{~}{B}`$ mass, $`M(\stackrel{~}{B})`$, from 3.5 to 4.5 GeV/$`c^2`$, has been explored. Also, the sneutrino mass, $`M(\stackrel{~}{\nu })`$, was varied within the kinematically allowed range. In the process $`e^+e^{}\stackrel{~}{b}\overline{\stackrel{~}{b}}`$, followed by hadronization, the supersymmetric $`\stackrel{~}{B}`$ mesons have energies $`E(\stackrel{~}{B})`$ distributed from $`M(\stackrel{~}{B})`$ to $`E_{\mathrm{beam}}`$ according to some unknown fragmentation function, which we approximate by a delta function $`\delta (E(\stackrel{~}{B})E_0)`$, with $`E_0`$ given a value between $`M(\stackrel{~}{B})`$ and $`E_{\mathrm{beam}}`$. We simulate both $`\stackrel{~}{B}`$ mesons and their daughters, as well as the additional hadrons that result from fragmentation. We determine the dependence of efficiency on $`M(\stackrel{~}{B})`$, $`M(\stackrel{~}{\nu })`$, and $`E_0`$. To simulate the decay $`\stackrel{~}{B}D\mathrm{}^{}X`$, we have taken $`X`$ to be a single pion and used three-body phase space for the decay. For $`\stackrel{~}{b}c\mathrm{}^{}\stackrel{~}{\nu }`$, we have used $`D\mathrm{}^{}\stackrel{~}{\nu }`$ three-body phase space.

These variations affect the $`\stackrel{~}{B}`$ meson detection efficiency primarily through the lepton momentum spectrum. Table I lists a representative set of detection efficiencies and the resultant 95% confidence level upper limits on the cross section, assuming $`M(\stackrel{~}{B})=4.0`$ GeV/$`c^2`$, $`M(\stackrel{~}{\nu })=0`$ GeV/$`c^2`$, and $`E_0=E_{\mathrm{beam}}`$. The efficiency exhibits weak dependence on $`M(\stackrel{~}{B})`$, but lower values of $`E_0`$ tend to soften the lepton momentum spectrum, resulting in reduced efficiency. For the smallest value of $`E_0`$ we considered (3.7 GeV), the efficiencies are one-fourth of those with $`E_0=E_{\mathrm{beam}}`$. We give the upper limit on the product of the $`\stackrel{~}{B}\overline{\stackrel{~}{B}}`$ production cross section and the $`\stackrel{~}{B}D^{()}\mathrm{}\{\pi ,\stackrel{~}{\nu }\}`$ branching fraction as a function of $`M(\stackrel{~}{B})`$ and $`E_0`$, assuming $`M(\stackrel{~}{\nu })=0`$ GeV/$`c^2`$. For all modes, the variation of the efficiency with $`M(\stackrel{~}{\nu })`$ is roughly linear between $`M(\stackrel{~}{\nu })=0`$ GeV/$`c^2`$ and $`M(\stackrel{~}{\nu })_{\mathrm{max}}0.6M(\stackrel{~}{B})0.8`$ GeV/$`c^2`$, where the efficiency vanishes. Hence, the upper limits for a finite $`M(\stackrel{~}{\nu })`$ would be those for $`M(\stackrel{~}{\nu })=0`$ GeV/$`c^2`$, scaled by a factor $`\left(1\frac{M(\stackrel{~}{\nu })}{M(\stackrel{~}{\nu })_{\mathrm{max}}}\right)^1`$.

The detection efficiency is also affected by the additional particles produced in fragmentation. From Monte Carlo studies, we have determined that for each fragmentation particle, the fractional decrease in efficiency is approximately 3% and depends on the momentum and species of the fragmentation particles. By considering Monte Carlo simulation of $`e^+e^{}c\overline{c}`$, we estimate roughly two fragmentation particles per GeV of fragmentation energy. The detection efficiency has been corrected accordingly. A systematic error per fragmentation particle of $`\pm 2.5\%`$ for the $`D`$ modes and $`\pm 5\%`$ for the $`D^{}`$ modes has been included in the upper limits to allow for the uncertainty in this correction.

In addition, the upper limits have been inflated to account for uncertainties in the detector simulation and in the $`D`$ and $`D^{}`$ branching fractions, which amount to systematic errors of 20%. We also assign an error of 10% to the expected Standard Model yields due to undertainties in the simulation of both fake and real leptons. The total systematic uncertainties are small compared to the statistical errors on the background-subtracted yields.

Figure 3 shows the 95% confidence level upper limit on the product of $`\sigma (e^+e^{}\stackrel{~}{B}\overline{\stackrel{~}{B}})`$ and $`(\stackrel{~}{B}D^{()}\mathrm{}\{\pi ,\stackrel{~}{\nu }\})`$ as a function of $`M(\stackrel{~}{B})`$ and $`E_0`$, assuming $`M(\stackrel{~}{\nu })=0`$ GeV/$`c^2`$. The strong rise of the upper limits with decreasing $`E_0/E_{\mathrm{beam}}`$ reflects the softening of the lepton spectrum. Because of the $`\beta ^3`$ dependence of the cross section near threshold, the cross section at $`\sqrt{s}=10.52`$ GeV is heavily dependent on the assumed $`M(\stackrel{~}{B})`$. These predicted cross sections are shown as the thick curves in Figure 3. For both the $`D`$-dilepton and $`D^{}`$-dilepton signatures, the upper limits in nearly all of the kinematically allowed $`M(\stackrel{~}{B})`$-$`E_0`$ parameter space with $`M(\stackrel{~}{\nu })=0`$ GeV/$`c^2`$ fall below the predicted $`\stackrel{~}{b}\overline{\stackrel{~}{b}}`$ production cross section. Below $`M(\stackrel{~}{B})=3.9`$ GeV/$`c^2`$, the $`D^{}`$-dilepton efficiency vanishes, so we exclude a $`\stackrel{~}{B}`$ that decays $`\stackrel{~}{B}D^{}\mathrm{}\{\pi ,\stackrel{~}{\nu }\}`$ only in the mass range 3.9–4.5 GeV/$`c^2`$. No portion of the parameter space is excluded if $`M(\stackrel{~}{\nu })`$ is greater than 1.2 GeV/$`c^2`$ and 1.3 GeV/$`c^2`$ for the $`D`$-dilepton and $`D^{}`$-dilepton signatures, respectively.

If the $`R`$-parity-violating decay $`\stackrel{~}{b}c\mathrm{}^{}`$ were to occur, then some fraction of the time it would appear as the dressed decay $`\stackrel{~}{B}D\mathrm{}^{}`$ or $`\stackrel{~}{B}D^{}\mathrm{}^{}`$. We have examined the $`D\mathrm{}^{}`$ and $`D^{}\mathrm{}^{}`$ invariant mass distributions in $`D^{()}\mathrm{}^+\mathrm{}^{}`$ events for evidence of a peak indicating such two-body decays. We find no evidence of a peak. Fitting the distributions to a polynomial background plus a Gaussian with width given by our experimental resolution (10 MeV/$`c^2`$), and stepping the Gaussian mean over the mass range, we obtain upper limits on the product of the $`\stackrel{~}{b}\overline{\stackrel{~}{b}}`$ production cross section and the $`\stackrel{~}{B}D^{()}\mathrm{}^{}`$ decay branching fraction, as a function of $`M(\stackrel{~}{B})`$. These upper limits, which include 15% systematic errors on the yields, are shown in Figure 3. For $`\stackrel{~}{B}D\mathrm{}^{}`$, the upper limit is less than 0.3 pb for all masses except near 4.36 GeV/$`c^2`$, where it is 1.0 pb. For $`\stackrel{~}{B}D^{}\mathrm{}^{}`$ the upper limit is weaker, typically 3 pb for masses around 4 GeV/$`c^2`$, dropping to 0.3 pb for masses near 4.5 GeV/$`c^2`$.

In conclusion, we have searched for associated $`\stackrel{~}{B}\overline{\stackrel{~}{B}}`$ production in $`e^+e^{}`$ collisions with center-of-mass energy below the $`B\overline{B}`$ production threshold. We assume a branching fraction of 100% for $`\stackrel{~}{B}D^{()}\mathrm{}^{}\pi `$ or $`\stackrel{~}{B}D^{()}\mathrm{}^{}\stackrel{~}{\nu }`$. Considering $`D`$-dilepton and $`D^{}`$-dilepton combinations, we find no evidence of a light scalar bottom quark produced at $`\sqrt{s}=10.52`$ GeV. Upper limits on the $`\stackrel{~}{B}\overline{\stackrel{~}{B}}`$ production cross section depend on the assumed mass and energy of the $`\stackrel{~}{B}`$ meson, as well as the mass of the sneutrino. For $`M(\stackrel{~}{\nu })`$ less than $`𝒪(1\mathrm{GeV}/c^2)`$, the existence of a light scalar bottom quark with mass between 3.5 GeV/$`c^2`$ and 4.5 GeV/$`c^2`$ has been excluded at the 95% confidence level. A light scalar quark decaying 100% of the time $`\stackrel{~}{b}c\mathrm{}\stackrel{~}{\nu }`$ and/or $`\stackrel{~}{b}c\mathrm{}`$ would have escaped our notice only if its decay matrix element results in a lepton spectrum much softer than three-body phase space.

###### Acknowledgements.

We thank Tung-Mow Yan and Matthias Neubert for useful discussions. We gratefully acknowledge the effort of the CESR staff in providing us with excellent luminosity and running conditions. I.P.J. Shipsey thanks the NYI program of the NSF, M. Selen thanks the PFF program of the NSF, A.H. Mahmood thanks the Texas Advanced Research Program, M. Selen and H. Yamamoto thank the OJI program of DOE, M. Selen thanks the Research Corporation, F. Blanc thanks the Swiss National Science Foundation, and H. Schwarthoff and E. von Toerne thank the Alexander von Humboldt Stiftung for support. This work was supported by the National Science Foundation, the U.S. Department of Energy, and the Natural Sciences and Engineering Research Council of Canada.


format: text

text:
# Discretized Diffusion Processes.

\[

## Abstract

We study the properties of the “Rigid Laplacian” operator, that is we consider solutions of the Laplacian equation in the presence of fixed truncation errors. The dynamics of convergence to the correct analytical solution displays the presence of a metastable set of numerical solutions, whose presence can be related to granularity. We provide some scaling analysis in order to determine the value of the exponents characterizing the process. We believe that this prototype model is also suitable to provide an explanation of the widespread presence of power-law in social and economic system where information and decision diffuse, with errors and delay from agent to agent.

\]

Most equations in science are continuous in value. While this is a good approximation to reality, natural processes often are discrete in microscopic details (e.g. atoms and molecules). We ask the following general questions: if the variables are subject to some small threshold such that no changes smaller than such a threshold are implemented, how do the continuous equations behave? Can we properly recover the continuous limit letting the threshold vanish? Are there new features emerging independent on the threshold value (and therefore robust down to the continuous limit)?

This problem is somewhat general and this is the reason of its importance. If we consider in particular the class of self-organized critical systems then the question becomes crucial. We know that for many SOC systems, such as the BTW  or the Zhang model , threshold and discretisation play a dramatic role in order to avoid the diffusive behaviour and to attain the critical one. By considering the continuous limit of the equations describing the microscopic dynamics of these cellular automata and by not preserving a finite threshold we get a trivial diffusive behaviour instead of the critical one.

In order to investigate this kind of subject, we study here one of the oldest and most important equations of physics: the Laplace equation $`^2\varphi (x)=0`$. Its time dependent version, describes the diffusion of heat and of particles, and the relaxation of incompressible, non-viscous fluids. In imaginary-time it corresponds to the Schroedinger equation of quantum mechanics for a free particle (the difference vanishes completely if we are interested in stationary solutions/eigenstates). Therefore, we do not overestimate its importance if we consider it to be ubiquitous in physics. More recently, it was also considered to play a major role in fractal growth processes such as Diffusion Limited Aggregation (DLA) and the Dielectric Breakdown Model (DBM) .

It is thus not a surprise that many methods have been developed to solve it, at least numerically. In general, both time and space are discretized on a lattice. The iterative method is the most commonly used method of solution, especially for problems such as DLA and DBM where the boundary conditions change in time. Essentially, the method resorts to iterating the equation

$$\varphi (i,t+1)=\frac{1}{2d}\underset{jn.n.i}{}\varphi (j,t),$$
(1)
whose fixed point correspond to the discretized version of $`^2\varphi (x)=0`$. Numerically convergence to the fixed values $`\varphi (i)=\varphi (i,t=\mathrm{})`$ has to be defined by some error: when two successive configurations differ by less than some small $`\delta `$ then the iteration stops. Of course $`\delta `$ has to be greater than the machine precision, $`\delta ^{}`$. This means that when the field $`\varphi `$ is of order of unity, it is impossible to change its value of quantities less than $`\delta ^{}`$.

Therefore it is of extreme importance to explore what happens if we explicitely set a precision in the definition of $`\varphi `$, that is, if we say that we can change its value, according to (1), only by integer multiples of some chosen $`\epsilon `$. This corresponds to consider the fields $`\stackrel{~}{\varphi }`$

$$\stackrel{~}{\varphi }(i)=Int_\epsilon \left[\frac{1}{2d}\underset{jn.n.i}{}\stackrel{~}{\varphi }(j)\right]$$
(2)
$`Int_\epsilon [x]`$ corresponds to taking the integer multiple of $`\epsilon `$ closest to $`x`$ (other possible definitions of the integer part, such as the lowest integer closest to $`x`$, have been explored and have been shown to give the same results). The relevance of this problem is not only related to the numerical solution of Laplacian equation, but it comes also from real problems: granular materials, for example, behave sometimes as fluids. As such, we can expect them to obey some of the laws typical of fluids (hence also to relax according to the time dependent Laplacian equation), yet their intrinsic granularity forbids movements of quantities smaller than a single grain. Under this respect we refer to the properties of $`\stackrel{~}{\varphi }`$ as properties of rigid diffusion processes. Discrete-discrete Laplacian fields $`\stackrel{~}{\varphi }`$ (discrete in space-time, and discrete in the field $`\stackrel{~}{\varphi }`$) can also play a role in modeling of economic and social systems. In particular, for the simplest case of the Laplacian “harmonic” operator the field is computed by averaging over the field neighbours. If this is viewed as a social system where opinions of players are formed by consulting people around, we would like to introduce a more realistic version, where decisions are considered only when a certain threshold of information is overcome and the diffusion is therefore discretized. This effect known as “rigidity” in the literature is indeed responsible of the onset of criticality. Since Self-Organized Criticality is believed to play a role both in economic and social system we suggest that our results can be viewed as a way to describe the appearance of criticality in such environments.

In this Letter we explore the behavior of this Rigid Laplacian Model (RLM). As we shall see, non trivial properties emerge already in one dimension.

The standard setup of our simulations is the following: we take a $`1d`$ lattice of $`L+1`$ sites, and we set the boundary conditions $`\varphi (0)=1`$ and $`\varphi (L)=0`$. Then we iterate the equation

$$\varphi (i,t+1)=Int_\epsilon \left[\frac{1}{2}(\varphi (i1,t)+\varphi (i+1,t))\right]$$
(3)
The iteration stops when there are no more possible rearrangements. The exact solution of $`^2\varphi =0`$ on the lattice with the above boundary conditions is the straight line $`\varphi (i)=i/L+1`$. As we show in Fig.1, the typical solution of (3) is instead a quadratic curve.

This can be easily explained by considering that as soon as a certain precision $`\epsilon `$ is introduced, the Laplace equation can be considered as a first approximation as

$$^2\varphi =\alpha \epsilon $$
(4)
In the case where the boundary condition are $`\varphi (L)=0`$ and $`\varphi (0)=1`$ we obtain the following equation

$$\varphi _\epsilon (x)=\frac{1}{2}\alpha \epsilon L^2\left(\frac{x}{L}\right)^2\left(1+\frac{1}{2}\alpha \epsilon L^2\right)\frac{x}{L}+1,$$
(5)
in terms of the reduced variables $`x^{}=x/L`$ and $`\epsilon ^{}=\frac{1}{2}\alpha \epsilon L^2`$ one obtains

$$\varphi _\epsilon ^{}(x^{})=\epsilon ^{}x^2(1+\epsilon ^{})x^{}+1.$$
(6)
The precise value of $`\alpha `$ depends on the initial condition and it has to be determined with some fitting procedure. We can study the dynamics of the system observing that the evolution (1) always minimizes the functional

$$E=\underset{i=0}{\overset{L1}{}}(\varphi (i+1)\varphi (i))^2$$
(7)
Eq.(7) maps the RLM into the discrete space-time Gaussian model. Physically, since $`\varphi `$ represents an electrostatic field, $`E`$ represents the total electrostatic energy of the sample, whose minimization in the case of fractal growth has been recently studied. Letting the field $`\varphi `$ take on only integer values leads to the discrete Gaussian model. The dynamics described by (3) is therefore a $`T=0`$ dynamics, always decreasing the energy, and the stability of solutions different from the constant slope implies the presence of local energy minima. Since this process of energy minimization stops in local minima, we can perform some sort of simulated annealing by perturbing it. We increase the value of the field $`\varphi `$ in some randomly chosen point by $`n\epsilon `$, then we let the system relax. Indeed, we find that the curvature $`\alpha `$ decreases in time (time is defined in terms of perturbation steps; we consider the relaxation process to be much faster, as usual in self-organized critical models). Each perturbation implies a rearrangement of the profile.

For the energy of the system one has from the eq.(5)

$$EL^3\alpha ^2\epsilon ^2.$$
(8)
We know that for ordinary diffusion the energy variation $`E`$ is proportional to $`t^{3/2}`$. We will show in the conclusion that based on that result, in this case one can expect $`Et^{1/2}`$ or more precisely $`EE_0(t/L^z)^{1/2}`$, where we explicitely introduced a dynamical exponent $`z`$ for the process. Assuming for the moment this result, one has

$$EE_0(t/L^z)^{1/2}(t/L^z)^{1/2}E_0L^3\alpha ^2(t/L^z)$$
(9)
By derivating this expression and by requiring that it shows no dependence on $`L`$ (unless for the ratio $`t/L^z`$) one obtains that the dynamical exponent must be $`z=3`$, in agreement with the value $`z=3.0\pm 0.1`$ from numerical simulations, and

$$\alpha (t,L)=\alpha _0(t/L^z)^{1/2}f(t/L^z).$$
(10)
the scaling function $`f(x)`$ is such that $`\alpha \alpha _{\mathrm{}}0`$ for $`tL^z`$. Such behaviour has been tested in Fig.2 We see that for this kind of systems the relaxation process to the equilibrium is very slow because of granularity. Indeed, for a continuous Laplacian field the relaxation to the stationary state is an exponential process with a characteristic time scale $`t_0L^2`$.

Such rearrangements proceed through power-law distributed avalanches, whose measure is the number of sites that changed their value of $`\varphi `$. We show in Fig.3 the probability $`p(s)`$ to deal with an avalanche of size $`s`$. This function can be fitted with a power law of the kind $`p(s)s^\tau `$. The value of the exponent $`\tau `$ measured from the data is $`3.8\pm 0.3`$. To take in account the finite size effects of this distribution (the size of an avalanche cannot be greater than the size of the system) we can write the universal form $`p(s,L)=s^\tau g(s/L)`$, where $`g(x)`$ is a scaling function such that $`p(s=L)=0`$. This scale invariant behavior describes the dynamics in a transient period where perturbations are accumulated into the system, until a limit profile (not necessarily the straight one) is reached. This transient state can be suitably delayed by considering a different dynamics where also negative perturbations are considered. To describe the evolution to the stationary state, we measured the average size $`s`$ of the avalanches with respect to the time of the simulation.

In Fig.4 we present the behaviour of $`s(t)`$ that also behaves as a power law with exponent $`1/3`$. A simple scaling argument can be presented to explain such behaviour. Indeed,

$$s(t,L)=Lt^\nu g(t/L^z)$$
(11)
where we explicitely considered the finite size effects in the average size. In the above expression, $`\nu `$ represents the expected exponent and $`z`$ is the dynamical exponent previously introduced and present also in the scaling function $`g`$.

The collapse shown in the inset of Fig.4 has been obtained with the following values: $`z=3.0\pm 0.1`$, $`\nu =0.34\pm 0.01`$. The values of these two exponents are related, indeed the meaning of the dynamical exponent is as usual to determine the scale of time after which the process stops. By imposing that at the maximum time the average size of the avalanche is independent on $`L`$, one obtains the scaling relation

$$1z\nu =0$$
(12)
from which $`\nu =1/z=1/3`$ as numerically found in the limit $`t<<L^3`$.

From this scaling relation it’s possible to provide an argument to determine the value of the exponent $`\tau `$. Since the number of avalanches is the same whether they are classified according to their sizes (i.e. by means of the $`p(s)ds`$) or just counted in time (i.e. by means of $`dt`$) then one obtains $`p(s)=\frac{dt}{ds}`$. From the above scaling laws one finally obtains

$$\tau =1+1/\nu .$$
(13)
We numerically tested the validity of these relationships as above mentioned. The “rigidity” of the field makes the process sub-diffusive as it is implied by the value $`z=3`$ of the dynamical exponent, as opposed to the usual $`z=2`$ value. A more intriguing feature of these processes can be obtained by a suitable rescaling of the microscopic time scale. As a matter of fact one can check the number $`n(t)`$ of microscopical sweeps needed in order to update the values of the field $`\varphi `$ in an avalanche. This number $`n(t)`$ must be proportional to the average size of an avalanche at time $`t`$ and then it varies as $`t^{1/3}`$, as it can be shown by a direct measurement (Fig.5). It means that the “physical” time of the process defined as $`_{t^{}=1}^tn(t^{})`$ behaves as $`t^{2/3}`$. By rescaling the scale of time with this characteristic time scale, the process loses any information on the granular nature of the medium and the dynamics behaves as a usual diffusive process. In this way we demonstrate the initial assumption for the energy $`E`$ in the system. In fact, since it varies for ordinary diffusion phenomena as $`Et^{3/2}`$ and in this case the time of the process $`t^{}`$ varies as $`t^{}t^{1/3}`$ we obtain the above assumption of $`Et^{1/2}`$.

We also checked the properties of this model in dimension larger than $`d=1`$. We found a similar scale invariant behaviour. In the case of $`d=2`$ we found a value of $`\tau =2.65\pm 0.03`$, for $`d=3`$ we have $`\tau =2.60\pm 0.05`$ and for $`d=4`$ $`\tau =2.55\pm 0.05`$.

As described in the introduction, it would be tempting then to link this randomness to the known effect of “rigidity” in Self-Organized systems. Self-Organization has often been invoked as a way to describe the invariance observed in various social and economic systems. Indeed we propose the basic mechanism of diffusion of money and/or information in a random environment (that can stop or reduce the information interchanged between players) as one of the reasons that could lead to scale invariance and to avalanche dynamics. As regards the economic systems, we believe that this study could be used also to describe the time behaviour of a quantized field with respect to the value of the truncation error. A specific example in this field is be presented by considering the spurious effect in the price dynamics that could happen by changing the tick at which they are traded (i.e. from $`1/16`$ of dollar to $`1/100`$).

In conclusion, we introduced a model of diffusion, where the redistribution process is hindered by some sort of rigidity induced by the granularity of the fields. Randomness is present through a cutoff in the different possible values the field can assume. By truncating to a certain precision the value of the field a new and rich scenario of metastable states appears. The evolution of this system through the metastable states is driven by an avalanche dynamics with no particular time or length scale. By suitably rescaling the microscopic time of the evolution, one can use the scaling relations of diffusive process to give a theoretical ansatz for the quantities describing the system. From computer simulations we find a good agreement between numerical data and theoretical ones. We believe this sort of process could be responsible for the ubiquitous presence of power-law relations in the every day life. Clearly our work is readily generalised to other systems, notably Naviers-Stokes, KPZ, wave equations. One may wonder if there is some sort of super-universality in their scaling behavior. Indeed a long list of similar problems awaits us to study.

P. De Los Rios and Y-C. Zhang thank the Dipartimento di Fisica in the University of Rome “La Sapienza”, for its kind hospitality. G. Caldarelli thanks for the same reason the Institut de Physique Théorique in Fribourg. This work has been supported by the European Network contract FMRXCT980183.


text:
# Acknowledgments

## Acknowledgments

This work was partially financed by CNPq, BRAZIL.

## Tables

## Figure Captions

Fig. 1: The logarithmic $`Q^2`$ slope of $`F_2^A`$ as function of the variable $`Q^2`$ at different values of $`A`$ and $`W`$. See text.

Fig. 2: The logarithmic $`Q^2`$ slope of $`F_2^A`$ as function of the variable $`x`$ at different values of $`A`$ and $`W`$. See text.


text:
# Scattering properties of a cut-circle billiard waveguide with two conical leads

## I Introduction

The scattering properties of ballistic electron waveguides have been studied by a number of authors, but primarily in the semi-classical regime, where resonances strongly overlap , . Several experiments  have shown that the shape of the waveguide cavity can strongly affect the conductance, depending on whether the shape induces integrable or chaotic motion . In this paper, we investigate the scattering properties of an electron waveguide whose scattering cavity strongly resembles that of a quantum dot. We investigate this system in a regime where the resonances are well separated in energy, and we look at the effect of distortions in shape and changes of lead placement, which change the symmetry of the scattering cavity, on the behavior of the scattering resonances.

The geometry of the waveguide we consider in this paper is shown in Fig. 1. The scattering system consists of a cavity with two cone shaped leads attached, which we will treat as infinitely long. The walls of both the cavity and the leads are assumed to be infinitely hard (infinite potential) walls. The cavity has the shape of a circle of radius $`R`$, but with a segment of the wall, subtended by an angle $`\omega `$, replaced by a flat segment which we call the cut. The opening angles of the leads are $`\mathrm{\Delta }\theta _1`$ and $`\mathrm{\Delta }\theta _2`$. The angular positions of the center of the cut and the center of the lead openings are denoted by $`\zeta `$, $`\theta _1`$, and $`\theta _2`$, respectively.

The effect of conical leads has already been studied in earlier works (e.g., Ree and Reichl , Persson et al. , and Berggren, Ji, and Lundberg ). Compared to straight leads (which have most commonly been used in publications on waveguides), conical leads are more similar to the electron sea which provides the source of electrons in most experiments. As shown in Ref. , conical leads allow tunneling resonances through the waveguide cavity for energy regimes in which conduction is prohibited for straight leads.

A classical particle moving in a closed billiard whose shape is that of a circle with a cut (the shape of our cavity) will display a rich range of chaotic behavior as the size of the cut is varied . Quantum signatures of the classical chaos in this billiard have been studied in Ref. . One question we will be interested in is how well the scattering process can probe this dynamics of the closed cavity. In the low-energy regime we will consider, we expect that the waveguide will exhibit well-separated Fano resonances at energies close to the energy eigenvalues of the closed billiard.

Beyond simply allowing to probe eigenenergies of the closed billiard as scattering resonances, lead placement itself can alter the symmetry properties of the closed billiard. We will study the effect of symmetries by investigating how the positions of the leads affect the behavior of the Fano resonances.

In Sec. 2, we will discuss cavity states and the scattering states that are used to obtain the scattering S-matrix for a system with conical leads. In Sec. 3, quantum properties of the closed circular cavity with a cut will be reviewed. In Secs. 4 and 5, we will look at the behavior of conductance and Wigner delay times for waveguides with full-circle and cut-circle cavities. Leads will be attached to cavities symmetrically or asymmetrically. Finally, in Sec. 6, we will summarize our results. A short discussion of the numerical method used for our computations is given in the Appendix.

## II Scattering wavefunction and S-matrix

Matter waves inside the waveguide pictured in Fig. 1 are governed by the Schrödinger equation, which, in polar coordinates, takes the form

$$(^2+k^2)\mathrm{\Psi }(\rho ,\phi )=\left(\frac{^2}{\rho ^2}+\frac{1}{\rho }\frac{}{\rho }+\frac{1}{\rho ^2}\frac{^2}{\phi ^2}+\frac{2mE}{\mathrm{}^2}\right)\mathrm{\Psi }(\rho ,\phi )=0.$$
(1)
Here, $`k=\sqrt{2mE/\mathrm{}^2}`$ is the wavevector of the particle wave, $`m`$ is the particle mass, $`\mathrm{}`$ is Planck’s constant, and $`E`$ is the energy. The geometry of the system, which determines its scattering properties, is taken into account by the boundary condition that the scattering wavefunction $`\mathrm{\Psi }`$ vanish on the walls.

The solution $`\mathrm{\Psi }`$ is computed numerically by expanding it into suitable basis functions, and determining the expansion coefficients from the Schrödinger equation and its boundary conditions by the method described in the Appendix. Below, we will discuss the expressions for $`\mathrm{\Psi }`$ in the different regions of our system. In Sec. II A, we look at the wavefunction inside the cavity. The wavefunction inside the leads is considered in Sec. II B, where we will also define the scattering S-matrix of the waveguide.

### A Cavity states

In the interior region of the cavity, the solution of Eq. (1) can be written as

$$\mathrm{\Psi }=\underset{\gamma =n_\gamma +1}{\overset{n_\gamma }{}}c_\gamma f_\gamma .$$
(2)
Here, the $`f_\gamma `$’s are suitable expansion functions. The $`c_\gamma `$’s are the expansion coefficients, which must be determined from the boundary conditions, by the procedure described in the Appendix. The cutoff number $`n_\gamma `$ has to be chosen in a proper way to achieve optimal numerical accuracy within a reasonable computation time. We used $`n_\gamma =90`$ for all computations presented here.

For the functions $`f_\gamma `$, several choices are possible in principle. For our geometry, it proved to be numerically most efficient to use solutions of the Helmholtz equation \[$`(^2+k^2)f_\gamma =0`$, without boundary conditions\] which can be separated into polar coordinates, ($`\rho ,\phi `$). These functions take the form

$$f_\gamma (\rho ,\phi )=J_{|\gamma |}(k\rho )e^{i\gamma \phi },$$
(3)
where $`J_{|\gamma |}`$ denotes the Bessel function of order $`|\gamma |`$.

### B Lead states

For waveguide scattering problems, it is usual to express the properties of the system in terms of an S-matrix, which contains reflection and transmission amplitudes from incoming channels to outgoing channels in the leads. Before we can properly define the S-matrix for our system, we first have to specify what we mean by incoming and outgoing channels (also called “modes”) in this particular case.

Following the conventions from scattering problems with straight leads, we require that an incoming (outgoing) mode $`\chi ^{}`$ ($`\chi ^+`$) is a solution of

$$(^2+k^2)\chi ^\pm =0$$
(4)
which is separable into coordinates longitudinal and transversal to the lead walls. As a further condition, the transverse part has to be a standing wave vanishing on the walls and the longitudinal part describes a wave propagating towards (away from) the cavity.

For the conical shape of our leads, the longitudinal and transversal coordinates are just the polar coordinates $`\rho `$ and $`\phi `$, respectively. Thus, the separation ansatz with a standing wave in the transverse direction ($`\alpha =1,2,\mathrm{}`$) takes the form

$$\chi _\alpha ^{i\pm }=Nu_\alpha ^{i\pm }(\rho )\mathrm{sin}\left[\frac{\alpha \pi \left(\phi \theta _i+\frac{\mathrm{\Delta }\theta _i}{2}\right)}{\mathrm{\Delta }\theta _i}\right],$$
(5)
where $`\theta _i`$ locates the angular position of the center of the $`i`$th lead, $`R\rho <\mathrm{}`$, and $`\theta _i\mathrm{\Delta }\theta _i/2\phi <\theta _i+\mathrm{\Delta }\theta _i/2`$ ($`i=`$1, 2 for leads 1 and 2, respectively). $`N`$ is a (yet unknown) constant factor.

Inserting this into Eq. (4) leads to an equation for the radial part of the wavefunction,

$$\left[\rho ^2\frac{^2}{\rho ^2}+\rho \frac{}{\rho }\left(\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}\right)^2+k^2\rho ^2\right]u_\alpha ^{i\pm }(\rho )=0.$$
(6)
This is the Bessel differential equation and therefore $`u_\alpha ^{i\pm }(\rho )`$ can be expressed as a linear combination of Bessel functions, $`J_{\alpha \pi /\mathrm{\Delta }\theta _i}(k\rho )`$ and $`Y_{\alpha \pi /\mathrm{\Delta }\theta _i}(k\rho )`$. Since we required that the longitudinal part represent a propagating wave, we write the radial solution in terms of Hankel functions,

$$u_\alpha ^\pm (\rho )=H_{\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}}^\pm (k\rho )=J_{\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}}(k\rho )\pm iY_{\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}}(k\rho ),$$
(7)
which approach exponential functions at infinity:

$$H_{\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}}^\pm (k\rho )\stackrel{k\rho \mathrm{}}{}\sqrt{\frac{2}{\pi k\rho }}\mathrm{exp}\left\{\pm i\left[k\rho \frac{\pi }{4}\left(2\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}+1\right)\right]\right\}+𝒪\left[\frac{1}{(k\rho )^{3/2}}\right].$$
(8)
Here, we used the notation $`H_\gamma ^+`$ ($`H_\gamma ^{}`$) for the Hankel functions of first (second) kind, of order $`\gamma `$.

The normalization constant, $`|N|`$, is calculated from the condition that each channel carry unit current, i.e., we require that

$$\left|_{\theta _i\mathrm{\Delta }\theta _i/2}^{\theta _i+\mathrm{\Delta }\theta _i/2}\rho 𝑑\phi j_\rho ^{i\pm }\right|=1.$$
(9)
Here, the longitudinal component of the probability current density, $`𝐣^{i\pm }=j_\rho ^{i\pm }𝐞_\rho +j_\phi ^{i\pm }𝐞_\phi `$, is given by

$$j_\rho ^{i\pm }=\frac{i\mathrm{}}{2m}\left[\left(\frac{\chi ^{i\pm }}{\rho }\right)^{}\chi ^{i\pm }\chi ^{i\pm ^{}}\left(\frac{\chi ^{i\pm }}{\rho }\right)\right].$$
(10)
This yields

$$|N|=\sqrt{\frac{\pi m}{\mathrm{}\mathrm{\Delta }\theta _i}}.$$
(11)
With these results, the expression for the incoming/outgoing channel $`\alpha =1,2,\mathrm{}`$, corresponding to the $`\alpha `$th standing wave in transverse direction, becomes

$$\chi _\alpha ^\pm =\frac{N}{|N|}\sqrt{\frac{m\pi }{\mathrm{}\mathrm{\Delta }\theta _i}}\mathrm{sin}\left[\frac{\alpha \pi \left(\phi \theta _i+\frac{\mathrm{\Delta }\theta _i}{2}\right)}{\mathrm{\Delta }\theta _i}\right]H_{\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}}^\pm (k\rho ).$$
(12)
To fix the overall phase factor, $`N/|N|`$, in Eq. (12), we require that the solution be real at the interface between the leads and the cavity. Then, if a mode is completely reflected back at this position, the corresponding reflection coefficient will be 1. This can be achieved by multiplying the modes with a complex phase factor. At $`\rho =R`$, the Hankel functions (the only complex factor in $`\chi _\alpha ^{i\pm }`$) can be written as

$$H_{\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}}^\pm (kR)=\left|H_{\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}}^\pm (kR)\right|\mathrm{exp}\left\{\pm i\mathrm{arctan}\left[\frac{Y_{\alpha \pi /\mathrm{\Delta }\theta _i}(kR)}{J_{\alpha \pi /\mathrm{\Delta }\theta _i}(kR)}\right]\right\}.$$
(13)
Because we want the probability amplitude at this position to be real, we will use the complex phase factor

$$\frac{N}{|N|}=\mathrm{exp}\left\{i\mathrm{arctan}\left[\frac{Y_{\alpha \pi /\mathrm{\Delta }\theta _i}(kR)}{J_{\alpha \pi /\mathrm{\Delta }\theta _i}(kR)}\right]\right\}.$$
(14)
Then, the final form of the incoming and outgoing propagating modes in channel $`\alpha `$ becomes

$$\chi _\alpha ^{i\pm }(\rho ,\phi )=\sqrt{\frac{m\pi }{\mathrm{}\mathrm{\Delta }\theta _i}}\mathrm{exp}\left\{i\mathrm{arctan}\left[\frac{Y_{\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}}(kR)}{J_{\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}}(kR)}\right]\right\}\mathrm{sin}\left[\frac{\alpha \pi \left(\phi \theta _i+\frac{\mathrm{\Delta }\theta _i}{2}\right)}{\mathrm{\Delta }\theta _i}\right]H_{\frac{\alpha \pi }{\mathrm{\Delta }\theta _i}}^\pm (k\rho ).$$
(15)
As we now have derived the expression for incoming and outgoing channels, we are finally in the position to define the S-matrix of our system. The matter wave $`\mathrm{\Psi }`$ inside the leads is expanded in terms of the incoming and outgoing modes, $`\chi _\alpha ^i(\rho ,\phi )`$ and $`\chi _\alpha ^{i+}(\rho ,\phi )`$, respectively. Then the expansion coefficients describe the probability amplitudes of reflection or transmission of the matter wave from an incoming channel to the outgoing channels. More precisely, if the electron is incident in lead 1 and channel $`\beta `$, the wavefunction in the leads takes the form

$$\mathrm{\Psi }_{\beta ,1}=\chi _\beta ^1+\underset{\alpha }{}r_{\alpha \beta }\chi _\alpha ^{1+}\mathrm{and}\mathrm{\Psi }_{\beta ,2}=\underset{\alpha }{}t_{\alpha \beta }\chi _\alpha ^{2+}.$$
(16)
Here, $`\mathrm{\Psi }_{\beta ,1}`$ is the state in lead 1, $`\mathrm{\Psi }_{\beta ,2}`$ is the state in lead 2, and $`r_{\alpha \beta }(t_{\alpha \beta })`$ is the probability amplitude that the electron enters in channel $`\beta `$ and is reflected (transmitted) into channel $`\alpha `$. Likewise, if the electron is incident in lead 2 and channel $`\beta `$, the wavefunction in the leads takes the form

$$\mathrm{\Psi }_{\beta ,2}=\chi _\beta ^2+\underset{\alpha }{}r_{\alpha \beta }^{^{}}\chi _\alpha ^{2+}\mathrm{and}\mathrm{\Psi }_{\beta ,1}=\underset{\alpha }{}t_{\alpha \beta }^{^{}}\chi _\alpha ^{1+}.$$
(17)
Here $`\mathrm{\Psi }_{\beta ,2}`$ is the state in lead 2, $`\mathrm{\Psi }_{\beta ,1}`$ is the state in lead 1, and $`r_{\alpha \beta }^{^{}}`$ and $`t_{\alpha \beta }^{^{}}`$ are the reflection and transmission probability amplitudes, respectively. Note that the solution $`\mathrm{\Psi }`$ of Eq. (1) is not unique, but we obtain different solutions $`\mathrm{\Psi }_\beta `$ depending on which channel $`\beta `$ is assumed to carry the incoming part of the matter wave.

The S-matrix connects the amplitudes of the incoming propagating electrons to the amplitudes of the outgoing propagating electrons and it is therefore constructed out of all the reflection and transmission amplitudes. Thus, the S-matrix can be written

$$𝐒=\left(\begin{array}{cc}𝐑& 𝐓^{}\\ 𝐓& 𝐑^{}\end{array}\right),$$
(18)
where $`𝐑`$ and $`𝐑^{}`$ are square matrices composed of reflection amplitudes, $`r_{\alpha \beta }`$ and $`r_{\alpha \beta }^{}`$, and $`𝐓`$ and $`𝐓^{}`$ are square matrices composed of transmission amplitudes, $`t_{\alpha \beta }`$ and $`t_{\alpha \beta }^{}`$. The reflection and transmission amplitudes must be computed numerically (see Appendix).

## III Closed cut-circle billiard

One aspect of the scattering properties of the waveguide that we wish to explore is how close a correspondence there is between the Fano scattering resonances and the eigenstates of a closed version of the waveguide cavity. In this section, we therefore describe some features of these eigenstates as the cut size on the cut-circle billiard is changed.

The Schrödinger equation \[Eq. (1)\] for a quantum particle in a closed full-circle billiard is separable, with one part describing the radial and the other part the azimuthal motion. The (discrete) energy spectrum can be labeled with two good quantum numbers. One of them, $`n=1,2,\mathrm{}`$, is associated with radial motion, and the other, $`l=0,\pm 1,\pm 2,\mathrm{}`$, with the angular motion. The eigenfunctions are

$$\mathrm{\Psi }_{l,n}=J_{|l|}\left(k_{|l|,n}\rho \right)e^{il\phi },$$
(19)
where $`J_{|l|}`$ is the Bessel function of order $`|l|`$ and $`(k_{|l|,n}R)`$ is its $`n`$th zero (i.e., $`J_{|l|}`$ vanishes at $`\rho =R`$). The corresponding eigenenergies $`E_{|l|,n}`$ are then calculated via $`E=\mathrm{}^2k^2/(2m)`$. The energy eigenstates of the circle billiard with angular momentum $`|l|1`$, corresponding to clockwise or counter clockwise rotation of the matter wave, produce two-fold degenerate eigenenergies. Only the eigenenergies with $`l=0`$ are non-degenerate.

If we apply a cut to the closed system, the classical dynamics of the particle exhibits hard chaos , as long as $`\omega 180^{}`$. Looking at the quantum properties, we observe a splitting of the degenerate eigenenergies as the cut is inserted. The absence of degenerate eigenvalues (level repulsion) is a well-known quantum feature of classically chaotic systems (see, e.g., Ref. ).

Examining the quantum properties of the cut-circle billiard in more detail, we show how five of the energy levels of the cut billiard behave as the cut size is varied from $`\omega =50^{}`$ to $`\omega =70^{}`$ in Fig. 2. The states shown are continuations of the circle cavity degenerate pairs $`(l=3,n=2)`$, $`(l=6,n=1)`$, and one of the pair, $`(l=1,n=3)`$.

Here and in later numerical computations, we are using dimensionless variables. The dimensionless energy $`ϵ`$ is given by

$$ϵ=\frac{mR^2}{\mathrm{}^2}E=\frac{(kR)^2}{2}=\frac{\kappa ^2}{2},$$
(20)
where $`\kappa =kR`$ is a dimensionless wavevector. To get an idea of how this translates to the physical dimensions of our system, consider an electron density of $`n\mathrm{\hspace{0.25em}4}\times 10^{11}\mathrm{cm}^2`$, which is a typical value measured for a two dimensional electron gas in a GaAs/Al0.3Ga0.7As heterostructure . The corresponding Fermi wavevector is $`k_f=\sqrt{2\pi n}\mathrm{\hspace{0.25em}1.6}\times 10^8\mathrm{m}^1`$, the Fermi energy $`E_f=\pi \mathrm{}^2n/m0.014\mathrm{eV}`$ (for an effective electron mass in GaAs of $`m=0.067m_e`$). This yields for $`ϵ=100`$, which is the maximum energy we used in the computations, $`R\mathrm{\hspace{0.25em}90}\mathrm{nm}`$.

In Fig. 2, we observe a (distant) avoided crossing, which occurs between states 2b and 3a in the interval $`\omega =53^{}`$ to $`\omega =63^{}`$. We also see an actual crossing in this plot. This occurs between the states labeled 2a and 2b at about $`\omega =65^{}`$.

However, as discussed in Ref. , separate eigenvalues of a generic Hamiltonian can in general only be brought to coincide if at least three parameters are varied. Thus, at first sight, we would expect another level repulsion instead of the actual crossing, since we only vary $`\omega `$. Yet we ignored that the Hamiltonian of the cut-circle is not completely “generic”. Even if the cut is inserted, the system retains parity symmetry with respect to an axis perpendicular to the cut, i.e.,

$$P_\zeta \mathrm{\Psi }(\rho ,\phi )\mathrm{\Psi }(\rho ,\zeta \phi )=\mathrm{\Psi }(\rho ,\phi ),$$
(21)
where $`P_\zeta `$ is the parity operator for a symmetry axis through the center of the cavity, at angle $`\zeta `$ (i.e., the angle at which the center of the cut is placed). Because of parity symmetry, we have two different symmetry classes, namely the classes formed by even and odd states. Since parity is a discrete symmetry, states of one class cannot interact with states of the other, thus both classes are completely independent of each other. Therefore, crossings between states of different parity may occur. Only within each parity class, i.e., among states with the same parity, crossings remain impossible.

In Fig. 3, we show energy eigenstates corresponding to the eigenvalues 2a, 2b, and 3a in Fig. 2, for three different cut sizes taken before the avoided crossing, between the avoided and actual crossings, and after the actual crossing. It can be seen that 2a has even parity, and 2b and 3a both have odd parities. This explains why 2a and 2b can undergo the crossing, whereas a crossing between 2b and 3a is avoided.

We also see that the state 2a, which does not participate in the avoided crossing and undergoes the crossing with 2b, remains unchanged. However, states 2b and 3a, which do undergo an avoided crossing, become mixed and loose their original identities. The mutual parity of states 2b and 3a allows them to couple, whereas this is impossible for states 2a and 2b with differing parities.

## IV Waveguide with a circle cavity

If we attach leads to the circle billiard, we have a matter waveguide and the eigenstates of the billiard become unstable and may decay. These unstable states strongly affect the scattering dynamics of the waveguide. They appear as resonances in the transmission probabilities through the waveguide,

$$T_\beta =\underset{\alpha }{}|t_{\alpha \beta }|^2,$$
(22)
for incoming mode $`\beta `$ in lead 1, where $`\alpha `$ runs over the outgoing modes in lead 2. We can examine the individual $`T_\beta `$’s or consider their sum over all possible incoming modes in lead 1, which is proportional to the Landauer-Büttiker conductance ,

$$G=\frac{2e^2}{h}\underset{\beta }{}T_\beta .$$
(23)
As shown in Ref. , there is an energy range where, for conical leads, mode $`\beta `$ has nonzero transmission, whereas for straight leads with the same opening angles it would still be evanescent. This occurs for energies too small for the $`\beta `$th transverse standing wave to fit into the lead opening at the junction to the cavity, such that the wave has to “tunnel” into the cavity. The condition for this tunneling (see Ref. ) is $`k<\beta \pi /d_i`$, if $`d_i=2R\mathrm{sin}(\mathrm{\Delta }\theta _i/2)`$ is the width of the lead opening (for opening angle $`\mathrm{\Delta }\theta _i`$). Thus, tunneling occurs for energies such that

$$ϵ<\stackrel{~}{ϵ}_\beta =\frac{\beta ^2\pi ^2R^2}{2d_i^2}=\frac{\beta ^2\pi ^2}{8\mathrm{sin}^2(\mathrm{\Delta }\theta _i/2)},$$
(24)
where $`ϵ`$ is the dimensionless energy defined in Eq. (20).

In addition to the transmission probabilities and Landauer-Büttiker conductance, it is useful to look at the Wigner delay times,

$$\tau _n^W=\mathrm{}\frac{\theta _n}{E},$$
(25)
where $`\theta _n`$ are the phases of the S-matrix eigenvalues $`s_n=e^{i\theta _n}`$. The Wigner delay times are a measure of the time delay of the electron due to the presence of the cavity. In the following plots, we will use dimensionless delay times

$$\tau _n=\frac{\theta _n}{ϵ}=\frac{\mathrm{}\tau _n^W}{m_e^{}R^2}.$$
(26)
(For $`R100\mathrm{nm}`$, $`\tau _n20`$, and $`m_e^{}=0.067m_e`$, the Wigner delay time is of the order of 0.1 ns.) From Ref.  it is clear that peaks in the Wigner delay time correspond to resonances and poles of the scattering matrix in the complex energy plane. Therefore we can use peaks in the delay time spectra to indicate the resonances of our system.

### A Symmetric placement of leads

In Fig. 4, we show numerical data obtained for an electron waveguide with a circle cavity with symmetrically placed leads, one at $`\theta _1=0^{}`$ and the other at $`\theta _2=180^{}`$. Both leads subtend the same angle $`\mathrm{\Delta }\theta _1=\mathrm{\Delta }\theta _2=20^{}`$. The dimensionless energy $`ϵ`$ is varied in steps of 0.5.

The transmission probability in Fig. 4(a) shows the quantities $`T_\beta `$ for incoming modes $`\beta =1,2`$ in lead 1. The first mode ($`\beta =1`$) starts transmitting (i.e., becomes significantly larger than zero) at about $`ϵ20`$. The second mode starts transmitting too closely below 100 to be visible in the plot. The transmission for higher modes is zero everywhere in the energy interval below 100, and therefore these modes were not included in the computations or the plot.

The threshold energies $`\stackrel{~}{ϵ}_\beta `$, for which the modes would start propagating in the case of straight leads with $`\mathrm{\Delta }\theta =20^{}`$, are $`\stackrel{~}{ϵ}_1=40.91`$, and $`\stackrel{~}{ϵ}_2=163.66`$ for the first two channels. This means that for the conical leads, we observe tunneling of the first mode in the energy range between $`ϵ20`$ and 40.

Examining the correspondence between eigenenergies of the closed circle (indicated by filled circles for degenerate levels and open circles for non-degenerate levels) and resonances of the open system, one finds that there is obviously a strong association. In energy intervals without eigenenergies, the transmission changes smoothly, whereas there are sharp changes close to eigenenergies. In the tunneling regime, we see peaks close to the eigenenergies. The peaks tend to be slightly shifted to the left. Above the tunneling threshold, we also observe valleys, e.g., at $`ϵ=67.5,74.7`$. In this regime, most peaks and valleys in the vicinity of an eigenenergy are too close together to tell which is associated with the eigenenergy. Since the lead openings themselves alter the geometry of the system, we cannot expect the resonances to be found exactly at the same positions as the eigenenergies.

In Fig. 4(b), we plot the dimensionless delay times $`\tau _n`$. The energy derivatives at each point are approximated by carrying out two computations for each point, with energy spacing $`\mathrm{\Delta }ϵ=10^5`$. Since the transmission probabilities $`T_\beta `$ vanish for $`\beta 3`$, we only have to consider the two channels $`\beta =1`$ and $`\beta =2`$. Thus, the S-matrix is a 4x4 matrix and we get four eigenphases and four delay time curves.

Here, unlike for the transmission probabilities, we only need to look for peaks in the Wigner delay times at the resonance energies, in the tunneling as well as in the normal regime. The first two time delays display their first peak around $`ϵ15`$, the others at $`ϵ75`$. At these energies, the corresponding transmission probabilities are still too close to zero to produce a visible peak in Fig. 4(a). This shows that the delay times are a more sensitive indicator of resonances than the $`T_\beta `$’s.

We see that the degenerate eigenenergies of the circle billiard do not produce double resonances when leads are attached at the special angles $`\theta _1=0^{}`$ and $`\theta _2=180^{}`$. Every eigenenergy of the circle billiard is associated only with a single delay time peak. Although placing leads to the circle cavity destroys its radial symmetry, the breaking of the degeneracy is not observed in this case.

### B Asymmetric placement of leads

In Fig. 5, we show the same waveguide, but with asymmetric placement of leads instead of the symmetric placement we had in Fig. 4. Lead 1 is placed at $`\theta _1=0^{}`$ and lead 2 at $`\theta _2=125^{}`$.

Now we not only see again the close correspondence between delay time peaks and eigenenergies of the circle, but we also observe double peaks at most of the degenerate billiard eigenenergies, most evidently at $`ϵ=`$24.6, 28.8, 35.4, 38.5, and 49.4. At higher energies, the observation of double peaks is more difficult because the second mode starts transmitting. (Some of its first peaks are too small to be seen in the plot, but can be found by looking at the data file.) Many of the double peaks show a small energy spacing: One peak is slightly shifted to smaller energies, the other one to higher energies. Here, the degeneracy is obviously broken by the addition of the leads. For the non-degenerate eigenenergies on the other hand, no splitting can be observed. Especially at $`ϵ=37.4`$, this becomes very obvious.

The difference between this case and the previous case with symmetric lead placement ($`\theta _2=180^{}`$) is that, for symmetric placement, we had two discrete symmetries, namely invariance under the parity transformations $`P_0^{}`$ and $`P_{90^{}}`$ \[where the parity operators are defined as in Eq. (21)\]. When the leads are placed asymmetrically, the symmetry under $`P_0^{}`$, with respect to the axis bisecting both leads, is destroyed, and we are only left with one parity symmetry under $`P_{62.5^{}}`$. The breaking of one parity symmetry seems to be the reason for the splitting of resonances seen in Fig. 5(b).

### C Numerical accuracy

To get an estimate of the numerical error, the sum of the absolute values squared of all reflection and transmission amplitudes was calculated. Due to the unitarity of the S-matrix, this value has to be the dimension of the S-matrix, which is 4 here. The relative error,

$$1\frac{1}{4}\underset{\alpha ,\beta =1}{\overset{2}{}}(|r_{\alpha \beta }|^2+|t_{\alpha \beta }|^2+|r_{\alpha \beta }^{^{}}|^2+|t_{\alpha \beta }^{^{}}|^2),$$
(27)
was usually smaller than 0.2% for the individual data points. The errors of the absolute values of the S-matrix eigenvalues (which should be 1) were of the same order of magnitude. For the delay times, there is no easy way to estimate the error.

The most crucial factor determining numerical accuracy was found to be the number $`2n_\gamma `$ of parameters (expansion coefficients $`r_{\alpha \beta },t_{\alpha \beta },b_{\nu \beta }`$) used to approximate the total wavefunction $`\mathrm{\Psi }`$ of our system \[see Eqs. (16), (A10), and (A19)\]. An estimate for a suitable choice of this number can be found by comparing the Fermi wavelength, $`\lambda _f=2\pi /k_f`$, to a typical length scale of the system, e.g., the length of the boundary $`2\pi R`$ (for zero or small cut size). If we assume that the total wavefunction of our system varies rather smoothly on a length scale of $`\lambda _f`$, it seems plausible that a “small” number of parameters, say about 10, should be sufficient to approximate $`\mathrm{\Psi }`$ on an interval of length $`\lambda _f`$. Therefore, we expect that such a “small” multiple of $`\kappa =k_fR=\sqrt{2ϵ}`$=“the number of $`\lambda _f`$-intervals needed to cover the cavity boundary” should give a reasonable order of magnitude estimate for $`2n_\gamma `$.

In Fig. 6, we compare the transmission spectra \[Fig. 6(a)\] and the spectra of the time delay sums \[Fig. 6(b)\] for different values of $`n_\gamma `$. As deviations are expected to be most pronounced for large cut size, the calculations shown were done exemplarily for the symmetric cut-circle waveguide to be examined in Sec. V A, at $`\omega =80^{}`$.

We find noticeable deviations for $`n_\gamma =30`$ and, to a small extent, also for $`n_\gamma =60`$. The values for $`n_\gamma =90`$ practically coincide with those for $`n_\gamma =120`$. Therefore, in view of computation time, we used $`n_\gamma =90`$ in all computations presented in this paper.

## V Waveguide with a cut-circle cavity

We now consider a waveguide with a cut-circle cavity. Fig. 7 shows the sum of the partial delay times, $`\tau =_{n=1}^4\tau _n`$, and the Landauer-Büttiker conductance $`G`$. For this computation, leads and cut are placed asymmetrically. The leads are centered at $`\theta _1=0^{}`$ and $`\theta _2=166^{}`$, and the cut at $`\zeta =225^{}`$. The cut size $`\omega `$ varies from $`0^{}`$ (circle) to $`80^{}`$, with step size $`5^{}`$. The dimensionless energy $`ϵ`$ is varied in steps of 0.25.

In the delay time spectrum (upper plot), we observe chains of separated double peaks, corresponding to the split-up branches of degenerate circle eigenenergies. If we compare the positions of maxima in the Wigner delay time to the positions of eigenenergies of the closed system, we find a good agreement between (open system) resonances and (closed cavity) eigenenergies.

The conductance (lower plot) shows similar features. However, the first visible chains of resonances occur at significantly higher energies than in the delay time plot. For energies between 20 and 40, both plots show basically the same structures. Above the first tunneling threshold at $`\stackrel{~}{ϵ}_1=40.91`$, the conductance reflects the onset of “normal” (compared to tunneling) transmission of the first mode. In this energy regime, we cannot tell any more whether resonances are indicated by peaks or dips in the conductance.

We will now compare the properties of the cut-circle cavity waveguide for symmetric and asymmetric lead placements. Note, however, that although this may seem to be analogous to the analysis from the previous section (for the full-circle cavity), the insertion of the cut already leaves us with only one parity symmetry (under $`P_\zeta `$). As a result, the “symmetric” cut-circle waveguide is analogous to the asymmetric full-circle waveguide rather than to the symmetric full-circle waveguide. Finally, in the “asymmetric” cut-circle waveguide, there are no geometric symmetries present at all. Since this could not be achieved by merely placing two leads to a circle, this case has no parallel in Sec. IV.

### A Symmetric placement of leads and cut

Let us examine the Wigner delay times for the energy range where we found an actual crossing and an avoided crossing in Fig. 2 for the closed cut-circle billiard. We first consider a system in which the leads are placed at $`\theta _1=0^{}`$ and $`\theta _2=120^{}`$, symmetrically with respect to the cut at $`\zeta =240^{}`$. The lead openings are $`\mathrm{\Delta }\theta _1=\mathrm{\Delta }\theta _2=16^{}`$. The plot of the sum $`\tau `$ of the Wigner delay times is shown in Fig. 8.

We observe five delay time peaks, which are Fano resonances resulting from the closed cut billiard eigenstates 1a, 1b, 2a, 2b, and 3a shown in Fig. 2. Starting from the left, we label these delay time peaks in the same manner. For cut size $`\omega =20^{}`$, we designate the delay time peaks, going from left to right, as 1b, 2a, 2b, 3a, and 3b. 1a is the resonance chain in the upper left corner. The first two delay time peaks, 1a and 1b, are Fano resonances resulting from the circle eigenstate with $`(l=3,n=2)`$, the following two from $`(l=6,n=1)`$, and the last two, 3a and 3b, from $`(l=1,n=3)`$.

Again, we observe an avoided crossing between 2b and 3a, and the two delay time peaks 2a and 2b actually cross. Like for the closed cavity, the crossing is possible due to the parity symmetry of the system, which has been preserved by adding the leads symmetrically with respect to the cut. Since 2a and 2b have different parities, these eigenvalue chains are allowed to cross.

The symmetric placement of leads has not broken any additional symmetries. Only the energies and cut sizes at which the avoided crossing and actual crossing take place are slightly shifted from the case of the closed cut billiard.

### B Asymmetric placement of leads and cut

In Fig. 9 we show the same Fano resonances as in Fig. 8, but for asymmetric placement of leads. Now we place the leads at $`\theta _1=0^{}`$ and $`\theta _2=166^{}`$ and the cut at $`\zeta =225^{}`$.

The crossing of resonances 2a and 2b, which occurred for the symmetric leads and for the closed cut-circle billiard, has now become an avoided crossing. (In fact we now have a sequence of two closely spaced avoided crossings, similar to the case studied in Ref. .) Here, the parity symmetry of the closed system is destroyed by placing the leads asymmetrically. Therefore, in the asymmetric open system, we no longer have two separate classes of states (namely states of even and odd parity) and the crossing, which was possible in the closed system and in the symmetric open system, is avoided here.

### C Cavity wavefunctions

We now look at the distribution of electron probability in the cavity at the resonance energies. This is shown in Fig. 10 for the asymmetric system. A three-way avoided crossing has mixed the resonance states, and has destroyed the identity they had for smaller cut size. The upper row displays the four states before the avoided crossing, which still reflect the quantum numbers of the associated circle states very well. For the middle row, the character of states 2a, 2b, and 3a has been changed. Only the resonance state 1b, which never comes really close to the crossing region, has preserved its characteristic structure (two rings with six peaks each) quite well. In the bottom row, after the second avoided crossing, 2a and 2b have mixed even further.

## VI Conclusion

In this paper, we investigated a waveguide with a cut-circle cavity by calculating conductances and Wigner delay times. We attached conic leads to the cavity in two ways, symmetrically and asymmetrically. We observed two kinds of avoided crossings when the leads were placed asymmetrically, breaking the parity symmetry of the closed cavity. One avoided crossing was due to the chaos of the closed cavity, i.e., due to level repulsions between states in the same symmetry class, which are known to occur for fully chaotic systems. The other avoided crossing was due to the breaking of the discrete symmetry (the parity in our case) by placing leads asymmetrically (see Fig. 9).

There is now great interest in exploring the effects of underlying chaos on the dynamics of open quantum systems. However, open quantum systems must make contact with the outside world. This contact itself may induce some of the signatures of chaos. This fact has also been studied extensively by Jung and Seligman  for classical scattering systems. The avoided crossing between 2a and 2b, shown in Fig. 9, provides another example of this effect.

## A Numerical method

### 1 Basic concept

In order to obtain an explicit expression for the S-matrix, we need to find the stationary solutions $`\mathrm{\Psi }`$ of Eq. (1). This can only be done numerically for our system. For that purpose, we use the boundary integral method described in Ref. . Below, we will give an outline of this method adapted to our particular problem.

The boundary integral method is based on the use of Green’s Theorem, which, for two functions $`f`$ and $`g`$, states that

$$_\mathrm{\Gamma }𝑑l\left(f\frac{g}{n}\frac{f}{n}g\right)=_A𝑑A^{}(f^2gg^2f),$$
(A1)
if $`\mathrm{\Gamma }`$ is a closed contour confining an area $`A`$, $`l`$ is the longitudinal coordinate along $`\mathrm{\Gamma }`$, and $`n`$ is the normal coordinate. For the boundary integral method, we choose for the functions $`f`$ and $`g`$ the wavefunction $`\mathrm{\Psi }`$ and a weight function $`\mathrm{\Phi }`$ that satisfies the Helmholtz equation

$$^2\mathrm{\Phi }=k^2\mathrm{\Phi }.$$
(A2)
Then, since $`\mathrm{\Psi }`$ is a solution of the Schrödinger equation \[Eq. (1)\], i.e., $`\mathrm{\Psi }`$ also satisfies $`^2\mathrm{\Psi }=k^2\mathrm{\Psi }`$, the right hand side of Eq. (A1) is identically zero. Therefore we get

$$_\mathrm{\Gamma }𝑑l\left(\mathrm{\Psi }\frac{\mathrm{\Phi }}{n}\frac{\mathrm{\Psi }}{n}\mathrm{\Phi }\right)=0$$
(A3)
as the basic equation of this numerical approach. (Note that the difference between $`\mathrm{\Psi }`$ and $`\mathrm{\Phi }`$ is that there are no boundary conditions for $`\mathrm{\Phi }`$, whereas $`\mathrm{\Psi }`$ has to satisfy the boundary conditions imposed by the geometry of the waveguide.)

For our problem, we choose $`\mathrm{\Gamma }`$ to follow the walls of the cavity. For the places where the leads come in, $`\mathrm{\Gamma }`$ is taken to be the continuation of the circular wall segments in order to construct a closed contour. Then, Eq. (A3) allows us to obtain the S-matrix and the probability distribution inside the cavity for any incoming particle energy. This is done by inserting (approximate) expansions of $`\mathrm{\Psi }`$ with, say, $`n_{\mathrm{cutoff}}`$ yet unknown expansion coefficients. As we need $`n_{\mathrm{cutoff}}`$ equations to solve for these expansion coefficients, we apply Eq. (A3) with sufficiently many different weight functions $`\mathrm{\Phi }_n`$ ($`n=1,\mathrm{},n_{\mathrm{cutoff}}`$).

In principle, there are several choices possible for the $`\mathrm{\Phi }_n`$’s, as long as they satisfy Eq. (A2). For numerical convenience, we used the same functions as for the expansion of $`\mathrm{\Psi }`$ inside the cavity, namely $`\mathrm{\Phi }_\gamma =f_\gamma `$ with the $`f_\gamma `$’s given by Eq. (3), where $`\gamma =n_\gamma +1,\mathrm{},n_\gamma `$ (i.e., $`n_{\mathrm{cutoff}}=2n_\gamma `$). For $`n_\gamma `$, we also used the same value as in Sec. II A, $`n_\gamma =90`$, which allows us to determine $`n_{\mathrm{cutoff}}=180`$ expansion coefficients.

### 2 S-matrix

In order to compute the S-matrix, we split up $`\mathrm{\Gamma }`$ into parts $`C_1,C_2`$ across the lead openings, $`P_1,P_2,P_3`$ along the circular walls, and $`P_0`$ along the cut (see Fig. 11).

We first regard a part $`I_C`$ of the integral in Eq. (A3), for which we only integrate over the lead openings $`C_1`$ and $`C_2`$ instead of $`\mathrm{\Gamma }`$. Inserting for $`\mathrm{\Phi }`$ the $`f_\gamma `$’s and for $`\mathrm{\Psi }`$ its expansion into lead channels from Eq. (16) (for an electron incident in channel $`\beta `$ in lead 1), we obtain

$`I_C`$ $`=`$ $`{\displaystyle _{C_1}}𝑑l\left(\chi _\beta ^{}{\displaystyle \frac{f_\gamma }{n}}{\displaystyle \frac{\chi _\beta ^{}}{n}}f_\gamma \right)+{\displaystyle \underset{\alpha }{}}r_{\alpha \beta }{\displaystyle _{C_1}}𝑑l\left(\chi _\alpha ^+{\displaystyle \frac{f_\gamma }{n}}{\displaystyle \frac{\chi _\alpha ^+}{n}}f_\gamma \right)`$ (A5)
$`+{\displaystyle \underset{\alpha }{}}t_{\alpha \beta }{\displaystyle _{C_2}}𝑑l\left(\chi _\alpha ^+{\displaystyle \frac{f_\gamma }{n}}{\displaystyle \frac{\chi _\alpha ^+}{n}}f_\gamma \right)`$
$`=`$ $`(z_0^1)_{\gamma \beta }(z_L^1)_{\gamma \beta }+{\displaystyle \underset{\alpha }{}}r_{\alpha \beta }\left[(z_0^{1+})_{\gamma \alpha }(z_L^{1+})_{\gamma \alpha }\right]`$ (A7)
$`+{\displaystyle \underset{\alpha }{}}t_{\alpha \beta }\left[(z_0^{2+})_{\gamma \alpha }(z_L^{2+})_{\gamma \alpha }\right],`$
where we introduced the abbreviations

$`(z_L^{i\pm })_{\gamma \alpha }`$ $``$ $`{\displaystyle _{C_i}}𝑑l{\displaystyle \frac{\chi _\alpha ^{i\pm }}{n}}f_\gamma ,`$ (A8)
$`(z_0^{i\pm })_{\gamma \alpha }`$ $``$ $`{\displaystyle _{C_i}}𝑑l\chi _\alpha ^{i\pm }{\displaystyle \frac{f_\gamma }{n}},`$ (A9)
with $`i=`$1, 2 denoting lead 1 and 2. (For an electron incident in lead 2, the expressions look similar.)

Along the walls of the cavity, the electron wavefunction is zero, but it may have a finite slope perpendicular to the wall. We can expand this slope in terms of a complete orthonormal basis along the wall. Therefore, we can write

$`\mathrm{\Psi }_{\beta ,P_i}`$ $`=`$ $`0,`$ (A10)
$`\left({\displaystyle \frac{\mathrm{\Psi }_\beta }{n}}\right)_{P_i}`$ $`=`$ $`{\displaystyle \underset{\nu }{}}b_{\nu \beta }^{(i)}\xi _\nu ^{(i)}(l),`$ (A11)
with some set of basis functions $`\xi _\nu ^{(i)}`$ in variable $`l`$. The index $`i=0`$ is used to denote the cut, and $`i=1,2,3`$ denote the three different circular wall segments. We are free to choose any complete basis, and for the circular wall segments we choose a Fourier basis,

$$\xi _\nu ^{(i)}(l)=\sqrt{\frac{2}{\mathrm{\Delta }_i}}\mathrm{cos}\left[\frac{\pi \nu \left(l\delta _i+\frac{\mathrm{\Delta }_i}{2}\right)}{\mathrm{\Delta }_i}\right],$$
(A12)
where $`\delta _i`$ is the angular position of the center of wall segment $`i=1,2,3`$, $`\mathrm{\Delta }_i`$ is its opening angle, $`\nu =0,1,\mathrm{}`$, and $`\delta _i\frac{\mathrm{\Delta }_i}{2}l<\delta _i+\frac{\mathrm{\Delta }_i}{2}`$. For the cut, we could use the Fourier basis, too. However, as this results in very long computation times, a basis of small triangles proved to be more efficient. Thus, along the straight wall we use

$$\xi _\nu ^{(0)}(l)=\frac{2(n_{\nu _0}+1)}{l_0}\{\begin{array}{ccc}l\frac{\nu }{n_{\nu _0}+1}l_0\hfill & \mathrm{for}& \hfill \frac{\nu }{n_{\nu _0}+1}l_0l<\frac{\nu +1/2}{n_{\nu _0}+1}l_0\\ \frac{\nu +1}{n_{\nu _0}+1}l_0l\hfill & \mathrm{for}& \hfill \frac{\nu +1/2}{n_{\nu _0}+1}l_0l<\frac{\nu +1}{n_{\nu _0}+1}l_0\\ 0\hfill & \mathrm{else}& \end{array},$$
(A13)
for $`0l<l_0`$, where $`l_0=2R\mathrm{sin}(\omega /2)`$ is the length of the cut, $`\omega `$ is its opening angle, and $`\nu =0,1,\mathrm{},n_{\nu _0}`$. Now we can write down the part $`I_P`$ of the integral in Eq. (A3), for which we only integrate over the wall segments $`P_0,\mathrm{},P_3`$ instead of $`\mathrm{\Gamma }`$. Inserting the expansions from Eq. (A10), we get

$$I_P=\underset{i=0}{\overset{3}{}}\underset{\nu }{}b_{\nu \beta }^{(i)}_{P_i}𝑑l\xi _\nu ^{(i)}f_\gamma =\underset{i=0}{\overset{3}{}}\underset{\nu }{}b_{\nu \beta }^{(i)}(x^i)_{\gamma \nu },$$
(A14)
with the abbreviation

$$(x^i)_{\gamma \nu }_{P_i}𝑑l\xi _\nu ^{(i)}(l)f_\gamma .$$
(A15)
Now we can combine these results to obtain a matrix equation for the unknown expansion coefficients, which can easily be solved numerically. With the integrals $`I_C`$ and $`I_P`$, Eq. (A3) becomes $`I_C+I_P=0`$. Inserting Eqs. (A5) and (A14) yields (after a slight reordering of terms)

$`{\displaystyle \underset{\alpha }{}}\left[(z_0^{1+})_{\gamma \alpha }(z_L^{1+})_{\gamma \alpha }\right]r_{\alpha \beta }+{\displaystyle \underset{\alpha }{}}\left[(z_0^{2+})_{\gamma \alpha }(z_L^{2+})_{\gamma \alpha }\right]t_{\alpha \beta }`$ (A16)
$`+{\displaystyle \underset{i=0}{\overset{3}{}}}{\displaystyle \underset{\nu }{}}(x^i)_{\gamma \nu }b_{\nu \beta }^{(i)}=\left[(z_0^1)_{\gamma \beta }(z_L^1)_{\gamma \beta }\right],`$ (A17)
or, in compact matrix form:

$$\left(\begin{array}{cc}𝐙_0^+𝐙_L^+,\hfill & \hfill 𝐗\end{array}\right)\left(\begin{array}{c}𝐫\\ 𝐭\\ 𝐛\end{array}\right)_\beta =\left(𝐳_L^{}𝐳_0^{}\right)_\beta .$$
(A18)
The matrices $`𝐙_0^+`$, $`𝐙_L^+`$, and $`𝐗`$, as well as the column matrix $`\left(𝐳_L^{}𝐳_0^{}\right)_\beta `$ can be computed numerically from Eqs. (A8) and (A15). Note that in order to include both leads ($`i=1,2`$) in $`𝐙_0^+`$ and $`𝐙_L^+`$, the first columns of these matrices correspond to lead 1, and the columns for lead 2 are simply appended. (Similarly, $`𝐗`$ is composed of the different wall parts, $`i=0,\mathrm{},3`$, also columnwise.)

The quantity $`\left(\begin{array}{c}𝐫\\ 𝐭\\ 𝐛\end{array}\right)_\beta `$ is a column matrix containing the unknown transmission and reflection coefficients for a particle incident in channel $`\beta `$ in lead 1, as well as the coefficients $`b_{\nu \beta }^{(i)}`$, which describe the normal derivative of $`\mathrm{\Psi }_\beta `$ on the walls. As this is the quantity which we want to solve for, we need to invert the matrix $`\left(\begin{array}{cc}𝐙_0^+𝐙_L^+,\hfill & \hfill 𝐗\end{array}\right)`$. Therefore, this matrix has to be square. This can be achieved by using appropriate cutoff numbers $`n_{\alpha _i}`$ ($`i=1,2`$) in the lead expansions \[Eq. (16)\] and $`n_{\nu _i}`$ ($`i=0,\mathrm{},3`$) in the wall expansions \[Eq. (A10)\]. Since the matrix consists of $`2n_\gamma `$ rows, these numbers have to satisfy the relation

$$\underset{i=1}{\overset{2}{}}n_{\alpha _i}+\underset{i=0}{\overset{3}{}}(n_{\nu _i}+1)=2n_\gamma .$$
(A19)
(Note that $`\nu =0,\mathrm{},n_{\nu _i}`$, therefore we actually sum over $`n_{\nu _i}+1`$ terms for the wall parts.) For best numerical accuracy, the relations of the cutoff numbers are chosen approximately equal to the length relations of the corresponding integration paths.

### 3 Cavity probability amplitude

Once the elements of the S-matrix have been obtained by the method described above, it is an easy task to calculate the probability amplitude for the electron state inside the cavity. We know that the wavefunction, which, inside the cavity, is described by Eq. (2), is zero along all the walls and that is continuous at the interface between the leads and the cavity. Thus we require that along the walls,

$$\underset{\gamma }{}c_{\gamma \beta }f_\gamma (\rho ,\phi )=0,$$
(A20)
and along the interface $`C_1`$ and $`C_2`$ between the leads and cavity

$$\underset{\gamma }{}c_{\gamma \beta }f_\gamma (\rho ,\phi )=\mathrm{\Psi }_\beta (\rho ,\phi ),$$
(A21)
where $`\mathrm{\Psi }_\beta `$ (in the leads, for incoming mode $`\beta `$) is known from the previous section. It is straightforward to solve for the $`c_{\gamma \beta }`$’s, e.g., by evaluating the above expressions at $`2n_\gamma `$ different points on $`\mathrm{\Gamma }`$, and inverting the resulting system of linear equations.


text:
# Parity Violating Measurements of Neutron Densities and Nuclear Structure

## Introduction

The size of a heavy nucleus is one of its most basic properties. However, because of a neutron skin of uncertain thickness, the size does not follow from measured charge radii and is relatively poorly known. For example, the root mean square neutron radius in 208Pb, $`R_n`$ is thought to be about 0.25 fm larger then the proton radius $`R_p5.45`$ fm. An accurate measurement of $`R_n`$ would provide the first clean observation of the neutron skin in a stable heavy nucleus. This is thought to be an important feature of all heavy nuclei.

Ground state charge densities have been determined from elastic electron scattering, see for example ref.Sick . Because the densities are both accurate and model independent they have had a great and lasting impact on nuclear physics. They are, quite literally, our modern picture of the nucleus.

In this paper we discuss future parity violating measurements of neutron densities. These purely electro-weak experiments follow in the same tradition and can be both accurate and model independent. Neutron density measurements have implications for nuclear structure, atomic parity nonconservation (PNC) experiments, isovector interactions, the structure of neutron rich radioactive beams, and neutron rich matter in astrophysics. It is remarkable that a single measurement has so many applications in atomic, nuclear and astrophysics.

Donnelly, Dubach and Sickdonnelly  suggested that parity violating electron scattering can measure neutron densities. This is because the $`Z`$boson couples primarily to the neutron at low $`Q^2`$. Therefore one can deduce the weak-charge density and the closely related neutron density from measurements of the parity-violating asymmetry in polarized elastic scattering.

Of course the parity violating asymmetry is very small, of order a part per million. Therefore measurements were very difficult. However, a great deal of experimental progress has been made since the Donnelly et. al. suggestion, and since the early SLAC experimentSLAC . This includes the Bates 12C experimentcarbon12 , Mainz 9Be experimentHeil , SAMPLEsample1  and HAPPEXhappex1 . The relative speed of the HAPPEX result and the very good helicity correlated beam properties of CEBAF show that very accurate parity violation measurements are possible. Parity violation is now an established and powerful tool.

It is important to test the Standard Model at low energies with atomic parity nonconservation (PNC), see for example the Colorado measurement in CsWieman98 ; Wieman99 . These experiments can be sensitive to new parity violating interactions such as additional heavy $`Z`$bosons. Furthermore, by comparing atomic PNC to higher $`Q^2`$ measurements, for example at the $`Z`$ pole, one can study the momentum dependence of Standard model radiative corrections. However, as the accuracy of atomic PNC experiments improves they will require increasingly precise information on neutron densitiesPollock ; Chen . This is because the parity violating interaction is proportional to the overlap between electrons and neutrons. In the future the most precise low energy Standard Model test may involve the combination of an atomic PNC measurement and parity violating electron scattering to constrain the neutron density.

There have been many measurements of neutron densities with strongly interacting probes such as pion or proton elastic scattering, see for example ref. Ray . Unfortunately, all such measurements suffer from potentially serious theoretical systematic errors. As a result no hadronic measurement of neutron densities has been generally accepted by the field. Because of the uncertain systematic errors, modern mean field interactions are typically fit without using any neutron density information.

Finally, there is an important complementarity between neutron radius measurements in a finite nucleus and measurements of the neutron radius of a neutron star. Both provide information on the equation of state of dense matter. In a nucleus, $`R_n`$ is sensitive to the density dependence of the symmetry energy. Likewise the neutron star radius depends also on the density dependence of the symmetry energy at normal and somewhat higher densities. In the future, we expect a number of improving radius measurements for nearby isolated neutron stars such as GemingaGeminga  and RX J185635-3754Nature .

We now present general considerations for neutron density measurements, discusses possible theoretical corrections and then conclude.

## General Considerations

In this section we illustrate how parity violating electron scattering measures the neutron density and discuss the effects of Coulomb distortions and other corrections. These corrections are either small or well known so the interpritation of a measurement is clean.

### Born Approximation Assymetry

The effect of the parity-violating part of the weak interaction may be isolated by measuring the parity-violating asymmetry in the cross section for the scattering of left(right) handed electrons. In Born approximation the parity-violating asymmetry is,

$$A_{LR}=\frac{G_FQ^2}{4\pi \alpha \sqrt{2}}\left[4\mathrm{sin}^2\theta _W1+\frac{F_n(Q^2)}{F_p(Q^2)}\right],$$
(1)
with $`G_F`$ the Fermi constant and $`\theta _W`$ the weak mixing angle. The Fourier transform of the proton distribution is $`F_p(Q^2)`$ while that of the neutron distribution is $`F_n(Q^2)`$ and $`Q^2`$ is the momentum transfer squared. The asymmetry is proportional to $`G_FQ^2/\alpha `$ which is just the ratio of $`Z^0`$ to photon propagators. Since 1-4sin$`{}_{}{}^{2}\theta _{W}^{}`$ is small and $`F_p(Q^2)`$ is known we see that $`A_{LR}`$ directly measures $`F_n(Q^2)`$. Therefore, $`A_{LR}`$ provides a practical method to cleanly measure the neutron form factor and hence $`R_n`$.

### Coulomb distortions

By far the largest known correction to the asymmetry comes from coulomb distortions. By coulomb distortions we mean repeated electromagnetic interactions with the nucleus remaining in its ground state. All of the $`Z`$ protons in a nucleus can contribute coherently so distortion corrections are expected to be of order $`Z\alpha /\pi `$. This is 20 % for 208Pb.

Distortion corrections have been accurately calculated in ref.cjh . Here the Dirac equation was numerically solved for an electron moving in a coulomb and axial-vector weak potentials. From the phase shifts, all of the elastic scattering observables including the asymmetry can be calculated.

Other theoretical corrections from meson exchange currents, parity admixtures in the ground state, dispersion corrections, the neutron electric form factor, strange quarks, the dependence of the extracted radius on the surface shape, etc. are discussed in reference bigpaper . These are all small. Therefore the interpretation of a parity violating measurement is very clean.

## Conclusion

With the advent of high quality electron beam facilities such as CEBAF, experiments for accurately measuring the weak density in nuclei through parity violating electron scattering (PVES) are feasible. The measurements are cleanly interpretable, analogous to electromagnetic scattering for measuring the charge distributions in elastic scattering. From parity violating asymmetry measurements in elastic scattering, one can extract the weak charge density in nuclei and from this the neutron density.

By a direct comparison to theory, these measurements test mean field theories and other models of the size and shape of nuclei. They therefore can have a fundamental and lasting impact on nuclear physics. Furthermore, PVES measurements have important implications for atomic parity nonconservation (PNC) experiments. In the future it may be possible to combine atomic PNC experiments and PVES to provide a precise test of the Standard Model at low energies.

## Acknowledgments

This work was done in collaboration with Robert Michaels, Steven Pollock and Paul Souder. It was supported in part by DOE grant: DE-FG02-87ER40365.