Modeling of open quantum dots and wave billiards using imaginary potentials for the source and the sink

A heuristic model for particle states and current flow in open ballistic two-dimensional (2D) quantum dots/wave billiards is proposed. The model makes use of complex potentials first introduced in phenomenological nuclear inelastic scattering theory (the optical model). Here we assume that external input and output leads connecting the system to the source and the drain regions may be represented by complex potentials. In this way, a current may be set up between the two ‘pseudo-leads’. Probability densities and current flows for an open quantum dot are analyzed here numerically and the results are compared with the microwave measurements used to emulate the system. The model is of conceptual as well as practical interest. In addition to quantum billiards, it may be used as a tool per se to analyze transport in classical wave analogues, such as microwave resonators, acoustic resonators, effects of leakage on such systems, etc.


Introduction
There are a number of techniques to fabricate various planar enclosures in which wave motion becomes quantized, i.e. the wave length is of the order of the dimensions of the confinement. Typical examples are semiconductor quantum dots, which are currently of great scientific and technological interest (see, for example, [1]). Figure 1 shows an ultra-small rectangular quantum dot fabricated from a layered GaAs/AlGaAs semiconductor heterostructure material with a patterned metallic gate [2]. The dot in the middle region is connected to the right and left reservoirs by narrow leads ('quantum wires' or alternatively 'quantum point contacts'). The electrically active region is the GaAs/AlGaAs interface underneath the patterned top gate [1].
In another celebrated example, one makes use of iron atoms placed on the (111) surface of copper [3]. By means of a scanning tunneling microscope (STM), these atoms may be moved around on the surface and arranged into various closed geometric patterns. Copper surface states are then trapped inside such 'corrals' and the corresponding probability distributions may be observed by STM techniques. For a square enclosure, for example, one recognizes the standing waves of a particle in a hard-wall box, for the circular corral there are Bessel functions and so on. There are also quantum billiards with ultracold atoms in arbitrarily shaped two-dimensional (2D) boundaries achieved by optical dipole potentials. This has allowed the study of various chaotic and regular systems in quite flexible ways [4]. Quantization may not only be achieved on the micro level. For example, macroscopic planar electromagnetic resonators display similar phenomena. In a number of cases, they may serve as classical analogues onto which quantum systems may be mapped [5].
There is a great deal of interest in modeling systems as above. The main purpose of the present paper is to elaborate an elementary heuristic model for electron states in quantum dots (alt. 'cavities', 'billiards', 'corrals', 'resonators'), and to show how transport through such systems may be modeled by introducing complex potentials to simulate the source and the sink for the carriers. Our model originates from early phenomenological nuclear scattering theory [6,7] in which absorption of scattered particles is accounted for by a fictitious imaginary scattering potential. In our case, such a potential can act as a sink. Here, we extend this model by introducing an additional imaginary potential of opposite sign. In regions where this one is nonzero, it acts, not unexpectedly, as a source of particles. Having both the source and the sink Scanning electron micrograph of a gated 2D square quantum dot with dimensions 1 × 1 µm 2 (middle dark region). The dot is connected to the surrounding electron reservoirs (outer dark areas) by two openings (leads), through which an electric current may be passed from one reservoir to the other by applying a driving voltage between the reservoirs. Regions in light gray are patterned metallic gates under which the electron gas is depleted. (Adapted from [2], courtesy J P Bird and with permission from IOP.) embedded in a system, we may thus arrange it so that particles flow from one lead to the other, i.e. we may simulate quantum transport through, for example, a quantum dot. Here, we show how the model works in practice by applying it to the open 2D quantum dot in figure 1 and how this case may be validated by a mapping onto a microwave resonator of the same shape. The model has been addressed previously for studies of quantum transport in single molecules and nanotubes [8,9]. Here, our focus is, however, quite different. We will rather highlight the delicate features of probability density distributions, current flows, the formation of vortices, etc. We also propose that the model per se should be a useful tool for analyzing related features in analogous classical wave systems, such as microwave and acoustic resonators, effects of leakage and probe antennas, etc.
The paper is organized in the following way. In section 2, we introduce our model and show how the generalized complex potential may be used in the single-particle Schrödinger equation. In particular, we discuss the physical interpretation of the imaginary part of the complex potential. Numerical methods for closed as well as open 2D cavities are outlined in section 3, and explicit results for probability densities and current flow are given in the same section for the device in figure 1. In section 4, the similarity between our 2D quantum system and microwave billiards [5] is discussed and a comparison with measurements is presented. A summary and some concluding remarks are given in section 5.

Quantum mechanics with complex potentials
Consider first the usual time-dependent Schrödinger equation for a single particle where V is the potential, which is real, and m is the mass of the particle. The continuity equation is where (r, t) = | (r, t)| 2 is the probability density and J is the probability current density For stationary states, equation (1) turns into its time-independent form in the usual way where E is the energy of the particle. For clarity, we introduce the notations (r, t) → ψ(r)exp(−iEt/h), J(r, t) → j(r) and (r, t) → ρ(r) for stationary states. According to equation (2) there is a differential conservation relation between (r, t) and current J(r, t) when V is real. Consider now the general case for which V may be complex. As mentioned, the idea of using a complex potential was introduced long ago in a phenomenological model for absorption in nuclear scattering [6,7], i.e. in situations when particle conservation is not mandatory. The model with complex V , also referred to as the optical model, is indeed a heuristic one and is, as such, frequently used in nuclear physics as well as in other fields, such as atomic and molecular physics, nanotransport, etc to describe inelastic losses. We specifically mention the modeling of electron wave patterns in quantum corrals and the reduced reflectivity of the surrounding walls of iron atoms [3] and studies of electron transport of quantal trajectories and dephasing in 2D nano-devices [10]- [12]. Let the complex potential be where V R and V I are real. By taking the time derivative of (r, t) and using the Schrödinger equation, we arrive at the generalized version of the continuity relation Since (r, t) is positive, the term on the right-hand side of this equation acts as the source or the sink for positive or negative values of V I , respectively. This is seen more clearly if we integrate equation (6) over a volume bounded by the surface A. Following [9], for example, we then have where n is the unit vector normal to the surface. If the system is well contained in , the surface term vanishes, implying that the normalization integral decreases with time if V I is negative, i.e. particles are being absorbed. For V I > 0, on the other hand, the normalization increases as particles are created, in this case [7]. In general, the eigenvalues are now complex and there is either an exponential growth or decay, depending on the sign of V I .
To study spatial distributions of probability currents and densities, vortices, etc in 2D quantum dots and billiards, one may combine the two cases, as in [8,9], by letting V I have Modeling of a 'two leads' rectangular cavity for which the leads are replaced by stubs. Inside the cavity (including the stubs) the real part V R = 0; in the exterior shaded areas, V R should be given a value that is much higher than the levels that one wishes to compute. The dark areas at the two ends of the stubs indicate the location for the source and the sink (±iV I ).
different signs in different parts of the system. For example, to simulate the probability current density flowing between input and output leads in a cavity, as in figure 1, we may simply replace the extended leads by two regions in which V I is nonzero and of opposite signs. The condition for a steady, well-balanced and lossless flow between the two regions is Therefore, the expectation value of the Hamiltonian operator, H , is real in spite of the potential V = V R + iV I being complex. Furthermore, the corresponding eigenvalues are real when V R is independent of time. As we will see in the following sections, our elementary approach of using ±iV I for the source and the sink to emulate probability currents, vortices, etc within a 2D cavity performs remarkably well. We refer to this case as an open 2D cavity.

Particles and currents in an open two-dimensional quantum dot
As an application of the present model, consider the open quantum dot in figure 1. We choose this particular system because of previous joint experimental/modeling work [2] and because there are matching microwave measurements [13]- [17] (for a recent overview, see, for example, [17]). Rather than dealing with the entire system consisting of the cavity itself, the two leads and the reservoirs, we do away with all of that by introducing two stubs, as in figure 2, with imaginary potentials of opposite signs at the two ends of the stubs. Hence, by this step we emulate both emission and absorption of particles and a current flow between the 'input' and 'output' leads. The effective confinement in figure 2 is smoother as compared to the lithographic pattern in figure 1. This is typical of the electrostatic confinement at the active interface region in 6 which the actual dot is embedded [2]. In the following, we assume hard confining walls and a flat potential V R = 0 in the interior region of the dot. When also V I = 0, the Schrödinger equation reads where k = 2m E/h 2 is the usual wave number for real values of E.

Numerical solutions: the finite difference method
To solve the Schrödinger equation numerically, we use the finite difference method (FDM) for the states of interest [18]. Here we prefer to use FDM as it is robust and straightforward to implement and, once a wave function is known, it is easy to evaluate related quantities like the probability current density, as well as other features, like statistical distributions and correlations. FDM may also be extended into the time domain (finite-difference time domain, FDTD) for studies of time-dependent processes. The present geometries are, however, quite simple. For more complex structures and/or much higher states, other computational methods could be used, e.g. FEM, Green's functions and surface integral methods.
In FDM, one normally generates a computational square grid (i, j) (row, column) of size N × M with grid points (i, j) covering the system. Denote the value of the wave function ψ(x, y) in grid point (i, j) by ψ i, j . Usually only neighboring grid points are considered. We may therefore approximate first-and second-order derivatives of ψ(x, y) in grid point and similarly for the y-coordinate: a is the equidistance between grid points and (i, j) are integers. The discretized stationary Schrödinger equation is then where V i, j is the value of the potential V (x, y) at grid point (i, j) and E is the energy eigenvalue (which may be complex when the potential is complex). The elements ψ i, j may be collected into a column vector C = (ψ 1,1 , ψ 1,2 , ψ 1,3 , . . . , ψ 2,1 , ψ 2,2 , . . .) and we may therefore rewrite equation (12) as a matrix equation where H is tridiagonal. The eigenvalue equation (13) may be solved by standard routines.
Once the solutions C are known, the probability current density for the different states may be evaluated using equations (3) and (10). To solve this eigenvalue equation for the systems presented in the following sections, we have used standard MATLAB eigenvalue solvers, but standard LAPACK routines would have been an equally good choice (http://www.netlib.org/lapack/).

Results of numerical simulations
The structure in figure 2 is embedded within a rectangular computational area. A straightforward computational way is to discretize the entire rectangle as above and assign, for example, a large value to V R in the exterior shaded regions. The solutions ψ i, j therefore effectively vanish in those regions and the Dirichlet boundary condition ψ = 0 is satisfied automatically at the boundary of the cavity. This approach works for isolated cavities. It is less convenient, however, if one wants to implement the Neumann boundary conditions ∂ψ(x, y)/∂ x = 0, for example, at the two vertical ends of the stubs in an attempt to mimic infinite leads. In that case, we exclude the shaded areas from the computational grid. Additional steps are then taken to ensure that boundary conditions, either Dirichlet or Neumann, are properly included. Thus, the Dirichlet condition ψ(x B i , y B j ) = 0 applies for points at the boundary B, and the vector C is preconditioned accordingly before solving the matrix equation (13). The Neumann condition may be implemented in a similar way. To a good approximation, ∂ψ(x, y)/∂ x (ψ(x i+1 , y j ) − ψ(x i−1 , y j ))/2a, as above. We may now precondition C such that ψ(x i+1 , y j ) = ψ(x i−1 , y j ) at the two vertical ends of the stubs.

Closed quantum cavities
To test the general accuracy of the numerical procedure, we have solved equation (13) for some low-lying states for a closed cavity (V I = 0 in figure 2). The calculations were performed with the GaAs conduction band effective mass m * = 0.067m e , where m e is the bare electron mass [1]. The width and height of the inner cavity were chosen as 380 and 500 nm, respectively, and the radius of the upper rounded corners was set equal to 100 nm. The width of the two stubs was 70 nm and the radius 50 nm. These numbers are typical for the type of quantum dots we deal with here. Dirichlet boundary conditions were thus imposed throughout. To achieve an acceptable numerical accuracy, the number of grid points may be chosen as (75 × 75) or approximately so for the 300 lowest states. This is found by varying the number of grid points (N × M). The accuracy may also be estimated by comparing numerical and analytic solutions for a rectangular box of similar size. The relative error for the lowest state is then ∼0.01%; there is a linear increase with energy; for the 300th level the relative error is ∼3%. The numerical accuracy is thus quite sufficient for the purposes here.

Open quantum cavities and probability current density
To simulate the case of open quantum dots with interior currents, we introduce, as discussed above, imaginary potentials ±iV I in the two stubs, as indicated in figure 2. By this step our 'embryonic leads' thus become functional leads with a current flowing from one to the other, and we may think of the 'potential drop' 2V I as representing an applied voltage bias between the source and the drain. As above, we choose Dirichlet boundary conditions except for the two vertical ends of the leads, for which it appears more realistic to use the Neumann boundary conditions ∂ψ/∂ x = 0 as we wish to emulate extended leads. The parameters are the same as already listed above.   1.251 meV). Dark and light areas in the density plot correspond to low and high densities, respectively. The state is reminiscent of a 'particle-in-a-box' state with n x = 2 and n y = 7 that is distorted by the leads and the rounding of corners. The accumulation of density in the leads shows that the cavity is well open. mixing of adjacent states as V I is turned on, one should choose |V I | less than the spacing E between the energy levels of the closed unperturbed system, but besides this restriction the actual value of V I is not crucial in the present context. As a good estimate for E, one may use  , y). Therefore, V I in equation (8) in practice vanishes because (x, y) is in practice symmetric for our type of systems, as long as V I is chosen small. As a consequence, the eigenvalues become real to a very high degree of accuracy and the corresponding states may therefore be regarded as stationary. (For stronger sources, i.e. higher values of |V I |, the solutions become increasingly asymmetric. To satisfy equation (8), one would then resort to a self-consistent calculational procedure in which the left and right imaginary potentials V L and V R would be given different magnitudes.) Because we choose V I to be small, the densities displayed in figures 3-6 are very similar to the case of closed cavities. These figures also show that a current from one lead to the other is set. Because of the smallness of V I , the net current is also small, in particular for the case of the vertical bouncing mode in figure 6, as implied by the smallness of in the outermost regions of the leads. We will return to this point in the next section.
At the same time, as net currents may be small, local currents may be significant in any part of the cavities, as in the cases displayed. Typically, there are left-and right-handed vortices at the nodal points ψ(x, y) [17], [19]- [21] (cf figure 4). Because of the left-and righthandedness, the vortex currents tend to integrate to zero and, if any at all, normally give only small contributions to the net current through the system (for example, the integrated current passing a cross-sectional area in one of the leads). The complex current patterns, including the vortices, in figures 3-6 may be seen as the result of interference-the wave function is basically a scattering state that 'interferes with itself' because of the reflecting walls of the cavity.
As already indicated above, vortices of the kind discussed here are of a general nature, as they are related to the interference of a complex wave function. As will be seen in the following section, vortices may be observed in 2D microwave cavities in which the in-plane electromagnetic Poynting vector plays the role of the probability current density [5,17].  Figure 6. The same as in figure 3 but for the fifth vertical bouncing mode with E = 0.8117 meV, similar to 'particle-in-a-box' mode with n x = 1 and n y = 6. In this case, the cavity is only weakly open, as indicated by the low densities in the two leads.

Experimental studies with microwave billiards
The ultra-small structures we have discussed above are all in the nano-regime and therefore dominated by quantum-mechanical effects. Particle densities may be measured for special cases, like quantum corrals [3] on a Cu substrate, while at present it would be impractical to make similar observations for embedded ultra-small semiconductor structures, as in figure 1. It would be even more forbidding to measure the corresponding current distributions j without severely perturbing the quantum states. Measurements of the wave functions themselves are in principle impossible because of their quantum properties. As already indicated in the previous section, there are classical analogues that may come to help.
Fortunately, macroscopic planar microwave cavities may be used to emulate our 2D quantum dots provided that the geometries of the systems are chosen to match each other [5,17]. As explained in detail in [4], there is a one-to-one correspondence between the TM modes of the electromagnetic field and the wave functions of the corresponding quantum billiard. The z-component of the electromagnetic field E z , perpendicular to the planar billiard, then corresponds to the quantum-mechanical wave function ψ and obeys the Helmholtz equation, which in this case is identical to equation (9) for a particle in a hard-wall box. As mentioned, the Poynting vector is the analogue to the probability current density. The wave number for the microwave is k = ω/c, whee ω is the angular frequency of the TM mode and c the speed of light.
Microwave measurements are available for the rectangular cavity in figure 2 with rounded corners [13]- [15], [17]. The dimensions of the cavity itself are in this case (16 × 21) cm 2 , and the width of the two attached leads is 3 cm, i.e. there is a scaleup of our quantum billiard to the macroscopic world. Antennas placed in the leads, as in figure 2, act as the source and the drain for the microwaves. The field distribution inside the cavity has been obtained via a  probe antenna moved on the grid indicated in figures 7-9. The absorption through this third antenna may normally be made small, and for such cases it may be disregarded in the data analysis. In contrast to the quantum-mechanical case, one may thus measure the emulated 'wave function', its amplitude as well as its phase. Having this information one may evaluate the density distribution and the emulated probability current density. Concerning the densities, figures 7-9 show that there is an overall positive agreement between the microwave data and the simulations for the open quantum billiards in figures 3-6. The flow patterns show overall structural agreement but details differ, in part depending on the lower resolution in the measured data. There is, however, a principal difference between the two figures 7 and 8 and figure 9. In the first two cases, the wave functions are well connected to the leads. Consequently, there is a substantial current flow in and out of these regions. This is consistent with the theoretical simulations in figures 3 and 5. On the other hand, for the bouncing mode in figure 9, the coupling to both of the leads is much less, which is also evident in the simulated data in figure 6. As mentioned, the microwave technique relies on the fact that the probe antenna introduces a weak perturbation only. According to Stöckmann [22], 'this condition is violated if the transport through the billiard is blocked by anti-resonances which occur generically, or if there are patterns with strong standing wave contributions such as bouncing ball states. In these cases the probe antenna means a severe source of leakage spoiling the measurement of the current.' We recall that a probe antenna is not included in the theoretical model, which apparently gives rise to the difference in the current patterns in figures 6 and 9. In principle, it would, of course, be possible to extend the present model by introducing a third imaginary potential to represent absorption by the probe antenna.

Concluding remarks
We have introduced a basic heuristic model for numerical simulations of wave functions in open 2D quantum dots and cavities. Once wave functions are known, flow patterns, phases, nodal points and other topological features may be evaluated [17,21]. The model, which is an extension of the optical model for inelastic nuclear scattering [6,7], makes use of complex potentials, which act as input and output leads by feeding and absorbing particles. To avoid level mixing, the modulus of the imaginary potential should be chosen much less than the level spacing of the unperturbed system. To ensure that the eigenvalues are in practice real, the geometry of the quantum system should be symmetric; otherwise the imaginary potential V I must be designed such that V I = 0 in equation (8). The model discussed so far is a 2D model for non-interacting particles. However, it may be extended to 3D in a straightforward way. It should also be possible to extend the model to interacting particles and to include magnetic fields because the derivations in previous sections are quite general. For example, Varga and Pantelides [9] have simulated nonlinear transport in single molecules using the Kohn-Sham self-consistent equations in combination with imaginary potentials for the source and the drain, as here. In the case of nonlinear transport, which implies increased values for |V I |, we recall that the imaginary potential may have to be readjusted in each step of the self-consistent calculations to maintain real eigenvalues.
Finally, we have compared the results of the theoretical simulations for a weakly open quantum dot with measurements for an analogous microwave billiard. The overall agreement between theory and measurements is satisfactory, particularly for intensities, as long as the perturbations from the probe antenna may be neglected. For this reason it would be of interest to extend the present modeling to include absorption by the probe antenna. In the same way, leakage through, for example, the walls and other kinds of absorption may be incorporated into the model. Likewise, it would be of interest to see how different geometries of a billiard influence wave function forms and current flows, i.e. one may inspect the crossover from regular to chaotic behavior, the evolution of different statistical distributions, etc. It would also be possible to extend this type of modeling to other classical wave analogues, for example, classical acoustic reverberation rooms [23].
In this paper, we have used the finite-difference method to solve the Schrödinger/ Helmholtz differential equation. In a way, this is a brute force method but an effective one to get a quick computational start and early focus on the physics involved. However, it would be of obvious interest to improve on the computational aspects by implementing more efficient methods, like FEM, Green's function methods, surface integral methods, etc, which would be needed when dealing with more complex structures than here. So far, we have investigated weakly open billiards, a limit in which even second-order perturbation theory may be used.