Analysis of tunneling probability in heavily doped 4H-SiC Schottky barrier diodes based on complex band structure considering barrier potential

To achieve a carbon-neutral society, further developments in semiconductor power devices are needed. Although silicon (Si) is now most widely used for power devices, the technology for Si power devices is relatively mature. ¹⁾ In this context, the use of wide-bandgap semiconductors ^2–6) has been considered. Among them, silicon carbide (SiC), ⁷⁾ especially 4H-SiC, is attracting attention due to both its excellent physical properties and mature device fabrication technologies. ⁸⁾ It is now being studied not only as a next-generation power device material, ^9–15) but also for applications other than the power devices. ^16–19)

Further development of 4H-SiC devices requires the understanding of the properties of 4H-SiC. In particular, a better understanding of the tunneling effects in 4H-SiC is needed to obtain 4H-SiC ohmic contacts with lower resistance. Therefore, in this study, we consider heavily doped 4H-SiC Schottky barrier diodes (SBDs), which have a higher doping concentration than those practically used, to make the tunneling currents prominent, and perform a theoretical analysis of the tunneling effects in SBDs. Hara et al. ²⁰⁾ reported that smaller values of the effective mass and the barrier height than the experimental values had to be used in order to reproduce the experimental tunneling currents in SBDs. In previous studies including Ref. 20, the parabolic band approximation (PBA) is widely used, i.e. the complex bands for the tunneling probability are assumed to be parabolic. However, complex bands are inherently non-parabolic, which means that the parabolic assumption can lead to be inaccuracy of the tunneling probability.

In our previous report, ²¹⁾ we calculated non-parabolic complex bands by the empirical pseudopotential method (EPM) and analyze the tunneling effects within the WKB approximation. ²²⁾ Though we found that considering the non-parabolicity of the complex bands is important when calculating the tunneling probability, there were still some concerns. One of them is the use of the WKB approximation when calculating the tunneling probability. The WKB approximation is widely used to calculate tunneling probabilities. However, the WKB approximation is derived assuming a slowly varying potential. ²³⁾ Therefore, its use for steeply varying potentials is controversial. ²⁴⁾ Since the barrier potential of a heavily doped 4H-SiC SBD as a function of position z, V_b(z), changes steeply especially when the reverse voltage is large, the validity of using the WKB approximation for the tunneling probability in the heavily doped 4H-SiC SBDs needs to be confirmed before calculating the reverse voltage dependence of the tunneling currents. Another is that V_b is not taken into account when calculating the complex bands; we calculated the tunneling probability using the electronic states without V_b. In this study, we investigate the validity of the previous method ²¹⁾ by calculating the EPM electronic states considering V_b in the Schr$\ddot{{\rm{o}}}$dinger equation and the tunneling probability without using the WKB approximation. We then compare the calculated tunneling probabilities with the previous and experimental results.

This paper is divided into four sections. In Sect. 2, we describe the calculation methods used in this study. In the extended abstract, ²⁵⁾ we calculated the electronic states considering V_b using a periodic superlattice model. In the present study, we extend the method and derive the tunneling probability by combining the electronic states of the superlattice structure and bulk 4H-SiC. In Sect. 3, we show and analyze the calculated results. Finally, a conclusion is given in Sect. 4.

We first introduce the models to calculate the tunneling probability without using the WKB approximation in Sect. 2.1. In Sect. 2.2, we describe the calculation method of the electronic states with V_b. Then, the calculation methods of the tunneling probability without and with the WKB approximation are explained in Sects. 2.3–2.5.

2.1. Model

In this section, we explain the details of the calculation models for calculating the tunneling probability, T(E_z ), as a function of the kinetic energy along the transport direction and the complex bands.

Figure 1(a) shows the band diagram of a SBD we consider in this study. The vertical axis represents the energy E and the horizontal axis represents the z direction. E is defined with reference to the conduction band minimum of 4H-SiC without V_b and the origin in the z direction is the Schottky interface. The barrier potential V_b(z) includes the image force effects ²⁶⁾ and is assumed to be ²⁷⁾

$$\begin{eqnarray}&&{V}_{{\rm{b}}}(z)=\displaystyle \frac{{e}^{2}{N}_{{\rm{d}}}}{2{\epsilon }_{{\rm{s}}}{\epsilon }_{0}}{\left(d-z\right)}^{2}-\displaystyle \frac{{e}^{2}}{16\pi {\epsilon }_{{\rm{s}}}{\epsilon }_{0}z}-e\left({V}_{{\rm{d}}}+{V}_{{\rm{r}}}\right),\end{eqnarray} \tag{ 1 }$$

where V_r is the reverse voltage, V_d the diffusion potential, N_d the donor density, d the length of the depletion layer, and _s the relative dielectric constant of 4H-SiC. The barrier height of a SBD including the image force lowering, Φ, is expressed as ${\rm{\Phi }}={V}_{{\rm{b}}}({z}_{\max })-{E}_{{\rm{F}}}^{{\rm{m}}},$ where ${z}_{\max }$ is the z-coordinate at which V_b(z) has its maximum value and ${E}_{{\rm{F}}}^{{\rm{m}}}$ is the Fermi level of the metal.

Zoom In Zoom Out Reset image size

The schematic diagram of the model to obtain T(E_z ) is shown in Fig. 1(b). We divide the system into three regions; the left bulk 4H-SiC region, the central Schottky barrier region with V_b(z), and the right bulk 4H-SiC region. We call the boundary between the left (right) bulk region and the central Schottky barrier region the left (right) boundary. The central region consists of several primitive unit cells of 4H-SiC along the z direction, which is the transport direction. Note that the barrier potential V_b(z) includes a constant region near each end of the central region as shown in Fig. 1. The reason for this is explained in Sect. 3.

When N_d and V_r are assumed to be 2.6 × 10¹⁸ cm⁻³ and 10 V, respectively, the tunneling currents at energies near the Fermi level of the metal are dominant in the whole tunneling current, and the tunneling length is at most 6 nm. Since the length of one primitive unit cell of 4H-SiC along the 〈0001〉 direction, L_cell, is about 1 nm, we can take the tunneling effects into account by cosidering the superlattice whose one period consists of 6 primitive unit cells. If the form of V_b(z) changes, the calculation can be performed in the same way by changing the period of the superlattice accordingly. We use the EPM using the plane wave expansion when calculating the complex band structures. We, therefore, consider a periodic superlattice model when calculating the electronic states with V_b(z) by periodically arranging the central Schottky barrier region of Fig. 1(b). The schematic diagram of the resulting superlattice model is shown in Fig. 1(c).

2.2. Complex band structure

The EPM electronic states ²⁸⁾ can be obtained by solving the following equation: ^29,30)

$$\begin{eqnarray}&&\left[-\frac{{{\hslash }}^{2}}{2{{bm}}_{0}}{{\rm{\nabla }}}^{2}+{V}_{\mathrm{ps}}({\boldsymbol{r}})+{V}_{{\rm{b}}}(z)-{E}_{n,{\boldsymbol{k}}}\right]{\varphi }_{n,{\boldsymbol{k}}}({\boldsymbol{r}})=0.\end{eqnarray} \tag{ 2 }$$

Here V_ps( r ) is a pseudopotential, V_b(z) the barrier potential, b a factor to take the non-local effects into account, φ( r ) the wavefunction, and E_{n, k} the energy eigenvalue corresponding to φ_{n, k} ( r ). Note that the subscript n, k represents the n-th smallest eigenvalue with a wave vector k . V_b(z) is ignored when calculating electronic states of bulk 4H-SiC. The potential V_ps( r ) is expressed as the sum of the pseudopotentials of each Si and C atom as follows:

$$\begin{eqnarray}\begin{array}{rcl}{V}_{{\rm{ps}}}({\boldsymbol{r}}) & = & \displaystyle \sum _{{\boldsymbol{G}}}\displaystyle \sum _{j=1}^{4}\left[\displaystyle \frac{\left.{V}^{{\rm{Si}}}({\boldsymbol{G}}\right)}{2}\exp (-i{\boldsymbol{G}}\cdot {{\boldsymbol{r}}}_{j}^{{\rm{Si}}})\right.\\ & & +\left.\displaystyle \frac{\left.{V}^{{\rm{C}}}({\boldsymbol{G}}\right)}{2}\exp (-i{\boldsymbol{G}}\cdot {{\boldsymbol{r}}}_{j}^{{\rm{C}}})\right]\exp (i{\boldsymbol{G}}\cdot {\boldsymbol{r}}),\end{array}\end{eqnarray} \tag{ 3 }$$

where G = ( G _xy , G_z ) = (G_x , G_y , G_z ) is the reciprocal lattice vector, V^Si( G ) (V^C( G )) the pseudopotential of Si (C) atom, and ${{\boldsymbol{r}}}_{j}^{\mathrm{Si}}$ (${{\boldsymbol{r}}}_{j}^{{\rm{C}}}$) (j = 1, 2, 3, 4) the vectors representing the positions of the Si (C) atoms in a primitive unit cell. The wavefunction φ_{n, k} ( r ) is expressed as

$$\begin{eqnarray}&&{\varphi }_{n,{\boldsymbol{k}}}({\boldsymbol{r}})=\displaystyle \sum _{{\boldsymbol{G}}}{f}_{n,{\boldsymbol{k}}}({\boldsymbol{G}})\exp (i({\boldsymbol{k}}+{\boldsymbol{G}})\cdot {\boldsymbol{r}}).\end{eqnarray} \tag{ 4 }$$

From Eqs. (3) and (4), Eq. (2) can be transformed to

$$\begin{eqnarray}\begin{array}{l}\left[\begin{array}{cc}0 & {\boldsymbol{I}}\\ -{{\boldsymbol{H}}}^{0} & -\displaystyle \frac{\sqrt{2{{bm}}_{0}}}{{\hslash }}{{\boldsymbol{H}}}^{1}\end{array}\right]\qquad \left[\begin{array}{c}{\boldsymbol{F}}\\ \displaystyle \frac{{\hslash }}{\sqrt{2{{bm}}_{0}}}{k}_{z}{\boldsymbol{F}}\end{array}\right]\\ \quad =\,\displaystyle \frac{{\hslash }}{\sqrt{2{{bm}}_{0}}}{k}_{z}\left[\begin{array}{c}{\boldsymbol{F}}\\ \displaystyle \frac{{\hslash }}{\sqrt{2{{bm}}_{0}}}{k}_{z}{\boldsymbol{F}}\end{array}\right],\end{array}\end{eqnarray} \tag{ 5 }$$

where I is an identity matrix and F is given by

$$\begin{eqnarray}{\boldsymbol{F}}={\left[\begin{array}{cccc}{f}_{n,{\boldsymbol{k}}}({{\boldsymbol{G}}}_{1}) & {f}_{n,{\boldsymbol{k}}}({{\boldsymbol{G}}}_{2}) & \cdots & {f}_{n,{\boldsymbol{k}}}({{\boldsymbol{G}}}_{p})\end{array}\right]}^{{\rm{T}}}.\end{eqnarray} \tag{ 6 }$$

Here p is the number of plane wave bases for φ_{n, k} ( r ). The (l, j) elements of H ⁰ and H ¹ are expressed as follows:

$$\begin{eqnarray}\begin{array}{rcl}{H}_{l,j}^{0} & = & \frac{{\hslash }^{2}}{2{{bm}}_{0}}\left[{\left({{\boldsymbol{k}}}_{{xy}}+{{\boldsymbol{G}}}_{\text{l},{xy}}\right)}^{2}+{G}_{l,z}^{2}-E\right]{\delta }_{{{\boldsymbol{G}}}_{l},{{\boldsymbol{G}}}_{j}}\\ & & +\displaystyle \sum _{j=1}^{4}\left[\frac{\left.{V}^{{\rm{Si}}}({\boldsymbol{G}}\right)}{2}\exp (-i{\boldsymbol{G}}\cdot {{\boldsymbol{r}}}_{j}^{{\rm{Si}}})\right.\\ & & \left.+\frac{\left.{V}^{{\rm{C}}}({\boldsymbol{G}}\right)}{2}\exp \left(-i{\boldsymbol{G}}\cdot {{\boldsymbol{r}}}_{j}^{{\rm{C}}}\right)\right]\\ & & +{\delta }_{{{\boldsymbol{G}}}_{l,{xy}},{{\boldsymbol{G}}}_{j,{xy}}}{\int }_{0}^{{L}_{{\rm{sl}}}}\frac{\left.{V}_{{\rm{b}}}(z\right)}{{L}_{{\rm{sl}}}}\exp \left[i\left({G}_{j,z}-{G}_{l,z}\right)\cdot z\right]{dz},\\ {H}_{l,j}^{1} & = & \frac{{\hslash }^{2}}{2{{bm}}_{0}}2{G}_{l,z}{\delta }_{{{\boldsymbol{G}}}_{l},{{\boldsymbol{G}}}_{j}},\end{array}\end{eqnarray} \tag{ 7 }$$

where k _xy = (k_x , k_y ) is the in-plane wave vector. By solving the eigenvalue problem of Eq. (5) for a given k _xy and E, we obtain a complex waevenumber k_z ≡ k + i κ (both k and κ are real), from which we can construct the complex band structures. Hereafter, from the wavefunction φ_{n, k} , we omit the band index n and add the subscript E to explicitly indicate that the state is at energy E, and we refer to it as φ _{k ,E} .

2.3. Scattering states

We calculate the scattering states by connecting the incident, transmitted, and reflected EPM wavefunctions in the model of Fig. 1(b). We first calculate the bulk wavefunction ${\varphi }_{{{\boldsymbol{k}}}^{{\rm{L}}},E}$ (${\varphi }_{{{\boldsymbol{k}}}^{{\rm{R}}},E}$) for the left (right) region; the wavefunction for the central region, ${\varphi }_{{{\boldsymbol{k}}}^{{\rm{C}}},E}$, is calculated within the superlattice model of Fig. 1(c). The scattering wavefunctions associated with an incident wave ${\varphi }_{{{\boldsymbol{k}}}^{\mathrm{inc}},E}^{\mathrm{inc}}$ having a real wave vector k ^inc from the left are then expanded by ${\varphi }_{{{\boldsymbol{k}}}^{\alpha },E}$ (α = L, C, R):

$$\begin{eqnarray}&&{\psi }^{{\rm{L}}}={\varphi }_{{{\boldsymbol{k}}}^{\mathrm{inc}},E}^{\mathrm{inc}}+\displaystyle \sum _{{k}_{z}^{{\rm{L}}}}{C}_{{k}_{z}^{{\rm{L}}}}{\varphi }_{{{\boldsymbol{k}}}^{{\rm{L}}},E},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\psi }^{{\rm{C}}}=\displaystyle \sum _{{k}_{z}^{{\rm{C}}}}{C}_{{k}_{z}^{{\rm{C}}}}{\varphi }_{{{\boldsymbol{k}}}^{{\rm{C}}},E},\,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\psi }^{{\rm{R}}}=\displaystyle \sum _{{k}_{z}^{{\rm{R}}}}{C}_{{k}_{z}^{{\rm{R}}}}{\varphi }_{{{\boldsymbol{k}}}^{{\rm{R}}},E},\,\end{eqnarray} \tag{ 10 }$$

where ${C}_{{k}_{z}}$ is a expansion coefficient. Now, considering the condition that ψ^L and ψ^C are continuous at the left boundary, the following equations have to be satisfied for any x and y,

$$\begin{eqnarray}&&{\left.{\psi }^{{\rm{L}}}\right|}_{z={z}_{{\rm{L}}}}={\left.{\psi }^{{\rm{C}}}\right|}_{z={z}_{{\rm{L}}}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial {\psi }^{{\rm{L}}}}{\partial z}\right|}_{z={z}_{{\rm{L}}}}={\left.\displaystyle \frac{\partial {\psi }^{{\rm{C}}}}{\partial z}\right|}_{z={z}_{{\rm{L}}}},\end{eqnarray} \tag{ 12 }$$

where z = z_L (=0) is the boundary between the left and the central region. By Fourier transforming in xy-plane, these equations can be written as

$$\begin{eqnarray}&&{\left.{\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{L}}}\right|}_{z={z}_{{\rm{L}}}}={\left.{\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{C}}}\right|}_{z={z}_{{\rm{L}}}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial {\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{L}}}}{\partial z}\right|}_{z={z}_{{\rm{L}}}}={\left.\displaystyle \frac{\partial {\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{C}}}}{\partial z}\right|}_{z={z}_{{\rm{L}}}},\end{eqnarray} \tag{ 14 }$$

where ${\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{L}}}$ and ${\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{C}}}$ are the Fourier components of ψ^L and ψ^C, respectively. From Eqs. (8) and (9), ${\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{L}}}$ and ${\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{C}}}$ can be written as

$$\begin{eqnarray}\begin{array}{rcl}\,{\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{L}}} & = & \displaystyle \sum _{{G}_{z}}{f}_{{{\boldsymbol{k}}}^{\mathrm{inc}}}({\boldsymbol{G}})\exp \left[i({k}_{z}^{\mathrm{inc}}+{G}_{z})\cdot z\right],\\ & & +\displaystyle \sum _{{k}_{z}^{{\rm{L}}}}\displaystyle \sum _{{G}_{z}}{C}_{{k}_{z}^{{\rm{L}}}}{f}_{{{\boldsymbol{k}}}^{{\rm{L}}}}({\boldsymbol{G}})\exp \left[i({k}_{z}^{{\rm{L}}}+{G}_{z})\cdot z\right],\end{array}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{C}}}=\displaystyle \sum _{{k}_{z}^{{\rm{C}}}}\displaystyle \sum _{{G}_{z}}{C}_{{k}_{z}^{{\rm{C}}}}{f}_{{{\boldsymbol{k}}}^{{\rm{C}}}}({\boldsymbol{G}})\exp \left[i({k}_{z}^{{\rm{C}}}+{G}_{z})\cdot z\right].\end{eqnarray} \tag{ 16 }$$

In the same way, for ψ^C and ψ^R to be continuous at the right boundary at z = z_R (=L_sl), the following boundary conditions have to be satisfied for any G _xy ,

$$\begin{eqnarray}&&{\left.{\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{C}}}\right|}_{z={z}_{{\rm{R}}}}={\left.{\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{R}}}\right|}_{z={z}_{{\rm{R}}}},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial {\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{C}}}}{\partial z}\right|}_{z={z}_{{\rm{R}}}}={\left.\displaystyle \frac{\partial {\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{R}}}}{\partial z}\right|}_{z={z}_{{\rm{R}}}},\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}&&{\phi }_{{{\boldsymbol{G}}}_{{xy}}}^{{\rm{R}}}=\displaystyle \sum _{{k}_{z}^{{\rm{R}}}}\displaystyle \sum _{{G}_{z}}{C}_{{k}_{z}^{{\rm{R}}}}{f}_{{{\boldsymbol{k}}}^{{\rm{R}}}}({\boldsymbol{G}})\exp \left[i({k}_{z}^{{\rm{R}}}+{G}_{z})\cdot z\right].\end{eqnarray} \tag{ 19 }$$

Solving Eqs. (13), (14), (17), and (18) for every G _xy results in making the wavefunction continuous in the whole region. As a result, we can determine the transmission wavefunction ψ^R corresponding to an incident wavefunction ${\varphi }_{{{\boldsymbol{k}}}^{\mathrm{inc}},E}^{\mathrm{inc}}$ in the left bulk region.

2.4. Probability current density

Probability current density J_z is expressed as ^31,32)

$$\begin{eqnarray}&&{J}_{z}=\displaystyle \frac{{\hslash }}{2{mi}}\displaystyle \frac{1}{{S}_{{xy}}}{\int }_{{S}_{{xy}}}\left({\varphi }^{* }\displaystyle \frac{\partial \varphi }{\partial z}-\varphi \displaystyle \frac{\partial {\varphi }^{* }}{\partial z}\right)d{{\boldsymbol{r}}}_{{xy}},\end{eqnarray} \tag{ 20 }$$

where S_xy is the area in the plane perpendicular to the transport direction, m = bm₀, and r _xy = (x, y). The wavefunction of an incident wave in bulk 4H-SiC ${\varphi }_{{{\boldsymbol{k}}}^{\mathrm{inc}}}^{\mathrm{inc}}$ is expressed as

$$\begin{eqnarray}&&{\varphi }_{{{\boldsymbol{k}}}^{\mathrm{inc}},E}^{\mathrm{inc}}=\displaystyle \sum _{{\boldsymbol{G}}}{f}_{{{\boldsymbol{k}}}^{\mathrm{inc}}}({\boldsymbol{G}})\exp \left[i({{\boldsymbol{k}}}^{\mathrm{inc}}+{\boldsymbol{G}})\cdot {\boldsymbol{r}}\right].\end{eqnarray} \tag{ 21 }$$

Substituting Eq. (21) into Eq. (20), the probability current density associated with ${\varphi }_{{{\boldsymbol{k}}}^{\mathrm{inc}},E}^{\mathrm{inc}}$ is given by

$$\begin{eqnarray}\begin{array}{rcl}{J}_{z}^{{\rm{inc}}} & = & \displaystyle \frac{\hslash }{2m}\displaystyle \frac{1}{{S}_{{xy}}}{\int }_{{S}_{{xy}}}\exp \left[i({{\boldsymbol{G}}}_{{xy}}-{\boldsymbol{G}}{{\rm{{\prime} }}}_{{xy}})\cdot {{\boldsymbol{r}}}_{{xy}}\right]d{{\boldsymbol{r}}}_{{xy}}\\ & & \times \displaystyle \sum _{{{\boldsymbol{G}}}^{{\rm{{\prime} }}}}\displaystyle \sum _{{\boldsymbol{G}}}{f}_{{{\boldsymbol{k}}}^{{\rm{inc}}}}^{\ast }({{\boldsymbol{G}}}^{{\rm{{\prime} }}}){f}_{{{\boldsymbol{k}}}^{{\rm{inc}}}}({\boldsymbol{G}})(2{k}_{z}^{{\rm{inc}}}+{G}_{z}+G{{\rm{{\prime} }}}_{z})\\ & & \times \,\exp \left[i({G}_{z}-G{{\rm{{\prime} }}}_{z})\cdot z\right].\\ & = & \displaystyle \frac{\hslash }{2m}\sum _{{{\boldsymbol{G}}}^{{\rm{{\prime} }}}}\sum _{{\boldsymbol{G}}}{\delta }_{{{\boldsymbol{G}}}_{{xy}},{\boldsymbol{G}}{{\rm{{\prime} }}}_{{xy}}}{f}_{{{\boldsymbol{k}}}^{{\rm{inc}}}}^{\ast }({{\boldsymbol{G}}}^{{\rm{{\prime} }}}){f}_{{{\boldsymbol{k}}}^{{\rm{inc}}}}({\boldsymbol{G}})\\ & & \times \exp \left[i({G}_{z}-G{{\rm{{\prime} }}}_{z})\cdot z\right](2{k}_{z}^{{\rm{inc}}}+{G}_{z}+G{{\rm{{\prime} }}}_{z}).\end{array}\end{eqnarray} \tag{ 22 }$$

Using Eq. (10), the probability current density associated with ψ^R is given by

$$\begin{eqnarray}\begin{array}{rcl}{J}_{z}^{{\rm{tra}}} & = & \displaystyle \frac{\hslash }{2m}\sum _{{{\boldsymbol{G}}}^{{\rm{{\prime} }}}}\sum _{{\boldsymbol{G}}}{\delta }_{{{\boldsymbol{G}}}_{{xy}},{\boldsymbol{G}}{{\rm{{\prime} }}}_{{xy}}}\sum _{{k}_{z}^{{\rm{R}}{\text{'}}^{}}}\sum _{{k}_{z}^{{\rm{R}}}}{C}_{{k}_{z}^{{\rm{R}}{\text{'}}^{}}}^{\ast }{C}_{{k}_{z}^{{\rm{R}}}}{f}_{{{\boldsymbol{k}}}_{{xy}}^{{\rm{R}}},{k}_{z}^{{\rm{R}}{\text{'}}^{}}}^{\ast }({{\boldsymbol{G}}}^{{\rm{{\prime} }}})\\ & & \times \,{f}_{{{\boldsymbol{k}}}_{{xy}}^{{\rm{R}}},{k}_{z}^{{\rm{R}}}}({\boldsymbol{G}})\exp \left[i({k}_{z}^{{\rm{R}}}-{k}_{z}^{{\rm{R}}{\text{'}}^{}\ast }+{G}_{z}-G{{\rm{{\prime} }}}_{z})\cdot z\right]\\ & & \times \,({k}_{z}^{{\rm{R}}}+{k}_{z}^{{\rm{R}}{\text{'}}^{}\ast }+{G}_{z}+G{{\rm{{\prime} }}}_{z}).\end{array}\end{eqnarray} \tag{ 23 }$$

The tunneling probability for the kinetic energy along the transport direction is then obtained from the ratio of the probability current densities of the incident and transmitted waves as

$$\begin{eqnarray}&&{T}_{{\rm{C}}-\mathrm{EPM}}({E}_{z})=\displaystyle \frac{{J}_{z}^{\mathrm{tra}}}{{J}_{z}^{\mathrm{inc}}}.\end{eqnarray} \tag{ 24 }$$

Here, E_z = E − E_xy with E_xy being the kinetic energy correspondng to k _xy in the considered band structure. We call this method connecting wavefunction approach in the C-EPM. This method does not use the WKB approximation.

2.5. WKB approximation

For comparison, we calculate the tunneling probabilities within the WKB approximation. The WKB tunneling probability using the electronic states in bulk 4H-SiC by the B-EPM is given by

$$\begin{eqnarray}&&{T}_{{\rm{B}}-\mathrm{EPM}}({E}_{z})=\exp \left[-2{\int }_{{z}_{1}}^{{z}_{2}}\kappa ({E}_{z}-{V}_{{\rm{b}}}(z)){dz}\right],\end{eqnarray} \tag{ 25 }$$

where z₁ and z₂ are the classical turning points. The damping factor κ is obtained from the complex band structure of bulk 4H-SiC. When using the PBA, the tunneling probability is expressed as

$$\begin{eqnarray}&&{T}_{{\rm{B}}-\mathrm{PBA}}({E}_{z})=\exp \left[-2{\int }_{{z}_{1}}^{{z}_{2}}\sqrt{\displaystyle \frac{-2{m}_{z}({E}_{z}-{V}_{{\rm{b}}}(z))}{{\hslash }^{2}}}{dz}\right],\end{eqnarray} \tag{ 26 }$$

where m_z is the effective mass of bulk 4H-SiC along the transport direction. We call this method bulk PBA (B-PBA).

In this section, the calculated complex band structures and tunneling probabilities are presented. As in our previous study, ²¹⁾ we consider a SBD whose donor density N_d and barrier height without the image force effects are 2.6 × 10¹⁸ cm⁻³ and 1.46 eV, respectively. The parameters for the pseudopotentials and the material parameters of 4H-SiC are also the same as those in our previous study. ²¹⁾ In the present study, k _xy is assumed to be that at the M point in bulk 4H-SiC. Note that we set the cut-off of G to $| {\boldsymbol{G}}{| }^{2}\leqslant (33+1/3){\left(2\pi /{a}_{0}\right)}^{2}$ for φ_{n, k} ( r ) in the C-EPM, $| {\boldsymbol{G}}{| }^{2}\leqslant (25+1/3){\left(2\pi /{a}_{0}\right)}^{2}$ for φ_{n, k} ( r ) in the B-EPM, and $| {\boldsymbol{G}}{| }^{2}\leqslant 11{\left(2\pi /{a}_{0}\right)}^{2}$ for V_ps( r ) , where a₀ is the a-axis lattice constant of 4H-SiC (a₀ = 0.308 051 nm ³³⁾). The length of one period of the superlattice L_sl is set to be 6L_cell, where the c-axis lattice constant L_cell is 1.008 480 nm. ³³⁾ The minimum barrier potential, ${V}_{{\rm{b}}}^{\mathrm{low}}$, [see Fig. 1(c)] is set to ${V}_{{\rm{b}}}^{\mathrm{low}}=-1.65\,\mathrm{eV}$.

Figure 2 shows the complex band structures of bulk 4H-SiC by the EPM (blue dotted line) and those by the PBA (red line). As shown in Fig. 2, there is a forbidden region within the conduction band of bulk 4H-SiC along the transport direction. To prevent the electron injection to this forbidden region, we introduce the constant potential region of V_b(z) as shown in Fig. 1.

Zoom In Zoom Out Reset image size

Figure 3 shows the energy dependence of the tunneling probabilities by the B-EPM, the B-PBA, and the C-EPM. The calculation results obtained by fitting to reproduce the experimental tunneling currents ²⁰⁾ with the B-PBA are also shown with a black dotted line as pseudo-experimental results. Note that T_C-EPM(z) corresponds to the tunneling probabilities for the lower conduction band passing through k = E = 0 in Fig. 2. We see that T_C-EPM(z) as well as T_B-EPM(z) overestimates the experimental results ²⁰⁾ because of considering the non-parabolicity of the complex bands as we mentioned in our previous study. ²¹⁾ It should be noted that considering the effect of k _xy on T_B-EPM(z) leads to improving the reproduction of the experimental result ²⁰⁾ by the EPM. The C-EPM results generally agree with the B-EPM results, but there are some large deviations at a few energy points such as E ≃ −1.46 eV and −1.37 eV. To understand the origin of this phenomena, we plot the complex band structure of the superlattice in Fig. 4. We find some flat real bands appearing, which implies that the resonant tunneling occurs at these energies. The peaks in T_C-EPM(z) coincide with the energies where the resonant tunneling occurs in the superlattice structure. These peaks may be a numerical artifact and could be removed. Apart from these peaks, we see that there is not significant difference between T_C-EPM(E_z ) and T_B-EPM(E_z ). This suggests that we can safely use the WKB approximation with electronic states without considering the barrier potential when calculating the tunneling current in the heavily doped 4H-SiC SBDs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here, according to Ref. 24, the condition for the approximate tunneling probability by the WKB approximation and the exact probability to have the similar functional forms is that ${\rm{\Delta }}z\cdot {\bar{{\rm{\Phi }}}}^{\tfrac{1}{2}}$ is sufficiently large, where Δz and $\bar{{\rm{\Phi }}}$ are the length of the barrier potential and the mean barrier height, respectively. This corresponds to the tunneling probability T ≪ 1. In the present model, Δz is at least about 4 nm and Φ is larger than 1 eV, leading to T of at most 10⁻⁶. This means that T ≪ 1 holds, i.e. the condition is considered to be satisfied. Hence, the present result that it is reasonable to use the WKB approximation for the considered barrier potential is consistent with the discussion in Ref. 24.

Finally, we comment on the possible causes of the discrepancy between our calculation and the experimental results. First, in our model, the electronic states in bulk 4H-SiC are used for those in metals, which could lead to the overestimation of the tunneling probability. Moreover, the barrier potential V_b(z) is simply kept constant at both ends of the Schottky barrier, which may lead to deviations from the actual electronic states. The scatterings, such as electron-phonon scattering, may also influence both the electronic states and transport phenomena. To obtain a definite answer, we have to take all these effects into account.

The tunneling probability of heavily doped 4H-SiC SBDs is calculated without using the WKB approximation. First, we calculate the complex bands of bulk 4H-SiC and those of a periodic superlattice model with the barrier potential by the EPM. Then, by making the wavefunctions continuous in the whole region, we determine the transmission wave for a given incident wave and obtain the tunneling probability using the probability current densities for the transmission and incident waves. Comparing the calculated tunneling probability with that by the WKB approximation, we find that there is not much difference between them. These results suggest that we can safely use the WKB approximation with the bulk 4H-SiC electronic states when calculating the tunneling current in the heavily doped 4H-SiC SBDs.

This work was supported by JSPS KAKENHI Grant No. JP21H05003.