A Numeric Approach for Investigating Electron Dynamics in Zinc‐Blende Semiconductor Heterostructures

Charge transport via interfaces in materials is crucial in comprehending the constraints of many electronic devices. Gaining the knowledge and understanding of the processes entail when two materials form contact between them may allow designing better performance devices. Nonetheless, simulating experimental length scale structures puts a computational challenge. This work investigates several heterostructure systems of typical semiconductors with zinc‐blende unit cells forming an interface in the contact plane. The purpose is to understand the trends in the electronic properties of such heterostructures regarding the strain profile, as it occurs near the interface and far away in deeper atomic planes. The local band‐gap values of each layer of atoms are calculated to demonstrate how the gaps have changed with respect to the distance from the interface. Finally, the potential energy at the interface region is parameterized. The correlation between the transmission coefficient of the one‐electron wave function is analyzed to pass through the interface and interface characteristic parameters. The primary conclusion is that the system's electrostatic potential embeds the heterostructure's unique properties. It is demonstrated how one can predict the transmission coefficient of the electronic wave function at the interface by examining the material's potential parameters at the interface exclusively.


Introduction
A fundamental aspect of a material's microstructure evolution occurs as two materials featuring different crystallographic structures are combined.A better understanding of the interface nature will facilitate controlling material properties from the atomic scale.When it comes to investigating systems for brand-new DOI: 10.1002/adts.202300155applications, there is insufficient information about the preferred materials for the applied system, making this task even more complex.Modern technology has evolved dramatically, developing modern electronic and optoelectronic devices utilizing quantum phenomena that occur when an interface between two semiconductors is formed.In a famous paper published in 2001, the catchphrase "the interface is the device" was coined, emphasizing the central role of interfaces in device functionality. [1]herefore, enhancing the performance of the heterostructure-based device can be achieved by better understanding their interfaces' properties.
Interface properties include the alignment of the electronic bands affecting the charge carriers' transport as well as the lattice defects and impurities at the materials' crystal lattice.These factors affect the mobility of both electrons and holes.Those also influence the recombination processes inhibiting charge transport along the materials.The formation of an interface at the contact between two phases might also result in a strain that alters the electronic structure of both materials. [2]escribing interfacial processes in heterostructures more accurately is becoming possible with advanced algorithms based on the quantum mechanics interactions between the nuclei and the electrons at the nanoscale.
First-principle modeling methods for atomistic level descriptions have gained a lot of interest and have become frequently used for investigating interfacial phenomena.Such techniques require accurately predicting the actual modeled interface's atomic geometry.However, the prescription for obtaining the atomic geometry for an interface is not trivial.Since these calculations deal with the crystalline system, the structure configuration should satisfy periodicity.Therefore, searching for a common supercell that contains the constructed interface between the materials is needed.This may turn out to be a challenge since it is almost impossible to fit two semiconductors having two distinctive crystal structures perfectly.When combining two semiconductors to form a unified supercell, straining one of the materials with respect to the elastic constraint the other component imposes is inevitable.Some approximations have been suggested to take the average between the two semiconductors' lattice constants. [3]Still, the mismatch in the lattice constants results in strain profiles in both materials.The one with the longer lattice parameter experiences biaxial compression strain; alternatively, the other is biaxially tensed.The strain profiles occur in the directions parallel to the interface's plane, where strain conditions are imposed on the system.Many considerations may be accounted for when modeling interfaces.Extensive atomistic models are required to address a more realistic description for investigating materials' interfaces.An automated method for determining the optimal supercell representation of interfaces was proposed for searching for the optimal orientation, relative angle, lattice matching, and elastic energy.The introduced algorithm provides a robust scanning tool for finding common supercells between materials' surfaces. [4]nother atomistic modeling approach has been suggested to minimize a heterostructure's total energy, thereby optimizing the interface configuration.Then the generated relaxed heterostructure interfacial strain can be evaluated while relaxation is performed with the Metropolis Monte Carlo scheme. [5]ur study involves computational and numerical analysis of conventional zinc-blende semiconductors.To address the impact of interface strain on charge transport, we created a supercell where two semiconductors form an interface.One of the cubic zinc-blende type semiconductors advantages is the broad range of band-gap values that characterizes them.The majority of this group are direct band-gap materials which is a desired feature when referring to energy conversion in terms of efficiency.Besides the variety of interests this materials family attracts, they also feature a simple crystal structure, making them the perfect candidates for examining the model since using them is convenient for constructing heterostructures.In this work, we investigated the effect of strain on the electronic structure and charge dynamics properties of zinc-blende heterostructures.To assess charge transport through the interfaces, we employed a technique that involves performing density functional theory (DFT) calculations on small systems and subsequently using the results for larger-scale supercells as we used a wave-packet propagation.

Computational Methodology
DFT calculations were conducted using VASP, [6,7] version 5.4.4.The generalized gradient approximation (GGA) functional of Perdew-Burke-Ernzerhof (PBE) [8] was used as it provides accurate lattice constants, which agree with experimental results at a relatively low computational cost.Nonetheless, it is known that GGAs tend to underestimate materials band gap values.Since the interest was in studying only trends in the band gap values concerning the interface implications, while insisting on precise band gap values might be superfluous, the choice of GGA functional of PBE served well.Justification for this approach can be found in reports indicating that the band alignment obtained with GGA demonstrated only a slight difference compared to the results obtained using hybrid functionals, which typically provide more precise band gap values, albeit at a higher computational cost.
The good agreements between band alignment results from GGA-based calculations and hybrid functional-based calculations supported this practice. [9]he bulk materials' core electrons were described by projector augmented-wave (PAW) pseudopotentials provided by VASP.The valence electron configurations for each atom can be found in Table 1.For each bulk structure, the relevant computational information can be found in Table 2 regarding the implemented lattice parameter, k-grid density, energy cut-off values, and the lattice unit cell symmetry.The tolerance value for both plane wave cut-off energy and k-point grid size convergence test was set to 1 meV atom −1 .The heterostructures were built from two different semiconductors posing two characteristic computational parameters, as shown in Table 2: the k-points mesh density and the energy cut-off value.Whenever a conflict between the input parameters arose, the stricter one was chosen to produce the calculations.A denser k-point mesh was sampled for the density of states (DOS) calculations to produce more accurate and detailed plots, usually set to be 11 × 11 × 1 mesh size.For bulk band structure calculations, the KPOINTS file was configured to describe the highsymmetry path along the first Brillouin zone (BZ).It was autogenerated for the zinc-blende unit cell in real space, using the Python package PYMATGEN, [10] applying 40 k-points between every two adjacent high-symmetry points.

Constructing Semiconductors Heterostructures
After sampling some well-known zinc-blende semiconductor candidates for the case studies, a superlattice structure was constructed for each combination.An example of a zinc-blende cubic unit cell of AlAs is shown in Figure 1.Several considerations were taken while building the heterostructures.First, structures that had only a minor lattice mismatch were pursued.For this reason, it was essential to take care of finding the optimal common lattice constant to minimize the lattice mismatch related to the multiplications needed to be applied for each unit cell.Second, the final form of the heterostructure should have had an equal number of atomic layers for each one of the materials forming the interface.In this way, the correct stoichiometry could be preserved and a symmetric interface could be ensured, thereby avoiding polarizations at the interface.In Figure 2, a 2D representation of such heterostructure as presented in AlAs-GaAs heterostructure is demonstrated.
Accordingly, regions within the superlattice bounded between two interfaces should be sufficiently broad to investigate the influence of interface-induced strains at different distances from the interface locations.A supercell must be constructed by a large number of atoms in order to validate such a situation.Nonetheless, it was crucial to maintain a proper balance between the computational cost, which dramatically rises as the number of atoms increased, and an accurate description of the investigated model.Therefore, it was decided to build heterostructures consisting of ten unit cells from each semiconductor, making it 20 unit cells overall (each supercell contains 160 atoms).In the case study of GaAs-AlAs heterostructure, it was also investigated how the number of atomic layers affects the local band gap.This analysis was conducted considering that each group of atoms had different distances from the interface.These findings revealed that by reducing the length of the system, the local band gaps approached the average values observed in bulk structures.The calculation for bulk GaAs resulted in 0.153 eV band gap and for bulk AlAs 1.499 eV band gap.When only a single unit cell of each semiconductor from both sides of the interface was constructed, the resulting supercell exhibited a consistent band gap of approximately 0.9 eV across all atomic layers.As the number of the unit cells of each semiconductor composing the heterostructure was extended, an increasing variation in the band gaps for each atomic layer was observed.The results regarding this part are presented in the Results section, demonstrating the largest constructed supercells.
Last, the orientation relationships between the two crystal surface planes were selected.A popular epitaxy growth direction documented in the literature for such zinc-blende semiconductors is [001], perpendicular to the crystallographic c direction.This crystal growth direction yields high-quality structures, avoiding piezoelectric fields and providing enhanced control of the growth procedure. [3]The primary motivation for modeling these heterostructures having interface contact on the (001) planes came from that.The mentioned orientation relation with no relative angle between the contact planes of the two phases was kept, simply forming an abrupt interface.

Wave Function Time Evolution
A primary objective of this work was to numerically calculate the transmission coefficient of an electron passing through a semiconductor heterostructures interface.The strategy was to propagate the one-electron Schrodinger equation solution in time and space, applying the Kohn-Sham (KS) potential outputs from the DFT calculations.
The initial one-electron wave function was in the form of a Gaussian wave packet and was first centered at an atomic site, relatively close to the interface.The initial electron position was about one and a half atomic layer distance from the interface at a potential curve minimum point.Locating the initial one-electron wave function at the specified place not too far from the interface was done to save simulation time of propagating toward the interface.As the initial wave function evolved in time and propagated throughout space, the accumulative probability for the electron to be located after the interface position could be calculated.
The time evolution of the wave function was simulated by applying a time-propagating operator.The simulation was then initiated by advancing the wave function by a sequence of single time-steps.The total simulation time would be summed up through the complete series of the single time steps, equals to the time step size multiplied by the number of the simulated time steps.
To achieve the transmission coefficients of the one-electron wave function to get transferred through the interface, one should calculate the probability of current flux at the interface position, integrated over the entire simulation time [11,12] ⃗ where ⃗ J(z, t) is the probability current flux, which exhibits dependencies on time, denoted with t, and space, denoted with z. m e is the electron mass, Ψ is the one-electron wave function and Ψ * is the complex conjugate of the wave function.For calculating the electron transmission coefficient, the probability current flux, as given in Equation ( 1), should be evaluated at the interface position, ⃗ J(z = z interface , t), and then be integrated with respect to simulation time where Φ is the transmission coefficient, t final is the total simulation time, and t 0 is the initial time, regularly set to be t 0 = 0.Alternatively, the accumulative probability of the electron to be found after the interface might be obtained by integrating the squared absolute of the wave function over the space ranging from the interface position to the far system edge.Then, for a long enough time, the calculated accumulative probability should stabilize and reach a plateau whose value represented the transmission coefficient where Φ is the electron transmission coefficient, t → t final marks the limit when the simulation time reaches a plateau, t 0 is the initial time, z_interface is the interface position, and z = system_length is the coordinate that describes the end of the system.
Nevertheless, it was crucial to consider the impact of edge effects on the calculation results as the wave function was transferred or reflected to the other side of the system.To mitigate this issue, avoiding excessively long simulation times was necessary.Two convergence tests were conducted to prevent the wave function from reaching the edges and to obtain reliable results.The first test ensured that the overall simulation time and the supercell's length were converged, which was achieved by observing a plateau in the accumulative probability value.This approach ensured that the simulation was not too long but sufficiently long to obtain accurate estimates of the charge transmission extents.
The difference in the total energy between consecutive iterations was enforced to be less than 0.005 eV.A tolerance value of 0.005 as the difference in the transmission coefficients of two successive iterations was also determined.The converged time steps and the system lengths achieved under those constraints varied for each system since the convergence tests were carried out repeatedly in each case (reported in Table S3.0.0,Supporting Information).
It was chosen to utilize the split-step operator approximation method for implementing the time evolution simulations of the one-electron wave function.The solution to the Schrodinger solution was propagated by applying the propagator to the initial wave function The time variable is t, and the initial time is denoted as t 0 .The time propagator is a unitary operator assuming the Hamiltonian Ĥ is hermitian, and it acts on the wave function in the following way The time-propagator has in its exponent term the Hamiltonian, which in turn can be expressed as the sum of the kinetic and potential energy operators.An exponent of this form, having a sum of operators, can be expressed as the multiplication of the operators in an individual exponent only if the operators commute with respect to each other Consequently, the mathematical procedure becomes complicated as the two energetic operators in the exponent do not commute and can not be separated.Instead, by expanding the expression for the time-propagator using a Taylor series, neglecting the Δt 3 related terms, one can get an approximated expression that held only for small enough time steps Δt, and has a general error in the order of O(Δt 3 ) [13][14][15] e Then, the Hamiltonian operator in the exponent term is decomposed into its ingredients, the sum of the kinetic and potential energy operators.
A practical and efficient operational scheme can be established by applying the fast-Fourier transform (FFT) algorithm, thus converting the function into the frequency domain where the derivatives can be taken more easily.The mathematical formalism for this action is shown in the following equation In order to clarify the algorithm operations, grouping of some terms is marked in Equation ( 8) The recipe for calculating |Ψ(t + dt)⟩, which stands for each time step, goes in the order presented by groups in the last formula.First, one should calculate the multiplication of the wave function and e − i 2ℏ Vdt as shown in the term I. Applying the FFT algorithm getting term II, which result was then multiplied by e Tdt founding III term.Afterward, the result gained up till this stage should be plugged into the inverse Fourier transform denoted as FFT −1 , bringing the calculation back to the time domain; this would yield the IV term.The entire process is completed by computing the V term, which is determined by multiplying the previous result by e − i 2ℏ Vdt .The last point to be elucidated is the initial wave function definition.The initial wave function is a superposition of the oneelectron Schrodinger equation eigenfunctions represented by a wave packet of a Gaussian form function, centered around z 0 , and characterized by a width of , featuring the standard deviation of the function as well.The Gaussian wave packet is a well-known initial wave-function for the time evolution of a quantum system One can notice that the Gaussian wave packet is characterized by a kinetic energy magnitude manifested in k 0 , through the subsequent relations: where m e is the electron mass, p is the momentum, and k is the wave vector.The initial k 0 is determined via the analytical solution for the following equations, derived from the expectation value of the kinetic energy operator [12] T 0 where erf is the error function, L is the system length, z 0 is the Gaussian wave-packet center,  is the Gaussian width, and m e is the electron mass.T 0 (k 0 ) is the initial kinetic energy whose momentum is extracted by equating Equations ( 14) and ( 15) where there is a subtraction of the potential energy expectation value V(t 0 ) given in Equation ( 12) from the initial total energy of the electron E tot (t 0 ), delivered from the DFT calculation as the energy of an electron located at the edge of the conduction band.

1D Electrostatic Potential
A typical potential energy output file called the LOCPOT output file from VASP is a 3D quantity describing a 3D system.The simulation was designed to represent transport along the interface axis, and as such, a 1D vector was solely concerned.To achieve this, it was necessary to integrate over the perpendicular direction to the heterostructure length, which was typically the z-direction Accordingly, the 1D special grid was divided into N grid samples.Thus, the electrostatic potential vector depended upon only one spatial coordinate and had the exact discretization as the sampled spatial grid.
It is well known that DFT has a system size limitation due to computational cost.To simulate carrier dynamics at a larger scale, such as several atomic layers, it was necessary to extend the length of the electrostatic potential vector.Areas with bulklike properties within the electrostatic potential range were first located.These areas were then multiplied to create an extended electrostatic potential vector.This procedure allowed to stitch together several electrostatic potential regimes for larger-scale simulations.

Bulk Material Calculations
We chose to construct well-studied zinc-blende semiconductors heterostructures as our case studies.The bulk system calculations included results for two imposed strain profiles, contraction and expansion.In addition, we obtained the corresponding DOS and band structure plots.Finally, the local potential for the strained systems, together with the unstrained one, was drawn to understand the impact of the strain on the materials' potential.
All the structures were imposed to biaxial tensile and compression in the direction of [010] and [100] perpendicularly to the interface plane (001).Figure 3 does not reveal a clear and consistent trend.As previously mentioned, the influence of strain on a material is contingent on its electronic structure and chemical properties.Our show different groups of materials that exhibit distinct behavior.The first group comprises materials whose band gap widens when subjected to compression.These were GaAs, Ge, GaSb, and GaN.Except for Ge, all these materials are a compound of Ga atom, which may hint at its chemical bonding characteristics.We can also see that AlAs and AlN increase their band gap values when relating to the tensed structure.
The other group of materials was of Si and AlP that, under biaxial compression, showed narrowed band gap values concerning the unstrained form.Finally, the tensed structure of bulk InAs exhibited higher band-gap values when compared to both the unstrained and compressed structures.
There are several approaches to quantitatively analyze the influence of different types of strain on the material's electronic structure.It can be investigated through the perturbation theory under the KP framework (the k•p perturbation theory, or 'K dot P'), the tight binding method, and potential deformation. [16]ur analysis of the strain-induced effects was based on chemical intuition and some chemical bonding considerations.We have suggested some explanations that were approximated to the zincblende sp 3 hybridization way of bonding, which does not accurately hold for several of the materials studied in this work.In most instances, this approach is suitable for accurately describing the valence band.However, using the same method to represent the conduction band states may lead to inaccuracies.
Since the higher energy bands also contain more delocalized orbitals, they cannot be reproduced successfully.An additional s orbital from a higher electronic shell, or d orbitals, is often required to describe the chemical bonding better.Due to the fact that 3d orbitals are spherically asymmetrical, their interaction with the 3s and 3p states is not trivial.Hence, some deviations re-garding the position of the conduction band valleys may appear.Furthermore, this approach does not take into consideration the electronic spin that interacts with the electronic orbitals' motion which results in band degeneration removal. [16]n example of such analysis can be demonstrated by examining the case of bulk AlAs.
In Figure 6, we observed that the unstrained state of the bulk AlAs yielded the highest band gap, while the compressed and tensed structures resulted in lower values of band gaps.In direct band gap and zinc-blende semiconductors such as AlAs, the valence band edge comprises mainly bonding p orbitals, while the conduction band edge is primarily s orbitals.When a biaxial tensile strain is applied in the plane of (001), where the Al atoms are positioned, the Al atoms are taken apart with increased distance between them, and the As atoms bonded to them on the next atomic layer are stretched out in the same manner.As a result, the bond angles parallel to the strained plane increase, while the angles beside them decrease to retain 180°of a straight line, distorting the tetrahedra bonding.This manner is illustrated in Figure 7.
This distortion reduces the overlapping between the LCAO (linear combination of atomic orbitals).The band gap modifications are related to the coupling between the atomic orbitals or the overlap integrals of the orbitals wave function.Zinc-blende lattice arrangement is characterized by the tetrahedral configuration of the chemical bond between the crystal atoms that usually results from sp 3 orbital hybridization.Therefore, when a tensile or compression strain is applied, the bonds rotate toward the plane or out of the plane alternatively, decreasing the coupling between the atomic orbitals; thus, the separation of bonding and antibonding states is reduced and yields a narrower band gap.
A further discussion is brought here to elucidate the connection between the local potential peaks' parameters and the simulated one-electron transmission coefficient.
Throughout the bulk material analysis, each system was described by a three-curve plot of the local potential curves for each strain mode, as shown in    be seen in Figure 5.We aimed to correlate the electrostatic potential parameters and the calculated transmission coefficients.Our analysis was grounded in chemical principles related to the motion of electrons under potential barriers.For example, we anticipated that higher barriers would restrict electron motion, resulting in a greater distance between global minima and maxima.Furthermore, the barriers' width may affect the transport probabilities if it is sufficiently thin with respect to the Gaussian wave function spatial spread, .When the electron wave function widths were around 0.7-0.5 Å, the narrower the potential peaks, the enhanced the electron transport.Finally, it was also reasonable that whenever the electron was provided with higher kinetic energy, it overcame the barriers more easily and favored being transmitted.We tested these trends via the correlation procedure given in Figure 8, by determining the relative weight, the barriers' height, the barriers' width, and the electron initial energy contribute to the electron conduction.Hence, the correlation pa-  rameter in Equation ( 13) is constructed with each parameter's influence on electron transport.
After a repetitive trial and error process, guided by the analytical expression for a single barrier potential transmission coefficient, [17] we did not manage to find a satisfactory and clear correlation between the potential peaks' parameters and the transmission coefficients (reaching a ∼ 0.4 R 2 ).The comparison to the rectangular potential barrier addresses only a limited description of our examined fluctuated potential.In addition, the potential curves exhibit complicated patterns, making it even more challenging to model such potential behavior in terms of the curve heights and widths.Nonetheless, we constructed a general correlation parameter in which the expected relationships between the potential curve parameters are taken part where E is the initial energy of an electron, V 0 is the potential height from the minima to the maxima points, and d is the approximated peaks' width., , and  are the exponent values acting on each of the above, respectively.The power of , , and  determines the relative weight each parameter has on the correlation parameter.In Figure 8, the fitted values for each were as the following (they are unitless)  = 0.5,  = 0.25, and  = 0.25 There is no evident physical meaning for the magnitudes of the exponents but rather a result of the numerical fitting.When examining the correlation in Figure 8, one can see an apparent trend that holds well for the entire tested scope, except for two clusters of sample points that deviate from the trend line.Meaning, that even if the linearity is not entirely kept, the observation on the correlation curve tells us that, generally, the proportionalities that have been suggested are the right ones, and further mathematical investigation is needed to refine this analysis.

Heterostructure Calculations
We tried to reproduce the same changes in the band gap influenced by strain for the bulk structures when an interface is formed (see Table S7.0.0 in Supporting Information for more detailed results).For each heterostructure case study, we examined the local energy gap modifications with respect to the distance from the interface atomic layer by observing the projected DOS.An example of GaAs/AlAs heterostructure analysis can be observed in Figure 9.The same comment occurred for three different heterostructures of the same composition, each constructed from different characteristic lattice constants.Each heterostructure had three other supercells built from the longer lattice parameter among the two materials, the smaller one, and the averaged lattice parameter among the two materials.
When a longer lattice parameter is enforced on a semiconductor with a shorter one, it experiences biaxial tensile stress parallel to the interface plane.Similarly, when a semiconductor is constrained into a smaller lattice parameter that characterizes the lattice geometry, it experiences biaxial compression stress, again, parallel to the plane of the interface, which imposes the constraint.We anticipated to find a gradual strain relaxation as the atomic layer gets far away from the interface position.In contrast, the nearest atomic plane to the interface will experience the most considerable effects stemming from the interface presence.
Next, in order to find out how well the interfaces allow electronic transport, we investigated the transmission coefficient and the probability of an electron passing through the present interface.We approached this by integrating the electron flux over the entire simulation time, as in formulas (1) and ( 2).Each examined system was obtained from our heterostructure systems, as the primary input for the simulation was the electrostatic potential vector.
Each simulation was preceded by a series of numeric convergence tests to minimize the errors.Then, the supercell's local potential was extended to allow the wave function to propagate in space without reaching the edges, avoiding self-interactions and interference that might affect the results.Afterward, the initial wave function was centered near the interface position (about 3-4 Å distance from the interface) and had the highest probability of finding the electron.With the initialization of the simulation, the wave function spread in space, and as the simulation time went on, some probability amplitudes were passing the interface.
After the simulation was completed, the transmission coefficient for each heterostructure could be calculated.Simulating the opposite direction was applied, and the heterostructures were flipped to initialize the wave function at the other material region.This habit allowed us to double the obtained data for each system.
In Figure 10, we present three frames of the generated simulation animation of the propagating wave function as it passes through the interface.We can see in Figure 10 Eventually, all the transmission coefficients were calculated to understand composition-function relations, and the interfaces' parameters were gathered to correlate with the electron transmission probability.Each local potential curve was parameterized in terms of its peaks' characteristics, as they are shown in Figure 11.Our focus was on elucidating the connection between potential-related variables and the electronic conductivity of the heterostructure by modeling the interfacial contact between the two semiconductors.
In Figure 11, we demonstrate some meaningful curve characteristic parameters.First, when examining the left-hand side of the interface, E 3 is the difference between the left-side averaged height of the maximum and minimum peaks, E 4 is the difference between the left side absolute maximum value with averaged maximum peaks height, d( 001) is the d-spacing of the relevant crystal lattice plane, ϕ is the difference between the conduction band edge and the left side absolute maximum value,  is the interface region width, Δ is the difference between the left side averaged minimum peaks height and the height of the last minimum peak at the interface from its left, representing bend bending, E 2 is the difference between the right side averaged height of the maximum and minimum peaks, and E 1 is the difference between the left side averaged maximum peak height and the right side averaged maximum peak height.
To develop a fundamental understanding, we employed the basic quantum mechanics framework for the transmission coefficient of electron motion in a step-like potential.We reached good correlations, as shown in Figure 12 for the direct propagation and Figure 13 for the electron coming to the interface in the opposite direction of propagation (namely, the wave function is initialized at the other material, then propagates until it crosses the interface from the other side) ) The direct step-like correlation parameter takes the form of Equation ( 14), taken from ref. [17] and it is  1 , where Θ is the transmission coefficient.E 3 is the difference between the average of the maximum points' height and the average of the minimum points' height for the left-hand material, and E 2 is the difference between the average of the maximum points' height and the average of the minimum points' height of the right-hand material. is the energy difference between the average of the maximum points' height and the conduction band minima energy, corresponding to the initial total energy the electron is provided.The same quantities also appear in Equation ( 15) for the opposite direction correlation parameter  2 [17]    Both tests show good correlations with the quantum mechanical models of the step-like potentials, encouraging us to believe that modeling the shape of the interface is the first step to assess the electronic conductivity of more complex interface features.In addition, by accepting the relations between the interface parameters, one can predict the electron transmission coefficient in future works done for other systems.Accordingly, it alleviates heterostructure-based systems' design by suggesting optimal pairing of semiconductors.
The correlation expressions prove that the average height differences between the potential curves around the interface, e.g., E 3 and E 2 (see Figure 11), have the most dominant impact on the transmission coefficients.In contrast, it reveals that the interface potential width does not play a significant role in those cases.
We also investigated the propagation of the wave function in the opposite direction to validate the proposed method for linking the electrostatic potential parameters and the transmission coefficient at the interface.The resulting correlation was also satisfactory, with a high R 2 value of ≈0.92.However, the correlation quality was slightly poorer than that of the direct propagation method.This could be due to the step-like potential model used to derive the correlation parameter, which may deviate from the actual oscillating nature of the electrostatic potential near the interface.Hence, the deviation may be a bit larger for the incoming left electron with the lower energy close to the barrier.
Our numerical investigation of the electrostatic potential at the interfaces has been validated by the good agreement between the correlated models and the calculated transmission coefficients.This confirms the effectiveness of our approach for studying heterostructures.One of the main advantages of this method is its low computational cost, which creates new opportunities for materials research.

Conclusions and Summary
In the present work, some binary heterostructures composed of two typical semiconductors were under investigation to understand trends induced by lattice-mismatched interfaces.The local DOS calculations were carried out to obtain the band gap values of each group of atoms at different distances from the interface position.A better understanding of the trends has been achieved by reproducing the same calculations for the bulk material systems, thus having a degree of comparison.The main output from this part of the work was that the central regions within the heterostructure systems, far from the interfaces, retained the band gap changing trends similar to what we observed in the bulk systems.The gradual changes in the band gap along the interface area suggest that structural changes induced by strain may impact the modifications of densities at the band gap energy range.Alternatively, the introduction of chemical bonds to an additional material results in the alteration of atomic orbitals that form the energy bands of the crystal.Consequently, the band gaps of the new supercell differ from those of the bulk materials that constitute it.
Numerical simulations followed the exploration of the heterostructure systems for charge transport through the interface.Accomplishing such a numerical approach was done by evolving in time the wave pocket as a solution for the one-electron nonrelativistic and time-dependent Schrodinger equation.The time propagation of the initial wave function was made by approximating the propagator terms acted on the wave function using the split-step approximation.The transmission probability of an electron passing through the interface was calculated for each system on both sides of the interface.After parameterizing those, we successfully correlated the transmission coefficients of the heterostructure with the descriptor parameters for the near interface potential.Taking a fundamental quantum mechanics model of a step-like potential and adapting the actual parameters of the interface potential to have the same functional relations as the step-like potential, we reached a good fitting between the logarithm of the transmission coefficients and the fitting model quantity.
The main conclusions of our work are listed below: 1) The band gap values exhibit changes that resembles the values of the bulk structures as we move further away from the interface, as anticipated.2) No clear rule allows predicting a material response to biaxial strain.Each semiconductor has its own characteristics that determine how its electronic structure changes under imposed stress.3) Compressed bulk structures showed enhanced transmission coefficients in most cases.Deviations are attributed to the contribution of the supplied initial energy of the electron.4) Each interface has a unique nature that can be modeled and then correlated to predict electron transfer probability through the interface.We saw from the correlation with the step-like potential pattern that the most significant measures of the interface electrostatic potential were differences between the averaged heights of the maximum and minimum peaks of both sides of the interface.We also observed that the relative energy difference between the initial total energy of the electron and the highest potential energetic position significantly affected the interface transmission coefficient ( in Figure 11).This work demonstrated how materials studies could be made using numerical tools based on a relatively simple entity achieved from DFT calculation, as we did with the electrostatic potential vector of the material.We showed that we could fit the transmission coefficient of the electron passing through the material from the electrostatic potential vector.By altering the electrostatic potential of the material numerically, one can expand the range of possibilities of materials phenomena that can be tackled operating at a low computational cost.

Figure 2 .
Figure 2. 2D projection of AlAs-GaAs heterostructure exhibiting an interface between two different materials.

Figure 3 .
Figure 3. Calculated band-gap values for strained and unstrained values.The bulk materials with zero band-gap have no bars.

Figure 4 .
The relative position of the electrostatic potential curve regarding the atoms' positions can

Figure 4 .
Figure 4. Local potential curve with respect to the 2D projection of the AlAs atoms' arrangement.

Figure 5 .
Figure 5. Bulk AlAs local-potential curves.Three curves represent three strain profiles applied.

Figure 7 .
Figure 7. Illustration of the stretching mode in which two adjacent atomic layers of AlAs experience as a result of biaxial tensile strain.

Figure 8 .
Figure 8. Linear correlation between the calculated transmitted coefficients and the suggested correlation parameter, .

Figure 9 .
Figure 9. GaAs/AlAs heterostructure constructed from the averaged lattice constant of the two semiconductors.From above, the supercell atomic representation was drawn to provide a reference corresponding to the matching DOS plots below.From the left side of the interface is the GaAs region and from the right side is the AlAs region.Ga atoms are in red, Al atoms are in blue, and As atoms are in green.

Figure 10 .
Figure 10.Three simulation frames, taken from: a) initial time t = 0, b) t ≈ 48 s, demonstrating first probability current transfer, c) t ≈ 360 s, representing the spreading probability throughout space.
(a) the initial state of the wave function when it is highly localized on initialization, centered at a specific position.The next frame in Figure 10(b) demonstrates how the probability current represented by the squared absolute of the electron wave function passes the interface layer marked by a straight blue line.The last frame in Figure 10(c) shows how the probability of finding the electron spreads throughout the entire system scope for sufficient long durations.

Figure 12 .
Figure 12. a) Correlation curve for directly propagating the one-electron wave function via the interface.Displaying the connection between the logarithm of the transmission coefficient Θ and the correlation parameter  1 .b) Schematic illustration of the direct step like potential form.

Figure 13 .
Figure 13.a) Correlation curve for the opposite direction propagating oneelectron wave function via the interface.Displaying the connection between the logarithm of the transmission coefficient Θ and the correlation parameter  2 .b) Schematic illustration of the opposite direction step like potential form.

Table 2 .
Bulk structures computational details.