Expanded boundary integral method and chaotic time-reversal doublets in quantum billiards

We present the expanded boundary integral method for solving the planar Helmholtz problem, which combines the ideas of the boundary integral method and the scaling method and is applicable to arbitrary shapes. We apply the method to a chaotic billiard with unidirectional transport, where we demonstrate the existence of doublets of chaotic eigenstates, which are quasi-degenerate due to time-reversal symmetry, and a very particular level spacing distribution that attains a chaotic Shnirelman peak at short energy ranges and exhibits Gaussian Unitary Ensemble (GUE) like statistics for large energy ranges. We show that, as a consequence of such particular level statistics or algebraic tunnelling between disjoint chaotic components connected by time-reversal operation, the system exhibits quantum current reversals.


Introduction
The solution of the planar Dirichlet problem for the Helmholtz equation has served as one of the principal numerical set-ups for verifying and demonstrating the main ideas of quantum chaos [1]- [7]. The reason is that this problem, usually referred to as the quantum billiard, is one of the simplest quantum Hamiltonian stationary problems for which good quality spectra of highly excited eigenstates can be obtained for various, even non-integrable, geometries.
In the literature there exists a good account of numerical techniques to tackle the problem [8]. Apart from many ingenious ideas developed for specific geometric set-ups, there are two competitive general purpose approaches: (i) the boundary integral method (reviewed in [9]) or (ii) Heller's plane wave decomposition method [2]. While the boundary integral method is really a general rigorously based method, the plane wave decomposition method is a rather heuristic approach which can fail in certain important general cases, e.g., of non-convex geometries [10]. Still, there exists an extremely efficient implementation of plane wave decomposition method due to Vergini and Saraceno [11], which makes this method an attractive option in spite of potential mathematical problems. The crucial new idea of Vergini and Saraceno was to look for a minimum of a boundary norm (an integrated norm of the wavefunction or the field amplitude along the boundary of a two-dimensional (2D) domain) in terms of a solution of a generalized eigenvalue problem, the gain being that in a single computational step a number of accurate eigenvalues are obtained in a constant proportion to the dimension of the matrices scaling as O(k) where k = 2π/λ is a referential wavenumber, in contrast to the boundary integral method, where a number of matrix computations have to be performed in order to locate a single eigenvalue. In both cases each matrix computation is of the order O(k 3 ).
In this paper we propose to use a similar (scaling) idea in conjunction with the rigorous boundary integral method. We develop a completely general numerical boundary integral technique which produces a constant fraction of k eigenvalues per single matrix operation involving O(k 3 ) scalar operations, thus being asymptotically orders of magnitude faster than traditional implementations of boundary integral method.
The details of the method, which comprises the first part of the paper, are elaborated in section 2.
The second main idea of the paper is to apply our method in a chaotic quantum billiard whose phase space structure is not time-reversal invariant, in a sense that its chaotic phase space components are not mapped onto themselves upon the time-reversal operation. 4 Namely we study 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT the billiard in the domain of non-convex and non-simply connected shape, the so-called Monza billiard. Our choice of the model system also represents a good benchmark for numerical methods since it is likely one of the most difficult types of the billiard shape that one can think of. The corresponding classical billiard-being a curved closed corridor with parallel walls-possesses the property of unidirectional motion, namely the classical phase space separates into two disjoint ergodic and chaotic components of clockwise and counter-clockwise motions. The absence of any geometric (point) symmetries of the model, and the fact that the time-reversal symmetry does not preserve the invariant ergodic components, has an interesting consequence, namely the existence of algebraically (in effective Planck constanth) split quasi-degenerate pairs of energy levels. The effective Planck's constant of the system is defined as the ratio of the physical Planck's constant to a typical action of a system. In billiard systems, the effective Planck's constanth ∝ 1/k where k is the wavenumber. The existence of quasi-degenerate pairs has a dramatic effect on energy level statistics, in particular since the level splitting is of the same order as the mean level spacing. The effect can be interpreted as a chaotic analogue of the Shnirelman peak [12,13] known for a long time for systems with non-time reversible KAM islands.
Another and perhaps even more dramatic effect concerns the long range correlations among levels which behave according to the Gaussian Unitary Ensemble (GUE) of random matrices in spite of the fact that the system as a whole possesses time reversal symmetry. This is explained by the fact that even though the eigenfunctions of the Hamiltonian are strictly real, the doublets of functions exchanged by the time-reversal operation, namely Tψ L,R = ψ R,L , which become preserved in the classical limit even though they are not exactly eigenfunctions of the Hamiltonian, are complex and hence GUE should be used for modelling spectral statistics at large energy ranges. 5 This is a new effect, which to our knowledge has not been predicted or known before.
Detailed discussion of the Monza billiard, its spectral statistics and dynamics, is given in section 3. In section 4, we discuss our main result and conclude.

Expanded boundary integral method
The scaling method as proposed by Vergini and Saraceno [11] allows one to compute a number of states of a quantum billiard, or any other system described by the Helmholtz equation with the Dirichlet boundary condition, lying close to a chosen reference value of the wavenumber k without losing any state in a chosen interval. One of its main advantages is that the matrices to be diagonalized are not of the dimension of the order of ∝ k 2 , as would be the case when using the usual diagonalization techniques, but only of the order ∝ k. While this fact makes it one of the most efficient approaches to solve for eigenvalues and eigenstates of quantum billiard problems, especially for high values of k, experience shows that its domain of applicability is rather limited. As shown by Gutkin in [10], the plane wave decomposition method [2] on which the scaling method is typically based can in general only be applied to a convex billiard. The scaling method does not necessarily involve a plane wave decomposition and a different basis set can be used, which enables the method to solve problems inaccessible by the plane wave decomposition. One such example is the use of Bessel functions of fractional order to solve the mushroom shaped billiard problem [15]. The choice of the basis in such cases must, however, be specifically tailored to the system. The boundary integral method, (see [9] for a review of the method), has a similar demand on matrix size as the scaling method and can in principle be applied to arbitrary shapes as it is based on the Green's theorem that shows how the wavefunction inside the billiard can be represented with its normal derivative on the boundary (in the case of Dirichlet boundary conditions). From this representation the integral equation for the normal derivative at the boundary u = ∂ ∂n ψ follows as [9] u(s) = −2 where ∂ denotes the boundary of the domain, s is the arc length parametrization of the boundary. The coordinates of the boundary are given as q(s) and ∂ ∂n G k (q, q ) is the inward normal derivative of the Green's function for the unbound problem at wavenumber k. Following reference [9] we may discretize equation (1) on a finite set of points as determined by their arc length parameters s i , yielding where u i is the value of the normal derivative at the point q(s i ), τ ij = |q(s i ) − q(s j )| is the distance between points q(s i ) and q(s j ), cos φ ij = n(s i ) · (q(s i ) − q(s j ))/τ ij gives the angle between the inward boundary normal n at the point q(s i ) and the distance vector between the two points and κ i is the boundary curvature at the point q(s i ), where positive curvature is taken for concave parts of the boundary. The arc-length of the boundary section centred at q(s i ) is denoted by l i . The discretization is characterized by the parameter b = 2πN p /(kL b ), where N p is the number of discretization points and L b is the total length of the billiard boundary. The parameter b measures the number of discretization points per wavelength. Equation (2) only has solutions for those values of k where the matrix A becomes singular. While in principle any Green's function could be used, the choice of the complex Hankel function of the first kind is necessary in order to avoid spurious solutions that occur because equation (1) only then becomes a sufficient condition for determining the normal derivative. The usual approach towards solving this equation is based on solving for zeros of the Fredholm determinant. It therefore yields not more than a single level (if at all) per determinant calculation and is typically plagued with loss of levels and accuracy when close to a degeneracy. The latter problem has been resolved by Bäcker [9] by calculating not the determinant directly but the singular value decomposition (SVD) of the matrix A and following the behaviour of individual singular values.
At this point, we introduce a similar procedure to the one presented for the scaling method, namely trying to determine the solution close to a chosen reference value of the wavenumber k 0 by varying k = k 0 + δk. We may write the equation (2) as a Taylor expansion where A (k 0 ) and A (k 0 ) are the first and the second derivative of the matrix A with respect to k at the point k 0 , obtained by taking the derivatives of each matrix element. Let us expand δk = δk 0 + 2 δk 1 + . . . and u = u 0 + u 1 + 2 u 2 + . . . in terms of the order of the formal perturbation parameter , whose powers give the order of the terms with respect to the small parameter δk 0 and should be set to one in the final result. This also shows that the procedure can only be applied to those eigenvectors u whose wavenumbers lie close to the reference value k 0 . By inserting these expressions into the equation (4) and ordering the terms with respect to we obtain (dropping the implied argument k 0 ) We will now show that taking the first two terms of the above equation gives a consistent starting point for solving the perturbative problem. We may do so as there is a degree of freedom involved as to how to define the terms u 0 and u 1 in the perturbation series, and the equation (6) therefore defines u 0 . 6 This equation also shows that, while A and u 0 are themselves of the zero-order in the expansion, their product is of the order O( ). This is due to the fact that we are only slightly shifted in terms of wavenumber k from the exact solution of the equation (2), where this product is exactly zero.
Looking again at all the terms of the equation (5) with up to the linear order in and using equation (6), we obtain This statement implies that the first-order corrections to the eigenvector should belong to the subspace of the matrix A for which Au 1 = O( ). Observing the previous argument regarding the product of A and u 0 , this subspace is the space of those solutions of our problem that reside close to the reference value k 0 . The first-order correction to an eigenstate will therefore mostly contain contributions from the other states in its vicinity. We may improve the accuracy of the levels by multiplying the equation (5) with the adjoint of the left eigenvector v 0 , defined as Due to the equations (6) and (8), the first four terms are exactly eliminated in pairs. Although A and v 0 are both of the zero order in terms of , the product v † 0 A = O( ) is of a higher order, as can be seen from equation (8). This means that the last term in the equation (5) when multiplied by the left eigenvector also becomes of the order O( 3 ). Only two terms remain, giving and therefore k = k 0 + δk 0 By calculating the equations (6) and (9) for various values of k 0 , the chosen spacing between the values depending on the desired accuracy, we may obtain all levels of a system within some 6 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT chosen interval. As the average density of states varies approximately as ∝ k 0 , by taking the above error estimate into account, the spacing between various reference wavenumbers k 0 at which the individual calculations are performed must vary as k 0 ∝ k if the error of calculating an individual level is to remain a constant percentage of the mean level spacing. The method itself is applicable even for calculating the ground state of the system.
For simply connected domains the above procedure is found to be robust with respect to the shape of the boundary, even in the case of arbitrary corners. The domains with holes require an additional step as spurious solutions are found in the spectrum. These solutions appear because the transpose of the matrix in equation (2), which is here used to solve the interior Dirichlet problem, gives the solution to the exterior Neumann problem as well. As the holes of the domain represent the exterior of our problem and are themselves compact, the left eigenvectors of our spurious solutions are simply the boundary values of the wavefunctions of the Neumann problem within the holes. We may therefore eliminate the spurious solutions by solving the Neumann problem for each of the holes by performing the same computation for each hole, where the boundary taken is now just the boundary of the hole in question, and then discard the levels obtained for each hole from the spectrum of the total problem.

Monza billiard
The method was applied to a billiard that was named the Monza billiard due to its similarity with the famous Italian racetrack. It comes from a family of unidirectional billiards, as defined in [16]. These billiards are shaped as channels with parallel walls and have the property that the classical trajectories going in one direction along the channel cannot reverse this direction. The phase space of such systems is therefore split into two disjoint regions, and the motion within each separate region can be fully chaotic. The two regions are separated by a 1D family of marginally stable bouncing ball orbits. The shape of the system as shown in figure 1 is one of the simplest nontrivial closed shapes without any geometric symmetry that belongs to this class of systems. For all the subsequent examples shown, the system parameters as defined in the figure were chosen to be q = 1/12, a = 1/2, b = 1/3, r = 1/3, α = 1. The billiard was numerically tested to be classically fully chaotic and ergodic within each of the two invariant phase space components. 7 We used the expanded boundary integral method to compute all levels of the Monza billiard up to the wavenumber k = 70. We used the Weyl formula [17] to test that we did not lose any levels in the chosen interval. The discretization parameter was chosen to be b = 12. The spacing of the reference wavenumber k 0 between individual runs of the method was chosen to be k 0 = 0.05k The principle of uniform semiclassical condensation (PUSC, [7,18,19]) states that in the semiclassical limit of effectiveh → 0 the individual eigenstates of a system correspond to invariant objects in the classical phase space. One would perhaps naively expect the eigenstates of a Monza billiard to correspond to either one of the two chaotic components. This, however, cannot be the case since states corresponding to either chaotic component would necessarily have a nonzero probability current [14]. The system, on the other hand, possesses the time reversal symmetry and therefore its non-degenerate eigenfunctions must be fully real and as such can have no probability current. This seeming contradiction is resolved by noting that in the Monza billiard most eigenenergies are to be found in nearly degenerate doublets. Only the states corresponding to the bouncing ball modes are singlets, with their relative measure vanishing in the semiclassical limit, as they correspond to a classical invariant separatrix between the two chaotic components which itself is the family of bouncing ball orbits corresponding to the particle bouncing perpendicularly between the two walls and is of measure zero in the classical phase space. The doublets are a chaotic manifestation of the Shnirelman peak in the level spacing distribution [12] that is generally found in all systems with a time reversal symmetry but having no point symmetries. Typically in such systems, the Shnirelman peak is due to the states that correspond to the classical invariant tori which do not map back on to themselves upon the operation of time reversal. There is usually exponentially small tunnelling (in terms of effectiveh) between the tori that are connected via the time reversal operation and this is reflected in an exponentially small energy splitting between the states that correspond to the torus pair. In the Monza billiard, however, the tunnelling does not occur between tori but between two chaotic components that are separated by a bouncing ball manifold of measure zero in phase space. This means that the expected tunnelling effects will not be exponentially small but will rather scale as a power law in terms of effectiveh as the two chaotic components are a distance zero apart. In fact, it appears that the average splitting of the pair remains constant as a fraction of the mean level spacing. Heuristically, this can be explained by noting that the tunnelling amplitude between the states localized in classically invariant chaotic phase space components can be estimated in terms of the phase space overlap of the corresponding Wigner functions which can be semiclassically estimated to scale as a linear function of an effective Planck constant.

8
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT To demonstrate this, let us fix an arbitrary point in the billiard q and locally represent the wavefunction as a combination of random plane waves where p ϕ = p 0 (cos (ϕ + β(q)), sin (ϕ + β(q))). Due to the unidirectionality, only the waves in one half of the wavevector space are taken, where β(q) denotes the polar angle of the bouncing ball trajectory going through the position q. We assume a distribution of the stochastic variable f such that with ρ = 1 πA where A is the area of the billiard and is the Heaviside step function. Inserting the expression (11) into the definition of the Wigner function we obtain We define the 2D domain The expected value of such a Wigner function is then Performing the integration with respect to x the exponential turns into a wide delta function such that The typical width σ of the delta function is nonzero due to the finite size of the system. Its value is σ ∝h/ where is some linear-dimension of the system independent ofh. To get an estimate of the overlap in the equation (10) we also need to introduce the Wigner function of the state travelling in the opposite direction 9 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT with the integration boundaries complementing those in the expression (18). If we neglect the correlations we may use the expected values for the Wigner functions in the equation (10). The expected values of the two Wigner functions then overlap due to the width σ, and the overlap itself is proportional to that width, thus demonstrating the statement (10). 8 The spectrum can therefore be thought of as a composition of two nearly identical sequences, where their difference remains significant on the scale of the mean level spacing. It is important to stress that this near degeneracy cannot be removed by any means such as desymmetrizing of the billiard as the system chosen does not possess any point symmetries. From this, we can deduce that the unidirectional family of billiards will generally exhibit a non-universal spectral statistics. In the level spacing distribution P(S ), which is the distribution of spacings S between consecutive levels in an unfolded spectrum where the mean level spacing is equal to one, we expect two main contributions, namely one from the spacings within individual pairs and another from the spacings between consecutive pairs. The first contribution is a widened delta distribution with the weight 1/2 close to the spacings S = 0, and the other contribution would be expected to be a stretched Gaussian Orthogonal Ensemble (GOE) [20] contribution. The widening of the second contribution by a factor of two is due to the unfolding procedure which requires that the mean level spacing is equal to one. In a sequence of pairs, the average spacing between centres of different pairs is therefore twice the mean level spacing.
The cumulative level spacing distribution for the system with the parameters q = 1/12, a = 1/2, b = 1/3, r = 1/3, α = 1 is given in figure 2. The number of levels used in the figure was N = 2403. For comparison, the same statistics using only the last 500 levels is shown (dashed red) as well in order to demonstrate that the level spacing distribution indeed does not change with energy (orh). Apart from the initial contribution that is mostly due to the spacings within individual pairs, the level spacing distribution then follows the stretched GOE prediction where is the Wigner surmise [20]. This prediction holds well for small and intermediate spacings, S < 2 (shown dotted green), yet at larger spacings the distribution starts to deviate from the GOE prediction but rather approaches the (stretched) GUE prediction (dot-dashed blue). This can be obtained from equation (21) but replacing P GOE with While not exactly following this prediction, a clear trend does emerge for the W(S) to move from the GOE to the GUE behaviour at larger spacings. The longer range spectral statistics are investigated by using the spectral rigidity [20], which is defined as where N(E) is the spectral staircase function whose value increases by one at each (unfolded) energy level. The rigidity measures the average deviations of the spectral staircase function from the best fitting straight line of length L. It is given for the spectrum of the Monza billiard in the figure 3. We are comparing the numerical results to the GO(U)E predictions for the random matrix ensembles as defined in [20]. Since most levels come in pairs, the predictions need to be scaled as˜ The stretching in L is necessary because the average level spacing between level pairs is twice the average level spacing. The factor of four occurs since each step of the spectral staircase function is twice as large for each pair as it would be for an individual level, and the 3 statistics is quadratic in spectral staircase fluctuations. The results for the 3 statistics of the Monza billiard are given in figure 3. At low values of L, the results of both predictions as well as the numerical results agree. The numerics undershoots both predictions because the pairs are not exactly degenerate, with the small spacings within the pairs reducing the fluctuations of the spectral staircase function in comparison to the assumed exact degeneracy in the GOE and GUE predictions. At larger distances L, the numerics clearly follows the GUE rather than the GOE prediction. At distances larger than about L = 10 the numerical 3 statistics starts to increase faster than any of the predictions. This effect could be expected because of the inherent presence of the large bouncing ball mode contribution to the spectrum [21]. Excluding this effect, however, the behaviour of the spectrum at larger distances appears to follow the GUE prediction.
To understand such behaviour, we need to look into the properties of state pairs. In figure 4, we show two states corresponding to a neighbouring pair of states ψ a and ψ b . Both states are of course fully real, with the phase switching between zero (shown red) and π (shown cyan), but a linear combination ψ L,R = (ψ a ± iψ b )/ √ 2 of two real states can be shown to have the largest current amplitude. In our case, the direction of the phase change, which itself corresponds to the direction of the current probability, shows that this current flow is clearly seen to be predominantly in a single direction. The complex conjugate of the function is orthogonal to the function itself and corresponds to a current in the opposite direction. The current behaviour is even more clearly exposed by looking at the Husimi plots of the 1D wavefunction cross-section along the coordinate z taken at w = 1/2 (see figure 1 for coordinate definitions). While for the two eigenstates the Husimi plots are nearly identical, examining the Husimi plot for the complex combination of the two states clearly yields a state that corresponds to a current in a single direction. This is typical for most pairs of states and, while PUSC cannot be strictly applied to individual eigenstates, taking a nearly degenerate pair and performing the above complex rotation in this subspace yields states that support the PUSC conjecture. 9 The nature of the spectrum naturally poses the question of the dynamics of such a quantum system. Classically, a particle traversing the billiard in one direction will keep that direction permanently. Starting with a quantum state such as shown in the bottom part of the figure 4, as this state is not an eigenstate of the system, it will exactly reverse its character every t = πh/ E where E is the energy splitting of the pair of the two eigenstates forming the initial state. Let us now define the operator for the current along the billiard direction aŝ 9 A similar analysis applies to (almost) degenerate tori in a KAM-like scenario.  figure 1) are shown to the right of each corresponding wavefunction. The darkness of the wavefunction corresponds to the probability density, whereas the hue going circularly from red to green to blue to red again represents the change of the phase of the wavefunction. In the Husimi plots there are 10 contours of constant value shown, starting from zero to the maximum value, where on the x-axis the coordinate that runs along the length of the Monza billiard z is given, while on the y-axis the momentum along this direction is shown in the units of the maximum attainable momentum at the corresponding classical energy. where e z is the unit vector in the direction of the z coordinate in the billiard (see figure 1). The initial state is chosen to be a Gaussian wavepacket with p 0 = (10, 0), r 0 = (0, 1/2) and σ = 1/5 with the coordinate frame as shown in figure 4. We expand this initial state in the basis of all eigenstates up to the wave-vector k = 20, where these states were obtained using the boundary discretization parameter b = 18. The evolution equation is taken to be the dimensionless Schrödinger equation The expected value of the current operator J as a function of time is shown in figure 5. As can be seen, for a very short time the initial wavepacket has a large current that up to the time t = 1 relaxes to the value v 0 2/π (shown dashed green), where v 0 = 20 is the classically expected velocity for the wavepacket with the chosen average energy. Classically, the current is expected to remain at this plateau, but the quantum system experiences a decrease of current that actually turns to fluctuating reversals of the current as it does not stabilize at zero. This behaviour remains even for times much larger than the ones shown in the figure almost uniformly covers the billiard and therefore corresponds to the wavepacket as spread over the whole unidirectional chaotic component. The phase of the wavefunction clearly indicates that the direction of the quantum current is indeed directed in a single direction along the billiard for all the times shown. In figure 8, we show a still from the video representing the evolution of the wavefunction at the time t = 175 until the time t = 175.3, where the current as shown in the figure 5 experiences its strongest reversal. From the phase representation the direction of the current is not easily read. Figure 5, however, shows that the global unidirectional current contribution is quite sizable but the local fluctuations seem to obscure it.

Conclusions
The goal of the present paper was twofold.

Presentation of a novel efficient method for numerical solution of the 2D Helmholtz equation
with Dirichlet boundary condition (e.g. a quantum billiard) which is obtained by combining rigorously founded boundary integral method with the ideas of Vergini and Saraceno [11] of expanding the boundary norm as a function of energy. 2. Application of the new technique to compute the energy spectra, their statistical properties, and time evolution in a quantum billiard with a chaotic classical limit whose phase space structure is not time reversal invariant. Here we mean that the system itself is (globally) timereversal invariant but different disjoint chaotic components are not mapped onto themselves upon time reversal operation. Such situation is possible in the so-called Monza billiard due to unidirectionality of the classical motion [16]. We show the presence of energy level statistics converging to a well defined distribution in the semiclassical limit, which is non-universal and exhibits a Shirelman-type peak at small energy ranges and GUE fluctuations at large energy ranges. We have also computed the time evolution of the expectation value of the tangential momentum (or particle current) which exhibits interesting irregular oscillations due to algebraic quantum tunnelling phenomena. The time-scale of the current reversal and consequent oscillations is a constant factor (typically much larger than unity) of the Heisenberg time. The results have been qualitatively explained in terms of a heuristic semiclassical picture, namely the PUSC combined with the nontrivial effects due to time reversal symmetry.
The phenomenon of GUE-like fluctuations in a time-reversal invariant system, which we discovered at large energy ranges, is similar to the behaviour observed by Leyvraz et al [22] in systems with point symmetries lacking real representations. Indeed the two phenomena can be understood to have some formal similarities. In future explorations it would be interesting also to consider systems which live in both situations simultaneously, having classically unidirectional motion and point symmetry with complex representations.
As the quantum billiard models can nowadays be easily realized in various experimental set-ups (e.g. quantum dots, microwave cavities, optics etc), we expect that predicted effects should be experimentally observable.
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT function ρ(x) = 4x(1 − x). Firstly, we calculate which gives the estimate of the (weighted) number of levels in the interval (k 0 , k 0 + k 0 ). If this weighed estimate is smaller than, say 1/2, we know there are no levels in the neighbourhood of the middle point of the interval and we may just merge the two subsequences by trimming them in the middle and joining them. Otherwise we trim the first sequence by discarding all the levels that are higher than the middle point of the interval and then vary the trimming point k c for the second sequence in such a way that differs as little as possible from σ. One should also always demand that k c > k 0 (preferably even k c > k 0 + ( k 0 /4) to avoid small contributions to σ from the edges of the interval), otherwise possible levels outside the chosen interval may be inadvertently added.