Structure and electron affinity of the 4H–SiC (0001) surfaces: a methodological approach for polar systems

The ability to accurately and consistently determine the surface electronic properties of polar materials is of great importance for device applications. Polar surface modelling is fundamentally limited by the spontaneous polarisation of these materials in a periodic boundary condition scheme. Surface data are sensitive to supercell parameters, including slab and vacuum thicknesses, as well as the non-equivalence of surface adsorbates on opposite surfaces. Using 4H–SiC as a specific case, this study explores calculation of electron affinities (EAs) of (000 1̄ ) and (0001) surfaces varying chemical termination as a function of computational parameters. We report the impact in terms of band-gap, electric fields across the vacuum and slab for single and double cell slab models, where the latter is constructed with inversional symmetry to eliminate the electric field in the vacuum regions. We find that single cells are sensitive to both slab and vacuum thickness. The band-gap narrows with slab thickness, ultimately vanishing and inducing charge transfer between opposite surfaces. This has a consequence for predicted EAs. Adsorbate species are found to play a crucial role in the rate of narrowing. Back to back cells with inversional symmetry have larger electric fields present across the slab than the single slab cases, resulting in a greater band-gap narrowing effect, but the vacuum thickness dependence is completely removed. We discuss the relative merits of the two approaches.


Introduction
Silicon carbide, and in particular its 4H-polytype, exhibits an array of attractive physical and electrical properties, including high thermal-conductivity, large breakdown-field and high carrier-mobility [1,2]. The wide band-gap and radiation resistance makes applications in hostile environment, such as high T, feasible [3][4][5]. * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Surfaces and interfaces are of great significance in technological applications of SiC. For field-emission and electrochemical applications, work functions, ionization potentials (IPs) and electron affinities (EAs) are of great importance, and an approximation for band-offsets in heterojunctions is the difference in either IPs or EAs [6]. The EA is the energy difference between an electron in vacuum and in the conduction band of a non-metal, so it is affected by the surface, including its chemical termination, with a layer of surface dipoles producing a step up or down in electron-energy, potentially by several volts. We have recently reported the EAs of nonpolar surfaces of 4H-SiC [7], but there is a paucity of data relating to EAs of polar SiC surfaces. This is partly because The C and Si labels represent the positions of the carbon and silicon atoms in each bilayer respectively. The large electronegativity of carbon compared to silicon (2.55 and 1.90 respectively [15]) results in the charge density peaking closer to the carbon atom. An array of dipoles in the [0001] surface normal direction within a PBC scheme explains the origin of spontaneous bulk polarization. hexagonal polytypes of SiC display spontaneous electric polarisation (figure 1), which affects surface geometry, electronic structure and interaction with adsorbates. Computationally the intrinsic polarity leads to a significant degree of challenge in atomistic modelling. In particular, it is common practice when employing quantum-mechanically based methods such as density functional theory (DFT), to represent surfaces using periodic boundary conditions (PBCs) in three directions. The bulk electric polarization requires a cancelling electric field across the vacuum separating slabs of continuous crystalline material. The thicker the slab, the bigger the potential difference across it, and for a fixed vacuum layer, the bigger the cancelling electric field.
Some previous studies aiming to address this issue have introduced additional difficulties. For example, references [8] and [9] propose the use of group-IV layers in the middle of the slab to artificially enforce inversional symmetry not present in polar materials. Lattice matching the polar material and the group-IV layer is not generally possible, and mismatched lattices introduce strain affecting electronic structure and geometry.
An existing intrinsic-dipole correction method based on the theory proposed by Neugebauer and Scheffler [10] and improved upon by Bengtsson [11] has been implemented in some DFT packages [12]. The electric-field artifact due to PBCs is compensated for by the inclusion of a cancelling electric field: a field-free vacuum region is established for each surface, with a discontinuity in the potential in the middle of the vacuum. A similar approach [13,14] utilizes pseudo-hydrogen with fractional charges and/or altered surface-adatom distances in conjunction with a self-consistent Laplace correction scheme [13] to eliminate the electric-field in the vacuum. Such additional chemical components introduce challenges in systematically modelling different surface terminations.
Although both approaches have merit, they are not universally available for all circumstances and it remains important to explore the impact of modelling parameters such as slab thickness when determining observable surface properties that are directly linked to the electrostatic potential in vacuum, such as is the case with the electron affinity. Due to the bulk polarization, increasing the slab thickness increases the potential difference between the surfaces, reducing and ultimately eliminating the band-gap. This effect is not solely a property of the PBC. It is not clear that the choice of slab thickness is always considered or justified in this respect. To our knowledge, no study so far has addressed a large data set composed of structures many lattice constants thick, with particular emphasis on adsorbates and how they may influence choice of underlying slab size.
Ideally a consistent and broadly applicable method independent of cell parameters and choice of bulk polarization cancellation technique would permit accurate and comparable surface species analysis with experiment. In the absence of a built-in dipole correction scheme, this work attempts to assess the viability of and directly compare two cell construction types, shedding light on key structural choices. Finally, we demonstrate the utility of a modified electron affinity calculation with the potential for broad applicability among polar surfaces. We report the results of density-functional modelling of polar 4H-SiC surfaces, presenting computed values of EAs for H-, F-and Cl-terminated surfaces.

Methodology
Calculations are performed using the AIMPRO [16,17] modelling software package, using the local density approximation [18]. Atoms were modelled using pseudo potentials [19], with 1, 4, 4, 7 and 7 electrons in the valence sets for H, C, Si, F and Cl, respectively. Kohn-Sham functions were expanded using atom-centred, independent s-and p-type Gaussian orbital functions [20] of four widths for H, with the addition of two sets of d-type orbitals for Si, C, Cl and F to account for polarization. This amounts to 16 functions per H atom and 28 per atom in all other cases. A plane-wave expansion of the charge-density [16] and Kohn-Sham potential was used to create the Hamiltonian matrix-elements, with a 175 Ha cut-off yielding well converged total energies (better than 0.1 meV atom −1 ) with respect to this parameter. Integration over the Brillouin-zone for calculation of the total energy follows the Monkhorst-Pack (MP) sampling scheme [21]. (0001) and (0001) surfaces were made up from non-primitive, orthorhombic cells with in-plane lattice vectors along [1120] and [1100]. All surface cells have two-dimensional Brillouin zones due to confinement along the c-axes, and were sampled using 10 × 6 in-plane MP-sampling.
Cells ranging in thickness from 3c to 16c with 70 a.u. of vacuum were constructed. For the 6c-thick slab vacuum thickness was also varied in the range 70-260 a.u. In-plane lattice vectors were fixed to those of bulk, and all atoms allowed to move during structural optimization. The energy was deemed  The projected density of states (pDoS) were obtained using Gaussian broadening of 0.1 eV and an MP sampling-grid of 100 × 60 k-points.

Electron affinity calculations
The standard method for determining a surface EA based upon aligned electrostatic potentials [25] suitable for nonpolar cases is problematic for polar cases where there are electric fields within each slab and vacuum. Figure 3 illustrates the underlying gradient in the bulk and slab cases. As discussed above, there are several dipolar correction schemes employed in literature. Despite this, to our knowledge there is no widely agreed method for calculating the electron affinity of a given polar surface post-correction. As the nature of the system is similar to that of non-polar material, it is evident that some modified form of the accepted alignment procedure can be employed. In an uncorrected scheme, there is a cancelling electric field in the vacuum, the periodic potential varies with an underlying triangular waveform. For triangular potentials where the vacuum potential is not well defined, one could approximate an EA as the difference in energy of the conduction band just inside the SiC and the vacuum just outside. In practice one might post hoc add a constant electric field cancelling the vacuum field exactly to obtain a constant potential in the vacuum region. Addition of a cancelling field everywhere results is a series of vacuum potentials separated by twice the potential difference across the original SiC slab.
Alternatively, rather than cancelling the field present in the vacuum, a field can be added to cancel the field inside the SiC, as illustrated in figure 4. The valence-band maximum and conduction-band minimum inside the slabs are position dependent, but their values close to the surfaces can be obtained by comparison with the bulk electrostatic potential in a fashion similar to that for non-polar cases [25]. The vertical lines in figure 4 indicate the points where the charge density has decayed to negligible value. The energy difference between the conduction band minimum and where these vertical lines cross the electrostatic potential is the estimate of the local EA. We term this approach the triangular potential method.
The electric-field in the vacuum can also be eliminated by construction of periodic systems with pairs of slabs alternating in polar orientation between [0001] and [0001] along z, similar to the symmetric slabs constructed to investigate twodimensional electron gases localized at perovskite interfaces [26,27]. We refer to these as back-to-back (b2b) systems. This approach has the disadvantage of doubling the system size, significantly increasing computational load, but eliminates the electric field in the vacuum by construction. Although there is no net dipole, since each slab has an internal electric field the pairs form an electric quadropole so there are alternating positive and negative electrostatic potential steps between constant vacuum potentials along [0001] (figure 5). The EAs can  be obtained from b2b calculations by fitting the bulk electrostatic potentials in a similar fashion to the triangular potential cases. Additionally, this method has quite general applicability for the calculation the EAs of a polar system.

Results
To explore the use of triangular potential and b2b geometries in determination of EAs, including the roles of the slab and vacuum thicknesses, we have analyzed electrostatic potentials for a wide range of conditions. We present the results and analysis first for a triangular potential system for each termination type, and then for b2b.

Si and C termination
Si and C terminated surfaces are constrained in our calculations to 1 × 1 structures, resulting in one dangling bond per surface atom. Figure 6 illustrates the optimized surface geometries, showing a significant relaxation, particularly on the (0001) face. Bondlengths and displacements are consistent with previous theoretical work on unterminated 4H-SiC surfaces [28,29].
The band structure of the unterminated 6c-thick slab is shown in figure 7. The narrowing of the band-gap for the surfaces is a consequence of the built-in electric field, so the conduction-band minimum for the slab corresponds to the opposite surface for the valence-band maximum. In the 0.0-0.7 eV range we find two pairs of surface states. The electron chemical potential crosses the surface states in this energy range, so the 1 × 1 surfaces are metallic.
The pDoS for layers in the vicinity of the two surfaces are plotted in figure 8. The peak in the Si-face pDoS is 0.1 eV higher than that of the carbon face, and slightly broader. A stronger electronegativity on the carbon surface atoms and  more spatially diffuse orbitals, particularly in the directions parallel to the surface normal, explain these differences, with similar states observed for 6H-SiC (0001)-(1 × 1) surfaces [31]. Analysis of the valence-band maximum orbitals indicates a significant overlap between surface and bulk states centred on carbon atoms, consistent with previous studies [32]. While both faces exhibit surface states dominated by p-contributions, we find that the surface states extend relatively deeper into the  C-face, consistent with these bands lying closer to the valenceband maximum in the adjacent region of the slab than for the Si-face.
Modelling of unterminated 4H-SiC forms a baseline for the investigation of chemically satisfied cases, and the impact on the slab and vacuum variables, and it is to this that we next turn.

Vacuum thickness
3.2.1. Triangular potential. The space between periodic images has to be sufficient to ensure no significant overlap of the evanescent surface states and produce a linear potential. However, an excessively wide vacuum impacts computation time. PBCs (section 1) applied to the non-centrosymmetric polar systems leads to electrostatic interactions between images, irrespective of vacuum size. Figure 9 shows the electric field across the slab, E slab , converges to different values depending upon the terminating species. This shows significant image-image electrostatic contributions to E slab for vacuum thicknesses.
Additionally much larger vacuum thicknesses were also examined, and the results are listed in table 1. It can be seen that from 500-1500 a.u. data, the values of the field across the slab are far from converged even at 250 a.u.
The magnitude of the electric field across the vacuum is such that it cancels the potential difference across the slab. Naturally, for a fixed slab the electric field across the vacuum, E vac , must tend to zero as a function of vacuum size. Figure 10 shows this trend.   Figure 11 shows a plot of the potential across the vacuum, ΔV vac , as a function of vacuum thickness. It reveals a modest increase over the range, consistent with the changes in E slab .

Back to back.
We now turn to the results from b2b structures. Figures 12 and 13 show ΔV vac and E slab , respectively, for the b2b simulations. In contrast to the triangular potential the values are essentially independent of vacuum width.
The b2b values are consistent with the data in table 2. H-termination produces the largest electric field of the cases examined. This is consistent with the surface dipoles in each case: for H-termination they point in the [0001] direction due to H having an electronegativity between those of C and Si, whereas for F and Cl the surface dipoles are all directed towards the vacuum, so partially cancel out. The trend of the triangular potentials to converge to those of the vacuumthickness independent values of the b2b model (figures 12 and 13) strongly advocate the use of the b2b approach to mitigate vacuum thickness effects.  We now move to the dependence of derived properties on the thickness of the polar material.

Slab thickness
3.3.1. Triangular potential method. One might seek to converge properties of surface by increasing the thickness of the underlying crystalline material. For non-polar materials such an approach is limited only by the computational effort required. For polar materials we generally still seek structural models that converge, for example, surface geometry, and the underlying crystal structure should closely match the equilibrium geometry of the bulk material. However, polar materials potentially possess an upper limit of slab thickness, since the band-gap narrows and ultimately vanishes with thickness. When the valence band maximum at one surface lies higher in energy than the conduction band minimum at the other, charge transfer will occur. Figure 15 shows this narrowing of  the band-gap, and in for all terminations there is a reduction in the band-gap with thickness. The magnitude of the band-gap for a given slab-thickness is termination-species dependent. For Cl-termination the bandgap vanishes for a relatively thin slab, being around 11c for the triangular potential method. In comparison, for F-termination this occurs at around 16c, and for H-termination a band-gap is present for all thicknesses examined. The physical origin of the termination-dependence lies in the existence of surface states in the case of halogen termination, which are absent for H-termination [33].
The pDoS for the Cl-termination (figure 14) reveal a number of key features. The Cl pDoS on each surface exhibit a sharp peak followed by a hard shoulder and drop-off to 0 eV cell −1 at higher energies. The shape is similar in nature to that observed in 2D materials such as graphene [34] and hexagonal boron-nitride [35]. The prominent surface-state peaks are predominantly lone pair Cl-centred p-orbitals with small contributions from the top-most SiC layer ( figure 16).
It is useful to examine the correlation of the host bands with the surface states. When the (0001) surface-states are fixed at zero ( figure 14(a)) the peaks for the (0001) surface-states shown in shades of blue between around −0.2 eV to 1.0 eV, that track with where the valence band DoS vanishes (shades of grey) which co-incidently cover the same energy range. The (0001) surface-states (blue peaks) increase in energy in approximately equal increments until the slab thickness coincides with disappearance of the band-gap (figure 15) at 11c (ninth peak). For greater slab thicknesses the charge transfer between the two surfaces offsets any further shift in the relative position of the surface states.  Figure 15 shows the band-gap as approximately constant up to 5c. For the relatively thin slabs the (0001) surfacestates are lower in energy than those from the (0001) face. In figure 14(a) the (blue) (0001) surface state high-energy limit only exceeds that of the (0001) surface states at around 1.2 eV for thicknesses exceeding 5c. Once this happens, as the thickness increases the (0001) surfaces states continue to increase in energy relative to those of the (0001) surface, corresponding to the most rapid phase of band-gap reduction in figure 15. The same effect can be seen in figure 14(b) relative to the (0001) surface states. The transition can also be seen in the wave functions (figure 16); the highest energy surface states are p-orbitals centred on Cl atoms of the C face for 3 and 4 lattice constants, transitioning to Cl atoms of the Si face for all Cl-SiC slabs 5 lattice constants and above.
The general trend is repeated for fluorine termination. The highest energy occupied states are associated with F atoms on the (0001) surface for slabs up to 4c in thickness, being overtaken in energy by the (0001) surface states for thicker slabs (the initial flat region in figure 15). This effect is absent in Hterminated surface, as there are no lone-pair orbitals in this case to form the characteristic surface states.
The three terminating species examined for this study illustrate the complex dependence of the band-gap on slab thickness and different surface states associated with the (0001) and (0001) surfaces.
We now focus upon the regime under which the band-gap vanishes. Our calculations show that where the band-gap is eliminated, the extent of charge transfer is that which generates an electric-field in opposition of the bulk polarization. Its magnitude is such that the associated electrostatic potential cancels out the increase that corresponding to the bulk polarization. This effect is shown in the change of gradient at around 10c for Cl-termination in figures 17-19. Focusing on figure 17, a consistent trend is observed from 3 to 9 lattice constants for all terminations. As the slabs increase in thickness, the proportionate contribution of the surface dipoles decreases, and all terminations show the same general trend. For Cl termination the loss of a band-gap increases the rate of decrease in electric field, and fixes the potential across the slab (and hence the vacuum) to a constant ( figure 19).
The absolute values of the electric fields and the potential drops across the slabs, ΔV, can be further illustrated in terms of the dipole size and direction on each surface, touched up previously. Table 2 shows the direction of the electric-dipoles for each termination. The dipoles for hydrogen termination point in the same direction, and in the same net direction as the bulk polarization. The sum of the two dipoles increase the electric field. In contrast, the fluorine and chlorine cases have  surface dipoles that will tend to cancel, leading to smaller electric fields relative to H-termination, consistent with figures 18 and 19. The magnitudes of dipole of F and Cl must be similar as they lead to very similar profiles of electric field strength and potential difference.

3.3.2.
Back to back systems. Some elements are common between b2b and triangular potential approaches. Although there is no electric field in the vacuum, there is inside the slabs, so there remains a slab-thickness dependence of the band-gap. We find the same general trends in band-gap as for the triangular potential approach ( figure 20 and 15), but b2b systems show a more rapid reduction in band-gap with slab thickness. This is a consequence of the higher electric field magnitudes in the b2b slabs, as shown in figure 21. At all thicknesses and surface-species the electric field in the slab is higher in the b2b model than the single-slab approach with a triangular potential. Taking 6c-slabs as an example, E slab for the H, F and Cl terminations with the b2b model are 0.014, 0.013 and 0.013 eV a.u. −1 respectively, compared to 0.009, 0.011 and 0.010 eV a.u. −1 for the triangular potentials. The consequences  of a larger slab electric field are that the potential difference between states is higher, and their dispersion through space is larger.
The confinement effect can also be seen in E slab values with a direct comparison between methods. In section 3.3.1, E slab values in the region of 4-9 lattice constants for figure 17 decrease linearly for all terminations. H-and F-termination E slab values in figure 21 are both higher and more consistent with thickness compared to the same variables for the triangular method.
Taking H-surface as a specific example, figure 17 shows that with varying thickness H-termination has the lowest E slab values for all thicknesses except for when Cl becomes metallic at 10 lattice constants. However this trend is the opposite for the back to back case where in figure 21 H-surfaces have higher E slab values for all thicknesses.

EA dependence
To illustrate the impact of the use of triangular potential or b2b methods, we have calculated the EA for the systems studies for this paper. The results are listed in table 3.  Table 3 lists estimated EA values as a function of termination, slab thickness and vacuum thickness. The ranges in values indicated provide the limits accumulated over the various slab and vacuum thicknesses. For triangular potential calculation the impact of slab-thickness is modest (0.02-0.06 eV). Vacuum thickness has a greater impact (0.04-0.08 eV). For b2b calculations, slab thickness yields a more significant variation than the triangular potential approach (0.09-0.19 eV). However, the vacuum-width dependence is completely removed. The electric field in b2b slabs where band-gaps remain is much larger than in the triangular potential scheme, reducing the range of viable slab thicknesses for surface-property modelling.

Discussion and conclusions
Reliable methods to determine surface properties for polar materials, such as structure and electron affinity, are problematic when using PBCs. In this paper we have made a comparison of two schemes to extract electrostatic data for 4H-SiC surfaces with different terminations. Under the assumption that complete saturation of (0001) and (0001) surfaces is possible, we have shown that relatively narrow ranges of values of electron affinity can be obtained from both a single and double slab model.
The existence of an electric field across the polar material in either model necessitates a maximum slab thickness that can be achieved before the band-gap is eliminated with the conduction band minimum at one face to lie below the valence band top at the other. The loss of band gap makes it energetically favourable for charge to move from one surface to the other, generating an addition electric field in opposition to the bulk polarization. In the event that the surfaces are metallic, such as in the unterminated (1 × 1) (0001) surfaces, the dependence of the charge accumulation at the surface will be a combination of rearrangement of charge at the metallic surface that would be expected to occur in nature, and a transfer between opposing surfaces that is simulation-model dependent. Although completely unterminated surfaces are of largely academic interest, less than 100% chemical termination is quite probable. For example, hydrogen termination is a species that has historically proved challenging to achieve full surface coverage [28,36,37]. Incomplete coverage both alters the surface dipole density affecting the extent to which electron affinity is changed compared to an unterminated surface, but also generates surface states that may be non-physically populated by charge transfer between surfaces in the event that the slab thickness is sufficiently large to eliminate the band-gap.
However, there are some general points that can be made regarding a preference between the triangular potential and b2b approaches. The latter completely eliminates the dependence of the vacuum width on the slab properties, such as electron affinity. However, the b2b method shows a more rapid decrease in band-gap with slab-thickness, so that the range of slabs for which charge transfer between the two surfaces can be safely avoided is smaller.
Furthermore, the computational cost of b2b calculations is significantly higher than the triangular potential approach. Conversely, single layer calculations that exhibit the triangular potential show very slow convergence with slab thickness of properties such as the electric field inside the slabs.
So far as the 4H-SiC specific results are concerned, we have estimated EAs for H, F and Cl termination of both (0001) and (0001) surfaces. The largest values are found for Ftermination, being around 5 eV for both surfaces. The smallest EA is obtained for H-termination on the (0001) face at around 1 eV. These values and their associated ranges would be anticipated and indeed are similar to previous 4H-SiC non-polar terminated surfaces [7] with matching surface coverage and species.