---
abstract: 'A detailed analysis of the wave-mode structure in a bend and its incorporation into a stable algorithm for calculation of the scattering matrix of the bend is presented. The calculations are based on the modal approach. The stability and precision of the algorithm is numerically and analytically analysed. The algorithm enables precise numerical calculations of scattering across the bend. The reflection is a purely quantum phenomenon and is discussed in more detail over a larger energy interval. The behaviour of the reflection is explained partially by a one-dimensional scattering model and heuristic calculations of the scattering matrix for narrow bends. In the same spirit we explain the numerical results for the Wigner-Smith delay time in the bend.'
address: 'Physics Department, Faculty of Mathematics and Physics, University of Ljubljana, Slovenia'
author:
- Martin Horvat and Tomaž Prosen
title: 'The bends on a quantum waveguide and cross-products of Bessel functions'
---

\#1\#2 \#1[\#1|]{} \#1[|\#1]{} \#1\#2[\#1|\#2]{} \#1[\#1 ]{} \#1\#2 \#1\#2[\#1|\#2]{}

\#1\#2\#3\#4 \#1\#2

Introduction
============

The wave propagation in bent waveguides has a long and rich history of research that dates back to Lord Rayleigh [@reyleigh:phil_mag:1897] and continues to the present days. Initially bends have been investigated in the framework of the electromagnetic theory, but more recently also quantum mechanical aspects attracted a lot of attention. The bends are popular subject of investigation, because they are typical elements incorporated into designs of waveguides. The computation of their properties to sufficiently high precision seems to be a difficult problem in the regimes of high energies and high curvatures even today.

We are discussing a bend as a scatterer of non-relativistic quantum waves on a two-dimensional ideal straight waveguide as shown in figure \[pic:schema\]. Such a structure is referred to as [*an open billiard*]{}. The past research of quantum aspects of bends can be separated into two branches. These are studies of bound states, their existence [@jensen:ann:71; @exner:jmp:89; @exner:jmp:93] and spectra [@lin:prb:96] and the scattering properties, which are both reviewed in reference [@londergan:book:99]. In order to describe quantum phenomena over our open billiard several approaches have been used in the past: Green function approach [@spivack:wrm:02], finite difference mesh calculations [@lent:apl:90] and mode-matching techniques (MMT) using natural modes i.e. eigenfunctions of the Laplacian in the bend [@cochran:radsci:66; @accatino:proc:90; @sols:prb:90; @sprung:jap:92; @lin:prb:96; @rashid:proc:02] and other bases [@amari:proc:00]. The work [@lin:prb:96] is particularly interesting as it raises the question on how to stabilise the calculations and gives a MMT method that is stable, but unfortunately a bit ambiguous. The MMT based on natural modes is called [*the modal approach*]{} and is the main topic of discussion in the present paper. The modal approach looks the most promising to deal with bends, because of its simplicity, power of interpretation and precision of results. But it also hides some problems that we examine here in detail.

![A schematic picture of a finite bend of the inner radius $r=q$ and the outer radius $r=1$ on a straight waveguide of width $a=1-q$.[]{data-label="pic:schema"}](finite_bend.eps){width="8cm"}

Let us introduce the modal approach in our open billiard composed of a bend, with the inner radius $r=q$ and the outer radius $r=1$, and a straight waveguide of width $a=1-q$ as shown in figure \[pic:schema\]. The area of the billiard, denoted by $\Omega$, can be separated into three sections: A – the left lead, B – the bend and C – the right lead. We are searching for the wave function $\psi({\bf r}\in\Omega)$, which solves the stationary Schrödinger (Helmholtz) equation on $\Omega$ with Dirichlet boundary conditions reading -([**r**]{}) = k\^2 ([**r**]{}), |\_[[**r**]{} ]{} = 0, \[eq:eigen\_problem\] where $E = k^2$ is the energy and $k$ the corresponding wavenumber. The Helmholtz equation (\[eq:eigen\_problem\]) is written in Cartesian coordinates ${\bf r} =(x,y)$ in the asymptotic regions A and C as -( + ) = k\^2 , |\_[y=0,1-q]{} = 0, \[eq:eigen\_straight\] and in polar coordinates ${\bf r} =(r,\phi)$ across the bend B as -( + + )= k\^2 , |\_[r=q,1]{} = 0. \[eq:eigen\_bend\] The modal approach suggests that we first solve Helmholtz equation (\[eq:eigen\_problem\]) on each region of the open billiard separately. Thereby we obtain partial solutions called [*mode functions*]{}, which are then used to describe the solution across the whole billiard. In the infinite straight waveguide and in the bend the mode functions are given by the ansätze $\psi \propto u(y) \exp(\ii g x)$ and $\psi \propto  U (r) \exp(\ii \nu \phi)$, respectively. The ansätze in equations (\[eq:eigen\_straight\]) and (\[eq:eigen\_bend\]) give the following equations for the mode functions: & + (k\^2 -g\^2) u = 0, &u|\_[y=0,a]{}= 0,y, \[eq:straight\_cross\_eq0\]\
& + + (k\^2-)U = 0, &U|\_[r=q,1]{}= 0, r, \[eq:bend\_cross\_eq0\] where $a=1-q$ is the channel’s width and the scalars $g,\mu\in\bC$ are called [*mode numbers*]{}. The equations (\[eq:straight\_cross\_eq0\]) and (\[eq:bend\_cross\_eq0\]) have a discrete set of solutions i.e. the mode functions and the corresponding mode numbers, which we refer to as [*modes*]{}.

The modes in the straight waveguide and in the bend are denoted by pairs $(g_n, u_n(y))$ and $(\nu_p, U_p(r))$, respectively, where $p,n\in\bN$. The modes in the straight waveguide are explicitly written as k\^2 = g\_n\^2 + ()\^2, u\_n(x) = ( n x ), \[eq:straight\_mode\] whereas modes in the bend are more complicated. They are discussed in Section \[sec:bess\], where we also show that mode numbers are either real or imaginary. The mode functions $u_n(y)$ and $U_p(r)$ are then used in the bases of functions in which we expand waves over parts of the waveguide. The basis in the straight waveguide is given by e\_n\^([**r**]{})= u\_n(y) , \[eq:mode\_fun\_straight\] and in the bend by f\_p\^([**r**]{}) = U\_p (r) , \[eq:mode\_fun\_bend\] where the sign $\pm$ labels the two directions of phase (probability) flux propagation. We define the square-root of a complex number $z=|z|\exp(\ii \phi)$, $\phi\in [0,2\pi)$ as $\sqrt{z} = |z|^{1/2} \exp (\ii \phi/2)$. The basis functions are called [*wave modes*]{} or [*modes of the Laplacian*]{}. We distinguish two types of wave modes. The wave modes corresponding to real and imaginary mode numbers are called [*open modes or travelling waves*]{} and [*closed modes or decaying (evanescent) waves*]{}, respectively. The wave-function $\psi({\bf r})$ (\[eq:eigen\_problem\]) in the entire open billiard region $\Omega$ is expressed in terms of the wave modes as ([**r**]{}) &= \_n a\^+\_n e\_n\^+ ([**r**]{})+ a\^-\_n e\_n\^- ([**r**]{}), &[**r**]{} \_[A]{}, \[eq:param\_A\]\
([**r**]{}) &= \_p \^+\_p f\_p\^+([**r**]{})+\^-\_p f\_p\^-([**r**]{}), &[**r**]{} \_[B]{}, \[eq:param\_B\]\
([**r**]{}) &= \_n b\^+\_n e\_n\^+([**r**]{})+ b\^-\_n e\_n\^-([**r**]{}), &[**r**]{} \_[C]{}, \[eq:param\_C\] where $\Omega_{A, B, C}$ are regions corresponding to sections A, B and C, respectively. The expansion coefficients $a^\pm_n$, $\lambda^\pm_p$ and $b^\pm_n$ are determined by the condition that the wave function $\psi(\vec r)$ is smooth everywhere in $\Omega$, in particular on the boundaries between different regions $\Omega_{A,B,C}$. The solution of the presented problem will be discussed in Section \[sec:scatt\].

The paper is organised as follows. In Section \[sec:bess\] we present a detailed study of the mode structure in the bend, which is closely related to the work of Cochran [@cochran:jciam:64; @cochran:pcps:66; @cochran:qjmam:66]. In comparison to the work of others, ours is directed more towards the application of the mode structure to scattering calculations. In addition, we write explicit formulae for the mode functions in the bend, where we give special attention to the closed modes. In Section \[sec:scatt\] we outline a numerically stable MMT for calculation of the scattering matrix [@newton:book:02] of a single bend. The section \[sec:scatt\] is concluded with the presentation of numerical results obtained by our method and compared to analytic estimates of the quantum transport properties of the bend. By considering the analogy between the quantum theory and EM theory we can connect our work to the EM wave propagation of longitudinal magnetic waves [@cochran:radsci:66].

The cross-product of Bessel functions {#sec:bess}
=====================================

In this section we analyse the properties of the mode numbers and the corresponding mode functions for a given wavenumber $k$ and inner radius $q$. The mode functions in the bend $U_p(r)$ are proportional to well known [*cross-products of Bessel functions*]{} [@cochran:jciam:64] of the first kind, t $J_\nu$, and Bessel functions of the second kind, $Y_\nu$, [@olver:72] written as Z\_(k,r) = J\_(kr)Y\_(k) - Y\_(kr) J\_(k), \[eq:bcp1\] or Z\_[,k]{} (r) = , , \[eq:bcp2\] where the allowed values of mode numbers $\nu$ are determined by the Dirichlet boundary conditions $Z_{\nu,k}(q) = 0$. In equation (\[eq:bcp2\]) we have used the relation $Y_\nu (z) = (J_\nu(z)\cos(\nu\pi)- J_{-\nu}(z))/\sin(\nu\pi)$ valid for orders $\nu\notin\bZ$. The understanding of the mode structure in the bend is essential for calculations of the scattering over our open billiard in the modal approach.

The properties of mode numbers
------------------------------

The set of mode numbers at a given wave-number $k\in \bR,\, k>0$ and inner radius $q\in (0,1)$ is denoted by ${\cal M}_{k,q}=\{\nu\in\bC: Z_{\nu,k}(q) = 0\}$. The functions $Z_{\nu,k}(r)$ are even $Z_{-\nu,k}(r) = Z_{\nu,k}(r)$ and analytic in the order $\nu$ [@cochran:pcps:66]. These properties yield the following symmetry of the set of mode numbers: \_[k,q]{} = - [M]{}\_[k,q]{},\_[k,q]{} \^\* = [M]{}\_[k,q]{}. \[eq:bcp\_symm0\] In addition we conclude that mode numbers are either purely real or purely imaginary \_[k,q]{} . \[eq:bcp\_symm1\] The number of real modes is finite, whereas the number of imaginary modes is infinite at a finite wavenumber $k$. The proof of the later is given in \[sec:mode\_symm\]. The properties (\[eq:bcp\_symm0\]) and (\[eq:bcp\_symm1\]) enable a decomposition of ${\cal M}_{k,q}$ into two disjoint subsets of mode numbers laying on the positive ${\cal M}_{k,q,+}$ and the negative ${\cal M}_{k,q,-}$ real and imaginary axes: \_[k,q,+]{} &=& {\_[k,q]{}: 0 0 },\
[M]{}\_[k,q,-]{} &=& {\_[k,q]{}: &lt; 0 &lt; 0 } -[M]{}\_[k,q,+]{}. \[eq:mode\_set\_minus\] It easy to see that ${\cal M}_{k,q} = {\cal M}_{k,q,+} \cup  {\cal M}_{k,q,-}$. We call $\Re\{ {\cal M}_{k,q,+}\}$ and $\Im \{{\cal M}_{k,q,+}\}$ the set of [*real modes*]{} and [*imaginary modes*]{}, respectively. The number of real modes in the bend $N_{\rm b}=\card \Re\{{\cal M}_{k,q,+}\}$ is equal or one more than the number of real modes in the straight waveguide $N_{\rm s}= \lfloor ka/\pi\rfloor$: 0 N\_[b]{} - N\_[s]{} 1 . \[eq:mode\_ineq\] where $\lfloor x\rfloor$ denotes the largest integer smaller than $x$. Taking into account the analyticity of $Z_{\nu,k}(r)$ in the order $\nu$ and in the wavenumber $k$ [@cochran:pcps:66] we find that $N_{\rm b}$ can be computed in the semi-classical limit, for $q\neq 0$, as N\_[b]{} = + + O(k\^[-2]{}), \[eq:number\_of\_modes\] and $N_{\rm s}$ is asymptotically, as k$\to\infty$, close to $N_{\rm b}$. The expressions (\[eq:mode\_ineq\]) and (\[eq:number\_of\_modes\]) are explained in \[sec:mode\_diff\]. The asymptotic form of $Z_{\nu,k}(r)$ in the order parameter [@olver:62] reads as Z\_[,k]{}(r) = \[r\^O(1) - r\^[-]{} O(1)\]. || 1. \[eq:tay\_asim\] From equation (\[eq:tay\_asim\]) we learn that $Z_{\nu,k}(q)$ diverges exponentially with increasing order parameter on the real axis as $O(q^{-|\nu|})$ and oscillates along the imaginary axis. This bounds real mode numbers $\Re\{{\cal M}_{k,r,+}\}$ from above and indicates that there is an infinite number of almost periodic imaginary mode numbers.

The mode numbers $\nu\in{\cal M}_{k,r}$ and consequently the mode functions can be in general obtained only numerically. We use different approximations of mode numbers to improve their numerical computation. By using the Debye approximation of Bessel functions for imaginary orders $\nu=\ii y$ ($y\in\bR$) valid for $y^2 + (qk)^2 \gg 1$ and the Dirichlet condition $Z_{\nu,k}(q)=0$ we obtain the following relation - + y = n, n. \[eq:WKB\_imag\] The solution of the equation (\[eq:WKB\_imag\]) in variable $y\gg k$ represents an asymptotic approximation of imaginary mode numbers $\nu_n=\ii y_n$ and is written as y\_n = - + O(n\^[-3]{}), n1. \[eq:asym\_imag\] The first term in equation (\[eq:asym\_imag\]) is already well known and can be also obtained from equation (\[eq:tay\_asim\]), see [@cochran:jciam:64]. The divergence of $Z_{\nu,k}(q)$ for $\nu\to\infty$ makes the finding of high real mode numbers, especially at large $k$, extremely difficult. We stabilize the search by using an analytic approximation of the highest real mode number $\nu_{\rm max} (k, q) = \max\{\nu: Z_{\nu,k}(q) =0\}$ for $q$ sufficiently far away from 0. This is achieved by using the asymptotic expansion of Bessel functions [@olver:72] in the transitional regime yielding \_[max]{}(k,q) = k - 3 + O(k\^[-[13]{}]{}), \[eq:nu\_bess\_asym\] where a\_1 = - 0.323685, a\_0 2.3381074. The constant $a_0$ is the negative first zero of the Airy function, ${\rm Ai}(-a_0)=0$. The exact implicit formula $Z_{\nu,k}(q) = 0$ for mode numbers at given $k$ and $q\neq 0$ has an interesting simple first order approximation reading ()\^2 - ()\^2 1, \[eq:useful\_bess\] which is asymptotically exact in two independent limits: $\nu={\rm fixed}$, $k\to\infty$, and $k={\rm fixed}$, $|\nu|\to\infty$. The relation (\[eq:useful\_bess\]) represents a useful approximation of mode numbers and is to our knowledge a new uniform approximation of modes in a bend. The expression $k_0(n,q)$ is the $n$-th zero of $Z_{0,k}(q)$, which can be easily found numerically. In the limit of large $n$, where we can use $k_0 (n,q) \approx \pi n/a$, the relation (\[eq:useful\_bess\]) in simplified to (k a)\^2 - (q)\^2 (n)\^2. \[eq:useful\_bess\_apr\] The validity of this formula is illustrated in figure \[pic:bess\_value\], where we compare mode numbers obtained from the approximate relation (\[eq:useful\_bess\_apr\]) with the exacts ones.

![Two different representations of mode numbers $\nu$ at the corresponding wavenumbers $k$ in the bend with $q=0.6$. In (a) we plot $k$ vs. $\nu$, where we separately discuss real and imaginary mode numbers depicted on the right and left side of the abscissa, respectively, and in (b) we plot $k^2$ vs. $\nu^2$. The dashed lines represent solutions of equation (\[eq:useful\_bess\_apr\]).[]{data-label="pic:bess_value"}](modes_value.eps "fig:"){width="7.5cm"}![Two different representations of mode numbers $\nu$ at the corresponding wavenumbers $k$ in the bend with $q=0.6$. In (a) we plot $k$ vs. $\nu$, where we separately discuss real and imaginary mode numbers depicted on the right and left side of the abscissa, respectively, and in (b) we plot $k^2$ vs. $\nu^2$. The dashed lines represent solutions of equation (\[eq:useful\_bess\_apr\]).[]{data-label="pic:bess_value"}](modes_value2.eps "fig:"){width="7.5cm"} to15cm[(a)(b)]{}

The highest real mode $\nu_{\rm max}$ for small wave-numbers $k < \nu$ can be approximated using equation (\[eq:useful\_bess\]) as \_[max]{}(k,q) . \[eq:nu\_bess\_min\] In practical applications it is important that below the wave-number $k_{\rm low}(q)= k_0(1,q)$ there are no real modes. In wide bends with $q \approx 0$ one can expand the cross-product of Bessel functions around $q=0$ and obtain the formula k\_[low]{}(q) = b\_0 + + O(q\^2), \[eq:k\_low\_qmin\] where $b_0$ is the smallest zero of the Bessel function $J_0(x)$, $J_0(b_0)=0$ and b\_02.404825558, b\_1 = - 1.54288974. In narrow bends, where $a (=1-q) \to 0$, we can use standard stationary perturbation theory [@morse:book:53], see equation (\[eq:schrod\_1d\]) in \[sec:mode\_diff\], to approximate $k_{\rm low}(q)$. By introducing the matrix elements V\_[nm]{} = -\_0\^1 x , = , we can express the lowest wave-number as k\_[low]{}(q) = \^[1/2]{}. \[eq:k\_low\_qmax\] The formula (\[eq:k\_low\_qmax\]) has a simple first order expansion in $a=1-q$ k\_[low]{}(q) = - + O(a\^2). \[eq:k\_low\_qmax\_exp\] From the expression (\[eq:k\_low\_qmax\_exp\]) we learn that the lowest wave-number at which real modes exists increases with increasing $q$ and converges to $\frac{\pi}{a}$.

Numerical evaluation of mode functions in a bend
------------------------------------------------

The mode functions in the bend at a given wave-number $k$ and inner radius $q$ are proportional to cross products of Bessel functions $Z_{\nu,k}(r)$ (\[eq:bcp1\]), where the order parameter $\nu$ takes values from the set ${\cal M}_{k,q}$. Because of the symmetry $Z_{-\nu,k}(r) = Z_{\nu,k}(r)$, we only consider mode numbers from the set ${\cal M}_{k,q,+}=\{\nu_n : \nu_n^2 > \nu_{n+1}^2\>,\; n \in\bN\}$, which are ordered by decreasing square. To illustrate the basic properties of mode functions, we plot in figure \[pic:bess\_fun\] the functions $Z_{\nu,k}(r)$ for real and first few imaginary mode numbers for $q=0.6$ and some low wave-number $k$. We see that the first mode function $Z_{\nu_1,k}(r)$ has no zeroes on the interval $r\in(q,1)$ and each consecutive mode function has one additional zero. In the following, we present formulae and numerical recipes for stable calculation of mode functions, where we assume that the mode numbers are given.

![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic1.eps "fig:"){height="4cm"}![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic2.eps "fig:"){height="4cm"}![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic3.eps "fig:"){height="4cm"}\
![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic4.eps "fig:"){height="4cm"}![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic5.eps "fig:"){height="4cm"}![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic6.eps "fig:"){height="4cm"}\
![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic7.eps "fig:"){height="4cm"}![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic8.eps "fig:"){height="4cm"}![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic9.eps "fig:"){height="4cm"}\
![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic10.eps "fig:"){height="4.25cm"}![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic11.eps "fig:"){height="4.25cm"}![Plots of $Z_{\nu_n,k}(r)$ for first 12 mode-numbers $\nu_n \in {\cal M}_{k,q,+}$ $(\nu_n^2 > \nu_{n+1}^2)$ calculated at the wave-number $k = 78.5398$ $(N_{\rm b}=10)$ and inner radius $q=0.6$.[]{data-label="pic:bess_fun"}](pic12.eps "fig:"){height="4.25cm"}

The Bessel functions for real orders are well implemented in the currently available numerical libraries, i.e. [@slatec:lib:93]. From the definitions of Bessel functions it is not unexpected that we encounter problems at evaluating $Z_{\nu,k}$ for large wavenumbers $k\gg 1$ and high real orders $\nu\gg 1$. In order to overcome these problems we apply the stable forward recursions in the order parameter [@olver:72] written as p\_[+ 1]{} &=& p\_[- 1]{} - 2(x q\_+ y r\_),\
q\_[+ 1]{} &=& (2(+1)y\^2 -1) r\_+ 2(+1) x y q\_- (+1)y p\_[-1]{} + x p\_,\
r\_[+ 1]{} &=& (2(+1)x\^2 -1) q\_+ 2(+1) x y r\_- (+1)x p\_[-1]{} + y p\_, where we write $x=(kr)^{-1}$ and $y=k^{-1}$ and define the following symbols p\_&=& J\_(kr) Y\_(k) - J\_(k) Y\_(kr),\
q\_&=& J\_(kr) Y’\_(k) - J’\_(k) Y\_(kr),\
r\_&=& J’\_(kr) Y\_(k) - J\_(k) Y’\_(kr). The initial conditions for the recursion, at low orders, are calculated using standard routines and the expressions for the derivatives of Bessel functions, e.g. a relation valid for any Cylindrical function: $ C'_\nu = ( C_{\nu-1} - C_{\nu+1} ) / \nu$. However, we encountered a problem at high wavenumbers and low orders due to the lack of precision in routines. Therefore in that regime we use the Hankel approximation [@olver:72] to evaluate $Z_{\nu,k}$ Z\_[,k]{}(r) = , \[eq:bcp\_hankel\_real\] where we write $A = P_\nu(kr)$, $B = Q_\nu(kr)$, $C = P_\nu(k)$ and $D = Q_\nu(k)$, which are expressed in terms of the asymptotic series: P\_(z) = 1 - + + …\[eq:bess\_hank\_P\] and Q\_(z) = - + …,= 4\^2, \[eq:bess\_hank\_Q\] which we sum up to $l_{\rm max} =\lfloor\nu/2\rfloor+1$ terms. There is another difficulty occurring at high wave-numbers, which can not be corrected. The first few real modes $\nu_p\in {\cal M}_{k,q,+}$ scale linearly with the wavenumber $k$ and functions $Z_{\nu,k}(r)$ diverge with increasing order $\nu$. Consequently, the values $Z_{\nu,k}(q)$ are exponentially sensitive on the precision of the first few mode numbers $\nu_p$: |Z\_[\_p + ,k]{}|(q) \~q\^[-|\_p|]{} (|q| ). \[eq:bcp\_nu\_asym\] This problem can not be solved completely, but only partially corrected by manually setting the values of $Z_{\nu_p,k}(r)$ to zero around the inner radius. This can be done without any real loss of precision, because the mode functions are localised near outer radius. In practice we calculate the left side of equation (\[eq:bcp\_nu\_asym\]) in a finite range arithmetic $[-m,m]$ with the maximal number $m$ (e.g. $m\approx 10^{308}$ in double precision). By considering that together with the known property $\nu_{\rm max}\sim k$ (\[eq:nu\_bess\_asym\]), we find that in practice our modal approach breaks down above some wave-number $k_{\rm break} \sim \log m /|\log q|$ and consequently bounding the number of open modes in our numerical analysis below N\_ \~ , \[eq:max\_open\_modes\] where we used the relation $\rNo\sim ka/\pi$ valid for narrow bends and high wave-numbers.

The numerical evaluation of $Z_{\nu_p,k}(r)$ at imaginary orders $\nu\in\ii \bR$ is almost unsupported in currently available numerical libraries. Therefore we have developed procedures for their evaluation ourselves and give here a summary of our work. For each regime of order parameters $\nu=\ii y$ and wavenumbers $k$ we use a different strategy to evaluate $Z_{\ii y,k}(r)$ in order to achieve an optimal precision control and a CPU time consumption. The formula for the cross-products of Bessel functions (\[eq:bcp1\]) takes for imaginary orders a simple form Z\_[y, k]{}(r) = {J\_[y]{}(k)J\_[y]{}\^\*(k r)}. \[eq:bcp\_imag\] At small wave-numbers $k$ or more generally for $k\ll y$, we use the Taylor expansion of the Bessel function [@olver:72] and rewrite equation (\[eq:bcp\_imag\]) into Z\_[y,k]{}(r) = {r\^[-y]{} u\_[y]{}(k) u\_[y]{}\^\*(kr)}, \[eq:bcp\_taylor\_imag\] where we use the series u\_(z) = \_[l=0]{}\^ , (x,y) = . \[eq:bcp\_taylor\_imag\_series\] The series (\[eq:bcp\_taylor\_imag\_series\]) is summed up to the index $l_{\rm max}= \lfloor\sqrt{z^2/\epsilon + |\nu|^2/4}\rfloor$, where $\epsilon$ is the desired accuracy of the expression.

At higher wave-numbers $k$ and orders $y \lesssim k$ we combine the backwards recursion valid for Cylindrical functions a\_[l-1]{}(z) = a\_[l]{}(z) - a\_[l + 1]{}(z), \[eq:bess\_it\] with an appropriate normalisation formula for $a_l(z)$ [@olver:72] and thereby obtain an expression for $u_{\ii y}(k)$ (\[eq:bcp\_taylor\_imag\]) given by u\_[y]{}(z) = a\_0(z) \^[-1]{}. \[eq:bcp\_rec\_imag\] The terms $a_{2l}$ in the series (\[eq:bcp\_rec\_imag\]) are given with the recursion (\[eq:bess\_it\]) started at the index $l_{\rm max}=2\lfloor(x + 1)/2\rfloor$ with initial conditions $a_{l_{\rm max}} = \eps$ and $a_{l_{\rm max}+1} = 0$, where constant $\eps$ is the smallest number supported by the CPU architecture and $x$ is determined by the equation || -1 + (1+y\^2) - y y =\
x (x\^2 + y\^2) - x ( +1 ) - x . \[eq:bcp\_rec\_eq\] The later equation (\[eq:bcp\_rec\_eq\]) is meaningful only if the right side is positive, yielding that presented approach with the iteration formula is valid only for orders $y$ below some value scaling as $O(|\log\eps|)$. By increasing the wavenumber $k$ further up and keeping orders small $y \ll k$ we can use the Hankel approximation (\[eq:bcp\_hankel\_real\]) with the order parameter $\nu=\ii y$. The asymptotic series (\[eq:bess\_hank\_P\]) and (\[eq:bess\_hank\_Q\]) in the expression (\[eq:bcp\_hankel\_real\]) are summed at least up to $l_{\rm max}=\lfloor z/2 + \sqrt{2 z^2 - y^2/2}\rfloor$ terms. At large enough wavenumbers $k$ and imaginary orders $y>k$ we can make use of the Debye approximation of Bessel functions [@erdelyi:book:55] and write Z\_[y]{}(k,r) = , \[eq:bcp\_debye\_imag\] with substitutions = , = , \^[-1]{} = 1+ , \^[-1]{} = 1+ . The expression $G(x,y)$ in formula (\[eq:bcp\_debye\_imag\]) is given in the form of an asymptotic series G(x,y) = \_[m=0]{}\^, \[eq:bcp\_debye\_G\] where polynomials $v_m(t)$ are generated by the following recursion v\_[k+1]{}(t) = (1-t) (k v\_k(t) + 2t [v’]{}\_k(t) ) + t\^[-]{} \_0\^t (1-5) \^ v\_k () . The first few $v_m(t)$ read as v\_0(t) = 1, v\_1(t) = - t, v\_2(t) = - t + t\^2,…. The formulae (\[eq:bcp\_hankel\_real\], \[eq:bcp\_taylor\_imag\], \[eq:bcp\_rec\_imag\], \[eq:bcp\_debye\_imag\]) enable a stable high precision calculation of the mode functions in the bend at imaginary mode numbers.

The overlap of mode-functions in different geometries
-----------------------------------------------------

The main ingredient in the modal description of the scattering are the overlap integrals of mode functions in the straight waveguide and in the bend. These overlap integrals “tell” about the compatibility of both scattering regions and are discussed in the following.

The cross products of Bessel functions with order parameter $\nu\in {\cal M}_{k,q,+}$ at given wavenumber $k$ and inner radius $q$ form a set of functions \_[k,q]{} = {Z\_[, k]{} (r):\_[k,q,+]{}, r}, which is complete in $L_2[q,1]$ and orthogonal w.r.t. the weight function $w(r) = r^{-1}$. The later is derived in \[sec:mode\_symm\]. The orthogonality relation for $Z_{\nu,k}\in{\cal Z}_{k,q}$ reads [@luke:book:62] \_q\^1 r w(r) Z\_[,k]{} (r) Z\_[,k]{} (r) = \_[,]{} . \[eq:ortho\_bess\] The line separating the bend and the straight wave-guide will be called [*the cross-section*]{} of our open-billiard. On the cross-section we define two different scalar products denoted by $(\cdot,\cdot)$ and $\ave{\cdot,\cdot}$, and written as (a,b) = \_q\^1 a(r) b(r), = \_0\^a y a(y) b(y). \[eq:def\_scalar\] Let us now introduce modes at some fixed wave-number $k$ and inner radius $q$ for different regions of the open billiard. In the bend, mode numbers $\nu_p$ and normalised mode functions $U_p(r)$ read as \_p \_[k,q,+]{}, U\_p(r) = , p, \[eq:def\_modes\_bend\] where we order the mode numbers so that $\nu^2_p > \nu_{p+1}^2$, and in the straight leads connected to the bend we have mode numbers $g_n$ and mode functions $u_n(x)$ defined by g\_n=, u\_n(y) = ( n y), n. \[eq:def\_modes\_straight\] The modes with real and imaginary mode numbers are called open and closed modes, respectively. The number of open modes in some geometry is denoted by $\rNo$. The overlap integrals of mode functions are given by A\_[np]{} = ,B\_[np]{} = (u\_n, U\_p), \[eq:mat\_elem\] where we use the relation $r = q + y$ between the coordinates. In figure \[pic:AB\] we show a density plot of the matrix elements, in log scale, namely $\log|A_{np}|$ and $\log|B_{np}|$, at two values of inner radii with the same number of open modes $\rNo$ in both geometries.

\

to15cm[ ]{}

We see that the matrices $A_{np}$ and $B_{np}$ have a similar form for all $q$ and $k$. This is starting at small indices with a wide area of high values of matrix elements that squeezes to almost a single intensified point at $n,p \approx \rNo$ again spreading in a triangular shape with increasing indices. The parameter $q$ has a strong influence on the shape of the area with high intensities in $A$ and $B$. In the case of small values of $q$ in contrast to larger $q$, the area of high values in matrices $A$ and $B$ covers almost the whole open-open block of indices and with crossing of the narrowing at $n,p \approx \rNo$ spreads faster with increasing indices. The shape of matrices $A$ and $B$ is similar therefore in the following we only show results for the matrix $A_{np}$. We found numerically that the area of high intensities in matrices $A$ and $B$ scales with the number of open modes $\rNo$ as $A_{np}, B_{np} \sim F (n/\rNo, p/\rNo)$, where $F$ is some a well behaved function. We demonstrate this by plotting the matrix elements $|A_{np}|$ in relative indices $n/\rNo$ and $p/\rNo$ for different numbers of open modes $\rNo$ shown in figure \[pic:AB\_amp\_k\].

to 15cm

In addition, we find numerical evidence that intensities of matrix elements $A_{np}$ and $B_{np}$ in the region of open modes and in the region of closed modes scale differently with $\rNo$ &:& |A\_[np]{}|, |B\_[np]{}| F\_[max]{} (, ), \[eq:AB\_scaling\_max\]\
&:& |A\_[np]{}|, |B\_[np]{}| F\_[min]{} ( , ). \[eq:AB\_scaling\_tails\] where $F_{\rm max}$ and $F_{\rm tail}$ are some well behaved functions. Taking into account these phenomenological findings enables a better precision control of scattering calculations. In figure \[pic:AB\_dif\_k\] we see that $F_{\rm max}$ is an envelope function for maximal values of $\rNo^{1\over 2} |A_{np}|$ and that $F_{\rm min}$ can be chosen to fit the tails of $\rNo^{2} |A_{np}|$.

We prove the scaling relation (\[eq:AB\_scaling\_tails\]) by using an asymptotic approximation of the mode function in the bend U\_p(r) = ( r ), p1, \[eq:U\_p\_asym\] which yields the following asymptotic behaviour of matrix elements A\_[np]{} (q\^2(-1)\^p + (-1)\^[n+1]{}) n p\^[-3]{} + O (p\^[-5]{}) n=, p, \[eq:AB\_tails\_asym\_p\]\
A\_[np]{} (q\^[-2]{}(-1)\^[p+1]{} + (-1)\^n) p n\^[-3]{} + O (n\^[-5]{})p= 1, n. \[eq:AB\_tails\_asym\_n\] The off-diagonal diagonal elements $A_{np}$ and $B_{np}$ decay algebraically with increasing index. The pre-factor of the decay is decreasing with increasing $q$ and is singular at $q=0$. This means that for large enough $q$ it is possible to approximately express the open modes of the bend solely in terms of open modes of the straight waveguide and vice versa, as indicated by the relation (\[eq:AB\_scaling\_max\]).

From the definition of matrix elements $A_{np}$ and $B_{np}$ (\[eq:mat\_elem\]), and completeness of the mode functions at given $k$ and $q$, it follows that $A$ and $B$ are transition matrices between the sets of mode functions in different regions, u\_n(x) = \_[p]{} B\_[np]{} U\_p(r), U\_p(r) = \_[n]{} A\_[np]{} u\_n(x), \[eq:mat\_trans\] yielding the relation A B\^T = A\^T B = . \[eq:mat\_AB\_id\] In practice we work with finite sets of modes, where the identity (\[eq:mat\_AB\_id\]) cannot hold exactly.

to13.5cm[$q=0.2$$q=0.6$$q=0.9$]{}

In figure \[pic:AB\_id\] we plot $AB^T$ for different inner radii $q$ and fixed $\rNo=100$. The mismatch from the identity (\[eq:mat\_AB\_id\]) on some sub-set of indices starting at the origin $n=1$, or $p=1$, increases with decreasing inner radius $q$. This means that the numerical calculation of scattering for smaller $q$ should be less accurate at finite dimensions. The discrepancy between $AB^T$ and the identity on (truncated) finite dimensional spaces is strongly non-uniform in indices. Before going into practical aspects of this problem, we examine the convergence of matrix elements $(AB^T)_{nm}$ to $\delta_{nm}$ with increasing number of all considered modes $N=\rNo + \rNc$ at fixed $\rNo$, where $\rNc$ is the number of closed modes. An example of such convergence is shown in figure \[pic:AB\_id\_testing\].

to13.5cm[$\rNc=1$$\rNc=100$$\rNc=500$]{}

The convergence proceeds by the standard scenario, where the agreement between $[A B^T]_{np}$ and $\delta_{np}$ propagates block-wise from low to higher indices with increasing $N$. The propagation is slow due to the triangular shaped area of high intensities in the closed-closed modes block of $A$ and $B$. The speed of propagation of accuracy to higher indices increases with increasing inner radius $q$.

The SVD decompositions [@demmel:book:97] of matrices $A$ and $B$ is useful for improving and stabilising the scattering calculations and will be used in the next section. From definitions of transition matrices (\[eq:mat\_elem\]) and completeness of mode functions we obtain ( A A\^T )\_[nn’]{} &= , ( A\^T A )\_[pp’]{} &= (U\_p, rU\_[p’]{}), \[eq:rel\_AA\]\
( B B\^T )\_[nn’]{} &= , ( B\^T B )\_[pp’]{} &= (U\_p, r\^[-1]{}U\_[p’]{}). \[eq:rel\_BB\] which we use to bound the image , a\_2 0. \[eq:AB\_norm\] These results together with $AB^T=\id$ can be used to determine the form of the SVD decomposition A = U V\^T, B = U \^[-1]{} V\^T, = {\_i}\_[i]{}, \[eq:AB\_svd\] where $U$ and $V$ are orthogonal matrices. We show here an example of the SVD decomposition of finite dimensional matrices $A$ and $B$ at $q=0.2$ and $\rNo=100$. In figure \[pic:AB\_svd\_spectra\] we show singular values and in figure \[pic:AB\_svd\_vecs\] we show density plots of the corresponding matrices $U$ and $V$, where the inner indices are ordered by decreasing magnitude of singular values.

![Singular values of matrices $A$ (left) and $B$ (right) for different number of closed modes $\rNc$ = 100, 200, 300, 400 and 500 at $q$ =0.2 and $\rNo$ = 100.[]{data-label="pic:AB_svd_spectra"}](singAq02.eps "fig:"){width="7cm"} ![Singular values of matrices $A$ (left) and $B$ (right) for different number of closed modes $\rNc$ = 100, 200, 300, 400 and 500 at $q$ =0.2 and $\rNo$ = 100.[]{data-label="pic:AB_svd_spectra"}](singBq02.eps "fig:"){width="7cm"}

We see that the relative dimension of the space, which violates the bounds of the singular spectra (\[eq:AB\_norm\]), converges with increasing space dimension $N=\rNo + \rNc$, where $\rNc$ is the number of imaginary modes. In the presented case $q=0.2$ the relative dimension is around $20\%$. It is important to see that vectors in $V$ and $U$ corresponding to singular values, which violate the bounds, have non-zero components only at closed modes. We conclude that due to a finite dimensional representation of matrices $A$ and $B$ we have deviations from the infinitely dimensional case only at high laying decaying modes that span a space of almost fixed relative dimension for some value of $q$.

The scattering across a bend {#sec:scatt}
============================

In this section we solve the on-shell scattering of a non-relativistic particle across our open billiard using the modal approach, initiated in the introduction. The scattering is discussed at fixed wave-number $k$ and inner radius $q$. The control of the precision of scattering calculations is studied in detail. In the second part of this section we investigate some interesting physical scattering properties of the bend.

The scattering matrix of a bend
-------------------------------

The scattering matrix $S$ [@newton:book:02] is a linear mapping between the incoming and outgoing “waves” with respect to our scatterer. We reorganise the expansion coefficients of the wave function over different regions (\[eq:param\_A\]), (\[eq:param\_B\]), (\[eq:param\_C\]) into the incoming contributions denoted as $ v_{\rm in} = (a^+_1,a^+_2,\ldots,b^-_1,b^-_2,\ldots)$ and outgoing contributions $v_{\rm out} = (a^-_1,a^-_2,\ldots, b^+_1, b^+_2,\ldots)$. Then the scattering matrix $S$ can be defined as S v\_[in ]{} = v\_[out]{}. \[eq:S\_symb\] The S-matrix in our case has a simple symmetric block form S = , \[eq:bend\_scatt\_def\] with $R$ and $T$ being the reflection and the transmission matrix, respectively. By reordering of rows and columns in $S$ so that the matrix elements concerning open (subscript $\rm o$) and closed (subscript $\rm c$) modes are separated and grouped together we obtain the matrix ${\cal S}$, reading S = . \[eq:S\_reorder\] We find that blocks of matrix ${\cal S}$ (\[eq:S\_reorder\]) obey the [*generalised unitarity*]{} [@prosen:jpa:95] defined by the following relations \_\_\^&=& [S]{}\_\^\_= , \[eq:S\_sim\_1\]\
[S]{}\_\_[co]{}\^&=& [S]{}\_[oc]{}, \[eq:S\_sim\_2\]\
[S]{}\_[oc]{}\^\_&=& [S]{}\_[co]{}, \[eq:S\_sim\_3\]\
[S]{}\_[co]{} [S]{}\_[co]{}\^&=& [S]{}\_[oc]{}\^\_[oc]{} = [S]{}\_[cc]{} - [S]{}\^\_[cc]{}. \[eq:S\_sim\_4\] The relations (\[eq:S\_sim\_1\] - \[eq:S\_sim\_4\]) result from the probability current conservation, which is also equivalent to the condition $AB^T=\id$. Due to the time-reversal symmetry of the physical problem the scattering matrices $S$ and ${\cal S}$ are symmetric S\^T = S,\^T = [S]{}. \[eq:S\_sim\_T\] The block symmetry of the scattering matrix (\[eq:bend\_scatt\_def\]) simplifies its calculation. We may consider individual incoming waves $e_n^+$ represented by the following wave function ansatz ([**r**]{}) &= e\^+\_n([**r**]{}) + \_m e\^-\_m([**r**]{}) R\_[mn]{}, &[**r**]{} \_[A]{}, \[eq:scatt\_A\_rt\]\
([**r**]{}) &= \_p f\^+\_p([**r**]{})\^+\_[pn]{} +f\^-\_p([**r**]{})\^-\_[pn]{}, &[**r**]{} \_[B]{}, \[eq:scatt\_B\_rt\]\
([**r**]{}) &= \_m e\^+\_m([**r**]{}) T\_[mn]{}. &[**r**]{} \_[C]{}, \[eq:scatt\_C\_rt\] The continuity of $\psi({\bf r})$ and its normal derivative on the connecting cross-sections between regions $\Omega_{\rm A,B,C}$ determines the matrix elements $R_{mn}$, $T_{mn}$ and $\Lambda^{\pm}_{pn}$ and yields the following system of matrix equations + R = M( \^+ + \^-), T = M ( [F]{} \^+ + [F]{}\^[-1]{} \^- ), \[eq:scatt\_sys\_1\]\
- R = N ( \^+ - \^- ), T = N ( [F]{} \^+ - [F]{}\^[-1]{} \^- ), \[eq:scatt\_sys\_2\]\
M = G\^ A V\^[-]{}, N = G\^[-]{} B V\^, which we write using diagonal matrices $V= \diag\{\nu_n\}_{n\in\bN}$, $G = \diag\{g_n\}_{n\in\bN}$ and ${\cal F} = \exp(\ii V)$, and transition matrices $A$ and $B$ (\[eq:mat\_elem\]). The elimination of matrices $\Lambda^\pm$ from equations (\[eq:scatt\_sys\_1\]) and (\[eq:scatt\_sys\_2\]) yields the blocks of the scattering matrix $S$, reading T = (C - [F]{} D C\^[-1]{} [F]{} D)\^[-1]{} [F]{} (C - DC\^[-1]{}D), \[eq:scatt\_T\]\
R = (C - [F]{} D C\^[-1]{} [F]{} D)\^[-1]{} ([F]{} DC\^[-1]{} [F]{}C - D), \[eq:scatt\_R\] that we express by using the following auxiliary matrices C = M\^[-1]{} + N\^[-1]{} = V\^[12]{} B\^T G\^[-[12]{}]{} + V\^[-[12]{}]{} A\^T G\^[12]{}, \[eq:mat\_C\]\
D = M\^[-1]{} - N\^[-1]{} = V\^[12]{} B\^T G\^[-[12]{}]{} - V\^[-[12]{}]{} A\^T G\^[12]{}, \[eq:mat\_D\] where we take into account the relation $AB^T=\id$. The presented form of the matrix $T$ (\[eq:scatt\_T\]) and $R$ (\[eq:scatt\_R\]) is chosen in order to increase its numerical stability i.e. minimising the use of inverses and avoiding direct computation of ${\cal F}^{-1}$.

Numerically stable scheme for scattering matrix calculation
-----------------------------------------------------------

A bend on a straight waveguide is a paradigmatic example for testing numerical schemes and ideas on how to accurately calculate the scattering matrix. In particular, the high curvature case $q\to 0$ turns to be highly non-trivial. Here we give a simple and stable procedure to obtain the scattering matrix with a clear precision control for practically all curvatures.

The scattering across a bend of angle $\beta$ and inner radius $q$ back to asymptotic region at some wavenumber $k$ is described by the scattering matrix $S(\beta)$ (\[eq:bend\_scatt\_def\]), which is composed of the reflection matrix $R(\beta)$ (\[eq:scatt\_R\]) and the transmission matrix $T(\beta)$ (\[eq:scatt\_T\]). In practice, we work with finite dimensional matrix approximations, denoted by R\_N(), T\_N() \^[NN]{}, S\_N() \^[2N2N]{}, A\_N, B\_N \^[NN]{}. where $N \ge \rNo$ is the number of modes used in the asymptotic regions. The main objective is to construct these finite dimensional matrices $R_N(\beta)$ and $T_N(\beta)$ so that:

1.  calculations are numerically stable and precise,

2.  $S_N(\beta)$ satisfies the time reversal symmetry (\[eq:S\_sim\_T\]) and the generalised unitarity relation (\[eq:S\_sim\_1\] - \[eq:S\_sim\_3\]),

3.  the sub-block of $R_N(\beta)$ and $T_N(\beta)$ of dimension $N'\in (\rNo,N]$ is calculated with controllable accuracy, where $\rNo = \lfloor k a/\pi\rfloor$ is the number of open modes in the asymptotic region.

The recipe to achieve these assumptions may be separated into two parts. In the first part we cure the numerical instability caused by the maximal element $r(\beta)= \exp(\beta\nu_N)$ in ${\cal F}_N(\beta)$, which are exponentially diverging with increasing $N$. This is achieved by separating the bend into $2^{n_0}$ identical subsections of angle $\beta'= 2^{-n_0}\beta$ so that $r(\beta')\eps_{\rm num}\approx 1$, where $\eps_{\rm num}$ is the numerical precision e.g. $\eps_{\rm num} = 2^{-52}$ in double precision floating point arithmetic. The number $n_0$ is calculated as n\_0 = {0, 1 + \_2 ( )}. \[eq:comb\_n\_steps\] The scattering matrix $S_N(\beta')$ of a small subsection of the bend can be calculated in a very stable way and up to a high precision. By concatenating scattering matrices of subsections together with the recursion S\_N(2\^[-m+1]{}) = S\_N(2\^[-m]{}) S\_N(2\^[-m]{}). \[eq:comb\_recurs\] we obtain the scattering matrix of the whole bend $S_N(\beta)$. The symbol $\odot$ denotes the [*operation for concatenating scattering matrices*]{} associated to scatterers on the waveguide and is defined in \[sec:comb\_scatt\].

In the second part we are discussing the problem that $S_N(\beta')$ diverges with increasing $N$, which is due to violation of the identity $AB^T = \id$ (\[eq:mat\_AB\_id\]) for finite truncated transition matrices $A_N$ and $B_N$. We eliminate the problem by deforming $A_N$ and $B_N$ so that they are non-singular and exactly fulfil the condition $A_N B_N^T = \id$. We make the SVD decomposition of the truncated matrix $A_N = U_N \Sigma_N V^T_N$, modify its singular values $\Sigma_N = \diag\{\sigma_i\}_{i=1}^N$ to $\tilde \Sigma_N = \diag\{\tilde\sigma_i\}_{i=1}^N$ so that they fit in the bounds obtained for infinitely dimensional case (\[eq:AB\_svd\]) \_i = {

[lll]{} \_i&:& \_i \^\*&:& [otherwise]{}

. , \[eq:correct\_sigma\] and again generate both matrices A\_N = U\_N \_N V\^T\_N, B\_N = U\_N \_N\^[-1]{} V\^T\_N. \[eq:correct\_A\] The same procedure can also be done using SVD decomposition of the matrix $B_N$ as a base for generation of both deformed matrices $\widetilde A_N$ and $\widetilde B_N$. The value of $\sigma^*>0$ can be chosen arbitrarily, but the most elegant choice is $\sigma^*=1$. By using matrices $\widetilde A_N$ and $\widetilde B_N$ (\[eq:correct\_A\]) instead of $A_N$ and $B_N$ in $S_N(\beta')$ and consequently in $S_N(\beta)$ these become generalised unitary with a well behaved and physically precise limit $N\to\infty$ at least on the sub-space of dimension $N'$. We first check this by discussing the precision of transition between the modes in the bend and in the straight waveguide on the sub-space of dimension $N'<N$. The [*error of transition*]{} from the asymptotic region (infinite waveguide) into the bend is quantified by \_[[s]{}]{}(N,N’) = \_[n, m L]{} | | = \_[n, m L]{} |\_[p=1]{}\^N A\_[np]{} B\_[mp]{} - \_[nm]{}|, \[eq:id\_eps\_sb\] and for transition in the opposite direction by \_[[b]{}]{}(N,N’) = \_[p, r L]{} |(U\_p| |U\_r)| = \_[p, r L]{} |\_[n=1]{}\^N A\_[nr]{} B\_[np]{} - \_[pr]{}|. \[eq:id\_eps\_bs\] with $L= \{1,2,\ldots N' \}$. The introduced transition errors (\[eq:id\_eps\_sb\]) and (\[eq:id\_eps\_bs\]) measure the violation of the identity $A B^T =\id$ on the subspace $N'$, when working with the finite number of modes $N$. It is expected and supported by our numerical studies that the errors vanish in the limit $N\to\infty$ at fixed $N'$ and other parameters. In figure \[pic:id\_eps\] we show transition errors as a function of $N$ at a fixed $N'$ for two values of $\rNo$.

![The transition errors $\eps_{{\rm s}\to {\rm b}}(\rNo+\rNc,N')$ (top) and $\eps_{{\rm b}\to {\rm s}}(\rNo+\rNc,N')$ (bottom) calculated on the sub-space of dimension $N'=\rNo+10$ vs. the number of closed modes $\rNc$ for $\rNo=10,100$ (left, right).[]{data-label="pic:id_eps"}](ids_keys.eps){width="1.5cm"}

The errors decrease with increasing $\rNc$ down to a certain plateau, which is determined by the precision of mode functions. The most problematic in the precision are the highest few open modes in the bend. A conservative estimate for the plateau is around $10^{-8}$. The transition errors (\[eq:id\_eps\_sb\]) and (\[eq:id\_eps\_bs\]) increase with decreasing $q$ indicating that we need more modes to achieve equally small error as for higher $q$. We did not find any analytic approximation for transition errors. Therefore we numerically estimate the minimal dimension of the functional space $N_{\rm t} = \rNc + \rNo$ needed for transition errors to be smaller than some $\eps$, defined as N\_[t]{} (N’, ) = { N : \_[[s]{}]{} (N,N’) &lt; \_[[b]{}]{} (N,N’) &lt; }, \[eq:N\_geom\] where $N'= \rNc' + \rNo$ is the dimension of the observed sub-space. In our numerical analysis we set $\eps \approx 10^{-7}$. We check the convergence of $R_N(\beta)$ and $T_N(\beta)$ with increasing $N$ on some fixed sub-space of dimension $N'<N$. The matrices $R_N$ and $T_N$ are calculated using the method of dividing the bend into sub-sections together with deforming the transition matrices. The convergence is measured through the relative difference of matrices at subsequent changes of the dimension $N$ \_[R]{} (N,N’) = , \_[T]{} (N,N’) = , \[eq:conv\_cauchy\_measure\] where we introduce a matrix norm $\| A \|_M = \max_{i,j\in [1, M]} |A_{ij}|$ on the sub-space of dimension $M$. The expressions (\[eq:conv\_cauchy\_measure\]) give upper bounds for the deviations of matrices from their asymptotic forms R\_() - R\_N() \_[N’]{} C\_1(N’) \_[M=N+1]{}\^\_[R]{}(M,N’), \[eq:conv\_abs\_measure\_1\]\
T\_() - T\_N() \_[N’]{} C\_2(N’) \_[M=N+1]{}\^\_[T]{}(M,N’), \[eq:conv\_abs\_measure\_2\] with expressions $C_1 (N')=\max_{N\ge N'} \|R_N(\beta)\|_{N'}$ and $C_2 (N')=\max_{N\ge N'} \|T_N(\beta)\|_{N'}$, which are of the order of magnitude 1. We have numerically studied quantities $\eps_{\rm R, T}(N,N')$ as functions of $N$ at fixed $N'$ and the results are shown in figure \[pic:RT\_convergeX\].

![The measures of convergence $\eps_{\rm R,T}(\rNo+\rNc,\rNo+10)$ (top, bottom) as functions of $\rNc$ for $\rNo = 10,100$ (left, right) at various inner radii $q$ as indicated in the figure. The bend is of angle $\beta=\pi$.[]{data-label="pic:RT_convergeX"}](RT_conv_key2.eps){width="1.5cm"}

In the case $N = N' = \rNo$ we are talking about the open-open block of the scattering matrix, which is sometimes called the semi-quantal approximation [@bogomolny:nonlin:92]. We see that $\eps_{\rm R}$ and $\eps_{\rm T}$ decrease with increasing total number of modes $N=\rNo+\rNc$ in a similar fashion as transition errors down to some plateau around $10^{-12}$. The plateau is almost equal to the machine precision, which is surprisingly better than the transition errors.

The presented method for the calculation of the scattering matrix and its accuracy control works well for all $q>0$. Nevertheless, a treatment of high curvature cases $q\to 0$ are difficult as we need to consider a large number of closed modes to reach a sufficient precision of the scattering matrix. However, this is feasible to achieve by our method in a stable and controlled way.

Quantum transport across the bend
---------------------------------

The scattering matrix of the open billiard describes the stationary quantum transport of a particle over the bend. We discuss here the most important and obvious measures of the transport, which are the reflection probability and the [*Wigner-Smith delay time*]{}. Since the classical particle cannot scatter back [@horvat:jpa:unidir04] we are particularly interested in the reflection probability as a genuine quantum (wave) property. This is a measure of quantum tunnelling between the two classically invariant components of the phase space corresponding to right and left going waves. The scattering of an incoming Gaussian ray over the bend is illustrated in figure \[eq:qm\_bend\_figs\]. The ray follows the classical trajectories, but as it is of finite width, its parts are scattered differently when hitting the curved wall. The parts of the ray travel different lengths and interfere among themselves.

We are discussing the scattering over the bend at some fixed wavenumber $k$ and inner radius $q$. The scattering properties are contained in the transmission matrix $T$ (\[eq:scatt\_T\]) and the reflection matrix $R$ (\[eq:scatt\_R\]). The wavefunction over the asymptotic region is described in $N$ modes, where $N \ge \rNo = \lfloor ka/\pi \rfloor$. We consider an incoming wave coming to the bend from the left side written in the asymptotic region as \_[in]{}([**r**]{}) = \_[n=1]{}\^N a\_n e\_n\^+([**r**]{}). \[eq:incoming\] By introducing a vector of complex coefficients $a = \{a_n\}_{n=1}^{\rNo}$ we can write the transmitted and reflected probability flux, $j_{\rm T}$ and $j_{\rm R}$ respectively, in an elegant form j\_[R]{} = a\^a, j\_[T]{} = a\^a, j\_0 = a\^a = j\_[R]{} + j\_[T]{}, \[eq:ref\_trans\_prop\_current\] where we introduce matrices $\Pi$ and $\Sigma$ calculated from open-open mode blocks of $R$ and $T$: = R\_\^R\_, = T\_\^T\_, + = . \[eq:ref\_trans\_amp\_rel\] The average transport properties are given by the first and the second moment of the probability currents $j_{\rm R, T}$ averaged over an ensemble of incoming states $a$. The ensemble represents states (vectors $a$) uniformly distributed over the 2N-dimensional sphere of radius $j_0$ [@newton:book:02; @prosen:jpa:02]. The [*average probability currents*]{} are given by = = {}, = 1 - [R]{}, and the [*standard deviations of probability currents*]{} (giving fluctuations within an ensemble) are written as \^2\_[R]{} = = = \^2\_[T]{}. In the following we thus discuss the average reflection $\cal R$ and the dispersion of reflection $\sigma^2_{\rm R}$. An approximation of the transmission matrix $T$ can be determined from the semi-classical calculations, whereas for the reflection matrix $R$ it can not, as the reflection in the bend is a purely quantum phenomenon. The gross structure of matrices $R$ and $T$, similarly as of matrices $A$ and $B$, does not change significantly with increasing wavenumber $k$. In figure \[pic:S\_matrix\] we show the density plot of matrices $R$ and $T$ with $\rNo=100$ open modes. The high probabilities in the matrix $T$ have a classical correspondence, which is revealed through the calculation of the classical scattering matrix $T_{\rm class}$ [@mendez:pre:02]. Both, the classical and the quantum transmission matrices feature similar patterns, but due to the quantum interference, we can not establish a clear correspondence. In the matrix $T$ we have a large area of high values so we can expect that transmission probability of individual modes should be high. It is important to notice the area in the reflection matrix of high intensity is concentrated around the last open mode with the index $\rNo$.

At the so called resonant wavenumbers $k_m = \frac{\pi}{a} m$ ($m\in\bN$) a new open mode appears in the asymptotic region and causes a strong increase in the reflection matrix elements at open modes with high indices. This is demonstrated in figure \[pic:S\_matrix\_res\], where we show the scattering matrices $R$ and $T$ around the highest open mode calculated calculated at $k\approx k_{100}$ and at $q=0.6$. The changes are centred around the index $\rNo$ and significantly influence the average transport.

In order to clarify the contributions to the total reflection we plot in figure \[pic:R\_mat\_elem\] the reflection probability of individual modes $\Pi_{nn}$ for wavenumbers near and far from the resonance. We see that the highest open mode has the strongest reflection and the reflection probability of “all” modes increases at resonant wavenumbers $k=k_{\rNo}$. In particular the reflection of highest open mode is almost perfect $\Pi_{\rNo, \rNo}\approx 1$. In the vicinity of the resonance we could effectively approximate the average reflection as ${\cal R}\approx \Pi_{\rNo,\rNo}/\rNo$ (see figure \[pic:reflect\_res\]).

![The diagonal matrix elements $\Pi_{nn}$ at wave-numbers $k=r\frac{\pi}{a}$ as indicated in the figure, and fixed inner radius $q=0.6$. Note that $\rNo = \lfloor r\rfloor$.[]{data-label="pic:R_mat_elem"}](tr_mat_elem.eps){width="7.5cm"}

The resonant wavenumbers $k_{\rNo}$ are important markers for anomalously strong reflection. This is illustrated in figure \[pic:reflect\_res\], where we plot the average reflection $\cal R$ as a function of the wave-number around $k=k_{100}$. We see the average reflection $\cal R$ has a strong sharp maximum at resonant wavenumbers and decreases in an irregular oscillating manner with increasing wavenumber until crossing the next resonant wavenumber. The frequency of irregular oscillations increases with increasing wavenumber. From numerical results we see that $\cal R$ decreases with increasing $q$. In narrow channels $a\to 0$ at large wave-numbers we showed with a perturbative approach (see \[sec:pert\_scatt\]) that \~, which is confirmed numerically. The resonant behaviour around the resonant wavenumber can be partially explained by neglecting all open modes expect the one with the highest index $\rNo$. Such system can be treated as an independent 1d scatterer ($d=1$) with the reflection and transmission matrix elements reading R\_[1d]{} = -, T\_[1d]{} = , \[eq:R\_T\_resonance\] with the phase shift \_ = 2( K\_ ), K\_p = . \[eq:R\_T\_phase\_shift\] At $k= k_{\rNo}$, the mode number $g_{\rNo} \sim \mu_{\rNo}$ and the phase shift $\mu_{\rNo}$ become zero yielding a perfect reflection in a 1d scattering model, $R_{\rm 1d} = -1$ and $T_{\rm 1d}= 0$. This treatment is meaningful, because the matrices $A$ and $B$ are approximately diagonal at $(\rNo,\rNo)$ with an algebraic decay of matrix elements when we move away from the diagonal. If the modes were strictly independent we would have $\Pi_{\rNo, \rNo}=\|R_{\rm 1d}\|^2$, but the algebraic tails in matrices $A$ and $B$ make this solution to hold only as a rough approximation as can be seen in figure \[pic:reflect\_res\].

to16cm[(a)(b)]{}

In figure \[pic:R\_gross\] we study measures of reflection $\cal R$ and $\sigma_{\rm R}$ over a larger range of wave-numbers $k$. From figure \[pic:R\_gross\].a we see that $\cal R$ strongly oscillates with peaks at resonant wave-number $k_n$ and its upper bound decreases proportionally to $k^{-1}$, as predicted. The numerical results in figure \[pic:R\_gross\].b indicate that $\sigma_{\rm R} < {\cal R}$ and ${\cal R} \sim \sigma_{\rm R}$ as $k$ goes to infinity.

![The average reflection $\cal R$ and the relative deviation of reflection $\sigma_{\rm R}$ over a larger interval of wave-numbers $k$ at $q=0.6$. []{data-label="pic:R_gross"}](tr_bent_R.eps "fig:"){width="7cm"}![The average reflection $\cal R$ and the relative deviation of reflection $\sigma_{\rm R}$ over a larger interval of wave-numbers $k$ at $q=0.6$. []{data-label="pic:R_gross"}](tr_bent_reldev.eps "fig:"){width="7cm"}

To get rid of oscillations and get an overall average behaviour of $\cal R$ and $\sigma_{\rm R}$ we calculate their cumulative integrals with respect to the wave-number. The results are shown in figure \[pic:R\_int\_gross\] and yield the following dependence \_[k\_0]{}\^k [R]{} \_[k\_0]{}\^k \_[R]{} = O((k)). \[eq:int\_measures\_asymp\] This indicates together with previous conclusions that the reflection measures, averaged over small wave-number ranges, indeed scale as \~\_[R]{} = O(k\^[-1]{}),k. \[eq:measures\_asymp\] It seems that this relation (\[eq:measures\_asymp\]) is valid for an arbitrary inner radius $q$ and represents a new and very useful information for the study of wave-guides and general billiards that include bends.

![The integral of the average reflection $\cal R$ and deviation $\sigma_{\rm R}$ over larger interval of wave-numbers $k$ of a bend with a radius $q=0.6$ and an angle $\beta=\pi$. The inserted lines are $f_1(k) = 0.01103 + 0.01022 \log k$ and $f_2(k) = -0.01553 + 0.006924 \log k$.[]{data-label="pic:R_int_gross"}](tr_bent_intRdev.eps){width="7.5cm"}

Another insight into the scattering properties gives the Wigner-Smith delay time $\tau_{\rm ws}$ [@wigner:phys_rev:55; @smith:phys_rev:60], which is the quantum analog of the geometric length travelled by a wave. In the semi-classical limit, where we could apply geometric optics, $\tau_{\rm ws}$ is equal to the average geometric length of classical trajectories over the open billiard. By using the hermitian variant of the [*lifetime matrix*]{} Q = S\_\^, \[eq:lifetime\_matrix\] the [*Wigner-Smith delay time*]{} is defined as \_[ws]{} = {Q} = {R\_\^} + {T\_\^} , where we have used block symmetries of our matrix $S$ (\[eq:bend\_scatt\_def\]). $\tau_{\rm ws}$ can be thought of as an average delay time corresponding to particular modes, which are defined as \^n\_[ws]{} = \_[m=1]{}\^ {\_[nm]{} \_[mn]{} + \_[nm]{} \_[mn]{}}, \[eq:time\_part\_mode\] with derivative defined as $(\bullet)' = \dd/ \dd k$. Numerical results shown in figure \[pic:wigner\_time\] point to a similar dependence of Wigner-Smith delay time $\tau_{\rm ws}$ on the wave-number $k$ as the average reflection $\cal R$, just the oscillations are smoothed out. The time $\tau_{\rm ws}$ strongly increases near the resonance wave-number $k=k_{\rNo}$ due to intense changes in the scattering matrices $R$ and $T$ in the area around the index of the newly open mode. Qualitatively we can explain the singular behaviour by treating the highest open mode in the resonance regions within 1d scattering model in which the delay time is given by \^[1d]{}\_[ws]{} = {T\_[1d]{}\^\* T\_[1d]{}’ + R\^\*\_[1d]{} R\_[1d]{}’} = . \[eq:time\_one\_mode\] The first term in the numerator of equation (\[eq:time\_one\_mode\]) corresponds to the transmission and the second term to the reflection. By slowly increasing the wavenumber across the region of the resonance we can notice three different regimes: before, in the vicinity and after the resonant wave-number. Slightly before the resonance $k < k_{\rNo}$, a new real mode appears in the bend (see formula (\[eq:Nb\_asym\])) making the propagation across the bend very slow. From the formula (\[eq:time\_one\_mode\]) we learn that this results in a large transmission time and consequently in a large time delay $\tau^{\rm 1d}_{\rm ws}$. In the instance of crossing the reflection resonance a new mode appears in the straight wave-guide, which causes a square-root singularity $\tau^{\rm 1d}_{\rm ws}\sim (k-k_{\rNo})^{-{1\over 2}}$ for $k> k_{\rNo}$ and its sign is determined by $2\beta \nu_{\rNo}$. This reflection term has a short-scale influence to the behaviour of the time delay and can enhance or reduce its size. Obviously, this is a very non-classical situation. By going further away from the resonance wave-numbers the reflection contribution to the time delay is levelled by an increasing transmission term due to a very slow propagation of the mode in the asymptotic region, which again increases the transition time. So we can experience one or two peaks of the time delay in the vicinity of the reflection resonance. Away from the reflection resonance the time delay drops even below the classical time. The latter we assume is due to reflection phenomena which reduces the classically expected phase shift. The presented 1d scattering model has only an instructive purpose and does not represent any useful quantitative approximation, similarly as was the case in the discussion of the reflection.

Conclusions
===========

We present mathematical, numerical and physical background of the non-relativistic 2D scattering of a quantum particle on a circular bend connected to infinite straight waveguides. We discuss mathematical properties and derive numerical recipes for accurate and reliable calculation of the mode functions and the corresponding mode numbers in a bend. We take a special care of closed (evanescent) modes in the bend. The obtained modal structure and its properties are incorporated in a robust and stable numerical scheme for computing the scattering matrix with a controllable precision. Our numerical apparatus is applied to the study of transport properties. We focus mainly on the reflection, which is a purely quantum (wave) phenomenon since the back-reflection of classical rays is not possible. Our study is particularly focused on the possibility to investigate the (semi-classical) regime of very large wave-numbers. Some of the obtained physical properties of the scattering problem can be explained analytically. In addition, we present results on the Wigner-Smith delay in the bend. The obtained transport properties can be useful in discussing and predicting properties of open billiards (or wave guides) composed of arbitrary combination of bends and straight segments.

Acknowledgements {#acknowledgements .unnumbered}
================

MH would like to thank Prof. Dr. Nico Temme for references on literature considering cross-products of Bessel functions at imaginary orders. Useful discussions with M. Žnidarič as well as the financial support by Slovenian Research Agency, grant J1-7347 and programme P1-0044, are gratefully acknowledged.

References {#references .unnumbered}
==========

[10]{}

Strutt J W 1897 On the passage of electric waves through tubes or the vibrations of dielectric cylinders [*Phil. Mag.*]{} [**(Ser. 5) 53**]{} 125–132

Jensen H and Koppe H 1971 Quantum mechanics with constraints [*Ann. Phys.*]{} [**63**]{} 586–91

Exner P and Seba P 1989 Bound states in curved quantum waveguides [*J. Math. Phys.*]{} [**30**]{} 2574–80

Exner P 1993 Bound states in quantum waveguides of slowly decaying curvature [*J. Math. Phys.*]{} [**34**]{} 23–28

Lin K and Jaffe R L 1996 Bound states and threshold resonances in quantum wires with circular bends [*Phys. Rev. B*]{} [**54**]{} 5750–62

Londergan J T and Carini J Pand Murdock D P 1999 [*Binding and scattering in two-dimensional systems: application to quantum wires, waveguides and photonic crystals, Lecture Notes in Physics, volume [**60**]{}*]{} (Berlin \[etc\]: Springer Verlag)

Spivack M, Ogilvy J and Sillence C 2002 Electromagnetic propagation in the curved two-dimensional waveguide [*Waves Random Media*]{} [**12**]{} 47–62

Lent C S 1990 Transmission through a bend in an electron waveguide [*Appl. Phys. Lett.*]{} [**56**]{} 2554–6

Cochran J A and Pecina R G 1966 Mode propagation in continuously curved waveguides [*Radio science*]{} [**1**]{} 679–96

Accation L and Bertin G 1990 Modal analysis of curved waveguides [*in Proc.20th Eur, Microwave Conf. Budapest sep. 1990*]{}

Sols F and Macucci M 1990 Circular bends in electron waveguides [*Phys. Rev. B*]{} [**41**]{} 11887–91

Sprung D W L and Wu H 1992 Understanding quantum wires with circular bends [ *J. Appl. Phys.*]{} [**71**]{} 515–7

Rashid M A and Kodama M 2002 Analysis of propagation properties in junctions between straight and bend waveguides using cylindrical functions of complex orders [*[*in Proc. ITC-CSCC-2002 Conference July, 2002, Phuket, Thailand*]{}*]{}

Amari S and Bornemannm J 2000 Modelling of propagation and scattering in waveguide bends [*in Proc. 30th European Microwave Conf. Paris, France, Oct. 2000*]{} [**Vol. 2**]{} 353–6

Cochran J A 1964 Remarks on the zeros of cross-product [B]{}essel functions [ *J. Soc. Indust. Appl. Math.*]{} [**12**]{} 580–7

Cochran J A 1966 The analyticity of cross-product [B]{}essel function zeros [ *Proc. Camb. Phil. Soc*]{} [**62**]{} 215–56

Cochran J A 1966 The asymptotic nature of zeros of cross-product [B]{}essel function [*Quart. Journ. Mech. and Applied Math.*]{} [**19**]{} 511–22

Newton R G 2002 [*Scattering Theory of Waves and Particles*]{} (Mineola, New York: Dover Publications, inc.)

Olver F W J 1972 [B]{}essel functions of integer order [*in Handbook of Mathematical Functions 10th ed., eds. Abramowitz M and Stegun I A (New York: Dover)*]{} pp 355–389

Olver F W J 1962 Tables for [B]{}essel functions of moderate or large orders [ *in Mathematical Tables vol. 6 (Her Majesty’s Stationary office, London, England)*]{}

Morse P M and Feshbach H 1953 [*Methods of Theoretical physics, vol. 1*]{} (New York: McGraw Hill)

Fong K W, Jefferson T H, Suyehiro T and Walton L 1993 [S]{}[L]{}[A]{}[T]{}[E]{}[C]{} common mathematical library, version 4.1

Erdélyi A 1955 [*Higher transcendental functions, vol. [I]{}[I]{}*]{} (New York, Toronto, London: McGraw-Hill)

Luke Y L 1962 [*Integrals of [B]{}essel functions*]{} (New York, Toronto, London: McGraw-Hill, cop.)

Demmel J W 1997 [*Applied numerical linear algebra*]{} (SIAM)

Prosen T 1995 General quantum surface-of-section method [*J. Phys. A: Math. Gen*]{} [**28**]{} 4133–55

Bogomolny E B 1992 Semiclassical quantization of multidimensional systems [ *Nonlinarity*]{} [**5**]{} 805–66

Horvat M and Prosen T 2004 Uni-directional transport properties of a serpent billiard [*J. Phys. A: Math. Gen*]{} [**37**]{} 3133–3145 ([*Preprint nlin.CD/0601055*]{})

Prosen T and Žnidarič M 2002 Stability of quantum motion and correlation decay [*J. Phys. A: Math. Gen.*]{} [**35**]{} 1455–81

Méndez-Bermudéz J A e a 2002 Understanding quantum scattering properties in terms of purely classical dynamics: Two-dimensional open chaotic billiards [*Phys. Rev. E*]{} [**66**]{} 046207

igner W P 1955 Lower limit for the energy derivative of the scattering phase shift [*Phys. Rev.*]{} [**98**]{} 145–7

Smith F T 1960 Lifetime matrix and the collision theory [*Phys. Rev.*]{} [**118**]{} 349–56

Mayer A and Vigneron J P 1999 Accuracy-control techniques applied to stable transfer-matrix computations [*Phys. Rev. E*]{} [**59**]{} 4659–65

The symmetry of mode numbers {#sec:mode_symm}
============================

We prove the symmetry (\[eq:bcp\_symm1\]), by changing the variable $r = e^{-x}$ and transforming the Bessel equation (\[eq:eigen\_bend\]) and the corresponding boundary condition into the equation + (k\^2 e\^[-2x]{} -\^2)Z = 0, Z|\_[x=0,q]{} = 0, which can be interpreted as a one-dimensional quantum mechanical eigenvalue problem, with the Hamiltonian $\hat H$ and potential $V'$: H Z = -\^2 Z, H = - + V’(x), V’(x) = {

[lll]{} -k\^2 e\^[-2x]{} &:& x &:&

.. \[eq:schrod\] Because the Hamiltonian $\hat H$ is a Hermitian operator, the eigenvalues are real $-\nu^2\in\bR$ yielding $\nu\in\bR\cup\ii\bR$. From the form of the potential $V'(x)$, depicted in figure \[pic:bess\_pot\], we see that there are only finite number of real mode numbers $\nu^2>0$ and infinite number of imaginary mode numbers $\nu^2<0$. The independent solution of equation (\[eq:schrod\]) are orthogonal with respect to the measure $\dd x = r^{-1}\,\dd r$ and so the weight function between mode functions $U_p(r)$ in the bend is $w(r) = r^{-1}$.

![The analog of the quantum potential in the eigenvalue equation for the mode functions in the bend.[]{data-label="pic:bess_pot"}](qpoten.eps){width="7cm"}

The number of modes in the straight and the bent waveguide {#sec:mode_diff}
==========================================================

We discuss the number of open modes $N_{\rm b}(k,q)=\card \Re\{ {\cal M}_{k,q,+}\}$ in the bend and its deviation from the number of modes in the straight waveguide $N_{\rm s} = \lfloor ka/\pi\rfloor$. The mode numbers in ${\cal M}_{k,q,+}$ continuously slide with increasing $k$ and fixed $q$ from the imaginary to the real axis by crossing the point $\nu=0$. This dynamics is depicted in figure \[pic:bess\_value\]. This means that $N_{\rm b}(k,q)$ is equal to the number of zeros $x$ of $Z_{0,x}(q)$ up to the value $k$ N\_[b]{}(k,q) = {[M]{}\_[k,q,+]{} } = { x k: Z\_[0,x]{}(q) = 0 }. By using substitutions $U = r^{-{1\over 2}} \varphi$ and $r = q + a x$ the mode problem in the bend (\[eq:bend\_cross\_eq0\]) is in the case $\nu=0$ transformed into a 1d stationary Schrödinger equation $(\gamma=q/(1-q))$ - + V(x) = e , V(x) = -\[2 (x + )\]\^[-2]{},x, \[eq:schrod\_1d\] with the eigen-energy dentoed by $e=(ka)^2$. The discrete set of eigen-energies is ordered as $e_{n+1} > e_n$, $n\in\bN$. By setting $V=0$ in expression (\[eq:schrod\_1d\]) we obtain the mode problem for appropriately rescaled straight waveguide. In the eigen-energies in this case se $(\pi n)^2$, $n\in\bN$. Then taking into account the (empirical) fact $e_n \le (\pi n)^2 < e_{n+1}$ we can conclude 0 N\_[b]{}(k,q) - N\_[s]{}(k,q) 1, q (0,1). This means that at certain $k$ and $q$ we can have in the bend one open mode more, but not less than in the straight waveguide. In the semi-classical limit $k\to\infty$ the eigen-energies $e_n$ can be obtained using the Debye approximation valid for $qk\gg 1$. In this way we get a relation between the eigenvalues $e$ and its counting number $N_{\rm b}$ 2 N\_[b]{} &=& (4 k\^2 + 1)\^[12]{} - (4(q k)\^2 + 1)\^[12]{}\
&-&((4k\^2 + 1)\^[-[12]{}]{} ) + ((4(q k)\^2 + 1)\^[-[12]{}]{}), \[eq:WKB\_real\_0\] which yields with asymptotic expansion in $k$ the expression N\_[b]{}(k,q) = + + + O(a(qk)\^[-3]{}). \[eq:Nb\_asym\] We see that $N_{\rm b}$ and $N_{\rm s}$ are close to each other for high wave-number and and not too small inner radius $q$.

The method of concatenating scattering matrices {#sec:comb_scatt}
===============================================

Here we outline a method to concatenate the scattering matrices [@mayer:phys_rev:99] associated to scatterers on sectioned wave guides. Let us assume to have two scatterers labelled by A and B and with scattering matrices $S_{\rm A}$ and $S_{\rm B}$, respectively. S\_[A,B]{} = \^[2N2N]{}. By combining both scatterers A and B in the order AB we build a “larger” scatterer with the scattering matrix $S$. The matrix $S$ is calculated from matrices $S_{\rm A, B}$ by a nonlinear operation $\odot: \bC^{2N\times2N}\times \bC^{2N\times2N}\to \bC^{2N\times2N}$ defined as S = S\_[A]{} S\_[B]{} = \^[2N2N]{}, \[eq:S\_merge\] which explicitly reads &r\^= r\_[A]{}\^+ t\_[A]{}\^r\_[B]{}\^L\^[-1]{} t\_[A]{}\^, &t\^= t\_[B]{}\^L\^[-1]{}t\_[A]{}\^, \[eq:S\_merge\_1\]\
&r\^= r\_[A]{}\^+ t\_[A]{}\^r\_[B]{}\^\^[-1]{} t\_[A]{}\^, &t\^= t\_[B]{}\^\^[-1]{}t\_[A]{}\^,\[eq:S\_merge\_2\] where we define $L = 1 - r_{\rm A}^\rR r_{\rm B}^\rL$ and $L' = 1 - r_{\rm B}^\rL r_{\rm A}^\rR$. Note that a bend on a straight waveguide can be treated as a scatterer. By combining bends of angles $\gamma'$ and $\delta'$ with scattering matrices $S(\gamma')$ and $S(\delta')$, respectively, we get a bend of angle $\gamma'+\delta'$ with the scattering matrix $S(\gamma'+\delta')$. The latter matrix can be obtained from matrices $S(\gamma')$ and $S(\delta')$ by the formula S(’+’) = S(’) S(’) = S(’) S(’). \[eq:S\_bend\_merge\]

Perturbative calculation of the scattering matrix for narrow bent wave-guide {#sec:pert_scatt}
============================================================================

We present a semi-classical approximation, for $k\gg 1$, of a scattering matrix corresponding to a single bend on a straight wave-guide of width $a$, as one shown in figure \[pic:schema\]. Here we are discussing only narrow channels $a \ll 1$, where the influence of closed modes on the scattering diminishes. Therefore closed modes are neglected in our calculations. We are working at wavenumbers $k$, where in all regions of the open-billiard the number of open modes is equal. This enables us to write the reflection and the transmission matrix in the following simpler form R = -\[C\_+ - S\_+\]\^[-1]{} (C\_- + S\_-), T = 2 \[C\_+ - S\_+\]\^[-1]{}, \[eq:semiclaas\_scatt\] where we use the diagonal matrices ${\cal C} = \Re \{{\cal F}\}$, ${\cal S}= \Im \{{\cal F}\}$ and $G =  \diag \{g_n\}_{n=1}^{\rNo}$ to express the introduced matrices C\_ &=& G\^[12]{} A [C]{} B\^T G\^[-[12]{}]{} G\^[-[12]{}]{} B [C]{} A\^T G\^[12]{}, \[eq:semiclaas\_C\]\
S\_ &=& G\^[12]{} A [S]{} A\^T G\^[1 2]{} G\^[-[12]{}]{} B[S]{} B\^T G\^[-[12]{}]{}. \[eq:semiclaas\_S\] We proceed by rescaling the variables to dimensionless form by the following substitutions y = a , r = q + a, = a k, h\_n = a g\_n =, \_p = v\_p, \[eq:semiclass\_substitutions\] with a new transverse coordinate $\xi\in[0,1]$, and geometric properties being described by the parameter $\alpha = a/q \ll 1$. The transition matrices are then expressed as A\_[np]{} = q\^[12]{} Q\_[np]{}, Q\_[np]{} = , \[eq:trans\_mat\_asym\_Q\]\
B\_[np]{} = q\^[-[12]{}]{} P\_[np]{}, P\_[np]{} = . \[eq:trans\_mat\_asym\_P\] with $b_n(\xi)=\sqrt{2} \sin(\pi n \xi)$. The eigen-pairs $(v_p,\phi_p(\xi))$ are defined by the following differential equation and the boundary condition: + (\^2 - )\_p = 0 , \_p(0)= \_p(1) = 0. \[eq:mode\_equation\_dimless\] We can easily recognize that the solutions of equation (\[eq:mode\_equation\_dimless\]) converge in the limit $\alpha\to 0$ to $v_p = h_p$ and $\phi_p(\xi) = b_p(\xi)$. We assume that the solutions can be expanded in a power series of variable $\alpha$. The eigen-pairs can then be obtained using the standard perturbation theory with the perturbation parameter $\alpha$. The rescaled mode numbers are written as v\_p\^2 = h\_p\^2 + \^2 ( + O(h\_p\^4) ) + O(\^3), \[eq:mode\_num\_asym\] and the rescaled mode functions read as \_[p]{}() = \_n V\_[np]{} b\_n (), V\_[np]{} = \_[np]{} + \_[n+p]{}\^[odd]{} + O(\^2), \[eq:mode\_fun\_asym\] where we use the symbol $\delta_n^{\rm odd} = (1: n \textrm{ is odd};\;0 : \textrm{otherwise})$. By plugging the mode functions $\phi_p (\xi)$ (\[eq:mode\_num\_asym\]) into transition matrices $Q$ (\[eq:trans\_mat\_asym\_Q\]) and $P$ (\[eq:trans\_mat\_asym\_P\]) we obtain Q\_[np]{} = \_[np]{} + (F\_[np]{} + \_[np]{}) + O(\^2), \[eq:trans\_mat\_asym\_pert\_Q\]\
P\_[np]{} = \_[np]{} - (F\_[pn]{} + \_[np]{}) + O(\^2), \[eq:trans\_mat\_asym\_pert\_P\]\
F\_[np]{} = ( - ) \_[n+p]{}\^[odd]{}. \[eq:trans\_mat\_asym\_pert\_F\] Note that the rescaled transition matrices $Q$ and $P$ satisfy the known identity $Q P^T = P Q^T = \id$. We insert the expressions for $Q$ (\[eq:trans\_mat\_asym\_pert\_Q\]) and $P$ (\[eq:trans\_mat\_asym\_pert\_P\]) back into $S_\pm$ (\[eq:semiclaas\_S\]) and $C_\pm$ (\[eq:semiclaas\_C\]) and write the reflection and the transmission matrix as R = - + O(\^2), \[eq:R\_asym\_pert\]\
T = [F]{} - (H\^\[[F]{},F\]H\^[-]{} + H\^[-]{}\[F,[F]{}\]H\^) + O(\^2). \[eq:T\_asym\_pert\] where we have introduced the symbol $(A)^{\rm s} = A + A^T$ and the diagonal matrix $H=\diag\{h_n\}_{n=1}^{\rNo}$. The approximations of the reflection matrix $R$ (\[eq:R\_asym\_pert\]) and the transmission matrix $T$ (\[eq:T\_asym\_pert\]) are valid far away from the resonant condition $ak = \pi n$, because we assumed that $|g_m| > \alpha $ for all $m \le \rNo$. We conclude that the strength of reflection scales as ${\cal R}\sim \alpha^2$ and that narrow channels can be treated as perturbed straight wave-guides.
