Exploring dielectric properties in atomistic models of amorphous boron nitride

We report a theoretical study of dielectric properties of models of amorphous Boron Nitride, using interatomic potentials generated by machine learning. We first perform first-principles simulations on small (about 100 atoms in the periodic cell) sample sizes to explore the emergence of mid-gap states and its correlation with structural features. Next, by using a simplified tight-binding electronic model, we analyse the dielectric functions for complex three dimensional models (containing about 10.000 atoms) embedding varying concentrations of sp1, sp2 and sp3 bonds between B and N atoms. Within the limits of these methodologies, the resulting value of the zero-frequency dielectric constant is shown to be influenced by the population density of such mid-gap states and their localization characteristics. We observe nontrivial correlations between the structure-induced electronic fluctuations and the resulting dielectric constant values. Our findings are however just a first step in the quest of accessing fully accurate dielectric properties of as-grown amorphous BN of relevance for interconnect technologies and beyond.


Introduction
As the need for larger information storage and processing explodes, improvements of computing device performances such as the lateral scaling, and integration of active devices in the back-end-of-line are crucially needed [1].In particular, the optimization of interconnect technologies becomes increasingly important as the signal delay of the interconnect increases rapidly compared to gate delay of the transistor [2,3].It is therefore critically important to develop new materials or improve existing ones to decrease interconnect energy loss, through dielectric constant and metal resistivity reduction (resistance-capacitance (RC) delay).Since the early 2000s, SiCOH, with a dielectric constant of about 3.0, has been the material of choice instead of SiO 2 , which has a dielectric constant of 4.0.Then, by embedding pore structures with a dielectric constant of 1.0, the dielectric constant of porous-SiCOH (p-SiCOH) could be linearly reduced, although unfortunately this decreases its Young's modulus exponentially, making it difficult to perform post deposition processes such as chemical mechanical processing and packaging ports.Ultimately, this effectively limits the dielectric constant of p-SiCOH in such applications to about 2.5 [4,5].
Recent reports of amorphous Boron Nitride (aBN) have shown unprecedented ultralow dielectric constant of 1.8, together with a robust Young's modulus of over 50 GPa in absence of pore morphology, sparking a great interest in view of the long sought-after improved interconnects technologies [6,7].Additionally, aBN has been found to display good mechanical properties overall as well as excellent thermal properties and diffusion barrier against metal migration.All such properties turn out to be perfectly suited for the development of next-generation interconnect technologies [5,8,9].
Within this context, we present a theoretical study of electronic and dielectric properties of complex models of aBN, and we attempt to correlate them with the system's atomic structure.Indeed, because aBN is an amorphous system, its structure, and therefore its properties of merit, depend strongly on its fabrication parameters (growth rate, temperature, B-N stoichiometry,. ..) [7,10].To progress in material optimization, proper and accurate simulation studies are thus highly desirable.
However, unfortunately, the modeling of highly disordered (large-scale) materials is generally a very challenging task for predictive theoretical investigation.Indeed, realistic modeling of nanoscale electronic properties in such complex systems benefits strongly from the predictive power of first-principles techniques, which can be quite accurate.But at the same time, such methods are also limited in accessible sample sizes (typically several thousands of atoms) and are generally unable to capture disorder effects such as multiple scattering, impurity states, quantum interference and localization in large-scale (realistic) models.
To cope with the complexity of the modeling problem, we employ a two-pronged approach.First, we generate a dataset of small aBN samples (lattices with ∼10 2 atoms in their periodic cell) through simulated annealing with machine learned interatomic potentials [11] and compute their electronic and dielectric properties using density functional theory (DFT).This allows us to explore the influence of the system's microstructure on its electronic and dielectric properties with ab initio precision.In this framework, we show the emergence of some correlations between the nature of the structural disorder and bonding character statistics inside the sample, and the resulting formation of mid-gap states and behavior of the dielectric constant (section 2).Second, we construct large scale atomic structural models (∼10 4 atoms) using molecular dynamics with such machine learning potentials, which enjoy ab-initio accuracy [11,12], and evaluate the electronic and dielectric properties of these large systems through a simple Slater-Koster tight-binding model.Although this latter model presents significant limitations, it allows us to get a first picture of the electronic and dielectric properties in large, highly complex three-dimensional geometries mixing boron and nitrogen atoms, linked together by sp 1 , sp 2 and sp 3 bonds.The obtained results are not aimed at quantitative predictions, but we expect a certain validity concerning some qualitative trends.
Finally, by summarizing the limitations of employed methodologies, we put in perspective the obtained results and discuss the need for more sophisticated modeling strategies to achieve better future quantitative predictions, able to guide experiments and research at the industrial level.

Generation of a dataset of small structures
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][14][15][16][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 −6 eV.A k-point grid density of 100 k-points per Å −3 of the reciprocal cell was used for almost all calculations in this dataset, with only 1 data point with 80 k-points per Å −3 of the reciprocal cell.With these calculation parameters, the band gap of hexagonal BN is 4.48 eV, while the dielectric constants are 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 ⩽10 −3 eV Å −1 ) 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][22][23][24][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 Å −1 .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 −3 (the density of hexagonal BN and approximately that of experimentally grown low-k a-BN [6,7]), 3.46 g cm −3 (the density of cubic BN), and 2.8 g cm −3 (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 −1 or 13.5 K ps −1 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.

Dataset exploration
Using the 209 aBN samples in the dataset, we try to reveal the relationship between their microstructure and their static dielectric constant ϵ 1 .We first investigate the effect of B/N imbalance, density, and hybridization of atoms (sp 1 , sp 2 , and 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 ∆c = nB−nN nB+nN 9 .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.
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, (ϵ 1 − 1) ∝ Ω −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 ϵ 1 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 (∆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 10 .As can be seen on figures 1(b) and 2, typical gaps for individual samples are found to vary between 0 and ∼3 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 9 We prefer this metric of B-N imbalance over the B/N ratio because it is symmetric in B and N. 10 We choose this set of samples because they are of comparable density with that of hBN.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).
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 c = 1.9 Å.
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 c , with the idea that 2/3/4-coordinated atoms typically correspond to sp 1 /sp 2 /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 tot : f i = n i /n tot .In this picture, structures with f 3 → 1 tend to be ordered sp 2 dominated structures, close to 'hBN order' , while structures with f 4 → 1 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 3 of sp 2 coordinated atoms.Samples with the lowest dielectric constant ϵ 1 occur for f 3 = 1, i.e. fully sp 2 -coordinated samples.It can be seen that, as f 3 decreases away from 1, ϵ 1 tends to increase steadily: in effect, non-sp 2 coordinated atoms in an sp 2 -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 3 -rich structures, although the latter are seen to have a higher dielectric constant compared to the sp 2 -rich case.This last point can likely be understood from the fact that the sp 3 -dominated structures tend to be markedly denser, as expected structurally.In fact, since the samples in the dataset  Radial distribution function for small aBN samples averaged over the whole dataset.A marked peak corresponding to the first nearest neighbor shell can be observed, which allows the definition of a cutoff radius for first-nearest neighbor interaction at rc = 1.9 Å.A weaker second peak is also visible, which allows for the estimation of a second nearest neighbor shell cutoff r mostly have sp 2 and sp 3 coordinated atoms (with a small fraction of sp 1 ), f 3 + f 4 ≈ 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 1 -dominated, sp 1 -coordinated atoms tend to always act as defects detrimental to ϵ 1 .We should note, however, that such sp 1 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: where n i (A) is the number of atoms of type A in atom j's i th shell, N i 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 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 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 − 1 cA (−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: 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 (2) c = 3.1 Å.Then, for each atom j, the average distance between its closest and furthest neighbors within r (2) c 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.

Partial conclusions
To summarize, we have so far observed that departure from an sp 2 -or sp 3 -rich phase tended to produce samples with comparatively higher dielectric constants compared to the 'pure' structure (with sp 2 -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.

Molecular dynamics and tight-binding investigation of large size models
We now move to the study of larger systems, using a combination of molecular dynamics to generate large structural atomic models, and a simple tight-binding model to investigate the electronic and dielectric response of these structures.

Generation of large structural atomic models
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 sp 1 -hybridized atoms and homonuclear bonds (5).
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 sp 2 -hybridized atoms, other aBN1−x samples have a low amount of N-N bonds (<5%).The amount of sp 3 -hybridized atoms are 6.0% for each aBN1−x sample and the ratio of sp 1 and sp 2 -hybridized atoms are between 17.1%-20.2%and 73.8%-76.9%,respectively.While the ratio of 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 ± 0.01.The composition of sample aBN2 is 7.2% sp 1 , 84.4% sp 2 and 8.4% sp 3 , and nearly all bonds are B-N bonds (>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 function
The dielectric response of the system is described by its dielectric function, ϵ(q, ω).In the single particle approximation and within the long wavelength limit (q → 0), the electronic contribution to the latter is given by [35]: where Ω is the system's volume, m e is the electron mass,⃗ e = q |q| is the electric field's polarization direction, p is the momentum operator, the (E α , |α⟩) and (E β , |β⟩) are the system's electronic eigenpairs and η → 0 + 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 ϵ(ω) for ϵ(0, ω).
It is fruitful to consider separately the real and imaginary part of ϵ(ω), denoted here respectively ϵ 1 (ω) and ϵ 2 (ω).ϵ 1 (ω) describes the system's dielectric screening, and the low-frequency dielectric constant, which is the figure of merit in this work, is given by ϵ 1 (ω = 0).At zero temperature, and in the η → 0 + limit, we have: 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: ϵ 2 (ω) 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, ϵ 2 (ω) nevertheless provides insight on the contributions to ϵ 1 (ω = 0).Indeed, it can be seen from equation ( 5) that it provides an energy-resolved description of the optical matrix elements |⟨c|⃗ e•p|v⟩| 2 (Ec−Ev) 2 that make up ϵ 1 (ω = 0).This can also be seen from the Kramers-Krönig relation [35]:

A simple tight-binding model
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 11 .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 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.
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: ) 2 (7) where the first sum is taken over atomic positions and the second over the corresponding orbitals, with |µ, n⟩ being the orbital µ ∈ Heuristically, if a state is delocalized over the whole system, i.e. if it has a Bloch-wave type behavior, |⟨µ, n|Ψ⟩| 2 ≈ 1 N where N is the number of orbitals in the system, and therefore IPR(Ψ) ≈ N × 1 N 2 → 0. However, for localized states, which in the extreme case have probability densities of the form |⟨µ, n|Ψ⟩| 2 ≈ δ n,n0 for some position n 0 , 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 ϵ 1 , 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 |⟨c|⃗ e • p|v⟩| 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 ϵ 1 integrates transitions from all energies, the presence of mid-gap states, even localized, can increase it.6).The light green curve (hBN PSL) corresponds to the response of periodically stacked hBN single layer (see text).Inset (a') displays a zoom on the low-frequency part of ϵ1, showing how a longer annealing or a shorter quenching time may lead to a lower dielectric constant.All panels share their legend.The higher energy behavior of ϵ2(ω) and a more detailed analysis of transitions are given in E.

Dielectric response
To investigate these questions more directly, we now turn to the computation of ϵ(ω).To this end, we use a common tight-binding approximation to express the momentum matrix elements in terms of the Hamiltonian matrix elements [38][39][40][41]: 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) 12 .As expected from the tendency of large systems to self-average, we find that at these scales (∼10 4 atoms) the dielectric properties of the system are essentially isotropic 13 .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 14 .In effect, this 'stacking' thus only enters the computation of ϵ through the interlayer distance, which affects the system's volume Ω (see equation ( 3)).
Panel (a) displays ϵ 1 (ω).The figure of merit, here, is the static dielectric constant ϵ 1 (ω = 0).We must stress that in the context of aBN, the low-frequency dielectric constant in capacitance experiments is typically 12 See appendix D for computational details, and caveats regarding the computation of ϵ1(ω) for aBN2. 13It should be noted that our sample generation procedure does not include the presence of a substrate.In real systems, and for thin enough films, substrate interaction could induce an anisotropic response, e.g. by favoring alignment of BN rings. 14As a consequence, the averaged dielectric function which we display for hBN PSL for consistency is given by ε = 1 + 2 3 (ϵ ∥ − 1), where ϵ ∥ is its dielectric constant for in-plane polarized fields.measured at frequencies of the order of hundreds of kHZ [6,7], which correspond to energies hω of the order of 10 −9 eV.Variations of ϵ 1 (ω) on these energy scales likely have origins in vibrational (or other non-directly electronic) contributions to ϵ 1 , 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 ϵ 1 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), ϵ 2 (ω) 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 ϵ 2 (ω) 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.

Beyond current modeling
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. ..).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 ϵ(ω > 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 ϵ 1 ≈ 2 and below have been reported in the litterature [6,7].

Conclusions and perspectives
We report on some exploratory modeling studies towards the understanding of electronic and dielectric properties in models of aBN.We combined the analysis of a dataset of small (∼10 2 atoms) samples at the DFT level with the study of large (∼10 4 atoms) samples generated with machine learning interatomic potentials whose electronic properties have been explored using a simple tight-binding model.In both cases, we have observed the formation of mid-gap states inside the otherwise large ionic gap (as observed in clean hBN for example).The density of these states, and the dielectric constant of the samples are tentatively related to the system's atomic density, and to the nature of its atomic bonding characteristics.For small samples, we highlight in particular that, in sp 2 dominated systems, the presence of sp 1 and sp 3 coordinated atoms appears to act as a source of extra disorder which increases the dielectric constant.Likewise, departure from a B-N alternating structure or a B/N = 1 stoichiometric ratio are also observed to have detrimental effects on the system's dielectric constant.However, the disordered nature of these materials demands to consider very large scale effects, which are out of scope for first-principles simulations.We have thus analyzed large systems through the prism of a simple tight-binding Hamiltonian, and, in this model, found them to display an overall weakening of their transitions's oscillator strengths due to their amorphous nature, which in turn lowers their dielectric constant.We also observe real space localization effects for the mid-gap states, which may lead to a partial suppression of their contribution in low-energy transitions and therefore a lower contribution to dielectric properties.To obtain a more realistic and complete view of the competing phenomena at play, efforts are therefore further needed to develop more refined tight-binding Hamiltonians, retaining the physics related to the bonding effects in the various observed structural conformations.It will demand adequate adjustment of parameters with ab-initio data on smaller system sizes, eventually harnessing the power of machine learning techniques to automatize the extraction of tight-binding parameters in an arbitrary complex disordered structure [45,46].Such models seem particularly appropriate in the context of aBN, as tight-binding naturally accounts for the geometric structure of the system, while benefiting from well-established linear scaling methods [47][48][49], which allow the study of extremely large samples (>10 6 atoms).The use of Artificial Intelligence-based techniques, meanwhile, has become key to obtain realistic electronic models able to deliver quantitative predictions [50,51].This will be fundamental to reach a predictive modeling of dielectric constants or other properties in amorphous forms of boron-nitride composites, enabling further materials optimization and performance improvement [52][53][54] 15 .
Center (TACC) at The University of Texas at Austin.O K is supported by the REDI Program, a project that has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 101034328.S R acknowledges support from PID2019-106684GB-I00 funded by MICIU/AEI /10.13039/501100011033 and by ERDF and EU and 2021 SGR 00997, funded by Generalitat de Catalunya.

Appendix A. Structure generation using molecular dynamics with machine-learned potentials
A.1.Pertinence of machine-learning models for the generation of aBN structures 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.

A.2. Details of the DeePMD training procedure
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 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 Å −1 ⩽ ϵ < 0.35 eV Å −1 .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 Å −1 .The interatomic potential is capable of >1 ns of simulation time per hour of wall time, thereby enabling very fast structure generation for the dataset used in this work.

Appendix B. Slater-Koster tight-binding model
In this study, a simple tight-binding model is used to gain qualitative insights into the electronic and dielectric response of aBN.The model relies on atomic orbitals from s and p subshells (i.e. 4 orbitals per atomic site) to expand the low energy electronic manifold.An effective onsite energy h I θ , with θ ∈ {s, p} and I ∈ {B, N}, is associated with each atomic subshell so as to parameterize the diagonal Hamiltonian matrix elements.The non-diagonal Hamiltonian matrix elements are described as two-center integrals and parameterized according to the Slater-Koster formalism, where t γ ∈ {t ssσ , t spσ , t psσ , t ppσ , t ppπ } are the Slater-Koster parameters associated with electron hoping between boron and nitrogen atoms, r is the atomic distance (r 0 = 1.57Å), and β γ ∈ {β ssσ , β spσ , β psσ , β ppσ , β ppπ } are dimensionless parameters introducing functional dependence of energy integrals with respect to radial distance.The cutoff radius (r c = 1.7 Å) restricts effectively the direct interaction to nearest-neighbour atoms.Note that t spσ ̸ = t psσ as the left hand side is associated with hoping between B s − N p orbitals, while the right hand side relates to B p − N s hoping.The model parameters have been fitted with respect to reference electronic structure calculations by minimizing the following loss function, built from the difference between the reference and tight-binding eigenvalues (ϵ tb mk − ϵ ref mk ) across the first Brillouin zone (k ∈ BZ).Here, m is the band index and the sum runs over all bands.The weight function w(ϵ) is introduced to restrict the energy range.The weight function used here is depicted in figure B1.The integral is evaluated for different atomic configurations, here labelled Γ, that contribute equally to the loss function.The set of atomic configurations considered in this work is composed of various allotropes of Boron Nitride, namely the cubic, wurtzite and monolayer hexagonal allotropes, both at equilibrium volume and under an isotropic dilation (∆ = 10%).The reference data have been computed within the framework of DFT.The equilibrium structural parameters were obtained upon geometry optimization with VASP with the same functional and cutoff energy as in section 2 and regular k-points grids of 6 × 6 × 6, 6 × 6 × 6 and 16 × 16 × 1, respectively for the cubic, wurtzite and monolayer hexagonal allotropes.The set of reference eigenvalues used for the fitting were then obtained by self-consistently expanding the electronic manifold of the equilibrium and dilated allotropes upon a basis sets of numerical atomic orbitals (double-ζ + polarization) as implemented in the Siesta simulation package [56].Integration in reciprocal space was performed on regular grids characterized by an effective cutoff ⩾10 Å.The plane wave cutoff for the real-space grid was set to 800 Ry.The parameters obtained after fitting are reported in table B1.The band-structures computed with our tight-binding model are depicted in figure B1.Comparison with the reference data shows a surprisingly good accordance for the low energy valence eigenstates indicating a high degree of transferability for this manifold.The agreement is less pronounced for the conduction eigenstates.While some characteristic energy dispersion of the low energy conduction bands are well reproduced, other features associated with longer ranged interaction or with significant contribution from atomic subshells of higher angular momenta are missing.Nevertheless, our tight-binding model, despite its simple form, captures qualitatively well the energy dispersion range of the first few valence and conduction bands.
In the allotropes of BN considered above, only B-N bonds are present.This is not the case in aBN, where B-B and N-N bonds are also possible and have been observed both in simulation and experimentally [57]  16 .While we observe that, in large samples at least, such bonds are typically not dominant in number (<5% of 16 To our knowledge, while B-B bonds are observed experimentally, N-N bonds typically are not, although they are seen in simulations.nearest neighbor bonds for the samples of section 3.1), we still require a prescription to incorporate them in our tight-binding framework.Since our fit to the crystalline allotropes does not provide this information, we employ an ad hoc parametrization with the same functional forms.For hoppings between same subshell (i.e.s-s and p-p) we use the same hopping and distance parameters for the B-B and N-N bonds as for the corresponding B-N bonds.For the s-p hopping terms, we use the arithmetic average of the hopping parameters s-p and p-s derived for BN with distance parameters set to the Harrison scaling value of 2 [58].

Appendix C. One-orbital tight-binding toy model
In this appendix, we briefly discuss a very simple tight-binding toy model which nevertheless captures many of the features discussed in section 3. We consider, generally, a binary alloy with two atoms of very different electronegativities, with one effective isotropic orbital per atom.For concreteness, let the two species under consideration be B and N, so that we work with a basis of localized orbitals {|B, n⟩, |N, m⟩} n,m where n and m run over the Nitrogen and Boron atomic positions respectively.The effective Hamiltonian is then given by: In other words: Here h B and h N are fixed onsite energies for all atoms of the given species whose difference h B − h N > 0 accounts for the difference in electronegativity between B and N, t is the hopping element for a reference interatomic distance which corresponds to a typical equilibrium bond length r 0 , and g(|n ′ − n|) is a decreasing function which describes the evolution of the hopping strength with distance.In keeping with the simplicity of the model, we choose here an exponential decay with a first nearest neighbors cutoff r c : For consistency, we choose values for the parameters which correspond to the π orbitals of the Slater-Koster model used in the main text (see appendix B), viz.: h B = h B p = 6.1 eV, h N = h N p = 0 eV, t = t ppπ = −1.8eV and β = β ppπ = 3.0 with an equilibrium bond length of r 0 = 1.57Å and a cutoff radius r c = 1.9 Å.The crucial parameter here is the ratio between the typical hoppings and the difference in onsite energies, | t h B −h N |.We have selected here specifically the parameters of p orbitals, as they are typically the ones involved near the band gap for BN materials [59], and the pp π hopping parameters because such bonds make up the π/π * bands of hBN and we mostly consider sp 2 -dominated samples.
Figure C1 presents the DOS and IPR in this simple model for sample aBN1−0 and aBN2 (see figure 6 displays for their DOS in the full Slater-Koster model).Except for the model, we used the same calculation parameters as in the main text (see appendix D).Again, localized mid-gap states are observed, although they are here completely suppressed for the slower cooled sample.We also notice that the valence and conduction DOS are essentially symmetric, except for small difference in the mid-gap states.This is at variance with the situation in the full Slater-Koster model, as seen in figure 6.In particular, we notice there a markedly higher amount of mid-gap states on the conduction side as opposed to the valence side, while this is not the case in the one-orbital model.We hypothesize that this may be due to the relatively flat conduction bands produced by the Slater-Koster model (see figure B1), whose states can thus be localized more easily under the action of disorder.
In this simple model, the presence of mid-gap states appears strongly connected to the presence of B-B and N-N bonds.In many samples, we observed a strong diminution or even a suppression of mid-gap states denominators of equations ( 4) and ( 6)) and avoid the associated increase in ϵ 1 , as seems here to happen for hω ≳ 1 − 2 eV.
Finally, as a consistency check and an additional visualization of the above, we display in the bottom panel of figure E1 the quantity 1 + 2 π ´hω 0 ϵ2(ω ′ ) hω ′ dhω ′ (note the replacement of +∞ by hω in the integral bound), which is thus the contribution to the static dielectric constant of transitions up to a cutoff energy of hω.We verify in particular that ϵ 2 (ω) was indeed computed for energies hω high enough that the static value of ϵ 1 is recovered over the chosen range.

Figure 1 .
Figure 1.Distribution of dielectric constant, ϵ1, of small aBN samples with various B/N compositions (a), densities (b) and different hybridizations of atoms (c).In panel (a), average dielectric constants were estimated by averaging over all dataset samples with a given B-N composition, while error bars represent the associated standard deviation.Sixteen out of the 209 samples were found to exhibit high dielectric constants ϵ1 > 10 and are excluded from the above plots for lisibility.Out of these sixteen samples, seven have ϵ1 > 20 and up to a few hundreds, and these outliers are excluded from the averages in panel (a).

Figure 2 .
Figure 2. Left panel: averaged DOS of low-density, equal B and N concentration aBN samples (densities comparable to hBN-48 samples).Right panel: zoom on the region of the electronic gap for these samples.The red curve depicts one of the samples averaged over, for reference, while the DOS of other samples of that set are represented as translucent curves.Individual DOS are computed from the Kohn-Sham energies with a Lorentzian broadening of 26 meV.

Figure 3 .
Figure 3.Radial distribution function for small aBN samples averaged over the whole dataset.A marked peak corresponding to the first nearest neighbor shell can be observed, which allows the definition of a cutoff radius for first-nearest neighbor interaction at rc = 1.9 Å.A weaker second peak is also visible, which allows for the estimation of a second nearest neighbor shell cutoff r

Figure 4 .
Figure 4. Short range order of aBN samples based on first (a) and second nearest neighbor shells (b).

Figure 5 .
Figure 5. Melt-quenching protocol of aBN samples used GAP-MD simulations (temperatures and times are not to scale).Insets: aBN sample with 500 atoms where sp 1 , sp 2 , and sp 3 -hybridized atoms are shown as blue, green, and red, respectively.Bonding of B and N atoms in an aBN sample.

Table 1 .
Cooling rates for the MD aBN samples studied in this section.See text and figure5for the meaning of the different rates.

Figure 6 .
Figure 6.Tight-binding density of states (DOS, in black) and Inverse Participation Ratio (IPR, blue and red dots for valence and conduction states respectively) for two ∼ 10 4 atoms aBN samples with equal concentration of Boron and Nitrogen.States are present within the electronic gap (between about 0 and 5 eV in this simple model), but their high IPR suggests that they are localized in real space.Panels (a) and (b) show respectively samples aBN1−0 and aBN2 (aBN2 is quenched faster than aBN1); note how a faster quenching step increases the amount of mid-gap states.Insets: zoom on the DOS of the corresponding samples in the mid-gap region.The dashed line displays the energy of the last occupied state.

Figure 7 .
Figure 7. Real (panel (a)) and imaginary (panel (b)) parts of the dielectric function for the aBN samples generated with different quenching and annealing time (see text and DOS displayed in figure6).The light green curve (hBN PSL) corresponds to the response of periodically stacked hBN single layer (see text).Inset (a') displays a zoom on the low-frequency part of ϵ1, showing how a longer annealing or a shorter quenching time may lead to a lower dielectric constant.All panels share their legend.The higher energy behavior of ϵ2(ω) and a more detailed analysis of transitions are given in E.

Figure B1 .
Figure B1.Band structure of allotropes of Boron Nitride used to fit the Slater-Koster parameters.Cubic, monolayer hexagonal and wurtzite allotropes (from left to right) have been considered both at equilibrium volume (top row) and under an isotropic dilation of ∆ = 10% (bottom row).The reference DFT and tight-binding calculations are depicted in blue and orange respectively.The weight function (w(ϵ)) used to restrict the energy range of the fit is depicted on the right.

ĤFigure C1 .
Figure C1.Single orbital tight-binding density of states (DOS, in black) and Inverse Participation Ratio (IPR, blue and red dots for valence and conduction states respectively) for samples aBN1−0 and aBN2 (aBN2 is cooled faster than aBN1−0) in the one-orbital tight-binding toy model.Mid-gap states are observed again, and are suppressed by slower cooling rates.The dashed line denotes the energy of the highest occupied state.Figure 6 presents the DOS of the same samples using the full Slater-Koster model.

Figure E1 .
Figure E1.Top panel: imaginary part ϵ2(ω) of the dielectric function for the samples introduced in section 3 over the model's full transition energy range.Middle panel: Kramers-Krönig contribution of each transition energy to the static dielectric constant (see equation (E.1)).Bottom panel: cumulative energy integration of the former.Dotted horizontal lines: static dielectric constants of the samples (colors are matched with the corresponding solid lines).The legend is shared by all panels.