Modelling structure and properties of amorphous silicon boron nitride ceramics

Silicon boron nitride is the parent compound of a new class of high-temperature stable amorphous ceramics constituted of silicon, boron, nitrogen, and carbon, featuring a set of properties that is without precedent, and represents a prototypical random network based on chemical bonds of predominantly covalent character. In contrast to many other amorphous materials of technological interest, a-Si3B3N7 is not produced via glass formation, i.e. by quenching from a melt, the reason being that the binary components, BN and Si3N4, melt incongruently under standard conditions. Neither has it been possible to employ sintering of μm-size powders consisting of binary nitrides BN and Si3N4. Instead, one employs the so-called sol-gel route starting from single component precursors such as TADB ((SiCl3)NH(BCl2)). In order to determine the atomic structure of this material, it has proven necessary to simulate the actual synthesis route. Many of the exciting properties of these ceramics are closely connected to the details of their amorphous structure. To clarify this structure, it is necessary to employ not only experimental probes on many length scales (X-ray, neutronand electron scattering; complex NMR experiments; IRand Raman scattering), but also theoretical approaches. These address the actual synthesis route to a-Si3B3N7, the structural properties, the elastic and vibrational properties, aging and coarsening behaviour, thermal conductivity and the metastable phase diagram both for a-Si3B3N7 and possible silicon boron nitride phases with compositions different from Si3N4 : BN = 1 : 3. Here, we present a short comprehensive overview over the insights gained using molecular dynamics and Monte Carlo simulations to explore the energy landscape of a-Si3B3N7, model the actual synthesis route and compute static and transport properties of a-Si3B3N7.


I. Introduction
One of the most fascinating new classes of hightechnology materials are the amorphous nitridic ceramics of the composition a-Si x B y N l C z [1−11] such as a-Si 3 B 3 N 7 and a-SiBN 3 C.These compounds are synthesized via the sol-gel route and exhibit a very high stability against crystallization up to 1900 K and 2100 K for a-Si 3 B 3 N 7 and a-SiBN 3 C, respectively.Furthermore, e.g.a-SiBN 3 C is stable against oxidation up to 1700 K, and exhibits a high bulk modulus B = 200−300 GPa, while at the same time the density of the materials is very low, e.g.ρ ≈ 1.9 g/cm 3 for a-Si 3 B 3 N 7 , compared to the weighted average of the binary endphases Si 3 N 4 and BN, ρ = 2.8 g/cm 3 .
Understanding the reason for this surprising stability is of great importance for the future design of such new materials.Since these materials are amorphous, experimental data [11−13] cannot fully reveal their structure, and thus it is necessary to combine theoretical and experimental studies.We have focussed our theoreti-cal work on the prototypical material, a-Si 3 B 3 N 7 , where we have studied the actual synthesis route to a-Si 3 B 3 N 7 [14,15], the structural properties [16,17], the elastic and vibrational properties [17], aging and coarsening behaviour [18,19], thermal conductivity [20], and the metastable phase diagram both for a-Si 3 B 3 N 7 [21] and possible crystalline silicon boron nitride phases with compositions different from Si 3 N 4 : BN = 1 : 3 [22].Here, we present a short overview over the insights gained using molecular dynamics and Monte Carlo simulations to explore the energy landscape of a-Si 3 B 3 N 7 , model the actual synthesis route and compute static and transport properties of a-Si 3 B 3 N 7 .

II. Methods
For our simulations, we have employed a number of numerical techniques, ranging from straightforward molecular dynamics and Monte Carlo simulations to multistep modelling, including statistical and analytical models.We employed a classical potential taken from the literature [23] and used system sizes ranging from 351 to 22464 atoms; however, most of the work was performed with simulation cells of about 20×20×20 Å 3 containing N atom = 702 atoms (162 Si,162 B,378 N) or with cells containing N atom = 1300 atoms.We found that once we went beyond the shorter time scales, the results of the molecular dynamics and Monte Carlo simulations were quite compatible.Thus, we have used the computation of the diffusion coefficients with these two methods to establish a simulation time scale for the Monte Carlo simulations where 1 Monte Carlo cycle (i.e.N atom individual Monte Carlo steps) corresponds to about 1/2 fs.The step size of the moves in the Monte Carlo simulations was set such that for a given temperature about 50% of the moves were accepted, and the time step in the molecular dynamics calculations was usually 1/2 fs.We performed both constant volume and constant pressure simulations depending on the issue being addressed.For more details on particular calculations, we refer to the literature cited in the results section.

Modelling the sol-gel synthesis of a-Si 3 B 3 N 7
As mentioned above, and as we will note in the next subsection, the only way to generate realistic structure models of a-Si 3 B 3 N 7 is to simulate the actual synthesis route, i.e. the polymer-precursor or sol-gel route.Since the complete polymer-precursor route extends over many time and length scales, we have developed a stepping stone approach in order to simulate the synthesis as faithfully as possible [14,15].Initially, precursor molecules are in solution with an excess of NH 3 , which react with each other based on the number of available reactive sites per molecule and their density within the solution.During this first linking stage, we take the like-lihood that a given reaction attempt is successful to depend on whether a Si-N or a B-N bond is formed, with B-N bonds being energetically preferred.The diffusion rates of the individual molecules being relatively high, local depletion effects in the precursor-concentration do not play a big role and thus the precursor density can be taken to be spatially homogeneous.The first reaction stage can therefore be modelled by generating lists of linked precursor molecules (n TADB < 10) according to these probabilities.
However, once several precursor molecules have linked up to form a larger aggregate, this oligomer will move much more slowly compared to the remaining original precursor and NH 3 molecules.Thus, the oligomers will become essentially stationary and serve as condensation centres for the still mobile reactants.This allows us to model this latter phase of the syn thesis as a multiple condensation process of individual oligomers and monomers.At the end of this stage, we are left with many isolated oligomers, which are now beginning to cross-link to form the polymer stage.This is modelled by placing the oligomers randomly on a lattice and shrinking the average distance until they can interact and form bonds. Finally, the pyrolysis stage is simulated by a Monte Carlo simulation at T = 1200 K, which stays well below the computational melting temperature (T melt (Si 3 B 3 N 7 cryst ) ≈ 2500 K and T melt/glass (Si 3 B 3 N 7 amorph ) ≈ 2000 K) of the system but nevertheless eliminates most of the many dangling bonds still present in the polymer.During the pyrolysis simulations, the density of the polymer increases from ρ ≈ 1.5 g/cm 3 to the final value of ρ E = 1.8−2.0g/cm 3 .Figure 1 shows the resulting amorphous compound.

Structural modelling of a-Si 3 B 3 N 7 and computation of bulk properties
In principle, there exist a large number of ways to generate an amorphous ternary compound: quenching from the melt, amorphization under high pressure, etc.In order to gain additional understanding of a-Si 3 B 3 N 7 , we have generated five different classes of models representing five such synthesis routes [17]: (A) quenching from a melt via Monte-Carlo / molecular dynamics simulations; (B) sintering of BN and Si 3 N 4 nanocrystals [16], with a diameter of each fragment ≈ 1/2 nm; (C) sintering of amorphous Si 3 B 3 N 7 -clusters, where each cluster contains several hundred atoms and is generated via a locally optimized "build-up-build-down" algorithm [24,25]; (D) deposition from the vapour phase, modelled via a random close packing algorithm for the N 3− anions where the B 3+ and Si 4+ cations are distributed over the holes in the packing, followed by a tempering simulation [26]; (E) the actual sol-gel route described in the previous subsection.
The most surprising observation is that the radial distribution functions (RDF) of the random network struc-tures generated via these different routes are distinctly different: The variance for the many samples generated for one class of models is smaller than the average distance in "RDF-space" between models belonging to different classes.The best agreement with the experimental RDF is found for the models of class (B) and class (E).As far as the first local coordination of Si and B by N and N by Si and B is concerned, all models give very similar results, with SiN 4 -tetrahedra and planar BN-triangles being the most dominant coordinations by far, in agreement with experiment.However, the models differ regarding the second coordination sphere, i.e. the first Si-Si/B and B-Si/B coordination.Here, the models belonging to classes (A), (C) and (D) show a close-to-random statistical distribution of Si/B neighbours to a given B or Si atom.In contrast, the class (B) models generated by sintering binary nanocrystallites exhibit a very high preference for Si-Si and B-B neighbours over a mixed Si-B coordination.Best agreement with NMR-measurements, however, is found for the models of class (E), where a distinct but not overwhelming preference for Si-Si and B-B neighbours is observed (Fig. 1).
Regarding the overall density of the structures generated, models of classes (A), (C) and (D) show densities in the range of 2.7−2.8g/cm 3 , and those belonging to classes (B) and (E) have densities of 2.2 g/cm 3 and 1.8−2.0g/cm 3 , respectively.For comparison, the density of a hypothetical crystalline compound [27] Si 3 B 3 N 7 is about 2.9 g/cm 3 .Again, the structures generated by following the experimental synthesis route show the best agreement with experiment (ρ exp ≈ 1.8−1.9g/cm 3 ) [6].Thus, we conclude that the models produced by simulating the actual synthesis are most appropriate as far as structural properties are concerned.
When analyzing the models of class (E) in more detail, we note that they contain a large number of tiny voids with diameters ranging from 0.2 to 0.7 nm [19].This is again in good agreement with the experiment, since TEM-measurements did not indicate any voids larger than 2 nm [28].In contrast, e.g. the models generated by quenching from the melt show no voids with diameters larger than 0.3 nm [19].However, we note that as far as the total energy is concerned, the models with the lower overall density are higher in energy than those with high densities.Thus, one would expect that at sufficiently high temperatures a volume shrinkage of the ceramics might take place (unless the temperatures required are so high as to lead to the decomposition of the ceramics with loss of N 2 and/or the formation of Si 3 N 4 / BN crystallites).
In order to investigate this question, we have simulated the possible coarsening of the material generated by the sol-gel-route via very long constant pressure Monte Carlo simulations at several temperatures for systems with up to 5200 atoms [19].We find that for temperatures below the melting / glass transition temperature of the amorphous phase (2000−2250 K) the average size of the voids increases approximately like a power law in time while the total number of voids decreases in agreement with expectations from coarsening theory.The overall density increase can best be approximated by a logarithmic law.
Repeating these calculations at elevated pressures [19], the material behaves like an elastic solid (with a bulk modulus of about 150 GPa) at low temperatures up to about 5 GPa.But for pressures above this value, irreversible shrinkage takes place, and after releasing the pressure, the overall density has increased to about 2.8 g/cm 3 .Similarly, already relatively small pressures are sufficient to cause an irreversible shrinkage during simulations at temperatures of 1500 K and above.It would clearly be interesting to see, whether a highpressure-high-temperature experiment would produce such a high-density amorphous phase of a-Si 3 B 3 N 7 , which should have very good mechanical properties.We note that in real systems, after a small shrinkage in the volume, a crust would be formed at the surface that is essentially impenetrable to voids.As a consequence, while coarsening also takes place inside the real material, this happens at a much lower rate than the one observed in our constant pressure simulations, where the next "outside surface" is effectively only one simulation cell length (here: 2 nm) away!For all the models belonging to all five classes, we have computed their bulk moduli [17].We find that the bulk modulus computed from the E(V)-curve on average increases with the density of the model, with values ranging from B = 50 GPa to about B = 250 GPa, in good agreement with the result from the simulations above and experimental data (B exp ≈ 200 GPa) [29].
Similarly, the computed phonon spectra show the same shape for all the model classes [17,30] (Fig. 2): A large broad peak at about 0.05 eV (≈ 400 cm −1 ), from which the phonon density of states slowly decreases to a minimum at about 0.15 eV (≈ 1200 cm −1 ), followed by a smaller peak at around 0.17 eV (≈ 1400 cm −1 ), and no contributions beyond 0.25 eV (≈ 2000 cm −1 ).This second peak is less pronounced for the models of class (B) and (E).A comparison with calculated phonon densities of states for crystalline β-Si 3 N 4 and hexagonal BN shows that the high-frequency peak can be associated with vibrations of BN while the lower frequency bulge should stem mostly from Si-N lattice vibrations.
When we analyze the localization of the computed vibrational modes using the average participation ratio [31,32], ρ(X )(ω i ) = (N Σ i (e j (i) e j (i) ) 2 ) -1 (where e j (i) is the (3-dimensional) contribution of atom j to the normalized eigenvector e (i) ), we find that we can assign certain frequency ranges to particular atom and building unit contributions.The vibrational modes in the main peak at 0.05 eV are all non-localized, while those with energies greater than 0.1 eV are localized and only very few atoms (≈ 3−7 atoms) are involved.This is supported by repeating the calculations for a system of twice the size (2600 instead of 1300 atoms), where the participation ratio scales by a factor of 2 for the localized modes (as it should for such modes), but does not clearly scale for modes in the low-frequency range.If one considers the building unit projections of the eigenmodes, one finds that the vibrations centered on Si and B atoms are dominated by vibrations of SiN 4 -and BN 3 -units, respectively.In addition, one finds in the low-energy modes also contributions of under-coordinated SiN 3 -and BN 2 -units, and in the high-frequency region the peaks at about 0.2 eV are clearly due to BN 2 -building units.
Where the models belonging to different classes can be most clearly distinguished is in the amount of contributions of units of the type NSi 3−x B x (x = 0,1,2,3), due to the different probabilities to find such units in the different models.The main peak in the nitrogen contribution to the vibrational densities of states is dominated by either NSi 2 B 1 units (in particular for models from class A) or NSi 3 -units.However, in general one would not expect that one could distinguish between the model classes just based on the measured phonon densities of states.The above assignments are in qualitative agreement with the experimental IR-spectra [33], which show broad peaks at 0.11 eV (≈ 930 cm −1 ) and 0.166 eV (≈ 1340 cm −1 ); the frequencies above 0.2 eV (≈ 1620 cm −1 ) which are missing in the experiment, are due to defects in the local coordination of the simulated structures which would be expected to have been annealed out during the experimental synthesis.
Regarding transport properties, the ceramics are insulators electronically, but their thermal conductivity is clearly of interest.Thus, we have performed a computer experiment [20], where we select two planes parallel  to the y-z-plane of the orthorhombic simulation cell that are 1/2L apart, where L is the side length of the simulation cell in the x-direction.After equilibrating the system at the simulation temperature T (T = 20, 50, 100, ..., 500 K), we raise the temperature of the atoms in one plane by a small amount ΔT and lower the temperature at the other plane by the same amount, and keep the two planes at these temperatures for the remainder of the molecular dynamics simulation.After equilibration of the system in this configuration, a temperature gradient has been established between the two planes, and we measure the amount of heat that needs to be added to the hotter plane and removed from the colder plane, in order to keep the temperatures of the two planes constant.Taking the ratio of heat flow to thermal gradient yields the is the average heat flow through a cross section A, and | T | = 4ΔT/L is the temperature gradient, for the given length L of the simulation cell.Repeatedly doubling the length of the cell in the x-direction and again performing the simulation yields κ(L), and by plotting κ(L) versus 1/L, one can derive the thermal conductivity for an infinite system.We find κ = 4 W/m/K, essentially constant for the whole temperature range investigated (Fig. 3).Since κ has not yet been measured for a-Si 3 B 3 N 7 , we estimate the reliability of the procedure by studying the corresponding binary systems.We find that the experimental value exceeds the one from the simulations by a factor of 3−5, similar to what we have found for Selen [34].Thus, one would conclude that the true value of κ should be about κ ≈ 15 W/K/m.Similarly, we can compute the dissipation of a heat pulse in a-Si 3 B 3 N 7 , which allows us to measure the speed with which energy is transported through the system, υ S ≈ 3−4×10 3 m/sec, a value which lies in the range of the speed of sound found experimentally in crystalline Si 3 N 4 (υ S ≈ 1.1 × 10 4 m/sec) and BN (υ S ≈ 7 × 10 3 m/sec).For comparison, we have performed simulations of energy transport in Si 3 N 4 and BN, which yield velocities υ S ≈ 7.5 × 10 3 m/sec and υ S ≈ 8 × 10 3 m/sec, respectively.We also find the dissipation time for the heat pulse τ diss ≈ 50 fs; this time competes with the time for energy transport through the simulation cell, τ trans = L/ υ S ≈ 0.2 ps.3.3 Phase diagrams in the Si/B/N system: a-Si 3 B 3 N 7 and the quasi-binary Si 3 N 4 /BN system When we consider the energies of formation of the various (meta) stable compounds in the Si/B/N system, we find the following ranking by energy [17,21]: For the overall composition Si 3 B 3 N 7 , the lowest energy is found for the corresponding combination of the binary end phases, BN : Si 3 N 4 = 3 : 1, followed by a number of hypothetical crystalline compounds of this composition [27,35].The next higher in energy are the dense amorphous phases we have generated in models (A), (C) and (D), then come the sintered nanocrystallites (1/2 nm diameter), and the ones highest in energy are the amorphous structures found via the sol-gel route.We note that the larger the nanocrystallites selected for the annealing are, the lower the final energy of the sintered product will be: Once their size exceeds ≈ 3 nm, the energy lies below the dense amorphous phase, and above ≈ 15 nm diameter, the sintered crystallites are lower in energy than the best hypothetical ternary compounds.Of course, such sinters of large crystallites no longer form true ternary phases but instead constitute a nanocrystalline array of Si 3 N 4 and BN crystallites.However, this analysis explains, why all attempts have failed to generate amorphous silicon boron nitride via sintering of μm size crystallites [6].

Δ Δ
Nevertheless, there are two important questions that should be considered: how would the metastable phase diagram of Si 3 B 3 N 7 look like as function of temperature and pressure, and what is the standard phase diagram for the quasi-binary Si 3 N 4 /BN system, with a particular focus on the possible existence of thermodynamically stable ternary compounds?We have addressed the first issue [21] by performing long constant volume Monte Carlo simulations for a large number of volume-temperature combinations ranging from 250 K to 7000 K and densities ρ N from 0.12 atoms/Å 3 to 10 −5 atoms/Å 3 , where the pressure was evaluated using the standard virial formulation [36].
A central part of the analysis is to determine, whether for the given pressure/volume and temperature conditions the system is in the solid, liquid or gaseous state.Several criteria are used, which are based on one-(time)point properties such as potential energy and size distribution of connected atom clusters as well as on two-point properties such as the diffusion coefficient and the bond-survival probability (corresponding to the likelihood that the same two atoms are nearest neighbours at two substantially different points in time).We first assign the system to the gaseous or "condensed" (solid / liquid) state of matter according to whether the cluster size distribution favours small units containing only a few atoms or large clusters with hundreds of atoms, respectively.Next, the bond-survival probabilities (BSP) are employed to distinguish between the liquid and solid state, since the liquid state exhibits low values for the BSP. Figure 4 shows the resulting phase diagram.We find that the condensed amorphous phase stays solid for all temperatures up to 1750 K for all densities, beyond which softening and a glass transition occurs.Similarly, the condensed crystalline phase remains solid up to the melting temperature of 2500 K. Up to about 3000 K a liquid phase is present, which at higher temperatures behaves like a non-ideal gas with a liquid-gas-coexistence region and a critical point at T cr ≈ 8000 K, p cr ≈ 1.3 GPa and ρ cr ≈ 0.032 atoms/Å 3 , which is in good agreement with a fit of the p(V,T)-curves to a van-der-Waals gas model.We also note that we can estimate the N 2 pressure required to keep the ternary liquid from decomposing to be about 4 GPa at 2000 K.
Regarding the search for possible ternary compounds in the Si/B/N-system, we have employed our modular approach to structure prediction [27,37−39].This method is based on the search for local minima and locally ergodic regions on the energy landscape of the chemical system, which correspond to the (meta) stable compounds capable of existence [40].For compositions Si 3 N 4 : BN = 1 : 1, 2 : 3, 1 : 2, 1 : 3, and 1 : 4, we have performed global searches for local minima on the empirical potential energy landscape, with up to 6 formula units per simulation cell [22].As search algorithms, we use simulated annealing [41], the threshold algorithm [42] and the prescribed-path method [43], where both atom positions (≈ 80% of the moves) and cell parameters (≈ 20% of the moves) are freely varied.Subsequently, we rank all structure candidates using ab initio energy calculations on Hartree-Fock level.Figure 5 shows the phase diagram found at 0 K, where only the best structure candidate for each composition investigated is shown.We find that for all compositions the thermodynamically stable state is achieved by a phase separation into the two binary end phases BN and Si 3 N 4 .Among the possible metastable ternary compounds, the "best" candidates are found for the composition Si 3 B 3 N 7 .Since all these compounds constitute quite intricate nitridic networks with many possible side minima on the energy landscape favouring the formation of amorphous phases, the actual synthesis is expected to be rather difficult, though all these compounds should be highly stable once they have been synthesized.An overview over the most promising structure candidates is given in the Appendix (Tables 1 and 2). .G: gaseous phase; L: liquid phase; S: solid phase; L-G: gas-liquid co-existence region; S-L: gas-solid co-existence region; S-L: crystalline solid / super-cooled liquid.Note that each data point corresponds to an average of several simulations performed with a finite simulation cell at the given values of T and ρ.

Properties of the energy landscape of a-Si 3 B 3 N 7 : relaxation and aging phenomena
Underlying the properties of a-Si 3 B 3 N 7 is the highly complex energy landscape of the Si 3 B 3 N 7 system, with a very large part belonging to the amorphous phase, leading to many non-ergodic features in the dynamics of the system.We determine the temperature below which non-ergodicity prevails by investigating, whether the fluctuation-dissipation theorem (FDT) holds for long simulations at a given temperature.As indicator variable, we use the potential energy and check, whether the derivative of the average energy with respect to temperature equals the fluctuations in the energy during the simulation: where ‹ E › tω< t < tω+tobs (T ) is the average potential energy of the configurations generated at temperature T and averaged over time after a waiting time t ω and σ 2 (E) is the variance of the energy during the simulation run.
A crucial feature of non-ergodic behaviour is the dependence of many quantities on both the starting and the end point of the measurement, t ω and t ω + t obs , in contrast to ergodic situations where only the elapsed time t obs matters.We find that below 2000 K, the FDT is violated, which manifests itself in a freezing-out of the structure, i.e. a strong increase in the bond-survival probability obeying a Kohlrausch-Williams-Watts law, BSP(t) = exp(−(t/τ BSP (T)) β ) (β ≈ 0.3−0.5)[18], a peak in the specific heat around 2000−2500 K[18], and a rapid decrease in the diffusion constant D = A(T − T X ) γ , with T X ≈ 2000 K.In the liquid phase above 2500 K, the diffusion coefficient found, D ≈ 2 × 10 −4 cm 2 /s, is in good agreement with the one obtained for amorphous silica, while in the solid phase below 1750 K, D becomes very small.
One particularly fascinating aspect of non-ergodicity is the aging behaviour of the complex system, i.e. if one lets the system relax on a time scale of t w in a nonergodic state, then for subsequent observation times t obs small compared to t ω , the FDT appears to hold, while for times t obs > t ω the FDT is violated.This can be studied explicitly by investigating the two-time energy-energy av- where the subscript "ens" always denotes taking an average over all trajectories in the ensemble of simulations.One finds that for t obs < t ω , φ ≈ 1, while for t obs > t ω φ increases logarithmically with t obs , as one would expect from non-ergodic systems.A similar aging behaviour is observed when computing the bond-survival probability and the diffusion coefficient as function of t ω and t obs [18].
Finally, one can connect the non-equilibrium behaviour to the properties of the energy landscape of a-Si 3 B 3 N 7 .For this, we explore the landscape via many long constant temperature Monte Carlo simulations, where we periodically stop the simulation and perform many stochastic quenches to explore the minima lying  56 below the trajectory.A semi-log plot of the average energy of the minima x min observed vs. the time along the trajectory for temperatures below 2000 K shows that ‹ E(t;T,x min ) › ens , decreases logarithmically with time, similar to what is seen in other complex systems [44].
A fit of the logarithmic slope as function of temperature yields A(T) = 76T − 134T 2 , which qualitatively agrees with results from so-called LS-tree models [45], suggesting that the landscape of a-Si 3 B 3 N 7 might possess some hierarchical aspects in the low-energy range relevant for the dynamics at temperatures below the glass transition temperature.This agrees with the observation that if one lowers the temperature below 1500 K, the average energy of the minima encountered starts to increase again with decreasing simulation temperature, in contrast to what one would expect from equilibrium-like behaviour [18].On the other hand, for temperatures above 2000 K in the (super-cooled) liquid state, the average energy of the minima below the trajectory increases with temperature as expected in equilibrium situations.Furthermore, we note that the distance in configuration space between the holding point and the corresponding quench points is about an order of magnitude larger than the distances between the quench points belonging to the holding point.This indicates that the holding points are typically associated with only one large basin on the landscape containing many similar individual minima and saddle points, and not with a transition region connecting two such basins.
Next, the vibrational spectra of the holding points (Fig. 6) and the minima (Fig. 2) are computed, in or-der to determine, to what degree the walker hops between basins or moves smoothly above them [30].For T < 1000 K, the fraction of imaginary modes is close to zero, i.e. the trajectory of the system spends most of its time close to local minima, hopping once in a while between them.Between 1000 K and 2000 K, the harmonic approximation breaks down and the number of imaginary modes increases, with a maximum in the increase between 2000 K and 2500 K when the system moves from the glassy to the liquid state.At the same time, the character of the imaginary modes changes from localized below 2000 K to extended motion for T > 2000 K. Here, Si and N atoms contribute the most to these modes indicating that Si-N bond breaking is the dominant mechanism for the structural re-arrangements in the amorphous phase.However, even at 7000 K, only 20% of the modes at the holding points are imaginary, in nice agreement with our observation during the phase diagram simulations that relatively long-lived clusters exist in the liquid phase even at high temperatures.
It is interesting to note that at high temperatures many stopping points of the quench runs are actually saddle points (upon a small shift away from such a quench point a subsequent gradient minimization led to the nearest local minimum).In particular, for T ≥ 2000 K, the energy differences between the holding points and these quench points are close to the vibrational energy at the quench points, suggesting that the description of dynamics above the glass transition temperature as a hopping dynamics on saddle points [46] has some justification.Typical vibrational spectrum of a-Si 3 B 3 N 7 for holding points along tr fferent temperatures between 25 K and 4000 K. 'Negative' frequencies co ary modes.
The reason for this might be connected to the trapping behaviour of the landscape, since the local density of states sampled during the simulations appears to grow approximately exponentially in the low-energy regions of the landscape, leading to the so-called exponential trapping [47,48]: Due to the local exponential growth by a growth factor 1/T trap in the number of states N(E) = N(E 0 ) exp((E − E 0 )/T trap ), the walker cannot enter this region for temperatures T > T trap .Thus he spends all the time in the flat parts of a band of the landscape of width k B T below the top of the exponentially growing region (these are precisely the saddle points in the accessible energy range), and the system can easily move between basins on the landscape, since the effective barriers that need to be crossed are very small (i.e. of order k B T).This is consistent with the earlier observation that the walkers spend most of their time in individual basins and not in transition regions connecting them: if the effective barriers are very small, the walkers are highly mobile even though the actual transition regions are relatively small.On the other hand, for T < T trap , the walker drops into the deep exponential valleys and thus must cross rather high barriers for further structural re-arrangements.As a consequence, the system requires considerably longer for moving among basins on the landscapes, and this is reflected in the rather sudden freezing-in of the amorphous structure associated with the glass transition.
Taken together, these results suggest that a-Si 3 B 3 N 7 has a rugged energy landscape with a "basin"-based hierarchy, but where for each basin both the local density of states and the local density of minima grow very rapidly, possibly exponentially for the accessible density of states 1 .In addition, the energy differences of the saddle points to their "associated minima" also show a steady increase with the energy of the minima, but probably not so fast as the density of states itself.Finally, the associated flat transition regions between the basins slowly change in complexity with increasing energy.But even here, most of the movement is not diffusive but some kind of oscillatory motion, in agreement with the fact that up to about 4000 K the system is still in a condensed (liquid) phase.Only for even higher temperatures do we explore the part of the landscape that describes a fluid of smaller clusters of sizes below 50 atoms.

IV. Conclusions
In order to understand the properties of amorphous silicon boron nitride in greater depth, we have modeled its structure, computed its properties and analyzed its underlying energy landscape.We find that it is necessary to follow the actual synthesis route when generating a realistic structural model of a-Si 3 B 3 N 7 .This consists of a random network of SiN 4 -tetrahedra and BN 3 -and NSi x B 3-x -triangles, where Si and B show a moderate preference (due to the synthesis route) for Si and B as first cationic neighbours, respectively.The presence of a large number of subnanometersize voids in the structure accounts for the low density of the material.The physical properties of these models show good agreement with experimental data.From the computed phase diagram of the system one sees that the thermodynamically stable state is a phase separation into BN and Si 3 N 4 crystallites that is preferable even over the best metastable ternary crystalline Si 3 B 3 N 7 -modifications once the crystallite size exceeds a few dozen nanometers, explaining the failure of synthesizing a-Si 3 B 3 N 7 by sintering of BN and Si 3 N 4 micro-crystals.As a consequence, at sufficiently high temperatures and long enough time scales, formation of such crystallites is expected to occur in the amorphous phase.But already before this happens, aging will take place in the ceramics resulting in a very slow coarsening of the distribution of voids.This is in good agreement with the complex energy landscape of a-Si 3 B 3 N 7 , which most likely consists of a hierarchy of minima basins, with a quasi-exponential growth of the local densities of states in the low-energy range of the landscape.

Figure 1 .
Figure 1.Structure of a-Si 3 B 3 N 7 generated via simulation of the polymer-precursor route depicted as ball models with an edge length of ≈ 4.5 nm.[50] Si-atoms: red; B-atoms: blue; N-atoms: green.Note the inhomogeneous distribution of Si and B atoms.

Figure 2 .
Figure 2. Typical vibrational spectrum of a-Si 3 B 3 N 7 for local minima of the energy landscape observed along trajectories for five different temperatures between 25 K and 4000 K

Figure 3 .
Figure 3. Heat conductivity κ of a-Si 3 B 3 N 7 as function of temperature T.

Figure 3 :
Figure 3: Heat conductivity κ of a-Si 3 B 3 N 7 as function of temperature T .

Figure 4 .Figure 4 :
Figure 4. Metastable phase diagram of (amorphous) silicon boron nitride [21].G: gaseous phase; L: liquid phase; S: solid phase; L-G: gas-liquid co-existence region; S-L: gas-solid co-existence region; S-L: crystalline solid / super-cooled liquid.Note that each data point corresponds to an average of several simulations performed with a finite simulation cell at the given values of T and ρ.

Figure 5 .
Figure 5. Excess enthalpies for the quasi-binary Si 3 N 4 /BN system as function of composition, representing the (metastable) phase diagram at 0 K.The vertical axis depicts the excess enthalpy at standard pressure of the metastable crystalline structure candidates for various compositions compared to the weighted average of the enthalpy of Si 3 N 4 and BN.Only candidates with lowest energy are shown for each composition.

Figure 5 :
Figure 5: Excess enthalpies for the quasi-binary Si 3 N 4 /BN system as function of composition, representing the (metastable) phase diagram at 0 K.The vertical axis depicts the excess enthalpy at standard pressure of the metastable crystalline structure candidates for various compositions compared to the weighted average of the enthalpy of Si 3 N 4 and BN.Only candidates with lowest energy are shown for each composition.

Figure 6 .
Figure 6.Typical vibrational spectrum of a-Si 3 B 3 N 7 for holding points along trajectories for five different temperatures between 25 K and 4000 K. 'Negative' frequencies correspond to imaginary modes