PAMELA: An Open-Source Software Package for Calculating Nonlocal Exact Exchange Effects on Electron Gases in Core-Shell Nanowires

We present a new pseudospectral approach for incorporating many-body, nonlocal exact exchange interactions to understand the formation of electron gases in core-shell nanowires. Our approach is efficiently implemented in the open-source software package PAMELA (Pseudospectral Analysis Method with Exchange&Local Approximations) that can calculate electronic energies, densities, wavefunctions, and band-bending diagrams within a self-consistent Schrodinger-Poisson formalism. The implementation of both local and nonlocal electronic effects using pseudospectral methods is key to PAMELA's efficiency, resulting in significantly reduced computational effort compared to finite-element methods. In contrast to the new nonlocal exchange formalism implemented in this work, we find that the simple, conventional Schrodinger-Poisson approaches commonly used in the literature (1) considerably overestimate the number of occupied electron levels, (2) overdelocalize electrons in nanowires, and (3) significantly underestimate the relative energy separation between electronic subbands. In addition, we perform several calculations in the high-doping regime that show a critical tunneling depth exists in these nanosystems where tunneling from the core-shell interface to the nanowire edge becomes the dominant mechanism of electron gas formation. Finally, in order to present a general-purpose set of tools that both experimentalists and theorists can easily use to predict electron gas formation in core-shell nanowires, we document and provide our efficient and user-friendly PAMELA source code that is freely available at http://alum.mit.edu/www/usagi


I. Introduction
The unique electronic properties of semiconductor nanowires continue to be an area of intense interest for diverse applications in optoelectronic, electromechanical, and actuation nanosystems. 1 Because of their small cross sections, quantum confinement effects in nanowires are dominant, and electrons occupy discrete energy levels that are considerably different than the continuum of energy levels found in bulk materials. An area of immense technological interest is the creation of core-shell nanowires (Fig. 1a) where the material composition and/or doping concentration are modulated in the radial direction of the heterostructure nanowire. 2 The core-shell nanowire structure provides a unique advantage compared to homogenous nanowires since it allows confinement of two-dimensional electron gases (2DEGs) at the semiconductor-semiconductor heterojunction interface. [3][4][5][6][7][8][9] For the specific case of a GaN/AlGaN core/shell nanowire as shown in Figs. 1a and 1b, carriers are confined within the GaN core by the larger bandgap of the AlGaN shell. In order to optimize performance of these nanowires, it is essential to know which set of material properties (i.e., bandgap alignment, material composition, cross-sectional size, doping density, or many-body electronic effects) favor the spontaneous formation of a mobile electron gas in these nanosystems. Since the experimental parameter space for controlling The most commonly used approach for predicting electronic properties in these nanowires is a simplistic Schrödinger-Poisson formalism 10 4 is solved, where m * is the effective mass, i  is the electron wavefunction for state i, E i is its energy, V n is the spatially dependent free carrier concentration due to the n-type doping density, V CB is the conduction band edge, and V e-e is the (local) Poisson potential due to electron-electron interactions.
Within the conventional Schrödinger-Poisson formalism, the Poisson potential, V e-e , is given by  ( 2) where N(E F ,E j ,T) is the integrated one-dimensional density of states along the nanowire, and the summation in Eq. ( 2) is over all occupied orbitals, j  . However, since the summation is over all j, an electron feels the electrostatic potential from the other electrons as well as from itself. This unphysical self-interaction error is well-known from electronic structure theory and catastrophically leads to overdelocalized electrons, underestimated band gaps, and unstable charge states. 13 where V e-e is now an orbital-dependent potential. In other words, the difference between Eq. ( 2) and ( 3) is the exclusion of the "diagonal term," j = i, which can be interpreted as the interaction of the electron with itself. However, the use of Eq. ( 3) (which, again, is still not quite correct as it ignores antisymmetry of the electronic wavefunction) poses severe computational difficulties in the conventional Schrödinger-Poisson approach since this orbital-dependent operator is different for every orbital (because of the restricted summation over j ≠ i), and one would have to solve a different effective Schrödinger equation for each individual orbital in the self-consistent approach.
Herein, we provide a simple and efficient approach for rigorously eliminating self-interaction errors and assessing the effects of many-body, electronic interactions on the formation of electron gases in core-shell nanowires. Our approach utilizes a pseudospectral numerical method 18 Approximations) available online. 20 We specifically focus on the effects of nonlocal exact exchange in nanowire systems for three important reasons: (1) Very recent work by us and others have shown that electron gas formation is insensitive to local correlation effects. 11,12 While results obtained with a simplistic local density approximation (LDA) were nearly identical to uncorrelated calculations, the effects of nonlocal exact exchange on large core-shell nanowires are not known; (2) We and other researchers have also previously shown that exchange-correlation effects are inherently nonlocal, and electronic properties in molecules and bulk systems can be significantly improved with nonlocal exchange, [21][22][23][24][25][26][27][28] and (3) Electronic wavefunctions with exact exchange are usually the initial reference state for including higher-level treatments of electron correlation, 29-32 so the formalism used in this work serves as a starting point for even more extensive treatments of strong correlation effects in nanowires.
While the incorporation of nonlocal exchange in atomistic calculations of molecules and materials is well-known, we are unaware of any previous work assessing the effect of exact exchange on large coreshell nanowires.
The work presented in this paper is divided into two separate parts: In the first part (Sections I-IV), we give a brief derivation of the exact-exchange equations for electrons in a core-shell nanowire, and we demonstrate how to reduce the integro-differential exact-exchange equations to simultaneous (coupled) partial differential equations that can be more easily solved numerically. Using this approach, we present specific situations where exact exchange effects become dominant and give several examples showing electron-occupancy levels, band-bending diagrams, electron density delocalization patterns, and electron subband energies. In the second part of the paper (Appendix Sections 1-9), we give a more detailed description of the various approaches and numerical routines which comprise the PAMELA software package that is available online. 20 The availability of the efficient and user-friendly PAMELA source code (< 400 lines, including several lines of comments) is intended to reduce the technical burden of writing computer code from scratch and to present a general-purpose set of tools 6 that both experimentalists and theorists can easily use to predict electron gas formation in core-shell nanowires. With these new approaches and predictive methods, we highlight several areas where manybody quantum effects play a significant role in both the formation and distribution of electrons in these low-dimensional nanostructures.

II. Theory and Methodology
The Hamiltonian for N electrons in a nanowire is given by ˆn V is the spatially dependent free carrier concentration due to the n-type doping density, and where an overbar denotes that a spatial orbital, In Eqs. (8) and (9), N(E F ,E j ,T) is the integrated electron concentration given by where g(E,E j ) is the one-dimensional density of states, E j is the energy of the jth wavefunction, E F is the Fermi energy, T is the temperature, and f(E,E F ,T) is the Fermi-Dirac distribution (all described further in Section 4 of the Appendix). Note that the exchange interaction is a nonlocal operator (where i  is inside the integral) and involves an "exchange" of electron i and j to the right of 1 ij rr  in Eq. (9) relative to Eq. (8). Also note that the summation in Eq. (7) is unrestricted and includes terms where i = j.
The unrestricted summation is a unique property of using an antisymmetric wavefunction since the direct and exchange interactions for a single electron are equal but opposite in sign, and the selfinteraction cancels out exactly 33 (this is easily seen by inspecting Eqs. (8) and (9), and showing that Schrödinger-Poisson approaches. [10][11][12] Using the same techniques as in the treatment of the direct interaction term, 19,34,35 the exchange interaction between electrons is and acting on Finally, the potential energy due to the n-type doping density is given by solving where ρ D is the charge density arising from the ionized dopants. It is important to note that we have derived an approach where the integro-differential exact-exchange equations have been reduced to simultaneous (coupled) partial differential equations without any nonlocal integral terms.
Eqs. (11), (13), (15), and (16) define our zero reference of energy to be at the shell rr  conduction band edge). In addition, we determine the position of the Fermi level via a charge-neutrality condition (described in Section 4 of the Appendix) that adjusts the Fermi level at each iteration step.
The specific material system we consider in this study is a cylindrical nanowire with a GaN core and an Al 0.3 Ga 0.7 N shell. We assume the GaN/Al 0.3 Ga 0.7 N layers are in the Wurtzite structure oriented with the [0001] direction along the axis of the nanowire (see Fig. 1a), which we also assume to be defect-free. 5 While we focus on this particular material system, we expect our findings to hold for other similar core-shell nanowires with Type I (straddling) heterojunctions. The bulk electronic properties and material parameters used in our self-consistent calculations are given in Table 1.

III. Results and Discussion
Numerous simulations with the PAMELA program 20 were carried out on cylindrical core-shell nanowires with GaN core radii ranging from 10 nm to 50 nm with a fixed 20-nm-thick shell layer of Al 0.3 Ga 0.7 N. We investigated N-type doping densities ranging from 10 16 to 1.2×10 18 cm -3 , assuming a uniform doping density throughout the nanowire. In order to numerically solve the coupled pseudospectral differential equations, we chose 79 Chebyshev points to discretize the radial direction, and 24 equispaced points were chosen for the θ discretization (see Sections 6-9 of the Appendix), resulting in matrices of rank 936 used throughout our calculations.
We first discuss the effects of nonlocal exact exchange on the number of occupied wavefunctions in our nanowire, as demonstrated in the contour plots in Fig. 2. Calculations without exact exchange consistently overestimate the number of occupied levels, which become even more noticeable at higher doping densities and at higher core radii. At smaller doping densities and small core-radii, exact exchange has a noticeably larger effect on the electron density profile within the cross section of the nanowire, as shown in Figs. 3a-d. For a small nanowire with a 20-nm-wide core, a 20-nm-11 thick shell, and a doping density of n D = 0.5×10 17 cm -3 , we find that the electrons stay localized at the center of the core when nonlocal exact exchange effects are incorporated. In contrast, we find that calculations without exchange (using the same doping density as before) predict a charge density that is more delocalized, showing a pronounced annular distribution. In particular, we obtain an electron density with a peak height that is reduced by almost half when exchange interactions are not incorporated.  We see similar localization effects in nanowires with larger core sizes and higher doping densities. In Figs. 4a and b, we use the PAMELA program to plot the free electron gas distribution for nanowires with an 80-nm-wide core, a 20-nm-thick shell, and a doping density of n D = 5.5×10 17 cm -3 . In both cases we see small variations in the electron density in the middle of the core region; however, the localization of the electron density at the core-shell interface is sharper with exact exchange, whereas the electron gas without exact exchange decays slowly towards the nanowire center. The doping density used for both cases is n D = 5.5×10 17 cm -3 . With exact exchange, the free electron gas is more localized at the heterojunction interface, whereas the electron gas without exact exchange decays slowly towards the nanowire center.
In order to explain the electron delocalization patterns in Figs. 3-4, we must take a closer look at the effective electron-electron potential energy in a nanowire. As derived in Eq. (11), the effective potential energy due to electron-electron repulsion in a nanowire is given by a difference of two terms:  (11)) and compensates for the exaggeration due to the D, j V term, leading to a more localized electronic density. As a result, the inclusion of quantummechanical exchange plays a crucial role in correcting the over-delocalization of electrons in nanowire systems that arises from including only the bare Coulombic interaction.
We can further analyze the effects of exact exchange on the energy eigenvalues (or quasiparticle energies) of the electronic subbands outputted by the PAMELA program. As shown in Fig. 5, we plot the lowest 20 electron subband eigenvalues as a function of doping density for calculations with and without nonlocal exact exchange. Due to the cylindrical symmetry of the nanowire, some eigenvalues are either singly or doubly degenerate and are shown as blue or red lines, respectively. In general we see that the energy spacings between subbands are smaller in the absence of exchange, while exchange effects dramatically increase these energy separations (i.e., the subbands are more spread out in Fig. 5a).
This widening of energy spacings is consistent with many-body electronic effects where it is wellknown that nonlocal exact exchange widens quasiparticle gaps in both molecules and bulk materials. 23,24,26 In contrast, this nonlocal many-body electronic effect is completely absent when solving the conventional Schrödinger-Poisson equations without exchange, leading to energy differences between subbands that are considerably underestimated. It is also important to note that the relative energy spacings between subbands further explain the overestimation in the number of occupied levels shown in Fig. 2b. Specifically, the total electron density, n e , is proportional to the difference between the Fermi energy and the occupied electronic subbands: However, in order to satisfy the charge neutrality constraint, we require the total number of positive and negative charges over the entire nanowire to be equal: De dA n dA n   (where A is the cross-sectional area; see Section 4 of the Appendix). If the energy separation between subbands is small (as it is for the case without exchange), then the various Fj EE  terms will also be small, and we will require considerably more occupied levels to satisfy the charge neutrality constraint. For the same reason, since 15 the energy spacings between subbands are considerably larger when exact exchange is incorporated, the Fj EE  terms are also larger, and we do not need as many occupied levels to satisfy charge neutrality (cf. Fig. 2a).
where A is the cross-sectional area of the nanowire. By construction, 2 r measures the relative spread of the electron density and is bounded between 0 and 22 shell core rr . From its definition, values of 2 r less than 1 indicate localization within the core, and values greater than 1 represent confinement within the shell layer. In Fig. 6a, we plot the core-normalized variance as a function of doping density for a set of 16 representative core radii. At low doping density (< 2×10 17 cm -3 ), we find that 2 r rises rapidly and eventually saturates to ~0.7 at moderate doping. This variation in 2 r corresponds to a shift in electron localization from the center of the nanowire to a spreading of electron density within the core as the doping density is increased. At a critical doping density (~ 7×10 17 cm -3 for these particular core radii), Fig. 6a shows a dramatic, abrupt change as 2 r rapidly increases towards 1. This sudden change in behavior corresponds to a transition where the electron gas in the middle of the core becomes strongly localized at (or even beyond) the heterojunction core-shell interface. It is important to note that these trends for 2 r still persist whether exact exchange is incorporated or not; however, the transition to an interface-localized state occurs at higher doping densities when exact exchange is included. This interesting behavior in 2 r has a direct correspondence with the position of the Fermi level in the coreshell nanowire. Fig. 6b shows the Fermi energy as a function of doping density for the same set of core radii considered previously. As the doping density is increased, the Fermi level rises linearly as we require a larger number of occupied electrons to satisfy charge neutrality. However, at exactly the same critical doping densities shown in Fig. 6a, we find that the Fermi energy in Fig. 6b has reached a saturation point such that further increases in doping density do not dramatically affect the position of the Fermi level. We can understand these surprising trends in both 2 r and the Fermi Energy by using the PAMELA program to plot the conduction band and electron density for the nanowire near this critical density. Fig. 7a depicts the band-bending, Fermi energy, and electron density for a nanowire with a 30nm-radius core and a 20-nm-thick shell with a doping density of 7.5×10 17 cm -3 . This specific doping density corresponds to the critical density for the 30 nm-radius core nanowire when exact exchange is included (blue crosses in Fig. 6a). In Fig. 7a, the band diagram shows that the Fermi energy is positioned well above the conduction band at the outer edge of the nanowire. For doping densities beyond this critical value as shown in Fig. 7b, tunneling into the shell layer occurs, and we observe some localization of the electron gas at the outer edge of the AlGaN shell. It is at this doping density that a sharp transition in 2 r begins to develop in Fig. 6a, and as the doping density is increased further, 2 r even exceeds 1 as the electron gas moves to the outer shell edge of the nanowire. It is interesting to note that we only observe these effects when the Fermi energy has reached a level where the shell barrier becomes thin enough for tunneling to occur. This further implies that there is a critical tunneling width (of the shell), that once reached, allows the electron gas to tunnel from the core into the shell layer. In all of our self-consistent calculations, we observe the onset of tunneling to occur when the width of the barrier is approximately two-thirds of the shell thickness, regardless of core radius.

IV. Conclusions
We have both demonstrated and provided an efficient approach for calculating the effects of nonlocal exact exchange on heterojunction electron gases in cylindrical core-shell nanowires. Our approach utilizes a pseudospectral numerical method that is efficiently implemented in the user-friendly PAMELA software package available online. 20 The extensive documentation and availability of our software package is a unique feature of this work since it provides a self-contained set of tools that both experimentalists and theorists can use to explore the effects of bandgap alignment, material composition, cross-sectional size, doping density, and many-body nonlocal exchange on electron gas formation in these core-shell nanowires. Furthermore, the availability of the fully open-source PAMELA computer code encourages researchers to get a detailed "look under the hood" to obtain a deeper understanding of how these effects are numerically incorporated in practice. Using the PAMELA software package, we find that the inclusion of exact exchange has a crucial effect on the localization of electrons in a GaN/AlGaN core/shell nanowire and leads to a much larger energy separation between individual electronic subbands compared to simplistic treatments without exact exchange. In addition, we also investigated nanowire systems in the high doping density regime and found that tunneling from the core-shell interface to the nanowire edge becomes the dominant mechanism of electron gas formation. With the ability to efficiently carry out calculations with the PAMELA program, this work plays an important role in understanding which geometries and material parameters favor the spontaneous formation of electron gases in these nanoscale systems. Looking forward, it would be extremely interesting to understand the effects of exact exchange on both different cross-sectional geometries (i.e., hexagonal and triangular cross-sections 11 ) and excited-state excitonic properties. Since electron/hole energies are extremely sensitive to the energy separation between electronic subbands, nonlocal exchange effects are expected to be essential for accurately predicting these excitations. These studies are currently underway in our group and will be presented in due course.

Schrödinger Equation
Assuming translational invariance along the nanowire z-axis, the two-dimensional Schrödinger equation in the effective-mass approximation (c.f. Eq. (11)) is given by where m * is the effective mass, i  is the electron wavefunction for state i, E i is its energy, V n is the spatially dependent free carrier concentration due to the n-type doping density, V CB is the conduction band edge, V D is the (local) potential of the direct interaction, and V EXX is the (nonlocal) potential due to exact exchange. In the PAMELA source code for the Schrödinger equation, we make use of the dimensionless variables 2 Cm  , m 0 is the mass of an electron, and 0 is a natural length scale that we define to be 1 nm. In terms of these dimensionless variables, Eq. (A1) reduces to: where the kinetic energy term in (dimensionless) cylindrical polar coordinates is As discussed briefly in Chapter 9 of Ref. 18, the pseudospectral method computes eigenfunctions that are linearly independent but not necessarily orthogonal to each other. Since we require the wavefunctions to be mutually orthogonal, we use Löwdin's procedure (also known as symmetric orthogonalization 33 ) in the PAMELA program to orthogonalize our wavefunctions. This procedure is carried out by constructing the following N/2 × N/2 overlap matrix: where  ψ and ψ are column-matrices of orthogonal and non-orthogonal wavefunctions, respectively.

V n and the Direct Interaction
Within the PAMELA source code, both V n and V D are solved simultaneously in a single equation (both V n and V D are local potentials): where j   are orthogonal wavefunctions (see end of Section 1 in the Appendix), and we have made use of Eqs. (16) and (13). The integrated electron concentration is given by

The Fermi Energy and Charge Neutrality
We determine the Fermi level, E F , by requiring the total number of positive and negative charges over the entire nanowire to be equal: where shell r is the radius of the outer shell region. Since n e is a function of E F (see Eq. (A24)), the Fermi energy is determined by solving Eq. (A25) with a standard root finding algorithm within PAMELA.

Boundary Conditions and Iteration Procedure
As discussed in the Theory and Methodology Section, we assume Dirichlet boundary conditions  in the Schrödinger equation, and the exchange interaction is a nonlocal, integral operator (the inverse matrix operator is mathematically equivalent to an integral operator).
The coupled set of Schrödinger, direct, and exchange interaction equations in Eqs. (A2), (A27), and (A29) are solved iteratively until self-consistency is achieved. For our starting guess to initialize the self-consistent procedure in the PAMELA source code, the bare conduction band offset (Eq. (A5)) for the core-shell system is used as a seed potential for the Schrödinger equation to generate a set of normalized wavefunctions and eigenenergies. The wavefunctions are orthogonalized using the Löwdin orthogonalization procedure, and these solutions are then subsequently used in the calculation of the Poisson and Exchange interaction potentials. The resulting potentials are then incorporated into a new Schrödinger equation for the next iteration, and the entire process iterates until convergence. To maintain a stable cycle for iteration to self-consistency, we use an under-relaxation technique on the matrix representation of the total interaction potential with an under-relaxation parameter ω (set to 0.05 in the PAMELA source code) such that after iteration i, the matrix representation of the total input potential (V input,i ) is calculated as   input, input, 1 CB D, EXX, input, 1 2.
These iterations are continued to self-consistency, which we have defined as the maximum change among all the matrix elements of the full interaction potential between successive iterations to fall below 0.001 eV.

Pseudospectral Methods
As discussed extensively in Ref. 18, pseudospectral methods are an efficient and accurate approach for numerically solving partial differential equations and computing integrals. In contrast to finite difference or finite element methods that approximate a differential operator by a large and sparse