Solitons in binary compounds with stacked two-dimensional honeycomb lattices

We model the electronic properties of thin films of binary compounds with stacked layers where each layer is a two-dimensional honeycomb lattice with two atoms per unit cell. The two atoms per cell are assigned different onsite energies in order to consider six different stacking orders: ABC, ABA, AA, ABC$^{\prime}$, ABA$^{\prime}$, and AA$^{\prime}$. Using a minimal tight-binding model with nearest-neighbor hopping, we consider whether a fault in the texture of onsite energies in the vertical, stacking direction supports localized states, and we find localized states within the bulk band gap for ABC, ABA, and AA$^{\prime}$ stacking. Depending on the stacking type, parameter values, and whether the soliton is atomically sharp or a smooth texture, there are a range of different band structures including soliton bands that are either isolated or that hybridize with other states, such as surface states, and soliton bands that are either dispersive or flat, the latter yielding narrow features in the density of states. We discuss the relevance of our results to specific materials including graphene, hexagonal boron nitride and other binary compounds.

Recently, stacking faults (in the vertical, out-of-plane direction) have been considered for thin films of rhombohedral graphite [8,[20][21][22][23][24].As with solitons in the SSH model [13,14,25], stacking faults in thin graphitic films support localized states, manifested as bands in quasi-2D.Such bands tend to appear in pairs at about the same energy because the sequence of intra-and interlayer bonding either side of the fault is fixed, i.e., the texture of bonding strengths away from the fault cannot be changed by the fault.Hence, such faults are effectively coupled soliton-antisoliton pairs [23].
In this paper, we investigate the possible existence of a single, isolated band localized on a stacking fault in quasi-2D materials consisting of layers of atoms on a honeycomb lattice.We consider materials with two different chemical elements such as hexagonal boron nitride [26,27], using a minimal tight-binding model to provide a generic description, modeling stacking faults for different stacking orders.Specifically, for each layer, we consider two nonequivalent atoms per unit cell, labeled A and B, each with a single orbital which is isotropic within the plane, e.g. a p z orbital for sp 2 hybridization in graphene or hexagonal boron nitride.The two atoms have different onsite energies ϵ A = U and ϵ B = −U such that U = (ϵ A − ϵ B )/2 and ϵ A + ϵ B = 0. We consider naturally occurring stackings, i.e., those whereby every atom is either directly above or below another atom, or above or below the center of a hexagon.There are six of them [28][29][30][31], as shown in the top row of Fig. 1, consisting of rhombohedral stacking (ABC), Bernal stacking (ABA), and AA stacking.For each of these, we also include a primed version (ABC ′ , ABA ′ , and AA ′ ) whereby every other layer has the sign of the onsite energy U reversed [28,29].
For rhombohedral stacking (ABC) with alternating onsite energies, dimensional reduction [17], by considering the in-plane wave vector as a fixed parameter, relates the Hamiltonian to that of the Rice-Mele model [32], which is a generalization of the SSH model with alternating onsite energies as well as alternating hopping in one dimension.This is generally not topological, but, at certain wave vectors (k x = k c in Fig. 1), the intra-and interlayer hoppings are effectively equal, and ABC stacking is related to the charge density wave (CDW) model [15,16,[33][34][35][36] which is a one-dimensional Z 2 topological insulator [37] with constant hoppings and alternating onsite energies.Stacking faults for ABC stacking are analogous to solitons in the CDW model [33,34,36], albeit only at a certain wave vector, and our aim is to explore this analogy.
We focus on electronic properties near the corner of the first Brillouin zone (K point) and Fig. 1 summarizes our results.The top row shows the six faultless systems with their electronic band structures.Three of them (ABC, ABA, and AA ′ stacking) exhibit band gaps, and these are the configurations which support localized states on a stacking fault that lie within the band gap.A sharp stacking fault consists of an inversion of the signs of the onsite energy U between two layers without any change in the intra-or interlayer hopping, as shown in the bottom row of Fig. 1; we refer to this as a soliton because it creates a change in texture of the onsite energies either side of the fault.We find a range of different behaviors, depending on the stacking; arrows in Fig. 1 (bottom row) show the energies of states localized on the soliton exactly at the K point.For ABC stacking, there are two local-   The second row (g)-(l) shows the lattices for N = 8 layers with an atomically-sharp fault in the onsite energies occurring between the middle two layers and the bottom of each panel shows the corresponding band structures.Horizontal arrows show the energies at kx = 0 of the states localized on the soliton.In all plots, parameter values are γ0 = 3.16 eV, γ1 = 0.381 eV [40], a = 2.46 Å [41], and the magnitude of the alternating onsite energies is u = U/γ1 = 0.4.ized states.One of them is within the valence band, but another generally lies within the bulk band gap.Depending on parameter values, it may hybridize with surface states, as shown in Fig. 1(g).For ABC ′ , there are also two localized states, but they don't lie near zero energy; the two states within the band gap in Fig. 1(h) are surface states.For ABA stacking, there are two localized states.One of them is within the conduction band, but another generally lies within the bulk band gap, Fig. 1(i), and there are no surface states to hybridize with.As this band eventually moves into the conduction band, we refer to the corresponding texture as an antisoliton whereas the band shown in Fig. 1(g) for ABC stacking moves into the valence band, and its texture is called a soliton.In general, both of these stackings (ABC and ABA) can support either a soliton or an antisoliton, depending on the position of the fault.
For both ABA ′ and AA stacking, there are no states localized on the fault, although the presence of the fault does have an impact on the band structure, Fig. 1(j) and Fig. 1(k), respectively.For AA ′ , a single fault supports four localized states: Two of them are in the bulk bands, but two of them lie within the bulk band gap and, depending on parameter values, they may hybridize with each other as shown in Fig. 1(l).
In the following, we describe the cases shown in Fig. 1, focusing on the systems with a single, isolated soliton band, and modeling their parameter dependence.We describe the methodology used in Section II.We describe rhombohedral stacking (ABC and ABC ′ ) in Section III, Bernal stacking (ABA and ABA ′ ) in Section IV, and AA stacking (AA and AA ′ ) in Section V.The properties of smooth solitons with a finite width are described in Section VI, and defects which consist of a reversal of the signs of the onsite energies on a single layer only [38] are described in Section VII.Section VIII briefly describes the robustness of soliton features in the density of states to the presence of interlayer disorder.In Section IX, we describe ballistic, coherent electronic transport for ABC stacking, comparing the energy dependence of conductivity for systems with and without a soliton.Section X describes the influence of additional tight-binding parameters beyond the minimal model, which only has nearestneighbor intra-and interlayer coupling.Relevance to particular materials and the stability of structural defects is discussed in Section XI.

II. METHODOLOGY
We use a minimal tight-binding model with nearestneighbor intralayer hopping γ 0 and nearest-neighbor interlayer hopping γ 1 .Assuming translational invariance within each layer, we Fourier transform to reciprocal space where q = (q x , q y ) is the in-plane wave vector measured with respect to the center of the Brillouin zone (the Γ point).In a basis of a single orbital on each atomic site written in terms of 2 × 2 blocks.Intralayer blocks are where ±U i are onsite energies and intralayer hopping is described [39] by where a is the in-plane lattice constant.For faultless systems with non-primed stackings, e.g.Fig. 1(a), (c), (e), then U i = U for all i, whereas, for faultless systems with primed stackings, e.g.Fig. 1(b), (d), (f), then U i = U for odd i and U i = −U for even i.
The form of the interlayer block V j in the Hamiltonian (1) depends on the stacking type: ABA (Bernal) : AA : without any dependence on whether the stacking is primed or not.
A system of N layers with an atomically-sharp fault in the onsite energies is labelled (m, n) where m is the number of layers below the fault, n is the number of layers above the fault, and N = m + n.Such a fault is modeled by reversing the signs of U i on the layers above the fault, as shown in Fig. 1(g)-(l).Note that the two sites, A 1 and B 1 , on the first layer are always considered to have fixed energies of ϵ A1 = U and ϵ B1 = −U .If one were to reverse this convention, there would be the same outcomes with electron-hole inversion of the band structures so that, for example, solitons would become antisolitons.
For the band structure plots we numerically diagonalize the Hamiltonian (1).We use tight-binding parameters experimentally measured [40] for Bernal-stacked bilayer graphene γ 0 = 3.16 eV and γ 1 = 0.381 eV, and lattice constant a = 2.46 Å [41].Band structures are calculated in the vicinity of valley K + with wave vector K + = (4π/(3a), 0) by shifting the wave vector as q = K + + k, with k = (k x , k y ).Specifically, we determine the energy eigenvalues E n on a square grid of points centered on K + .The band structures have valley degeneracy (under a suitable rotation) and they are approximately isotropic around each valley at the energies we consider, so we plot them as a function of component k x (with k y = 0), normalized by the characteristic wave vector k c = γ 1 /ℏv [42] where v = √ 3aγ 0 /(2ℏ) is the velocity related to intralayer hopping.By plotting k x as a function of k c and energy E as a function of γ 1 , there are only two free parameters in our model, namely the ratios γ 0 /γ 1 and U/γ 1 .The density of states g(E) per unit energy per unit area (L 2 ) is determined numerically by approximation using a Lorentzian with a finite width δ, States localized on solitons may or may not appear in the bulk band gap, but we focus on those that do.Such states are extended states in the in-plane direction, but are localized on the soliton in the out-of-plane (z) direction.We identify them by observing the form of the state in the atomic basis of the Hamiltonian (1); an example is shown in Fig. 4(a).
For analytic calculations, we expand the function f (q), Eq. ( 3), for |k| ≪ k c , as f (K ξ + k) ≈ − √ 3a(ξk x − ik y )/2 so that the intralayer hopping matrix element for layer i is H AiBi ≈ ℏv(ξk x − ik y ) where ξ = ±1 is a valley index for the valley centers at wave vectors K ξ = ξ(4π/(3a), 0).With this approximation, the phase of H AiBi may be gauged away (using a diagonal unitary transformation) so that the intralayer matrix element may be replaced by H AiBi ≈ ℏvk where k = |k| = (k 2 x + k 2 y ) 1/2 , and the intralayer block is In the interpretation of numerical results, we use dimensional reduction [17] whereby we consider the inplane wave vector k to be a fixed parameter, typically either k = 0 or k = k c .Then, we can relate the model of stacked two-dimensional layers to a one-dimensional tight-binding model.In particular, we use analogy with the Rice-Mele model [32] which is a one-dimensional tight-binding model with two orbitals per unit cell, alternating nearest-neighbor hopping, and alternating onsite energies.It has two phases with chiral symmetry: The Su-Schrieffer-Heeger (SSH) model [13,43,44] with alternating nearest-neighbor hopping but constant onsite energies, and the charge density wave (CDW) model [15,16,[33][34][35][36] with alternating onsite energies but constant nearest-neighbor hopping.In position space with open boundary conditions, the CDW model for a system with J orbitals may be written in the atomic basis as where ±U is the alternating onsite energy and t is the constant nearest-neighbor hopping parameter.This is a two band insulator with a band gap of 2U .For a soliton in the center of the system, i.e. with onsite energies of . . ., +U, −U, +U, −U, −U, +U, −U, +U, . . . it was recently shown [36] that a state is localised on the soliton with an energy E sol that isn't exactly at zero energy, but within the band gap for 0 < t/U < J/2, i.e., where J is the total number of orbitals.For U/t ≫ 1, perturbation theory to second order in t shows [36] that the energy of this level E sol is This energy E sol is not at zero because the chiral symmetry of the CDW model is nonsymmorphic [16,34,36,37,45,46] meaning it involves a translation by half a unit cell, and the ends of the system and a soliton of finite width break the chiral symmetry.The nonsymmorphic chiral symmetry may be written as T a/2 S z [16,36] where Here T a/2 describes translation by half a unit cell, and S z describes sublattice chiral symmetry for the SSH model, The expectation value of T a/2 S z is a generalization of electric polarization, where subscript 'y' is used because operator T a/2 S z is represented as σ y in reciprocal space.In one-dimensional topological insulators with chiral symmetry, the polarization related to a chiral operator takes values ±1 for zeroenergy topological states [14].However, since the nonsymmorphic chiral symmetry is broken in a finite system, the corresponding polarization takes values with magnitude less than one.For U/t ≫ 1, perturbation theory [36] gives for a soliton state, so that p sol y = 1/2 for t/U = 0 and p sol y > 1/2 for t/U > 0. For an antisoliton, p asol y = −p sol y .

III. RHOMBOHEDRAL STACKING
A. ABC

Faultless system
For bulk, faultless ABC stacking, the unit cell consists of two atoms A i , B i on the ith layer, and the Hamiltonian may be written as using the approximation for the intralayer hopping (8) with k < k c , where q z is the z component of the wave vector and d is the layer separation.The two corresponding bulk bands are which have a band gap of 2U .If we consider the inplane wave vector k to be a fixed parameter, then the Hamiltonian ( 16) is that of the Rice-Mele model [32].
For U = 0, the system is rhombohedral graphite which has a bulk energy gap and two flat bands at zero energy ABC stacking (8,8) FIG. 2. ABC stacking with N = 16 layers and a sharp soliton at the center (8,8).The top row shows band structures, the bottom row shows the corresponding density of states g(E) for different values of the magnitude of the alternating onsite energies as (a), (f) U/γ1 = 0.0, (b), (g) U/γ1 = 0.2, (c), (h) U/γ1 = 0.4, (d), (i) U/γ1 = 0.6, and (e), (j) U/γ1 = 2.0.In the band structure plots, horizontal arrows show the energies at kx = 0 of the states localized on the soliton.Note that the axis scales for (e), (j) differ from the rest.In all plots, parameter values are γ0 = 3.16 eV, γ1 = 0.381 eV [40], a = 2.46 Å [41].For the band structures, ky = 0, and, for the density of states, δ = 0.01γ1.
arising from surface states (localised on sites A 1 and B N ), as shown in Fig. 2(a) and (f).The flat bands give rise to a narrow peak in the density of states at zero energy.This band structure may be understood by dimensional reduction [17][18][19] as being related to the SSH model, whereby, on treating the in-plane wave vector k as a fixed parameter, the Hamiltonian of rhombohedral graphite is the same as that of the SSH model.The point in the band structure where k x is equal to the characteristic wave vector k c = γ 1 /(ℏv) is equivalent to the metallic phase of the SSH model, and, thus, k c defines the extent of the flat bands in k space.

Single soliton
For finite U , we consider a sharp fault consisting of two adjacent sites with negative onsite energies as shown in Fig. 1(g).Band structures and density of states for N = 16 layers and a sharp fault occurring between the middle two layers, denoted (8,8), are shown in Fig. 2.There are two bands localized on the fault, as indicated by the arrows in the band structure plots of Fig. 2. Both of these bands move into the valence band for k x ≫ k c , hence we refer to this fault as a soliton [47].We focus on the higher energy band of the two because it is near zero energy.For U ≤ γ 1 /2, this soliton band hybridizes with the flat bands arising from the surface states, Fig. 2(b) and Fig. 2(c).Note that the soliton changes the texture of onsite energies, flipping the sign of the onsite energy on the top surface (site B N ).This can be seen in Fig. 1(g) where the soliton (between sites B 4 and A 5 ) changes the order of the onsite energies for higher layers (with layer index j ≥ 5) so that the final orbital, B 8 , has the same onsite energy as the first, A 1 , in contrast to the faultless case, Fig. 1(a).As a result, both surface states, localized on A 1 and B 8 , are at positive energy and move into the conduction band for k x ≫ k c .
For U > γ 1 /2, Fig. 2(d) and Fig. 2(e), the soliton band and the surface bands are always separate, although the soliton band remains within the band gap.Thus, for U > γ 1 /2, there is a band gap, but it is smaller than the bulk band gap of 2U , Fig. 2(j).The numerical plots may be understood by examining the form of the Hamiltonian (1) at k = 0 because the intralayer hopping −γ 0 f (q) is zero there, and the Hamiltonian (1) becomes block diagonal.The two surface states, on A 1 and B N , are completely disconnected at k = 0 and they give two degenerate states at energy E = +U .There are N − 2 dimers with onsite energies ±U which each give two states of energies E = ± U 2 + γ 2 1 .Finally, the soliton is a dimer of sites B m and A m+1 with onsite energies both at −U , yielding two states with energies E = −U ± γ 1 .The soliton level is the highest energy of the pair, Given that the soliton band begins at energy −U + γ 1 for k = 0, and joins the valence band for k ≫ k c , there will be a part of this band within the bulk band gap for all non-zero values of U .In the regime U ≫ γ 1 and U ≫ ℏvk, it is possible to estimate the energy of the soliton level using perturbation theory in the hopping strength.Using the analytic approximation (8) for the intralayer hopping matrix element gives Note that expressions (17,18) for the energy E sol of the band localized on the soliton are independent of the position of the soliton and the number of layers N in the system.Soliton states are extended states in the in-plane di-rection, but are localized in the out-of-plane (z) direction.At k = 0, the soliton state is localized on the two dimer sites B m and A m+1 , i.e. it is given by ψ sol = (ψ Bm + ψ Am+1 )/ √ 2 where ψ j denotes the atomic orbital on site j.For U ≈ γ 1 /2 and k ≪ k c we derive an effective 3 × 3 Hamiltonian describing the hybridization of the soliton state, ψ sol = (ψ Bm + ψ Am+1 )/ √ 2, with the two surface states ψ A1 and ψ B N .This is done using the linear expansion f (K ξ + k) ≈ − √ 3a(ξk x − ik y )/2, eliminating the other orbitals [39,48] and performing an expansion for k/k c ≪ 1 and |U/γ 1 − 1/2| ≪ 1.In basis ψ A1 , (ψ Bm + ψ Am+1 )/ √ 2, ψ B N , the effective Hamiltonian is where ∆ is a dimensionless contribution to the onsite energies, U/γ 1 = 1/2 + ∆ with |∆| ≪ 1, κ is a dimensionless complex wave vector κ = (ξk For each matrix element in the effective Hamiltonian (19), we retain the leading order term in |κ|.
An antisoliton (m, n) has the same structure as a soli-ton (m, n) but with the signs of all of the onsite energies reversed.At k = 0, the antisoliton is a dimer of sites B m and A m+1 with onsite energies both at U , yielding two states with energies E = U ± γ 1 .The antisoliton state is the lowest state of the pair, and it is localized on the two dimer sites B m and A m+1 as where To analyze the hybridization of the soliton band with the surface states, we simplify the soliton effective Hamiltonian ( 19) by considering a soliton at the center of a system with an even number of layers, m = n = N/2, and we consider a large number of layers where α = (4/5)|κ| 2 and β = (1/ √ 2)(−4/5) m−1 κ m .The eigenvalues of the three states are E/γ 1 = 1/2 + ∆ and 2 or, in terms of physical parameters, there is a flat band E = U and two hybridized bands, (k/k c ) N .The analytic solutions for the flat band E = U and the hybridized bands Eq. ( 24) are plotted in Fig. 3 for energy values near U = γ 1 /2.For k = k c , the Hamiltonian ( 16) is approximately equal to that of the CDW model (this equivalence is only approximate because k/k c is not small when k = k c ) with the replacement t ≡ γ 1 .For k = 0, the soliton state is localized on the two adjacent dimer sites as ψ sol = (ψ Bm + ψ Am+1 )/ √ 2. For k = k c , it remains localized on the soliton, but has a broader extent in position space as shown in Fig. 4(a) for U/γ 1 = 0.6.Fig. 4(b) shows the polarization p y = ⟨ψ|T a/2 S z |ψ⟩, Eq. ( 14), as a function of k x (with k y = 0) for U/γ 1 = 0.6 (solid line) and U/γ 1 = 1.8 (dashed line).For U/γ 1 ≳ 1, the value of p y at k x = k c is in good agreement with the analytic expectation (15) with τ = γ 1 /(2U ) (dashed line), and the value of p y increases as U/γ 1 decreases (solid line).The polarization p y never reaches unity because the nonsymmorphic chiral symmetry is broken by the ends of the system and the finite width of the soliton.

Soliton-antisoliton pair
A soliton-antisoliton pair is denoted (ℓ, m, n) meaning a soliton (with consecutive onsite energies of −U, −U ) after ℓ layers and an antisoliton (with consecutive onsite energies of +U, +U ) after ℓ+m layers, with N = ℓ+m+n.For U ≈ γ 1 /2, the soliton state will strongly hybridize with the bottom surface state (on A 1 ) at energy +γ 1 /2, and the antisoliton state will strongly hybridize with the top surface state (on B N ) at energy −γ 1 /2.We derive an effective 4 × 4 Hamiltonian in the basis ψ A1 , (ψ Here we neglect the contributions (H ), and f n,m is a numerical factor, for n = 1 and m = 1 2/5 for n = 1 and m > 1 28/45 for n > 1 and m = 1 4/5 for n > 1 and m > 1 For U ≈ γ 1 , the soliton and antisoliton state are both near zero energy at k = 0 and they hybridize for k ̸ = 0 (the surface bands are distant, at energies ±U ).We use a two-component basis (ψ where U/γ 1 = 1 + δ with |δ| ≪ 1, |κ| ≪ 1, and g m,ℓ is a Hamiltonian (26) has eigenvalues given by describing hybridization of the soliton and antisoliton states for U ≈ γ 1 and k ≪ k c .

Faultless system
For bulk, faultless ABC ′ stacking, the unit cell consists of four atoms such as A 1 , B 1 , A 2 , B 2 (in Fig. 1(b)), and the Hamiltonian may be written as using the approximation for the intralayer hopping (8) with k < k c , where q z is the z component of the wave vector and d is the layer separation (the length of the unit cell in the z direction is 2d).The four corresponding bulk bands are given by For a thin film of finite thickness, we begin by considering the faultless system for an even number of layers N , Fig. 6 (top row).At k = 0, the intralayer hopping is zero, and the system simplifies as a collection of separate dimers plus two isolated surface states.Thus, the N bands (for N ≥ 4) are at five distinct energies at k = 0 with E = +U (twice) arising from the surface states, and from the dimers: , and E = −U − γ 1 (N/2 times).The surface bands appear at lower energy (E = +U ) than the bulk bands for 0 ≤ U < γ 1 /2, and the bulk bands touch at k = 0 for U = γ 1 .Considering the bands for all k values, the numerical plots in Fig. 6 show that there is always a band within the nominal band gap for U < γ 1 , whereas there is a bandgap for U > γ 1 .

Single soliton
For even N and a single fault, Fig. 6 (second row), we consider k = 0 again.In this case, the N bands (for N ≥ 4) are at eight distinct energies at k = 0 with E = +U and E = −U levels arising from the surface states; the presence of the fault switches the texture of onsite energies, causing the energy of the top surface to flip sign as compared to the faultless case.For dimers at k = 0, the energies are E = +U +γ 1 , E = +U −γ 1 , E = −U +γ 1 , and E = −U − γ 1 , each of these occurring N/2 − 1 times.In addition, there is the dimer consisting of the stacking fault giving energies E = ± U 2 + γ 2 1 .For nonzero k, for 0 < U < γ 1 , there is a small anticrossing between two low-energy bands, then, for U = γ 1 , the bands touch, and, for U > γ 1 , there is a band gap E g = 2(U − γ 1 ).However, the states localized on the soliton are not near low energy, as indicated by the arrows in Fig. 6 (second ABC stacking (5,6,5) row).They are at energies E = ± U 2 + γ 2  1 for k = 0, and then move further away from zero for nonzero k.
The band structure for ABC ′ stacking depends on whether there is an even or odd number of layers, even for a faultless system.The band structure for a faultless system with an odd number of layers N is shown in Fig. 6 (third row).For N ≥ 3, at k = 0, the N bands are at six distinct energies with E = +U and E = −U levels arising from the surface states, and from the dimers: and E = −U −γ 1 , each of these occurring (N −1)/2 times.As with the faultless even N case, the surface bands appear at lower energy (E = ±U ) than the bulk bands for 0 ≤ U < γ 1 /2, and the bulk bands touch at k = 0 and U = γ 1 .Nevertheless, the band structure is more similar to that of the even N case with a fault.This can be seen by considering the bands for nonzero k values: For 0 < U < γ 1 , there is a small anticrossing between two low-energy bands, then, for U = γ 1 , the bands touch, and, for U > γ 1 , there is a band gap E g = 2(U − γ 1 ).
The band structure for an odd number of layers N and a fault near the center is shown in Fig. 6 (bottom row).In this case, the N bands (for N ≥ 5) are at seven distinct energies at k = 0 with E = +U (twice) arising from the surface states, and from the dimers: , and E = −U − γ 1 ((N − 1)/2 times).In addition, there is the dimer consisting of the stacking fault giving energies E = ± U 2 + γ 2 1 .Again, the surface bands appear at lower energy (E = ±U ) than the bulk bands for 0 ≤ U < γ 1 /2, and the bulk bands touch at k = 0 and U = γ 1 .Near low energy, the band structure is similar to that of the faultless even N case.Considering the bands for nonzero k values, there is always a band within the nominal band gap for U < γ 1 , whereas the bandgap, for U > γ 1 , is given by E g = 2(U − γ 1 ).As with the even N case, the states localized on the soliton are not near low energy, as indicated by the arrows in Fig. 6 (bottom row).They are at energies E = ± U 2 + γ 2  1 for k = 0, and then move further away from zero for nonzero k.

IV. BERNAL STACKING
A. ABA

Faultless system
For bulk, faultless ABA stacking, the unit cell consists of four atoms such as A 1 , B 1 , A 2 , B 2 (in Fig. 1(c)), and the Hamiltonian may be written as using the approximation for the intralayer hopping (8) with k < k c , where q z is the z component of the wave vector and d is the layer separation (the length of the unit cell in the z direction is 2d).The four corresponding bulk N=16 faultless N=16 (8,8) fault N=15 faultless N=15 (8,7) fault ABC' stacking FIG. 6. Low-energy band structures for ABC ′ stacking.Each column shows a different value of the onsite energies U , each row shows a different system.The top row shows a faultless system with N = 16 layers, whereas the second row shows N = 16 layers with a sharp fault at its center (8,8).The third row shows a faultless system with N = 15 layers, and the final row shows N = 15 layers with a sharp fault near its center (8,7).Horizontal arrows show the energies at kx = 0 of the states localized on the soliton.Note that the axis scales for the final column (with U = 2.00γ1) differ from the rest.In all plots, parameter values are γ0 = 3.16 eV, γ1 = 0.381 eV [40], a = 2.46 Å [41].
bands are given by which shows that the band gap is E g = 2U .The two conduction bands are degenerate at the Brillouin zone edge (q z d = ±π/2) because cos q z d = 0 there giving a block diagonal Hamiltonian (30), likewise the two valence bands.
For a thin film of finite thickness, a faultless system is gapped with bandgap E g ≈ 2U , Fig. 1(c).This can be understood by considering k = 0 where the intralayer hopping is zero and the system separates into isolated atoms with energies E = ±U plus a N -mer which is a single chain of N atoms (N = 8 in Fig. 1(c)).This chain corresponds to the CDW model [15,16,[33][34][35][36] which is a one-dimensional chain with constant hopping γ 1 and alternating onsite energies ±U , having a band gap of 2U , as described in the methodology, Section II.

Single soliton
For a system with a single fault, there are two bands localized on the fault, one of which generally lies within the bulk band gap, Fig. 1(i) and Fig. 7.The nature of this band depends on whether the fault occurs after an even or odd number of layers [47] as shown in Fig. 7.When the fault occurs after an even number of layers, Fig. 1(i), Fig. 7(a), and Fig. 7(b), the fault consists of two consecutive onsite energies of +U (on atoms A 4 and B 5 in Fig. 1(i)) and we refer to this as an antisoliton; the energy band moves into the conduction band for k ≫ k c .However, when the fault occurs after an odd number of layers, Fig. 7(c) and Fig. 7(d), the fault consists of two consecutive onsite energies of −U , a soliton texture, and the energy band moves into the valence band for k ≫ k c .
We consider k = 0 where the intralayer hopping is zero.Once again, the system separates into isolated atoms with energies E = ±U plus a N -mer (N = 8 in Fig. 1(i)).But now the N -mer, which corresponds to the ABA stacking Fault at center (8,8) CDW model [15,16,[33][34][35][36], has a soliton or antisoliton in it [34,36].It was shown [36] that the energy of the state localized on such a fault lies within the bulk band gap as long as 0 < γ 1 /U < N , when the fault lies at the center of the system of N layers.
For large U (U ≫ γ 1 and U ≫ ℏvk), it is possible to estimate the energy of the soliton level using perturbation theory in the hopping strength.Using the analytic approximation (8) for the intralayer hopping matrix element gives for {γ 1 , ℏvk} ≪ U , assuming the soliton is not at the surface of the system.For an antisoliton, E asol (k) = −E sol (k).The approximation (31) is shown as the dashed line in Fig. 7(b) and Fig. 7(d).
Soliton states are extended states in the in-plane direction, but are localized in the out-of-plane (z) direction.For k = 0, the soliton state is localized on the soliton [34,36], and, for k x = k c (and k y = 0), the state is still localized near the soliton as shown in Fig. 4(c).Since the system is equivalent to the CDW model at k = 0, the soliton state is polarized, with maximum polarization at k = 0, Fig. 4(d).In this case, the equivalence with the CDW model only applies to the N -mer (the atomic sites connected by interlayer hopping γ 1 ), excluding the isolated non-N -mer atoms.Therefore, we modify the definition of the translation operator T a/2 , Eq. ( 13), in order to describe translation between the atoms on the N -mer only, which is shown for an even number of layers N (for N ≥ 4).The translation Ta/2 for an odd number of layers is the same with the replacements ( Ta/2 ) N,1 = ( Ta/2 ) N −1,2 = 0 and ( Ta/2 ) N −1,1 = ( Ta/2 ) N,2 = 1.
Fig. 4(d) shows the polarization py , Eq. ( 34), as a function of k x (with k y = 0) for U/γ 1 = 0.6 (solid line) and U/γ 1 = 1.8 (dashed line).The polarization is a maximum at k = 0, and, at k = 0, it is larger for smaller U/γ 1 (solid line).Note that the values of py at k = 0 in Fig. 4(d) are approximately the same as the values of p y at k x = k c in Fig. 4(b) because, at these points, both systems are equivalent to the CDW model.The polarization py never reaches unity because the nonsymmorphic chiral symmetry is broken by the ends of the system and the finite width of the soliton.

Soliton-antisoliton pair
For ABA stacking, the band structure and density of states of a soliton-antisoliton pair (5,6,5) with N = 16 layers is shown in Fig. 8 for different values of U .There are two localized states associated with the soliton and two with the antisoliton, of which one soliton state and one antisoliton lie within the bulk bandgap of E g = 2U .
For small values of U/γ 1 , the low-energy soliton and antisoliton states hybridize with an anticrossing, Fig. 8(a) and Fig. 8(g).For larger values of U/γ 1 , the soliton and antisoliton bands separate.For large U (U ≫ γ 1 and U ≫ ℏvk), we can use the perturbation expression (31) for the energy of a soliton to estimate the separation of these bands as 2U − 2γ 1 + γ 2 1 /U .This separation appears as an extra feature within the bulk band gap in the density of states, Fig. 8 (bottom row).

Faultless system
For bulk, faultless ABA ′ stacking, using the approximation for the intralayer hopping (8) with k < k c , we can use a unit cell of two atoms, either A i , B i on an odd layer or B i , A i on an even layer, to give the Hamiltonian N=16 faultless N=16 (8,8) fault ABA' stacking where q z is the z component of the wave vector and d is the layer separation.The two corresponding bulk bands are This shows that there is no band gap for U ≤ γ 1 , but that the band gap is For a thin film of finite thickness, we begin by considering the faultless system, Fig. 9 (top row).At k = 0, the system separates into isolated atoms with energy +U [47] and a N -mer monoatomic chain with onsite energies −U and nearest neighbor hopping γ 1 , Fig. 1(d).Such a chain has a bulk band given by E(q z ) = −U + 2γ 1 cos(q z d) where d is the layer separation and q z is the z component of the wave vector; this band is centered on energy −U and it has a bandwidth of 4γ 1 .Thus, this band overlaps with the isolated atoms of energy +U for U ≤ γ 1 .Hence, Fig. 9(c) for U = γ 1 shows the bands due to the isolated atoms of energy +U just touching the lower bands due to the N -mer chain whereas Fig. 9(d) for U > γ 1 shows that they are separated.

Single soliton
Band structures for ABA ′ stacking with a sharp fault at its center are shown in Fig. 9 (bottom row).At k = 0, the system separates into isolated atoms with energies ±U plus the N -mer chain.Now the N -mer chain contains a fault, so it consists of a monoatomic chain with onsite energies −U and nearest neighbor hopping γ 1 connected to a monoatomic chain with onsite energies +U and nearest neighbor hopping γ 1 .The former has a bulk band E(q z ) = −U + 2γ 1 cos(q z d), the latter has a bulk band E(q z ) = U + 2γ 1 cos(q z d), i.e., the system consists of one part with a band centered on energy +U and of bandwidth 4γ 1 , connected to another part with a band centered on energy −U and of bandwidth 4γ 1 .For U ≤ 2γ 1 , the two bands overlap leading, in a system of a finite number of layers, to overlapping bands with anticrossings at low energy as in Fig. 9(h) to Fig. 9(j).For U > 2γ 1 , the two bands do not overlap, giving an overall band gap as in Fig. 9(l).
Although the presence of a fault for ABA ′ stacking creates changes in the band structure, Fig. 9, there are no states localized on the fault.This can be understood by considering the system at k = 0: The fault divides two sections of monoatomic chains with different onsite energies (+U or −U ) and, thus, the two sections are not degenerate, nor do they have a Dirac-like low-energy continuum description [25].
V. AA STACKING A. AA

Faultless system
For bulk, faultless AA stacking, the unit cell consists of two atoms A i , B i on the ith layer, and the Hamiltonian may be written as using the approximation for the intralayer hopping (8) with k < k c , where q z is the z component of the wave vector and d is the layer separation.The two corresponding bulk bands are At k = 0, the bands are centered on energy +U or −U and each of them has a bandwidth of 4γ 1 .This shows that there is no band gap for U ≤ 2γ 1 , but that the band gap is For a thin film of finite thickness, band structures for the faultless system are plotted in Fig. 10 (top row).The two bands overlap for U ≤ 2γ 1 as shown in Fig. 10 (top row).It is also noteable that there are no anticrossings when the bands overlap.This is because the system has reflection symmetry (swap A 1 with A N , B 1 with B N , etc.) which may be used to block diagonalize the position space Hamiltonian (1) into two separate blocks with either even or odd parity eigenstates.

Single soliton
In the presence of an atomically-sharp fault, there are different ladders either side of the fault, Fig. 1(k), resulting in band structures as plotted in Fig. 10 (bottom row).As with the faultless case, there are two bands overall, centered on energy +U or −U and each of them with a bandwidth of 4γ 1 .Thus, the two bands overlap for U ≤ 2γ 1 as shown in Fig. 10 (bottom row).Unlike the faultless case, however, there are some anticrossings (and crossings) of the bands in the region when the bands overlap.The situation shown in Fig. 10 (bottom row) is somewhat special because the fault is at the center so that the system has inversion symmetry (swap A 1 with B N , B 1 with A N , etc.) and this may be used to block diagonalize the position space Hamiltonian (1) into two separate blocks with either even or odd parity eigenstates.Hence there are level crossings and all bands are doubly degenerate at k = 0.However, within each even or odd parity block, anticrossings can arise giving the band structure in Fig. 10(h) and Fig. 10(i).
As with ABA ′ stacking, although the presence of a fault for AA stacking creates changes in the band structure, Fig. 10, there are no states localized on the fault.The term describing hopping along the ladder, 2γ 1 cos(q z d), appears on the main diagonal in the ladder Hamiltonian (36) and it breaks chiral symmetry.The trivial nature of the fault may be further understood by considering the limit k = 0 when the system separates into monoatomic chains which do not have a Dirac-like low-energy continuum description [25].

Faultless system
For bulk, faultless AA ′ stacking, using the approximation for the intralayer hopping (8) with k < k c , we can use a unit cell of two atoms, either A i , B i on an odd layer or B i , A i on an even layer, to give the Hamiltonian where q z is the z component of the wave vector and d is the layer separation.The two corresponding bulk bands are which have a band gap E g = 2U .
For a thin film of finite thickness, there is indeed a band gap of E g ≈ 2U in the numerically-derived band structure, Fig. 1(f).Note that the system at k = 0 consists of two separated CDW chains and, for k ̸ = 0, the chains are coupled by intralayer hopping to form a ladder.By dimensional reduction, treating k as a fixed parameter, this ladder is equivalent to a one-dimensional model in class CI [16] which is a topologically trivial insulator.

Single soliton
In the presence of a single atomically-sharp fault, there are four states localized on the fault, Fig. 1(l) and Fig. 11.Of these, one is at the top of the conduction band, one is at the bottom of the valence band, and the other two generally lie within the band gap.At low U values, they can hybridize together as in Fig. 11(b) and Fig. 11(g).
At k = 0, the system consists of two separated CDW chains, each of them has either a soliton or an antisoliton.This is why there are two low-energy states.The analogy with the CDW model means that, at k = 0, we can say that the energies of the two low-energy states are within the band gap for 0 < γ 1 /U < N for a fault at the center of a system with N layers.For U/γ 1 ≫ 1, these energies are at k = 0, as indicated by the two middle arrows in Fig. 11(e).The solitons give rise to a step-like nonzero density of states within the bulk band gap.For U ≫ γ 1 , the nonzero density of states due to the solitons is visible at energies −U to −(U − γ 1 + γ 2 1 /(2U )) and at (U − γ 1 + γ 2  1 /(2U )) to U , Eq. ( 38), see Fig. 11(j).In addition to the soliton feature, there are also sharp peaks in the density of states at the bulk band edges, E = ±U .
For small values of U/γ 1 , there are crossings between levels, including the two low-energy levels, Fig. 11(b), whereby the energies appear to oscillate as a function of k.Using spatial inversion symmetry, it is possible to block diagonalize the Hamiltonian (1) into a block with even parity states and a block with odd parity states.The odd parity block has eigenstates which are the electron-hole reflection (E → −E) of the even parity eigenstates.Each individual block describes a system similar to AA ′ stacking with a spectrum having anticrossings at U ̸ = 0 and a single soliton level.For small U/γ 1 , the soliton level crosses zero energy, at which point its electron-hole partner (from the opposite parity block) crosses zero energy in the opposite direction.Oscillations are caused by anticrossings with levels from the same block.The anticrossings, and resulting oscillations, are reduced as U/γ 1 increases because the other levels (in the bulk conduction and valence bands) move further away in energy and the soliton-antisoliton pair move away from zero energy and separate.This picture holds in the presence of spatial inversion symmetry.When the soliton is off-center, the band structure is similar, but the level crossings between the levels with even or odd parity are replaced by anticrossings.

VI. SMOOTH SOLITONS
Atomically sharp solitons support localized states within the bulk band gap for ABC, ABA and AA ′ stacking, Fig. 1.Now we will describe the properties of smooth solitons with a finite width for these three stacking types.Band structures and density of states are determined as ABC stacking (8,8) U=0.25� 1 (a) 12. Band structures and density of states g(E) for smooth solitons with width ζ = 8 at the center (8, 8) of a system with N = 16 layers for ABA, ABC, and AA ′ stacking.In the band structure plots, horizontal arrows show the energies at kx = 0 of the states localized on the soliton.There is a different axis scale for AA ′ stacking.In all plots, parameter values are γ0 = 3.16 eV, γ1 = 0.381 eV [40], a = 2.46 Å [41].For the band structures, ky = 0, and, for the density of states, δ = 0.01γ1.described in the methodology, Section II, by numerically diagonalizing Hamiltonian (1) after replacing the diagonal elements with a smooth soliton texture.Specifically, for ABC and ABA stacking, we model the onsite energy of an A site on the jth layer, where j = 1, 2, . . ., N , as where U is the magnitude of the texture at infinity, m is the number of layers below the fault, and ζ is the soliton width in dimensionless units, i.e., measured in units of the interlayer spacing.For AA ′ stacking, For all solitons, we assume charge neutrality within each layer so that, for the onsite energy of a B site on the jth layer, ϵ B,j = −ϵ A,j .Band structures and density of states for the three stacking types are shown in Figure 12 for a smooth soliton of width ζ = 8 at the center (8, 8) of a system with N = 16 layers, and different values of U/γ 1 .

A. ABC stacking
For ABC stacking, Fig. 12 shows that the soliton level remains near the same energy E sol ≈ γ 1 at k = 0 for a range of U values, in contrast to a sharp soliton, Eq. (17).Since the soliton level moves into the valence band for large k, this means that the soliton state crosses the bulk band gap, leading to a nonzero density of states.
The behavior at k = 0 may be understand because the system separates into dimers and isolated atoms there, and the soliton state is localized on a dimer consisting of site B m and A m+1 connected by interlayer hopping γ 1 , as in Fig. 1(g).The onsite energy of these sites is identical, ϵ A,m+1 = ϵ B,m = −U tanh(1/2ζ) according to Eq. ( 39).The dimer yields two energy levels and the soliton level is the highest energy of the pair, This is in agreement with the sharp soliton energy ( 17), E sol = γ 1 − U , in the limit ζ → 0, as expected.For ζ ≫ 1/2, we expand the hyperbolic tangent to give E sol ≈ γ 1 − U/(2ζ) which shows that E sol ≈ γ 1 for a range of U values if the width ζ is large enough.As in the case of a sharp soliton, the two surface states, on A 1 and B N , are completely disconnected at k = 0 and they give two degenerate states at energy E ≈ U according to Eq. ( 39) (they are at E = U in a sufficiently large system N ≫ ζ).For N ≫ ζ ≫ 1, the soliton level moves from E sol ≈ γ 1 at k = 0 to the valence band at large k, crossing the surface states at E ≈ U provided that U < γ 1 .
At k = k c , the soliton level is at E sol ≈ 0 for a wide range of U values.This is because the system is approximately equivalent to the CDW model at this point and, for smooth solitons ζ ≫ 1, the system approaches the continuum limit which supports a zero-energy soliton level, as discussed in detail previously for the CDW model [34,36].
For U < γ 1 , all levels in the conduction and valence bands originate from near E = ±γ 1 at k = 0, as for rhombohedral graphite at U = 0.For U > γ 1 , the levels no longer coalesce at k = 0, but bifurcations are visible in the conduction and valence bands whereby levels that are doubly degenerate at k = 0 split at nonzero k.This is because of spatial inversion of the texture of onsite energies (39) about the center of the soliton.At k = 0, there are identical dimers either side of the soliton and equidistant from it, creating a degeneracy which is broken at nonzero k.
As stated above, the soliton state at k = 0 is localized on two adjacent sites near the center of the soliton as For k = k c , it remains localized on the soliton, but has a broader extent in position space as shown in Fig. 13(a) for U/γ 1 = 0.6.Fig. 13(b) shows the polarization p y = ⟨ψ|T a/2 S z |ψ⟩, Eq. ( 14), as a function of k x (with k y = 0) for U/γ 1 = 0.6 (solid line) and U/γ 1 = 1.8 (dashed line).There is a noteable dip in p y for U/γ 1 = 0.6 (solid line) at k x ≈ 0.6k c which arises when the soliton level crosses the surface states, Fig. 12(g).Otherwise, p y ≈ 1 for a range of k values.This is related to the fact that, for k = k c , the Hamiltonian ( 16) is approximately equal to that of the CDW model.However, the polarization p y is never exactly one because the nonsymmorphic chiral symmetry is broken by the ends of the system and the finite width of the soliton [34,36].

B. ABA stacking
For ABA stacking, Fig. 12 shows that the antisoliton level remains near the same energy E sol ≈ 0 at k = 0 for a wide range of U values, and it moves into the conduction band at large k.The level is near zero at k = 0 because the central N -mer of the system is separated from the individual non-dimer atoms at this point and, thus, the N -mer is equivalent to the CDW model.For smooth solitons ζ ≫ 1, the N -mer approaches the continuum limit with a soliton level approaching zero exponentially quickly as a function of ζ, as discussed in detail previously for the CDW model [34,36].
In addition to the soliton state, however, there are many levels in the conduction and valence bands, and also in the nominal bulk band gap, which bifurcate, being doubly degenerate at k = 0 and splitting at nonzero k.These levels arise from the individual non-N -mer atoms which are separated from the rest of the system at k = 0. Owing to spatial inversion of the texture of onsite energies (39) about the center of the soliton, there are single atoms with identical onsite energies either side of the soliton and equidistant from it, creating a degeneracy which is broken at nonzero k.The two atoms nearest to the center of the soliton have energies E = −U tanh(1/2ζ) at k = 0 according to Eq. ( 39), so, for ζ ≫ 1, they are very close to zero energy as E ≈ −U/(2ζ).
For k = 0, the soliton state is localized near the center of the soliton [34,36], and, for k x = k c (and k y = 0), the state is still localized near the soliton center as shown in Fig. 13(c).Since the system is equivalent to the CDW model at k = 0, the soliton state is polarized, with maximum polarization approaching one at k = 0, Fig. 13(d).Again, the polarization py is never exactly one because the nonsymmorphic chiral symmetry is broken by the ends of the system and the finite width of the soliton [34,36].

C. AA ′ stacking
For AA ′ stacking, Fig. 12 shows that there are a pair of bands (soliton and antisoliton) near zero energy, as for a sharp soliton, but the major difference is that these two bands do not separate as U/γ 1 increases, but remain near zero energy as flat bands with a large extent up to k ≈ 2k c .
At k = 0, the system consists of two separated CDW chains, each of them has either a soliton or an antisoliton.For smooth solitons ζ ≫ 1, the CDW chains approach the continuum limit with a soliton (or antisoliton) level approaching zero exponentially quickly as a function of ζ, as discussed previously [34,36].
For small values of U/γ 1 , Fig. 12(e) and (k), there are also level crossings between the soliton and antisoliton bands which appear to oscillate as a function of k, as described for sharp solitons in Section V B 2. The oscillations, due to anticrossings with the other conduction and valence bands, are reduced as U/γ 1 increases and the other bands move away from zero energy, separating from the soliton-antisoliton pair.The flat bands due to the soliton-antisoliton pair result in a large peak in the density of states near zero energy, Figure 12(r) and  (x).We consider the effect of disorder on this peak in Section VIII.

VII. SINGLE LAYER DEFECTS
So far, we have considered solitons which consist of a change in the texture of onsite energies without changing the interatomic hopping, and they may be sharp on the atomic scale or smooth with a nonzero width ζ.There are many other types of stacking faults one could consider.An example, which also involves onsite energies, is to reverse the signs of the onsite energies on a single layer only.This doesn't change the texture of the onsite energies either side of the defect, and we refer to this as a 'single layer defect', not a soliton.Such defects were modeled for AA, AA ′ , ABA and ABA ′ stacking using density functional theory (DFT) in the context of hexagonal boron nitride [38], and we will only discuss them briefly here, for ABC, ABA and AA ′ stacking.The lattice structures, with the defects, are shown schematically in the top panels of Fig. 14 for a system we denote as (3, 1, 4) indicating a system of N = 8 layers consisting of 3 layers (at the bottom) with regular stacking, followed by a single layer defect, followed by 4 layers (at the top) with regular stacking.Band structures and density of states for ABC, ABA, and AA ′ stacking, are shown in the remaining panels in Fig. 14 for a system with a larger number of layers N = 16 and a defect near its center (7,1,8).

A. ABC stacking
For ABC stacking, Fig. 14 shows that there are generally two defect states within the band gap, one merges with the conduction band for k ≫ k c , the other merges with the valence band.For U < γ 1 /2, Fig. 14(a), the two surface states are closer to zero energy than the defect states whereas, for U > γ 1 /2, the defect states are closest to zero energy and, for U = γ 1 , Fig. 14(m), they touch at E = 0 and there is no band gap at this point.For U = γ 1 /2, Fig. 14(g), the defect states hybridize strongly with the surface states.
The behavior at k = 0 may be understand because the system separates into dimers and isolated atoms there.There are the two isolated surface states with energies E = ±U , N − 3 dimers each with energies E = ± U 2 + γ 2 1 , and, near the defect, a dimer with onsite energies −U giving E = −U ± γ 1 and a dimer with onsite energies +U giving E = +U ± γ 1 .These four energies are shown with the horizontal arrows in Fig. 14.The presence of two defect states within the band gap is in stark contrast to a single soliton which only supports one state within the gap, Fig. 2. Instead, a single layer defect has a band structure very similar to that of a soliton-antisoliton pair, Fig. 5, because the latter can also be viewed as a defect, albeit of larger spatial extent than just one layer, which doesn't change the texture of the onsite energies either side of the defect.

B. ABA stacking
For ABA stacking, Fig. 14 shows that there is a single defect state within the bulk band gap.For small U/γ 1 , it is near the bottom of the conduction band for k = 0, crossing the band gap to move into the valence band at ABC stacking  14).The level is the one within the band gap.A single soliton also creates a single level within the band gap, Fig. 7, and the band structure and density of states of a single soliton and a defect are similar, particularly for U ≫ γ 1 .For U ≪ γ 1 , the defect state extends across most of the band gap, creating a wide energy region with a nonzero density of states, unlike the soliton state.

C. AA ′ stacking
For AA ′ stacking, Fig. 14 shows that there are generally two defect states within the band gap, one merges with the conduction band for k ≫ k c , the other merges with the valence band.The lattice structure, Fig. 14, shows the ladder structure consisting of two coupled Nmers.In one of them, there is a trimer of states with onsite energy E = +U near the defect, and, in the other, there is a trimer of states with onsite energy E = −U near the defect.Together, they contribute six states which are, at k = 0 and , and E = U ± √ 2γ 1 (as indicated by the horizontal arrows in Fig. 14).The two levels E = −U + √ 2γ 1 and E = U − √ 2γ 1 are within the band gap.A single soliton also creates two levels within the band gap, Fig. 11, and the band structure and density of states of a single soliton and a defect are similar.There are differences such as the hybridization of the two levels with the levels being close together and strongly hybridized for U ≪ γ 1 in the soliton and for U ≈ γ 1 in the defect.

VIII. THE ROLE OF INTERLAYER DISORDER
We have shown that solitons, both atomically-sharp and smooth in position space, are able to support localized states with energies within the bulk band gap for ABC, ABA and AA ′ stacking.The soliton bands are generally dispersive as a function of the in-plane wave vector k so that they give rise to a non-zero density of states without any particularly sharp features.The exception are smooth solitons for AA ′ stacking which give flat bands yielding a narrow peak at zero energy in the density of states, Fig. 12(x).
For certain values of the in-plane wave vector k, the lattice structure of these systems may be related by dimensional reduction to that of the CDW model [15,16,[33][34][35][36] which has nonsymmorphic chiral symmetry.The influence of symmetry-breaking disorder has been considered previously for the CDW model [36] so here we consider it only for the most interesting case of AA ′ stacking, for both atomically-sharp and smooth solitons.
To focus on the topology related to the CDW model, we consider random tight-binding parameters in the perpendicular-to-layer direction while preserving translational invariance in the in-plane direction.We consider four types of disorder: (i) "Sharp onsite disorder" (diagonal) which gives an additional contribution δϵ A,j to the onsite energy of an A site on the jth layer, where j = 1, 2, . . ., N , which is drawn randomly from a uniform distribution −W ≤ δϵ A,j ≤ W with disorder strength W .We maintain charge neutrality within each layer so that δϵ B,j = −δϵ A,j for all layers j. (ii) "Sharp hopping disorder" (off-diagonal) where the interlayer coupling takes random values γ 1 + δ n where n = 1, 2, . . ., N − 1 indexes the N − 1 interlayer spaces and each δ n takes a value drawn randomly from a uniform distribution −W ≤ δ n ≤ W .For AA ′ stacking, the values of coupling between A atoms and between B atoms for any pair of adjacent layers are identical, so that there are only N − 1 random values of δ n in total.(iii) "Smooth onsite disorder" (diagonal) with an additional contribution δϵ A,j to the onsite energy of an A site on the jth layer described by a Gaussian-correlated potential [49,50] as given by where η is the correlation length in dimensionless units, i.e., measured in units of the interlayer spacing.The summation is over all layers m = 1, 2, . . ., N with w m drawn randomly from a uniform distribution −W ≤ w m ≤ W with disorder strength W .We maintain charge neutrality within each layer so that δϵ B,j = −δϵ A,j for all layers j. (iv) "Smooth hopping disorder" (off-diagonal) where the interlayer coupling takes random values γ 1 +δ n where n = 1, 2, . . ., N −1 indexes the N −1 interlayer spaces and each δ n is described by a Gaussian-correlated weighting as in Eq. ( 42), the only difference being that now there are N − 1 bonds and N − 1 independent parameters w m .
Again, the values of coupling between A atoms and between B atoms for any pair of adjacent layers are identical, so that there are only N − 1 random values of δ n in total.
The normalisation with the square root factor in Eq. ( 42) is used so that smooth disorder interpolates between sharp disorder for η ≪ 1 and sample-to-sample parameter variations [36] for η ≫ N .The latter are parameter values that are spatially uniform across a single member of the ensemble, i.e., δϵ A,j = δϵ A for all j, but that differ from sample to sample with δϵ A drawn randomly from a uniform distribution −W ≤ δϵ A ≤ W . Smooth solitons in the CDW model are fairly robust to parameter variations [36] because the nonsymmorphic chiral symmetry holds in the continuum limit (of a smooth soliton with parameter variations).However, the equivalence of AA ′ stacking with the CDW model only holds at k = 0, so smooth solitons for AA ′ stacking are not expected to be as robust.As parameter variations were studied in detail for the CDW model [36], we consider instead the case of a smooth potential with a finite range 1 ≲ η ≲ N .
For AA ′ stacking and an atomically sharp soliton in the absence of disorder, Fig 11(j), the density of states has a step-like nonzero density of states within the bulk band gap −U ≤ E ≤ U due to the soliton bands, Fig 11(e).In addition, there are also sharp peaks in the density of states at the bulk band edges, E = ±U .Figure 15 shows the disorder-averaged density of states ⟨g(E)⟩ for an atomically-sharp soliton at the center (8,8) of a system with AA ′ stacking and N = 16 layers.We consider the example of onsite energy U = 2.0γ 1 and disorder strength W = 0.5γ 1 .Figure 15 shows that all types of disorder tend to smooth out the step-like feature due to the soliton states within the bulk band gap.The sharp peaks at E = ±U at the band edges are quite robust to disorder, as they are only smoothed out by atomicallysharp onsite disorder, Figure 15(a).For AA ′ stacking and a smooth soliton in the absence of disorder, Fig 12(x), the density of states shows a narrow peak at zero energy due to the flat bands related to the soliton states.Figure 16 shows the disorder-averaged density of states ⟨g(E)⟩ for a smooth soliton at the center (8, 8) of a system with AA ′ stacking and N = 16 layers.Sharp disorder, Fig. 16(a) and (b), tends to destroy the narrow peak at zero energy, leaving a small but finite density of states in the bulk band gap.The narrow peak is fairly robust to smooth disorder, however, and more robust to hopping disorder than onsite disorder, Fig. 16(c) and (d).This is possibly because the solitons are textures in the onsite energies, not the hopping parameters.

IX. ELECTRONIC TRANSPORT
Solitons, both atomically-sharp and smooth in position space, are able to support localized states with energies within the bulk band gap for ABC, ABA and AA ′ stacking.The contrasting electronic band structures and densities of states will impact transport and spectroscopic measurements [31,[51][52][53][54][55][56][57].As an example, we model coherent, ballistic, electronic transport for rhombohedral-stacking, generalizing a model developed for monolayer and bilayer graphene [42,58] in which the central sample is connected to leads of the same material that are heavily doped in order to provide a large density of states.The conductivity is determined using the Landauer-Büttiker formalism [59] where the transmission probability is found by wave-matching at the boundaries of the leads to the central sample, and Hamiltonians are taken in the continuum limit.Details of our numerical calculations are described in Appendix A.
Generalizing previous results for monolayer and bilayer graphene [42,58], we find that the minimal conductivity for rhombohedral stacking in the absence of onsite energies (U = 0) is given by σ min = 4N e 2 /(πh) for N layers with W ≳ L where W and L are the sample width and length, respectively.Figure 17 shows band structures and corresponding conductivities for a rhombohedrally-stacked system with N = 10 layers and onsite energy U = 0.1 eV.The band structures are plotted near the K point, with the direction of the Γ and M points indicated for negative and positive k x , respectively.Fig. 17

A. ABC stacking: rhombohedral graphite
We consider a tight-binding model of rhombohedral graphite containing skew interlayer hopping γ 3 and γ 4 , and next-nearest layer hopping γ 2 .In a basis of a single orbital on each atomic site (A 1 , B 1 , A 2 , B 2 , . . ., A N , B N ), the Hamiltonian (1) is modified [63] as where intralayer blocks D i are defined in Eq. ( 2) and we use U i = U for all i (for a faultless system).Interlayer blocks are and next-nearest layer blocks are The band structures and corresponding conductivities for this model are shown in Fig. 17   g) and (h) are for an atomicallysharp soliton at the center of the system.In comparison to the minimal model, it can be seen that the additional parameters have a generally small effect on the band structure and conductivity, without producing any qualitative changes.
B. AA ′ stacking: hexagonal boron nitride As a second example of the influence of additional tight-binding parameters, we use a tight-binding model of AA ′ -stacked hexagonal boron nitride (h-BN) with parameters fitted to calculations using density functional theory (DFT) and the GW approximation for a faultless system [64].In this model [64], tight-binding parameters are γ 0 = 2.33 eV, γ 1 = 0.5 eV, and U = 3.625 eV.Additional parameters, as compared to the minimal model, are in-plane next-nearest neighbor hopping γ n = −0.4eV and next-nearest layer hopping γ 2 = −0.1 eV.
In this tight-binding model, for AA ′ -stacking, the Hamiltonian (1) is modified [64] as where U i = U for odd i and U i = −U for even i (for a faultless system).Interlayer blocks are The band structure plotted across the whole Brillouin zone is shown for a faultless system in Fig. 18(a) and for an atomically-sharp soliton at the center in Fig. 18(b).For the faultless system, the band gap is near the K point but is indirect in the presence of additional parameters, Fig. 18 C. ABA stacking: hexagonal boron phosphide As a third example of the influence of additional tightbinding parameters, we use a tight-binding model of ABA-stacked hexagonal boron phosphide (h-BP) with parameters fitted to calculations using DFT for a faultless system [61].As compared to the minimal model, we add skew interlayer couplings γ 3 and γ 4 .In this tightbinding model, for ABA-stacking, the Hamiltonian (1) is modified [65] as where intralayer blocks D i are defined in Eq. (2) and we use U i = U for all i (for a faultless system).Interlayer blocks V are given by Eq. (44).We use parameter values fit to DFT [61] near the K point for bilayers, neglecting small differences in the hopping between boron and phosphorous atoms, giving γ 0 = 1.65 eV, γ 1 = 0.761 eV, γ 3 = −0.456eV, γ 4 = 0.255 eV, and U = 0.636 eV.The band structure plotted across the whole Brillouin zone is shown for a faultless system in Fig. 18(c) and for an atomically-sharp soliton at the center in Fig. 18(d).For the faultless system, the band gap is near the K point but is indirect in the presence of additional parameters, Fig. 18(c).For the soliton, there is a state localized on the soliton within the bulk band gap, above the bulk valence bands, Fig. 18(d).This is qualitatively similar to the minimal model, Fig. 7(c).Note that h-BP is in a regime of a relatively-small band gap with U/γ 1 ≈ 0.84.

XI. RELEVANCE TO SPECIFIC MATERIALS
Although graphite occurs with both Bernal and rhombohedral stacking, our results do not apply to it because it is practically impossible to induce different onsite energies in a controlled way.However, we note that it may also be possible to obtain isolated bands in graphite by considering rhombohedrally-stacked sections within Bernal-stacked graphite [66,67], by applying a large displacement field to a stacking fault near the surface of rhombohedral graphite [8], or with a large in-plane magnetic field in rhombohedral graphite [68].Localized topological states have also been considered in a different context, that of domain walls in the horizontal, in-plane direction [69,70].
The stacking types to be realized experimentally (AA ′ , ABA and ABC) all have different atoms connected in the vertical direction, i.e., boron to nitrogen, and a fault consisting of a connection between two identical atoms, i.e., boron to boron or nitrogen to nitrogen, will be energetically expensive.Single layer defects, as described in Section VII, were modeled using DFT [38] for AA, AA ′ , ABA and ABA ′ stacking, and it was found that the formation energy of such faults is very small.For example, an AAtype defect in AA ′ stacking was estimated to cost of the order of 50 meV.Even so, it is likely that such a defect is difficult to realize experimentally because a relative shift in the in-plane direction would yield the more favorable ABA stacking [29,78].Nevertheless, experiments [90,91] have recently demonstrated the possiblilty of engineering stacking of boron nitride (to create ferroelectric materi-als), and a relative twist between layers creates domains of ABA stacking separated, in the in-plane direction, by domain walls and regions of AA stacking.

XII. CONCLUSIONS
Using a minimal tight-binding model, we determined the electronic properties of thin films of binary compounds with two atoms per unit cell arranged as stacked two-dimensional honeycomb lattices.The two atoms per cell are assigned different onsite energies.We considered six different stacking orders to determine whether a fault in the texture of onsite energies in the stacking direction supports localized states.Faults for ABC, ABA, and AA ′ stacking support localized states within the band gap, whereas faults for ABC ′ , ABA ′ , and AA stacking do not.For ABC and ABA stacking, there is a single localized state within the bulk band gap and, for ABC stacking, this state may hybridize with surface states.For AA ′ stacking, there are two states within the bulk band gap which may hybridize with each other, depending on parameter values.
We consider smooth solitons where the texture of onsite energies changes over a length scale greater than the interlayer spacing, leading to different band structures as compared to atomically sharp solitons.In particular, a smooth soliton in AA ′ stacking results in a flat band at zero energy and a corresponding narrow peak in the density of states.We show that this feature is fairly robust to long-range correlated disorder.Finally, we also consider the band structure due to single layer defects [38] where the signs of the onsite energies are reversed on a single layer only without a subsequent change of texture.Overall, different stackings and types of fault produce a range of contrasting electronic band structures and densities of states which will manifest themselves in various transport and spectroscopic measurements [31,[51][52][53][54][55][56][57].
dispersion is approximately linear so that U ∞ ≈ ℏvQ (R) ℓ .For our numerical calculations, Fig. 17, we choose U ∞ = 50γ 1 so that Q The unknown coefficients are found by determining C = M −1 A. Then, for a given energy value, the con-ductivity [42,59] is given by where L is the length and W is the width of the sample.The trace may evaluated by summing over the eigenvalues of tt † for all k y values where, for periodic boundary conditions [42], k y = 2πn/W for integer n = 0, ±1, ±2, . ...

FIG. 1 .
FIG. 1. Summary of the main results for the six different stacking types.The first row (a)-(f) shows the faultless systems with a schematic side view of the stacking for N = 8 layers with red (yellow) circles indicating A (B) atoms, horizontal (vertical) solid lines indicating nearest-neighbor intralayer (interlayer) hopping.The bottom of each panel shows the corresponding band structure near low energy as a function of wave vector component kx (with ky = 0) measured from the center of the K+ valley, where kc = γ1/(ℏv).The second row (g)-(l) shows the lattices for N = 8 layers with an atomically-sharp fault in the onsite energies occurring between the middle two layers and the bottom of each panel shows the corresponding band structures.Horizontal arrows show the energies at kx = 0 of the states localized on the soliton.In all plots, parameter values are γ0 = 3.16 eV, γ1 = 0.381 eV [40], a = 2.46 Å [41], and the magnitude of the alternating onsite energies is u = U/γ1 = 0.4.

FIG. 3 .
FIG.3.ABC stacking with N = 16 layers and a sharp soliton at the center(8,8).Plots of the analytic band energies of bands hybridized between the soliton and surface states for (a) U/γ1 = 0.45 and (b) U/γ1 = 0.55.The black dotted line is the flat band E = U , the red dashed line is the solution corresponding to the '+' sign in Eq. (24), and the blue solid line corresponds to the '−' sign.
For example, (2, 2, 1) would indicate onsite energies of +U, −U, +U, −U, −U, +U, −U, +U, +U, −U .The band structure and density of states of a solitonantisoliton pair (5, 6, 5) with N = 16 layers is shown in Fig. 5 for different values of U .There are two localized states associated with the soliton, two with the antisoliton, and two surface states.The two surface states, one soliton state and one antisoliton state generally lie within the bulk band gap.At k = 0, the soliton has energy E sol = −U + γ 1 , the antisoliton E asol = U − γ 1 , and the surface states have energies ±U .For U ≤ γ 1 /2, the soliton and antisoliton state hybridize with the two surface states, Fig. 5(b) and Fig. 5(h).For γ 1 > U > γ 1 /2, the soliton and antisoliton are closer to zero energy than the surface states, and they tend to hybridize together with a tiny anticrossing, Fig. 5(c) and Fig. 5(d).For U > γ 1 , the soliton and antisoliton states separate, leaving an overall band gap E g = 2(U − γ 1 ), Fig. 5(f) and Fig. 5(l).

FIG. 7 .
FIG. 7. Low-energy band structures for ABA stacking with N = 16 layers and a sharp soliton near the center.The top row shows band structures, the bottom row shows the corresponding density of states g(E), with the left side showing a soliton at the center (8, 8) after an even number of layers, and the right side shows a soliton near the center (7, 9) after an odd number of layers.There are different values of the magnitude of the alternating onsite energies as (a), (c), (e), (g) U/γ1 = 0.5, (b), (d), (f), (h) U/γ1 = 2.0, with different scales on the axes.Horizontal arrows in the top row show the energies at kx = 0 of the states localized on the soliton.The dashed line in (b) shows the perturbative approximation for the antisoliton energy E asol (k) and that in (d) shows the soliton energy E sol (k) where E asol (k) = −E sol (k) and E sol (k) is given in Eq. (31).In all plots, parameter values are γ0 = 3.16 eV, γ1 = 0.381 eV[40], a = 2.46 Å[41].For the band structures, ky = 0, and, for the density of states, δ = 0.01γ1.

FIG. 9 .
FIG.9.Low-energy band structures for ABA ′ stacking.Each column shows a different value of the onsite energies U .The top row shows a faultless system with N = 16 layers, whereas the second row shows N = 16 layers with a sharp fault at its center(8,8).In all plots, parameter values are γ0 = 3.16 eV, γ1 = 0.381 eV[40], a = 2.46 Å[41].

FIG. 10 .
FIG.10.Low-energy band structures for AA stacking.Each column shows a different value of the onsite energies U .The top row shows a faultless system with N = 16 layers, whereas the second row shows N = 16 layers with a sharp fault at its center(8,8).In all plots, parameter values are γ0 = 3.16 eV, γ1 = 0.381 eV[40], a = 2.46 Å[41].

FIG. 14 .
FIG.14.Band structures and density of states g(E) for single layer defects near the center (7, 1, 8) of a system with N = 16 layers for ABA, ABC, and AA ′ stacking.In the band structure plots, horizontal arrows show the energies at kx = 0 of the states localized on the defect.There is a different axis scale for AA ′ stacking, and the other stackings at large values of U/γ1.In all plots, parameter values are γ0 = 3.16 eV, γ1 = 0.381 eV[40], a = 2.46 Å[41].For the band structures, ky = 0, and, for the density of states, δ = 0.01γ1.

FIG. 17 .
FIG. 17. Rhombohedral stacking with N = 10 layers and onsite energy U = 0.1 eV showing band structures (left) and electrical conductivity (right).(a), (b) show the faultless system and (c), (d) show a system with an atomicallysharp soliton at the center using the minimal model with γ2 = γ3 = γ4 = 0. (e), (f) show the faultless system and (g), (h) show a system with an atomically-sharp soliton at the center using the full parameter model.Parameter values are γ0 = 3.16 eV, γ1 = 0.381 eV, γ2 = −0.02eV, γ3 = 0.315 eV, γ4 = 0.044 eV, a = 2.46 Å.The conductivity calculation is performed for a system of length L = 1754a and width W = 1754a, where a is the lattice constant.Details are described in Appendix A.
(a) and (b) are for the minimal tight-binding model (with γ 0 , γ 1 , and U only, Eq. (1)) and the faultless system.The region of the band gap, Fig.17 (a), corresponds to zero conductivity, Fig.17(b), and, when the energy is equal to that of higher-energy bands, the presence of additional transport modes produces steps in the conductivity.Fig.17 (c) and (d) are for the minimal tight-binding model with an atomically-sharp soliton at the center.The state localized on the soliton creates a non-zero density of states resulting in non-zero conductivity, and the position of the strong hybridization of the soliton state with the surface states corresponds to a characteristic dip in the conductivity (at E ≈ U = 0.1 eV in Fig. 17 (d)).
(e) and (f) are for the faultless system, whereas Fig.17 ( of 2 × 2 blocks.Intralayer blocks are Di = (a).For the soliton, there are two states localized on the soliton within the bulk band gap, Fig.18(b), and this is qualitatively similar to the minimal model, Fig.11(e).However, h-BN is in the wide band gap regime with U/γ 1 ≈ 7.25, and, as a result, the soliton states do not extend very far into the band gap, in relative terms.