Exploring dielectric properties in atomistic models of amorphous boron nitride

Calculations of the dielectric constant of the small (i.e. typically 100 atoms in the unit cell) structures were done with the Vienna ab initio Simulation Package (VASP) and projector-augmented wave pseudopotentials [13–17] using density functional perturbation theory [18]. The exchange-correlation functional was calculated in the generalized gradient approximation [19], and the 1s electrons were treated as core states for both B and N atoms. We only consider the electronic contribution to the dielectric constant and only calculate the static value. For all calculations, we use a cutoff energy of 520 eV in the plane-wave expansion and converge the energy to 10⁻⁶ eV. A k-point grid density of 100 k-points per Å⁻³ of the reciprocal cell was used for almost all calculations in this dataset, with only 1 data point with 80 k-points per Å⁻³ of the reciprocal cell. With these calculation parameters, the band gap of hexagonal BN is 4.48 eV, while the dielectric constants are $\left(\epsilon_{xx} = \epsilon_{yy}, \epsilon_{zz}\right) = (4.58, 2.15)$.

Using atomistic simulation techniques, the generation of the amorphous structures needed as input for the electronic calculations can be challenging. This is due to their complexity, the diverse atomic environments they exhibit, and the need for large unit cells to accurately account for their disordered nature. In this work, to reduce the computational cost and accelerate the generation of the aBN structures (and allow the generation of the large structures used in section 3), we use molecular dynamics with two machine-learned interatomic potentials trained on DFT data. This allows us to efficiently generate amorphous aBN samples with first-principles accuracy. The first of these interatomic potentials is a Gaussian approximation potential trained by Kaya et al [11] which was used for all structures with B/N ratios not equal to 1, as well as some of the structures with B/N = 1. Since the potential of Kaya et al [11] was trained using first-principles data generated with Quantum Espresso, for consistency, the final structures were fully relaxed with VASP (until the atomic forces were ${\unicode{x2A7D}} 10^{-3}\,\mathrm{eV}$ Å⁻¹) before the dielectric constants were calculated. To avoid the need for structure optimization using VASP for later structures, a second interatomic potential was trained [20] on first principles data generated from VASP using DeePMD [21–25] and DP-GEN [26]. The training data for the DeePMD potential was generated using the same pseudopotentials and exchange-correlation functional approximation as used for the dielectric constant calculations (detailed above). The model was trained on 181 621 structures with 64, 96, and 100 atoms, and validated on 20 181 structures. The mean absolute error of the energy per atom predicted by the interatomic potential on the validation set was 27.3 meV/atom, and the mean absolute error of the force components was 0.355 eV Å⁻¹. Structures generated with this potential all had B/N = 1 and did not require additional relaxation before the dielectric constant was calculated. Further discussion regarding the pertinence of machine-learned interatomic potentials in this context and the construction of the DeePMD potential may be found in appendix A.

The same procedure was followed for the generation of the small structures regardless of the potential used. First, B and N atoms were randomly scattered using PACKMOL [27] in a box with volume chosen to give a target initial density for the structure at hand. Densities were centered around 2.1 g cm⁻³ (the density of hexagonal BN and approximately that of experimentally grown low-k a-BN [6, 7]), 3.46 g cm⁻³ (the density of cubic BN), and 2.8 g cm⁻³ (in between). To generate low-energy structures, simulated annealing was done using LAMMPS [28] on the initial configurations of randomly scattered atoms. Each structure is heated from 300 K to 3000 K over 200 ps of simulation time, held at 3000 K for 20 ps, and then cooled back to 300 K at either 100 K ps⁻¹ or 13.5 K ps⁻¹ before a final relaxation. Some structures were generated in the NPT ensemble, and some were generated in the NVT ensemble to keep the density fixed. When the NPT ensemble was used, the lattice in the final relaxation was also allowed to relax, and when it was not, the lattice was held fixed.

Using the 209 aBN samples in the dataset, we try to reveal the relationship between their microstructure and their static dielectric constant ε₁. We first investigate the effect of B/N imbalance, density, and hybridization of atoms (${\mathrm{sp^1, sp^2}}$, and ${\mathrm{sp^3}}$) on the dielectric constant and electronic band gap. Then, we use the Cowley short-range order parameter (SRO) [29] as a way to evaluate the effects of disorder in boron–nitrogen alternation.

Let us first examine the question of B–N imbalance. Figure 1(a) displays the average dielectric constant of dataset samples with a given difference in their Boron and Nitrogen concentrations $\Delta c = \frac{n_\mathrm{B} - n_\mathrm{N}}{n_\mathrm{B} + n_\mathrm{N}}$ ⁹ . In agreement with Lin and co-authors [7], we find that B–N imbalance leads to a significant increase in the dielectric constant, both in the B-rich and N-rich cases.

Zoom In Zoom Out Reset image size

Figure 1(b) displays the observed relationship between density, dielectric constant and band gap. We recover the usual trends that a lower density and a higher band gap favor a lower dielectric constant (as can be intuited from equation (4) below, ($\epsilon_1-1) \propto \Omega^{-1}$ where Ω is the cell volume, and the presence of transition energies E_c  - E_v in the denominators), but we attract the reader's attention to the strong variability in ε₁ even at fixed density. This variability, in turn, is partially explained by the variation of the gap (color on figure 1(b)), which is an electronic property related to the precise details of the system's microstructure.

To provide insight on these electronic properties, we display in figure 2 the ensemble averaged density of states (DOS) over all equal concentration ($\Delta c = 0$) low-density samples in the dataset, i.e. samples with a density around that of hBN or lower, as well as one representative sample of this set ¹⁰ . As can be seen on figures 1(b) and 2, typical gaps for individual samples are found to vary between 0 and ${\sim} 3 \ \textrm{eV}$, significantly lower than the DFT gap of hBN. The DOS of the representative sample visibly displays mid-gap states, which in these small samples effectively set the value of the electronic gap. Given the data of figure 1(b), this suggests that these mid-gap states may have a notable influence on the system's dielectric constant. The exact energy of these states depends on the specific structure of the sample, so that in a larger sample displaying many different local atomic configurations, one may expect to see a distribution of these states inside what would otherwise be the gap of a 'reference crystalline structure', in a manner akin to the tail states typically observed in other amorphous solids [30]. To test this idea, we have performed an ensemble average of the DOS in the subset of the database described above. From this, it can indeed be seen that these mid-gap states effectively fill the gap. In disordered systems, even in this case, such states often do not actually contribute to DC conduction due to localization effects [30]. However, likely since the samples of the dataset are small (a few hundreds of atoms), we could not convincingly demonstrate localization in this particular case. We will discuss this in section 3, where we consider larger samples using a simple tight-binding model. We still mention here that, even if these states are localized, they can still contribute to the dielectric constant of the system by inducing low-transitions to or from delocalized states (see section 3 for details).

Zoom In Zoom Out Reset image size

To get a better understanding of the samples' microstructures, we examine their radial distribution functions (RDF). Figure 3 depicts an ensemble average for the RDF of all samples in the dataset. There is a very clear first peak, allowing us to define a first nearest neighbor shell through a cutoff radius of $r_\mathrm{c} = 1.9$ Å.

Zoom In Zoom Out Reset image size

Using this definition for the first nearest neighbor shell, we display in figure 1(c) the correlations between the coordination number fractions of the samples and their dielectric constants. We define the coordination number of a given atom by counting its neighbors within the cutoff radius $r_\mathrm{c}$, with the idea that $2/3/4$-coordinated atoms typically correspond to $\mathrm{sp}^1/\mathrm{sp}^2/\mathrm{sp}^3$ environments. The coordination number fractions f_i are then the ratio of the number n_i of i-coordinated atoms in a sample over its total number of atoms $n_\mathrm{tot}$: $f_i = n_i/n_\mathrm{tot}$. In this picture, structures with $f_3\to1$ tend to be ordered $\mathrm{sp}^2$ dominated structures, close to 'hBN order', while structures with $f_4\to1$ tend to be close to an ordered 3D allotrope of BN. Let us first focus on the second panel of figure 1(c), depicting the fraction f₃ of sp² coordinated atoms. Samples with the lowest dielectric constant ε₁ occur for f₃ = 1, i.e. fully sp²-coordinated samples. It can be seen that, as f₃ decreases away from 1, ε₁ tends to increase steadily: in effect, non-sp² coordinated atoms in an sp²-rich structure appear to act as defects which degrade the dielectric properties. The third panel of figure 1(c) shows the same behavior for sp³-rich structures, although the latter are seen to have a higher dielectric constant compared to the sp²-rich case. This last point can likely be understood from the fact that the sp³-dominated structures tend to be markedly denser, as expected structurally. In fact, since the samples in the dataset mostly have sp² and sp³ coordinated atoms (with a small fraction of sp¹), $f_3+f_4\approx 1$ and so the second and third panels are effectively mirrors of each other. As can be seen in the first panel, since no structures are sp¹-dominated, sp¹-coordinated atoms tend to always act as defects detrimental to ε₁. We should note, however, that such sp¹ bonds are likely very reactive, and that in real samples we may observe different effects due to contaminants which are not modeled here.

To further this discussion on the local bonding character in the samples, we now take the view of aBN as a disordered binary (A-B) alloy and examine the boron–nitrogen alternation in the samples. To this end, we rely on the Cowley short range order (SRO) [29]. The SRO for the i^th nearest neighbors shell of a given atom j of type B is defined as follows:

$$\begin{equation} {\mathrm{SRO}}_{B_j}\left(i\right) = 1-\frac{1}{c_A}\frac{n_i^{\left(j\right)}\left(A\right)}{N_i^{\left(j\right)}}\end{equation} \tag{ 1 }$$

where $n_i(A)$ is the number of atoms of type A in atom j's i^th shell, $N_i^{(j)}$ the total number of atoms in that same shell, and c_A is the global concentration of A in the sample. The ratio $n_i^{(j)}(A)/N_i^{(j)}$ can be seen as the local concentration of A atoms in the i^th shell: if the attribution of atomic types in the alloy were random, then on average one would have $n_i^{(j)}(A)/N_i^{(j)} = c_A$ and therefore the SRO would vanish for all shells. Oppositely, if the alloy was perfectly alternating in the sense that consecutive shells alternate in composition as purely one type of atoms (i.e. if j is a B atom, its first shell contains only A atoms, its second shell only B atoms, etc), then the SRO would oscillate between $1-\frac{1}{c_A}$ (−1 for equal concentrations alloys) and 1 for odd and even shells respectively. The average SRO of order i relative to a given atomic type (say B) for a sample is then obtained by averaging over all atoms of that type:

$$\begin{equation} {\mathrm{SRO}}_{B}\left(i\right) = \sum_j {\mathrm{SRO}}_{B_j}\left(i\right).\end{equation} \tag{ 2 }$$

It remains to be seen how nearest neighbor shells are to be attributed in aBN, given the variability in the samples. To this end, we follow (and use) the Pyscal [31] implementation of the Cowley SRO. Using the averaged RDF displayed in figure 3, we define a cutoff for the 2nd nearest neighbor shell, $r_\mathrm{c}^{(2)} = 3.1$ Å. Then, for each atom j, the average distance between its closest and furthest neighbors within $r_\mathrm{c}^{(2)}$ is used to define the boundary between the first and second shells for the purpose of the SRO calculation.

Figure 4 shows the average SRO relative to boron atoms for the first and second shells for all structures in the dataset. Overall, we observe that structures with the lowest dielectric constants tend to have the lowest first-shell SRO, which, as per the previous discussion, corresponds to the more ordered structures in terms of alternation. A similar trend can be observed from the second shell SRO, for which the highest values correspond to the more ordered alternation. It has been suggested by several authors [7, 32] that boron clusters or boron–boron bonds could lead to a high dielectric constant. Glavin and co-authors [32] have specifically noted the possibility that B–B bonds could create mid-gap states and conductive pathways in their samples. This is consistent with our analysis thus far, and echoes with the idea that non-alternating samples in our dataset also tend to exhibit higher dielectric constants.

Zoom In Zoom Out Reset image size

To summarize, we have so far observed that departure from an sp²- or sp³-rich phase tended to produce samples with comparatively higher dielectric constants compared to the 'pure' structure (with sp²-rich phases producing the lowest dielectric constants in general), as did departing from a B–N alternating structure in the Cowley SRO sense (with equal B and N concentrations and alternating samples producing the lowest dielectric constants in general). Both of these indicators can be seen as measures of disorder in the samples. It thus follows that, with the important caveat that only small unit cells were considered in the dataset, more disordered aBN samples are expected to display higher dielectric constants. This is in agreement with the findings of Lin and coauthors [7], who found that slower growth rates led to films with lower dielectric constants, since slower growth rates tend to produce samples with less disorder.

However, taking this line of reasoning to the extreme may suggest that the systems with the lowest dielectric constants would be the most ordered ones, i.e. the crystals. This seems paradoxical, as many works [6, 7, 33] have shown that aBN can have a lower dielectric constant than its crystalline counterparts. For this reason, we turn to the investigation of large samples.

We employ a melt-quenching protocol to generate samples, as presented in figure 5, using the Large-scale Atomic/Molecular Massively Parallel Simulator and machine-learning based GAP model generated by Kaya et al [11]. Here, we first place an equal number of B and N atoms in a simulation box (1) and melt the samples at 5000 K (2). Later, the samples are cooled down to 2000 K fast, after a short equilibration run, and then quenched to 1000 K with a high cooling rate (3). The samples are later cooled down to 300 K by employing different cooling rates (4). Finally, we employ an annealing step at 500 K to reduce the amount of $\mathrm{sp^1}$-hybridized atoms and homonuclear bonds (5).

Zoom In Zoom Out Reset image size

In this section, we will focus our analysis on two different samples obtained with different cooling rates (sample aBN2 is quenched faster than sample aBN1) and annealing rates (the properties of sample aBN1 are explored for different annealing rates, from aBN1-0 having the slowest annealing to aBN1-3 the fastest). Detailed values of their cooling rates are given in table 1.

Finally, we visualize the samples as in figure 5 using the code OVITO [34]. Samples aBN1−x, which correspond to the same initial structure at different times in the annealing process, have similar morphological properties. While sample aBN1-0 (longest annealing time) has only B–N bonds and mostly $\mathrm{sp^2}$-hybridized atoms, other aBN1−x samples have a low amount of N–N bonds (${\lt} 5\%$). The amount of $\mathrm{sp^3}$-hybridized atoms are 6.0% for each aBN1−x sample and the ratio of $\mathrm{sp^1}$ and $\mathrm{sp^2}$-hybridized atoms are between 17.1%–20.2% and 73.8%–76.9%, respectively. While the ratio of $\mathrm{sp^1}$ is 20.2% for sample aBN1-3, it is reduced slightly with a higher cooling rate during annealing. Samples aBN1−x and aBN2 all have a density of $1.967 \pm 0.01$. The composition of sample aBN2 is 7.2% $\mathrm{sp^1}$, 84.4% $\mathrm{sp^2}$ and 8.4% $\mathrm{sp^3}$, and nearly all bonds are B–N bonds (${\gt} 95\%$). Even though this sample has the same density as the aBN1−x samples, there are several noticeably large microvoids within the structure.

The dielectric response of the system is described by its dielectric function, $\epsilon(\mathbf{q},\omega)$. In the single particle approximation and within the long wavelength limit ($\mathbf{q}\to\mathbf{0}$), the electronic contribution to the latter is given by [35]:

$$\begin{equation} \epsilon\left(\mathbf{0},\omega\right) = 1+\frac{8\pi}{\Omega}\frac{\hbar^2}{{m_e}^2}\sum_{\alpha,\beta} \frac{\vert{\langle{\beta}\vert{\vec{e}\cdot\hat{\mathbf{p}}}\vert{\alpha}\rangle}\vert^2}{\left(E_\beta-E_\alpha\right)^2}\frac{f\left(E_\alpha\right)-f\left(E_\beta\right)}{E_\beta-E_\alpha-\hbar\omega -i\eta}\end{equation} \tag{ 3 }$$

where Ω is the system's volume, m_e is the electron mass, $\vec{e} = \frac{\mathbf{q}}{\vert{\mathbf{q}}\vert}$ is the electric field's polarization direction, $\hat{\mathbf{p}}$ is the momentum operator, the $(E_\alpha, \vert{\alpha}\rangle)$ and $(E_\beta, \vert{\beta}\rangle)$ are the system's electronic eigenpairs and $\eta\to0^+$ acts as a phenomenological broadening. Since we remain in the long wavelength limit in this work, we will omit the q-dependence the dielectric function below and write simply $\epsilon(\omega)$ for $\epsilon(\mathbf{0},\omega)$.

It is fruitful to consider separately the real and imaginary part of $\epsilon(\omega)$, denoted here respectively $\epsilon_1(\omega)$ and $\epsilon_2(\omega)$. $\epsilon_1(\omega)$ describes the system's dielectric screening, and the low-frequency dielectric constant, which is the figure of merit in this work, is given by $\epsilon_1(\omega = 0)$. At zero temperature, and in the $\eta\to0^+$ limit, we have:

$$\begin{equation} \epsilon_1\left(\omega\right) = 1+\frac{16\pi}{\Omega}\frac{\hbar^2}{{m_e}^2}\sum_{v,c} \frac{\vert{\langle{c}\vert{\vec{e}\cdot\hat{\mathbf{p}}}\vert{v}\rangle}\vert^2}{\left(E_c-E_v\right)^2}\frac{E_c-E_v}{\left(E_c-E_v\right)^2-\left(\hbar\omega\right)^2}\end{equation} \tag{ 4 }$$

where $v,c$ refer respectively to the valence (filled) and conduction (unfilled) states in the system. In the same conditions, the imaginary part of the dielectric function is given by:

$$\begin{equation} \epsilon_2\left(\omega\right) = \frac{8\pi}{\Omega}\frac{\hbar^2}{{m_e}^2}\sum_{v,c} \frac{\vert{\langle{c}\vert{\vec{e}\cdot\hat{\mathbf{p}}}\vert{v}\rangle}\vert^2}{\left(E_c-E_v\right)^2}\delta\left(E_c-E_v-\hbar\omega\right)\end{equation} \tag{ 5 }$$

$\epsilon_2(\omega)$ is typically related to the system's absorption spectrum. While the methods discussed in this work are not necessarily suitable to precisely access optical properties, $\epsilon_2(\omega)$ nevertheless provides insight on the contributions to $\epsilon_1(\omega = 0)$. Indeed, it can be seen from equation (5) that it provides an energy-resolved description of the optical matrix elements $\frac{\vert{\langle{c}\vert{\vec{e}\cdot\hat{\mathbf{p}}}\vert{v}\rangle}\vert^2}{\left(E_c-E_v\right)^2}$ that make up $\epsilon_1(\omega = 0)$. This can also be seen from the Kramers–Krönig relation [35]:

$$\begin{equation} \epsilon_1\left(\omega\right) = 1+\frac{2}{\pi}\int_0^{+\infty}\frac{\omega^{\prime}}{{\omega^{\prime}}^2-\omega^2} \epsilon_2\left(\omega^{\prime}\right) \mathrm{d}\omega^{\prime}\end{equation} \tag{ 6 }$$

To access the electronic properties of aBN in a qualitative way, we rely on a simple tight-binding model. We employ a first nearest neighbor Slater–Koster model [36] fitted on the cubic, wurtzite and single-layer hexagonal allotropes of BN at equilibrium and under 10% isotropic dilation to determine the Slater–Koster parameters and their dependence on interatomic distances. B–B and N–N hoppings are parametrized ad hoc, using the equivalent B–N parameters as order of magnitude estimates. Details of the fitting procedure and parameters are discussed in appendix B. We must stress that, while the tight-binding methodology accounts correctly for the different local coordination numbers and geometries, this model does not include any energetic corrections for them apart from the aforementioned distance-dependent hoppings. It is therefore limited only to a qualitative exploration of the phenomena at play and will serve to guide intuition. In fact, we found that a minimal one-orbital toy model displayed similar trends. We briefly discuss this other model in appendix C.

Figure 6 presents the DOS of two $11\,520$ atoms structures generated using the GAP methodology (see section 3.1) with different quenching rates ¹¹ . Mid-gap states are present in both cases, although in manifestly lesser amounts than in the DFT results of section 2. This discrepancy may be due either to the better relaxation of these large structures, in the sense that they are less constrained by small-cell periodic boundary conditions, to the lesser statistical impact of single 'defects' in a comparatively large structure with a broad range of local environments, or by the inability of the tight-binding model to fully describe the energetic correction associated to them. Even with these limitations, we do observe markedly more mid-gap states in the sample with the faster quenching rate (aBN2): in fact, the electronic gap effectively vanishes in this sample, as suggested by the averaging procedure for small samples performed in 2. In contrast, the slowly quenched sample (aBN1) retains an electronic gap of about 4–5 eV, reduced from its value of $6.1\,\textrm{eV}$ for the reference crystalline hBN in our simple model (see hBN PSL in section 3.4). A similar effect occurs for samples aBN1-x, for which we observe a more marked population of mid-gap states for shorter annealing times, although these samples always retain an electronic gap.

Zoom In Zoom Out Reset image size

In section 2, we had suggested that the mid-gap states may be localized in real space, but the small size of the samples made it difficult to explore that hypothesis. Having now moved to larger samples, we investigate this point by computing the inverse participation ratio (IPR) for all eigenstates as an indicator of localization. For a given eigenstate Ψ, we compute its IPR as:

$$\begin{equation} \mathrm{IPR}\left(\Psi\right) = \sum_{\mathbf{n}} \left(\sum_\mu\vert{\langle{\mu,\mathbf{n}}\vert{\Psi}\rangle}\vert^2\right)^2\end{equation} \tag{ 7 }$$

where the first sum is taken over atomic positions and the second over the corresponding orbitals, with $\vert{\mu,\mathbf{n}}\rangle$ being the orbital $\mu\in\left\{B_{2s}, B_{2p_x},B_{2p_y},B_{2p_z}\right\}$ or $\left\{N_{2s},N_{2p_x}, N_{2p_y},N_{2p_z}\right\}$ at position n. Heuristically, if a state is delocalized over the whole system, i.e. if it has a Bloch-wave type behavior, $\vert{\langle{\mu,\mathbf{n}}\vert{\Psi}\rangle}\vert^2 \approx \frac{1}{N}$ where N is the number of orbitals in the system, and therefore $\mathrm{IPR}(\Psi) \approx N\times\frac{1}{N^2}\to 0$. However, for localized states, which in the extreme case have probability densities of the form $\vert{\langle{\mu,\mathbf{n}}\vert{\Psi}\rangle}\vert^2 \approx \delta_{n,n_0}$ for some position n₀, the IPR is finite (and of order unity).

The results are reported in figure 6. It can be seen in particular that, while IPR within 'bands' are almost zero, suggesting delocalized states, the IPR for states at the band-edges, and, notably, of the mid-gap states are of order unity. This suggests that these states are indeed localized in real space, and consequently should not or only weakly contribute to DC conduction in the system [30]. However, their impact on the system's dielectric constant is less clear. On the one hand, as can be recovered from equation (4), low-energy transitions such as those between mid-gap states can provide a strong contribution to ε₁, due to the presence of $(E_c-E_v)^3$ in the denominators. On the other hand, transitions between localized states are typically difficult because the optical matrix elements between them, the $\vert{\langle{c}\vert{\vec{e}\cdot\hat{\mathbf{p}}}\vert{v}\rangle}\vert^2$, are essentially local in real space and can be thus suppressed for the same reasons that the DC conductivity is suppressed (if the valence and localization states in a given (v, c) pair are localized at different positions in the solid, their overlap is small). In fact, it has been shown that, close to the Anderson transition, the dielectric constant of 3D systems should approximately correlate with the square of the wave function localization length [37]. Even in this case, however, (optical) transitions between localized mid-gap states and extended 'band-like' states are a priori still possible, and because ε₁ integrates transitions from all energies, the presence of mid-gap states, even localized, can increase it.

To investigate these questions more directly, we now turn to the computation of $\epsilon(\omega)$. To this end, we use a common tight-binding approximation to express the momentum matrix elements in terms of the Hamiltonian matrix elements [38–41]:

$$\begin{equation} \langle{\mu,\mathbf{n}}\vert{\hat{\mathbf{p}}}\vert{\mu^{\prime},\mathbf{n^{\prime}}}\rangle = -\frac{m_e}{i\hbar}\left(\mathbf{n^{\prime}}-\mathbf{n}\right)\langle{\mu,\mathbf{n}}\vert{\hat{H}}\vert{\mu^{\prime},\mathbf{n^{\prime}}}\rangle\end{equation} \tag{ 8 }$$

Equation (8) also substantiates our earlier remark about the local character of the momentum matrix elements, which in this tight-binding formulation can be seen to be a direct consequence of the local character of the tight-binding Hamiltonian (itself a consequence of the exponential decay of the atomic orbitals).

Figure 7 displays the dielectric function for samples aBN1−x and aBN2 (the DOS of aBN1−0 and aBN2 are also shown in figure 6) ¹² . As expected from the tendency of large systems to self-average, we find that at these scales (${\sim} 10^4$ atoms) the dielectric properties of the system are essentially isotropic ¹³ . We therefore represent the average of ε over the cartesian directions. To provide a basis for comparison, we also display the dielectric function for bulk hBN averaged in such a manner, using the same model. Because the hoppings in our tight-binding model have a cutoff radius that is inferior to the interlayer distance in hBN, interlayer hoppings are not taken into account for consistency. We thus call this reference system 'Periodically-stacked single layer hBN' (hBN PSL) for clarity. It also results from this that the dielectric constant of hBN PSL reduces to 1 along the stacking direction ¹⁴ . In effect, this 'stacking' thus only enters the computation of ε through the interlayer distance, which affects the system's volume Ω (see equation (3)).

Zoom In Zoom Out Reset image size

Panel (a) displays $\epsilon_1(\omega)$. The figure of merit, here, is the static dielectric constant $\epsilon_1(\omega = 0)$. We must stress that in the context of aBN, the low-frequency dielectric constant in capacitance experiments is typically measured at frequencies of the order of hundreds of kHZ [6, 7], which correspond to energies $\hbar\omega$ of the order of $10^{-9} \ \textrm{eV}$. Variations of $\epsilon_1(\omega)$ on these energy scales likely have origins in vibrational (or other non-directly electronic) contributions to ε₁, which are not included in our bulk frozen-atoms calculations. With this caveat, it can be seen that, in this model, aBN does have a lower average dielectric constant than our reference crystalline system. We also note that ε₁ decreases as the sample's annealing or quenching time is increased, as can be seen in the inset (a'), although the effect in this formulation is relatively weak.

To understand these effects we turn to the imaginary part of the dielectric function, which is displayed in panel (b). As per equation (5), $\epsilon_2(\omega)$ resolves in energy both the density of transitions and their strength (it is in effect the joint DOS weighted by oscillator strengths). Optically active transitions can be observed below the transition edge of the crystalline reference system, which by definition can be ascribed to the contribution of mid-gap states. We also observe the suppression of the van-Hove singularities, which is symptomatic of the loss of crystal order in the samples. In addition, and in contrast to what happens below the crystalline gap, there is a sizeable reduction of the values of $\epsilon_2(\omega)$ above the reference crystal gap, which, through equation (5), can be interpreted as an overall weakening of oscillator strengths in the system. Because the overall oscillator strength-density of mid-gap states is low in these models, even if they have lower transition energies, the net effect through equation (4) or (6) is a reduction of the static dielectric constant compared to the reference crystalline sample. This is further shown through the whole energy range of transitions by a Kramers–Krönig analysis in appendix E.

The above results suggest that there is a competition between two effects. On one side, as was discussed in section 2.3, disorder in the samples has a tendency to increase the dielectric constant, likely through the emergence of mid-gap states. This effect is also observed in this section, although only weakly in our simple tight-binding model. On the other side, as seen above, the amorphous structure of the sample yields an overall reduction in the oscillator strength of its transitions compared to our crystalline reference, which results in a lowering of its dielectric constant. In addition, disorder induces localization effects which are visible in large samples, and should contribute to a lowering of the optical matrix elements. To obtain a clearer understanding of these combined effects, it thus seems necessary to possess a good description of the mid-gap states, and more generally of the energetics of the numerous atomic configurations and local environments present in the samples (onsites, hoppings$\ldots$). This description is a priori available to ab initio techniques, but it is lacking in our simple tight-binding model. Conversely, the large-scale effects of disorder that we just discussed are very computationally demanding for ab initio methods. It thus appears desirable to move towards more sophisticated tight-binding (or generally second principles) models. Finally, it is believed that amorphous BN materials also contain hydrogen or oxygen atoms in non-negligible quantity, impacting their structural stability and mechanical properties [12], so that one should also consider proper tight-binding modeling of these additional atomic species for accessing the full dielectric response of measured samples.

For completeness, we make here a remark and a further caveat. Owing to the complexity and size of the amorphous systems considered here, we have remained at the theoretical level of (effectively) independent particles. It is well known, however that in hBN, quasiparticle corrections and electron–hole interaction effects (excitons) are far from negligible [42, 43]. The effect of such many-body corrections is typically larger if the Coulomb potential is poorly screened, making aBN a likely candidate to display such effects despite its 3D nature. While even in such systems, DFT methods have been used to compute quantities such as the static dielectric constant [44], they are typically insufficient to describe their optical properties. The same limitation applies to single-particle tight-binding methods such as the ones used here. Nevertheless, as per section 3.2, we believe that the ω > 0 behavior of the dielectric function, even at this level of theory, can be helpful both as a first step in this direction, and in understanding the static trends. Obtaining an accurate estimation of $\epsilon(\omega \gt 0)$ which could be compared with the results of optical experiments is however likely to necessitate some form of the aforementioned many-body corrections. This is however out of the scope of this work.

To conclude this section, we point out that significant work remains to be done to reach accurate modeling of aBN systems. In this context, it is not so surprising that the values of the static dielectric constants computed in this work remain above two, while experimental values of $\epsilon_1 \approx 2$ and below have been reported in the litterature [6, 7].

As underlined in the main text, generating atomistic structures of amorphous materials requires a model that can simultaneously handle many different atomic environments and a very large unit cell. While DFT can handle diverse atomic environments, it is limited to unit cells with a small number of atoms. On the other hand, Molecular Dynamics (MD) simulations, although they do not directly provide electronic information, can handle very large numbers of atoms thanks to the use of robust and quick empirical interatomic potentials. However, these potentials are often derived from experimental or high-level computational data of crystalline phases, in which amorphous phases may not be well represented. Therefore, empirical interatomic potentials cannot accurately capture the local structural variations and complex nature of amorphous structures. Machine learning-based interatomic potentials, such as Gaussian approximation potential (GAP) [55] and DeePMD [21], offer a bridge between the realms of very accurate DFT simulations and very robust MD simulations. By (machine) learning the local energies and forces associated to local atomic configurations from ab initio calculations on small structures without relying on predefined rules, these models can capture the atomic interactions in diverse local atomic environments with close to DFT accuracy, but at a fraction of the computational cost. This alleviates the cost of subsequently generating large sets of small structures, and makes it possible to generate large structures, both of which can later be used as input for electronic calculations.

In this work, we generated all samples used in DFT and tight-binding electronic calculations using an existing GAP model [11] and a new DeePMD model trained over large DFT-generated sets of structures to ensure the accuracy of MD simulations. The details of the training set and the training protocol for the new DeePMD model are presented below.

Several different descriptors are available in DeePMD-kit [21, 22, 25]. In this work, we use the DeepPot-SE descriptor with radial and angular information of the atomic configuration, referred to as $\mathrm{se\_e2\_a}$ in the code. We set a cutoff radius of 6 Å for the local atomic environment and set a maximum of 200 B and 200 N atoms within the cutoff radius, above the number of neighbors in any of our training configurations. We set the sizes of the embedding and fitting neural nets to (25, 50, 100) and (120, 120, 120), respectively.

Training data for the DeePMD potential was generated using the DP-GEN active learning framework [26]. DP-GEN utilizes an ensemble of DeePMD models trained on the same dataset but with different initial parameters to select samples for DFT calculations and addition to the training set as follows: first, the potential energy surface of the system is explored by running molecular dynamics simulations (with user defined parameters) on selected systems using a DeePMD interatomic potential (boot strapped on an initial set of training data), configurations are generated. Then, using the ensemble of models, the variance in the atomic forces of the configurations is used to calculate an error indicator ε, which is a measure of the deviation between models. Finally, configurations with model deviation within a range (very small deviations are configurations for which the models are all accurate, and very large deviations are likely unphysical) are deemed as good additions to the training set and selected for DFT calculation. After a new ensemble of models is trained with the newly added data, the process repeats.

To generate the training data for the DeePMD aBN interatomic potential used in this work, we used the following settings in the DP-GEN active learning process. The exploration steps involved molecular dynamics runs at temperatures ranging from 500 K to 3500 K in both the NVT and NPT ensembles. In the NPT ensemble, pressures ranging from 0 to 360 000 Bar were used. Several dozen structures were used as the initial structures for the molecular dynamics runs in exploration steps. The structures included supercells of crystalline hexagonal-, cubic-, and wurtzite-BN, as well as 100 atom amorphous structures generated using simulated annealing (initially done with ab initio molecular dynamics and, when the interatomic potential had trained for some time, preliminary versions of the interatomic potential). The range of model deviation in the forces we used to select structures for the training set is 0.05 eV Å⁻¹ $\unicode{x2A7D} \epsilon \lt 0.35$ eV Å⁻¹. 181 621 structures were generated for the training set in this manner. An additional 20 181 structures were used for validation

Once trained, the DeePMD interatomic potential has mean absolute error of 27.3 meV/atom in the total energy across systems in the validation set. The average mean absolute error of the force components is 0.355 eV Å⁻¹. The interatomic potential is capable of $\gt$1 ns of simulation time per hour of wall time, thereby enabling very fast structure generation for the dataset used in this work.