Directly Controlling the Transport Properties of All-Nitride Josephson Junctions by N-Vacancy Defects

All-nitride Josephson junctions are being actively explored for applications in superconducting quantum chips because of their unique advantages including their antioxidant chemical stability and high crystal quality. However, the theoretical research on their microstructure mechanism that determines transport properties is still absent, especially on the defects. In this paper, we apply the first principles and non-equilibrium Green’s function to calculate the electrical transport characteristics of the yellow preset model. It is first revealed that the N-vacancy defects play a crucial role in determining the conductivity of the NbN-based Josephson junctions, and demonstrate the importance for the uniformity of vacancy distribution. It is found that the uniform number of vacancies can effectively increase the conductance of Josephson junction, but the position distribution of vacancies has little effect on the conductance. The work clarifies the effect of the N-vacancy defects on the conductivity of the NbN-based Josephson junctions, which offers useful guidance for understanding the microscope mechanism of the NbN-based Josephson junction, thus showing a great prospect in the improvement of the yield of superconducting quantum chips in the future.


Introduction
Superconducting quantum chips have shown the unique application potential in superconducting quantum computers. The motivation in developing superconducting quantum chips comes from the increase of their decoherence time by the improvement in circuit design [1], new material choices [2][3][4][5][6][7][8], and the optimization of fabrication process [9,10]. As an important component of superconducting quantum chips, the choice of materials for Josephson junctions (JJs) affects the performance. In particular, epitaxially grown all-nitride Josephson junctions exhibit great potential to solve material-related problems, such as reducing microscopic two-level fluctuations because of their antioxidant chemical stability and high crystal quality [9]. Among them, niobium nitride (NbN)/aluminum nitride (AlN)/NbN JJ has been widely studied [3,4,11,12]. AlN, as the insulating layer, has the characteristics of a wide band gap and high resisitivity, while NbN, as the electrode material, has a high superconducting gap (~5.2 meV) and critical temperature (Tc = 17K), which suppresses the excitation of quasiparticles [9,13,14].
However, the properties of this all-nitride tunnel junction strongly depend on the crystal characteristics of the NbN and the interface around the tunnel barrier. Therefore, the NbN layers with a single crystal structure and a clean interface with the tunnel barrier are the necessary conditions for the development of high-quality all-nitride JJs [15]. To date, good experimental results have been achieved by optimizing the substrates and process conditions [9,12,[14][15][16][17][18][19][20]. For example, through increasing the buffer layers between the junction and the substrate to reduce the dielectric loss, the energy relaxation time T1 of the all-nitride junctions reaches 16 microseconds (ms), and the phase relaxation time T2 reaches 22 ms [13,21]. However, the results still fall short of expectations.
Recently, great efforts have been focused on exploring the simulation for the microstructure of all-nitride Josephson junctions owing to the close relationship between the physical properties and the electronic structure of the junction [22], and expected to help improve their performance. For example, John et al. [23] used density functional theory to study the microstructure and electrical properties of all-nitrogen interlayers, and the density of states revealed the causes of metallicity and rigidity in the interlayer. Zhang et al. [24] studied the growth mechanism of AlN by molecular dynamics. The effects of different N:Al fluxes ratios, temperatures, and stresses on the crystallinity of AlN were analyzed, and the environmental conditions of the best crystallization quality and the location of defects were obtained. However, the electrical characteristics of the microstructure in all-nitride JJs have not been described qualitatively/semi-quantitatively by computer simulation, especially with defects. Therefore, to make better use of the inherent advantages about all-nitride JJs, it is worth conducting further research.
In this paper, the NbN/AlN/NbN device models are created from the atomic level. As we know, it is difficult to fabricate single-crystal AlN thin films in the manufacture, which are primarily polycrystalline or amorphous, and have certain internal defects [25]. Among these defects, the N-vacancy is the most common defect type [25,26]. By studying the microstructure of the all-nitride JJs, the influence of N-vacancy defects on the electrical properties of the junction is explored. The critical factors that affect the conductivity of the models are illustrated from the two aspects, which are the uniformity and position of vacancies distribution in each atom layer. The results show that the uniformity of vacancies has a greater influence on the conductivity of the model. The study, in directly controlling the transport properties of the Josephson junction by N-vacancies, can guide the fabrication of ideal JJs and improve the yield of quantum chips.

Designing Defect Models with N Vacancies in the Barrier
In this study, we constructed a three-layer structure model of NbN/AlN/NbN, with NbN as the electrode and AlN as the central region. Considering that the lattice mismatch between the electrode and the center region should be as low as possible during modeling, three models, as shown in Figure A1, corresponding to different lattice constants, were selected for comparison. Model 3 (in Figure 1a), which has the lowest lattice mismatch degree and the highest binding energy (as shown in Table A1), was selected as the reference model, and its influence on electrical properties was studied by creating N vacancies.
To increase the comparability of the experiments, all defect models had the same total vacancy numbers and the same vacancy element type. The total number of vacancies was 8, and the vacancy element was nitrogen. Based on the reference model 3, all N vacancies were distributed at the interface layer to construct a defect model 3b, as shown in Figure 1b. The defect models 3c-3f, as shown in Figure 1c-f, respectively, were constructed by setting different N vacancies for each atomic layer of AlN in the reference model. Under the condition of ensuring that the number of N vacancies in each atomic layer of AlN in the defect model 3c was consistent, the position of N vacancies in the atomic layer was adjusted, and the models 3c-1com and 3c-com2 were constructed, as shown in Figure 1g,k. The three models were also called the group 3c. There were five groups of models, and we controlled the number of vacancies and changed the positions of vacancies in the same group model, to study the influence of the number and positions of vacancies on the electrical properties, respectively. The N atomic vacancies of model 3b were all on the upper interface, and the vacancy distribution of model 3c-3f was set in the different positions, as shown in Figure 1. The vertically aligned models were the same group comparison models, while the horizontally aligned models were different groups of models in the figure. (c,g,k) the vacancies distribution of the group 3c, including model 3c, 3c-com1, and 3c-com2. The number of atoms in each layer is 0, 3, 2, 3, respectively; (d,h,l) the vacancies distribution of the group 3d, including model 3d, 3d-com1, and 3d-com2. The number of vacancies of atoms in each layer is 1, 2, 3, and 2, respectively; (e,i,m) the vacancies distribution of the group 3e, including model 3e, 3e-com1, and 3e-com2. The number of vacancies of atoms in each layer is 2, 1, 1, and 4, respectively; (f,j,n) the vacancies distribution of the group 3f, including model 3f, 3f-com1, and 3f-com2. The number of vacancies of atoms in each layer is 2.

Interface Configuration
Miller Bravais indices are an essential parameter affecting interface architecture. The crystal properties of the substrate largely determine the structure of the deposited material. However, there is an interesting phenomenon that AlN shows a (001) orientation and high resistance in the lower layer of NbN sputtering deposition, but the reverse electrode NbN deposited on the AlN barrier layer has a polycrystalline structure [11]. The study indicates that the thickness of the barrier, such as AlN, affects the crystal orientation of the NbN layer at the top of the barrier [19]. Therefore, the thickness of AlN selected in this experiment is within 1 nm to ensure that the upper NbN still has a single crystal structure. The epitaxial growth of NbN on single crystal Si (200) substrate is mainly characterized by (111) and (100) cubic structure, and the AlN film is characterized by (111) hexagonal structure [20,27,28]. There is no primary orientation relationship between NbN and AlN grown under different process conditions. To better compare the orientation in the model, the lattice constant, Miller indices, and other parameters were selected as shown in Table A1.
The interface in the vacancy model was optimized after relaxation, when the model was constructed. The convergence criteria of total energy and force were consistent with the defect-free model (as shown in Appendix A.2). It was also necessary to carry out quantum transport calculation for the defect models, and the parameter setting of the calculation was consistent with that of the defect-free models (as shown in Appendix A.2), for which the overall energy of the model tends to be the lowest, and the overall structure of the model tends to be stable.

Model Calculation
After the device models were built, the quantum transport properties were calculated using the non-equilibrium Green's function (NEGF)-DFT quantum transport program Nanodcal. The density matrix ρ can be given by the intrinsic state ψ i of the KS equation: where β is the electron temperature, µ is the chemical potential, f EQ (E − µ; β) is the Fermi Dirac distribution function, expressed as: where E is the electron energy. To achieve self-consistency, the density matrix is iteratively solved by calculating the Green's function, which can be expressed as: where G R and G A indicate the retarded and advanced Green's functions, respectively. S is the overlap matrix caused by non-orthogonal orbits (if the orbits are orthogonal, S is the identity matrix), and H is the Hamiltonian matrix of the system. Σ R ,Σ A represents retarded and advanced self-energy respectively, and G A can be calculated by Considering an open system with several leads, the current can be calculated by the Landauer formula: where f l , f r is the distribution function of the electrons in the left or right of the electrodes, respectively. The factor of 2 comes from the spin degeneracy.
In addition, the conductance can be defined as where G 0 = (2e 2 )/h is the quantum conductance, h is Planck's constant, e is the charge of a single electron, and f l ( ) − f r ( ) is the difference between the Fermi distribution of the left and right electrodes, V r − V l is the difference in bias voltage applied to the left and right electrodes, and T( ) is the transmission coefficient under a certain energy. When the applied bias voltage approaches zero, the limit value obtained is the conductivity value in the equilibrium state. It can be seen from the formula that the conductivity is closely related to the applied bias voltage and the transmission coefficient. This theory, it should be noted, is based on the conventional ballistic transport theory and calculates the electrical properties of devices at 0 K . The tunneling effect of the Cooper pair at superconducting temperature and the influence of electroacoustic interactions are not considered, so the conductance generated by the tunneling effect of the JJs in the actual process is different. However, from the qualitative point of view, the calculation results have reference value, and in the study of the impact of defects on electrical properties, the first-principles simulation model can be more convenient to guide the process.

Electrical Properties of Defect Models
After optimization, the layer spacing at the interface was changed, as shown in Table 1. At the upper interface, the bond lengths between Al in AlN and N in NbN decreased from 2.09 Å to 1.99 Å. The same happened at the lower interface, where atoms converged near the vacancy [29]. The bond length between the N atom in the barriers and the Nb atom in buffers increased from 2.25 Å to 2.29 Å at the upper interface. Still, there was no obvious pattern at the lower interface. Due to the vacancies, the Al atoms in the interface layer attracted the N atoms in the NbN, and the distance between Al and N became smaller at the upper and lower interfaces. On the contrary, the N atom in AlN repulsed with the Nb atom in NbN, resulting in the bond between N and Nb getting longer, and the lower interface also showed such a rule. Table 1. Comparison results of the spacing between the layers at the interface before and after relaxation.

Model The Upper Interface
The Lower Interface Defects, including the quantity uniformity and position distribution, were observed by increasing N vacancies, and the effects of vacancies were analyzed on the electrical properties of the models. Since the conductance of the three defect-free models is in the same order of magnitude and the difference is slight (as shown in Figure A2), the vacancies were added to study the factors affecting the junction performance. To qualitatively characterize the differences in the conductivity of different vacancy models in the barrier layer, we calculated the conductivity under zero bias voltage, as shown in Figure 2a. It is found that vacancies have a significant influence on the conductance of the models. In Figure 2a, the same color represents the same group of models. It can be seen that vacancies lead to an increase in the conductance of the models. The conductance of model 3f is the largest, while that of model 3b is the smallest. The research also shows that, due to the existence of vacancies, additional electronic states will be formed in the insulator gap [30,31], which may cause an increase in conductance. The data are consistent with the variation of the transmission coefficient shown in Figure 2b, which reflects the transmission coefficient of different vacancy distribution models near the Fermi level in the equilibrium state. It can be seen from the figure that, near the Fermi level, the transmission coefficient of model 3f is the strongest, while that of model 3b is the weakest. Previous studies on oxide vacancies have shown that a single oxygen vacancy or oxygen vacancy chain would lead to delocalization, resulting in a wide transport range and high transmission coefficient [32,33], which was similar to our result.

Electron Density
Local density of states in all of vacancy models near the Fermi level were calculated, as shown in Figure 3, which is consistent with the vacancy distribution shown in Figure 1. In the defect models, the distribution of vacancies in each atomic layer has a great influence on the density of electron states in the barrier layers, as shown in the horizontal view of Figure 3. In each group of models, the number of vacancies in each layer is the same, but the positions are different, resulting in different degrees of change in density of states in the barrier, as shown in Figure 3c,g,k. In general, the density of states in the barrier is the highest in the group 3f, while there is basically no electron state in the barrier layer of the model 3b. It indicates that the possibility of electrons passing through the group 3f is the largest, so the conductivity will be better. This can also explain the transmission coefficient of model 3f being the highest, while that of model 3b is the lowest. On the whole, vacancies lead to an increase in the density of states. Not only the density of states near the Fermi level becomes prominent, but the peak values of the valence band and conduction band increase in varying degrees. It can be seen from Figure 3 that the higher the density of states in the barrier, the higher the corresponding transmission coefficient, which is also similar to the previous research results [29], thus leading to the increase of conductance.
The electronic structure of AlN is closely related to its crystal structure. The results show that the crystal system of zinc blende AlN in metastable phase is hexagonal, which belongs to an indirect band gap semiconductor, while the AlN with zinc blende structure belongs to a direct band gap semiconductor [34]. Model 3b and 3d were selected for the electron density comparison with the defect-free model. The vacancy distribution and conductance of these three models were different. Inspecting the total density of states in the device (as shown in Figure 4a,c,e, respectively), the density of states increases in both defect model 3b and 3d. The local density of states in Figure 3a,b,d can be related to the total density states shown in Figure 4a,c,e, respectively, which reveals that the device models of three structures are metallic. There is no state in the barrier layer at the Fermi surface (as shown in Figure 3a), while the metallicity of the devices in Figure 4a at the Fermi surface should be attributed to the role of NbN, which is consistent with the results of previous studies [23,33]. When the vacancies are made, AlN in the device model generates unsaturated bonds. The conduction band and valence band appear degeneracy and metallicity, increasing the density of states near the Fermi surface (as shown in Figure 4b,c).  It is found that the vacancies lead to the change of the electron density in the vacancy areas. From the diagrams (as shown in Figure 4b), the delta charge of the defect-free model in the bonding region is relatively large, evenly arranged, and has clear boundaries. Due to the absence of N atoms, the bond between Al and N breaks and the dangling bonds appear around Al atoms. The delta charge in the bonding region near the vacancy becomes smaller, and the local deformation occurs (as shown in Figure 3d,f). The delta charge at the N vacancy becomes smaller, indicating that the electrons gather on the Al ion near the N vacancy, causing the change of differential electron density, which is also responsible for enhancing the conductance [32].
Owing to the defects, fluctuations in local density and stoichiometry will lead to local transmission channels. Previous studies have shown that conduction mainly passes through metal channels in highly hypoxic structures [35]. The vacancies of the O atom help the electron to tunnel through the channel of the metal ion [32]. In our study, the absence of N atoms leads to an increase in local density of states near the Fermi level and the formation of local transport channels. The change of delta charge near the N vacancies indicates the aggregation orientation of electrons, which aggregate towards metallic Al atoms to achieve the channel tunneling of electrons through Al ions.

Analyzing Vacancies Distribution
To further analyze the influence of vacancies on the electrical properties of the models, we have made statistics on the vacancy distribution all models. The number distribution of N vacancies is counted in each barrier layer of each model, and the variance of the number distribution is calculated. Compared with the corresponding conductance, the average conductance growth amplitude of the vacancy model with that of the defect-free model is worked out as shown in Table 2. It is clear from the table that the uniformity of vacancy distribution affects the conductance of the model to a certain extent. When the variance of vacancy distribution is minor, that is, the number of vacancies is evenly distributed, and the conductance is relatively large. On the contrary, the conductance is small, as shown in Figure 5a. From the quantity distribution of vacancies, group 3f with the minimum variance has the most considerable conductance, while model 3b with the maximum variance has the smallest conductance. This is due to the uniform distribution of vacancies, which reduces the electron tunneling barrier, increasing the electron tunneling probability and the conductance correspondingly.  It can be seen that the positions of vacancies have little influence on conductance. By comparing models 3c, 3d, 3e, and 3f with their comparative experiments in the same group (as shown in Table 2), it is found that the conductance does not vary significantly, when the number of N vacancies in each layer remains unchanged and the positions are changed. To compare the influence of the position distribution of N vacancies on the conductance, the bond length between each Al atom of the barrier layer and N vacancies in each group of classification models was averaged. Then, the average bond length between Al atoms of the barrier layer and N vacancies in the whole model is obtained (as shown in Figure 5b). In general, the conductance is inversely proportional to the average bond length from the figure. By comparing the mean values of four groups, the mean values of group 3d and 3f become smaller, that is, the more uniform the vacancy distribution is, the average bond length tends to decrease. Meanwhile, the larger the conductance, the larger the average bond length fluctuation.
A large number of previous studies have shown that vacancy defects in oxides have a significant effect on the electrical properties of materials [32,33,[35][36][37][38][39]. In our study, vacancy defects enhance the electrical properties of the device model, indicating that nitrides and oxides have similar phenomenal-level defect characteristics. Compared with previous studies, the vacancy distribution was a single vacancy or continuous vacancy chain [32,33], and the distribution was very regular. Our study reveals the effect of the pseudorandom distribution of N vacancies on the electrical properties of materials, indicating that the uniformity of the vacancy distribution in the barrier layer has a significant influence on the electrical properties. At the same time, although AlN is a single crystal material and its structure is relatively simple, its microstructure is still extremely complex in the actual process. From the analysis of the distribution of vacancy positions, we can also see the complexity of its microscopic characteristics. The study of the microstructure of AlN is very worthwhile in reducing the two-level defects of superconducting materials and improving the performance of Josephson junctions.

Conclusions
The microstructure of the barrier layer in JJs plays an essential role in the performance of the junction. After calculating the conductance of three defect-free models, the model with the highest binding energy and the lowest mismatch rate is taken as the reference model. Considering the effect of vacancies on the junction performance, we characterize the conductance of different defect models at zero bias. Combining the DFT and NEGF, the results show that the vacancies increase the electrical transport properties of the junctions. When the total number of vacancies is controlled and the number of N vacancies is increased in each layer, the conductance is greatly affected. When the number of N vacancies in each layer is fixed and the positions of the vacancies are changed, the conductance is not obviously affected. By calculating the average bond length between each Al atom and N vacancies in the barrier layer and combining with the variance distribution of N vacancies, it is found that the uniformity of vacancy distribution has a great influence on the conductivity of the junction. We mainly guarantee the accuracy of calculation from the perspective of first principles. Although the model has undergone various relaxations and self-consistency of functional calculation, it will still be different from the actual fabrication of JJs. In addition, since the ballistic transport theory is used in the calculation and the electrical properties at 0K are calculated, the tunneling effect of Cooper pairs at the superconducting temperature and the influence of phonons are not considered, and we can only qualitatively discuss the impact of N-vacancy defects on the electrical transport properties of Josephson junctions. However, defects can be inevitable in the process of fabrication. The in-depth understanding of the defect mechanism can help us provide a more precise idea and guidance in the manufactures.

Patents
The Chinese patent "Method for regulating the transport properties of NbN/AlN/NbN Josephson Junction by using N-vacancy defects" has been accepted.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Acknowledgments: This work was undertaken with resources from the National Supercomputing Center in Zhengzhou.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A Appendix A.1. System Models
In this study, the models were based on the three-layer structure of JJs, and the stacking order was NbN/AlN/NbN, as shown in Figure A1a. The xand y-directions of all models were periodic, and the z-direction was the direction of electron transport. During modeling, the mismatch rate between the electrode regions and the central region should be kept as low as possible. Meanwhile, considering calculation cost and speed, the model scale should not be too large. This study was the first attempt to construct the theoretical model of the all-nitride junctions, and the model were limited to about 200 atoms. The choice of model was shown in Table A1. In model 1, AlN, as the barrier layer, select the hexagonal wurtzite structure, as shown in Figure A1b, which had stability in natural conditions. In model 2 and model 3, AlN were both tetra sphalerite structures, as shown in Figure A1d,f, respectively. Although it is difficult to obtain under natural conditions, there are technological conditions that can prepare the cubic structure of sphalerite AlN [34]. Considering the mismatch rate of the lattice, these two structures were chosen (as shown in Table A1).  NbN was chosen as the electrode region, which was approximately considered to be infinitely long. After the Hamiltonian of the electrode was obtained by DFT, the infinite periodic system was truncated into a semi-infinite electrode. The device model was built with different AlN structures as the central region, as shown in Figure A1a,c,e, respectively.

Appendix A.2. Optimizing Models
The structure of the atoms was optimized by using DFT based on the projector augmented wave (PAW) method. During relaxation, the Vienna Ab-initio Simulation Package (VASP) [21,28,40] was used for calculation. Perdew-Burke-Ernzerhof (PBE) [41] was used for the exchange-correlation functional between electrons. PAW was used for pseudopotential, and the sampling value of the Brillouin region was the Gamma center strategy. The plane wave cut-off energy of the NbN was set to 400 eV, and the k-point was carried out according to the solid model. The two sides of the AlN were added vacuum layers along the z-axis, and the position of atoms was optimized using k-point sampling of 8*8*1 according to the two-dimensional material. The total energy convergence criterion was at 10 −4 eV per atom, and the maximum force was less than 0.05 eV/Å.
To obtain a more stable and reasonable device model, we combined the optimized NbN and AlN. The suitable distance between the barrier layer and the superconducting layer interface can be obtained by DFT total energy relaxation. Then, we obtained the optimum distance between AlN and NbN of different models. In the best model construction, vacuum regions of 7-8 Å were added to both sides of the device in the z-axis direction to make the electric potential transition to the vacuum level. The model with the lowest energy was selected to establish the initial crystal model, which was subjected to interface relaxation to optimize the interface between NbN and AlN.
Nanodcal [42] software was used for quantum transport calculation, and the exchangecorrelation functional was the same as that used in structural optimization. The convergence criteria for the Hamiltonian and density matrices was set to 10 −5 eV, and the Kohn-Sham wave function was expanded using a linear combination basis of atomic orbitals. The cut-off energy of real space grid was 100 Hatree. The k-point was set to 5*5*100 and 5*5*1 for self-consistent calculation of the electrode part and the central region, respectively. The three defect-free models are shown in Figure A1a,c,e, respectively, and different atomic structures of the central region AlN are shown in Figure A1b,d,f, respectively, on the x-y plane.

Appendix A.3. Electrical Properties of Defect-Free Models
First, based on Nanodcal software, we calculated the conductivity of the three defectfree models without any bias (i.e., zero bias), as shown in Figure A2a. Because the barrier layer in the model is very thin (<1 nm), it becomes very sensitive to the effects of the interface. The differences in conductance among the different models are shown in the figure. The conductance of the three models is approximated in the same order of magnitude, and the I-V curve in Figure A2b shows the same relationship. Model 1 has the largest conductance and corresponding current, while model 2 has the smallest conductance and corresponding current. On the basis of zero bias, we calculated the relation between current and voltage in the range of 0.001-1.2 V. When the bias voltage is low in the sample, the current shows linear growth. When the bias voltage increases to a certain threshold, the current increases exponentially, which is consistent with the experimental results [28]. The current can be regarded as the integral of the transmission spectrum within a certain energy range. Figure A2c-e show the electron transmission spectra corresponding to the three models under different bias. The two ends of the black dashed lines in the figure correspond to the size of the bias window, and the zero energy point is set in the middle of the bias window. As can be seen in the figure, when the bias window is small, the transmission spectrum changes gently. While with the increase of the bias, the transmission coefficient increases, and a new transmission peak appears in the figure. This is consistent with the nonlinear growth of current shown in Figure A2b.